RStanで混合分布の推定

混合分布のパラメタをRstanで推定する。

1. 正規分布が二つの場合

library(rstan)

src <- "
data {
  int N;
  int M;
  real x[N];
}

parameters {
  real mu[M];
  real sigma[M];
  simplex[M] pi;
}

model {
  real ps[M];
  for(n in 1:N){
    ps[1] <- log(pi[1]) + normal_log(x[n], mu[1], sigma[1]);
    ps[2] <- log(pi[2]) + normal_log(x[n], mu[2], sigma[2]);
    increment_log_prob(log_sum_exp(ps));
  }
}
"

x <- c(rnorm(150, 1, 1), rnorm(850, 10, 3))
dat <- list(N=length(x), M=2, x=x)

out <- stan(
  model_code = src,
  data = dat,
  iter = 4000,
  warmup = 2000,
  thin = 1,
  chains = 2,
)

summary(out)

2. ポアソン分布が二つの場合

library(rstan)

src <- "
  data {
    int N;
    int M;
    int x[N];
  }

  parameters {
    real kappa[M];
    simplex[M] pi;
  }

  model {
    real ps[M];
    for(n in 1:N){
      ps[1] <- log(pi[1]) + poisson_log(x[n], kappa[1]);
      ps[2] <- log(pi[2]) + poisson_log(x[n], kappa[2]);
      increment_log_prob(log_sum_exp(ps));
    }
  }
"

x <- c(rpois(150, 1), rpois(850, 5))
dat <- list(N=length(x), M=2, x=x)

out <- stan(
  model_code = src,
  data = dat,
  iter = 4000,
  warmup = 2000,
  thin = 1,
  chains = 2
)

summary(out)

3. 一様分布の上限を推定

library(rstan) 

src <- " 
  data { 
    int<lower=1> N; 
    real<lower=0> x[N]; 
  } 

  parameters { 
    real<lower=0, upper=1000> high; 
  } 

  model { 
    high ~ uniform(0,1000); 
    x ~ uniform(0, high); 
  } 
" 

x <- c(sample(0:10, 1000, replace=T)) 
dat <- list(N=1000, x=x) 

out <- stan( 
  model_code = src, 
  data = dat, 
  iter = 4000, 
  warmup = 2000, 
  thin = 1, 
  chains = 2, 
) 

summary(out)