Plot logarithm(theta) vs theta

[homerjaysimpson]

I have created the below function for Em two binomial mixture and i want to plot the log(theta) versus theta of each iteration but i don't know how to do it?Any help i will appreciate it.

EMflipcoin = function(p,n,theta1,theta2,x) {
  
  tol=10^(-10)  
  maxiter=100
  diff=1
  iter=0
  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")}
  
}
x=c(5,9,8,4,7)
EMflipcoin(p=0.5,n=10,theta1=0.7,theta2=0.4,x=x)

[JossChavez]

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)