from numpy import array, exp, ones, mean
from random import uniform, gauss

class Gaussian_Likelihood_Updater(object):
    def __init__(self, standard_deviation, number_hypothesis=200):
        self.h, self.probabilies, self.standard_deviation, self.n = array(range(number_hypothesis+1))/float(number_hypothesis*0.5)-1.0, ones(number_hypothesis+1), standard_deviation, number_hypothesis + 1
        self.Normalize()
    def Normalize(self):
        self.probabilies /= self.probabilies.sum()
    def Likelihood(self, hypothesis, evidence):
        return exp(-1.0*(hypothesis - evidence)**2/(2.0*self.standard_deviation**2))
    def Update(self, evidence):
        for i, hypothesis, probability in zip(xrange(self.n), self.h, self.probabilies):
            self.probabilies[i] *= self.Likelihood(hypothesis, evidence)
        self.Normalize()
    def Mean_Likelihood(self):
        mean_like = 0.0
        for hypothesis, probability in zip(self.h, self.probabilies):
            mean_like += hypothesis * probability
        return mean_like

def main():
    standard_deviation = 0.5
    iterations = 100000
    number_of_updates = 5
    mean_likelihood_lists = []
    for i in xrange(number_of_updates+1):
        mean_likelihood_lists.append([])
    for i in xrange(iterations):
        mean_this_iteration = uniform(-1,1)
        our_hypothesis = Gaussian_Likelihood_Updater(standard_deviation)
        for i in xrange(number_of_updates):
            mean_likelihood_lists[i].append(our_hypothesis.Mean_Likelihood())
            our_hypothesis.Update(gauss(mean_this_iteration,standard_deviation))
        mean_likelihood_lists[-1].append(our_hypothesis.Mean_Likelihood())
    for i, mean_likelihood_list in enumerate(mean_likelihood_lists):
        print "Mean value of estimates following " + str(i) + " updates:", mean(mean_likelihood_list)
        print "Mean absolute value of estimates following " + str(i) + " updates:", mean([abs(val) for val in mean_likelihood_list])
    print "\n", "-"*60, "\n"
    for i in xrange(number_of_updates):
        shift_list = [posterior_mean_likelihood-prior_mean_likelihood for posterior_mean_likelihood, prior_mean_likelihood in zip(mean_likelihood_lists[i+1], mean_likelihood_lists[i])]
        print "Mean Value of Shift for update " + str(i+1) + ":", mean(shift_list)
        print "Mean Absolute Value of Shift for update " + str(i+1) + ":", mean([abs(val) for val in shift_list])

if __name__ == '__main__':
    main()