In [3]:
import numpy as np
import matplotlib.pyplot as plt
import pymc3 as pm
import tqdm
import spacepy.toolbox as tb

%matplotlib inline

Setup a convolution data set and try and doconvolve it


In [39]:
np.random.seed(8675309)

dat_len = 100
xval = np.arange(dat_len)
realdat = np.zeros(dat_len, dtype=int)
realdat[40:60] = 50
noisemean = 2
real_n = np.zeros_like(realdat)
for i in range(len(realdat)):
    real_n[i] = np.random.poisson(realdat[i]+noisemean)


# make a detector
# triangular with FWFM 5 and is square
det = np.array([1,1,1,1,1])

# the numpy convolve I don't understand the normalization
obs = np.convolve(real_n, det, mode='same')
obs = tb.normalize(obs)
obs *= real_n.max()
    
# Two subplots, the axes array is 1-d
f, axarr = plt.subplots(4, sharex=True)
axarr[0].plot(xval, realdat)
axarr[0].set_ylabel('Truth')
axarr[0].set_ylim((0,60))

axarr[1].plot(xval, real_n)
axarr[1].set_ylabel('T+N')

axarr[2].plot(np.arange(len(det)), det)
axarr[2].set_ylabel('Det')

axarr[3].plot(xval, obs)
axarr[3].set_ylabel('Obs')


Out[39]:
<matplotlib.text.Text at 0x117666668>

So the det provides a point spread function that is U(0,5)


In [44]:
# generate some data
with pm.Model() as model:
    truth_mc = pm.Uniform('truth', 0, 100, shape=dat_len)
    noisemean_mc = pm.Uniform('noisemean', 0, 100)
    noise_mc = pm.Poisson('noise', noisemean_mc, observed=obs[1:20])
    real_n_mc = pm.Poisson('real_n', truth_mc+noisemean_mc, shape=dat_len)
    psf = pm.Uniform('psf', 0, 5, observed=det)
    obs_mc = pm.Normal('obs', (truth_mc+noisemean_mc)*psf.max(), 1/5**2, observed=obs, shape=dat_len)
    
    trace = pm.sample(5000)


Applied interval-transform to truth and added transformed truth_interval_ to model.
Applied interval-transform to noisemean and added transformed noisemean_interval_ to model.
Assigned NUTS to truth_interval_
Assigned NUTS to noisemean_interval_
Assigned Metropolis to real_n
 [-----------------100%-----------------] 5000 of 5000 complete in 20.9 sec

In [45]:
pm.traceplot(trace)


Out[45]:
array([[<matplotlib.axes._subplots.AxesSubplot object at 0x186786e48>,
        <matplotlib.axes._subplots.AxesSubplot object at 0x1867c50f0>],
       [<matplotlib.axes._subplots.AxesSubplot object at 0x18679c128>,
        <matplotlib.axes._subplots.AxesSubplot object at 0x18839e438>],
       [<matplotlib.axes._subplots.AxesSubplot object at 0x1883d9588>,
        <matplotlib.axes._subplots.AxesSubplot object at 0x188513a58>]], dtype=object)

In [46]:
pm.summary(trace)


real_n:

  Mean             SD               MC Error         95% HPD interval
  -------------------------------------------------------------------
  
  22.053           24.462           2.442            [0.000, 81.000]
  14.330           20.935           2.088            [0.000, 72.000]
  13.760           20.308           2.024            [0.000, 69.000]
  34.145           26.765           2.672            [5.000, 98.000]
  35.526           28.548           2.851            [0.000, 89.000]
  35.915           31.123           3.109            [1.000, 92.000]
  39.125           23.811           2.376            [6.000, 94.000]
  26.210           29.557           2.951            [0.000, 98.000]
  32.137           26.551           2.651            [0.000, 89.000]
  51.525           25.033           2.499            [22.000, 98.000]
  28.203           24.434           2.438            [4.000, 86.000]
  29.116           28.007           2.797            [0.000, 88.000]
  23.918           22.783           2.273            [0.000, 88.000]
  17.391           23.780           2.374            [0.000, 73.000]
  30.831           23.291           2.324            [4.000, 75.000]
  18.621           19.696           1.962            [2.000, 71.000]
  18.845           23.206           2.315            [0.000, 80.000]
  38.192           28.921           2.888            [4.000, 99.000]
  21.425           22.726           2.268            [2.000, 86.000]
  15.395           26.530           2.649            [0.000, 87.000]
  44.221           31.486           3.145            [7.000, 104.000]
  24.683           30.797           3.075            [0.000, 95.000]
  28.238           28.203           2.816            [0.000, 82.000]
  15.051           25.734           2.568            [0.000, 92.000]
  19.440           26.940           2.689            [0.000, 85.000]
  25.200           32.855           3.283            [0.000, 97.000]
  36.478           23.500           2.345            [4.000, 80.000]
  84.666           14.513           1.445            [61.000, 108.000]
  49.683           23.997           2.395            [15.000, 97.000]
  30.249           26.686           2.664            [7.000, 98.000]
  16.726           26.158           2.611            [0.000, 95.000]
  32.023           28.994           2.895            [1.000, 86.000]
  13.603           22.480           2.243            [0.000, 79.000]
  24.622           27.304           2.726            [0.000, 85.000]
  25.601           27.385           2.733            [5.000, 99.000]
  24.019           24.766           2.472            [1.000, 84.000]
  16.596           28.983           2.893            [0.000, 88.000]
  21.455           30.883           3.084            [0.000, 96.000]
  40.954           27.545           2.750            [13.000, 98.000]
  56.690           20.007           1.997            [26.000, 101.000]
  81.299           11.336           1.123            [66.000, 104.000]
  56.195           15.371           1.530            [41.000, 94.000]
  70.756           12.306           1.221            [56.000, 102.000]
  72.599           7.509            0.738            [63.000, 86.000]
  82.313           13.401           1.332            [59.000, 104.000]
  64.521           11.567           1.148            [53.000, 92.000]
  63.750           11.377           1.126            [47.000, 89.000]
  97.534           12.343           1.226            [80.000, 122.000]
  86.530           9.382            0.927            [68.000, 102.000]
  63.018           14.112           1.402            [48.000, 99.000]
  56.993           14.807           1.474            [41.000, 94.000]
  67.760           14.903           1.482            [43.000, 105.000]
  56.418           11.878           1.177            [37.000, 85.000]
  54.140           11.229           1.112            [43.000, 83.000]
  80.555           9.300            0.918            [66.000, 100.000]
  68.405           7.509            0.737            [55.000, 82.000]
  65.314           16.588           1.652            [45.000, 105.000]
  47.815           13.592           1.351            [32.000, 77.000]
  97.344           13.530           1.344            [81.000, 123.000]
  42.584           25.191           2.515            [16.000, 96.000]
  57.487           24.832           2.478            [21.000, 101.000]
  40.131           23.701           2.364            [15.000, 82.000]
  34.427           24.539           2.449            [4.000, 83.000]
  29.584           27.764           2.772            [4.000, 94.000]
  35.686           29.037           2.901            [6.000, 98.000]
  42.929           28.687           2.864            [0.000, 87.000]
  31.472           30.703           3.066            [6.000, 99.000]
  22.394           30.721           3.069            [0.000, 94.000]
  27.532           31.583           3.155            [0.000, 96.000]
  26.002           29.122           2.909            [0.000, 97.000]
  38.490           36.835           3.680            [0.000, 107.000]
  20.623           27.526           2.749            [0.000, 96.000]
  26.437           26.662           2.661            [6.000, 97.000]
  30.634           30.040           3.001            [0.000, 87.000]
  33.572           27.086           2.704            [0.000, 75.000]
  34.271           23.348           2.330            [11.000, 82.000]
  19.073           24.490           2.445            [0.000, 83.000]
  29.633           30.516           3.048            [4.000, 107.000]
  37.839           28.506           2.846            [7.000, 100.000]
  30.362           27.145           2.711            [7.000, 88.000]
  17.108           27.037           2.699            [0.000, 91.000]
  26.405           29.852           2.982            [3.000, 100.000]
  36.351           27.248           2.720            [12.000, 104.000]
  35.192           28.412           2.836            [3.000, 94.000]
  37.813           24.689           2.464            [12.000, 92.000]
  28.583           35.801           3.577            [0.000, 103.000]
  19.426           28.539           2.849            [0.000, 84.000]
  36.352           28.256           2.822            [2.000, 86.000]
  22.870           23.092           2.304            [0.000, 76.000]
  30.740           24.681           2.463            [5.000, 88.000]
  14.257           24.240           2.420            [0.000, 82.000]
  43.280           33.272           3.324            [10.000, 107.000]
  30.923           22.272           2.222            [8.000, 86.000]
  36.913           26.332           2.629            [0.000, 91.000]
  15.735           28.512           2.846            [0.000, 99.000]
  44.611           23.196           2.315            [9.000, 81.000]
  43.327           28.983           2.895            [10.000, 98.000]
  25.525           24.447           2.440            [0.000, 75.000]
  18.636           23.068           2.301            [0.000, 79.000]
  22.473           26.082           2.604            [0.000, 76.000]

  Posterior quantiles:
  2.5            25             50             75             97.5
  |--------------|==============|==============|--------------|
  
  1.000          6.000          11.000         28.000         91.000
  1.000          3.000          6.000          14.000         85.000
  0.000          3.000          6.000          12.000         79.000
  6.000          14.000         29.000         45.000         101.000
  1.000          11.000         29.000         61.000         97.000
  2.000          8.000          27.000         67.000         95.000
  10.000         23.000         29.000         50.000         103.000
  2.000          7.000          13.000         30.000         103.000
  0.000          17.000         23.000         41.000         98.000
  23.000         30.000         44.000         73.000         100.000
  7.000          12.000         17.000         34.000         92.000
  1.000          9.000          14.000         49.000         92.000
  1.000          13.000         16.000         24.000         97.000
  0.000          2.000          5.000          24.000         83.000
  7.000          15.000         18.000         47.000         91.000
  6.000          9.000          12.000         16.000         84.000
  0.000          4.000          13.000         20.000         93.000
  4.000          14.000         31.000         58.000         101.000
  4.000          7.000          15.000         21.000         95.000
  0.000          1.000          3.000          16.000         92.000
  9.000          19.000         31.000         66.000         115.000
  0.000          5.000          9.000          35.000         100.000
  0.000          4.000          22.000         41.000         85.000
  0.000          3.000          6.000          10.000         96.000
  0.000          3.000          6.000          32.000         101.000
  0.000          4.000          9.000          32.000         99.000
  6.000          16.000         32.000         54.000         87.000
  62.000         74.000         81.000         97.000         111.000
  17.000         31.000         46.000         68.000         102.000
  9.000          14.000         19.000         30.000         102.000
  0.000          4.000          6.000          12.000         103.000
  1.000          6.000          21.000         58.000         88.000
  0.000          2.000          5.000          11.000         91.000
  0.000          4.000          8.000          41.000         95.000
  6.000          11.000         15.000         22.000         102.000
  2.000          5.000          11.000         34.000         90.000
  0.000          2.000          3.000          7.000          104.000
  0.000          3.000          6.000          28.000         101.000
  13.000         18.000         32.000         58.000         102.000
  27.000         46.000         52.000         61.000         106.000
  66.000         72.000         80.000         88.000         105.000
  42.000         47.000         51.000         58.000         97.000
  57.000         64.000         68.000         72.000         110.000
  64.000         68.000         71.000         74.000         100.000
  60.000         74.000         81.000         93.000         106.000
  54.000         58.000         60.000         64.000         99.000
  49.000         56.000         61.000         67.000         95.000
  80.000         90.000         96.000         102.000        123.000
  68.000         81.000         87.000         93.000         103.000
  50.000         55.000         58.000         63.000         102.000
  43.000         48.000         52.000         57.000         97.000
  47.000         60.000         64.000         72.000         110.000
  40.000         48.000         55.000         60.000         90.000
  45.000         48.000         50.000         55.000         88.000
  67.000         73.000         80.000         85.000         103.000
  57.000         63.000         67.000         73.000         88.000
  46.000         53.000         63.000         69.000         107.000
  35.000         39.000         42.000         53.000         89.000
  83.000         86.000         93.000         106.000        127.000
  18.000         25.000         34.000         48.000         104.000
  23.000         38.000         51.000         81.000         109.000
  15.000         20.000         32.000         59.000         85.000
  5.000          15.000         28.000         47.000         86.000
  5.000          11.000         18.000         34.000         98.000
  7.000          13.000         21.000         53.000         100.000
  1.000          17.000         42.000         73.000         94.000
  7.000          13.000         17.000         30.000         102.000
  0.000          1.000          5.000          37.000         98.000
  0.000          4.000          12.000         38.000         98.000
  0.000          4.000          18.000         39.000         99.000
  0.000          10.000         21.000         72.000         110.000
  0.000          3.000          7.000          26.000         101.000
  8.000          12.000         15.000         26.000         103.000
  1.000          5.000          16.000         60.000         99.000
  0.000          6.000          27.000         57.000         84.000
  13.000         18.000         22.000         42.000         91.000
  1.000          5.000          8.000          22.000         88.000
  5.000          9.000          12.000         47.000         111.000
  7.000          15.000         30.000         45.000         102.000
  9.000          12.000         16.000         52.000         96.000
  0.000          3.000          7.000          11.000         99.000
  2.000          6.000          14.000         31.000         100.000
  12.000         18.000         23.000         38.000         105.000
  4.000          7.000          30.000         57.000         97.000
  13.000         18.000         30.000         44.000         95.000
  0.000          4.000          12.000         48.000         109.000
  0.000          1.000          4.000          27.000         102.000
  3.000          11.000         33.000         59.000         96.000
  0.000          5.000          14.000         33.000         88.000
  6.000          11.000         24.000         41.000         91.000
  0.000          2.000          5.000          9.000          91.000
  11.000         19.000         24.000         78.000         110.000
  9.000          18.000         21.000         41.000         95.000
  1.000          18.000         35.000         49.000         94.000
  0.000          1.000          4.000          9.000          101.000
  12.000         25.000         43.000         64.000         99.000
  10.000         18.000         41.000         64.000         102.000
  0.000          4.000          22.000         38.000         97.000
  0.000          4.000          11.000         21.000         91.000
  0.000          2.000          9.000          45.000         82.000


