Two different models have the same log likllihood

songb5454 · June 17, 2024, 6:34pm

Hello, I am undergraduate students major in psychology.

As a current class assignment, we are making a model of existing research into a heiarchical bayesian model.

The problem is that we created stan codes for two different models, and the log likelihood value of the two models is the same.

I will attach the code below, I am always using your program well, thank you.

The first code is below.

// model with alpha and beta
data {
  int<lower=1> N; //31,
  int<lower=1> T; //231
  int<lower=1> option_number;
  int<lower=1> context_number;
  int<lower=1> stimulus_number;
  int<lower=1, upper=2> choice[T, N];
  int<lower=1, upper=8> stimulus[T, option_number, N];
  int<lower=0, upper=4> context[T,N];
  real Rt[T, 2, N];
  real first_Q;
}

transformed data {
  
}

parameters {
  real<lower=0 ,upper = 1> mu_p_alpha;
  real<lower=0, upper = 1> mu_p_beta;
  real<lower=log(10), upper= log(1000)> x_for_sigma_alpha;
  real<lower=log(10), upper = log(1000)> x_for_sigma_beta;
    
  // Subject-level raw parameters (for Matt trick)
  vector<lower=0, upper=1>[N] alpha_pr;    // learning rate [0, 1]
  vector<lower=0, upper=1>[N] beta_pr;  // inverse temperature [0, 20]
}

transformed parameters {
  // subject-level parameters
  vector<lower=0,upper=1>[N] alpha;
  vector<lower=0,upper=20>[N] beta;
	real<lower=10, upper=1000> sigma_alpha = exp(x_for_sigma_alpha);
	real<lower=10, upper=1000> sigma_beta = exp(x_for_sigma_beta);
  
  for (i in 1:N) {/
    alpha[i] = alpha_pr[i];
    beta[i] = beta_pr[i]*20;
  }
}
model {
  // Hyperparameters,
  mu_p_alpha  ~ normal(0.25, 0.1);
  mu_p_beta  ~ normal(0.25, 0.1); 
  x_for_sigma_alpha ~ uniform(log(10),log(1000));  
  x_for_sigma_beta ~ normal(log(10),log(1000));
  
  
  alpha_pr ~ beta(mu_p_alpha * sigma_alpha ,(1-mu_p_alpha) * sigma_alpha);
  beta_pr ~ beta(mu_p_beta * sigma_beta ,(1-mu_p_beta) * sigma_beta);
  
  vector[stimulus_number] Q;
	vector[option_number] v;
  
  
  
  for (n in 1:N) {
    Q = rep_vector(first_Q,stimulus_number); //(27,27,27,27,27,27,27,27)
    
    for (t in 1:T) {        
      
      choice[t,n] ~ categorical_logit(beta[n]*Q[stimulus[t, 1:option_number, n]]); 
      
      if(context[t,n] !=0){
        
        for(k in 1:option_number) {		
					  
						v[k] = Rt[t, k, n]; 
        }	
					Q[stimulus[t, 1:option_number, n]] = Q[stimulus[t, 1:option_number, n]] + alpha[n] * (v - Q[stimulus[t, 1:option_number, n]]);

      }
    }
  }
}


generated quantities {

  real<lower=0, upper=20> mu_b = mu_p_beta * 20; 
  vector[N] log_lik;
  int<lower=0> fake_learn[context_number, N];
  int<lower=0> fake_test[stimulus_number, N];
  { 
   int what_pick;
   vector[stimulus_number] Q;
   vector[option_number] v;
	
 	 for (n in 1:N) {
	   log_lik[n] = 0;
	   fake_learn[1:context_number, n] = rep_array(0, context_number);
	   fake_test[1:stimulus_number, n] = rep_array(0, stimulus_number);
	  
	   Q = rep_vector(first_Q, stimulus_number);
	  
	   for (t in 1:T) {
	     log_lik[n] = categorical_logit_lpmf(choice[t, n] | beta[n] * Q[stimulus[t, 1:option_number, n]]);
	    
	     if (context[T,N] != 0){
	       what_pick = categorical_rng(softmax(beta[n] * Q[stimulus[t, 1:option_number, n]]));
	       fake_learn[context[t, n], n] = fake_learn[context[t, n], n] + (what_pick == option_number); 
	      
	       for (k in 1:option_number) {
	        
	         v[k] = Rt[t, k, n];
	        
	       }
	      
				 Q[stimulus[t, 1:option_number, n]] = Q[stimulus[t, 1:option_number, n]] + alpha[n] * (v - Q[stimulus[t, 1:option_number, n]]);
	      
	     }
	     else { 
	       what_pick = categorical_rng(softmax(beta[n] * Q[stimulus[t, 1:option_number, n]]));
	       fake_test[stimulus[t, what_pick, n], n] = fake_test[stimulus[t, what_pick, n], n] + 1;
	     }

	   }
	  }
	}
 
}

