Instantly share code, notes, and snippets.

@WilsonMongwe WilsonMongwe/CoinTossing.R
Last active Aug 29, 2015

Embed
What would you like to do?
Coin Tossing
#### Author: Wilson Mongwe
#### Date: 21/06/2015
#### Website: www.wilsonmongwe.co.za
#### Title: Implementing the EM algorithm to the coin tossing example
####################### Simulating the population data #########################
N= 10000
numberOfCOins=2
components = sample(1:numberOfCOins,prob=c(0.45,0.55),size=N,replace=TRUE)
thetas= c(0.7,0.2) ### probablity of obtaining a head on Coin 2 and Coin 3 respectively
trials =10 ### number of times that each of Coin 2 and Coin 3 gets tossed.
x = rbinom(n=N,size=trials,prob=thetas[components]) ### sampling from a binomial distribution
plot(density(x),xlim=c(-2,12),ylim=c(0,0.26),xlab="Values",
ylab="Density",main="Approximate 'density' of binomial mixture",lwd=4)
par(new=TRUE)
########################### Setting up the EM algorithm #############################
#### Distribution of unboserved data given obseerved data #####
q_one=function(w)
{
lam = (w[1])
p_one=(w[2])
p_two=(w[3])
(dbinom(x,trials,p_one,log=FALSE)*lam)/
(dbinom(x,trials,p_one,log=FALSE)*lam+dbinom(x,trials,p_two,log=FALSE)*(1-lam))
}
q_two=function(w)
{
lam = (w[1])
p_one=(w[2])
p_two=(w[3])
(dbinom(x,trials,p_two,log=FALSE)*(1-lam))/
(dbinom(x,trials,p_one,log=FALSE)*lam+dbinom(x,trials,p_two,log=FALSE)*(1-lam))
}
#### Setting up the Likelihood ####
likelihood=function(w)
{
lam = (w[1])
p_one=(w[2])
p_two=(w[3])
result = log(dbinom(x,trials,p_one,log=FALSE)*lam/q_one(w))*q_one(w)+
log(dbinom(x,trials,p_two,log=FALSE)*(1-lam)/q_two(w))*q_two(w)
-sum(result)
}
##### Implementing the EM algorithm #####
plot(density(x),xlim=c(-2,12),ylim=c(0,0.26),xlab="Values",
ylab="Density",main="Approximate 'density' of binomial mixture",lwd=4)
par(new=TRUE)
start = c(runif(1,0.4,0.9),runif(1,0.4,0.9),runif(1,0.4,0.9))
tolerance = 20
i=1
functionValueCurrent=0
functionValueNew=100000
while(tolerance>10^-10)
{
if(i==1)
{
cat("Iteration number:\t", i,"Estimates:\t",c((start[1]),
(start[2]),
(start[3])), "\n")
}
functionValueCurrent=functionValueNew
p=optim(start,likelihood,method = "L-BFGS-B", lower=c(0.01,0.01,0.01),upper=c(0.99,0.99,0.99))
functionValueNew= p$value
tolerance= abs(functionValueNew-functionValueCurrent)
i=i+1
start=p$par
#cat("Current:\t", functionValueCurrent,"\n")
#cat("New:\t", functionValueNew,"\n")
#cat("Tolerance:\t", tolerance,"\n")
cat("Iteration number:\t", i,"Estimates:\t", c((start[1]),
(start[2]),
(start[3])), "\n")
components = sample(1:numberOfCOins,prob=c((start[1]),1-(start[1])),size=N,replace=TRUE)
thetas= c((start[2]),(start[3]))
samples = rbinom(n=N,size=trials,prob=thetas[components])
plot(density(samples),xlim=c(-2,12),ylim=c(0,0.26),xlab="",ylab="",main="",col="red",lwd=2)
par(new=TRUE)
}
plot(density(x),xlim=c(-2,12),ylim=c(0,0.26),xlab="Values",
ylab="Density",main="Approximate 'density' of binomial mixture",lwd=2)
################# GIF generation for website ###########################
ani.options(convert = 'C:\\Program Files\\ImageMagick-6.9.1-Q16\\convert.exe')
library(animation)
ani.options(interval=0.5)
start=c(runif(1,0.1,0.9),runif(1,0.1,0.9),runif(1,0.1,0.9))
saveGIF({
components = sample(1:numberOfCOins,prob=c((start[1]),1-(start[1])),size=N,replace=TRUE)
thetas= c((start[2]),(start[3]))
samples = rbinom(n=N,size=trials,prob=thetas[components])
plot(density(x),xlim=c(-2,12),ylim=c(0,0.4),xlab="0",
ylab="Density",main="Approximate 'density' of binomial mixture",col="black",lwd=4) ### plotting the density of the mixture
par(new=TRUE)
plot(density(samples),xlim=c(-2,12),ylim=c(0,0.4),xlab="",
ylab="",main="Approximate 'density' of binomial mixture",col="red",lwd=2)
legend('topright',c("Original", "EM Estimate"),lty=1, col=c('black', 'red'), bty='n', cex=1)
for( k in 1:10)
{
plot(density(x),xlim=c(-2,12),ylim=c(0,0.4),xlab=k,
ylab="Density",main="Approximate 'density' of binomial mixture",col="black",lwd=4) ### plotting the density of the mixture
par(new=TRUE)
p=optim(start,likelihood,method = "L-BFGS-B", lower=c(0.01,0.01,0.01),upper=c(0.99,0.99,0.99))
start=p$par
cat("Iteration number:\t", k,"Estimates:\t", c((start[1]),
(start[2]),
(start[3])), "\n")
components = sample(1:numberOfCOins,prob=c((start[1]),1-(start[1])),size=N,replace=TRUE)
thetas= c((start[2]),(start[3]))
samples = rbinom(n=N,size=trials,prob=thetas[components])
samples = rbinom(n=N,size=trials,prob=thetas[components])
plot(density(samples),xlim=c(-2,12),ylim=c(0,0.4),xlab=k,
ylab="",main="Approximate 'density' of binomial mixture",col="red",lwd=2)
legend('topright',c("Original", "EM Estimate"),lty=1, col=c('black', 'red'), bty='n', cex=1)
}
})
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment