Which one of the thetas you are interested in to plot? Anyway, you need to store the initial value in some variables and add the new values for each iteration. Then, add the plot instructions at the end of your function. Here is a proposal using theta1:
EMflipcoin = function(p,n,theta1,theta2,x) {
tol=10^(-10)
maxiter=100
diff=1
iter=0
value_theta = theta1 # save original values
value_log_theta = log(theta1)
while(diff>tol & iter<maxiter){
theta1 = theta1
theta2 = theta2
p=p
## E-step:
a = (theta1^x)*(1-theta1)^(n-x);a
b = (theta2^x)*(1-theta2)^(n-x);b
pa = a/(a+b);pa
pb = b/(a+b);pb
eah = pa*x;eah
eat = pa*(n-x);eat
ebh = pb*x;ebh
ebt = pb*(n-x);ebt
c = dbinom(x,n,p);c
c2 = dbinom(x,n,1-p);c2
pc = c/(c+c2);pc
## M-step:
theta1new = sum(eah)/(sum(eah)+sum(eat));theta1new
theta2new = sum(ebh)/(sum(ebh)+sum(ebt));theta2new
pnew = mean(pc);pnew
## calculate diff for tolerance in while condition
diff = (pnew-p)^2+(theta1new-theta1)^2+(theta2new-theta2)^2
iter =iter+1;
theta1 = theta1new
theta2 = theta2new
p =pnew
cat("Iter", iter,
", theta1=",theta1new,
", theta2=",theta2new,
", p =",pnew,
", diff=", diff, "\n")
value_theta = c(value_theta,theta1)
value_log_theta = c(value_log_theta, log(theta1) )
}
# plot theta vs log(theta)
df <- data.frame(value_log_theta, value_theta) #create a data frame with your variables
plot <- ggplot(df,aes(x = value_log_theta, y = value_theta)) +
geom_point() +
geom_line() +
labs(x = "log(theta)", y = "theta") +
theme_bw()
print(plot)
}
x=c(5,9,8,4,7)
EMflipcoin(p=0.5,n=10,theta1=0.7,theta2=0.4,x=x)