# 2012.11.08

library(DPpackage)

data(rolling)

y =rolling$y1

p.hat = mean(y)/9 ; p.hat

freq=table(y) ; freq

expected.if.1.binomial = 320 * dbinom(1:9,9,p.hat)

expected.if.1.binomial

plot(cumsum(expected.if.1.binomial), cumsum(freq) )

plot(1:9,expected.if.1.binomial,type="h")
points(1:9,freq,pch=19)

sum( (freq - expected.if.1.binomial)^2 /  expected.if.1.binomial )

pchisq(100,7,lower.tail=F)


# ML fit of 2 different probabilities

log.lik = function(par) {
  p1 = exp(par[1]) / (1+exp(par[1]));
  p2 = exp(par[2]) / (1+exp(par[2]));
  fitted.p = dbinom(1:9,9,p1) + 
             dbinom(1:9,9,p2)
  return( sum(freq * log(fitted.p)) )
}

fit= optim(c(0,1),log.lik,method="BFGS",hessian=T,
     control=list(fnscale=-1,parscale=c(0.1,0.1)))
fit     
fit$par

logit1 = fit$par[1] ;              p1.hat = exp(logit1) / (1+exp(logit1))
logit2 = fit$par[2] ; p2.hat = exp(logit2) / (1+exp(logit2))
round(c(p1.hat, p2.hat, (p1.hat+p2.hat)/2),2)

fitted.2.binomials  = 160 * dbinom(1:9,9,p1.hat) + 160 * dbinom(1:9,9,p2.hat)
points((1:9)+0.1, fitted.2.binomials, type="h",col="blue")

sum( (freq - fitted.2.binomials)^2 /  fitted.2.binomials )

##


#################################################### method of Moments

#    

var(y) # 3.463548

var.if.1.binomial = 9*p.hat*(1-p.hat) ; var.if.1.binomial # 2.050303

var(y) = V [E ] + E [V ]  

E [V ]  = (1/2)*9*p1*(1-p1) + (1/2)*9*p2*(1-p2)  = slightly less than 9*p.hat*(1-p.hat)


V[ E] = ( 1/2 distance b/w 9*p1 and 9*p2) ^ 2

est  of V[E] = 3.463548 - 2.050303 = 1.41

est of ( 1/2 distance b/w 9*p1 and 9*p2) =  sqrt(1.41) = 1.19

est of ( 1/2 distance b/w p1 and p2) =  1.19/9

est of distance b/w p1 and p2 =  2*1.19/9 = 0.26

ave of p1 and p2 =  0.65

so p1 and p2 estimated to be 0.52 and 0.68


#################################################### EM algorithm

UNDER CONSTRUCTION.. not correct 

fraction.from.1st.source = .4; p1=.45 ; p2=.85

for (iteration in 1:10) { 
	fitted.freq1 = 320 *  fraction.from.1st.source    * dbinom(1:9,9,p1)
	fitted.freq2 = 320 * (1-fraction.from.1st.source) * dbinom(1:9,9,p2)
  	p1 = sum((1:9)*fitted.freq1)/ sum(fitted.freq1)
  	p2 = sum((1:9),fitted.freq2)/ sum(fitted.freq2)
    fraction.from.1st.source = sum(fitted.freq1)/320
    print(c(fraction.from.1st.source,p1,p2))
  } 



#################################################### MCMC

n=length(y)

MixtureXModel=readLines(,11)
model
{
  for (i in 1:n) { 
  	z[i] ~ dbern(fraction.from.1st.source)
  	p[i]  <- p1*(1-z[i]) + p2*z[i] 
  	y[i]  ~ dbinom(p[i],9)
  } 
  p1 ~ dunif(0,1)
  p2 ~ dunif(0,1)
  fraction.from.1st.source ~ dunif(0,1)
}

writeLines(MixtureXModel,'MixtureXModel.txt')

INITS= list( p1 = .4, p2 = .9)


library('rjags')

jags <- jags.model("MixtureXModel.txt",
             data = list(n=n,y=y),
             inits=INITS,
             n.chains = 1,
             n.adapt = 4000)
 
update(jags, 10000)
 
results=jags.samples(jags, c('p1','p2','fraction.from.1st.source'), 10000)
str(results)

summary(results$p1[1,,1])
summary(results$p2[1,,1])
summary(results$fraction.from.1st.source[1,,1])