I hereby claim:
- I am pierredc on github.
- I am pierredc (https://keybase.io/pierredc) on keybase.
- I have a public key ASB9VnZNw9mP3CqUh3-Ck40NjzqcsaWezeVpNldjh57BLAo
To claim this, I am signing this object:
I hereby claim:
To claim this, I am signing this object:
x = np.array([ | |
[1,29,57,85,113,141,169,197,225,253,1,29,57,85,113,141,169,197,225,253], | |
[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], | |
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1] ]) | |
counts = np.array([ | |
[966,675,542,406,362,355,296,284,238,244,960,503,322,342,291,252,245,267,256,250], | |
[943,608,487,404,310,276,290,248,250,248,951,623,429,337,279,248,248,251,268,252], | |
[971,687,550,391,372,293,278,264,259,260,940,596,424,320,285,251,249,231,253,253], | |
[961,660,501,416,316,294,299,268,253,239,966,632,413,339,269,265,276,241,254,262], |
from scipy.optimize import curve_fit | |
def gsp(x, a, alpha, alpha_c, b): | |
return x[1] * (a * alpha ** x[0] + b) \ | |
+ x[2] * (a * (alpha * alpha_c) ** x[0] + b) |
# curve fit | |
popt, pcov = curve_fit(f = gsp, xdata = x, | |
ydata = np.mean(counts, axis = 0) / shots, | |
sigma = np.std(counts, axis = 0) / shots / np.sqrt(num_samples), | |
absolute_sigma = True, | |
bounds = ([scale - .1, .9, .9, (1. - scale) - .1], | |
[scale + .1, 1., 1., (1. - scale) + .1])) | |
# statistical error on pararameters | |
perr = np.sqrt(np.diag(pcov)) |
import pymc3 as pm | |
h_model = pm.Model() |
with h_model: | |
BoundedUniform = pm.Bound(pm.Uniform, | |
lower=np.fmax(popt-0.1, np.full(popt.shape, 1.e-9)), | |
upper=np.fmin(popt+0.1, np.full(popt.shape, 1.-1.e-9))) | |
pi = BoundedUniform("π", testval = popt, shape = popt.shape) | |
EPC = pm.Deterministic('EPC', scale * (1 - pi[2]) ) |
with h_model: | |
GSP = pi[0] * ( x[1] * pi[1] ** x[0] + \ | |
x[2] * (pi[1] * pi[2] ) ** x[0] ) + pi[3] |
with h_model: | |
BoundedGamma = pm.Bound(pm.Gamma, | |
lower= 0.0005, | |
upper= 0.004 + 0.001 * (N - 1) ) | |
sigma = BoundedGamma("σ_Beta", | |
alpha = 10 - 5 * (N - 1), | |
beta = 10000 - 8000 * (N - 1), | |
testval = 0.001 + 0.0015 * (N - 1), | |
shape = counts.shape[1]) |
with h_model: | |
theta = pm.Beta('θ', mu = GSP, sigma = sigma, | |
shape = counts.shape[1]) |
with h_model: | |
p = pm.Binomial("Counts", p = theta, observed = counts, n = shots) |