-
Notifications
You must be signed in to change notification settings - Fork 0
/
sequential_importance_sampling_with_resample.py
133 lines (127 loc) · 3.27 KB
/
sequential_importance_sampling_with_resample.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
from scipy.stats import beta
from scipy.stats import norm
import matplotlib.mlab as mlab
import numpy as np
import matplotlib.pyplot as plt
import random
import bisect
from math import *
import collections
from matplotlib.pyplot import *
def cdf(weights):
total = sum(weights)
result = []
cumsum = 0
for w in weights:
cumsum += w
result.append(cumsum / total)
return result
def choice(population, weights):
assert len(population) == len(weights)
cdf_vals = cdf(weights)
x = random.random()
idx = bisect.bisect(cdf_vals, x)
return population[idx]
#weights=[0.3, 0.4, 0.3]
#population = 'ABC'
#counts = collections.defaultdict(int)
#for i in range(10000):
#counts[choice(population, weights)] += 1
#print(counts)
####sample from the proposal distribution at time t=1
N=100
resample_particles=[]
particles=np.random.normal(3,3,N)## proposal distribution
target=mlab.normpdf(particles,-3,1)## calculate the target pdf
proposal=mlab.normpdf(particles,3,3)## calculate the proposal pdf
particles=particles.tolist()
w_t=target/proposal ## calculate the weights of the particles
#newparticles=np.random.choice(particles,N,w_t)
for i in range(100):
particles_index=choice(particles, w_t)
#print particles_index
resample_particles.append(particles_index)
#print index
#print newparticles
#print w_t
## at time t=1:50 random walk and sample from new distributions
#print len(w_t)
w_t=np.ones(N)*0.01
#print w_t
a_update=[]
a=w_t
#print len(resample_particles)
#print suma
#suma=w_t
#print len(a)
for t in range(6):
particles_update=[]
w=[]
tar=[]
#a=[]
a_new=[]
for i in range(N):
#regenerate particles
#print resample_particles[i]
#print i
newparticles=np.random.normal(resample_particles[i],3,1)## generate new particles
particles_update.extend(newparticles)
proposal_update=mlab.normpdf(newparticles,resample_particles[i],3)
target_update=mlab.normpdf(newparticles,-3,1)
tar.extend(target_update)
#print particles_update
a[i]=a[i]*(target_update/(target[i]*proposal_update))
#print a_update
#a_update=w_t[i]*a
#print a[i]
a_new.append(a[i])
#### reform the varibles
#target=tar
#resample_particles=particles_update
resample_particles=[]
#print len(resample_particles)
for i in range(100):
particles_index=[]
particles_index=choice(particles_update, a_new)
#print particles_index
resample_particles.append(particles_index)
target=tar
#a=
#w_t=w
#print a_update[1]
numBins=500
#x=np.linspace(-5,5,100)
#plt.plot(x,)
#plt.show()
print resample_particles
#plot
#print len(suma)
#print len(w_t)
#w=w.__iter__()
#w=w.tolist()
#sumw=w[99]
#print w
#w=w/sum(w)*1.0
#print w_t
#w=w/sumw
#print w
#a=a/sum(a)
#print len(a)
#mean=sum(a*particles)
#x=np.linspace(-20,20,100)
#plt.plot(x,mlab.normpdf(x,-6,1),'r',label='target distribution')
#plt.plot(x,mlab.normpdf(x,1,2),'b',label='proposal distribution')
#plt.legend()
#vlines(particles,[0],a)
#plt.show()
#samples=np.random.choice(particles,N,w_t)
#plt.hist(samples)
#list=[1.1412, 4.3453, 5.8709, 0.1314]
#k=w.index(min(w))
#print k
#print particles[k]
#l=[]
#for i in range(N):
#l.append(w[i]*particles[i])
#s=sum(l)/N*1.0
#print s