1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
******************************************************
***    ПОДГОНКА МОДЕЛЕЙ С ИЗБЫТОЧНОЙ ДИСПЕРСИЕЙ    ***
***  (распределение Пуассона с экстра-вариацией)   ***
******************************************************
* Подобные модели довольно полезны в биологии при изучении популяций паразитов  *
******************************************************

**************************************************************************************************
* По MJ Moroney 'Facts From Figures'                                                             *
* (London: Penguin Books, 1951; 2nd ed 1953)                                                     *
*                                                                                                *
* В пуассоновских моделях E(x)=V(x)=lambda (только один параметр)                                *
* Функция вероятности:                                                                           *
* P(x)=e^(-lambda)*(lambda^x)/x!                                                                 *
*                                                                                                *
* В распределениях с избыточной дисперсией V(X)>E(x)                                             *
* Отрицательный биномиальный закон распределения - прототип распределений с избыточной дисперсией*
*                                                                                                *
* В моделях с отрицательным биномиальным законом требуется дополнительный параметр (k):          *
* Параметр                     Чем оценивается:                                                  *
* E(x)=lambda ----------------> выборочное среднее                                               *
* V(x)=lambda*(1+lambda/k) ---> выборочная дисперсия                                             *
* k = EІ(x)/(V(x)-E(x))      -> выб. средн.^2/(выб. дисп. - выб. средн.) ('метод моментов')      *
* Комментарий: метод максимального правдоподобия лучше, но гораздо более трудоёмкий!             *
* (см. веб и программу ParaDis)                                                                  *
* Функция вероятности:                                                                           *
* P(x)=[(k/(k+lambda))^k]*[(lambda/(k+lambda))^x]*GAMMA(x+k)/(x!*GAMMA(k))                       *
**************************************************************************************************
 
* Пример набора данных (из руководства к программе ParaDis):
* (http://www.bondy.ird.fr/~pichon/paradis/parad2.html).

DATA LIST FREE/ number freq (2 F8.0).
BEGIN DATA
0 70
1 38
2 17
3 10
4  9
5  3
6  2
7  1
END DATA.
VARIABLE LABEL number 'Клещи/лист растения'.
 
* Добавим среднюю, дисперсию и число наблюдений к данным и проведём проверку Фишера на избыточность вариации *.
CACHE.
EXECUTE.
MATRIX.
PRINT /TITLE='ДАННЫЕ И СТАТИСТИКИ'.
GET x/VAR=number.
GET freq/VAR=freq.
COMPUTE ndata=NROW(x).
COMPUTE nt=MSUM(freq).
COMPUTE mean=MSUM(freq&*x)/nt.
COMPUTE variance=(MSUM(freq&*(x&**2))-nt*mean&**2)/(nt-1).
COMPUTE k=(mean&**2)/(variance-mean).
PRINT {x,freq}
 /FORMAT='F8.0'
 /TITLE='Данные выборки'
 /CLABELS='X','Част'.
PRINT nt
 /FORMAT='F8.0'
 /TITLE='Общий объём выборки'
 /RLABEL='Nt'.
PRINT {mean,variance,(variance/mean),k}
 /FORMAT='F8.2'
 /TITLE='Статистики'
 /CLABELS='Среднее','Дисперсия','DP (**)','K(*)'.
PRINT /TITLE='(*) Подсчитан по методу моментов.'.
PRINT /TITLE='(**) Отношение дисперсии к средней.'.
PRINT /TITLE='ПРОВЕРКА ФИШЕРА НА ИЗБЫТОЧНОСТЬ ДИСПЕРСИИ'.
COMPUTE fd=(MSUM(freq&*(x-mean)&**2))/mean.
COMPUTE fdsig=1-CHICDF(fd,MSUM(freq)-1).
PRINT {fd,fdsig}
 /FORMAT='F8.4'
 /TITLE='Статистика проверки (df=Nt-1)'
 /CLABELS='Хи^2','Знач.'.
DO IF fdsig LE 0.05.
- PRINT /TITLE='Присутствует избыточная дисперсия! Обнаружена экстра-вариация в распр. Пуассона'.
END IF.
COMPUTE means=MAKE(ndata,1,mean).
COMPUTE n=MAKE(ndata,1,nt).
COMPUTE kp=MAKE(ndata,1,k).
COMPUTE namevec={'number','n','mean','k'}.
SAVE {x,n,means,kp} /OUTFILE='c:\\temp\\temp.sav' /NAMES=namevec.
END MATRIX.
MATCH FILES /FILE=*
 /FILE='C:\\Temp\\temp.sav'
 /BY number.
EXECUTE.
 
* Подсчёт ожидаемых частот в предположении отрицательного биномиального распределения *
  (Для критеря хи-квадрат происходит потеря 3-х степеней свободы: используется N, среднее и дисперсия) *.
COMPUTE expect=n*((mean/(k+mean))**number)*((k/(k+mean))**k)*EXP(LNGAMMA(number+k))/(EXP(LNGAMMA(number+1))*EXP(LNGAMMA(k))).
COMPUTE sumexp=expect.
* Небольшая хитрость: делаем сумму ожидаемых частот равной N, корректируя ожидаемую частоту последней категории *.
DO IF $casenum GT 1.
- COMPUTE sumexp=sumexp+LAG(sumexp).
END IF.
SORT CASES BY number(D).
DO IF $casenum EQ 1.
- COMPUTE expect=expect+(n-sumexp).
END IF.
SORT CASES BY number(A).
 
* Табличный и графический вывод наблюдаемых и ожидаемых частот *.
VAR LABEL freq 'Наблюдаемые' /expect 'Ожидаемые'.
REPORT
 /FORMAT=LIST
 /TITLE='Наблюдаемые и ожидаемые частоты (с учётом избыточной дисперсии)'
 /VAR=number freq expect.
GRAPH /BAR(GROUPED)=VALUE(freq expect) BY number
 /TITLE='Наблюдаемые и ожидаемые(*) частоты'
 /Footnote='(*) С учётом избыточной дисперсии'.
 
* Проверяем, нет ли малых ожидаемых частот и объединяем ячейки, чтобы избавиться от таких частот *.
COMPUTE id=$casenum.

DO IF (expect LT 2).
- COMPUTE id=LAG(id).
END IF.
AGGREGATE
 /OUTFILE=*
 /BREAK=id
 /observed = SUM(freq)
 /expected = SUM(expect).
STRING groups (A2).
COMPUTE groups = STRING(id,F2.0) .
 
* Проверка согласия по критерию хи-квадрат *.
MATRIX.
PRINT /TITLE='ПРОВЕРКА СОГЛАСИЯ ХИ-КВАДРАТ'.
GET groups/VAR=groups.
GET obs/VAR=observed.
GET expect/VAR=expected.
PRINT {obs,expect,(obs-expect);MSUM(obs),MSUM(expect),0}
 /FORMAT='F8.2'
 /TITLE='Частоты (после объединения категорий)'
 /CLABELS='Наблюдаемые','Ожидаемые','Разность'
 /RNAMES=groups.
COMPUTE k=NROW(obs).
PRINT {k-3}
 /FORMAT='F8.0'
 /TITLE='Число степеней свободы (k-3)'.
COMPUTE chi2=MSUM((obs-expect)&**2/expect).
COMPUTE chisig=1-CHICDF(chi2,k-3).
PRINT {chi2,chisig}
 /FORMAT='F8.4'
 /TITLE='Статистика проверки'
 /CLABELS='Хи^2','Знач.'.
COMPUTE minexp=CMIN(expect).
COMPUTE flag=0.
LOOP i=1 TO k.
- DO IF expect(i) LT 5.
-  COMPUTE flag=flag+1.
- END IF.
END LOOP.
COMPUTE pflag=100*flag/k.
DO IF flag GT 0.
- PRINT pflag
 /FORMAT='F8.1'
 /TITLE='ВНИМАНИЕ: присутствуют ячейки с ожидаемой частотой, меньшей, чем 5.'
 /RLABEL='% таких ячеек='.
- PRINT minexp
 /FORMAT='F8.1'
 /TITLE='Минимальная ожидаемая частота ячейки:'
 /RLABEL='Ожид.='.
END IF.
END MATRIX.

**********************
* КОНЕЦ СИНТАКСИСА.