and below is the stan code for the second model.

// model with alpha and beta
data {
  int<lower=1> N; //31,
  int<lower=1> T; //231
  int<lower=1> option_number;
  int<lower=1> context_number;
  int<lower=1> stimulus_number;
  int<lower=1, upper=2> choice[T, N];
  int<lower=1, upper=8> stimulus[T, option_number, N];
  int<lower=0, upper=4> context[T,N];
  real Rt[T, 2, N];
  real first_Q_for_reference;
  real first_vs;
}

transformed data {
}

parameters {
  // Hyper(group)-parameters  
  real<lower=0, upper = 1> mu_p_alpha;
  real<lower=0, upper = 1> mu_p_beta;
  real<lower=0, upper = 1> mu_p_alpha_for_vs;
  real<lower=log(10), upper= log(1000)> x_for_sigma_alpha;
  real<lower=log(10), upper = log(1000)> x_for_sigma_beta;
  real<lower=log(10), upper = log(1000)> x_for_sigma_alpha_for_vs;  
  
  
    
  vector<lower=0, upper=1>[N] alpha_pr;    
  vector<lower=0, upper=1>[N] beta_pr;  // inverse temperature [0, 5]
  vector<lower=0, upper=1>[N] alpha_for_vs_pr; 
}

transformed parameters {
  vector<lower=0,upper=1>[N] alpha;
  vector<lower=0,upper=20>[N] beta;
  vector<lower=0,upper=1>[N] alpha_for_vs;
  
  for (i in 1:N) {
    alpha[i] = alpha_pr[i];
    beta[i] = beta_pr[i]* 20;
    alpha_for_vs[i] =alpha_for_vs_pr[i]; 
  }
  real<lower=10, upper=1000> sigma_alpha = exp(x_for_sigma_alpha);
  real<lower=10, upper=1000> sigma_beta = exp(x_for_sigma_beta);
  real<lower=10, upper=1000> sigma_alpha_for_vs = exp(x_for_sigma_alpha_for_vs);
}
model {
  mu_p_alpha  ~ normal(0.25,0.1);
  mu_p_beta  ~ normal(0.25,0.1);
  mu_p_alpha_for_vs ~ normal(0.25,0.1);
  
  alpha_pr ~ beta(mu_p_alpha * sigma_alpha ,(1-mu_p_alpha) * sigma_alpha);
  beta_pr ~ beta(mu_p_beta * sigma_beta ,(1-mu_p_beta) * sigma_beta);
  alpha_for_vs_pr ~ beta(mu_p_alpha_for_vs * sigma_alpha_for_vs ,(1-mu_p_alpha_for_vs) * sigma_alpha_for_vs);
  
  x_for_sigma_alpha ~ uniform(log(10),log(1000));  
  x_for_sigma_beta ~ normal(log(10),log(1000)); 
  x_for_sigma_alpha_for_vs ~ normal(log(10),log(1000)); 
  
  
  vector[stimulus_number] Q;
	vector[option_number] v;
	vector[context_number] vs;
  
  
  
  for (n in 1:N) {
    Q = rep_vector(first_Q_for_reference, stimulus_number);
    vs = rep_vector(first_vs,context_number);
    
    for (t in 1:T) {        
      
      choice[t,n] ~ categorical_logit( beta[n] * Q[stimulus[t, 1:option_number, n]] );
      
      if (context[t,n] != 0){
        
        for(k in 1:option_number) {		
          
					  vs[context[t,n]] = vs[context[t,n]] + alpha_for_vs[n] * ((Rt[t,1,n] + Rt[t,2,n])/2 - vs[context[t,n]]);

						v[k] = Rt[t, k, n]- vs[context[t,n]]; 
						
					
					
					}
			Q[stimulus[t, 1:option_number, n]] = Q[stimulus[t, 1:option_number, n]] + alpha[n] * (v - Q[stimulus[t, 1:option_number, n]]);
      }
    }
  }
}

