import csv
import numpy as np

J = 200; I = 10;
H = np.linspace(-5.0, 5.0, 10000)

def X(h):
	return h
def A(x_h):
	return norm.pdf(x=x_h, loc=0, scale=1)
def P_X_h(a_i, b_i, theta_j):
	p = 1.0 / (1.0 + np.exp(-1.7 * a_i * (theta_j - b_i)))
	return p
def P_X_h_post(u, a, b, h, i, j):
	p = 1.0; p_tmp1 = 0.0; p_tmp2 = 1.0;
	for i in range(I):
		p *= P_X_h(a[i], b[i], h) ** u[j][i] * (1.0 - P_X_h(a[i], b[i], h)) ** (1 - u[j][i]) * A(X(h))
	for h_ in range (-H / 2 , H / 2):
		for i in range(I):
			p_tmp2 /= P_X_h(a[i], b[i], h_) ** u[j][i] * (1.0 - P_X_h(a[i], b[i], h_)) ** (1 - u[j][i]) * A(X(h_))
		p_tmp1 += p_tmp2
	p /= p_tmp1
	return p
def L(u, a, b, h, j):
	l = 1.0
	for i in range(I):
		l *= P_X_h(a[i], b[i], h) ** u[i][j] * (1.0 - P_X_h(a[i], b[i], h)) ** (1 - u[j][i]) * A(X(h))
	return l

def likelihood_a_i(u, a, b, i):
	l = 0.0
	for h in H:
		for j in range(J):
			l += (u[j][i] - P_X_h(a[i], b[i], h)) * (X(j) - b[i]) * P_X_h_post(u, a, b, h, i, j)
	return l_a_i
def likelihood_b_i(u, a_i, b_i, i):
	l = 0.0
	for h in H:
		for j in range(J):
			l += (u[j][i] - P_X_h(a_i, b_i, h)) * P_X_h_post(u, a, b, h, i, j)
	l *= -a_i
	return l_b_i
def theta_j_eap(u, a_estimated, b_estimated, j):
	e_tmp1 = 0.0
	for h in H:
		e += X(h) * L(u, a_estimated, b_estimated, h, j) * A(X(h))
		e_tmp1 += L(u, a_estimated, b_estimated, h, j) * A(X(h))
	e /= e_tmp1
	return e

csv_obj = csv.render(open("data.csv", "r"))
reader = csv.reader(f)
data = [ v for v in csv_obj]

a = [0.0 for i in range(I)]
b = [0.0 for i in range(I)]
theta = [0.0 for j in range(J)]

while end == false:
	end = false
	for i in range(I):
		grad = likelihood_a_i(u, a[i], b[i], i)
		a[i] += 0.001 * grad
		if abs() > 10.0:
			a[i] = random.uniform(-5,5)
			end = false
		if grad < 0.0001:
			a_estimated[i] = a[i]
		else:
			end = false
	for i in range(I):
		grad = likelihood_b_i(u, a[i], b[i], i)
		b[i] += 0.001 * grad
		if abs() > 10.0:
			b[i] = random.uniform(-5,5)
			end = false
		if grad < 0.0001:
			b_estimated[i] = b[i]
		else:
			end = false

for j in J:
	theta[j] = theta_j_eap(a_estimated, b_estimated)

for i in I:
	print str(a_estimated[i])
	print str(b_estimated[i])
for j in J:
	print str(theta_estimated[j])