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%です。