generated quantities {
  real<lower=0, upper=20> mu_b = mu_p_beta * 20;
  vector[N] log_lik;
  int<lower=0> fake_learn[context_number, N];
  int<lower=0> fake_test[stimulus_number, N];
  
  
  {
   int what_pick;
   vector[stimulus_number] Q;
	 vector[option_number] v;
	 vector[context_number] vs;
	
	 for (n in 1:N) {
	   
	   log_lik[n] = 0;
	   fake_learn[1:context_number, n] = rep_array(0, context_number);
	   fake_test[1:stimulus_number, n] = rep_array(0, stimulus_number);
	  
	   Q = rep_vector(first_Q_for_reference, stimulus_number);
	   vs = rep_vector(first_vs,context_number);
	  
	   for (t in 1:T) {
	     
	     log_lik[n] = categorical_logit_lpmf(choice[t, n] | beta[n] * Q[stimulus[t, 1:option_number, n]]);
	    
	     if (context[T,N] != 0){
	       
	       // 1이면 가장 좋은 옵션 고른 거고 0이면 가장 안 좋은 옵션을 고른 거로 시뮬레이션함
	       what_pick = categorical_rng(softmax(beta[n] * Q[stimulus[t, 1:option_number, n]]));
	       fake_learn[context[t, n], n] = fake_learn[context[t, n], n] + (what_pick == option_number); 
	      
	       for (k in 1:option_number) {
	        
	         v[k] = Rt[t, k, n]- vs[context[t,n]]; 
	        
	       }
	      
				 Q[stimulus[t, 1:option_number, n]] = Q[stimulus[t, 1:option_number, n]] + alpha[n] * (v - Q[stimulus[t, 1:option_number, n]]);
	       vs[context[t,n]] = vs[context[t,n]] + alpha_for_vs_pr[n] * (mean(Rt[t, 1:option_umber, n]) - vs[context[t, n]]);  
	       
	     }
	     else { 
	       //transfer 상황 
	       what_pick = categorical_rng(softmax(beta[n] * Q[stimulus[t, 1:option_number, n]]));
	       fake_test[stimulus[t, what_pick, n], n] = fake_test[stimulus[t, what_pick, n], n] + 1;
	      
	       
	       }
	   }
	 }
  }
}

and this is the R code.

library(loo)  
library(rstan)
library(openxlsx)
library(dplyr)
setwd("C:/Users/손자연/OneDrive/바탕 화면/손경배 작업/인지모델링 기말 프로젝트/인지 기말프로젝트")
###########데이터 전처리####################
blocked_data = read.xlsx("blocked data.xlsx")

blocked_data$context[blocked_data$context == "NR"] <- 1
blocked_data$context[blocked_data$context == "NS"] <- 2
blocked_data$context[blocked_data$context == "PR"] <- 3
blocked_data$context[blocked_data$context == "PS"] <- 4

blocked_data$context = as.integer(blocked_data$context)

blocked_data
# 변수 설정
trial_number = 232  # 시행 횟수
subject_number = 31  # 참가자 수
stimulus_number = 8  # 자극 수
context_number = 4  # 컨텍스트 수 (baseline에는 필요 없음)
option_number = 2  # 옵션 수

# 배열 및 행렬 생성
Rt = array(dim = c(trial_number, option_number, subject_number))  # 실제 얻은 보상
stimulus = array(dim = c(trial_number, option_number, subject_number))  # 자극 배열 생성
context = matrix(nrow = trial_number, ncol = subject_number)  # 컨텍스트 행렬 생성
choice = matrix(nrow = trial_number, ncol = subject_number)  # 선택 행렬 생성

# 반복문을 사용하여 데이터를 피험자별로 처리
for (j in 1:subject_number) {
  # 컨텍스트 데이터를 해당 열에 할당
  context[, j] <- (blocked_data$context[blocked_data$PID == j])
  
  # 옵션 데이터를 해당 차원에 할당
  stimulus[, , j] <- cbind(blocked_data$option_1[blocked_data$PID == j], blocked_data$option_2[blocked_data$PID == j])
  
  # 보상 데이터를 해당 차원에 할당
  Rt[, , j] <- cbind(blocked_data$outcome_1[blocked_data$PID == j], blocked_data$outcome_2[blocked_data$PID == j])
  
  # 선택 데이터를 해당 열에 할당
  choice[, j] <- blocked_data$choice_recoded[blocked_data$PID == j]
}

