Nesting two loops for Expectation-Maximisation algorithm in r

Question

Hi I have a code that I would like to loop but I don't know how. You see I start with some prior values theta_A and theta_B, the code runs and spits new values which at the bottom of the code I named new theta_A and new theta_B. After the first iteration i would like the new theta_A to replace theta_A and similar new theta_B to replace old theta_B. Then I would like the code to keep running until the values converge to 1 decimal point.

I already have a for loop in my code and I am not really sure how to nest two loops together and if it is even possible.

n=10 #patients undergoing treatment
R=c(9,8,5,7,4)  #patients recovered in each hospital
NR=c(1,2,5,3,6) #patients not recovered in each hospital

theta_A=0.6 #prior probability of recovering using treatment A
theta_B=0.5 #prior probability of recovering using treatment B

#We now need to create empty vectors for our probability functions

f_A=matrix(0, 1, 5)
f_B=matrix(0, 1, 5)

#We need to create empty vectors for our normalisation

fsum=matrix(0, 1, 5)
A=matrix(0, 1, 5)
B=matrix(0, 1, 5)

#We need to create empty vectors for expected numeber of recoveries using each treatment

AR=matrix(0, 1, 5)
ANR=matrix(0, 1, 5)
BR=matrix(0, 1, 5)
BNR=matrix(0, 1, 5)

# Now we employ the expectation algorithm, the for loops run through all 5 hospitals

for (i in 1:5) {
  f_A[i]<-(theta_A)^(R[i])*(1-theta_A)^(n-R[i]) #Chances of exact number of recoveries using treatment A
  f_B[i]<-(theta_B)^(R[i])*(1-theta_B)^(n-R[i]) #Chances of exact numebr of recoveries using treatment B
  #Normalisation
  fsum[i]<-f_A[i]+f_B[i] # Sum of the two recoveries
  #print(totalf[i])
  A[i]<-f_A[i]/fsum[i] #Chances of using treatment A
  B[i]<-f_B[i]/fsum[i] #Chances of using treatment B
  AR[i]<-R[i]*A[i]     #Expected recovered patients using treatment A
  ANR[i]<-NR[i]*A[i]   #Expected non-recovered patients using treatment A
  BR[i]<-R[i]*B[i]     #Expected recovered patients using treatment B
  BNR[i]<-NR[i]*B[i]   #Expected non-recovered patients using treatment B


}

# Now employ maximaisation algorithm
total_recA=sum(AR)
total_nonrecA=sum(ANR)
total_recB=sum(BR)
total_nonrecB=sum(BNR)

# Posterior probability of recovery
new_theta_A=total_recA/(total_recA+total_nonrecA)
new_theta_B=total_recB/(total_recB+total_nonrecB)

Answer 1

I'd use a function then a repeat

The problem is that it doesn't appear to be converging for me. It goes to theta_A at 0.797 and theta_B at 0.52 and that repeats with the same result:

n = 10 #patients undergoing treatment
R = c(9, 8, 5, 7, 4)  #patients recovered in each hospital
NR = c(1, 2, 5, 3, 6) #patients not recovered in each hospital

theta_A = 0.6 #prior probability of recovering using treatment A
theta_B = 0.5 #prior probability of recovering using treatment B

thetaFunction <- function(theta_A, theta_B) {
    #We now need to create empty vectors for our probability functions

    f_A = matrix(0, 1, 5)
    f_B = matrix(0, 1, 5)

    #We need to create empty vectors for our normalisation

    fsum = matrix(0, 1, 5)
    A = matrix(0, 1, 5)
    B = matrix(0, 1, 5)

    #We need to create empty vectors for expected numeber of recoveries     using each treatment

    AR = matrix(0, 1, 5)
    ANR = matrix(0, 1, 5)
    BR = matrix(0, 1, 5)
    BNR = matrix(0, 1, 5)

    # Now we employ the expectation algorithm, the for loops run through all     5 hospitals

    for (i in 1:5) {
        f_A[i] <-
            (theta_A) ^ (R[i]) * (1 - theta_A) ^ (n - R[i]) #Chances of     exact number of recoveries using treatment A
        f_B[i] <-
            (theta_B) ^ (R[i]) * (1 - theta_B) ^ (n - R[i]) #Chances of     exact numebr of recoveries using treatment B
        #Normalisation
        fsum[i] <- f_A[i] + f_B[i] # Sum of the two recoveries
        #print(totalf[i])
        A[i] <- f_A[i] / fsum[i] #Chances of using treatment A
        B[i] <- f_B[i] / fsum[i] #Chances of using treatment B
        AR[i] <-
            R[i] * A[i]     #Expected recovered patients using treatment A
    ANR[i] <-
        NR[i] * A[i]   #Expected non-recovered patients using treatment A
        BR[i] <-
            R[i] * B[i]     #Expected recovered patients using treatment B
        BNR[i] <-
            NR[i] * B[i]   #Expected non-recovered patients using treatment B


    }

    # Now employ maximaisation algorithm
    total_recA = sum(AR)
    total_nonrecA = sum(ANR)
    total_recB = sum(BR)
    total_nonrecB = sum(BNR)

    # Posterior probability of recovery
    new_theta_A = total_recA / (total_recA + total_nonrecA)
    new_theta_B = total_recB / (total_recB + total_nonrecB)

    c(new_theta_A, new_theta_B)
}

result <- c(theta_A, theta_B)

repeat{
    priorResult <- result
    result <- thetaFunction(result[1], result[2])
    if (abs(result[1] - priorResult [1]) <=0.01 & abs(result[2] - priorResult [2]) <=0.01){
        break
    }
}