-
Notifications
You must be signed in to change notification settings - Fork 0
/
bayesian_inference_random_generate.py
45 lines (37 loc) · 1.3 KB
/
bayesian_inference_random_generate.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
import numpy as np
from scipy import io
from scipy.stats import wishart
import sys
if __name__ == '__main__':
if len(sys.argv) < 4:
print('usage: %s <inputfile> <N> <number of random precision>' % sys.argv[0])
sys.exit(1)
data = io.loadmat(sys.argv[1])
for key in data:
if key[0:2] != '__':
data = data[key]
break
N = int(sys.argv[2])
rand_p = int(sys.argv[3])
total, D = data.shape
true_cov = np.cov(data.T)
rand_data = data[np.random.choice(data.shape[0], N, replace=False), :]
w0 = np.identity(D)
v0 = 1
mean = np.array([1, -1])
post_v = v0+N
post_w = np.linalg.inv(w0) + np.dot((data[:N] - mean).T, data[:N] - mean)
post_w = np.linalg.inv(post_w)
test = wishart.rvs(df=post_v, scale=post_w, size=rand_p)
probability = np.empty(rand_p)
print '----testing----'
for i in range(rand_p):
probability[i] = wishart.pdf(test[i], df=post_v, scale=post_w)
pMAP = np.argmax(probability)
covMAP = np.linalg.inv(test[pMAP])
print 'the approximated MAP solution of covariance is'
print covMAP
print 'the true covariance is'
print true_cov
print 'error of approximated MAP solution is', np.mean(np.abs(covMAP - true_cov))
print 'the posterior probability is', probability[pMAP]