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 | 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.
|