<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Shots are fired isotropically from a point and hit a position sensitive detector
There is no scattering
y is fixed to be 1 away
Step1:  if both are unknown
Step2:  Now repeat all this updating the prior with the posterior one count at a time
Step3:  How does the posterior compare to just riunning N points?
<ASSISTANT_TASK:>
Python Code:
# generate some data
with pm.Model() as model:
    x = pm.Cauchy(name='x', alpha=0, beta=1)
    trace = pm.sample(10000, njobs=4)
    pm.traceplot(trace)
sampledat = trace['x']
trace.varnames, trace['x']
sns.distplot(sampledat, kde=False, norm_hist=True)
# plt.hist(sampledat, 200, normed=True);
plt.yscale('log');
np.random.randint(0, len(sampledat), 10)
# generate some data
bins = np.linspace(-4,4,100)
hists = {}
stats = {}
for npts in tqdm.tqdm_notebook(range(1,102,40)):
    d1 = sampledat[np.random.randint(0, len(sampledat), npts)]
    with pm.Model() as model:
        alpha = pm.Uniform('loc', -10, 10)
        #     beta = pm.Uniform('dist', 1, 1)
        x = pm.Cauchy(name='x', alpha=alpha, beta=1, observed=d1)
        trace = pm.sample(5000, njobs=4)
        hists[npts] = np.histogram(trace['loc'], bins)
        stats[npts] = np.percentile(trace['loc'], (1, 5, 25, 50, 75, 95, 99))
keys = sorted(list(hists.keys()))
for k in keys:
    p = plt.plot(tb.bin_edges_to_center(bins), hists[k][0]/np.max(hists[k][0]), 
                 drawstyle='steps', label=str(k), lw=1)
    c = p[0].get_color()
    plt.axvline(stats[k][3], lw=3, color=c)
    print(k, stats[k][2:5], stats[k][3]/(stats[k][4]-stats[k][2]), )
plt.legend()
plt.xlim((-2,2))
Explanation: Shots are fired isotropically from a point and hit a position sensitive detector
There is no scattering
y is fixed to be 1 away
End of explanation
# generate some data
bins = np.linspace(-4,4,100)
hists2 = {}
stats2 = {}
hists2d = {}
binsd = np.linspace(0.1,5,100)
for npts in tqdm.tqdm_notebook((1,2,5,10,20,40,60,80,200)):
    d1 = sampledat[np.random.randint(0, len(sampledat), npts)]
    with pm.Model() as model:
        alpha = pm.Uniform('loc', -10, 10)
        beta = pm.Uniform('dist', 0.1, 5)
        x = pm.Cauchy(name='x', alpha=alpha, beta=beta, observed=d1)
        trace = pm.sample(5000, njobs=4)
        hists2[npts] = np.histogram(trace['loc'], bins)
        stats2[npts] = np.percentile(trace['loc'], (1, 5, 25, 50, 75, 95, 99))
        hists2d[npts] = np.histogram2d(trace['loc'], trace['dist'], bins=(bins, binsd))
keys = sorted(list(hists2.keys()))
for k in keys:
    p = plt.plot(tb.bin_edges_to_center(bins), hists2[k][0]/np.max(hists2[k][0]), 
                 drawstyle='steps', label=str(k), lw=1)
    c = p[0].get_color()
    plt.axvline(stats2[k][3], lw=3, color=c)
    print(k, stats2[k][2:5], stats2[k][3]/(stats2[k][4]-stats2[k][2]), )
plt.legend()
plt.xlim((-2,2))
# plt.contour(hists2d[1][0], 5)
from matplotlib.colors import LogNorm
keys = sorted(list(hists2.keys()))
for k in keys:
    plt.figure()
    plt.pcolormesh(tb.bin_edges_to_center(binsd), 
                   tb.bin_edges_to_center(bins),
                   hists2d[k][0], 
                   norm=LogNorm())
    plt.title(str(k))
    plt.colorbar()
    plt.axvline(1, lw=0.5, c='k')
    plt.axhline(0, lw=0.5, c='k')
Explanation: if both are unknown
End of explanation
def from_posterior(param, samples):
    smin, smax = np.min(samples), np.max(samples)
    width = smax - smin
    x = np.linspace(smin, smax, 100)
    y = stats.gaussian_kde(samples)(x)
    
    # what was never sampled should have a small probability but not 0,
    # so we'll extend the domain and use linear approximation of density on it
    x = np.concatenate([[x[0] - 3 * width], x, [x[-1] + 3 * width]])
    y = np.concatenate([[0], y, [0]])
    return Interpolated(param, x, y)
# sampledat is the data
dat2 = sampledat.copy().tolist()
# traces = []
with pm.Model() as model:
    data = dat2.pop()
    data = [2.5, 2.5] # start kinda bad
    alpha = pm.Uniform('loc', -10, 10)
    #     beta = pm.Uniform('dist', 1, 1)
    x = pm.Cauchy(name='x', alpha=alpha, beta=1, observed=data)
    trace = pm.sample(5000, njobs=4)
    
traces= [trace]
sns.distplot(traces[-1]['loc'][1000:], kde=False)
plt.axvline(0, c='k', linestyle='--')
plt.axvline(data[0], c='r')
plt.axvline(data[1], c='r')
plt.title('{} points'.format(len(traces)))
plt.xlim((-4,4))
alpha
traces[-1].varnames
for _ in tqdm.tqdm_notebook(range(20)):
    with pm.Model() as model:
        data = [dat2.pop(), dat2.pop()]
        alpha = from_posterior('loc', trace['loc'])
        x = pm.Cauchy(name='x', alpha=alpha, beta=1, observed=data)
        trace = pm.sample(5000, njobs=4)
        traces.append(trace)
for ii, t in enumerate(traces):
    plt.figure()
    sns.distplot(t['loc'][1000:], kde=False, norm_hist=True)
    plt.axvline(0, c='k', linestyle='--')
    plt.title('{} points'.format((ii+1)*2))
    plt.xlim((-4,4))
    plt.show()
for ii, t in enumerate(traces):
    sns.distplot(t['loc'][1000:], kde=True, hist=False, color=plt.cm.OrRd(ii/len(traces)))
    plt.axvline(0, c='k', linestyle='--')
plt.title('{} points'.format(ii))
plt.xlim((-4,4))
plt.ylim((0, 1.3))
    
plt.figure()
sns.distplot(traces[-1]['loc'][1000:], kde=True, hist=False, color=plt.cm.OrRd(ii/len(traces)))
for t in traces:
    pm.plot_posterior(t)
for t in traces:
    print(pm.gelman_rubin(t))
Explanation: Now repeat all this updating the prior with the posterior one count at a time
End of explanation
# traces = []
with pm.Model() as model:
    alpha = pm.Uniform('loc', -10, 10)
    #     beta = pm.Uniform('dist', 1, 1)
    x = pm.Cauchy(name='x', alpha=alpha, beta=1, observed=sampledat[:42])
    trace = pm.sample(5000, njobs=4)
    
pm.plot_posterior(trace)
sns.distplot(trace['loc'][1000:], kde=True, hist=False, color=plt.cm.OrRd(ii/len(traces)))
plt.axvline(0, c='k', linestyle='--')
plt.xlim((-4,4))
plt.ylim((0, 1.3))
    
plt.figure()
sns.distplot(trace['loc'][1000::10], kde=True, hist=False, color=plt.cm.OrRd(ii/len(traces)))
plt.axvline(0, c='k', linestyle='--')
plt.xlim((-4,4))
plt.ylim((0, 1.3))
pm.forestplot(trace)
pm.autocorrplot(trace)
pm.traceplot(trace)
pm.gelman_rubin(trace)
pm.energyplot(trace)
energy = trace['loc']
energy_diff = np.diff(energy)
sns.distplot(energy - energy.mean(), label='energy')
sns.distplot(energy_diff, label='energy diff')
plt.legend()
Explanation: How does the posterior compare to just riunning N points?
End of explanation
<END_TASK>