Mixed Posterior Sampling with Known and Unknown Conditional Posteriors

Modeling
1

I have a problem，parameters ‘theta’& ‘tau’ have known condition posterior close form(theta&tau have conjugate prior)，but parameter ‘nu’ condition posterior is unknown，only ‘nu’ need to use MCMC，how to implement code in R & Stan?

data {
  int<lower=0> N;
  int<lower=0> M;
  matrix [N,M] y;
  matrix [N,M] deltat;
  //hyperparameters
  real<lower=0> alpha_tau;
  real<lower=0> beta_tau;
  real mu_theta;
  real sigma_theta;
  real<lower=0> alpha_nu;
  real<lower=0> beta_nu;
}


parameters {
  //real theta;
  real eta;
  real<lower=0> nu;
  //real<lower=0> sigma;
  real<lower=0>tau;
  matrix<lower=0>[N,M] g;
}
transformed parameters{
  //prior
  real theta;
  theta = mu_theta+sigma_theta*eta;

}

model {
  // prior
  nu ~ gamma(alpha_nu, beta_nu);
  tau ~ gamma(alpha_tau, beta_tau);
  for(i in 1:N){
    for(j in 1:M){
      //print("before: i=",i,", j=",j,", gamma[i,j]=",g[i,j]);
      g[i,j] ~ gamma(deltat[i,j]*nu,1/nu);
      //print("after: i=",i,", j=",j,", gamma[i,j]=",g[i,j]);
      y[i,j] ~ normal(theta*g[i,j],sqrt(g[i,j]*1/tau));
      //y[i,j] ~ normal(theta*g[i,j],sigma*sqrt(g[i,j]));
    }
  }
}

library(rstan)
library(ggplot2)
# library(extraDistr)
library(MASS)
library(statmod)
library(progress)

stan_file <- "VGP_conjugate.stan"
stan_model <- stan_model(file = stan_file)
stan_model <- normalizePath(paste0(getwd(),"/VGP_conjugate.stan"))

#VGP sample code structure
VGP_generator_BGSS <- function(tau, theta, nu, n, T, mu=1){
  # n <- 2^k
  shape <- mu^2/nu
  scale <- nu/mu
  GP <- rgamma(n, shape=shape*T/n, scale=scale)
  Z <- rnorm(n, mean=0, sd=1)
  VGP <- rep(0, times=n+1)
  
  for (i in 1:n) {
    #VGP[i+1] <- VGP[i] + theta*GP[i] + sigma*sqrt(GP[i])*Z[i]
    VGP[i+1] <- VGP[i] + theta*GP[i] + sqrt(1/tau)*sqrt(GP[i])*Z[i]
  }
  
  VGP_df <- data.frame(
    degradation = VGP,
    time = seq(from=0,to=T,by=T/n)
  )
  degradation_increment <- diff(VGP_df$degradation)  
  time_increment <- diff(VGP_df$time)              
  
  increments_df <- data.frame(
    degradation_increment = degradation_increment,
    time_increment = time_increment
  )
  
  return(list(VGP_df = VGP_df, increments_df = increments_df))
}

#####Generate data (generate data do not need prior)
#set.seed(123)
n <- 1000
#parameters (true parameter, answer, later we will estimate these two parameter)
tau = 1
theta = 1
nu = 1
# VGP 
samples <- VGP_generator_BGSS(tau, theta, nu, n, T=1000)
y = samples$increments_df['degradation_increment'] 
deltat = samples$increments_df['time_increment']
y_matrix <- as.matrix(samples$increments_df['degradation_increment'])
deltat_matrix <- as.matrix(samples$increments_df['time_increment'])

####Bayesian setting
#hyperparameters
#mu_sigma = 1
#sigma_sigma = 1
alpha_tau = 1
beta_tau = 1
mu_theta = 1
sigma_theta = 1
alpha_nu = 1
beta_nu = 1

# fit model
stan_data <- list(N = n, M = 1, y = y_matrix, deltat = deltat_matrix, alpha_tau = alpha_tau, beta_tau = beta_tau, mu_theta = mu_theta, sigma_theta = sigma_theta, alpha_nu = alpha_nu, beta_nu = beta_nu)

stan_fit <- stan(file = stan_model, data = stan_data, iter = 2000, chains = 4, warmup = 500, thin = 10)
2

Hi, @wilson24.mg12 and welcome to the Stan community.

You can just code the model directly in Stan—you don’t have to worry about the fact that some of the variables have conjugate priors. Stan doesn’t use conjugacy information.

Now if you want to accelerate the model, you can exploit the conjugacy. There are some examples in the efficiency chapter of the User’s Guide. The idea is that you just code up the posterior rather than the prior and likelihood (or whatever else you’re using if this is in the middle of a stack of layers).

In a simpler case, consider a binomial with a conjugate beta prior.

data {
  int<lower=0> N;
  int<lower=0, upper=N> n;
}
parameters {
  real<lower=0, upper=1> theta;
}
model {
  theta ~ beta(3.2, 9.7);
  n ~ binomial(N, theta);
}

This is going to have the same sampling properties as just coding the posterior directly, because

beta(theta | 3.2, 9.7) * binomial(n | N, theta) 
  = constant * beta(theta | 3.2 + n, 9.7 + N).

We don’t need to worry about constants in Stan programs, so the Stan code for the unfolded prior is as follows.

data {
  int<lower=0> N;
  int<lower=0, upper=N> n;
}
parameters {
  real<lower=0, upper=1> theta;
}
model {
  theta ~ beta(3.2 + n, 9.7 + N - n);
}

All I’ve done is coded the posterior for theta rather than a prior and likelihood. As long as the target log density is equal to the posterior log density plus a constant, Stan is OK with unfolding conjugacy manually. You can do the same with normals.

2 Likes