truth:

  Mean             SD               MC Error         95% HPD interval
  -------------------------------------------------------------------
  
  12.309           9.877            0.952            [0.019, 33.633]
  8.624            8.599            0.829            [0.033, 30.394]
  8.551            8.626            0.831            [0.001, 30.228]
  17.914           9.920            0.952            [2.969, 38.229]
  18.077           11.113           1.077            [0.121, 36.938]
  17.875           12.076           1.171            [1.236, 38.061]
  20.305           8.545            0.801            [5.875, 38.322]
  13.915           11.116           1.081            [0.559, 38.815]
  16.824           10.236           0.984            [0.012, 36.444]
  24.613           8.219            0.768            [12.641, 40.581]
  15.505           9.227            0.880            [3.480, 35.780]
  15.431           10.975           1.066            [0.017, 36.367]
  13.890           8.696            0.827            [0.018, 34.123]
  9.955            10.262           0.999            [0.004, 32.097]
  16.937           9.003            0.856            [4.019, 35.513]
  11.576           7.702            0.721            [2.729, 32.347]
  10.984           9.704            0.938            [0.003, 33.229]
  19.261           10.683           1.030            [2.408, 39.395]
  12.422           8.824            0.847            [0.938, 34.333]
  8.419            11.009           1.076            [0.004, 34.979]
  21.800           10.825           1.042            [5.910, 42.549]
  12.926           11.976           1.168            [0.028, 38.184]
  14.928           11.775           1.145            [0.005, 35.760]
  8.787            10.260           0.996            [0.008, 36.061]
  10.737           11.056           1.077            [0.018, 35.064]
  12.867           12.650           1.238            [0.000, 39.297]
  19.362           9.266            0.881            [3.627, 35.836]
  35.301           4.985            0.381            [25.987, 44.650]
  24.456           8.309            0.776            [10.768, 40.488]
  16.702           9.529            0.908            [5.536, 40.493]
  9.521            10.368           1.008            [0.007, 35.437]
  16.405           11.755           1.144            [0.010, 35.975]
  8.149            9.263            0.898            [0.001, 32.002]
  13.131           11.336           1.107            [0.026, 34.313]
  14.273           9.926            0.952            [2.571, 39.071]
  13.580           10.273           0.994            [0.688, 35.631]
  8.766            11.616           1.139            [0.000, 36.400]
  11.216           12.061           1.181            [0.011, 38.276]
  25.828           10.344           0.976            [10.651, 46.125]
  37.995           7.081            0.610            [24.929, 53.398]
  52.533           4.956            0.293            [43.706, 62.287]
  54.524           6.072            0.424            [44.630, 68.008]
  66.202           5.283            0.289            [56.142, 76.963]
  64.745           4.693            0.179            [56.069, 74.125]
  68.216           5.551            0.343            [57.711, 78.464]
  62.950           5.339            0.296            [52.770, 73.336]
  62.801           5.357            0.302            [52.596, 73.278]
  70.818           5.241            0.306            [60.670, 80.464]
  66.443           5.068            0.240            [57.187, 76.242]
  58.378           5.765            0.367            [47.914, 70.286]
  56.179           5.961            0.403            [45.433, 68.850]
  57.658           5.783            0.384            [47.176, 69.876]
  57.163           5.538            0.326            [47.430, 68.424]
  57.265           5.359            0.305            [47.385, 68.013]
  65.680           4.818            0.227            [56.735, 74.712]
  65.358           4.756            0.186            [56.408, 74.551]
  62.219           6.172            0.438            [51.078, 75.080]
  56.370           5.877            0.383            [45.759, 68.338]
  62.897           5.262            0.330            [53.185, 72.900]
  38.032           9.227            0.833            [24.265, 57.859]
  35.893           8.663            0.791            [20.648, 51.593]
  25.201           9.347            0.876            [11.362, 42.212]
  18.255           9.488            0.903            [3.201, 36.595]
  15.989           10.384           1.001            [2.836, 38.910]
  18.363           10.730           1.035            [4.027, 39.214]
  20.821           11.081           1.072            [0.040, 37.225]
  16.500           10.889           1.053            [4.114, 40.408]
  11.577           12.429           1.219            [0.001, 37.907]
  14.272           12.578           1.232            [0.005, 38.792]
  13.918           11.778           1.148            [0.012, 37.913]
  18.558           13.560           1.328            [0.004, 42.082]
  11.376           11.105           1.083            [0.000, 36.239]
  14.967           9.706            0.928            [4.189, 39.069]
  15.831           12.195           1.189            [0.006, 36.654]
  17.404           11.525           1.119            [0.009, 34.393]
  18.680           8.546            0.803            [7.176, 36.516]
  10.984           10.026           0.968            [0.241, 34.138]
  15.468           11.177           1.088            [2.569, 40.587]
  19.182           10.259           0.984            [3.980, 39.768]
  16.126           10.060           0.970            [4.577, 36.809]
  9.499            10.576           1.029            [0.001, 36.526]
  13.965           11.317           1.101            [0.010, 38.502]
  18.960           9.447            0.902            [6.606, 40.248]
  18.174           11.297           1.093            [1.959, 37.792]
  20.033           8.880            0.836            [7.318, 39.003]
  14.554           13.616           1.337            [0.001, 41.028]
  10.493           12.040           1.178            [0.003, 35.500]
  18.948           11.117           1.073            [1.877, 38.341]
  13.126           9.834            0.947            [0.001, 32.916]
  16.984           9.490            0.903            [2.836, 36.503]
  8.389            9.895            0.963            [0.004, 34.537]
  21.262           11.274           1.086            [6.940, 43.266]
  17.516           8.316            0.778            [5.823, 36.800]
  19.265           10.686           1.029            [0.008, 37.187]
  8.718            11.401           1.115            [0.002, 38.409]
  22.862           8.550            0.802            [7.810, 38.186]
  22.012           10.522           1.009            [6.144, 40.820]
  14.360           10.528           1.020            [0.006, 33.443]
  11.157           9.534            0.917            [0.004, 33.874]
  12.120           11.192           1.090            [0.007, 32.810]

  Posterior quantiles:
  2.5            25             50             75             97.5
  |--------------|==============|==============|--------------|
  
  1.045          5.260          8.796          16.703         36.619
  0.417          2.963          5.726          10.292         34.277
  0.352          2.884          5.595          10.094         33.328
  4.143          9.884          16.568         23.366         40.068
  0.947          8.648          17.279         27.021         39.149
  1.841          6.590          15.894         29.320         39.302
  6.960          14.006         18.741         25.635         39.938
  1.340          5.784          9.966          18.605         41.077
  0.432          10.916         14.854         21.800         38.874
  12.815         17.981         22.805         31.449         40.976
  4.528          8.687          12.460         19.718         37.312
  1.000          7.031          11.084         24.375         38.498
  1.229          8.753          11.681         16.040         37.724
  0.173          2.053          5.669          14.960         34.567
  4.696          9.989          13.969         23.641         36.565
  3.679          6.910          9.446          12.619         35.131
  0.264          3.459          8.913          14.109         36.185
  3.162          10.398         18.215         26.977         40.474
  2.308          6.348          10.194         14.526         37.477
  0.099          1.484          3.303          10.494         37.917
  6.816          13.039         18.509         30.531         43.701
  0.554          4.273          7.612          19.795         40.186
  0.375          4.207          13.529         22.737         37.484
  0.217          2.709          5.323          8.946          39.029
  0.524          2.842          5.669          17.235         38.001
  0.308          3.640          7.707          18.267         41.017
  4.430          11.674         18.955         26.396         36.937
  26.087         31.800         35.264         38.822         44.859
  10.820         17.793         23.747         30.669         40.715
  5.835          10.258         13.638         19.333         41.142
  0.329          3.112          5.601          10.147         39.534
  0.904          5.756          13.974         27.177         37.947
  0.113          2.327          4.975          9.107          35.609
  0.545          3.612          7.995          22.574         37.280
  3.976          8.087          10.929         15.517         41.013
  1.244          4.989          10.503         20.076         37.679
  0.125          1.765          3.641          7.750          39.718
  0.297          2.695          5.769          16.841         40.396
  11.371         17.416         23.408         33.336         47.416
  25.474         33.336         37.212         41.780         54.079
  43.571         49.172         52.404         55.755         62.265
  44.555         50.505         53.755         58.050         68.001
  56.825         62.642         65.874         69.321         77.935
  56.187         61.588         64.671         67.696         74.410
  57.733         64.553         68.295         71.988         78.497
  53.473         59.392         62.650         66.246         74.249
  52.783         59.197         62.651         66.155         73.558
  61.186         67.230         70.695         74.297         81.282
  57.022         63.113         66.475         69.731         76.136
  48.688         54.491         57.890         61.521         71.556
  46.033         52.194         55.411         59.473         69.918
  47.359         53.747         57.396         61.186         70.084
  47.430         53.516         56.891         60.618         68.424
  47.700         53.715         56.933         60.508         68.548
  56.735         62.423         65.717         68.913         74.712
  56.526         62.275         65.245         68.347         74.950
  51.289         58.096         61.657         66.005         75.464
  46.357         52.309         55.918         59.838         69.368
  53.533         59.430         62.660         66.352         73.504
  24.433         31.485         35.897         42.754         58.117
  21.372         29.210         34.534         43.106         52.427
  11.892         17.289         23.088         33.329         42.881
  3.746          11.068         16.857         24.366         37.372
  3.335          8.353          13.146         20.043         40.024
  4.795          9.405          14.499         26.572         40.369
  1.106          11.907         21.802         30.348         38.860
  4.549          9.161          12.472         18.862         41.244
  0.170          1.852          5.469          19.817         40.037
  0.292          3.575          9.845          21.878         40.453
  0.294          3.639          11.844         20.876         40.226
  0.567          7.498          14.771         31.007         43.712
  0.236          2.884          6.988          16.421         39.903
  5.098          8.754          11.506         16.790         41.048
  0.575          4.964          12.378         27.004         38.781
  0.419          5.938          17.902         27.684         36.342
  8.052          12.288         15.685         23.502         38.039
  1.084          4.122          6.897          13.818         36.557
  3.436          7.033          10.535         23.018         42.503
  4.731          10.794         17.750         24.629         41.212
  5.311          8.890          11.718         23.996         38.179
  0.256          2.863          5.933          9.955          39.124
  1.129          5.492          10.467         18.095         40.727
  7.112          12.079         16.082         22.908         41.050
  2.682          7.038          17.683         27.803         39.527
  7.878          12.924         18.311         24.528         39.791
  0.243          4.203          9.441          24.281         42.967
  0.125          1.519          4.376          17.001         38.714
  2.365          9.038          18.572         27.857         39.500
  0.596          4.862          10.816         19.702         35.887
  4.166          9.444          15.040         22.648         38.680
  0.258          2.370          4.845          8.494          37.579
  7.162          12.536         16.593         31.514         43.582
  6.822          11.854         14.901         21.187         38.167
  0.787          11.798         19.942         26.231         39.441
  0.127          1.717          3.946          8.822          40.853
  8.326          15.902         22.601         29.466         39.058
  7.064          12.581         20.895         30.440         42.230
  0.603          4.525          14.016         21.896         37.592
  0.499          4.229          8.771          14.294         37.056
  0.236          2.816          7.504          22.498         35.409


