Вложение двух циклов для алгоритма ожидания-максимизации в r

Question

Привет, у меня есть код, который я хотел бы зациклить, но я не знаю как. Вы видите, я начинаю с некоторых предыдущих значений theta_A и theta_B, код запускается и выдает новые значения, которые в нижней части кода, который я назвал новая тэта_A и новая тэта_B. После первой итерации я хотел бы, чтобы новый theta_A заменил theta_A и аналогичный новый theta_B, чтобы заменить старый theta_B. Затем я хотел бы, чтобы код продолжал работать, пока значения не сойдутся в 1 десятичную точку.

У меня уже есть цикл for в моем коде, и я не совсем уверен, как объединить два цикла вместе и возможно ли это.

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

Я бы использовал функцию, тогда repeat

Проблема в том, что она не кажется мне сходящейся. Это идет к theta_A в 0.797 и theta_B в 0.52, и это повторяется с тем же результатом:

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
    }
}