Navigation Menu

Skip to content

Instantly share code, notes, and snippets.

@MatsuuraKentaro
MatsuuraKentaro / model.stan
Last active March 5, 2024 06:08
Tweedie distribution in Stan
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
MatsuuraKentaro / model.stan
Last active April 19, 2017 01:08
cumulative_sum to reduce calculation
data {
int T;
vector[T] Y;
}
parameters {
real<lower=0, upper=2> mu_l;
real<lower=0, upper=2> mu_r;
real<lower=0> sigma;
}
@MatsuuraKentaro
MatsuuraKentaro / 0-replica-exchange-mcmc.R
Last active December 17, 2020 02:21
replica exchange MCMC
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
MatsuuraKentaro / model.stan
Last active October 21, 2017 02:29
correction of age heaping
data {
int A;
int Y[A];
int J;
int From[J];
int To[J];
}
parameters {
simplex[A] q;
@MatsuuraKentaro
MatsuuraKentaro / model.stan
Last active July 14, 2016 12:40
solving knapsack problem
functions {
real max_value(int I, int W, vector value, int[] weight) {
real dp[I+1,W+1];
for (w in 0:W) dp[1,w+1] = 0;
for (i in 1:I) {
for (w in 0:W) {
if (w < weight[i]) {
dp[i+1,w+1] = dp[i,w+1];
} else {
dp[i+1,w+1] = fmax(dp[i,w+1], dp[i,w-weight[i]+1] + value[i]);
@MatsuuraKentaro
MatsuuraKentaro / model.stan
Created July 14, 2016 12:46
calculation of expected time on random walk
functions {
real expected_time(int I, int J, int[,] can_goal, int[] di, int[] dj, vector p) {
matrix[I*J,I*J] A;
vector[I*J] b;
vector[I*J] res;
for (k1 in 1:(I*J)) {
b[k1] = 0;
for (k2 in 1:(I*J)) A[k1,k2] = 0;
}
functions {
vector slide_min(int I, int K, vector x) {
vector[I-K+1] res;
int deq[I];
int s;
int t;
s = 1;
t = 1;
for (i in 1:I) {
while (s < t && x[deq[t-1]] >= x[i])
@MatsuuraKentaro
MatsuuraKentaro / model-CMP.stan
Created August 12, 2016 02:29
underdispersed Poisson alternatives
functions {
real CMP_log_lik(int Y, real mu, real nu) {
return(nu * (Y * log(mu) - lgamma(Y+1)));
}
}
data {
int C;
int Y[C];
int Count[C];
@MatsuuraKentaro
MatsuuraKentaro / model1-add1-speedup.stan
Last active August 16, 2022 09:46
WAIC with hierarchical models
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
MatsuuraKentaro / model.stan
Created October 26, 2016 08:03
Bayesian GPLVM
data {
int N;
int K;
int D;
vector[N] Y[K];
}
transformed data {
vector[N] Mu;
Mu = rep_vector(0, N);