****************************************************** *** ПОДГОНКА МОДЕЛЕЙ С ИЗБЫТОЧНОЙ ДИСПЕРСИЕЙ *** *** (распределение Пуассона с экстра-вариацией) *** ****************************************************** * Подобные модели довольно полезны в биологии при изучении популяций паразитов * ****************************************************** ************************************************************************************************** * По 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. ********************** * КОНЕЦ СИНТАКСИСА.