DEFINE !BINCINT (N !TOKENS(1) / X !TOKENS(1) /CONF !TOKENS(1) ).
*.
* .
* По данным результата биномиального эксперимента, X = число успехов и N = число попыток,
* эта программа возвращает точный (неасимптотический) доверительный интервал.
* .
title 'Точный доверительный интервал для биномиального параметра'.
* .
* Пример использования:  !BINCINT N=10 X=3 CONF=95.0 .
*.
* Автор: Mike Lacy, декабрь 1996.
* Инициализируем файл данных, чтобы можно было производить расчёты в SPSS.
new file .
input program.
compute dummy = 1 .
end case.
end file.
end input program.
*  .
***********************************************.
* Константы, изменяемые пользователем .
compute toler = 1.0e-4.              /* желаемая относительная точность. 
compute maxiter = 40.                /* максимальное число итераций метода.
* .
string lconverg (a3) uconverg (a3) . /* флаги для индикации схождения алгоритма.
* .
* Принимаем параметры макроса и проводим необходимые вычисления.
compute N = !N .
compute x = !X .
compute alpha2 =  (1 - (!CONF/100))/2 .
* Некоторые дополнительные вычисления.
compute pyes = x/n.      
compute se = sqrt(pyes * (1-pyes)/N).
*.
*.
*******************************************.
*  Найдём границы доверительного интервала. Сначала - нижнюю, затем - верхнюю.
* .
loop dolow = 1 to 0  by -1.
do if dolow .
* Инициализируем нижнюю границу.
+  compute lowbound = 0.0 . /*ограничим снизу нижнюю границу.
+  compute upbound  = pyes . /*ограничим сверху нижнюю границу
*  Сгенерируем два приближения нижней границу.  
+  compute OldGuess = pyes.
+  compute guess = pyes - probit(alpha2) * se.
+  if (guess >= upbound) or (guess <= lowbound) 
      guess = (lowbound + upbound)/2.0.
*   Вычислим вероятности для начальных приближений.
+   compute p = 1 - cdf.binom(x-1,n,guess).
+   compute OldP = 1- cdf.binom(x-1,n,OldGuess).
+   compute lconverg = 'no'.
*.
else.  /* Верхняя граница .
*.
*  Инициализируем верхнюю границу.
+  compute lowbound = pyes . /*ограничим снизу верхнюю границу.
+  compute upbound  = 1.0 . /*ограничим сверху верхнюю границу
*  Сгенерируем два приближения верхней границы.  
+  compute OldGuess = pyes.
+  compute guess = pyes + probit(1-alpha2) * se.
+  if (guess >= upbound) or (guess <= lowbound) 
      guess = (lowbound + upbound)/2.0.
*  Вычислим вероятности для начальных приближений.
+  compute p = cdf.binom(x,n,guess).
+  compute OldP = cdf.binom(x,n,OldGuess).
+  compute uconverg = 'нет'.
* .
end if.   /* Конец инициализации .
*.
*.
*. Продолжаем искать границы, проверим для начала на частные случаи.
*.
do if (x = 0) and dolow .      /* Число успехов (X) минимально - интервал не строим.
+  compute lcl = 0.0.
+  compute lconverg = 'да'.
else if (x=N) and not dolow.   /* Число успехов (X) максимально - интервал не строим.
+  compute ucl = 1.0 .
+  compute uconverg = 'да' . 
else .
*  Используем итеративный метод для поиска границ.
* Используем метод секущей (эмпирическая производная) с поиском методом деления.
*. 
+  compute numiter = 0 .
+  compute diff = p - alpha2 .
*.
*          ЦИКЛ ИТЕРАЦИЙ .
+  loop .  
+    compute numiter = numiter + 1 .
*    Установка исходных границ для итерации: разная для нижней и верхней границ.
+    do if (diff < 0) and dolow .  
+      compute lowbound = guess.
+    else if (diff <0) and not dolow. 
+      compute upbound = guess .
+    else if (diff >= 0) and dolow .
+      compute upbound = guess .
+    else if (diff >= 0) and not dolow .
+      compute lowbound = guess .
+    end if .
* .
*    Формируем следующее приближение и его вероятность, отсеивая плохие приближения и
* большие наклоны секущей.
+    compute slope = (OldP - P)/(OldGuess - Guess) .
+    compute OldP = P .
+    compute OldGuess = Guess .
+    do if ((Abs(slope) < 1e-5 ) or (Abs(slope) > 1e+5)) .
+      compute guess = (upbound + lowbound)/2.0 .
+    else .
+      compute guess = guess - diff/slope .
+      if (guess >= upbound) or (guess <= lowbound) 
         guess = (lowbound + upbound)/2.0.
+    end if .
+    do if dolow.
+      compute p = 1- cdf.Binom(x-1,n,Guess) .
+    else if not dolow .
+      compute p = cdf.Binom(x,n,Guess) .
+    end if.
+    compute diff = p - alpha2 .
*.
+  end loop if ( (abs(diff) <= toler * alpha2)
         or    (numiter > maxiter) ) .
* .
end if . /* если удалось достигнуть границы в очередной итерации .
*.
do if dolow .
+   compute lcl = guess .
+   compute l_iter = numiter .
+   if (numiter < maxiter) lconverg = 'да'.
else if not dolow .
+   compute ucl = guess .
+   compute u_iter = numiter .
+   if (numiter < maxiter) uconverg = 'да'.
end if . /* конец верхнего условия dolow .
end loop . /* цикл по dolow = 1 до 0 для нижней и верхней границ.
*.
format pyes (f7.5) ucl (f7.5) lcl (f7.5)  
       x (f6.0)  n (f6.0)  l_iter (f6.0) u_iter (f6.0).
list x n lcl pyes ucl l_iter lconverg u_iter uconverg.
!ENDDEFINE .
* .
!BINCINT N=10 X=3 CONF=95.0 .
* .
* X = наблюдавшееся число успехов в эксперименте. 
* N = число попыток .
* lcl = нижняя граница доверительного интервала.
* pyes = наблюдаемая доля успехов (параметр биномиального распределения).
* ucl = верхняя граница доверительного интервала.
* l_iter = число итераций, потребовавшихся для нахождения последнего приближения нижней границы.
* l_converg = была ли достигнута требуемая точность для нижней границы?.
* 'u_iter = число итераций, потребовавшихся для нахождения последнего приближения верхней границы.
* 'u_converg = была ли достигнута требуемая точность для верхней границы?.
* .
* .
exe.