# Copyright (c) 2005 Andrew Fraser
# This file is part of HMM_DS_Code.
#
# HMM_DS_Code is free software; you can redistribute it and/or modify it under
# the terms of the GNU General Public License as published by the Free Software
# Foundation; either version 2 of the License, or (at your option) any later
# version.
#
# HMM_DS_Code is distributed in the hope that it will be useful, but
# WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
# General Public License for more details.
#
# You should have received a copy of the GNU General Public License along with
# this program; if not, write to the Free Software Foundation, Inc.,
# 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.


#import psyco
#psyco.full()

import random
from math import *
random.seed(1)
mu = [-2.0,2.0]
Y = []
T = 10
for t in xrange(T):
   s = random.randint(0,1)       # Choose randomly 0 or 1
   y = random.gauss(mu[s],1.0)   # Draw from a Gaussian with mean
                                 # mu[s] and varaince 1
   print '%5.2f'%y,
   Y.append(y)
print '\n'
alpha = [0.5]
mu0 = [-1.0]
mu1 = [1.0]
for i in xrange(200):
    sum0 = 0.0
    sum1 = 0.0
    sums = 0.0
    ps = []
    for t in xrange(T):
        p0t = (alpha[i]/sqrt(2*pi)) * exp((-(Y[t]-mu0[i])**2)/2)
        p1t = ((1-alpha[i])/sqrt(2*pi)) * exp((-(Y[t]-mu1[i])**2)/2)
        ps.append(p0t/(p0t + p1t))
        sums += ps[t]
        print '%5.2f'%ps[t],
    print '\n',
    for t in xrange(T):
        sum0 += Y[t]*ps[t]
        print '%5.2f'%(Y[t]*ps[t]),
    print '\n',
    for t in xrange(T):
        sum1 += Y[t]*(1.0-ps[t])
        print '%5.2f'%(Y[t]*(1.0-ps[t])),
    print '\n',
    alpha.append(sums/T)
    mu0.append(sum0/sums)
    mu1.append(sum1/(T-sums))
    print i, alpha[-1], mu0[-1], mu1[-1], '\n'
#Local Variables:
#mode:python
#End: