In [2]:
%matplotlib inline

In [3]:
from pprint import pprint

import datetime

import pymc
import numpy as np
import spacepy.plot as spp # for the style
import matplotlib.pyplot as plt
import spacepy.toolbox as tb
import spacepy.plot as spp


This unreleased version of SpacePy is not supported by the SpacePy team.
/Users/balarsen/miniconda3/envs/python3/lib/python3.5/site-packages/matplotlib/__init__.py:872: 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 [4]:
print(datetime.datetime.now().isoformat())


2016-09-08T15:24:42.779314

In [10]:
# test data


obsntot = 170
obsbluetot = 78
obsnbkg = 2861
obsbluebkg = 1625
C = 16.61

In [34]:
pymc.Categorical('t1', [0.5, 0.5]).value


Out[34]:
array(0)

In [14]:
nbkg = pymc.Uniform('nbkg', 1, 1e7, value=100)
nclus = pymc.Uniform('nclus', 1,1e7, value=100)
fbkg = pymc.Beta('fbkg', 1,1)
fclus = pymc.Beta('fclus', 1, 1)

pymc.Binomial?

obsnbkg = pymc.Poisson('obsnbkg', nbkg, observed=True, value=obsnbkg)
obsbluebkg = pymc.Binomial('obsbluebkg', fbkg, fbkg/obsnbkg, observed=True, value=obsbluebkg)
obsntot = pymc.Poisson('obsntot', nbkg/C+nclus, observed=True, value=obsntot)
f = (fbkg*nbkg/C + fclus*nclus)/(nbkg/C + nclus)
obsbluetot = pymc.Binomial('obsbluetot', f, obsntot, observed=True, value=obsbluetot)


---------------------------------------------------------------------------
ZeroProbability                           Traceback (most recent call last)
<ipython-input-14-3421a511abb6> in <module>()
      7 
      8 obsnbkg = pymc.Poisson('obsnbkg', nbkg, observed=True, value=obsnbkg)
----> 9 obsbluebkg = pymc.Binomial('obsbluebkg', fbkg, fbkg/obsnbkg, observed=True, value=obsbluebkg)
     10 obsntot = pymc.Poisson('obsntot', nbkg/C+nclus, observed=True, value=obsntot)
     11 f = (fbkg*nbkg/C + fclus*nclus)/(nbkg/C + nclus)

/Users/balarsen/miniconda3/envs/python3/lib/python3.5/site-packages/pymc/distributions.py in __init__(self, *args, **kwds)
    318                     logp_partial_gradients=logp_partial_gradients,
    319                     dtype=dtype,
--> 320                     **arg_dict_out)
    321 
    322     new_class.__name__ = name

/Users/balarsen/miniconda3/envs/python3/lib/python3.5/site-packages/pymc/PyMCObjects.py in __init__(self, logp, doc, name, parents, random, trace, value, dtype, rseed, observed, cache_depth, plot, verbose, isdata, check_logp, logp_partial_gradients)
    773         if check_logp:
    774             # Check initial value
--> 775             if not isinstance(self.logp, float):
    776                 raise ValueError(
    777                     "Stochastic " +

/Users/balarsen/miniconda3/envs/python3/lib/python3.5/site-packages/pymc/PyMCObjects.py in get_logp(self)
    930                     (self._value, self._parents.value))
    931             else:
--> 932                 raise ZeroProbability(self.errmsg)
    933 
    934         return logp

ZeroProbability: Stochastic obsbluebkg's value is outside its support,
 or it forbids its parents' current values.

In [6]:
pymc.Matplot.plot(model)


Plotting s
Plotting bkg

In [14]:
from matplotlib.colors import LogNorm
_ = spp.plt.hist2d(s.trace()[:], bkg.trace()[:], bins=50, normed=True, norm=LogNorm(), vmin=1e-4, vmax=1e-2)
spp.plt.colorbar()


Out[14]:
<matplotlib.colorbar.Colorbar at 0x1215aa128>

In [15]:
spp.plt.hist2d?

In [ ]: