Nested MCMC on Stan

Algorithms
1

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!

2

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(StanHeaders)
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);
}

3

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?

4

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

5

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.

6

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