Оценка максимального правдоподобия модели обучения с подкреплением

Question

Я пытаюсь смоделировать скорость обучения для задачи, которая похожа на задачу обучения по усилению бандитов. Единственная разница здесь заключается в том, что вознаграждение не является стохастическим (если сделан правильный выбор, агент награждается 1, в противном случае 0) и существует 6 возможных "подсказок".

По сути, агент должен узнать истинные значения для этих 6 сигналов. Возможные варианты: "Go" и "No-Go". Ниже показано, как я смоделировал задачу:

sim_PALP_M1 <- function(alpha,beta,N,palp){

 Q <- matrix(0.5, nrow=2, ncol=6)
 choice <- matrix(0, nrow=N, ncol=1)
 reward <- matrix(0, nrow=N, ncol=1)

 for(k in 1:N){

   # Get current stimulus of trial k
   currentstim <- palp[k,1]
   correct_ans <- palp[k,2]

   Pb <- 1 / (1 + exp(-beta*(Q[2,currentstim] - Q[1, currentstim]))) 


   # Now generate a random choice
   choice[k] <- sample(c(0,1),size=1, prob=c(1-Pb,Pb))

   reward[k] <- ifelse(choice[k] == correct_ans, 1, 0)

   # Update Q value of chosen option after seeing reward
   Q[choice[k]+1,currentstim] <- Q[choice[k]+1,currentstim] + 
   alpha*(reward[k] - Q[choice[k]+1,currentstim])
 }

 simdata <- list("c" = choice, "r" = reward, "Q" = Q)
}

Ниже приведена логарифмическая функция правдоподобия, используемая оптимизатором для максимизации:

LL_RLS <- function(param,choice,outcome, palp){


  N <- length(choice)
  Q <- matrix(0.5, nrow=2, ncol=6) # Initialise Q-values to 0.5 (arbitrary)
  LL <- 0

  for (k in 1:N) {

    currentstim <- palp[k,1]
    P2 <- 1 / (1 + exp(-param[2]*(Q[2,currentstim] - Q[1,currentstim])))

    if (choice[k] == 1) {
      LL <- LL + log(P2)
    }
    else{ LL <- LL + log(1-P2)}

    Q[choice[k]+1,currentstim] <- Q[choice[k]+1,currentstim] +  param[1]*(outcome[k] - Q[choice[k]+1,currentstim])

  }

  return(LL)
}

А вот как я использую Optim со случайными запусками и симулированными данными:

L=20
alpha = runif(L)
beta = runif(L,0,10)
N=dim(palpdat)[1]

# Matrix containing the simulated output
simres <- matrix(0, nrow=L, ncol=4)

for (k in 1:L) {

  simdat2 <- sim_PALP_M1(alpha[k],beta[k],0,0,N,palpdat)
  choice <- simdat2$c
  outcome <- simdat2$r
  M <- 150 # Number of starts for optim

  controlparams=list(fnscale=-1)

  limit <- length(choice) # Ignore this, was for other implementation
  results <- matrix(0, nrow=M, ncol=5)

  for (ii in 1:M) {
      alpha0 <- runif(1)
      beta0 <- runif(1,0.001,10)
      start_vals = c(alpha0,beta0)

      fit <- optim(par=start_vals, fn= LL_RLS, lower = c(0,0.001), upper = c(1,10), choice=choice[1:limit], outcome=outcome[1:limit],
                   palp=palpdat[1:limit,], method = "L-BFGS-B", control = controlparams)

      results[ii,1] <- fit$value
      results[ii,2:3] <- fit$par
      results[ii,4] <- alpha0
      results[ii,5] <- beta0
    }

 best <- results[results[,1] == min(results[,1]),]

 simres[k,1:2] <- c(beta[2],beta[3])
 simres[k,3:4] <- c(alpha[k], beta[k])

}

Данные задачи (palp) - это просто двумерная матрица с первым столбцом cues (цифры от 1 до 6) и вторым столбцом с правильным ответом (0 или 1). Чтобы создать поддельные данные, я просто передискретизировал строки матрицы столько раз, сколько необходимо. Я генерирую это просто с помощью:

stims <- c(1,2,3,4,5,6)
answers <- c(0,1,1,0,1,0)
v <- cbind(stims, answers)
palpdat <- create_PALP(v,200)

create_PALP <- function(v,n){
  # n is the number of trials with stims shuffled
  t <- matrix(v,ncol=2)
  palp <- t
  k = 0
  while (k<n) {
    palp <- rbind(palp, t[sample(nrow(t)),])
    k = k + 1
  }

  return(palp)

  }

Ниже приведен график зависимости моего фактического значения от фактического для обоих параметров (красный = установленный, зеленый = истинный). Я надеялся получить более точные оценки, так как это простая модель. Любая помощь приветствуется, спасибо большое! ( график альфа-оценок) ( график бета-оценок)