In [1]:
import numpy as np
import scipy.stats as ss
import matplotlib.pyplot as plt
import stan
import arviz as az
import nest_asyncio
nest_asyncio.apply()
x_names = ['Lat.', 'W_Long.']
data = [#City Temp. Lat. W_Long.
['Boston', 14, 42.4, 71.1],
['Washington', 18, 38.9, 77],
['Miami', 33, 25.8, 80.2],
['Detroit', 13, 42.3, 83],
['Atlanta', 22, 33.7, 84.4],
['Chicago', 15, 41.9, 87.6],
['Houston', 32, 29.8, 95.4],
['Oklahoma_City',21, 35.5, 97.5],
['Denver', 16, 39.7, 105],
['Los_Angeles', 23, 34.1, 118.2],
['San_Francisco',19, 37.8, 122.4],
['Seattle', 23, 47.6, 122.3]
]
f = open('Results.txt', 'w') # Output text file
data = np.array(data)
case_id = data.T[0]
y = np.array(data.T[1], dtype=float)
X = np.array(data[:,2:4], dtype=float)
print(data)
print(case_id)
print(y)
print(X)
y_mean = np.mean(y)
y_sd = np.std(y)
f.write(f'\n\ny_mean = {y_mean:.3f}, y_sd = {y_sd:.3f}\n')
X = np.array(X, dtype='float')
x_means = np.mean(X, axis=0)
print(x_means)
x_sds = np.std(X, axis=0)
print(x_sds)
for j in range(len(x_names)):
f.write(f'\nb{j+1}({x_names[j]}): ')
f.write(f'mean = {x_means[j]:.3f}, sd = {x_sds[j]:.3f}\n')
y = ss.zscore(np.array(y, dtype ='float'))
X = ss.zscore(np.array(X, dtype = 'float'))
print(y)
print(X)
N, M = np.shape(X)
print(N, M)
stan_code = """
data {
int N;
int M;
vector[N] y;
array[N] vector[M] X;
}
parameters {
real b0;
vector[M] b;
real<lower = 0.0> vsgm;
}
transformed parameters {
real sgm;
sgm = vsgm + 0.001;
}
model {
b0 ~ normal(0.0, 2.5);
b ~ normal(0.0, 2.5);
vsgm ~ exponential(1/1.0);
for (i in 1:N) {
y[i] ~ normal(b0 + dot_product(b,X[i]), sgm);
}
}
"""
sm = stan.build(stan_code, data = {'N':N, 'M':M, 'y':y, 'X':X})
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\n')
f.write(az.summary(i_data).__str__())
az.plot_trace(i_data)
plt.tight_layout()
plt.show()
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
b0_med, b0_mad_sd = med_mad_sd(d_frm['b0'])
print(f'\nb0_med = {b0_med:.3f}, b0_mad_sd = {b0_mad_sd:.3f}')
f.write(f'\n\nb0_med = {b0_med:.3f}, b0_mad_sd = {b0_mad_sd:.3f}\n')
az.plot_kde(d_frm['b0'])
plt.title('Posterior distribution for b0' +
f'\nMed. = {b0_med:.3f}, mad sd = {b0_mad_sd:.3f}')
plt.show()
b_med = np.empty(M)
b_mad_sd = np.empty(M)
for j in range(M):
b_med[j], b_mad_sd[j] = med_mad_sd(d_frm[f'b.{j+1}'])
print(x_names[j])
print(b_med[j], b_mad_sd[j])
f.write(f'\nb{j+1}({x_names[j]}):' +
f' med = {b_med[j]:.3f}, mad_sd = {b_mad_sd[j]:.3f}\n')
az.plot_kde(d_frm[f'b.{j+1}'])
plt.title(f'Posterir distribution for b{j+1}: ' + x_names[j] +
f'\nMed. = {b_med[j]:.3f}, mad sd = {b_mad_sd[j]:.3f}')
plt.show()
y_hat = np.empty(N)
for i in range(N):
y_hat[i] = b0_med + np.dot(b_med, X[i])
for i in range(N):
plt.plot(y_hat[i], y[i], 'o', c='b')
plt.text(y_hat[i], y[i], case_id[i])
plt.xlabel('$\hat{y}$(model)')
plt.ylabel('$y$(data)')
x_min = np.min(y_hat)
x_max = np.max(y_hat)
plt.plot([x_min, x_max], [x_min, x_max], label ='y(data) = $\hat{y}$(model)')
plt.legend()
plt.show()
print('Results.txt was saved.')
Building: found in cache, done. Sampling: 0%
[['Boston' '14' '42.4' '71.1'] ['Washington' '18' '38.9' '77'] ['Miami' '33' '25.8' '80.2'] ['Detroit' '13' '42.3' '83'] ['Atlanta' '22' '33.7' '84.4'] ['Chicago' '15' '41.9' '87.6'] ['Houston' '32' '29.8' '95.4'] ['Oklahoma_City' '21' '35.5' '97.5'] ['Denver' '16' '39.7' '105'] ['Los_Angeles' '23' '34.1' '118.2'] ['San_Francisco' '19' '37.8' '122.4'] ['Seattle' '23' '47.6' '122.3']] ['Boston' 'Washington' 'Miami' 'Detroit' 'Atlanta' 'Chicago' 'Houston' 'Oklahoma_City' 'Denver' 'Los_Angeles' 'San_Francisco' 'Seattle'] [14. 18. 33. 13. 22. 15. 32. 21. 16. 23. 19. 23.] [[ 42.4 71.1] [ 38.9 77. ] [ 25.8 80.2] [ 42.3 83. ] [ 33.7 84.4] [ 41.9 87.6] [ 29.8 95.4] [ 35.5 97.5] [ 39.7 105. ] [ 34.1 118.2] [ 37.8 122.4] [ 47.6 122.3]] [37.45833333 95.34166667] [ 5.78409577 17.23530574] [-1.08992722 -0.44404442 1.97801606 -1.25139791 0.20183837 -0.92845652 1.81654536 0.04036767 -0.76698582 0.36330907 -0.28257372 0.36330907] [[ 0.85435423 -1.40651213] [ 0.24924668 -1.06419155] [-2.01558442 -0.87852614] [ 0.83706544 -0.71606891] [-0.64977025 -0.6348403 ] [ 0.76791029 -0.4491749 ] [-1.32403294 0.00338453] [-0.33857208 0.12522745] [ 0.38755698 0.56038074] [-0.5806151 1.32625053] [ 0.05907002 1.56993637] [ 1.75337115 1.56413433]] 12 2
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.7e-05 seconds 1000 transitions using 10 leapfrog steps per transition would take 0.27 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.7e-05 seconds 1000 transitions using 10 leapfrog steps per transition would take 0.27 seconds. Adjust your expectations accordingly! Gradient evaluation took 2.7e-05 seconds 1000 transitions using 10 leapfrog steps per transition would take 0.27 seconds. Adjust your expectations accordingly!
mean sd hdi_3% hdi_97% mcse_mean mcse_sd ess_bulk ess_tail \
b0 -0.004 0.225 -0.443 0.393 0.004 0.004 3615.0 2420.0
b[0] -0.814 0.225 -1.247 -0.391 0.004 0.003 3296.0 2504.0
b[1] 0.265 0.228 -0.173 0.697 0.004 0.003 3516.0 2428.0
vsgm 0.745 0.198 0.427 1.104 0.004 0.003 2334.0 2398.0
sgm 0.746 0.198 0.428 1.105 0.004 0.003 2334.0 2398.0
r_hat
b0 1.0
b[0] 1.0
b[1] 1.0
vsgm 1.0
sgm 1.0
b0_med = 0.001, b0_mad_sd = 0.210
Lat. -0.8117544851498923 0.2088277078456016
W_Long. 0.26429376051216646 0.20934461760597944
Results.txt was saved.