noisemean:

  Mean             SD               MC Error         95% HPD interval
  -------------------------------------------------------------------
  
  0.328            1.127            0.077            [0.078, 0.464]

  Posterior quantiles:
  2.5            25             50             75             97.5
  |--------------|==============|==============|--------------|
  
  0.087          0.169          0.234          0.320          0.493


In [50]:
# plt.plot(trace['truth'][0:5,:].T)
trace['truth'].shape
iqr = np.zeros((dat_len,3))
for i in range(dat_len):
    iqr[i] = np.percentile(trace['truth'].T[i], (25,50,75), axis=0)
plt.plot(xval, iqr[:,1],  label='recovered')
plt.fill_between(xval, iqr[:,0], iqr[:,2], alpha=0.2)
plt.plot(xval, real_n, c='r', label='pre psf')
plt.plot(xval, realdat, c='g',  label='truth')
plt.plot(xval, obs, c='k',  label='observed', lw=3)

plt.legend()
plt.figure()
snr = iqr[:,1]/(iqr[:,2], iqr[:,0])
perixval.shape, snr.shape
plt.plot(xval, snr)


---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-50-5a4d4d0734e1> in <module>()
     13 plt.figure()
     14 snr = iqr[:,1]/(iqr[:,2], iqr[:,0])
---> 15 plt.plot(xval, snr)

/Users/balarsen/miniconda3/envs/python3/lib/python3.5/site-packages/matplotlib/pyplot.py in plot(*args, **kwargs)
   3159         ax.hold(hold)
   3160     try:
-> 3161         ret = ax.plot(*args, **kwargs)
   3162     finally:
   3163         ax.hold(washold)

/Users/balarsen/miniconda3/envs/python3/lib/python3.5/site-packages/matplotlib/__init__.py in inner(ax, *args, **kwargs)
   1817                     warnings.warn(msg % (label_namer, func.__name__),
   1818                                   RuntimeWarning, stacklevel=2)
-> 1819             return func(ax, *args, **kwargs)
   1820         pre_doc = inner.__doc__
   1821         if pre_doc is None:

/Users/balarsen/miniconda3/envs/python3/lib/python3.5/site-packages/matplotlib/axes/_axes.py in plot(self, *args, **kwargs)
   1380         kwargs = cbook.normalize_kwargs(kwargs, _alias_map)
   1381 
-> 1382         for line in self._get_lines(*args, **kwargs):
   1383             self.add_line(line)
   1384             lines.append(line)

/Users/balarsen/miniconda3/envs/python3/lib/python3.5/site-packages/matplotlib/axes/_base.py in _grab_next_args(self, *args, **kwargs)
    379                 return
    380             if len(remaining) <= 3:
--> 381                 for seg in self._plot_args(remaining, kwargs):
    382                     yield seg
    383                 return

/Users/balarsen/miniconda3/envs/python3/lib/python3.5/site-packages/matplotlib/axes/_base.py in _plot_args(self, tup, kwargs)
    357             x, y = index_of(tup[-1])
    358 
--> 359         x, y = self._xy_from_xy(x, y)
    360 
    361         if self.command == 'plot':

/Users/balarsen/miniconda3/envs/python3/lib/python3.5/site-packages/matplotlib/axes/_base.py in _xy_from_xy(self, x, y)
    217         y = _check_1d(y)
    218         if x.shape[0] != y.shape[0]:
--> 219             raise ValueError("x and y must have same first dimension")
    220         if x.ndim > 2 or y.ndim > 2:
    221             raise ValueError("x and y can be no greater than 2-D")

ValueError: x and y must have same first dimension

In [48]:
print(np.percentile(trace['noisemean'], (25,50,75)), noisemean)


[ 0.16863086  0.23445815  0.3197733 ] 2

In [49]:
obs[1:20]


Out[49]:
array([ 0.        ,  0.72490706,  0.48327138,  0.72490706,  0.72490706,
        0.72490706,  0.48327138,  0.24163569,  0.        ,  0.48327138,
        0.96654275,  0.96654275,  1.44981413,  1.20817844,  0.72490706,
        1.20817844,  0.72490706,  0.24163569,  1.20817844])

In [ ]: