ベイズ統計： 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%です。