ベイズ統計： pystanを使う 〜 指数分布を用いた推測

朝倉書店の「基礎からのベイズ統計学」7章末問題2)をStanを使って解きます。
Sさんは桜並木を散歩しています。桜並木には普通の桜に混じって時々八重桜が咲いていて、八重桜が1本観察されるまでに歩いた距離を10回記録すると、
189.9
71.3
243.1
94.1
204.9
122.0
163.0
467.1
41.8
36.4
となりました（単位はメートル）。八重桜が1本観察されるまでに歩いた距離が指数分布に従うとすると、八重桜は平均的に何mおきに生えているでしょうか。また、これから100m歩く間に八重桜を1本以上見る確率はどれくらいでしょうか。そして、その確率が50%より高い確率はどれくらいでしょうか。

コードはこちらになります。

# -*- coding: utf-8 -*-
import pandas as pd
import numpy as np
import pystan

def get_data():
  x = np.array(pd.read_csv("ch0702.data", header=None)[0])
  print(x)

  return {
    "n" : x.size,
    "x" : x
  }

def get_code():
  f = open("ch0702.stan")
  code = f.read()
  f.close()
  return code

def main():
  stan_data = get_data()
  stan_code = get_code()

  fit = pystan.stan(model_code=stan_code,
                    data=stan_data,
                    iter=11000,
                    warmup=1000,
                    thin=1,
                    chains=4)
  print(fit)
  fit.plot()


if __name__ == "__main__":
  main()

データです。

189.9
71.3
243.1
94.1
204.9
122.0
163.0
467.1
41.8
36.4

Stanのコードです。

  data {
    int<lower=0>  n;
    real<lower=0>  x[n];
  }
  parameters {
    real<lower=0>  lambda;
  }
  model {
    x ~ exponential(lambda);
  }
  generated quantities {
    real mu;
    real<lower=0, upper=1> p100;
    real<lower=0, upper=1> p100_over50;
    mu <- inv(lambda);
    p100 <- exponential_cdf(100, lambda);
    p100_over50 <- step(p100 - .5);
  }

出力結果になります。

              mean se_mean     sd   2.5%    25%    50%    75%  97.5%  n_eff   Rhat
lambda      6.7e-3  1.0e-5 2.0e-3 3.4e-3 5.3e-3 6.5e-3 7.9e-3   0.01  40000    1.0
mu          163.39    0.27  54.04  89.25  126.1 152.92 189.47 296.69  40000    1.0
p100          0.48  5.0e-4    0.1   0.29   0.41   0.48   0.55   0.67  40000    1.0
p100_over50   0.42  2.5e-3   0.49    0.0    0.0    0.0    1.0    1.0  40000    1.0
lp__        -66.51  3.6e-3   0.71 -68.58 -66.67 -66.23 -66.06 -66.01  40000    1.0

muのEAP推定値から、八重桜は平均163.4mおきに生えていると言えます(確信区間は89.25mから296.69m)。100m歩く間に八重桜を1本以上見る確率(p100)は、48%となりました。ほぼ五分五分の確率です。その確率が50%を超える確率は42%です。

ベイズ統計： rstanを使う 〜 正規分布の対応のある2群の平均値差に関する推測

朝倉書店の「基礎からのベイズ統計学」6章末問題5)をStanを使って解きます。あるダイエット方法を20人が試して、ダイエット前後に平均値に差があるか検討します。平均値差のEAP推定値と2群の平均値差が0kgより大きいという研究仮説と、2kgより大きいという研究仮説が正しい確率を求めます。加えて相関係数が0.7より大きい確率も求めます。

今回はpystanではなく、rstanを使います。pystanでは実行中にmemory errorが発生してプログラムが完了しませんでした。stan実行時のchainsパラメータを小さくしましたが変わりませんでした。

コードはこちらになります。

library(rstan)
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())

setwd("~/R/repository/bayes_hmc")

stan_file <- "ch0605.stan"
x=read.csv(file="./ch0605.data", header=T)
stan_data <- list(n=nrow(x), x=x)
par<-c("mu","stdev","rho","stdevsq","sigma",
       "delta","over0","over2", "rho_over0", "rho_over07")

fit <- stan(file=stan_file,
            data=stan_data,
            iter=11000,
            warmup=1000,
            thin=1,
            chains=4)

print(fit, digits_summary=3)
plot(fit, pars=par)
traceplot(fit, inc_warmup=F, pars=par)

stanのコード(ch0605.stanの内容)です。

  data {
    int<lower=0>  n;
    vector[2]     x[n];
  }
  parameters {
    vector[2]               mu;     # mean
    vector<lower=0>[2]      stdev;  # standard deviation
    real<lower=-1, upper=1> rho;    # correlation coefficient
  }
  transformed parameters {
    vector<lower=0>[2] stdevsq;   # variance
    matrix[2, 2] sigma;
    stdevsq[1] <- pow(stdev[1], 2);
    stdevsq[2] <- pow(stdev[2], 2);
    sigma[1, 1] <- stdevsq[1];
    sigma[2, 2] <- stdevsq[2];
    sigma[1, 2] <- stdev[1] * stdev[2] * rho;
    sigma[2, 1] <- stdev[1] * stdev[2] * rho;
  }
  model {
    x ~ multi_normal(mu, sigma);    # multi normal distribution model
  }
  generated quantities {
    real delta;
    real<lower=0> over0;
    real<lower=0> over2;
    real<lower=0> rho_over0;
    real<lower=0> rho_over07;
    delta <- mu[1] - mu[2];
    over0 <- step(delta);
    over2 <- step(delta-2);
    rho_over0 <- step(rho);
    rho_over07 <- step(rho-0.7);
  }

データです。

prior, post
53.1, 51.3
51.5, 48.2
45.5, 49.6
55.5, 59.6
49.6, 44.2
50.1, 47.6
59.2, 54.9
54.7, 56.5
53.0, 48.4
48.6, 50.6
55.3, 53.6
54.6, 57.5
51.7, 52.0
48.6, 46.9
56.4, 56.8
58.9, 50.0
53.3, 53.8
42.4, 38.3
51.9, 52.6
39.1, 41.0

実行結果はこうなりました。

              mean se_mean     sd    2.5%     25%     50%     75%   97.5% n_eff Rhat
mu[1]       51.625   0.011  1.228  49.197  50.818  51.627  52.423  54.074 13241    1
mu[2]       50.638   0.011  1.312  48.015  49.791  50.640  51.489  53.197 13317    1
stdev[1]     5.434   0.009  0.950   3.956   4.758   5.306   5.962   7.659 12272    1
stdev[2]     5.817   0.009  1.018   4.225   5.095   5.687   6.385   8.192 12752    1
rho          0.744   0.001  0.108   0.482   0.687   0.763   0.823   0.901 14120    1
stdevsq[1]  30.433   0.103 11.255  15.648  22.639  28.153  35.548  58.657 11871    1
stdevsq[2]  34.874   0.117 12.937  17.853  25.963  32.338  40.774  67.103 12215    1
sigma[1,1]  30.433   0.103 11.255  15.648  22.639  28.153  35.548  58.657 11871    1
sigma[1,2]  24.487   0.106 10.469  10.325  17.305  22.467  29.341  50.625  9806    1
sigma[2,1]  24.487   0.106 10.469  10.325  17.305  22.467  29.341  50.625  9806    1
sigma[2,2]  34.874   0.117 12.937  17.853  25.963  32.338  40.774  67.103 12215    1
delta        0.987   0.005  0.899  -0.784   0.400   0.979   1.570   2.765 36824    1
over0        0.872   0.002  0.334   0.000   1.000   1.000   1.000   1.000 29392    1
over2        0.124   0.002  0.329   0.000   0.000   0.000   0.000   1.000 28841    1
rho_over0    1.000   0.000  0.005   1.000   1.000   1.000   1.000   1.000 40000    1
rho_over07   0.714   0.003  0.452   0.000   0.000   1.000   1.000   1.000 18164    1
lp__       -76.711   0.015  1.666 -80.816 -77.580 -76.375 -75.479 -74.500 12078    1

平均値差(delta)は0.987で差があるように見えます。平均値差が0kgより大きい確率(over0)は87.2%で差があることは確かなようです。平均値差が2kgより大きい確率(over2)は12.4%となり、非常に効果があるダイエット方法とは言えません。
また、相関係数が0.7より大きい確率(rho_over07)は71.4%ですので、中程度から強い相関はありそうです。

ベイズ統計： pystanを使う 〜 正規分布の独立な2群の平均値差に関する推測

朝倉書店の「基礎からのベイズ統計学」6章末問題4)をStanを使って解きます。それぞれ20個のデータをもつ2群において、分散が同一であるという仮定のもとで、平均値に差があるか検討します。

コードはこちらになります。

# -*- coding: utf-8 -*-
import pystan

# karaoke score in last year contest
experiment = [30.86, 29.75, 31.55, 32.29, 29.90, 31.71, 31.35, 29.03, 30.37, 31.55,
              29.26, 32.29, 29.90, 30.18, 30.72, 32.28, 30.72, 29.90, 31.55, 31.55]
