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