**** КРИТЕРИЙ СОГЛАСИЯ ДЛЯ РАСПРЕДЕЛЕНИЯ ПУАССОНА ****

* Пример данных *.

DATA LIST FREE/ NUMBER(f8.0).
BEGIN DATA
0 2 0 0 2 2 0 0 1 1 3 0 0 1 0
0 1 0 1 4 0 0 1 4 2 0 0 1 0 0
2 2 1 1 0 6 0 5 1 3 0 1 0 1 8
END DATA.
VAR LABEL number 'Число иксодовых клещей/мышь'.

* Проверка по критерию Колмогорова-Смирнова для одной выборки (низкая чувствительность
* если вручную не ввести поправку Лиллефорса (Lilliefors correction) *.

NPAR TESTS
  /K-S(POISSON)= number
  /MISSING ANALYSIS.

* Обратите внимание на пустую категорию между 6 и 8; она должны быть заполнена, чтобы 
* ожидаемые значения были подсчитаны верно (данная часть приводится только
* из образовательных соображений, она может быть удалена из практически используемого кода *.

FREQUENCIES
  VARIABLES=number
  /ORDER=  ANALYSIS .

*** Заполняем пустые категории ***

* Во-первых, сгруппируем наблюдения по частоте встречаемости вариат *.
AGGREGATE
  /OUTFILE=*
  /BREAK=number
  /obs=N.
EXEC.

* Создадим вспомогательный файл со всевозможными значениями от 0 до максимального значения из файла данных *.
PRESERVE.
SET ERRORS=NONE RESULTS=NONE.
MATRIX.
GET number /VAR=number.
COMPUTE max=CMAX(number).
RELEASE number.
COMPUTE number=MAKE(max+1,1,1).
COMPUTE wgt=MAKE(max+1,1,0.0001).
LOOP i=1 TO max+1.
- COMPUTE number(i)=i-1.
END LOOP.
COMPUTE namevec={"number","wgt"}.
SAVE {number,wgt} /OUTFILE='c:\\temp\\temp.sav' /NAMES=namevec.
END MATRIX.

* Объединим его с рабочим файлом для заполнения пустых категорий *.
MATCH FILES /FILE=*
 /FILE='C:\\Temp\\temp.sav'
 /BY NUMBER.
EXECUTE.
IF MISSING(obs) obs=wgt.
MATCH FILES /FILE=*
 /DROP=wgt.
EXECUTE.
WEIGHT BY obs.
RESTORE.
* Теперь пустоты заполнены как бы нулевой частотой
*(на самом деле, пустая категория имеет малый вес (данные взвешены) - А.Б.) *.
FREQUENCIES
  VARIABLES=number
  /ORDER=  ANALYSIS .

*** Критерий согласия хи-квадрат (с большей чувствительностью, чем проверка Колмогорова-Смирнова) ***

* Рассчитаем среднее и количество наблюдений (N) *.
COMPUTE const = 1 .
EXECUTE .
AGGREGATE
  /OUTFILE='C:\\temp\\aggr.sav'
  /BREAK=const
  /mean = MEAN(number)
  /N=N.
WEIGHT OFF.
PRESERVE.
SET ERRORS=NONE.
MATCH FILES /FILE=*
 /FILE='C:\\temp\\aggr.sav'
 /BY const.
EXECUTE.
RESTORE.
DO IF MISSING(mean).
- COMPUTE mean=LAG(mean).
- COMPUTE n=LAG(n).
END IF.
EXECUTE.

* Рассчитаем ожидаемые частоты *.
DO IF $casenum=1.
- COMPUTE expect=n*CDF.POISSON(number,mean) .
END IF.
DO IF $casenum GT 1.
- COMPUTE expect=n*CDF.POISSON(number,mean)-n*CDF.POISSON(LAG(number),mean) .
END IF.
SORT CASES BY number(D).
DO IF $casenum=1.
- COMPUTE expect=n*(1-CDF.POISSON(number-1,mean)).
END IF.
EXEC.
SORT CASES BY number(A).

* Графическое сопоставление перед объединением категорий.
* Хотя оно не является частью проверки, бывает полезно для визуализации отклонений от
* распределения Пуассона. Должно производиться перед объединением категорий *.
GRAPH /BAR(GROUPED)=VALUE( obs expect ) BY number .

* Проверяем, какие категории имеют ожидаемые частоты, меньшие, чем 1.
* Объединяем соответствующие категории.
* (можно доработать для объединения категорий с ожидаемыми частотами, меньшими 5) *.
COMPUTE id=$casenum.
DO IF expect LT 1.
- COMPUTE id=LAG(id).
END IF.
EXECUTE.
AGGREGATE
  /OUTFILE=*
  /BREAK=id
  /observed = SUM(obs) /expected = SUM(expect).

* Статистика проверки и её значимость *.
STRING groups (A2).
COMPUTE groups = STRING(id,F2.0) .
EXECUTE .
MATRIX.
GET groups/VAR=groups.
GET obs/VAR=observed.
GET expect/VAR=expected.
PRINT {obs,expect;MSUM(obs),MSUM(expect)}
 /FORMAT='F8.2'
 /TITLE='Наблюдаемые и ожидаемые частоты'
 /CLABELS='Набл.','Ожид.'
 /RNAMES=groups.
COMPUTE k=NROW(obs).
PRINT {k-2}
 /FORMAT='F8.0'
 /TITLE='Число степеней свободы (k-2)'.
COMPUTE chi2=MSUM((obs-expect)&**2/expect).
COMPUTE chisig=1-CHICDF(chi2,k-2).
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.