control = [31.36, 33.34, 33.16, 31.36, 36.19, 29.80, 31.11, 35.23, 31.36, 31.27,
           31.63, 31.63, 32.00, 31.11, 31.63, 31.36, 31.81, 31.63, 29.21, 33.37]
stan_data = {
  "exN" : len(experiment),
  "ex" : experiment,
  "conN" : len(control),
  "con" : control
}

stan_code = """
  data {
    int<lower=0> exN;
    real<lower=0> ex[exN];
    int<lower=0> conN;
    real<lower=0> con[conN];
  }
  parameters {
    real<lower=0> mu_ex;    # mean for experiment
    real<lower=0> mu_con;   # mean for control
    real<lower=0> stdev;    # standard deviation
  }
  model {
    ex ~ normal(mu_ex, stdev);    # normal distribution model
    con ~ normal(mu_con, stdev);  # normal distribution model
  }
  generated quantities {
    real delta;
    real<lower=0> over0;
    real<lower=0> over1;
    delta <- mu_con - mu_ex;
    over0 <- step(delta);
    over1 <- step(delta-1);
  }
"""

def main():
  fit = pystan.stan(model_code=stan_code, 
                    data=stan_data,
                    iter=11000,
                    warmup=1000,
                    thin=1,
                    chains=4)
  print(fit)
  fit.plot()

if __name__ == "__main__":
  main()

実行結果はこうなりました。

         mean se_mean     sd   2.5%    25%    50%    75%  97.5%  n_eff   Rhat
mu_ex   30.84  1.6e-3   0.32  30.22  30.63  30.83  31.05  31.46  40000    1.0
mu_con  31.98  1.6e-3   0.31  31.36  31.77  31.97  32.18  32.59  40000    1.0
stdev     1.4  8.3e-4   0.17   1.12   1.28   1.38    1.5   1.77  40000    1.0
delta    1.14  2.2e-3   0.45   0.27   0.84   1.14   1.43   2.02  40000    1.0
over0    0.99  3.8e-4   0.08    1.0    1.0    1.0    1.0    1.0  40000    1.0
over1    0.62  2.4e-3   0.49    0.0    0.0    1.0    1.0    1.0  40000    1.0
lp__   -25.44  6.3e-3   1.26  -28.7 -26.02 -25.12 -24.52 -24.01  40000    1.0

平均の差(delta)のEAP推定値は1.14となりました。事後標準偏差は0.45です。95%確信区間は0.27から2.02です。差がある確率(over0)は99%で、ほぼ差があると見ていいようです。差が1以上ある確率(over1)は62%です。

ベイズ統計： pystanを使う 〜 正規分布の分位数に関する推測

朝倉書店の「基礎からのベイズ統計学」6章末問題3)をStanを使って解きます。30個の去年のカラオケ大会データから0.90分位数のEAP (Expected A Posteriori)推定値と平均87点、標準偏差5点の人が本選に進める確率を求めます。

コードはこちらになります。

# -*- coding: utf-8 -*-
import pystan

# karaoke score in last year contest
karaoke = [75, 82, 77, 86, 98, 91, 85, 84, 82, 79,
           88, 79, 84, 87, 69, 93, 84, 82, 89, 78,
           90, 74, 75, 84, 89, 81, 74, 79, 82, 84]
stan_data = {
  "n" : len(karaoke),
  "x" : karaoke
}

stan_code = """
  data {
    int<lower=0> n;       # number of data
    real<lower=0> x[n];   # vector of karaoke score
  }
  parameters {
    real<lower=0> mu;     # mean
    real<lower=0> stdev;  # standard deviation
  }
  model {
    x ~ normal(mu, stdev);  # normal distribution model
  }
  generated quantities {
    real<lower=0> borderline;             # 90% percentile
    real<lower=0, upper=1> qualifying;    # probability of qualifying karaoke contest
    borderline <- mu + 1.282 * stdev;
    qualifying <- 1 - normal_cdf(borderline, 87, 5);
  }
"""

def main():
  fit = pystan.stan(model_code=stan_code, 
                    data=stan_data,
                    iter=11000,
                    warmup=1000,
                    thin=1,
                    chains=4)
  print(fit)
  fit.plot()

if __name__ == "__main__":
  main()

実行結果はこうなりました。

             mean se_mean     sd   2.5%    25%    50%    75%  97.5%  n_eff   Rhat