choice
context[121:232,] = 0

dataList <- list(
  N =subject_number ,
  T=trial_number,
  stimulus_number=stimulus_number,
  context_number = context_number,
  option_number = option_number,
  stimulus = stimulus,
  context = context,
  Rt= Rt,
  choice  = choice,
  first_Q = 27,
  first_vs = 27,
  first_Q_for_reference = 0,
  first_rmin = 10, #global maximum
  first_rmax = 44 #global minimum
)
###########stan code 돌리기, baseline####################
output_baseline = stan("~/모델링프로젝트/baseline_stan코드.stan", data = dataList, pars = c("mu_p_alpha","mu_p_beta","sigma_alpha","sigma_beta","alpha", "beta","log_lik"),
                       iter = 1000, warmup=500, chains=4, cores=4)
parameters_baseline <- rstan::extract(output_baseline)
output_baseline
1+1
###########stan code 돌리기, reference####################
output_reference = stan("~/모델링프로젝트/reference_stan코드.stan", data = dataList, pars = c("mu_p_alpha","mu_p_beta","mu_p_alpha_for_vs","sigma_alpha","sigma_beta","sigma_alpha_for_vs","alpha", "beta","alpha_for_vs","log_lik"),
                        iter = 1000, warmup=500, chains=4, cores=4)
parameters_output_reference <- rstan::extract(output_reference)
output_reference

#############baseline LOOIC 계산##############
baseline_log_lik <- loo::extract_log_lik(stanfit = output_baseline,
                                         merge_chains = FALSE)

baseline_r_eff <- relative_eff(exp(baseline_log_lik), cores = 8)
baseline_looic <- loo(baseline_log_lik, r_eff = baseline_r_eff)#여기 있던 estimates[3, 1:2] 지움 

print(baseline_looic)

##############reference LOOIC 계산##############
reference_log_lik <- loo::extract_log_lik(stanfit = output_reference, 
                                          parameter_name = "log_lik",
                                          merge_chains = FALSE)
reference_log_lik
reference_r_eff <- relative_eff(exp(reference_log_lik), cores = 8)
reference_looic <- loo(reference_log_lik, r_eff = reference_r_eff)#여기 있던 estimates[3, 1:2] 지움

jsocolar · June 18, 2024, 12:22am

In the model block for the first model, this whole computation

for (t in 1:T) {        
      
      choice[t,n] ~ categorical_logit(beta[n]*Q[stimulus[t, 1:option_number, n]]); 
      
      if(context[t,n] !=0){
        
        for(k in 1:option_number) {		
					  
						v[k] = Rt[t, k, n]; 
        }	
					Q[stimulus[t, 1:option_number, n]] = Q[stimulus[t, 1:option_number, n]] + alpha[n] * (v - Q[stimulus[t, 1:option_number, n]]);

      }
    }

comes after all increments to the target, and therefore doesn’t affect lp__ at all. Something similar happens in the log likelihood compuation in generated quantities and in the second model.
Coding errors like these often show up when users are more familiar with JAGS, which infers a DAG from the model structure in such a way that the order of the statements doesn’t matter. In Stan, the order of the statements matters deeply. Where you have

log_lik[n] = categorical_logit_lpmf(choice[t, n] | beta[n] * Q[stimulus[t, 1:option_number, n]]);

you are incrementing the target using whatever values Q is populated with at this moment in the program, not using the values that you overwrite Q with when you later do

Q[stimulus[t, 1:option_number, n]] = Q[stimulus[t, 1:option_number, n]] + alpha[n] * (v - Q[stimulus[t, 1:option_number, n]]);

Hope that helps!

songb5454 · June 18, 2024, 1:41am

Really thank you for your kind answer!!!

Topic		Replies	Views
Newly Specified Bayesian Hierarchical Model with High Runtime Modeling techniques , specification	17	653	October 12, 2023
Write model for log-likelihood in Stan General	9	2145	December 11, 2020
Log-Likelihood with different parts Modeling	4	483	July 31, 2020
Bayesian multi-level (partial pooling) model -- multiple likelihood statements Modeling specification	2	526	July 1, 2022
Calculating the log likelihood for the Boxes model in Statistical Rethinking (2nd Ed) Modeling specification , rethinking	2	162	May 24, 2024

Two different models have the same log likllihood

Related Topics