$G=\frac{1}{2}\pi^2 \left[R_1^2+R_2^2+l^2 -\left\{\left(R_1^2+R_2^2+l^2\right)^2-4R_1^2R_2\right\}^{\frac{1}{2}} \right]$
$G \ge \frac{A_1A_2}{R_1^2+R_2^2+l^2}$
In [2]:
from pprint import pprint
import numpy as np
import pymc3 as pm
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
sns.set(font_scale=1.5)
sns.set_context("notebook", rc={"lines.linewidth": 3})
%matplotlib inline
In [4]:
def getBoundedNormal_dist(mean=None, FWHM=None, name=None, lower=0, upper=1e6):
"""
Make a bounded normal distribution
NOTE: https://github.com/pymc-devs/pymc3/issues/1672 bounded dist fail until 3.1 on
non array bounds!!
"""
assert mean is not None
assert FWHM is not None
assert name is not None
BoundedNormal = mc3.Bound(mc3.Normal, lower=lower, upper=upper)
return BoundedNormal('{0}'.format(name), mu=mean, sd=FWHMtoSD_Normal(mean * (FWHM / 100.)))
def Sullivan_Bound(R1, R2, l):
A1 = np.pi*R1**2
A2 = np.pi*R2**2
top = A1*A2
bottom = R1**2+R2**2+l**2
return top/bottom
def Sullivan(R1, R2, l):
f = 0.5*np.pi**2
t1 = R1**2+R2**2+l**2
t2 = 4*R1**2*R2**2
G = f*(t1 - (t1**2-t2)**0.5 )
return G
def frac_bounds(trace):
med = np.median(trace)
bounds = np.percentile(trace, (2.5, 97.5))
frac = (med-bounds[0])/med
return med, frac*100
In [5]:
Sullivan(4, 5, 10)
Out[5]:
In [28]:
R1 = 0.5
R2 = 0.5
l = np.linspace(2, 20, 100)
fig, ax = plt.subplots(ncols=2, figsize=(8,4))
ax[0].plot(l, Sullivan(R1, R2, l), lw=2)
ax[0].set_xlabel('Distance between [$cm$]')
ax[0].set_ylabel('GF [$cm^2sr$]');
ax[1].loglog(l, Sullivan(R1, R2, l), lw=2)
ax[1].set_xlabel('Distance between [$cm$]');
ax[0].grid(True, which='both')
ax[1.grid(True, which='both')
In [33]:
with pm.Model() as model1:
T1 = pm.Normal('T1', 1.0, 0.1e-2) # 1cm +/- 0.1mm
T2 = pm.Normal('T2', 1.0, 0.1e-2) # 1cm +/- 0.1mm
T3 = pm.Normal('T3', 1.0, 0.1e-2) # 1cm +/- 0.1mm
T4 = pm.Normal('T4', 1.0, 0.1e-2) # 1cm +/- 0.1mm
T5 = pm.Normal('T5', 1.0, 0.1e-2) # 1cm +/- 0.1mm
R1 = 0.5
R2 = 0.5
R3 = 0.5
G = pm.Deterministic('G', Sullivan(R1, R3, T1+T2+T3+T4+T5))
Gbound = pm.Deterministic('Gbound', Sullivan_Bound(R1, R3, T1+T2+T3+T4+T5))
trace = pm.sample(1000, chains=4, target_accept=0.9)
In [34]:
pm.summary(trace).round(3)
Out[34]:
In [35]:
pm.traceplot(trace, combined=False);
In [36]:
gf = frac_bounds(trace['G'])
gbf = frac_bounds(trace['Gbound'])
print("G={:.5f} +/- {:.2f}%".format(gf[0], gf[1]))
print("Gbound={:.5f} +/- {:.2f}%".format(gbf[0], gbf[1]))
In [9]:
pm = np.logspace(-3, -1, 10)
ans = {}
for ii, p in enumerate(pm):
print(p, ii+1, len(pm))
with mc3.Model() as model2:
T1 = mc3.Normal('T1', 1.0, tau=(p)**-2) # 1cm +/- 0.1mm
T2 = mc3.Normal('T2', 1.0, tau=(p)**-2) # 1cm +/- 0.1mm
T3 = mc3.Normal('T3', 1.0, tau=(p)**-2) # 1cm +/- 0.1mm
T4 = mc3.Normal('T4', 1.0, tau=(p)**-2) # 1cm +/- 0.1mm
T5 = mc3.Normal('T5', 1.0, tau=(p)**-2) # 1cm +/- 0.1mm
R1 = 0.5
R2 = 0.5
R3 = 0.5
G = mc3.Deterministic('G', Sullivan(R1, R3, T1+T2+T3+T4+T5))
Gbound = mc3.Deterministic('Gbound', Sullivan_Bound(R1, R3, T1+T2+T3+T4+T5))
start = mc3.find_MAP()
trace = mc3.sample(10000, start=start, jobs=2)
ans[p] = gf = frac_bounds(trace['G'])
In [10]:
pprint(ans)
In [11]:
vals = np.asarray(list(ans.keys()))
gs = np.asarray([ans[v][0] for v in ans ])
gse = np.asarray([ans[v][1] for v in ans ])
In [12]:
valsf = (vals/1.0)*100
plt.errorbar(valsf, gs, yerr=gse, elinewidth=1, capsize=2, barsabove=True)
plt.ylim([0,15])
# plt.xscale('log')
Out[12]:
In [13]:
plt.plot(valsf, gse)
Out[13]:
In [ ]: