In [3]:
import numpy as np
import matplotlib.pyplot as plt
import pymc3 as pm
import tqdm
import spacepy.toolbox as tb
import spacepy.plot as spp
%matplotlib inline
/Users/balarsen/miniconda3/envs/python3/lib/python3.5/site-packages/matplotlib/__init__.py:878: UserWarning: axes.color_cycle is deprecated and replaced with axes.prop_cycle; please use the latter.
warnings.warn(self.msg_depr % (key, alt_key))
In [38]:
np.random.seed(8675309)
sim_pa = np.arange(20,175)
sim_c = 890*np.sin(np.deg2rad(sim_pa))**0.8
# at each point draw a poisson variable with that mean
sim_c_n = np.asarray([np.random.poisson(v) for v in sim_c ])
prob=0.1
sim_c_n2 = np.asarray([np.random.negative_binomial((v*prob)/(1-prob), prob) for v in sim_c ])
# Two subplots, the axes array is 1-d
f, axarr = plt.subplots(3, sharex=True, sharey=True, figsize=(7,9))
axarr[0].plot(sim_pa, sim_c, lw=2)
axarr[0].set_ylabel('Truth')
axarr[0].set_xlim((0,180))
axarr[0].set_yscale('log')
axarr[1].plot(sim_pa, sim_c_n, lw=2)
axarr[1].set_ylabel('Measured\nPoisson')
axarr[1].set_xlim((0,180))
axarr[1].set_yscale('log')
axarr[2].plot(sim_pa, sim_c_n2, lw=2)
axarr[2].set_ylabel('Measured\nNegBin')
axarr[2].set_xlim((0,180))
axarr[2].set_yscale('log')
axarr[2].set_ylim(axarr[1].get_ylim())
Out[38]:
(100.0, 10000.0)
In [54]:
# generate some data
with pm.Model() as model:
bkg = pm.NegativeBinomial('bkg', mu=pm.Uniform('m_bkg', 0, 1e5,shape=len(sim_pa), testval=1e3),
alpha=0.1,
observed=sim_c_n2, shape=len(sim_pa))
# 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(50000)
Applied interval-transform to m_bkg and added transformed m_bkg_interval_ to model.
Assigned NUTS to m_bkg_interval_
[-----------------100%-----------------] 50000 of 50000 complete in 9.0 sec
In [55]:
pm.traceplot(trace)
Out[55]:
array([[<matplotlib.axes._subplots.AxesSubplot object at 0x127d8fcc0>,
<matplotlib.axes._subplots.AxesSubplot object at 0x12f4259b0>]], dtype=object)
In [56]:
pm.summary(trace)
m_bkg:
Mean SD MC Error 95% HPD interval
-------------------------------------------------------------------
24518.784 0.000 0.000 [24518.784, 24518.784]
89061.869 0.000 0.000 [89061.869, 89061.869]
26730.966 0.000 0.000 [26730.966, 26730.966]
91726.542 0.000 0.000 [91726.542, 91726.542]
19092.396 0.000 0.000 [19092.396, 19092.396]
15813.599 0.000 0.000 [15813.599, 15813.599]
190.210 0.000 0.000 [190.210, 190.210]
85883.446 0.000 0.000 [85883.446, 85883.446]
13052.674 0.000 0.000 [13052.674, 13052.674]
4280.376 0.000 0.000 [4280.376, 4280.376]
5725.985 0.000 0.000 [5725.985, 5725.985]
89896.479 0.000 0.000 [89896.479, 89896.479]
3329.918 0.000 0.000 [3329.918, 3329.918]
1066.652 0.000 0.000 [1066.652, 1066.652]
3979.871 0.000 0.000 [3979.871, 3979.871]
4229.269 0.000 0.000 [4229.269, 4229.269]
78290.615 0.000 0.000 [78290.615, 78290.615]
75513.295 0.000 0.000 [75513.295, 75513.295]
69715.093 0.000 0.000 [69715.093, 69715.093]
7963.500 0.000 0.000 [7963.500, 7963.500]
3896.962 0.000 0.000 [3896.962, 3896.962]
78281.682 0.000 0.000 [78281.682, 78281.682]
1982.039 0.000 0.000 [1982.039, 1982.039]
312.858 0.000 0.000 [312.858, 312.858]
13637.283 0.000 0.000 [13637.283, 13637.283]
6438.958 0.000 0.000 [6438.958, 6438.958]
3880.202 0.000 0.000 [3880.202, 3880.202]
26386.489 0.000 0.000 [26386.489, 26386.489]
22539.140 0.000 0.000 [22539.140, 22539.140]
59746.479 0.000 0.000 [59746.479, 59746.479]
305.683 0.000 0.000 [305.683, 305.683]
40817.168 0.000 0.000 [40817.168, 40817.168]
11543.797 0.000 0.000 [11543.797, 11543.797]
1525.997 0.000 0.000 [1525.997, 1525.997]
1903.233 0.000 0.000 [1903.233, 1903.233]
79777.161 0.000 0.000 [79777.161, 79777.161]
1105.875 0.000 0.000 [1105.875, 1105.875]
30164.122 0.000 0.000 [30164.122, 30164.122]
20861.868 0.000 0.000 [20861.868, 20861.868]
19589.601 0.000 0.000 [19589.601, 19589.601]
3509.641 0.000 0.000 [3509.641, 3509.641]
90907.452 0.000 0.000 [90907.452, 90907.452]
6280.000 0.000 0.000 [6280.000, 6280.000]
68166.685 0.000 0.000 [68166.685, 68166.685]
1276.186 0.000 0.000 [1276.186, 1276.186]
1278.551 0.000 0.000 [1278.551, 1278.551]
172.852 0.000 0.000 [172.852, 172.852]
30181.398 0.000 0.000 [30181.398, 30181.398]
9227.869 0.000 0.000 [9227.869, 9227.869]
863.983 0.000 0.000 [863.983, 863.983]
10503.640 0.000 0.000 [10503.640, 10503.640]
39675.286 0.000 0.000 [39675.286, 39675.286]
3802.292 0.000 0.000 [3802.292, 3802.292]
672.958 0.000 0.000 [672.958, 672.958]
10163.156 0.000 0.000 [10163.156, 10163.156]
54687.818 0.000 0.000 [54687.818, 54687.818]
883.713 0.000 0.000 [883.713, 883.713]
9241.529 0.000 0.000 [9241.529, 9241.529]
14127.678 0.000 0.000 [14127.678, 14127.678]
5184.912 0.000 0.000 [5184.912, 5184.912]
2137.659 0.000 0.000 [2137.659, 2137.659]
51641.383 0.000 0.000 [51641.383, 51641.383]
63335.690 0.000 0.000 [63335.690, 63335.690]
16535.406 0.000 0.000 [16535.406, 16535.406]
720.711 0.000 0.000 [720.711, 720.711]
4447.431 0.000 0.000 [4447.431, 4447.431]
3445.203 0.000 0.000 [3445.203, 3445.203]
2319.363 0.000 0.000 [2319.363, 2319.363]
18386.615 0.000 0.000 [18386.615, 18386.615]
685.609 0.000 0.000 [685.609, 685.609]
1542.480 0.000 0.000 [1542.480, 1542.480]
2115.588 0.000 0.000 [2115.588, 2115.588]
3386.827 0.000 0.000 [3386.827, 3386.827]
64.710 0.000 0.000 [64.710, 64.710]
2530.634 0.000 0.000 [2530.634, 2530.634]
175.140 0.000 0.000 [175.140, 175.140]
74691.840 0.000 0.000 [74691.840, 74691.840]
32474.993 0.000 0.000 [32474.993, 32474.993]
18585.848 0.000 0.000 [18585.848, 18585.848]
40743.202 0.000 0.000 [40743.202, 40743.202]
15555.904 0.000 0.000 [15555.904, 15555.904]
1516.434 0.000 0.000 [1516.434, 1516.434]
2005.291 0.000 0.000 [2005.291, 2005.291]
9489.615 0.000 0.000 [9489.615, 9489.615]
7476.945 0.000 0.000 [7476.945, 7476.945]
18977.847 0.000 0.000 [18977.847, 18977.847]
31546.312 0.000 0.000 [31546.312, 31546.312]
36895.759 0.000 0.000 [36895.759, 36895.759]
35681.658 0.000 0.000 [35681.658, 35681.658]
46002.335 0.000 0.000 [46002.335, 46002.335]
49355.674 0.000 0.000 [49355.674, 49355.674]
94783.385 0.000 0.000 [94783.385, 94783.385]
2672.509 0.000 0.000 [2672.509, 2672.509]
21198.563 0.000 0.000 [21198.563, 21198.563]
7397.935 0.000 0.000 [7397.935, 7397.935]
25102.139 0.000 0.000 [25102.139, 25102.139]
3645.336 0.000 0.000 [3645.336, 3645.336]
8157.216 0.000 0.000 [8157.216, 8157.216]
78145.828 0.000 0.000 [78145.828, 78145.828]
421.356 0.000 0.000 [421.356, 421.356]
5459.179 0.000 0.000 [5459.179, 5459.179]
264.621 0.000 0.000 [264.621, 264.621]
252.255 0.000 0.000 [252.255, 252.255]
3725.515 0.000 0.000 [3725.515, 3725.515]
67013.741 0.000 0.000 [67013.741, 67013.741]
6426.947 0.000 0.000 [6426.947, 6426.947]
6618.732 0.000 0.000 [6618.732, 6618.732]
69948.402 0.000 0.000 [69948.402, 69948.402]
41692.932 0.000 0.000 [41692.932, 41692.932]
14469.020 0.000 0.000 [14469.020, 14469.020]
43630.078 0.000 0.000 [43630.078, 43630.078]
39365.821 0.000 0.000 [39365.821, 39365.821]
2725.494 0.000 0.000 [2725.494, 2725.494]
255.275 0.000 0.000 [255.275, 255.275]
3173.358 0.000 0.000 [3173.358, 3173.358]
5239.520 0.000 0.000 [5239.520, 5239.520]
3804.238 0.000 0.000 [3804.238, 3804.238]
4884.278 0.000 0.000 [4884.278, 4884.278]
44653.785 0.000 0.000 [44653.785, 44653.785]
28685.632 0.000 0.000 [28685.632, 28685.632]
27912.764 0.000 0.000 [27912.764, 27912.764]
50977.905 0.000 0.000 [50977.905, 50977.905]
81231.890 0.000 0.000 [81231.890, 81231.890]
372.418 0.000 0.000 [372.418, 372.418]
25860.739 0.000 0.000 [25860.739, 25860.739]
203.856 0.000 0.000 [203.856, 203.856]
61456.722 0.000 0.000 [61456.722, 61456.722]
87344.989 0.000 0.000 [87344.989, 87344.989]
1468.944 0.000 0.000 [1468.944, 1468.944]
9354.191 0.000 0.000 [9354.191, 9354.191]
4248.658 0.000 0.000 [4248.658, 4248.658]
4657.324 0.000 0.000 [4657.324, 4657.324]
14117.866 0.000 0.000 [14117.866, 14117.866]
57662.544 0.000 0.000 [57662.544, 57662.544]
6784.044 0.000 0.000 [6784.044, 6784.044]
37984.081 0.000 0.000 [37984.081, 37984.081]
137.558 0.000 0.000 [137.558, 137.558]
111.340 0.000 0.000 [111.340, 111.340]
29230.965 0.000 0.000 [29230.965, 29230.965]
74861.551 0.000 0.000 [74861.551, 74861.551]
3347.438 0.000 0.000 [3347.438, 3347.438]
30749.199 0.000 0.000 [30749.199, 30749.199]
16681.520 0.000 0.000 [16681.520, 16681.520]
6462.548 0.000 0.000 [6462.548, 6462.548]
6219.847 0.000 0.000 [6219.847, 6219.847]
37732.809 0.000 0.000 [37732.809, 37732.809]
2351.100 0.000 0.000 [2351.100, 2351.100]
10277.297 0.000 0.000 [10277.297, 10277.297]
76102.142 0.000 0.000 [76102.142, 76102.142]
85579.688 0.000 0.000 [85579.688, 85579.688]
87186.563 0.000 0.000 [87186.563, 87186.563]
92971.385 0.000 0.000 [92971.385, 92971.385]
95632.522 0.000 0.000 [95632.522, 95632.522]
6317.607 0.000 0.000 [6317.607, 6317.607]
34510.372 0.000 0.000 [34510.372, 34510.372]
Posterior quantiles:
2.5 25 50 75 97.5
|--------------|==============|==============|--------------|
24518.784 24518.784 24518.784 24518.784 24518.784
89061.869 89061.869 89061.869 89061.869 89061.869
26730.966 26730.966 26730.966 26730.966 26730.966
91726.542 91726.542 91726.542 91726.542 91726.542
19092.396 19092.396 19092.396 19092.396 19092.396
15813.599 15813.599 15813.599 15813.599 15813.599
190.210 190.210 190.210 190.210 190.210
85883.446 85883.446 85883.446 85883.446 85883.446
13052.674 13052.674 13052.674 13052.674 13052.674
4280.376 4280.376 4280.376 4280.376 4280.376
5725.985 5725.985 5725.985 5725.985 5725.985
89896.479 89896.479 89896.479 89896.479 89896.479
3329.918 3329.918 3329.918 3329.918 3329.918
1066.652 1066.652 1066.652 1066.652 1066.652
3979.871 3979.871 3979.871 3979.871 3979.871
4229.269 4229.269 4229.269 4229.269 4229.269
78290.615 78290.615 78290.615 78290.615 78290.615
75513.295 75513.295 75513.295 75513.295 75513.295
69715.093 69715.093 69715.093 69715.093 69715.093
7963.500 7963.500 7963.500 7963.500 7963.500
3896.962 3896.962 3896.962 3896.962 3896.962
78281.682 78281.682 78281.682 78281.682 78281.682
1982.039 1982.039 1982.039 1982.039 1982.039
312.858 312.858 312.858 312.858 312.858
13637.283 13637.283 13637.283 13637.283 13637.283
6438.958 6438.958 6438.958 6438.958 6438.958
3880.202 3880.202 3880.202 3880.202 3880.202
26386.489 26386.489 26386.489 26386.489 26386.489
22539.140 22539.140 22539.140 22539.140 22539.140
59746.479 59746.479 59746.479 59746.479 59746.479
305.683 305.683 305.683 305.683 305.683
40817.168 40817.168 40817.168 40817.168 40817.168
11543.797 11543.797 11543.797 11543.797 11543.797
1525.997 1525.997 1525.997 1525.997 1525.997
1903.233 1903.233 1903.233 1903.233 1903.233
79777.161 79777.161 79777.161 79777.161 79777.161
1105.875 1105.875 1105.875 1105.875 1105.875
30164.122 30164.122 30164.122 30164.122 30164.122
20861.868 20861.868 20861.868 20861.868 20861.868
19589.601 19589.601 19589.601 19589.601 19589.601
3509.641 3509.641 3509.641 3509.641 3509.641
90907.452 90907.452 90907.452 90907.452 90907.452
6280.000 6280.000 6280.000 6280.000 6280.000
68166.685 68166.685 68166.685 68166.685 68166.685
1276.186 1276.186 1276.186 1276.186 1276.186
1278.551 1278.551 1278.551 1278.551 1278.551
172.852 172.852 172.852 172.852 172.852
30181.398 30181.398 30181.398 30181.398 30181.398
9227.869 9227.869 9227.869 9227.869 9227.869
863.983 863.983 863.983 863.983 863.983
10503.640 10503.640 10503.640 10503.640 10503.640
39675.286 39675.286 39675.286 39675.286 39675.286
3802.292 3802.292 3802.292 3802.292 3802.292
672.958 672.958 672.958 672.958 672.958
10163.156 10163.156 10163.156 10163.156 10163.156
54687.818 54687.818 54687.818 54687.818 54687.818
883.713 883.713 883.713 883.713 883.713
9241.529 9241.529 9241.529 9241.529 9241.529
14127.678 14127.678 14127.678 14127.678 14127.678
5184.912 5184.912 5184.912 5184.912 5184.912
2137.659 2137.659 2137.659 2137.659 2137.659
51641.383 51641.383 51641.383 51641.383 51641.383
63335.690 63335.690 63335.690 63335.690 63335.690
16535.406 16535.406 16535.406 16535.406 16535.406
720.711 720.711 720.711 720.711 720.711
4447.431 4447.431 4447.431 4447.431 4447.431
3445.203 3445.203 3445.203 3445.203 3445.203
2319.363 2319.363 2319.363 2319.363 2319.363
18386.615 18386.615 18386.615 18386.615 18386.615
685.609 685.609 685.609 685.609 685.609
1542.480 1542.480 1542.480 1542.480 1542.480
2115.588 2115.588 2115.588 2115.588 2115.588
3386.827 3386.827 3386.827 3386.827 3386.827
64.710 64.710 64.710 64.710 64.710
2530.634 2530.634 2530.634 2530.634 2530.634
175.140 175.140 175.140 175.140 175.140
74691.840 74691.840 74691.840 74691.840 74691.840
32474.993 32474.993 32474.993 32474.993 32474.993
18585.848 18585.848 18585.848 18585.848 18585.848
40743.202 40743.202 40743.202 40743.202 40743.202
15555.904 15555.904 15555.904 15555.904 15555.904
1516.434 1516.434 1516.434 1516.434 1516.434
2005.291 2005.291 2005.291 2005.291 2005.291
9489.615 9489.615 9489.615 9489.615 9489.615
7476.945 7476.945 7476.945 7476.945 7476.945
18977.847 18977.847 18977.847 18977.847 18977.847
31546.312 31546.312 31546.312 31546.312 31546.312
36895.759 36895.759 36895.759 36895.759 36895.759
35681.658 35681.658 35681.658 35681.658 35681.658
46002.335 46002.335 46002.335 46002.335 46002.335
49355.674 49355.674 49355.674 49355.674 49355.674
94783.385 94783.385 94783.385 94783.385 94783.385
2672.509 2672.509 2672.509 2672.509 2672.509
21198.563 21198.563 21198.563 21198.563 21198.563
7397.935 7397.935 7397.935 7397.935 7397.935
25102.139 25102.139 25102.139 25102.139 25102.139
3645.336 3645.336 3645.336 3645.336 3645.336
8157.216 8157.216 8157.216 8157.216 8157.216
78145.828 78145.828 78145.828 78145.828 78145.828
421.356 421.356 421.356 421.356 421.356
5459.179 5459.179 5459.179 5459.179 5459.179
264.621 264.621 264.621 264.621 264.621
252.255 252.255 252.255 252.255 252.255
3725.515 3725.515 3725.515 3725.515 3725.515
67013.741 67013.741 67013.741 67013.741 67013.741
6426.947 6426.947 6426.947 6426.947 6426.947
6618.732 6618.732 6618.732 6618.732 6618.732
69948.402 69948.402 69948.402 69948.402 69948.402
41692.932 41692.932 41692.932 41692.932 41692.932
14469.020 14469.020 14469.020 14469.020 14469.020
43630.078 43630.078 43630.078 43630.078 43630.078
39365.821 39365.821 39365.821 39365.821 39365.821
2725.494 2725.494 2725.494 2725.494 2725.494
255.275 255.275 255.275 255.275 255.275
3173.358 3173.358 3173.358 3173.358 3173.358
5239.520 5239.520 5239.520 5239.520 5239.520
3804.238 3804.238 3804.238 3804.238 3804.238
4884.278 4884.278 4884.278 4884.278 4884.278
44653.785 44653.785 44653.785 44653.785 44653.785
28685.632 28685.632 28685.632 28685.632 28685.632
27912.764 27912.764 27912.764 27912.764 27912.764
50977.905 50977.905 50977.905 50977.905 50977.905
81231.890 81231.890 81231.890 81231.890 81231.890
372.418 372.418 372.418 372.418 372.418
25860.739 25860.739 25860.739 25860.739 25860.739
203.856 203.856 203.856 203.856 203.856
61456.722 61456.722 61456.722 61456.722 61456.722
87344.989 87344.989 87344.989 87344.989 87344.989
1468.944 1468.944 1468.944 1468.944 1468.944
9354.191 9354.191 9354.191 9354.191 9354.191
4248.658 4248.658 4248.658 4248.658 4248.658
4657.324 4657.324 4657.324 4657.324 4657.324
14117.866 14117.866 14117.866 14117.866 14117.866
57662.544 57662.544 57662.544 57662.544 57662.544
6784.044 6784.044 6784.044 6784.044 6784.044
37984.081 37984.081 37984.081 37984.081 37984.081
137.558 137.558 137.558 137.558 137.558
111.340 111.340 111.340 111.340 111.340
29230.965 29230.965 29230.965 29230.965 29230.965
74861.551 74861.551 74861.551 74861.551 74861.551
3347.438 3347.438 3347.438 3347.438 3347.438
30749.199 30749.199 30749.199 30749.199 30749.199
16681.520 16681.520 16681.520 16681.520 16681.520
6462.548 6462.548 6462.548 6462.548 6462.548
6219.847 6219.847 6219.847 6219.847 6219.847
37732.809 37732.809 37732.809 37732.809 37732.809
2351.100 2351.100 2351.100 2351.100 2351.100
10277.297 10277.297 10277.297 10277.297 10277.297
76102.142 76102.142 76102.142 76102.142 76102.142
85579.688 85579.688 85579.688 85579.688 85579.688
87186.563 87186.563 87186.563 87186.563 87186.563
92971.385 92971.385 92971.385 92971.385 92971.385
95632.522 95632.522 95632.522 95632.522 95632.522
6317.607 6317.607 6317.607 6317.607 6317.607
34510.372 34510.372 34510.372 34510.372 34510.372
In [ ]:
Content source: balarsen/pymc_learning
Similar notebooks: