This is an old question, but when it popped up in the timeline I thought it would be a nice example of working with latent variables using Bayesian inference in stan:
library(rstan)
library(tidyverse)
a1 = 5
a2 = 10
b1 = 2
b2 = 3
e1 = .5
e2 = .7
n = 1000
x = rnorm(n)
y1 = a1 + b1 * x + rnorm(n, 0, e1)
y2 = a2 + b2 * x + rnorm(n, 0, e2)
data = data.frame(y1, y2)
code = '
data {
int n;
real y1[n];
real y2[n];
}
parameters {
real a1;
real a2;
real<lower=0> b1; // Make slopes positive
real<lower=0> b2;
real<lower=0> e1;
real<lower=0> e2;
real x[n];
}
model{
x ~ normal(0, 1); // Enforce a distribution on x
for(i in 1:n){
y1[i] ~ normal(a1 + x[i]*b1, e1);
y2[i] ~ normal(a2 + x[i]*b2, e2);
}
}
'
stan_data = list(
n = n, y1 = y1, y2 = y2
)
model = stan_model(model_code = code)
approx_model = vb(model, data = stan_data)
# Quick result using variational bayes
summary(approx_model, pars = c('a1', 'a2', 'b1', 'b2', 'e1', 'e2'))$summary %>%
round(digits = 2)
## mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff khat
## a1 5.00 NaN 0.02 4.96 4.99 5.00 5.02 5.04 NaN 3.98
## a2 9.99 NaN 0.03 9.93 9.97 9.99 10.01 10.06 NaN 3.99
## b1 1.82 NaN 0.02 1.78 1.81 1.82 1.83 1.85 NaN 3.99
## b2 2.72 NaN 0.03 2.67 2.71 2.73 2.74 2.78 NaN 3.99
## e1 0.51 NaN 0.01 0.49 0.50 0.51 0.52 0.53 NaN 3.99
## e2 0.77 NaN 0.02 0.73 0.76 0.77 0.78 0.80 NaN 3.98
# More precise, slower result using MCMC
# mcmc_model = sampling(model, chains = 2, cores = 2)
# summary(mcmc_model, pars = c('a1', 'a2', 'b1', 'b2', 'e1', 'e2'))$summary
# Estimated values of x
xhat = summary(approx_model)$summary %>%
data.frame() %>% rownames_to_column('parameter') %>%
filter(str_detect(parameter, 'x'))
plot(x, xhat$mean, xlab = 'True value', ylab = 'Estimated value')
