### Input the data I talked about in class.

data <- c(1.3,1.8,2.4,2.9,5.1,6.2,7.8,9.4,10.3)

### This is going to contain the "unobserved data"
zvec <- rep(.5,9)

### Some inityial parameter estimates. How important are starting values?
params <- c(1, 7, 0.4,2)


# A function to calculate a new Z.

newz <- function(y,mu1,mu2,theta,sigma){
return(theta * dnorm(y,mu1,sigma)/(
theta * dnorm(y,mu1,sigma) + (1 - theta) * dnorm(y,mu2,sigma))
)
}


### This is a crude looping function, that will iterate through the EM 
### Notice that I have not had to define the marginal likelihood to do this
### But I put it in as a check. This should decrease at each step as I
### defined it as -loglikelihood.

newvals <- function(){
for(i in 1:100){
zvec <- newz(data,params[1],params[2],params[3],params[4])
params[1] <- sum(zvec * data)/sum(zvec)
params[2] <- sum((1 - zvec) * data)/sum(1 - zvec)
params[3] <- mean(zvec)
print(params)
print(mixfunction(params,data))
}
return()
}



### This function calculates the mixture likelihood.

mixfunction <- function(params,thedata){
mu1 <- params[1]
mu2 <- params[2]
theta <- params[3]
y <- thedata
loglike <- -log(prod(theta * 1/sqrt(2 * pi * 2^2) * exp(-1/(2 * 2^2) *
(y - mu1)^2) +
(1 - theta) * 1/sqrt(2 * pi * 2^2) * exp(-1/(2 * 2^2) *
(y - mu2)^2)))
return(loglike)
}
 

### Try the non-linear minimizer, it should agree!

nlminb(start=c(2.7,8.3,.54),mixfunction,thedata=data,
lower=c(-Inf,-Inf,0),upper=c(Inf,Inf,1))

 
### What does the likelihood surface look like?

x <- y <- seq(-5,10,length=16)



calcgrid <- function(theta, sigma){
gridvals <- matrix(rep(0,16*16),ncol=16)
for(i in -5:10){
	for( j in -5:10){
	    params <- c(i,j,theta,sigma)
	gridvals[i+6,j+6] <- mixfunction(params,data)
}
}
            return(gridvals)
}

### motif()
### win.graph()
	   
par(mfrow=c(3,3))
for( k in 1:9)
{
persp(-calcgrid(.8,k/2))
}
	 

###### <------------------------>
###### What do these things look like?
###### What is the magic number?


meanmix <- function(x, mu1,mu2,sigma1,sigma2,theta) {
return(theta * dnorm(x,mu1,sigma1) + (1 - theta) * dnorm(x,mu2,sigma2))
}

par(mfrow=c(1,1))


plot(seq(-5,20,length=100),meanmix(seq(-5,20,length=100),2,6,.5,.5,.5),
type="n")
lines(seq(-5,20,length=100),meanmix(seq(-5,20,length=100),2,6,.5,.5,.5))


plot(seq(-5,20,length=100),meanmix(seq(-5,20,length=100),2,6,.75,.75,.5),
type="n")
lines(seq(-5,20,length=100),meanmix(seq(-5,20,length=100),2,6,.75,.75,.5))

plot(seq(-5,20,length=100),meanmix(seq(-5,20,length=100),2,6,1,1,.5),type="n")
lines(seq(-5,20,length=100),meanmix(seq(-5,20,length=100),2,6,1,1,.5))


plot(seq(-5,20,length=100),meanmix(seq(-5,20,length=100),2,6,1.5,1.5,.5),
type="n")
lines(seq(-5,20,length=100),meanmix(seq(-5,20,length=100),2,6,1.5,1.5,.5))


plot(seq(-5,20,length=100),meanmix(seq(-5,20,length=100),2,6,2,2,.5),
type="n")
lines(seq(-5,20,length=100),meanmix(seq(-5,20,length=100),2,6,2,2,.5))

 

plot(seq(-5,20,length=100),meanmix(seq(-5,20,length=100),2,6,2.5,2.5,.5),
type="n")
lines(seq(-5,20,length=100),meanmix(seq(-5,20,length=100),2,6,2.5,2.5,.5))