mu          82.79  6.1e-3   1.22   80.4  81.97  82.79  83.61  85.18  40000    1.0
stdev        6.64  4.6e-3   0.92   5.13   5.99   6.54   7.18   8.72  40000    1.0
borderline   91.3  8.5e-3   1.71  88.36  90.11  91.17  92.33  95.07  40000    1.0
qualifying   0.21  4.4e-4   0.09   0.05   0.14    0.2   0.27   0.39  40000    1.0
lp__       -64.73  5.0e-3   1.01 -67.42 -65.14 -64.43 -64.01 -63.73  40000    1.0

0.90分位数EAP推定値は91.3点となりました。事後標準偏差は1.71点です。95%確信区間は88.36点から95.07点です。
平均87点、標準偏差5点の人が本選に進める確率は21%です。

ベイズ統計： pystanを使う 〜 正規分布の分散に関する推測

朝倉書店の「基礎からのベイズ統計学」６章末問題2)をStanを使って解きます。20個のデータから分散のEAP (Expected A Posteriori)推定値と分散が1.5以下である確率、1.0以下である確率を求めます。

コードはこちらになります。

# -*- coding: utf-8 -*-
import pystan

juryo = [128.27, 130.00, 130.68, 128.82, 130.08, 129.24, 129.82,
         131.76, 128.19, 131.01, 128.81, 131.53, 129.36, 130.82,
         130.89, 130.35, 130.18, 131.00, 131.03, 128.84]
stan_data = {
  "n" : len(juryo),
  "x" : juryo
}

stan_code = """
  data {
    int<lower=0> n;       # number of data
    real<lower=0> x[n];   # vector of juryo
  }
  parameters {
    real<lower=0> stdev;  # standard deviation
  }
  transformed parameters {
    real<lower=0> stdevsq;  # variance
    stdevsq <- pow(stdev, 2);
  }
  model {
    x ~ normal(130.0, stdev);  # normal distribution model
  }
  generated quantities {
    real<lower=0, upper=1> under0;   # probability of EAP under 1.5
    real<lower=0, upper=1> under1;   # probability of EAP under 1.0
    under0 <- step(1.5-stdevsq);
    under1 <- step(1.0-stdevsq);
  }
"""

def main():
  fit = pystan.stan(model_code=stan_code, 
                    data=stan_data,
                    iter=11000,
                    warmup=1000,
                    thin=1,
                    chains=4)
  print(fit)
  fit.plot()

if __name__ == "__main__":
  main()

実行結果はこうなりました。

          mean se_mean     sd   2.5%    25%    50%    75%  97.5%  n_eff   Rhat
stdev     1.12  9.5e-4   0.19   0.82   0.98   1.09   1.23   1.56  40000    1.0
stdevsq   1.29  2.3e-3   0.46   0.67   0.97    1.2   1.51   2.43  40000    1.0
under0    0.75  2.2e-3   0.44    0.0    0.0    1.0    1.0    1.0  40000    1.0
under1    0.29  2.3e-3   0.45    0.0    0.0    0.0    1.0    1.0  40000    1.0
lp__    -11.41  3.5e-3    0.7 -13.42 -11.57 -11.14 -10.97 -10.92  40000    1.0

分散のEAP推定値は1.29です。この分散の事後標準偏差は0.46です。
1.5以下である確率は75%、1.0以下である確率は29%となりました。

ベイズ統計： pystanを使う 〜 正規分布の平均に関する推測

朝倉書店の「基礎からのベイズ統計学」６章末問題1)をStanを使って解きます。１２個のデータから平均のEAP (Expected A Posteriori)推定値と平均が70,000円を超える確率と75,000円を超える確率を求めます。

コードはこちらになります。

# -*- coding: utf-8 -*-
import pystan

uriage = [76230, 73550, 80750, 71500, 75420, 74840,
          71580, 76920, 68450, 76990, 64070, 76200]
stan_data = {
  "n" : len(uriage),
  "x" : uriage
}

stan_code = """
  data {
    int<lower=0> n;         # number of data
    real<lower=0> x[n];   # vector of uriage
  }
  parameters {
    real<lower=0> mu;     # mean
    real<lower=0> stdev;  # standard deviation
  }
  model {
    x ~ normal(mu, stdev);  # normal distribution model
  }
  generated quantities {
    real<lower=0, upper=1> over0;   # probability of EAP over 70k JPY
    real<lower=0, upper=1> over1;   # probability of EAP over 75k JPY
    over0 <- step(mu-70000);
    over1 <- step(mu-75000);
  }
"""

def main():
  fit = pystan.stan(model_code=stan_code, 
                    data=stan_data,
                    iter=11000,
                    warmup=1000,
                    thin=1,
                    chains=4)
  print(fit)
  fit.plot()

if __name__ == "__main__":
  main()

