In [2]:

    

# http://twiecki.github.io/blog/2013/09/12/bayesian-glms-1/


    

In [1]:

    

import pymc3 as pm

import numpy as np
import matplotlib.pyplot as plt
import spacepy.plot as spp


    

    




/Users/balarsen/miniconda3/envs/python3/lib/python3.6/site-packages/matplotlib/__init__.py:913: 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 [12]:

    

def PA(amp, angles, order=2):
    return amp*np.sin(np.deg2rad(angles))**order

fig, ax = plt.subplots(ncols=3, nrows=3, figsize=(10,10), sharex=True, sharey=True)


ang = np.linspace(0, 180)
anp=10
order=np.linspace(0.8, 5, 9)
ii = 0 
for row in range(3):
    for col in range(3):
        ax[row][col].plot(ang, PA(amp, ang, order[ii]))
        ax[row][col].set_title('{:0.3}'.format(order[ii]))
        ii += 1


    

In [83]:

    

# make some data that follows sin2

angles = np.linspace(5, 175, 7)
amp = 10
order = 2
# add some possion noise to each
np.random.seed(8675309)
nangles = np.random.normal(loc=angles, scale=7.0, size=len(angles))
vals = PA(amp, nangles, order)

plt.plot(np.linspace(0, 180, 100), PA(amp, np.linspace(0, 180, 100), order))
plt.plot(angles, vals, 'x')


    

    
Out[83]:





[<matplotlib.lines.Line2D at 0x129d4a358>]






    

In [84]:

    

noise, add


    

    
Out[84]:





(array([ 0,  4,  4, 10,  3,  0]), array([ 0, -1,  0, -1, -1,  0]))

In [85]:

    

data = {'vals':vals, 'angles':angles, 'x':angles}
with pm.Model() as mdl_fish_alt:
    pm.glm.GLM.from_formula('vals ~ x', data)
    trace = pm.sample(3000, njobs=4) # draw 3000 posterior samples using NUTS sampling


    

    




Auto-assigning NUTS sampler...
Initializing NUTS using ADVI...
Average Loss = 32.193:   9%|▉         | 18868/200000 [00:02<00:25, 7042.14it/s]
Convergence archived at 18900
Interrupted at 18,900 [9%]: Average Loss = 3,623
 99%|█████████▉| 3482/3500 [00:10<00:00, 340.57it/s]/Users/balarsen/miniconda3/envs/python3/lib/python3.6/site-packages/pymc3/step_methods/hmc/nuts.py:456: UserWarning: Chain 0 contains 1 diverging samples after tuning. If increasing `target_accept` does not help try to reparameterize.
  % (self._chain_id, n_diverging))
100%|██████████| 3500/3500 [00:10<00:00, 333.47it/s]
/Users/balarsen/miniconda3/envs/python3/lib/python3.6/site-packages/pymc3/step_methods/hmc/nuts.py:456: UserWarning: Chain 1 contains 1 diverging samples after tuning. If increasing `target_accept` does not help try to reparameterize.
  % (self._chain_id, n_diverging))

In [86]:

    

pm.traceplot(trace, combined=True)


    

    
Out[86]:





array([[<matplotlib.axes._subplots.AxesSubplot object at 0x12a5085f8>,
        <matplotlib.axes._subplots.AxesSubplot object at 0x12b4ffe10>],
       [<matplotlib.axes._subplots.AxesSubplot object at 0x129a78748>,
        <matplotlib.axes._subplots.AxesSubplot object at 0x12bbe5b38>],
       [<matplotlib.axes._subplots.AxesSubplot object at 0x12b354978>,
        <matplotlib.axes._subplots.AxesSubplot object at 0x12e12ee10>]], dtype=object)






    

In [87]:

    

plt.figure(figsize=(7, 7))
plt.plot(np.linspace(0, 180, 100), PA(amp, np.linspace(0, 180, 100), order))
plt.plot(angles, vals, 'x')

pm.plot_posterior_predictive_glm(trace, samples=100, eval=np.linspace(0,180,90),
                              label='posterior predictive regression lines', )
# plt.plot(x, true_regression_line, label='true regression line', lw=3., c='y')


    

In [88]:

    

data = {'vals':vals, 'angles':angles, 'x':angles}
with pm.Model() as mdl_fish_alt:
    sigma = pm.HalfCauchy('sigma', beta=10, testval=1.)
    x_coeff = pm.Uniform('amp', 1, 100, )
    order = pm.Uniform('order', 0.5, 8)
    
    # Define likelihood
    likelihood = pm.Normal('y', mu=x_coeff * np.sin(np.deg2rad(data['x']))**order,
                        sd=sigma, observed=data['vals'])

    trace = pm.sample(3000, njobs=4, tune=5000) # draw 3000 posterior samples using NUTS sampling

    
    
pm.traceplot(trace, combined=True)


    

    




