# Two different models have the same log likllihood

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\$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] ์ง์

``````

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!

1 Like