混合二項分布の推定

 1: mix_bin <- function(alpha, N, p1, p2) {
 2:     r <- runif(1,0,1)
 3:     if(r < alpha) {
 4:         l <- 0
 5:         y <- rbinom(1,N,p1)
 6:     }else {
 7:         l <- 1
 8:         y <- rbinom(1,N,p2)
 9:     }
10:     return(c(y,l))
11: }
12: N <- 100
13: alpha <- 0.6
14: p1 <- 0.55
15: p2 <- 0.82
16: 
17: Num <- 100
18: ys <- c()
19: ls <- c()
20: for(i in 1:Num) {
21:     ret <- mix_bin(alpha,N,p1,p2)
22:     y <- ret[1]
23:     l <- ret[2]
24:     ys <- c(ys,y)
25:     ls <- c(ls,l)
26: }
27: #print(ys)
28: 
29: df <- data.frame(score=ys, label=ls)
30: head(df)
31: 
32: write.table(df,'mix_biom_test.csv',sep=',', quote=F,row.names=F,col.names=T)

1: df <- read.csv('mix_biom_test.csv', header=TRUE)
2: library(ggplot2)
3: 
4: g <- ggplot(data=df,mapping=aes(x=score,y=..density..)) + geom_density(alpha=0.4) + geom_histogram(alpha=0.2,color='black') + xlim(0,100)
5: fname <- 'mix_biom_test.png'
6: ggsave(file=fname, plot=g, dpi=70)

data {
  int N;
  int Y[N];
}

parameters {
  real<lower=0,upper=1> a;
  vector<lower=0,upper=1>[2] p;
}

model {
  for(n in 1:N) {
    target += log_sum_exp(log(a)   + binomial_lpmf(Y[n] | 100, p[1]),
                          log1m(a) + binomial_lpmf(Y[n] | 100, p[2]));
  }
}

1: data {
2:   int N;
3:   array[N] int Y;
4: }

data {
  int N;
  int Y[N];
}

parameters {
  real<lower=0,upper=1> a;
  # ordered<lower=0,upper=1>[2] p;  # エラーとなってしまう！！
  ordered[2] p_inv;
}

transformed parameters {
  real<lower=0, upper=1> p[2];
  p[1] <- inv_logit(p_inv[1]);
  p[2] <- inv_logit(p_inv[2]);
}

model {
  for(n in 1:N) {
    target += log_sum_exp(log(a)   + binomial_lpmf(Y[n] | 100, p[1]),
                          log1m(a) + binomial_lpmf(Y[n] | 100, p[2]));
  }
}

data {
  int N;
  int Y[N];
}

parameters {
  real<lower=0,upper=1> a;
  ordered[2] p_inv;
}

transformed parameters {
  real<lower=0, upper=1> p[2];
  p[1] = inv_logit(p_inv[1]);
  p[2] = inv_logit(p_inv[2]);

  vector[2] lp[N];
  for(n in 1:N) {
    lp[n,1] = log(a)   + binomial_lpmf(Y[n] | 100, p[1]);
    lp[n,2] = log1m(a) + binomial_lpmf(Y[n] | 100, p[2]);
  }
}

model {
  for(n in 1:N) {
    target += log_sum_exp(lp[n]);
  }
}

generated quantities {
  // データの予測分布を作る場合
  //vector[N] pi;
  //for(n in 1:N) {
  //  pi[n] = softmax(lp[n])[1];
  //}
  // 1 - 100 点までを使う場合
  vector[2] lx[100];
  for(n in 1:100) {
    lx[n,1] = log(a)   + binomial_lpmf(n | 100, p[1]);
    lx[n,2] = log1m(a) + binomial_lpmf(n | 100, p[2]);
  }
  vector[100] px;
  for(n in 1:100) {
    px[n] = softmax(lx[n])[1];
  }
}

1: library(cmdstanr)
2: df <- read.csv('mix_biom_test.csv', header=TRUE)
3: data <- list(N=nrow(df), Y=df$score)
4: model <- cmdstan_model('mix_biom_test.stan')
5: fit <- model$sample(data=data, chains=4, parallel_chains=4, seed=1234)
6: fit$save_output_files(dir="./", basename='mix_biom_test', timestamp=FALSE, random=FALSE)

 1: #df <- read.csv('mix_biom_test.csv', header=TRUE)
 2: #head(df)
 3: library(cmdstanr)
 4: files <-c('mix_biom_test-1.csv',
 5:           'mix_biom_test-2.csv',
 6:           'mix_biom_test-3.csv',
 7:           'mix_biom_test-4.csv')
 8: fit <- as_cmdstan_fit(files)
 9: fit$summary(c("a","p"))
10: fit$summary(c("px[1]","px[2]","px[3]","px[4]","px[5]","px[6]"))

1: library(cmdstanr)
2: files <-c('mix_biom_test-1.csv',
3:           'mix_biom_test-2.csv',
4:           'mix_biom_test-3.csv',
5:           'mix_biom_test-4.csv')
6: fit <- as_cmdstan_fit(files)
7: q95 <- function(x) quantile(x, probs=c(0.025, 0.25, 0.5, 0.75, 0.975))  # 自作関数
8: # 自作関数 q95 を含むいくつかの代表値を表示させる
9: fit$summary(c("a","p"), 'mean', 'sd', 'median', q95, 'rhat')

 1: library(rstan)
 2: files <-c('mix_biom_test-1.csv',
 3:           'mix_biom_test-2.csv',
 4:           'mix_biom_test-3.csv',
 5:           'mix_biom_test-4.csv')
 6: rstanfit <- rstan::read_stan_csv(files)  # rstan の fit object が作成できた. これを使う.
 7: options(width=180)
 8: summary(rstanfit)$summary[c('px[60]','px[90]'),]
 9: # 書き出す確率を変えたいとき
10: # summary(rstanfit, probs=c(0.1,0.5,0.9))$summary[c('px[60]','px[90]'),]
11: summary(rstanfit, pars=c('a','p'))$summary

 1: library(cmdstanr)
 2: library(ggplot2)
 3: files <-c('mix_biom_test-1.csv',
 4:           'mix_biom_test-2.csv',
 5:           'mix_biom_test-3.csv',
 6:           'mix_biom_test-4.csv')
 7: fit <- as_cmdstan_fit(files)
 8: p <- bayesplot::mcmc_dens(fit$draws(c('a','p')))
 9: # p <- bayesplot::mcmc_hist(fit$draws(c('a','p')))  # ヒストグラムを描く場合
10: fname <- 'mix_biom_test_dens.png'
11: ggsave(file=fname, plot=p, dpi=70)

 1: library(cmdstanr)
 2: library(ggplot2)
 3: files <-c('mix_biom_test-1.csv',
 4:           'mix_biom_test-2.csv',
 5:           'mix_biom_test-3.csv',
 6:           'mix_biom_test-4.csv')
 7: fit <- as_cmdstan_fit(files)
 8: p <- bayesplot::mcmc_dens_overlay(fit$draws(c('a','p')))
 9: fname <- 'mix_biom_test_dens_overlay.png'
10: ggsave(file=fname, plot=p, dpi=70)

 1: library(cmdstanr)
 2: library(ggplot2)
 3: files <-c('mix_biom_test-1.csv',
 4:           'mix_biom_test-2.csv',
 5:           'mix_biom_test-3.csv',
 6:           'mix_biom_test-4.csv')
 7: fit <- as_cmdstan_fit(files)
 8: p <- bayesplot::mcmc_pairs(fit$draws(c('a','p')))
 9: fname <- 'mix_biom_test_pairs.png'
10: ggsave(file=fname, plot=p, dpi=70)

 1: library(cmdstanr)
 2: library(ggplot2)
 3: files <-c('mix_biom_test-1.csv',
 4:           'mix_biom_test-2.csv',
 5:           'mix_biom_test-3.csv',
 6:           'mix_biom_test-4.csv')
 7: fit <- as_cmdstan_fit(files)
 8: p <- bayesplot::mcmc_trace(fit$draws(c('a','p')))
 9: fname <- 'mix_biom_test_trace.png'
10: ggsave(file=fname, plot=p, dpi=70)

 1: library(ggplot2)
 2: d <- read.csv('mix_biom_test.csv',header=TRUE)
 3: 
 4: p <- ggplot(data=d,mapping=aes(x=score,y=..density..))
 5: p <- p + geom_histogram(fill='white',color='black')
 6: p <- p + geom_density(fill='black',alpha=0.3)
 7: p <- p + xlim(0,100)
 8: 
 9: x <- seq(0,100, by=1)
10: alpha <- 0.55
11: p1 <- 0.55
12: y1 <- alpha * dbinom(x, 100, p1)
13: dy1 <- data.frame(x=x,y1=y1)
14: p <- p + geom_line(data=dy1,mapping=aes(x=x,y=y1),colour='green',size=2, linetype='dashed')
15: 
16: p2 <- 0.82
17: y2 <- (1 - alpha) * dbinom(x, 100, p2)
18: dy2 <- data.frame(x=x,y2=y2)
19: p <- p + geom_line(data=dy2,mapping=aes(x=x,y=y2),colour='blue',size=2, linetype='dashed')
20: 
21: fname <- 'mix_binom_test_hist.png'
22: ggsave(plot=p,file=fname,dpi=70)

 1: library(cmdstanr)
 2: files <-c('mix_biom_test-1.csv',
 3:           'mix_biom_test-2.csv',
 4:           'mix_biom_test-3.csv',
 5:           'mix_biom_test-4.csv')
 6: fit <- as_cmdstan_fit(files)
 7: mat <- fit$summary(c("px"), "median","rhat")
 8: mat$score <- 1:100
 9: 
10: library(ggplot2)
11: p <- ggplot(data=mat) + geom_point(mapping=aes(x=score, y=median)) + xlim(0,100)
12: fname <- 'mix_biom_test_pred.png'
13: ggsave(file=fname, plot=p, dpi=70)