Доверительный интервал для дисперсии методом бутстреп

PRESERVE.
SET printback=off.
DEFINE !BCNON(             
               SAMPLES=!TOKENS(1) !DEFAULT(1000)
              /VAR    =!TOKENS(1)
              /ALPHA1 =!TOKENS(1) !DEFAULT(0.05)
              /ALPHA2 =!TOKENS(1) !DEFAULT(0.95)
              /DEBUG  = !CMDEND   !DEFAULT('N')   ).
****************************************************************************
* Этот макрос строит непараметрический дверительный интервал на основе бутстреп-метода
* для дисперсии переменной из открытого файла данных. 
* Этот метод описан в книге "An Introduction to the Bootstrap" (Введение в бутстреп-оценки) B Efron
* и R J Tibshirani, Chapman & Hall, 1993, С. 184-188, (BC method).
* Я признателен корп. SPSS за основу кода (см.
* http://www.spss.com/tech/answer/result.cfm?tech_tan_id=100000795 )
* и Raynald Levesque за предложенный метод включения числа наблюдений в счётчик цикла LOOP.
*
* Стандартное использование этого кода: 
* !BCNON SAMPLES=1000 VAR=hhinc ALPHA1=0.05 ALPHA2=0.95.
* где
* SAMPLES=1000 число бутстреп-выборок, которые необходимо сделать
* VAR имя переменной, для дисперсии которой необходимо построить доверительный интервал
* ALPHA1 и ALPHA2  -  доверительные уровни, задающие нижнюю и верхнюю границу доверительного
* интревала (есть возможность построения несимметричного интервала).
*
* VAR - единственный обязательный параметр. Если остальные папаметры опущены, их значения 
* будут приняты по умолчанию. Т.е.
* !BCNON VAR=hhinc.
* есть минимальная команда для вызова этого макроса.
*
* David Hitchin, Университет Суссекса, 7 марта 2001
* ccfg4@sussex.ac.uk
* 
****************************************************************************.
PRESERVE.
!IF ( !DEBUG !EQ 'N') !THEN
SET printback=off mprint off.                                                   
!ELSE
SET printback on mprint on.
!IFEND .
SET SEED=random.
*---------------------------------------------------------------------------.
* Сохраняем исходный файл для восстановления его после работы макроса.
*---------------------------------------------------------------------------.
!IF (!DEBUG !EQ 'N') !THEN 
SET RESULTS ON.
DO IF $CASENUM=1.
PRINT / "Внимание вывод всех результатов, включая сообщения об ошибках, временно подавлен" 
      / "Если есть проблемы, перезапустите макрос, установив параметр DEBUG в Y"  
      / "Перед этим загрузите исходный файл данных"
      / "Так вы будете получать отчёт о всех ошибках".
END IF.
!IFEND .
SAVE OUTFILE='C:\\Temp\\rr__tmp1.sav'.

SELECT IF (NOT MISSING(!VAR)).
COMPUTE sample=1.
COMPUTE id=$casenum.

* Небольшое мошенничество для того, чтобы загнать число наблюдений в файле данных в цикл.
RANK VARIABLES=!VAR  (A) /N INTO N.
EXECUTE.
DO IF $casenum=1.
WRITE OUTFILE='C:\\Temp\\syntax.sps'  /'DEFINE !NBCASES()'N '!ENDDEFINE.'.
END IF.
EXECUTE.
INCLUDE FILE='C:\\Temp\\syntax.sps' .

* DESCRIPTIVES  VARIABLES=!VAR /STATISTICS=VARIANCE .
SAVE OUTFILE='C:\\Temp\\temp1.sav'/keep=!VAR,sample,id.
GET FILE='C:\\Temp\\temp1.sav'.
AGGREGATE  /OUTFILE='C:\\Temp\\sd.sav'/BREAK=sample/s = SD(!VAR).

* Непосредственно метод бутстреп начинается здесь.
INPUT PROGRAM .
LOOP SAMPLE=1 to !SAMPLES.
  LOOP N = 1 to !NBCASES.
    COMPUTE ID=TRUNC(UNIFORM(!NBCASES)) + 1.
    END CASE.
    LEAVE SAMPLE.
    END LOOP.
END LOOP.
END FILE.
END INPUT PROGRAM .

FORMATS ID SAMPLE N (F6.0).
SORT CASES BY ID .
MATCH FILES / FILE = * / TABLE 'C:\\Temp\\temp1.sav' / BY ID.
AGGREGATE  /OUTFILE=*/BREAK=sample/z = SD(!VAR).
COMPUTE z=z*z.
COMPUTE sample=1.
MATCH FILES / FILE = * / TABLE 'C:\\Temp\\sd.sav' / BY sample.

COMPUTE theta=s*s.
COMPUTE thetals=(z < theta).
AGGREGATE  /OUTFILE='C:\\Temp\\theta.sav'/ BREAK=sample / stless=SUM(thetals) / thetabar=MEAN(z) / N_BREAK=N.
MATCH FILES / FILE = * / TABLE 'C:\\Temp\\theta.sav' / BY sample.

COMPUTE z0=stless/n_break.
COMPUTE z0hat=IDF.NORMAL(z0,0,1).

COMPUTE zalph=IDF.NORMAL(!ALPHA2,0,1).
COMPUTE z1malph=IDF.NORMAL(!ALPHA1,0,1).

COMPUTE t2=(thetabar-theta)*(thetabar-theta).
COMPUTE t3=t2*(thetabar-theta).
AGGREGATE  /OUTFILE='C:\\Temp\\tt.sav'/ BREAK=sample / tt2=SUM(t2) / tt3=SUM(t3).
MATCH FILES / FILE = * / TABLE 'C:\\Temp\\tt.sav' / BY sample.

COMPUTE alphahat=tt3/ (6*tt2)**1.5 .
COMPUTE a=z0hat+(z0hat+z1malph)/(1-alphahat*(z0hat+z1malph)).
COMPUTE a1=CDF.NORMAL(a,0,1).
COMPUTE a=z0hat+(z0hat+zalph)/(1-alphahat*(z0hat+zalph)).
COMPUTE a2=CDF.NORMAL(a,0,1).

COMPUTE pl=a1*100.
COMPUTE pu=a2*100.
COMPUTE cl=RND(n_break*a1).
COMPUTE cu=RND(n_break*a2).
PRINT FORMATS z to theta,thetabar,z0 to pu (f6.4)/
 thetals,stless,cl,cu,n_break (f6.0)/ sample (f2.0).

COMPUTE filter_$=($casenum=1).
VARIABLE LABEL filter_$ '$casenum=1 (FILTER)'.
VALUE LABELS filter_$  0 'Не выбран' 1 'Выбран'.
FORMAT filter_$ (f1.0).
FILTER BY filter_$.
*DESCRIPTIVES VARIABLES=sample z s theta thetals stless thetabar n_break z0 z0hat zalph
  z1malph t2 t3 tt2 tt3 alphahat a a1 a2 pl pu cl cu
  /STATISTICS=MEAN STDDEV MIN MAX .
FILTER OFF.
USE ALL.

SORT CASES BY  z (A).
IF ($casenum=cl) lci=z.
IF ($casenum=cu) uci=z.
AGGREGATE  /OUTFILE=* /BREAK=sample /lci_1 = MIN(lci) /uci_1 = MAX(uci).
RENAME VARIABLES  (lci_1 =lci)(uci_1=uci).
* SUMMARIZE  /TABLES=lci uci  /FORMAT=VALIDLIST NOCASENUM TOTAL LIMIT=100
*  /TITLE='Case Summaries'  /MISSING=VARIABLE  /CELLS=COUNT .
MATRIX.
GET m /VARIABLES=lci,uci.
GET s /FILE='C:\\Temp\\sd.sav' /VARIABLES=s.
COMPUTE s=s*s.
PRINT s /FORMAT="F14.6"/TITLE="Дисперсия"/CLABELS=" ".
PRINT m /FORMAT="F14.6"/TITLE="Доверительные границы для дисперсии"/CLABELS="Нижняя","Верхняя".
COMPUTE lower=!ALPHA1.
COMPUTE upper=!ALPHA2.
PRINT lower/FORMAT="F7.3"/TITLE="Alpha1"/CLABELS=" ".
PRINT upper/FORMAT="F7.3"/TITLE="Alpha2"/CLABELS=" ".
END MATRIX.
ERASE FILE='C:\\Temp\\temp1.sav'.
ERASE FILE='C:\\Temp\\sd.sav'.
ERASE FILE='C:\\Temp\\theta.sav'.
ERASE FILE='C:\\Temp\\tt.sav'.
ERASE FILE='C:\\Temp\\syntax.sps'.
GET FILE='C:\\Temp\\rr__tmp1.sav'.

RESTORE.
!ENDDEFINE.
RESTORE.

* Проверочные данные из книги Efron, стр. 180. Опыт не обязательно точно повторит результаты
* из книги, т.к. будут осуществлены разные случайные выборки.
Data list free/zz.
Begin data.
48 36 20 29 42 42 20 42 22 41 45 14 6 
0 33 28 34 4 32 24 47 41 24 26 30 
end data.
!BCNON SAMPLES=2000 VAR=zz ALPHA1=0.05 ALPHA2=0.95.