In [3]:
import numpy as np
import matplotlib.pyplot as plt
import pymc3 as pm
import tqdm
import spacepy.toolbox as tb
%matplotlib inline
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 [ ]:
Content source: balarsen/pymc_learning
Similar notebooks: