Stan と R でベイズ統計モデリング 第四章

 1: data {
 2:   int N;
 3:   real Y[N];
 4: }
 5: 
 6: parameters {
 7:   real mu;
 8: }
 9: 
10: model {
11:   for(i in 1:N) {
12:     Y[N] ~ normal(mu,1)
13:   }
14:   mu ~ normal(0,100);
15: }

1: library(ggplot2)
2: d <- read.csv(file='data-salary.txt')
3: p <- ggplot(data=d) + geom_point(mapping=aes(x=X,y=Y),shape=1, size=3)
4: fname <- '04_01_fig.png'
5: ggsave(file=fname, plot=p, dpi=50)   # dpi で図の大きさ(解像度)を指定

1: d <- read.csv(file='data-salary.txt')
2: res <- lm(Y ~ X, data=d)
3: summary(res)

 1: d <- read.csv(file='data-salary.txt')
 2: res <- lm(Y ~ X, data=d)
 3: 
 4: X_new <- data.frame(X=23:60)
 5: conf_95 <- predict(res, X_new, interval='confidence', level=0.95)
 6: conf_50 <- predict(res, X_new, interval='confidence', level=0.50)
 7: conf_95 <- data.frame(X_new,conf_95)  # 列名 = (X,fit,lwr,upr)
 8: conf_50 <- data.frame(X_new,conf_50)
 9: 
10: library(ggplot2)
11: p <- ggplot()
12: p <- p + geom_ribbon(data=conf_95, mapping=aes(x=X,ymin=lwr,ymax=upr), alpha=1/6)
13: p <- p + geom_ribbon(data=conf_50, mapping=aes(x=X,ymin=lwr,ymax=upr), alpha=3/6)
14: p <- p + geom_line(data=conf_50,aes(x=X,y=fit),size=1) # 線
15: p <- p + geom_point(data=d, aes(x=X, y=Y), shape=1, size=3) # まる
16: p <- p + labs(x='X', y='Y') + coord_cartesian(xlim=c(22, 61), ylim=c(200, 1400))
17: p <- p + scale_y_continuous(breaks=seq(from=200, to=1400, by=400))
18: fname <- 'fig4-3-left.png'
19: ggsave(file=fname, plot=p, dpi=50)

 1: d <- read.csv(file='data-salary.txt')
 2: res <- lm(Y ~ X, data=d)
 3: 
 4: X_new <- data.frame(X=23:60)
 5: pred_95 <- predict(res, X_new, interval='prediction', level=0.95)
 6: pred_50 <- predict(res, X_new, interval='prediction', level=0.50)
 7: pred_95 <- data.frame(X_new,pred_95)  # 列名 = (X,fit,lwr,upr)
 8: pred_50 <- data.frame(X_new,pred_50)
 9: 
10: library(ggplot2)
11: p <- ggplot()
12: p <- p + geom_ribbon(data=pred_95, mapping=aes(x=X,ymin=lwr,ymax=upr), alpha=1/6)
13: p <- p + geom_ribbon(data=pred_50, mapping=aes(x=X,ymin=lwr,ymax=upr), alpha=3/6)
14: p <- p + geom_line(data=pred_50,aes(x=X,y=fit),size=1) # 線
15: p <- p + geom_point(data=d, aes(x=X, y=Y), shape=1, size=3) # まる
16: p <- p + labs(x='X', y='Y') + coord_cartesian(xlim=c(22, 61), ylim=c(200, 1400))
17: p <- p + scale_y_continuous(breaks=seq(from=200, to=1400, by=400))
18: fname <- 'fig4-3-right.png'
19: ggsave(file=fname, plot=p, dpi=50)

 1: data {
 2:   int N;
 3:   real X[N];
 4:   real Y[N];
 5: }
 6: 
 7: parameters {
 8:   real a;
 9:   real b;
10:   real<lower=0> sigma;
11: }
12: 
13: model {
14:   for(i in 1:N) {
15:     Y[i] ~ normal(a + b*X[i], sigma);
16:   }
17: }

1: library(cmdstanr)
2: d <- read.csv(file='data-salary.txt')
3: data <- list(N=nrow(d), X=d$X, Y=d$Y)
4: model  <- cmdstan_model('model4-5.stan')
5: fit <- model$sample(data=data, chains=4, parallel_chains=4, seed=1234)
6: fit$save_output_files(dir="./", basename='model4-5', timestamp=FALSE, random=FALSE)

1: library(cmdstanr)
2: files <- c('model4-5-1.csv',
3:            'model4-5-2.csv',
4:            'model4-5-3.csv',
5:            'model4-5-4.csv')
6: fit <- as_cmdstan_fit(files)
7: fit$summary()

1: library(rstan)
2: files <-c('model4-5-1.csv',
3:           'model4-5-2.csv',
4:           'model4-5-3.csv',
5:           'model4-5-4.csv')
6: rstanfit <- rstan::read_stan_csv(files)  # rstan の fit object が作成できた. これを使う.
7: options(width=180)
8: summary(rstanfit)$summary

1: library(cmdstanr)
2: files <- c('model4-5-1.csv',
3:            'model4-5-2.csv',
4:            'model4-5-3.csv',
5:            'model4-5-4.csv')
6: fit <- as_cmdstan_fit(files)
7: ms <- fit$draws(format='matrix')  # chain を一つにまとめる
8: dim(ms)
9: head(ms)

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

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

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

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

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

 1: library(cmdstanr)
 2: files <- c('model4-5-1.csv',
 3:            'model4-5-2.csv',
 4:            'model4-5-3.csv',
 5:            'model4-5-4.csv')
 6: fit <- as_cmdstan_fit(files)
 7: ms <- fit$draws(format='df')
 8: head(ms)
 9: head(ms$a)
10: length(ms$a)
11: quantile(ms$b, probs=c(0.025,0.975))

 1: library(ggplot2)
 2: library(cmdstanr)
 3: files <- c('model4-5-1.csv',
 4:            'model4-5-2.csv',
 5:            'model4-5-3.csv',
 6:            'model4-5-4.csv')
 7: fit <- as_cmdstan_fit(files)
 8: df_mcmc <- fit$draws(format='df')
 9: # 散布図
10: p_xy <- ggplot(df_mcmc,aes(x=a,y=b)) +
11:     geom_point(alpha=1/4,size=2,shape=1) +
12:     scale_x_continuous(breaks=seq(-400,200,200), limits=c(-420,210)) +
13:     scale_y_continuous(breaks=seq(15,25,5),limits=c(14.5,29))
14: 
15: # a の分布図
16: p_x <- ggplot(df_mcmc,aes(x=a)) + theme_bw(base_size=18) + theme(legend.position='none') +
17:     geom_histogram(aes(y=..density..),color='black',fill='white') +
18:     geom_density(alpha=0.3,fill='gray20') +
19:     scale_x_continuous(breaks=seq(-400,200,200),limits=c(-420,210))+
20:     labs(x='',y='')
21: # b の分布図
22: p_y <- ggplot(df_mcmc, aes(x=b)) + theme_bw(base_size=18) + theme(legend.position='none') +
23:     coord_flip() +
24:     geom_histogram(aes(y=..density..), colour='black', fill='white') +
25:     geom_density(alpha=0.3, fill='gray20') +
26:     scale_x_continuous(breaks=seq(15, 25, 5), limits=c(14.5, 29)) +
27:     labs(x='', y='')
28: 
29: p_emp <- ggplot(data.frame(0,0)) + theme(rect=element_rect(fill='white'), panel.border=element_blank())
30: 
31: g_xy  <- ggplotGrob(p_xy)
32: g_x   <- ggplotGrob(p_x)
33: g_y   <- ggplotGrob(p_y)
34: g_emp <- ggplotGrob(p_emp)
35: 
36: g1 <- cbind(g_x,  g_emp, size='first')
37: g2 <- cbind(g_xy, g_y,   size='first')
38: g  <- rbind(g1,   g2,    size='first')
39: fname <- 'fig4-7.png'
40: ggsave(file=fname, plot=g, dpi=70)

 1: library(ggplot2)
 2: library(cmdstanr)
 3: files <- c('model4-5-1.csv',
 4:            'model4-5-2.csv',
 5:            'model4-5-3.csv',
 6:            'model4-5-4.csv')
 7: fit <- as_cmdstan_fit(files)
 8: ms <- fit$draws(format='df')
 9: 
10: X_new <- 23:60
11: N_X <- length(X_new)
12: N_mcmc <- length(ms$lp__)
13: 
14: set.seed(1234)
15: y_base_mcmc <- as.data.frame(matrix(nrow=N_mcmc,ncol=N_X))
16: y_mcmc      <- as.data.frame(matrix(nrow=N_mcmc,ncol=N_X))
17: 
18: # i 列目のデータを作る
19: for (i in 1:N_X) {
20:     y_base_mcmc[,i] <- ms$a + ms$b * X_new[i]
21:     y_mcmc[,i]      <- rnorm(n=N_mcmc,mean=y_base_mcmc[,i],sd=ms$sigma)
22: }
23: conf <- data.frame(X=X_new,t(apply(y_base_mcmc,2,quantile,probs=c(0.025,0.25,0.50,0.75,0.975))))
24: colnames(conf) <- c('X','p025','p250','p500','p750','p975')
25: 
26: d <- read.csv(file='data-salary.txt')
27: p <- ggplot()
28: p <- p + geom_ribbon(data=conf,aes(x=X,ymin=p025,ymax=p975),alpha=1/6)
29: p <- p + geom_ribbon(data=conf,aes(x=X,ymin=p250,ymax=p750),alpha=2/6)
30: p <- p + geom_line(  data=conf,aes(x=X,y=p500),size=0.2) # 線
31: p <- p + geom_point( data=d,   aes(x=X,y=Y), shape=1, size=3) # まる
32: p <- p + labs(x='X', y='Y') + coord_cartesian(xlim=c(22, 61), ylim=c(200, 1400))
33: p <- p + scale_y_continuous(breaks=seq(from=200, to=1400, by=400))
34: 
35: fname <- 'fig4-8-left.png'
36: ggsave(file=fname, plot=p, dpi=70)

 1: library(ggplot2)
 2: library(cmdstanr)
 3: files <- c('model4-5-1.csv',
 4:            'model4-5-2.csv',
 5:            'model4-5-3.csv',
 6:            'model4-5-4.csv')
 7: fit <- as_cmdstan_fit(files)
 8: ms <- fit$draws(format='df')
 9: 
10: X_new <- 23:60
11: N_X <- length(X_new)
12: N_mcmc <- length(ms$lp__)
13: 
14: set.seed(1234)
15: y_base_mcmc <- as.data.frame(matrix(nrow=N_mcmc,ncol=N_X))
16: y_mcmc      <- as.data.frame(matrix(nrow=N_mcmc,ncol=N_X))
17: 
18: # i 列目のデータを作る
19: for (i in 1:N_X) {
20:     y_base_mcmc[,i] <- ms$a + ms$b * X_new[i]
21:     y_mcmc[,i]      <- rnorm(n=N_mcmc,mean=y_base_mcmc[,i],sd=ms$sigma)
22: }
23: 
24: # y_base_mcmc --> y_mcmc に変更しただけ
25: conf <- data.frame(X=X_new,t(apply(y_mcmc,2,quantile,probs=c(0.025,0.25,0.50,0.75,0.975))))
26: colnames(conf) <- c('X','p025','p250','p500','p750','p975')
27: 
28: d <- read.csv(file='data-salary.txt')
29: p <- ggplot()
30: p <- p + geom_ribbon(data=conf,aes(x=X,ymin=p025,ymax=p975),alpha=1/6)
31: p <- p + geom_ribbon(data=conf,aes(x=X,ymin=p250,ymax=p750),alpha=2/6)
32: p <- p + geom_line(  data=conf,aes(x=X,y=p500),size=0.2) # 線
33: p <- p + geom_point( data=d,   aes(x=X,y=Y), shape=1, size=3) # まる
34: p <- p + labs(x='X', y='Y') + coord_cartesian(xlim=c(22, 61), ylim=c(200, 1400))
35: p <- p + scale_y_continuous(breaks=seq(from=200, to=1400, by=400))
36: 
37: fname <- 'fig4-8-right.png'
38: ggsave(file=fname, plot=p, dpi=70)

data {
  int N;
  real X[N];
  real Y[N];
  // 追加部分
  int N_new;
  real X_new[N_new];
}

parameters {
  real a;
  real b;
  real<lower=0> sigma;
}

// transformed parameters では,
// data, parameter を組み合わせて,
// 新しいサンプリングパラメータを作成出来る.
transformed parameters {
  real y_base[N];
  for (i in 1:N) {
    y_base[i] = a + b*X[i];
  }
}

model {
  for (i in 1:N) {
    Y[i] ~ normal(y_base[i],sigma);
  }
}

// generated quantities では,
// data, parameter, transformed parameter を使って,
// 新しいサンプリングパラメータを作成出来る.
// 分布に従う乱数作成の際には, ~ ではなく
// y = 分布名_rng(parameter) という形で = で代入する点に注意する.
generated quantities {
  real y_base_new[N_new];
  real y_new[N_new];
  for (i in 1:N_new) {
    y_base_new[i] = a + b*X_new[i];
    y_new[i]      = normal_rng(y_base_new[i],sigma);
  }
}

library(cmdstanr)
d     <- read.csv(file="data-salary.txt")
X_new <- 23:60
data  <- list(N=nrow(d),X=d$X,Y=d$Y,N_new=length(X_new),X_new=X_new)
model <- cmdstan_model('model4-4.stan')
fit   <- model$sample(data=data, chains=4, parallel_chains=4, seed=1234)
fit$save_output_files(dir="./", basename='model4-4', timestamp=FALSE, random=FALSE)

1: library(cmdstanr)
2: files <- c('model4-4-1.csv',
3:            'model4-4-2.csv',
4:            'model4-4-3.csv',
5:            'model4-4-4.csv')
6: fit <- as_cmdstan_fit(files)
7: df <- fit$summary(variables=c('y_base_new'),'mean','sd', 'quantile', 'rhat','ess_bulk','ess_tail')
8: head(df)

 1: library(cmdstanr)
 2: q95 <- function(x) quantile(x, probs=c(0.025, 0.25, 0.5, 0.75, 0.975))
 3: files <- c('model4-4-1.csv',
 4:            'model4-4-2.csv',
 5:            'model4-4-3.csv',
 6:            'model4-4-4.csv')
 7: fit <- as_cmdstan_fit(files)
 8: df <- fit$summary(variables=c('y_base_new'),'mean','sd', q95)
 9: X_new <- 23:60
10: df$X_new <- X_new
11: print(df)
12: print(class(df$"50%"))  # 書き出すときはこれでもOKだが.
13: print(df$`50%`)  # ggplot の aes で使うときは `50%` (back quote?) を使う
14: print(X_new)

 1: library(cmdstanr)
 2: library(ggplot2)
 3: 
 4: # 95% 信頼区間
 5: q95 <- function(x) quantile(x, probs=c(0.025, 0.25, 0.5, 0.75, 0.975))
 6: files <- c('model4-4-1.csv',
 7:            'model4-4-2.csv',
 8:            'model4-4-3.csv',
 9:            'model4-4-4.csv')
10: fit <- as_cmdstan_fit(files)
11: df <- fit$summary(variables=c('y_base_new'),'mean','sd', q95)
12: X_new <- 23:60
13: df$X_new <- X_new
14: 
15: # 元データ
16: d     <- read.csv(file="data-salary.txt")
17: 
18: # 図
19: p <- ggplot() + geom_point(data=d, mapping=aes(x=X,y=Y), shape=1, size=3)
20: p <- p + geom_line(data=df,mapping=aes(x=X_new,y=`50%`))
21: p <- p + geom_ribbon(data=df, mapping=aes(x=X_new, ymax=`97.5%`,ymin=`2.5%`),alpha=1/6)
22: p <- p + geom_ribbon(data=df, mapping=aes(x=X_new, ymax=`75%`  ,ymin=`25%` ),alpha=2/6)
23: 
24: fname <- 'fig4-8-left_cmdstan.png'
25: ggsave(file=fname, plot=p, dpi=70)

 1: library(cmdstanr)
 2: library(ggplot2)
 3: 
 4: # 95% 信頼区間
 5: q95 <- function(x) quantile(x, probs=c(0.025, 0.25, 0.5, 0.75, 0.975))
 6: files <- c('model4-4-1.csv',
 7:            'model4-4-2.csv',
 8:            'model4-4-3.csv',
 9:            'model4-4-4.csv')
10: fit <- as_cmdstan_fit(files)
11: df <- fit$summary(variables=c('y_new'),'mean','sd', q95)
12: X_new <- 23:60
13: df$X_new <- X_new
14: 
15: # 元データ
16: d     <- read.csv(file="data-salary.txt")
17: 
18: # 図
19: p <- ggplot() + geom_point(data=d, mapping=aes(x=X,y=Y), shape=1, size=3)
20: p <- p + geom_line(data=df,mapping=aes(x=X_new,y=`50%`))
21: p <- p + geom_ribbon(data=df, mapping=aes(x=X_new, ymax=`97.5%`,ymin=`2.5%`),alpha=1/6)
22: p <- p + geom_ribbon(data=df, mapping=aes(x=X_new, ymax=`75%`  ,ymin=`25%` ),alpha=2/6)
23: 
24: fname <- 'fig4-8-right_cmdstan.png'
25: ggsave(file=fname, plot=p, dpi=70)

 1: import pandas as pd
 2: from cmdstanpy import CmdStanModel
 3: 
 4: d = pd.read_csv("data-salary.txt",header=0)
 5: N = len(d)
 6: X = d.iloc[:,0]
 7: Y = d.iloc[:,1]
 8: X_new = range(23,61)
 9: N_new = len(X_new)
10: data = {'N': N, 'X': X,  'Y': Y, 'N_new': N_new, 'X_new': X_new}
11: model = CmdStanModel(stan_file='model4-4.stan',model_name='model4-4')
12: fit = model.sample(data=data, chains=4, parallel_chains=4, seed=1234)
13: # ファイル名を指定できない... timestamp とか要らないんだが.
14: fit.save_csvfiles(dir='./')
15: fpd = fit.summary()
16: print(fpd.head(50))

 1: import cmdstanpy
 2: # 何か面倒
 3: files=['model4-4-20220331192118_1.csv',
 4:        'model4-4-20220331192118_2.csv',
 5:        'model4-4-20220331192118_3.csv',
 6:        'model4-4-20220331192118_4.csv']
 7: fit = cmdstanpy.from_csv(path=files, method='sample')
 8: df = fit.summary(percentiles=(2.5, 50, 97.5))
 9: print(type(df))   # pandas dataframe として保存されている様子.
10: print(df.head(5))

1: library(cmdstanr)
2: files <-c('model4-4-20220331192118_1.csv',
3:           'model4-4-20220331192118_2.csv',
4:           'model4-4-20220331192118_3.csv',
5:           'model4-4-20220331192118_4.csv')
6: fit <- as_cmdstan_fit(files)
7: fit$summary()

 1: import numpy as np
 2: import cmdstanpy
 3: import matplotlib.pyplot as plt
 4: from scipy.stats import kde
 5: files=['model4-4-20220331192118_1.csv',
 6:        'model4-4-20220331192118_2.csv',
 7:        'model4-4-20220331192118_3.csv',
 8:        'model4-4-20220331192118_4.csv']
 9: fit = cmdstanpy.from_csv(path=files, method='sample')
10: 
11: la = fit.stan_variables()  # { a:[...], b:[...], y_base:[[....],[....],....] }
12: adraw = la['a']
13: bdraw = la['b']
14: sdraw = la['sigma']
15: 
16: adens = kde.gaussian_kde(adraw)
17: x = np.linspace(-350,150,100)
18: y = adens(x)
19: 
20: bdens = kde.gaussian_kde(bdraw)
21: bx = np.linspace(15,30,100)
22: by = bdens(bx)
23: 
24: sdens = kde.gaussian_kde(sdraw)
25: sx = np.linspace(50,170,100)
26: sy = sdens(sx)
27: 
28: fig = plt.figure(dpi=100, figsize=(7,3))
29: ax = fig.add_subplot(1,3,1)
30: ax.set_title('a')
31: ax.plot(x,y)
32: ax2 = fig.add_subplot(1,3,2)
33: ax2.set_title('b')
34: ax2.plot(bx,by)
35: ax3 = fig.add_subplot(1,3,3)
36: ax3.set_title('sigma')
37: ax3.plot(sx,sy)
38: 
39: # plt.rcParams['font.size'] = 5
40: plt.tight_layout()
41: 
42: fname = 'model4-5_param_dens_py.png'
43: fig.show()
44: fig.savefig(fname)

 1: import numpy as np
 2: import cmdstanpy
 3: import matplotlib.pyplot as plt
 4: from scipy.stats import kde
 5: files=['model4-4-20220331192118_1.csv',
 6:        'model4-4-20220331192118_2.csv',
 7:        'model4-4-20220331192118_3.csv',
 8:        'model4-4-20220331192118_4.csv']
 9: fit = cmdstanpy.from_csv(path=files, method='sample')
10: 
11: la = fit.stan_variables()  # { a:[...], b:[...], y_base:[[....],[....],....] }
12: adraw = la['a']
13: bdraw = la['b']
14: sdraw = la['sigma']
15: 
16: fig = plt.figure(dpi=100, figsize=(7,3))
17: ax = fig.add_subplot(1,3,1)
18: ax.set_title('a')
19: ax.hist(adraw)
20: ax2 = fig.add_subplot(1,3,2)
21: ax2.set_title('b')
22: ax2.hist(bdraw)
23: ax3 = fig.add_subplot(1,3,3)
24: ax3.set_title('sigma')
25: ax3.hist(sdraw)
26: 
27: # plt.rcParams['font.size'] = 5
28: plt.tight_layout()
29: 
30: fname = 'model4-5_param_hist_py.png'
31: fig.show()
32: fig.savefig(fname)

 1: # テストプログラム
 2: import numpy as np
 3: import cmdstanpy
 4: import matplotlib.pyplot as plt
 5: from scipy.stats import kde
 6: files=['model4-4-20220331192118_1.csv',
 7:        'model4-4-20220331192118_2.csv',
 8:        'model4-4-20220331192118_3.csv',
 9:        'model4-4-20220331192118_4.csv']
10: fit = cmdstanpy.from_csv(path=files, method='sample')
11: 
12: la = fit.stan_variables()  # { a:[...], b:[...], y_base:[[....],[....],....] }
13: adraw = la['a']
14: xmax = len(adraw)
15: print(xmax)
16: rng = range(0, xmax)
17: print(len(rng))

 1: import numpy as np
 2: import cmdstanpy
 3: import matplotlib.pyplot as plt
 4: from scipy.stats import kde
 5: files=['model4-4-20220331192118_1.csv',
 6:        'model4-4-20220331192118_2.csv',
 7:        'model4-4-20220331192118_3.csv',
 8:        'model4-4-20220331192118_4.csv']
 9: fit = cmdstanpy.from_csv(path=files, method='sample')
10: 
11: la = fit.stan_variables()  # { a:[...], b:[...], y_base:[[....],[....],....] }
12: adraw = la['a']
13: bdraw = la['b']
14: sdraw = la['sigma']
15: x = range(0, len(adraw))
16: 
17: fig = plt.figure(dpi=100, figsize=(7,3))
18: ax = fig.add_subplot(1,3,1)
19: ax.set_title('a')
20: ax.plot(x,adraw)
21: ax2 = fig.add_subplot(1,3,2)
22: ax2.set_title('b')
23: ax2.plot(x,bdraw)
24: ax3 = fig.add_subplot(1,3,3)
25: ax3.set_title('sigma')
26: ax3.plot(x,sdraw)
27: 
28: # plt.rcParams['font.size'] = 5
29: plt.tight_layout()
30: 
31: fname = 'model4-5_param_trace_py.png'
32: fig.show()
33: fig.savefig(fname)

 1: # テストプログラム
 2: import numpy as np
 3: import cmdstanpy
 4: import matplotlib.pyplot as plt
 5: import pandas as pd
 6: d = pd.read_csv("data-salary.txt",header=0)
 7: print(len(d['X']))
 8: files=['model4-4-20220331192118_1.csv',
 9:        'model4-4-20220331192118_2.csv',
10:        'model4-4-20220331192118_3.csv',
11:        'model4-4-20220331192118_4.csv']
12: fit = cmdstanpy.from_csv(path=files, method='sample')
13: 
14: la = fit.stan_variables()  # { a:[...], b:[...], y_base:[[....],[....],....] }
15: #print(type(la))   # dict
16: ybase = la['y_base']
17: 
18: p = np.percentile(ybase[:,0], q=[2.5, 25, 50, 75, 97.5]) # [329.275975 382.54375  407.4335   432.47275  484.563225]
19: #print(p)
20: df = pd.DataFrame()
21: i = 42
22: df[23] = np.percentile(ybase[:,0], q=[2.5, 25, 50, 75, 97.5])
23: df[i] = np.percentile(ybase[:,i-23], q=[2.5, 25, 50, 75, 97.5])
24: # 転置
25: df2 = df.T
26: # 列名の設定
27: df3 = df2.set_axis(['p2.5','p25','p50','p75','p97.5'], axis='columns')
28: print(df3)

 1: # テストプログラム2
 2: import numpy as np
 3: import cmdstanpy
 4: import matplotlib.pyplot as plt
 5: import pandas as pd
 6: d = pd.read_csv("data-salary.txt",header=0)
 7: X = d['X']
 8: files=['model4-4-20220331192118_1.csv',
 9:        'model4-4-20220331192118_2.csv',
10:        'model4-4-20220331192118_3.csv',
11:        'model4-4-20220331192118_4.csv']
12: fit = cmdstanpy.from_csv(path=files, method='sample')
13: 
14: la = fit.stan_variables()  # { a:[...], b:[...], y_base:[[....],[....],....] }
15: ybase = la['y_base']
16: df = pd.DataFrame()
17: for i in range(0,len(X)):
18:     df[i] = np.percentile(ybase[:,i], q=[2.5, 25, 50, 75, 97.5])
19: df2 = df.T
20: df3 = df2.set_axis(['p2.5','p25','p50','p75','p97.5'], axis='columns')
21: print(df3)

 1: import numpy as np
 2: import cmdstanpy
 3: import matplotlib.pyplot as plt
 4: import pandas as pd
 5: d = pd.read_csv("data-salary.txt",header=0)
 6: X = d['X']
 7: Y = d['Y']
 8: files=['model4-4-20220331192118_1.csv',
 9:        'model4-4-20220331192118_2.csv',
10:        'model4-4-20220331192118_3.csv',
11:        'model4-4-20220331192118_4.csv']
12: fit = cmdstanpy.from_csv(path=files, method='sample')
13: 
14: la = fit.stan_variables()  # { a:[...], b:[...], y_base:[[....],[....],....] }
15: X_new = range(23,61)
16: ybase = la['y_base_new']
17: df = pd.DataFrame()
18: for i in range(0,len(X_new)):
19:     df[i] = np.percentile(ybase[:,i], q=[2.5, 25, 50, 75, 97.5])
20: df2 = df.T
21: df3 = df2.set_axis(['p2.5','p25','p50','p75','p97.5'], axis='columns')
22: 
23: # d['X'] と df3 を使う
24: 
25: fig = plt.figure(dpi=100, figsize=(5,5))
26: ax = fig.add_subplot(1,1,1)
27: ax.scatter(X,Y, marker='o', s=50, facecolor='None', edgecolors='black')
28: ax.fill_between(X_new, df3['p2.5'], df3['p97.5'], color='gray',alpha=0.2)
29: ax.fill_between(X_new, df3['p25'], df3['p75'], color='gray',alpha=0.5)
30: plt.tight_layout()
31: 
32: fname = 'fig4-8-left_py.png'
33: fig.show()
34: fig.savefig(fname)

 1: import numpy as np
 2: import cmdstanpy
 3: import matplotlib.pyplot as plt
 4: import pandas as pd
 5: files=['model4-4-20220331192118_1.csv',
 6:        'model4-4-20220331192118_2.csv',
 7:        'model4-4-20220331192118_3.csv',
 8:        'model4-4-20220331192118_4.csv']
 9: fit = cmdstanpy.from_csv(path=files, method='sample')
10: 
11: la = fit.stan_variables()  # { a:[...], b:[...], y_base:[[....],[....],....] }
12: X_new = range(23,61)
13: 
14: y_new = la['y_new']
15: print(type(y_new))   # numpy.ndarray
16: print(len(y_new))    # 4000 列 (横に乱数が並んでいる感じ)
17: print(len(y_new[0])) # 38 行
18: df = pd.DataFrame()
19: for i in range(0,len(X_new)):
20:     df[i] = np.percentile(y_new[:,i], q=[2.5, 25, 50, 75, 97.5])
21: df2 = df.T
22: df3 = df2.set_axis(['p2.5','p25','p50','p75','p97.5'], axis='columns')
23: print(df3.head(5))

 1: import numpy as np
 2: import cmdstanpy
 3: import matplotlib.pyplot as plt
 4: import pandas as pd
 5: d = pd.read_csv("data-salary.txt",header=0)
 6: X = d['X']
 7: Y = d['Y']
 8: files=['model4-4-20220331192118_1.csv',
 9:        'model4-4-20220331192118_2.csv',
10:        'model4-4-20220331192118_3.csv',
11:        'model4-4-20220331192118_4.csv']
12: fit = cmdstanpy.from_csv(path=files, method='sample')
13: 
14: la = fit.stan_variables()  # { a:[...], b:[...], y_base:[[....],[....],....] }
15: X_new = range(23,61)
16: 
17: y_new = la['y_new']
18: df = pd.DataFrame()
19: for i in range(0,len(X_new)):
20:     df[i] = np.percentile(y_new[:,i], q=[2.5, 25, 50, 75, 97.5])
21: df2 = df.T
22: df3 = df2.set_axis(['p2.5','p25','p50','p75','p97.5'], axis='columns')
23: 
24: # X_new と df3 を使う
25: 
26: fig = plt.figure(dpi=100, figsize=(5,5))
27: ax = fig.add_subplot(1,1,1)
28: ax.scatter(X,Y, marker='o', s=50, facecolor='None', edgecolors='black')
29: ax.fill_between(X_new, df3['p2.5'], df3['p97.5'], color='gray',alpha=0.2)
30: ax.fill_between(X_new, df3['p25'], df3['p75'],    color='gray',alpha=0.5)
31: plt.tight_layout()
32: 
33: fname = 'fig4-8-right_py.png'
34: fig.show()
35: fig.savefig(fname)