|
1
|
|
|
# Copyright (c) 2014, Salesforce.com, Inc. All rights reserved. |
|
2
|
|
|
# |
|
3
|
|
|
# Redistribution and use in source and binary forms, with or without |
|
4
|
|
|
# modification, are permitted provided that the following conditions |
|
5
|
|
|
# are met: |
|
6
|
|
|
# |
|
7
|
|
|
# - Redistributions of source code must retain the above copyright |
|
8
|
|
|
# notice, this list of conditions and the following disclaimer. |
|
9
|
|
|
# - Redistributions in binary form must reproduce the above copyright |
|
10
|
|
|
# notice, this list of conditions and the following disclaimer in the |
|
11
|
|
|
# documentation and/or other materials provided with the distribution. |
|
12
|
|
|
# - Neither the name of Salesforce.com nor the names of its contributors |
|
13
|
|
|
# may be used to endorse or promote products derived from this |
|
14
|
|
|
# software without specific prior written permission. |
|
15
|
|
|
# |
|
16
|
|
|
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS |
|
17
|
|
|
# "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT |
|
18
|
|
|
# LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS |
|
19
|
|
|
# FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE |
|
20
|
|
|
# COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, |
|
21
|
|
|
# INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, |
|
22
|
|
|
# BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS |
|
23
|
|
|
# OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND |
|
24
|
|
|
# ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR |
|
25
|
|
|
# TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE |
|
26
|
|
|
# USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. |
|
27
|
|
|
|
|
28
|
|
|
import numpy as np |
|
29
|
|
|
from scipy.special import betaln, binom |
|
30
|
|
|
from matplotlib import pyplot |
|
31
|
|
|
|
|
32
|
|
|
|
|
33
|
|
|
def binomln(n, k): |
|
34
|
|
|
return np.log(binom(n, k)) |
|
35
|
|
|
|
|
36
|
|
|
|
|
37
|
|
|
def simone(alpha, beta, n, ITERS): |
|
38
|
|
|
ks = [] |
|
39
|
|
|
for i in range(ITERS): |
|
40
|
|
|
theta = np.random.beta(alpha, beta) |
|
41
|
|
|
ks.append(np.random.binomial(n, theta)) |
|
42
|
|
|
return np.array(ks) |
|
43
|
|
|
|
|
44
|
|
|
|
|
45
|
|
|
def simtwo(alpha, beta, n, ITERS): |
|
46
|
|
|
counts = np.zeros((ITERS, 2), dtype=np.uint32) |
|
47
|
|
|
thetas = np.random.beta(alpha, beta, size=ITERS) |
|
48
|
|
|
for i in range(ITERS): |
|
49
|
|
|
if i % 1000000 == 0: |
|
50
|
|
|
print i |
|
51
|
|
|
|
|
52
|
|
|
counts[i] = np.random.binomial(n, thetas[i], size=2) |
|
53
|
|
|
return counts |
|
54
|
|
|
|
|
55
|
|
|
|
|
56
|
|
|
class Hypers(object): |
|
57
|
|
|
|
|
58
|
|
|
def __init__(self): |
|
59
|
|
|
self.alpha = 0.0 |
|
60
|
|
|
self.beta = 0.0 |
|
61
|
|
|
self.N = 0 |
|
62
|
|
|
|
|
63
|
|
|
def __str__(self): |
|
64
|
|
|
return "alpha=%f, beta=%f, N=%d" % (self.alpha, |
|
65
|
|
|
self.beta, |
|
66
|
|
|
self.N) |
|
67
|
|
|
|
|
68
|
|
|
|
|
69
|
|
|
class SuffStats(object): |
|
70
|
|
|
|
|
71
|
|
|
def __init__(self): |
|
72
|
|
|
self.heads = 0 |
|
73
|
|
|
self.tails = 0 |
|
74
|
|
|
self.binomln_accum = 0 |
|
75
|
|
|
self.dpcount = 0 |
|
76
|
|
|
|
|
77
|
|
|
def __str__(self): |
|
78
|
|
|
return "heads=%d, tails=%d" % (self.heads, self.tails) |
|
79
|
|
|
|
|
80
|
|
|
|
|
81
|
|
|
def add_dp(k, hps, suffs): |
|
82
|
|
|
suffs.heads += k |
|
83
|
|
|
suffs.tails += hps.N - k |
|
84
|
|
|
suffs.binomln_accum += binomln(hps.N, k) |
|
85
|
|
|
|
|
86
|
|
|
suffs.dpcount += 1 |
|
87
|
|
|
|
|
88
|
|
|
|
|
89
|
|
|
def compute_post_pred(k, hps, suffs): |
|
90
|
|
|
|
|
91
|
|
|
a = hps.alpha + suffs.heads |
|
92
|
|
|
b = hps.beta + suffs.tails |
|
93
|
|
|
|
|
94
|
|
|
return binomln(hps.N, k) + betaln(a + k, b + hps.N - k) - betaln(a, b) |
|
95
|
|
|
|
|
96
|
|
|
|
|
97
|
|
|
def compute_total_likelihood(hps, suffs): |
|
98
|
|
|
|
|
99
|
|
|
a = hps.alpha + suffs.heads |
|
100
|
|
|
b = hps.beta + suffs.tails |
|
101
|
|
|
|
|
102
|
|
|
r = betaln(a, b) - betaln(hps.alpha, hps.beta) |
|
103
|
|
|
|
|
104
|
|
|
r += + suffs.binomln_accum |
|
105
|
|
|
|
|
106
|
|
|
return r |
|
107
|
|
|
|
|
108
|
|
|
|
|
109
|
|
|
def two_d_hist(alpha, beta, N, ITERS): |
|
110
|
|
|
|
|
111
|
|
|
r = simtwo(alpha, beta, N, ITERS) |
|
112
|
|
|
h, b1, b2 = np.histogram2d(r[:, 0], r[:, 1], bins=N + 1) |
|
113
|
|
|
h[:, :] += 1 |
|
114
|
|
|
pyplot.imshow(h) |
|
115
|
|
|
pyplot.show() |
|
116
|
|
|
|
|
117
|
|
|
logp = np.log(h.astype(float) / (ITERS)) |
|
118
|
|
|
|
|
119
|
|
|
return logp |
|
120
|
|
|
|
|
121
|
|
|
|
|
122
|
|
|
def test_post_pred_equal_likelihood(): |
|
123
|
|
|
hps = Hypers() |
|
124
|
|
|
hps.N = 10 |
|
125
|
|
|
|
|
126
|
|
|
for alpha in [5.0]: |
|
127
|
|
|
for beta in [5.0]: |
|
128
|
|
|
for ks in [[9, 1], [1, 9]]: |
|
129
|
|
|
print "-" * 60 |
|
130
|
|
|
hps.alpha = alpha |
|
131
|
|
|
hps.beta = beta |
|
132
|
|
|
|
|
133
|
|
|
ss = SuffStats() |
|
134
|
|
|
|
|
135
|
|
|
pred_score = 0.0 |
|
136
|
|
|
for ki, k in enumerate(ks): |
|
137
|
|
|
pred_score += compute_post_pred(k, hps, ss) |
|
138
|
|
|
add_dp(k, hps, ss) |
|
139
|
|
|
|
|
140
|
|
|
total_score = compute_total_likelihood(hps, ss) |
|
141
|
|
|
res = two_d_hist(alpha, beta, hps.N, 10000000) |
|
142
|
|
|
print res |
|
143
|
|
|
print res.shape |
|
144
|
|
|
emp_score = res[ks[0], ks[1]] |
|
145
|
|
|
d1 = abs(pred_score - emp_score) |
|
146
|
|
|
d2 = abs(emp_score - total_score) |
|
147
|
|
|
d3 = abs(total_score - pred_score) |
|
148
|
|
|
print "alpha=", alpha, "beta=", beta, "ks =", ks |
|
149
|
|
|
print pred_score, total_score, emp_score |
|
150
|
|
|
print "error=", d1 + d2 + d3 |
|
151
|
|
|
|
posterior /
distributions
