In [1]:
import stan
import pandas as pd
import numpy as np
import arviz as az
import matplotlib.pyplot as plt
import nest_asyncio
nest_asyncio.apply() # To use PyStan3 on Notebook
data = [ #City Temp. Lat.
['Sapporo', 15, 43.1],
['Sendai', 20, 38.3],
['Niigata', 17, 37.9],
['Kanazawa', 18, 36.6],
['Tokyo', 24, 35.7],
['Osaka', 22, 34.7],
['Fukuoka', 22, 33.6],
['Kochi', 24, 33.6],
['Kagoshima', 23, 31.6],
['Naha', 26, 26.2]
]
data = np.array(data)
ID = data.T[0]
Y = np.array(data.T[1], dtype=float)
X = np.array(data.T[2], dtype=float)
N = len(Y)
for v in zip(ID, Y, X):
print(v)
f = open('Results.txt', 'w') # Output text file
mean_X = np.mean(X)
sd_X = np.std(X)
print(f'X: mean = {mean_X:.2f} sd = {sd_X:.3f}')
f.write(f'X: mean = {mean_X:.2f} sd = {sd_X:.3f}\n')
X = (X - mean_X) / sd_X
mean_Y = np.mean(Y)
sd_Y = np.std(Y)
print(f'Y: mean = {mean_Y:.2f} sd = {sd_Y:.3f}')
f.write(f'\nY: mean = {mean_Y:.2f} sd = {sd_Y:.3f}\n')
Y = (Y - mean_Y) / sd_Y
Data = {'N': N, 'Y': Y, 'X': X}
stan_code = """
data{
int N;
vector[N] Y;
vector[N] X;
}
parameters{
real b0;
real b1;
real<lower = 0> sgm;
}
transformed parameters{
vector[N] mu;
mu = b0 + b1 * X;
}
model{
b0 ~ normal(0.0, 2.5);
b1 ~ normal(0.0, 2.5);
sgm ~ exponential(1/1.0);
Y ~ normal(mu, sgm);
}
"""
sm = stan.build(stan_code, data = Data)
fit = sm.sample()
d_frm = fit.to_frame()
i_data = az.from_pystan(posterior = fit, posterior_model = sm)
print(az.summary(i_data))
f.write('\n')
f.write(az.summary(i_data).__str__())
def med_mad_sd(x):
""" Calculation of median and mad_sd
Gelman et al. (2021), p.73
"""
med = np.median(x)
mad = np.median(np.abs(x - med))
mad_sd = 1.483 * mad
return med, mad_sd
med_b0, mad_sd_b0 = med_mad_sd(d_frm['b0'])
med_b1, mad_sd_b1 = med_mad_sd(d_frm['b1'])
med_sgm, mad_sd_sgm = med_mad_sd(d_frm['sgm'])
f.write('\n\n')
f.write(f'b0: med.= {med_b0:.3f}, mad sd = {mad_sd_b0:.3f}\n')
f.write(f'\nb1: med.= {med_b1:.3f}, mad sd = {mad_sd_b1:.3f}\n')
f.write(f'\nsgm: med.= {med_sgm:.3f}, mad sd = {mad_sd_sgm:.3f}\n')
az.plot_kde(d_frm['b0'])
plt.xlabel('b0', fontsize = 16)
plt.title(f'b0: med. = {med_b0:.3f}, mad sd = {mad_sd_b0:.3f}', fontsize = 16)
plt.show()
az.plot_kde(d_frm['b1'])
plt.xlabel('b1', fontsize = 16)
plt.title(f'b1: med. = {med_b1:.3f}, mad sd = {mad_sd_b1:.3f}', fontsize = 16)
plt.show()
az.plot_kde(d_frm['sgm'])
plt.xlabel('$\sigma$', fontsize = 14)
plt.title(f'$\sigma$: med. = {med_sgm:.3f}, mad sd = {mad_sd_sgm:.3f}',
fontsize = 16)
plt.show()
plt.plot(d_frm['b0'], d_frm['b1'], 'o', alpha = 0.3)
plt.xlabel('b0', fontsize = 16)
plt.ylabel('b1', fontsize = 16)
plt.tight_layout()
plt.show()
#
# Regression line and Scatter plot
#
for i in range(len(ID)):
plt.plot(X[i], Y[i], 'o', c='b')
plt.text(X[i], Y[i], ID[i])
min_x = np.min(X)
max_x = np.max(X)
plt.plot([min_x, max_x], [med_b0 + med_b1*min_x, med_b0 + med_b1*max_x],
label = f'Y = {med_b0:.3f} + {med_b1:.3f} * X')
plt.xlabel('X', fontsize = 14)
plt.ylabel('Y', fontsize = 14)
plt.legend()
plt.tight_layout()
plt.show()
f.close()
print('Results.txt was saved.')
Building: found in cache, done. Sampling: 0%
('Sapporo', 15.0, 43.1)
('Sendai', 20.0, 38.3)
('Niigata', 17.0, 37.9)
('Kanazawa', 18.0, 36.6)
('Tokyo', 24.0, 35.7)
('Osaka', 22.0, 34.7)
('Fukuoka', 22.0, 33.6)
('Kochi', 24.0, 33.6)
('Kagoshima', 23.0, 31.6)
('Naha', 26.0, 26.2)
X: mean = 35.13 sd = 4.252
Y: mean = 21.10 sd = 3.330
Sampling: 25% (2000/8000) Sampling: 50% (4000/8000) Sampling: 75% (6000/8000) Sampling: 100% (8000/8000) Sampling: 100% (8000/8000), done. Messages received during sampling: Gradient evaluation took 2.2e-05 seconds 1000 transitions using 10 leapfrog steps per transition would take 0.22 seconds. Adjust your expectations accordingly! Gradient evaluation took 2.5e-05 seconds 1000 transitions using 10 leapfrog steps per transition would take 0.25 seconds. Adjust your expectations accordingly! Gradient evaluation took 2.3e-05 seconds 1000 transitions using 10 leapfrog steps per transition would take 0.23 seconds. Adjust your expectations accordingly! Gradient evaluation took 2.1e-05 seconds 1000 transitions using 10 leapfrog steps per transition would take 0.21 seconds. Adjust your expectations accordingly!
mean sd hdi_3% hdi_97% mcse_mean mcse_sd ess_bulk ess_tail \
b0 0.003 0.214 -0.425 0.391 0.004 0.004 2707.0 1877.0
b1 -0.858 0.205 -1.243 -0.465 0.004 0.003 2838.0 2230.0
sgm 0.633 0.189 0.346 0.961 0.004 0.003 2108.0 1951.0
mu[0] -1.606 0.438 -2.458 -0.784 0.009 0.006 2699.0 2188.0
mu[1] -0.637 0.262 -1.152 -0.172 0.005 0.004 2659.0 1977.0
mu[2] -0.556 0.251 -1.017 -0.078 0.005 0.004 2656.0 2006.0
mu[3] -0.294 0.225 -0.732 0.115 0.005 0.003 2671.0 2038.0
mu[4] -0.112 0.215 -0.546 0.279 0.004 0.004 2690.0 1933.0
mu[5] 0.090 0.215 -0.336 0.478 0.004 0.004 2720.0 1887.0
mu[6] 0.312 0.226 -0.129 0.725 0.004 0.004 2767.0 2028.0
mu[7] 0.312 0.226 -0.129 0.725 0.004 0.004 2767.0 2028.0
mu[8] 0.715 0.274 0.205 1.246 0.005 0.004 2828.0 2064.0
mu[9] 1.805 0.482 0.878 2.708 0.009 0.007 2882.0 1954.0
r_hat
b0 1.0
b1 1.0
sgm 1.0
mu[0] 1.0
mu[1] 1.0
mu[2] 1.0
mu[3] 1.0
mu[4] 1.0
mu[5] 1.0
mu[6] 1.0
mu[7] 1.0
mu[8] 1.0
mu[9] 1.0
Results.txt was saved.