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
ПРИМЕР:
В выборке 400 респондентов имеется:
(1)  60 голосов за "красных"
(2) 120 голосов за "синих"
(3)  20 голосов за "зелёных"
(4) 200 голосов за "серых"
(5)   0 голосов за "чёрных"
т.е. никто не "проголосовал" за "чёрных", но мы точно знаем, что такая партия существует.
Тогда байесовские оценки долей и 99% доверительные интервалы  могут быть подсчитаны как:

bayes 60 120 20 200 0 / 99%.

В случае 95% доверительного интервала, последний аргумент может быть опущен:

bayes 60 120 20 200 0.

**************************************************
* Макрос: Bayes.sps                              *
* Автор:  Евгений Ивашкевич (Evgeny Ivashkevich) *
* E-mail: ivashkev@yandex.ru)                    *
* Дата  : 2005/07/19                             *
**************************************************

preserve.
set printback=off.
define bayes  (!positional !charend ('/')
              /!positional !charend ('%') !default(95)).

!let !n = 0.
!let !r = 0.
!do !k !in (!1).
!let !n = !length(!concat(!blanks(!n),!blanks(!k))).
!let !r = !length(!concat(!blanks(!r),!blanks(1))).
!doend.

input program.
numeric id (f2) pr se lb ub (f7.5).
vector cov(!r e10.3) #p(!r).
string info(A10).
variable labels  info ' '
 id '#'
 pr 'Доля'
 se 'Стандартная ошибка'
 lb 'Нижняя граница'
 ub 'Верхняя граница'.

compute #i = 0.
!do !k !in (!1).
  compute #i = #i+1.
  compute #p(#i) = (!k+1)/(!n+!r).
!doend.

loop #i=1 to !r.
  compute id = #i.
  compute pr = #p(#i).
  compute se = sqrt(#p(#i)*(1-#p(#i))/(!n+!r+1)).
  compute lb = max(#p(#i)+probit((100-!2)/200)*se,0).
  compute ub = min(#p(#i)+probit((100+!2)/200)*se,1).
  loop #j=1 to !r.
    if (#j =#i) cov(#j) =  #p(#i)*(1-#p(#i))/(!n+!r+1).
    if (#j<>#i) cov(#j) = -#p(#i)*#p(#j)/(!n+!r+1).
  end loop.
  if ($casenum=1) info = concat('P =',string(!2,f4),'%').
  if ($casenum=2) info = concat('N =',string(!n,f5)).
  if ($casenum>2) info = ' '.
end case.
end loop.
end file.
end input program.
execute.

report format = list automatic align(center)
 /variables = id pr se lb ub info
 /title 'БАЙЕСОВСКИЕ ОЦЕНКИ ДОЛЕЙ И ДОВЕРИТЕЛЬНЫЕ ИНТЕРВАЛЫ'.
!enddefine.
restore.

bayes 60 120 20 200 0 / 99%.