MatsuuraKentaro

data {
  int N;
  int M;
  real<lower=0> Y[N];
}

parameters {
  real<lower=0> mu;
  real<lower=0> phi;
  real<lower=1, upper=2> theta;

MatsuuraKentaro

#include <ostream>

static int iteration_index = 1;

inline void add_iter(std::ostream* pstream__) {
  iteration_index += 1;
}

inline int get_iter(std::ostream* pstream__) {
  return iteration_index;

MatsuuraKentaro

data {
  int I; // number of data
  int R; // number of Rcen
  int F; // number of Freq
  array[I] int Rcen; // value of Rcen
  array[I] int Freq; // value of Freq
  array[I] int N;
  array[I] int PV;
}

MatsuuraKentaro

functions {
  real f(int Y, real[] theta, real x) {
    return(gamma_lpdf(x | theta[2], theta[2]/theta[1]) + poisson_lpmf(Y | x));
  }

  real log_lik_Simpson(int Y, real[] theta, real a, real b, int M) {
    vector[M+1] lp;
    real h;
    h = (b-a)/M;
    lp[1] = f(Y, theta, a);

MatsuuraKentaro

library(cmdstanr)

waic <- function(log_likelihood) {
  training_error <- - mean(log(colMeans(exp(log_likelihood))))
  functional_variance_div_N <- mean(colMeans(log_likelihood^2) - colMeans(log_likelihood)^2)
  waic <- training_error + functional_variance_div_N
  return(waic)
}

model1_Si <- cmdstan_model('model1-addG-Simpson.stan')

MatsuuraKentaro

library(rstan)
library(doParallel)

replica.exchange.mcmc <- function (inv_T, n_ex, stanmodel, data, par_list, init, iter, warmup) {
  n_rep <- length(inv_T)
  len <- iter - warmup
  n_param <- sum(unlist(lapply(par_list, prod))) + 2  # number of parameters included E and lp__
  ms_T1 <- matrix(0, len*n_ex, n_param)               # MCMC samples at inv_T=1
  idx_tbl <- matrix(0, n_ex, n_rep)                   # index table of (exchange time, replica)
  E_tbl   <- matrix(0, n_ex, n_rep)                   # E table along idx_tbl

MatsuuraKentaro

functions {
  // calculating the convolutions
  vector convolution(vector x, vector y_rev) {
    int T = num_elements(x);
    vector[T-1] res;
    for (t in 2:T) {
      res[t-1] = dot_product(x[1:(t-1)], y_rev[(T-t+2):T]);
    }
    return res;
  }

MatsuuraKentaro

MatsuuraKentaro

MatsuuraKentaro

data {
  int A;
  int Y[A];
  int J;
  int From[J];
  int To[J];
}

parameters {
  simplex[A] q;