黒木玄
2017-10-16
Mambaパッケージを
http://mambajl.readthedocs.io/en/latest/tutorial.html#model-specification
に書いてある方法でモデルを定義して、mcmc() 函数でシミュレーションを実行すると、モデル全体の対数尤度函数を直接的に定義して sample!() 函数の繰り返しでシミュレーションを行うよりも計算が非常に遅くなる。計算速度が約10倍ほど違う。この違いは何が原因で生じているのだろうか?
このノートでは結果的にこの疑問を解決できなかったが、サンプラーに渡す対数尤度函数の書き方によってシミュレーション速度が大きく変わることは確認できた。
## Test Sample
using Distributions
import PyPlot; plt = PyPlot
srand(32)
N = 32
X = rand(Uniform(-5.0, 5.0), N)
Y = 2.0 + 1.0*X + rand(Normal(0.0, 2.0), N)
@time begin
b0_lse, b1_lse = linreg(X,Y)
d = fit_mle(Normal, Y - b0_lse - b1_lse*X)
s2_lse = d.σ^2
end
@show b0_lse, b1_lse, s2_lse
sleep(0.2)
plt.figure(figsize=(5,3.5))
plt.scatter(X, Y, label="Sample")
x = -5.0:5.0
plt.plot(x, b0_lse + b1_lse*x, color="red", label="LSE")
plt.grid("on", ls=":")
plt.legend()
plt.title("Sample and LSE")
0.218353 seconds (99.81 k allocations: 5.062 MiB, 4.08% gc time) (b0_lse, b1_lse, s2_lse) = (1.9733544847267817, 1.0864307623437985, 3.334179811356694)
PyObject <matplotlib.text.Text object at 0x000000002F895BE0>
################################################################################
## Linear Regression
## y ~ N(b0 + b1 * x, s2)
## b0, b1 ~ N(0, 1000)
## s2 ~ invgamma(0.001, 0.001)
################################################################################
using Mamba
## Input Data
data = Dict{Symbol, Any}(
:X => X,
:Y => Y,
:N => length(X),
:b00 => 0.0,
:b10 => 0.0,
:s20 => 1.0,
)
## Model Specification
model = Model(
y = Stochastic(1,
(b0, b1, s2) -> UnivariateDistribution[
Normal(b0 + b1*data[:X][i], sqrt(s2)) for i in 1:data[:N]
],
false
),
b0 = Stochastic(
() -> Normal(0.0, 1000.0),
true
),
b1 = Stochastic(
() -> Normal(0.0, 1000.0),
true
),
s2 = Stochastic(
() -> InverseGamma(0.001, 0.001),
true
),
)
## Number of chains
chains = 1
## Initial Values
inits = [
Dict{Symbol, Any}(
:y => data[:Y],
:b0 => data[:b00],
:b1 => data[:b10],
:s2 => data[:s20],
)
for k in 1:chains
]
## Sampling Scheme Assignment
scheme = [
NUTS([:b0, :b1, :s2])
]
setsamplers!(model, scheme)
## MCMC
iterations = 10000
burnin = 2000
thinning = 1
@time mc = mcmc(model, data, inits, iterations, burnin=burnin, thin=thinning, chains=chains)
## Summary of MCMC
print("\nSummary:\n\n")
display(describe(mc))
## Convergence Diagnostics
#print("Gelman Diagnostic:\n")
#show(gelmandiag(mc))
## Plot
sleep(0.2)
draw(plot(mc), width=10inch, height=4inch, nrow=1, ncol=2, ask=false);
MCMC Simulation of 10000 Iterations x 1 Chain... Chain 1: 0% [0:21:42 of 0:21:43 remaining] Chain 1: 10% [0:00:26 of 0:00:29 remaining] Chain 1: 20% [0:00:18 of 0:00:23 remaining] Chain 1: 30% [0:00:15 of 0:00:21 remaining] Chain 1: 40% [0:00:12 of 0:00:20 remaining] Chain 1: 50% [0:00:10 of 0:00:19 remaining] Chain 1: 60% [0:00:08 of 0:00:19 remaining] Chain 1: 70% [0:00:06 of 0:00:18 remaining] Chain 1: 80% [0:00:04 of 0:00:18 remaining] Chain 1: 90% [0:00:02 of 0:00:18 remaining] Chain 1: 100% [0:00:00 of 0:00:18 remaining] 20.385970 seconds (151.33 M allocations: 3.838 GiB, 3.68% gc time) Summary: Iterations = 2001:10000 Thinning interval = 1 Chains = 1 Samples per chain = 8000 Empirical Posterior Estimates:
nothing
Mean SD Naive SE MCSE ESS s2 3.8145955 1.062420967 0.0118782275 0.0191226553 3086.7174 b1 1.0848857 0.118329961 0.0013229692 0.0016422079 5191.9819 b0 1.9646058 0.341234798 0.0038151210 0.0065682292 2699.0444 Quantiles: 2.5% 25.0% 50.0% 75.0% 97.5% s2 2.2692860 3.0871734 3.6364907 4.3645196 6.3212915 b1 0.8497315 1.0069151 1.0848232 1.1620178 1.3212111 b0 1.2916790 1.7439487 1.9664853 2.1926454 2.6359444
Mambaの通常の使い方だと 10000 iterations で20秒程度かかる
次のセルは
http://mambajl.readthedocs.io/en/latest/samplers/nuts.html
の例のほぼコピー。サンプルの与え方だけを変えてある。
logfgrad() 函数を直接定義してMCMCしている。
################################################################################
## Linear Regression
## y ~ N(b0 + b1 * x, s2)
## b0, b1 ~ N(0, 1000)
## s2 ~ invgamma(0.001, 0.001)
################################################################################
using Mamba
## Data
data = Dict(
:x => X,
:y => Y
)
## Log-transformed Posterior(b0, b1, log(s2)) + Constant and Gradient Vector
function logfgrad(x::DenseVector)
b0 = x[1]
b1 = x[2]
logs2 = x[3]
r = data[:y] - b0 - b1 * data[:x]
logf = (-0.5 * length(data[:y]) - 0.001) * logs2 -
(0.5 * dot(r, r) + 0.001) / exp(logs2) -
0.5 * b0^2 / 1000 - 0.5 * b1^2 / 1000
grad = [
sum(r) / exp(logs2) - b0 / 1000,
sum(data[:x] .* r) / exp(logs2) - b1 / 1000,
-0.5 * length(data[:y]) - 0.001 + (0.5 * dot(r, r) + 0.001) / exp(logs2)
]
logf, grad
end
## MCMC Simulation with No-U-Turn Sampling
n = 10000
burnin = 2000
sim = Chains(n, 3, start = (burnin + 1), names = ["b0", "b1", "s2"])
theta = NUTSVariate([0.0, 0.0, 0.0], logfgrad)
@time for i in 1:n
sample!(theta, adapt = (i <= burnin))
if i > burnin
sim[i, :, 1] = [theta[1:2]; exp(theta[3])]
end
end
describe(sim)
sleep(0.2)
draw(plot(sim), width=10inch, height=4inch, nrow=1, ncol=2, ask=false);
1.705615 seconds (5.76 M allocations: 266.812 MiB, 4.14% gc time) Iterations = 2001:10000 Thinning interval = 1 Chains = 1 Samples per chain = 8000 Empirical Posterior Estimates: Mean SD Naive SE MCSE ESS b0 1.9726776 0.34755550 0.0038857886 0.0069661842 2489.1921 b1 1.0862246 0.11639038 0.0013012840 0.0013473754 7462.0291 s2 3.7796685 1.03901525 0.0116165437 0.0190110668 2986.9701 Quantiles: 2.5% 25.0% 50.0% 75.0% 97.5% b0 1.2857420 1.7365435 1.9657448 2.2080103 2.6569081 b1 0.8581004 1.0081835 1.0850679 1.1632342 1.3171053 s2 2.2529309 3.0545961 3.6073976 4.3315463 6.2464812
直接 sample!() 函数を使う方法だと、同じ推定に2秒かからない!
Mambaの通常の使い方よりも10倍以上速い!
Mambaの通常の使い方はどうして遅いのだろうか?
次のセルは同じモデルを別の方法で試してみた。
こちらは gradient を Calculus.gradient で計算させた。
################################################################################
## Linear Regression
## y ~ N(b0 + b1 * x, s2)
## b0, b1 ~ N(0, 1000)
## s2 ~ invgamma(0.001, 0.001)
################################################################################
using Mamba
import Calculus; calc = Calculus
## Data
data = Dict(
:x => X,
:y => Y
)
## Log-transformed Posterior(b0, b1, log(s2)) + Constant and Gradient Vector
function loglik(x::DenseVector)
b0 = x[1]
b1 = x[2]
logs2 = x[3]
r = data[:y] - b0 - b1 * data[:x]
logf = (-0.5 * length(data[:y]) - 0.001) * logs2 -
(0.5 * dot(r, r) + 0.001) / exp(logs2) -
0.5 * b0^2 / 1000 - 0.5 * b1^2 / 1000
return logf
end
function loglikgrad(x::DenseVector)
logf = loglik(x)
grad = calc.gradient(loglik, x)
logf, grad
end
## MCMC Simulation with No-U-Turn Sampling
n = 10000
burnin = 2000
sim2 = Chains(n, 3, start = (burnin + 1), names = ["b0", "b1", "s2"])
theta = NUTSVariate([0.0, 0.0, 0.0], loglikgrad)
@time for i in 1:n
sample!(theta, adapt = (i <= burnin))
if i > burnin
sim2[i, :, 1] = [theta[1:2]; exp(theta[3])]
end
end
describe(sim2)
sleep(0.2)
draw(plot(sim2), width=10inch, height=4inch, nrow=1, ncol=2, ask=false);
0.791933 seconds (9.82 M allocations: 620.316 MiB, 10.87% gc time) Iterations = 2001:10000 Thinning interval = 1 Chains = 1 Samples per chain = 8000 Empirical Posterior Estimates: Mean SD Naive SE MCSE ESS b0 1.9653758 0.34017628 0.0038032864 0.0069244726 2413.4297 b1 1.0862501 0.11801214 0.0013194158 0.0014769405 6384.5068 s2 3.8004730 1.02505623 0.0114604771 0.0174218049 3461.8567 Quantiles: 2.5% 25.0% 50.0% 75.0% 97.5% b0 1.30641644 1.7418303 1.9656969 2.1879719 2.6336029 b1 0.85362785 1.0069135 1.0855996 1.1641017 1.3209748 s2 2.27995089 3.0790760 3.6309155 4.3448489 6.2384370
gradient を Calculus.gradient で計算させたら、倍以上速くなった!
gradient を手で計算した結果の式を使うよりも、Calculus.gradient で計算してもらった方が計算が速くなるようだ! これは意外な結果でちょっと驚き。
Distributionsパッケージの logpdf 函数を一部で使って見た。
################################################################################
## Linear Regression
## y ~ N(b0 + b1 * x, s2)
## b0, b1 ~ N(0, 1000)
## s2 ~ invgamma(0.001, 0.001)
################################################################################
using Mamba
import Calculus; calc = Calculus
## Data
data = Dict(
:x => X,
:y => Y
)
## Log-transformed Posterior(b0, b1, log(s2)) + Constant and Gradient Vector
function loglik(x::DenseVector)
b0 = x[1]
b1 = x[2]
logs2 = x[3]
r = data[:y] - b0 - b1 * data[:x]
logf = -0.5 * dot(r, r)*exp(-logs2) +
(-0.5 * length(data[:y]) * logs2) +
logpdf(Normal(0.0, 1000.0), b0) +
logpdf(Normal(0.0, 1000.0), b1) +
logpdf(InverseGamma(0.001, 0.001), exp(logs2))
# logf = (-0.5 * length(data[:y]) - 0.001) * logs2 -
# (0.5 * dot(r, r) + 0.001) / exp(logs2) -
# 0.5 * b0^2 / 1000 - 0.5 * b1^2 / 1000
return logf
end
function loglikgrad(x::DenseVector)
logf = loglik(x)
grad = calc.gradient(loglik, x)
return logf, grad
end
## MCMC Simulation with No-U-Turn Sampling
n = 10000
burnin = 2000
sim3 = Chains(n, 3, start = (burnin + 1), names = ["b0", "b1", "s2"])
theta = NUTSVariate([0.0, 0.0, 0.0], loglikgrad)
@time for i in 1:n
sample!(theta, adapt = (i <= burnin))
if i > burnin
sim3[i, :, 1] = [theta[1:2]; exp(theta[3])]
end
end
describe(sim3)
sleep(0.2)
draw(plot(sim3), width=10inch, height=4inch, nrow=1, ncol=2, ask=false);
0.851075 seconds (9.85 M allocations: 622.102 MiB, 10.06% gc time) Iterations = 2001:10000 Thinning interval = 1 Chains = 1 Samples per chain = 8000 Empirical Posterior Estimates: Mean SD Naive SE MCSE ESS b0 1.9782926 0.33809500 0.0037800170 0.0073908965 2092.5857 b1 1.0861015 0.11514201 0.0012873268 0.0014818940 6037.1688 s2 3.5505039 0.93607160 0.0104655987 0.0145592966 4133.6848 Quantiles: 2.5% 25.0% 50.0% 75.0% 97.5% b0 1.31737309 1.7646719 1.9736532 2.1943349 2.6438895 b1 0.85680974 1.0102578 1.0869083 1.1623297 1.3080492 s2 2.14473919 2.8814953 3.4071107 4.0660816 5.7536538
速さに変わりはない。
-0.5 * dot(r, r)*exp(-logs2) +
(-0.5 * length(data[:y]) * logs2) +
の部分を
logpdf(MvNormal(fill(exp(0.5*logs2), length(r))), r) +
に書き換えてみた。
################################################################################
## Linear Regression
## y ~ N(b0 + b1 * x, s2)
## b0, b1 ~ N(0, 1000)
## s2 ~ invgamma(0.001, 0.001)
################################################################################
using Mamba
import Calculus; calc = Calculus
## Data
data = Dict(
:x => X,
:y => Y
)
## Log-transformed Posterior(b0, b1, log(s2)) + Constant and Gradient Vector
function loglik(x::DenseVector)
b0 = x[1]
b1 = x[2]
logs2 = x[3]
r = data[:y] - b0 - b1 * data[:x]
logf = logpdf(MvNormal(fill(exp(0.5*logs2), length(r))), r) +
logpdf(Normal(0.0, 1000.0), b0) +
logpdf(Normal(0.0, 1000.0), b1) +
logpdf(InverseGamma(0.001, 0.001), exp(logs2))
# logf = -0.5 * dot(r, r)*exp(-logs2) +
# (-0.5 * length(data[:y]) * logs2) +
# logpdf(Normal(0.0, 1000.0), b0) +
# logpdf(Normal(0.0, 1000.0), b1) +
# logpdf(InverseGamma(0.001, 0.001), exp(logs2))
# logf = (-0.5 * length(data[:y]) - 0.001) * logs2 -
# (0.5 * dot(r, r) + 0.001) / exp(logs2) -
# 0.5 * b0^2 / 1000 - 0.5 * b1^2 / 1000
return logf
end
function loglikgrad(x::DenseVector)
logf = loglik(x)
grad = calc.gradient(loglik, x)
return logf, grad
end
## MCMC Simulation with No-U-Turn Sampling
n = 10000
burnin = 2000
sim4 = Chains(n, 3, start = (burnin + 1), names = ["b0", "b1", "s2"])
theta = NUTSVariate([0.0, 0.0, 0.0], loglikgrad)
@time for i in 1:n
sample!(theta, adapt = (i <= burnin))
if i > burnin
sim4[i, :, 1] = [theta[1:2]; exp(theta[3])]
end
end
describe(sim4)
sleep(0.2)
draw(plot(sim4), width=10inch, height=4inch, nrow=1, ncol=2, ask=false);
1.174281 seconds (9.91 M allocations: 1.164 GiB, 12.04% gc time) Iterations = 2001:10000 Thinning interval = 1 Chains = 1 Samples per chain = 8000 Empirical Posterior Estimates: Mean SD Naive SE MCSE ESS b0 1.9873032 0.33402085 0.0037344667 0.0070002962 2276.7447 b1 1.0839021 0.11823097 0.0013218624 0.0013226702 7990.2313 s2 3.5711500 0.95351640 0.0106606374 0.0183488432 2700.4675 Quantiles: 2.5% 25.0% 50.0% 75.0% 97.5% b0 1.34426011 1.7639035 1.9828091 2.1965159 2.6784113 b1 0.85017333 1.0078761 1.0847434 1.1601784 1.3211359 s2 2.15787938 2.8964387 3.4167551 4.0710239 5.8846329
MvNormal分布のlogpdfを使ったら少し遅くなった。
上で MvNormal を使っていた
logpdf(MvNormal(fill(exp(0.5*logs2), length(data[:y]))), r) +
の部分を次のように書き直した。
sum(logpdf.(Normal(0.0, exp(0.5*logs2)), r[i]) for i in 1:length(r)) +
################################################################################
## Linear Regression
## y ~ N(b0 + b1 * x, s2)
## b0, b1 ~ N(0, 1000)
## s2 ~ invgamma(0.001, 0.001)
################################################################################
using Mamba
import Calculus; calc = Calculus
## Data
data = Dict(
:x => X,
:y => Y
)
## Log-transformed Posterior(b0, b1, log(s2)) + Constant and Gradient Vector
function loglik(x::DenseVector)
b0 = x[1]
b1 = x[2]
logs2 = x[3]
r = data[:y] - b0 - b1 * data[:x]
logf = sum(logpdf.(Normal(0.0, exp(0.5*logs2)), r[i]) for i in 1:length(r)) +
logpdf(Normal(0.0, 1000.0), b0) +
logpdf(Normal(0.0, 1000.0), b1) +
logpdf(InverseGamma(0.001, 0.001), exp(logs2))
return logf
end
function loglikgrad(x::DenseVector)
logf = loglik(x)
grad = calc.gradient(loglik, x)
return logf, grad
end
## MCMC Simulation with No-U-Turn Sampling
n = 10000
burnin = 2000
sim5 = Chains(n, 3, start = (burnin + 1), names = ["b0", "b1", "s2"])
theta = NUTSVariate([0.0, 0.0, 0.0], loglikgrad)
@time for i in 1:n
sample!(theta, adapt = (i <= burnin))
if i > burnin
sim5[i, :, 1] = [theta[1:2]; exp(theta[3])]
end
end
describe(sim5)
sleep(0.2)
draw(plot(sim5), width=10inch, height=4inch, nrow=1, ncol=2, ask=false);
1.780044 seconds (10.87 M allocations: 714.573 MiB, 7.20% gc time) Iterations = 2001:10000 Thinning interval = 1 Chains = 1 Samples per chain = 8000 Empirical Posterior Estimates: Mean SD Naive SE MCSE ESS b0 1.9752781 0.340927228 0.0038116823 0.0065644924 2697.2492 b1 1.0863133 0.115562093 0.0012920235 0.0012355233 8000.0000 s2 3.5474988 0.929502108 0.0103921495 0.0181222109 2630.7428 Quantiles: 2.5% 25.0% 50.0% 75.0% 97.5% b0 1.2971724 1.7498005 1.9704395 2.2027267 2.6511457 b1 0.8567451 1.0097570 1.0888926 1.1616243 1.3057893 s2 2.1489168 2.8888104 3.3868456 4.0466380 5.7859321
MvNormalによる尤度函数の評価を sum と Normal による評価に変えたら少し遅くなった。
sum(logpdf.(Normal(0.0, exp(0.5*logs2)), r[i]) for i in 1:length(r)) +
の部分を
sum(logpdf.(Normal(0.0, exp(0.5*logs2)), r)) +
に書き直したら、以下のように非常に遅くなった。2秒を切っていたのが、35秒を超えてしまう。
################################################################################
## Linear Regression
## y ~ N(b0 + b1 * x, s2)
## b0, b1 ~ N(0, 1000)
## s2 ~ invgamma(0.001, 0.001)
################################################################################
using Mamba
import Calculus; calc = Calculus
## Data
data = Dict(
:x => X,
:y => Y
)
## Log-transformed Posterior(b0, b1, log(s2)) + Constant and Gradient Vector
function loglik(x::DenseVector)
b0 = x[1]
b1 = x[2]
logs2 = x[3]
r = data[:y] - b0 - b1 * data[:x]
logf = sum(logpdf.(Normal(0.0, exp(0.5*logs2)), r)) +
logpdf(Normal(0.0, 1000.0), b0) +
logpdf(Normal(0.0, 1000.0), b1) +
logpdf(InverseGamma(0.001, 0.001), exp(logs2))
return logf
end
function loglikgrad(x::DenseVector)
logf = loglik(x)
grad = calc.gradient(loglik, x)
return logf, grad
end
## MCMC Simulation with No-U-Turn Sampling
n = 10000
burnin = 2000
sim6 = Chains(n, 3, start = (burnin + 1), names = ["b0", "b1", "s2"])
theta = NUTSVariate([0.0, 0.0, 0.0], loglikgrad)
@time for i in 1:n
sample!(theta, adapt = (i <= burnin))
if i > burnin
sim6[i, :, 1] = [theta[1:2]; exp(theta[3])]
end
end
describe(sim6)
sleep(0.2)
draw(plot(sim6), width=10inch, height=4inch, nrow=1, ncol=2, ask=false);
39.154306 seconds (23.10 M allocations: 1.257 GiB, 0.72% gc time) Iterations = 2001:10000 Thinning interval = 1 Chains = 1 Samples per chain = 8000 Empirical Posterior Estimates: Mean SD Naive SE MCSE ESS b0 1.9701685 0.32089077 0.0035876679 0.0070264348 2085.6643 b1 1.0869933 0.11099434 0.0012409544 0.0013327149 6936.2882 s2 3.5018410 0.92357676 0.0103259021 0.0182594902 2558.4015 Quantiles: 2.5% 25.0% 50.0% 75.0% 97.5% b0 1.33580265 1.7554155 1.9707346 2.1863044 2.5932734 b1 0.86919504 1.0137925 1.0862727 1.1628614 1.2992822 s2 2.13516750 2.8458839 3.3479617 3.9812698 5.6973021
y = Stochastic(1,
(b0, b1, s2) -> UnivariateDistribution[
Normal(b0 + b1*data[:X][i], sqrt(s2)) for i in 1:data[:N]
],
false
),
の部分を
y = Stochastic(1,
(b0, b1, s2) -> MvNormal(b0 + b1*data[:X], sqrt(s2)),
false
),
に変えたら、20秒程度だったのに、10秒程度に速くなった。
これでもまだ遅い。どこに律速段階があるのだろうか?
################################################################################
## Linear Regression
## y ~ N(b0 + b1 * x, s2)
## b0, b1 ~ N(0, 1000)
## s2 ~ invgamma(0.001, 0.001)
################################################################################
using Mamba
## Input Data
data = Dict{Symbol, Any}(
:X => X,
:Y => Y,
:N => length(X),
:b00 => 0.0,
:b10 => 0.0,
:s20 => 1.0,
)
## Model Specification
model = Model(
y = Stochastic(1,
(b0, b1, s2) -> MvNormal(b0 + b1*data[:X], sqrt(s2)),
false
),
b0 = Stochastic(
() -> Normal(0.0, 1000.0),
true
),
b1 = Stochastic(
() -> Normal(0.0, 1000.0),
true
),
s2 = Stochastic(
() -> InverseGamma(0.001, 0.001),
true
),
)
## Number of chains
chains = 1
## Initial Values
inits = [
Dict{Symbol, Any}(
:y => data[:Y],
:b0 => data[:b00],
:b1 => data[:b10],
:s2 => data[:s20],
)
for k in 1:chains
]
## Sampling Scheme Assignment
scheme = [
NUTS([:b0, :b1, :s2])
]
setsamplers!(model, scheme)
## MCMC
iterations = 10000
burnin = 2000
thinning = 1
@time mc2 = mcmc(model, data, inits, iterations, burnin=burnin, thin=thinning, chains=chains)
## Summary of MCMC
print("\nSummary:\n\n")
display(describe(mc2))
## Convergence Diagnostics
#print("Gelman Diagnostic:\n")
#show(gelmandiag(mc2))
## Plot
sleep(0.2)
draw(plot(mc2), width=10inch, height=4inch, nrow=1, ncol=2, ask=false);
MCMC Simulation of 10000 Iterations x 1 Chain... Chain 1: 0% [0:00:58 of 0:00:58 remaining] Chain 1: 10% [0:00:10 of 0:00:11 remaining] Chain 1: 20% [0:00:08 of 0:00:10 remaining] Chain 1: 30% [0:00:07 of 0:00:10 remaining] Chain 1: 40% [0:00:06 of 0:00:10 remaining] Chain 1: 50% [0:00:05 of 0:00:10 remaining] Chain 1: 60% [0:00:04 of 0:00:10 remaining] Chain 1: 70% [0:00:03 of 0:00:10 remaining] Chain 1: 80% [0:00:02 of 0:00:10 remaining] Chain 1: 90% [0:00:01 of 0:00:10 remaining]
nothing
Chain 1: 100% [0:00:00 of 0:00:10 remaining] 10.344650 seconds (41.88 M allocations: 1.874 GiB, 4.71% gc time) Summary: Iterations = 2001:10000 Thinning interval = 1 Chains = 1 Samples per chain = 8000 Empirical Posterior Estimates: Mean SD Naive SE MCSE ESS s2 3.7809102 1.022606496 0.0114330882 0.0191838210 2841.4948 b1 1.0867164 0.119475718 0.0013357791 0.0014421200 6863.6713 b0 1.9703375 0.338640133 0.0037861118 0.0064969599 2716.7926 Quantiles: 2.5% 25.0% 50.0% 75.0% 97.5% s2 2.26738242 3.0553677 3.6205297 4.3143896 6.2202623 b1 0.85465384 1.0045591 1.0872332 1.1662116 1.3226845 b0 1.30661956 1.7427829 1.9716513 2.1955837 2.6265877
println("LSE:")
println("b0_lse = $b0_lse")
println("b1_lse = $b1_lse")
println("s2_lse = $s2_lse")
function printresult(mc)
for n in ["b0", "b1", "s2"]
k = findfirst(x -> x == n, mc.names)
m = mean(mc.value[:,k])
println("mean($n) = $m")
end
end
println()
println("mc:")
printresult(mc)
println()
println("sim2:")
printresult(sim2)
println()
println("sim3:")
printresult(sim3)
println()
println("sim4:")
printresult(sim4)
println()
println("sim5:")
printresult(sim5)
println()
println("sim6:")
printresult(sim6)
println()
println("mc2:")
printresult(mc2)
;
LSE: b0_lse = 1.9733544847267817 b1_lse = 1.0864307623437985 s2_lse = 3.334179811356694 mc: mean(b0) = 1.964605798870186 mean(b1) = 1.0848856816760002 mean(s2) = 3.814595473903992 sim2: mean(b0) = 1.965375812441522 mean(b1) = 1.0862501411085563 mean(s2) = 3.8004729628230827 sim3: mean(b0) = 1.978292560107559 mean(b1) = 1.0861015202379718 mean(s2) = 3.5505039327125942 sim4: mean(b0) = 1.9873032421548582 mean(b1) = 1.0839020582953462 mean(s2) = 3.5711500222303 sim5: mean(b0) = 1.9752781070252512 mean(b1) = 1.0863132980457235 mean(s2) = 3.5474988198089004 sim6: mean(b0) = 1.970168478479773 mean(b1) = 1.0869932970920726 mean(s2) = 3.501840973177938 mc2: mean(b0) = 1.9703375019120009 mean(b1) = 1.0867164406924 mean(s2) = 3.7809101839656822