In [1]:
import stan
import matplotlib.pyplot as plt
import numpy as np
import arviz as az
import seaborn as sb
import nest_asyncio
nest_asyncio.apply()
RawData = [[1, 1, 0, 0, 1, 0, 0],
[1, 1, 0, 0, 0, 1, 0],
[10, 1, 0, 0, 0, 0, 1],
[3, 0, 1, 0, 1, 0, 0],
[2, 0, 1, 0, 0, 1, 0],
[10, 0, 1, 0, 0, 0, 1],
[6, 0, 0, 1, 1, 0, 0],
[2, 0, 0, 1, 0, 1, 0],
[4, 0, 0, 1, 0, 0, 1]]
F = []
Xr1 = []
Xr2 = []
Xr3 = []
Xc1 = []
Xc2 = []
Xc3 = []
for v in RawData:
F.append(v[0])
Xr1.append(v[1])
Xr2.append(v[2])
Xr3.append(v[3])
Xc1.append(v[4])
Xc2.append(v[5])
Xc3.append(v[6])
for i in range(9):
print(' ', F[i], ' ', Xr1[i], ' ', Xr2[i], ' ', Xr3[i],
' ', Xc1[i], ' ', Xc2[i], ' ', Xc3[i])
N = 9
Data = {'N': N, 'F': F, 'Xr1': Xr1, 'Xr2': Xr2, 'Xr3': Xr3,
'Xc1': Xc1, 'Xc2': Xc2, 'Xc3': Xc3}
stan_code = """
data{
int N;
array[N] int F;
array[N] int Xr1;
array[N] int Xr2;
array[N] int Xr3;
array[N] int Xc1;
array[N] int Xc2;
array[N] int Xc3;
}
parameters{
real mu;
real a2;
real a3;
real b2;
real b3;
real g22;
real g23;
real g32;
real g33;
}
transformed parameters{
array[N] real Lmbd;
array[N] real LogLmbd;
for (i in 1:N)
LogLmbd[i] = mu + a2 * Xr2[i] + a3 * Xr3[i]
+ b2 * Xc2[i] + b3 * Xc3[i]
+ g22 * Xr2[i] * Xc2[i] + g23 * Xr2[i] * Xc3[i]
+ g32 * Xr3[i] * Xc2[i] + g33 * Xr3[i] * Xc3[i];
for (i in 1:N)
Lmbd[i] = exp(LogLmbd[i]);
}
model{
mu ~ normal(0.0, 10.0);
a2 ~ normal(0.0, 10.0);
a3 ~ normal(0.0, 10.0);
b2 ~ normal(0.0, 10.0);
b3 ~ normal(0.0, 10.0);
g22 ~ normal(0.0, 10.0);
g23 ~ normal(0.0, 10.0);
g32 ~ normal(0.0, 10.0);
g33 ~ normal(0.0, 10.0);
for (i in 1:N)
F[i] ~ poisson(Lmbd[i]);
}
"""
sm = stan.build(stan_code, data = Data)
fit = sm.sample()
f = open('Results.txt', 'w')
i_data = az.from_pystan(posterior = fit, posterior_model = sm)
smry = az.summary(i_data)
print('Summary =\n',smry)
f.write('\nSummary...\n')
f.write(smry.__str__())
f.write('\n')
az.plot_trace(i_data, var_names = ('a2', 'a3', 'b2', 'b3'))
plt.show()
d_fm = fit.to_frame()
a3 = d_fm['a3']
cnt_a3_pos = (a3 > 0.0).sum()
print('\nP(a3 > 0.0) = {0:.3f}'.format(cnt_a3_pos / len(a3)))
f.write('\nP(a3 > 0.0) = {0:.3f}\n'.format(cnt_a3_pos / len(a3)))
b3 = d_fm['b3']
cnt_b3_pos = (b3 > 0.0).sum()
print('P(b3 > 0.0) = {0:.3f}'.format(cnt_b3_pos / len(b3)))
f.write('\nP(b3 > 0.0) = {0:.3f}\n'.format(cnt_b3_pos / len(b3)))
az.plot_kde(d_fm['a2'], label = 'a2')
az.plot_kde(d_fm['a3'], label = 'a3')
az.plot_kde(d_fm['b2'], label = 'b2')
az.plot_kde(d_fm['b3'], label = 'b3')
plt.legend()
plt.show()
g33 = d_fm['g33']
print(len(g33))
print(np.shape(g33))
cnt_g33_neg = (g33 < 0.0).sum()
print('P(g33 < 0.0) = {0:.3f}'.format(cnt_g33_neg / len(g33)))
f.write('\nP(g33 < 0.0) = {0:.3f}\n'.format(cnt_g33_neg / len(g33)))
az.plot_kde(d_fm['g22'], label = 'g22')
az.plot_kde(d_fm['g23'], label = 'g23')
az.plot_kde(d_fm['g32'], label = 'g32')
az.plot_kde(d_fm['g33'], label = 'g33')
plt.legend()
plt.show()
lmbd_hat = np.empty(N)
for i in range(N):
lmbd_hat[i] = np.median(d_fm[f'Lmbd.{i+1}'])
plt.plot(lmbd_hat[i], F[i], 'o', c='b')
min_lmbd_hat = np.min(lmbd_hat)
max_lmbd_hat = np.max(lmbd_hat)
plt.plot([min_lmbd_hat, max_lmbd_hat], [min_lmbd_hat, max_lmbd_hat],
label = '$F = \hat{\lambda}$')
plt.xlabel('$\hat{\lambda}$', fontsize = 16)
plt.ylabel('F', fontsize = 14)
plt.legend(fontsize = 16)
plt.tight_layout()
plt.show()
f.close()
print('Results.txt was saved.')
Building: found in cache, done. Sampling: 0%
1 1 0 0 1 0 0 1 1 0 0 0 1 0 10 1 0 0 0 0 1 3 0 1 0 1 0 0 2 0 1 0 0 1 0 10 0 1 0 0 0 1 6 0 0 1 1 0 0 2 0 0 1 0 1 0 4 0 0 1 0 0 1
Sampling: 8% (600/8000) Sampling: 32% (2600/8000) Sampling: 58% (4600/8000) Sampling: 82% (6600/8000) Sampling: 100% (8000/8000) Sampling: 100% (8000/8000), done. Messages received during sampling: Gradient evaluation took 3.9e-05 seconds 1000 transitions using 10 leapfrog steps per transition would take 0.39 seconds. Adjust your expectations accordingly! Gradient evaluation took 4.4e-05 seconds 1000 transitions using 10 leapfrog steps per transition would take 0.44 seconds. Adjust your expectations accordingly! Gradient evaluation took 4.3e-05 seconds 1000 transitions using 10 leapfrog steps per transition would take 0.43 seconds. Adjust your expectations accordingly! Gradient evaluation took 4.6e-05 seconds 1000 transitions using 10 leapfrog steps per transition would take 0.46 seconds. Adjust your expectations accordingly!
Summary =
mean sd hdi_3% hdi_97% mcse_mean mcse_sd ess_bulk \
mu -0.289 1.075 -2.312 1.470 0.043 0.039 763.0
a2 1.186 1.256 -1.113 3.532 0.046 0.041 805.0
a3 1.986 1.152 -0.046 4.110 0.044 0.036 796.0
b2 -0.324 1.608 -3.183 2.820 0.059 0.041 770.0
b3 2.530 1.120 0.656 4.637 0.043 0.035 769.0
g22 -0.138 1.896 -3.576 3.491 0.063 0.045 922.0
g23 -1.177 1.333 -3.638 1.288 0.048 0.042 831.0
g32 -0.947 1.811 -4.405 2.382 0.062 0.050 849.0
g33 -2.960 1.303 -5.376 -0.596 0.046 0.036 888.0
Lmbd[0] 1.172 1.062 0.002 3.072 0.032 0.023 763.0
Lmbd[1] 0.953 0.962 0.001 2.769 0.021 0.015 1440.0
Lmbd[2] 9.890 3.102 4.205 15.593 0.045 0.032 4547.0
Lmbd[3] 2.950 1.743 0.225 6.076 0.031 0.022 2812.0
Lmbd[4] 2.018 1.412 0.092 4.596 0.021 0.016 4531.0
Lmbd[5] 9.991 3.166 4.275 15.814 0.047 0.034 4605.0
Lmbd[6] 5.963 2.499 1.617 10.591 0.038 0.027 4317.0
Lmbd[7] 1.981 1.376 0.084 4.476 0.020 0.015 4553.0
Lmbd[8] 4.045 1.982 0.841 7.856 0.029 0.022 4528.0
LogLmbd[0] -0.289 1.075 -2.312 1.470 0.043 0.039 763.0
LogLmbd[1] -0.613 1.231 -2.959 1.355 0.035 0.027 1440.0
LogLmbd[2] 2.241 0.322 1.651 2.857 0.005 0.003 4547.0
LogLmbd[3] 0.898 0.644 -0.355 2.029 0.013 0.009 2812.0
LogLmbd[4] 0.435 0.803 -1.023 1.833 0.012 0.010 4531.0
LogLmbd[5] 2.251 0.323 1.635 2.841 0.005 0.003 4605.0
LogLmbd[6] 1.697 0.432 0.924 2.549 0.007 0.005 4317.0
LogLmbd[7] 0.425 0.779 -1.060 1.738 0.012 0.010 4553.0
LogLmbd[8] 1.267 0.539 0.259 2.200 0.008 0.006 4528.0
ess_tail r_hat
mu 717.0 1.01
a2 756.0 1.01
a3 727.0 1.01
b2 731.0 1.00
b3 721.0 1.01
g22 1094.0 1.01
g23 808.0 1.01
g32 1009.0 1.00
g33 920.0 1.00
Lmbd[0] 717.0 1.01
Lmbd[1] 1389.0 1.00
Lmbd[2] 3376.0 1.00
Lmbd[3] 1962.0 1.00
Lmbd[4] 3529.0 1.00
Lmbd[5] 3605.0 1.00
Lmbd[6] 3332.0 1.00
Lmbd[7] 3414.0 1.00
Lmbd[8] 3326.0 1.00
LogLmbd[0] 717.0 1.01
LogLmbd[1] 1389.0 1.00
LogLmbd[2] 3376.0 1.00
LogLmbd[3] 1962.0 1.00
LogLmbd[4] 3529.0 1.00
LogLmbd[5] 3605.0 1.00
LogLmbd[6] 3332.0 1.00
LogLmbd[7] 3414.0 1.00
LogLmbd[8] 3326.0 1.00
P(a3 > 0.0) = 0.978 P(b3 > 0.0) = 0.998
4000 (4000,) P(g33 < 0.0) = 0.997
Results.txt was saved.