Auto-assigning NUTS sampler...
Initializing NUTS using ADVI...
Average Loss = 21.01:   6%|▌         | 12261/200000 [00:01<00:23, 7886.79it/s] 
Convergence archived at 12600
Interrupted at 12,600 [6%]: Average Loss = 469.16
 99%|█████████▉| 7920/8000 [00:21<00:00, 359.80it/s]/Users/balarsen/miniconda3/envs/python3/lib/python3.6/site-packages/pymc3/step_methods/hmc/nuts.py:456: UserWarning: Chain 1 contains 5 diverging samples after tuning. If increasing `target_accept` does not help try to reparameterize.
  % (self._chain_id, n_diverging))
100%|█████████▉| 7998/8000 [00:22<00:00, 373.82it/s]/Users/balarsen/miniconda3/envs/python3/lib/python3.6/site-packages/pymc3/step_methods/hmc/nuts.py:456: UserWarning: Chain 0 contains 6 diverging samples after tuning. If increasing `target_accept` does not help try to reparameterize.
  % (self._chain_id, n_diverging))
100%|██████████| 8000/8000 [00:22<00:00, 363.09it/s]
/Users/balarsen/miniconda3/envs/python3/lib/python3.6/site-packages/pymc3/step_methods/hmc/nuts.py:456: UserWarning: Chain 2 contains 9 diverging samples after tuning. If increasing `target_accept` does not help try to reparameterize.
  % (self._chain_id, n_diverging))
/Users/balarsen/miniconda3/envs/python3/lib/python3.6/site-packages/pymc3/step_methods/hmc/nuts.py:456: UserWarning: Chain 3 contains 2 diverging samples after tuning. If increasing `target_accept` does not help try to reparameterize.
  % (self._chain_id, n_diverging))






    




---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-88-e202e53a8b06> in <module>()
     16 
     17 plt.figure(figsize=(7, 7))
---> 18 plt.plot(np.linspace(0, 180, 100), PA(amp, np.linspace(0, 180, 100), order))
     19 plt.plot(angles, vals, 'x')
     20 

~/miniconda3/envs/python3/lib/python3.6/site-packages/matplotlib/pyplot.py in plot(*args, **kwargs)
   3315                       mplDeprecation)
   3316     try:
-> 3317         ret = ax.plot(*args, **kwargs)
   3318     finally:
   3319         ax._hold = washold

~/miniconda3/envs/python3/lib/python3.6/site-packages/matplotlib/__init__.py in inner(ax, *args, **kwargs)
   1896                     warnings.warn(msg % (label_namer, func.__name__),
   1897                                   RuntimeWarning, stacklevel=2)
-> 1898             return func(ax, *args, **kwargs)
   1899         pre_doc = inner.__doc__
   1900         if pre_doc is None:

~/miniconda3/envs/python3/lib/python3.6/site-packages/matplotlib/axes/_axes.py in plot(self, *args, **kwargs)
   1404         kwargs = cbook.normalize_kwargs(kwargs, _alias_map)
   1405 
-> 1406         for line in self._get_lines(*args, **kwargs):
   1407             self.add_line(line)
   1408             lines.append(line)

~/miniconda3/envs/python3/lib/python3.6/site-packages/matplotlib/axes/_base.py in _grab_next_args(self, *args, **kwargs)
    405                 return
    406             if len(remaining) <= 3:
--> 407                 for seg in self._plot_args(remaining, kwargs):
    408                     yield seg
    409                 return

~/miniconda3/envs/python3/lib/python3.6/site-packages/matplotlib/axes/_base.py in _plot_args(self, tup, kwargs)
    379         if len(tup) == 2:
    380             x = _check_1d(tup[0])
--> 381             y = _check_1d(tup[-1])
    382         else:
    383             x, y = index_of(tup[-1])

~/miniconda3/envs/python3/lib/python3.6/site-packages/matplotlib/cbook.py in _check_1d(x)
   2218     dimension; leaves everything else untouched.
   2219     '''
-> 2220     if not hasattr(x, 'shape') or len(x.shape) < 1:
   2221         return np.atleast_1d(x)
   2222     else:

TypeError: object of type 'TensorVariable' has no len()





    













    

In [89]:

    

pm.traceplot(trace[::2], combined=True)


    

    
Out[89]:





array([[<matplotlib.axes._subplots.AxesSubplot object at 0x10b11a400>,
        <matplotlib.axes._subplots.AxesSubplot object at 0x12939c7f0>],
       [<matplotlib.axes._subplots.AxesSubplot object at 0x124fbf208>,
        <matplotlib.axes._subplots.AxesSubplot object at 0x12e1aada0>],
       [<matplotlib.axes._subplots.AxesSubplot object at 0x12b081b38>,
        <matplotlib.axes._subplots.AxesSubplot object at 0x12cffcb38>]], dtype=object)






    

In [90]:

    

# grab 100 random samples
sami = np.random.random_integers?


    

In [ ]:

    

sami = np.random.random_integers


    

In [91]:

    

sami = np.random.randint(0, 3000*4, 200)


    

In [92]:

    

data


    

    
Out[92]:





{'angles': array([   5.        ,   33.33333333,   61.66666667,   90.        ,
         118.33333333,  146.66666667,  175.        ]),
 'vals': array([ 0.25140412,  3.86932124,  6.46708935,  9.95379338,  8.48258157,
         3.20993372,  0.18964863]),
 'x': array([   5.        ,   33.33333333,   61.66666667,   90.        ,
         118.33333333,  146.66666667,  175.        ])}

In [96]:

    

plt.figure(figsize=(10,10))

ampp = trace['amp'][sami]
angp = np.linspace(0,180,90)
orderp = trace['order'][sami]

for aa, oo in zip(ampp, orderp):
    plt.plot(angp, PA(aa, angp, oo), 'grey', alpha=0.5)

plt.errorbar(data['angles'], data['vals'], marker='o', c='r', markersize=5, yerr=np.sqrt(data['vals']))


    

    
Out[96]:





<Container object of 3 artists>






    

In [94]:

    

plt.errorbar?


    

In [ ]:

    




    

In [97]:

    

data = {'vals':vals, 'angles':angles, 'x':angles}
with pm.Model() as mdl_fish_alt:
    sigma = pm.HalfCauchy('sigma', beta=10, testval=1.)
    x_coeff = pm.Uniform('amp', 1, 100, )
    order = pm.Uniform('order', 0.5, 8)
    offset = pm.Normal('offset', 0.0, sd=5)
    
    # Define likelihood
    likelihood = pm.Normal('y', mu=x_coeff * np.sin(np.deg2rad(data['x']+np.deg2rad(offset)))**order,
                        sd=sigma, observed=data['vals'])

    trace = pm.sample(3000, njobs=4, tune=5000) # draw 3000 posterior samples using NUTS sampling

    
    
pm.traceplot(trace, combined=True)


    

    




Auto-assigning NUTS sampler...
Initializing NUTS using ADVI...
Average Loss = 20.284:   8%|▊         | 15576/200000 [00:02<00:27, 6756.26it/s]
Convergence archived at 16200
Interrupted at 16,200 [8%]: Average Loss = 369.92
 99%|█████████▉| 7912/8000 [00:34<00:00, 238.25it/s]/Users/balarsen/miniconda3/envs/python3/lib/python3.6/site-packages/pymc3/step_methods/hmc/nuts.py:456: UserWarning: Chain 1 contains 7 diverging samples after tuning. If increasing `target_accept` does not help try to reparameterize.
  % (self._chain_id, n_diverging))
100%|█████████▉| 7985/8000 [00:35<00:00, 229.94it/s]/Users/balarsen/miniconda3/envs/python3/lib/python3.6/site-packages/pymc3/step_methods/hmc/nuts.py:456: UserWarning: Chain 0 contains 3 diverging samples after tuning. If increasing `target_accept` does not help try to reparameterize.
  % (self._chain_id, n_diverging))
100%|██████████| 8000/8000 [00:35<00:00, 227.35it/s]
/Users/balarsen/miniconda3/envs/python3/lib/python3.6/site-packages/pymc3/step_methods/hmc/nuts.py:456: UserWarning: Chain 3 contains 2 diverging samples after tuning. If increasing `target_accept` does not help try to reparameterize.
  % (self._chain_id, n_diverging))






    
Out[97]:





array([[<matplotlib.axes._subplots.AxesSubplot object at 0x12c3534a8>,
        <matplotlib.axes._subplots.AxesSubplot object at 0x12c7cc390>],
       [<matplotlib.axes._subplots.AxesSubplot object at 0x12bfa4550>,
        <matplotlib.axes._subplots.AxesSubplot object at 0x12d202fd0>],
       [<matplotlib.axes._subplots.AxesSubplot object at 0x1257473c8>,
        <matplotlib.axes._subplots.AxesSubplot object at 0x12c856978>],
       [<matplotlib.axes._subplots.AxesSubplot object at 0x128cddc88>,
        <matplotlib.axes._subplots.AxesSubplot object at 0x130431cf8>]], dtype=object)






    

In [118]:

    

plt.figure(figsize=(10,10))

ampp = trace['amp'][sami]
angp = np.linspace(0,180,90)
orderp = trace['order'][sami]
offsetp = trace['offset'][sami]

for aa, oo, of in zip(ampp, orderp, offsetp):
    plt.plot(angp, PA(aa, angp+of, oo), 'grey', alpha=0.2)

plt.scatter(data['angles'], data['vals'], marker='o', c='r',  )


    

    
Out[118]:





<matplotlib.collections.PathCollection at 0x137c98eb8>






    

In [109]:

    

plt.hist2d(angp, PA(ampp[5:], angp+offsetp[5:], orderp[5:]), bins=20)


    

    




---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-109-453bc74d7adc> in <module>()
----> 1 plt.hist2d(angp, PA(ampp[5:], angp+offsetp[5:], orderp[5:]), bins=20)

ValueError: operands could not be broadcast together with shapes (90,) (195,) 

In [108]:

    

aa


    

    
Out[108]:





()

In [ ]:

    




    

In [ ]:

    




    

In [ ]:

    




    

In [ ]: