-
Notifications
You must be signed in to change notification settings - Fork 0
/
simpMixGaussNz.py
151 lines (125 loc) · 4.7 KB
/
simpMixGaussNz.py
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
import stats_util as s_u
import numpy as _N
import pickle
import matplotlib.pyplot as _plt
import scipy.cluster.vq as scv
import scipy.stats as _ss
import fitutil as _fu
import time as _tm
mvn = _N.random.multivariate_normal
class simpMixGauss:
################## hyper parameters - do not change during Gibbs
# HYPER PARAMS for prior covariance: nu, PSI
# how many clusters do I think there are
def fit(self, ITERS, M, x, _f_u=None, _f_q2=None, _q2_a=None, _q2_B=None, _ms_alp=None, f_0=None, q2_0=None, ms_0=None):
"""
Fit, with the inverting done in blocks
"""
oo = self
nSpks = x.shape[0]
xLo = _N.min(x)
xHi = _N.max(x)
f = _N.empty((ITERS, M, 1))
q2 = _N.empty((ITERS, M, 1))
ms = _N.empty((ITERS, M, 1))
xs = _N.sort(x)
if ms_0 is None:
rats = _N.ones(M)/M
rats += (1./M)*0.1*_N.random.randn(M)
rats /= _N.sum(rats)
ms[0, :, 0] = rats
else:
ms[0, :, 0] = ms_0
if f_0 is None:
ii = 0
for m in xrange(M):
f[0, m, 0] = _N.mean(xs[ii:ii+int(ms[0, m, 0]*nSpks)])
ii += int(ms[0, m, 0]*nSpks)
else:
f[0, :, 0] = f_0
if q2_0 is None:
ii = 0
for m in xrange(M):
q2[0, m, 0] = _N.std(xs[ii:ii+int(ms[0, m, 0]*nSpks)])**2
ii += int(ms[0, m, 0]*nSpks)
else:
q2[0, :, 0] = q2_0
xr = x.reshape((1, nSpks))
gz = _N.zeros((ITERS, nSpks, M), dtype=_N.int)
if _f_u is None:
_f_u = _N.empty(M)
ii = 0
for m in xrange(M):
_f_u[m] = _N.mean(xs[ii:ii+int(ms[0, m, 0]*nSpks)])
ii += int(ms[0, m, 0]*nSpks)
if _f_q2 is None:
_f_q2 = _N.ones(M)*(5**2)
if _q2_a is None:
_q2_a = _N.ones(M)*1e-4
if _q2_B is None:
_q2_B = _N.ones(M)*1e-3
if _ms_alp is None:
_ms_alp = _N.ones(M)*(0.001 / M)
# # termporary containers
econt = _N.empty((M, nSpks))
rat = _N.zeros((M+1, nSpks))
alp_ = _N.empty(M)
# print ms[0, :, 0]
# print f[0, :, 0]
# print q2[0, :, 0]
# print "^^^^^^^^"
oo.f = f
oo.q2 = q2
oo.ms = ms
oo.gz = gz
# initial values given for it == 0
for it in xrange(1, ITERS):
#################### STOCHASTIC ASSIGNMENT
norms = 1/_N.sqrt(2*_N.pi*q2[it-1])
zrs = _N.where(ms[it-1] == 0)[0]
ms[it-1, zrs, 0] = 1e-10
lms = _N.log(ms[it-1])
iq2 = 1./q2[it-1]
rnds = _N.random.rand(nSpks)
qdrSPC = (f[it-1] - xr)*(f[it-1] - xr)*iq2 # M x nSpks
cont = lms + norms - 0.5*qdrSPC
_N.max(cont, axis=0)
mcontr = _N.max(cont, axis=0).reshape((1, nSpks))
cont -= mcontr
_N.exp(cont, out=econt)
for m in xrange(M):
rat[m+1] = rat[m] + econt[m]
rat /= rat[M]
# print rat
M1 = rat[1:] >= rnds
M2 = rat[0:-1] <= rnds
gz[it] = (M1&M2).T
# prior for weights + likelihood used to sample
# _alp (prior) and alp_ posterior hyper
_N.add(_ms_alp, _N.sum(gz[it], axis=0), out=alp_)
############## SAMPLE WEIGHTS
ms[it, :, 0] = _N.random.dirichlet(alp_)
for m in xrange(M):
thisgr = _N.where(gz[it, :, m] == 1)[0]
nSpksC = len(thisgr)
if nSpksC > 0:
#### sample component MEANS
xbar = _N.mean(x[thisgr])
Mu = (_f_q2[m] * xbar + _f_u[m]*q2[it-1, m, 0]) / (_f_q2[m] + q2[it-1, m, 0])
S = (_f_q2[m] * q2[it-1, m, 0]) / (_f_q2[m] + q2[it-1, m, 0])
f[it, m, 0] = Mu + _N.sqrt(S)*_N.random.randn()
# print "----"
# print xbar
# print Mu
# print S
#### sample component variancs
a_ = 0.5*(nSpksC + 2*_q2_a[m])
B_ = 0.5*_N.sum((x[thisgr] - f[it, m, 0])*(x[thisgr] - f[it, m, 0])) + _q2_B[m]
q2[it, m, 0] = _ss.invgamma.rvs(a_, scale=B_)
else:
f[it, m, 0] = f[it-1, m, 0]
q2[it, m, 0] = q2[it-1, m, 0]
# print ms[it, :, 0]
# print f[it, :, 0]
# print q2[it, :, 0]
# print "^^^^^^^^"