# Nested MCMC on Stan

Hi Everyone!

I am trying to implement this Nested MCMC on RStan, i am not sure if it is possible!

1. for s = 1, . . . , N1,
Sample θ1(s) ~ p(y1|θ1) p(θ1)

2. For s = 1,…,N1,
Sample θ2(s) ~ p(θ2|θ1(s) ,Y2)

Return (θ1(s) * θ2(s))

These are the various likelihood distributions for θ1, θ2, Y1, Y2

Y1 = (Y1.1,Y1.2 . . . , Y1.n1).
Y1.i ∼ N (θ1, 1) for all 1 ≤ i ≤ n1.
A prior distribution for θ1 is N (0, λ1^-1).
Y1 = (Y2.1,Y2.2 . . . , Y2.n2). Y2.i~ N (θ1 + θ2, 1) for all 1 ≤ i ≤ n2.
The prior distribution on θ2 is N (0, λ2^-1).
The value of θ1* is set as 0, the value of θ2* is set as 1.
n1 is set as 100 and n2 = 1000,
Furthermore, we use λ1 = 1 and λ2 = 100.

I am able to do Part 1 on Stan, I am just not sure how to use the θ1 samples from Step 1 and Implement it to Step 2.

Sorry to inconvenience anyone!

Hi everyone, this is what i have attempted so far but it doesnt seem to work!

//model1.stan

data {
int<lower=0> N1;
vector[N1] Y1;
real<lower=0> lambda1;
}

parameters {
real theta1;
}

model {
for(i in 1:N1)
Y1[i] ~ normal(theta1, 1);
theta1 ~ normal(0,lambda1^-1);
}

library(Rcpp)
library(ggplot2)
library(rstan)
library(withr)
options(mc.cores = parallel::detectCores())
rstan_options(auto_write = TRUE)

N1 <-100
N2 <-1000
Y1 <- rnorm(N1, mean = 0, sd = 1)
Y2 <- rnorm(N2, mean = 1, sd = 1)
lambda1 <- 1
lambda2 <- 100

model_data1 <- list(N1=N1,Y1 = Y1,lambda1 = lambda1, lambda2 = lambda2)

fit2 <-stan(file = ‘model1.stan’, data = model_data1)

my_samples2 <- extract(fit2)

theta1samples <- my_samples2\$theta1[2001:2100]
Y2_2 <-rnorm(N1, mean = 1, sd = 1)
model_data4<-list(N1=N1,Y2 = Y2_2,theta1samples=theta1samples, lambda2 = lambda2)
fit4 <- stan(file = ‘model4.stan’, data = model_data4)
my_samples4 <- extract(fit4)

//model4.stan

data {
int<lower=0> N1;
vector[N1] Y2;
real<lower=0> lambda2;
vector[N1] theta1samples;
}

parameters {
real theta1;
real theta2;
}

model {
theta1 = theta1samples;
theta2 ~ normal(0, sqrt(lambda2^-1));
Y2 ~ normal(theta1 + theta2, 1);
}

It’s never (at least, I’ve never seen it) going to make sense to attempt to use one model’s output as the input to the next model. You need to think of a single generative model for all the data. And for your example, there is a very simple/straightforward one.

Is this a homework question for a course of some sort?

Hello Mike,

Thank you for your reply! This is for my thesis, my professor wanted me to work on this example and use Stan to code it up. I haven’t learnt Stan in school so i am learning it from scratch. Still working on it in the meanwhile

Ok, why don’t you take a look at the first few lectures here (Vids #8 through 12; feel free to skip 9). Your model will be a minor deviation from the two-groups example I cover pretty thoroughly by the end of 12.

Okay i will have a look at the lectures, thank you so much!