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')

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, data = stan_data, 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')