実行結果はこうなりました。

        mean se_mean     sd   2.5%    25%    50%    75%  97.5%  n_eff   Rhat
mu     7.4e4    7.53 1506.7  7.1e4  7.3e4  7.4e4  7.5e4  7.7e4  40000    1.0
stdev 5041.7    6.26 1252.8 3274.0 4166.7 4829.6 5655.0 8089.9  40000    1.0
over0   0.99  4.2e-4   0.08    1.0    1.0    1.0    1.0    1.0  40000    1.0
over1   0.21  2.0e-3   0.41    0.0    0.0    0.0    0.0    1.0  40000    1.0
lp__  -87.76  5.4e-3   1.08 -90.63 -88.18 -87.43 -86.99 -86.71  40000    1.0

EAP推定値は74,000円となりました。事後標準偏差は1,506円です。95%確信区間は71,000円から77,000円です。
70,000円を超える確率は99%です。
75,000円を超える確率は21%です。

Pythonでデータサイエンス入門〜クロス集計：どの属性の顧客が離脱しているのか

ビジネス活用事例で学ぶデータサイエンス入門 (SBクリエイティブ）では、Rでコーディングされたプログラムを使ってデータ分析を行っています。これをPythonでやるとどうなるのか、Python初心者、データサイエンス初心者がやってみました。
コーディングの環境はminiconda3でインストールしたSpyderを使ってます。

今回は第4章「どの属性の顧客が離脱しているのか？」(クロス集計）をPythonでコーディングします。

# -*- coding: utf-8 -*-
import os
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

"""
ファイル読込
"""
current_dir = os.getcwd()
data_dir = "data"
dau_file = "section4-dau.csv"
userinfo_file = "section4-user_info.csv"
# DAU
dau = pd.read_csv(os.path.join(current_dir, data_dir, dau_file))
print(dau.head())
# User Info
userinfo = pd.read_csv(os.path.join(current_dir, data_dir, userinfo_file))
print(userinfo.head())
"""
ファイルを結合 : merge
"""
dau_userinfo = pd.merge(dau, userinfo,
                        how="left",
                        on=["user_id", "app_name"])
print(dau_userinfo.head())
"""
セグメント分析 : 性別でクロス集計
"""
dau_userinfo["log_month"] = dau_userinfo.log_date.str[0:7]
cross_gender = pd.crosstab(dau_userinfo.log_month, dau_userinfo.gender)
print(cross_gender)
"""
セグメント分析 : 年代でクロス集計
"""
cross_generation = pd.crosstab(dau_userinfo.log_month, dau_userinfo.generation)
print(cross_generation)
"""
セグメント分析 : 性別 x 年代で集計
"""
pivot_gender_generation = pd.pivot_table(dau_userinfo,
                                        values="user_id",
                                        index=["log_month"],
                                        columns=["gender", "generation"],
                                        aggfunc=np.count_nonzero)
print(pivot_gender_generation)
"""
セグメント分析 : デバイスのクロス集計
"""
cross_device = pd.crosstab(dau_userinfo.log_month, dau_userinfo.device_type)
print(cross_device)
"""
セグメント分析の結果を可視化する
"""
# 日付、デバイス別にユーザー数を算出する
dau_userinfo_device_summary = dau_userinfo.groupby(
                                    by=["log_date", "device_type"],
                                    as_index=False
                                    ).user_id.count()
dau_userinfo_device_summary.columns = ["log_date", "device_type", "dau"]
# log_dateを日付型に変換する
dau_userinfo_device_summary.log_date = pd.to_datetime(dau_userinfo_device_summary.log_date)
print(dau_userinfo_device_summary)
# 時系列のトレンドグラフを描画する
xy_android = dau_userinfo_device_summary[
                    dau_userinfo_device_summary.device_type =="Android"
                    ][["log_date", "dau"]]
xy_ios = dau_userinfo_device_summary[
                    dau_userinfo_device_summary.device_type =="iOS"
                    ][["log_date", "dau"]]
plt.plot(
        xy_android.log_date,
        xy_android.dau,
        linestyle="--",
        marker="o",
        markersize=10,
        color="pink",
        label="Android"
        )
plt.plot(
        xy_ios.log_date,
        xy_ios.dau,
        linestyle="-",
        marker="^",
        markersize=10,
        color="skyblue",
        label="iOS"
        )
plt.axis(facecolor="black")
plt.ylim(ymin=0)
plt.xlabel("log date")
plt.ylabel("dau")
plt.subplots_adjust(right=0.86)
plt.legend(title="device type", loc="right", bbox_to_anchor=(1.32,.5))
plt.show()

セグメント分析：時系列