Masked TruncatedNormal

[mamacneil]

I can't seem to pass a masked_array for mu into the TruncatedNormal; Normal works fine - is this a bug?

E.G.

from pymc import *
from numpy import ma, array

x = array([1,2,-999])
x1 = ma.masked_array(x,mask=x==-999)

y = array([-999,2,4])
y1 = ma.masked_array(y,mask=y==-999)

Normal('Yi',mu=x1,tau=0.01,value=y1,observed=T)

Out[10]: <pymc.distributions.Normal 'Yi' at 0x1100ebfd0>

VS:

TruncatedNormal('Yi',mu=x1,tau=0.01,a=0,b=100,value=y1,observed=T)

"""

In [11]: TruncatedNormal('Yi',mu=x1,tau=0.01,a=0,b=100,value=y1,observed=T)

ZeroProbability Traceback (most recent call last)
/Volumes/ma_macneil/Documents/Current_work/AIMS_survey_impacts/rates_pymc/ in ()
----> 1 TruncatedNormal('Yi',mu=x1,tau=0.01,a=0,b=100,value=y1,observed=T)

/Users/amac/Downloads/pymc-devs-pymc-5b2d924/pymc/distributions.pyc in init(self, _args, *_kwds)
266 random = debug_wrapper(random)
267 else:
--> 268 Stochastic.init(self, logp=logp, random=random, logp_partial_gradients = logp_partial_gradients, dtype=dtype, **arg_dict_out)
269
270 new_class.name = name

/Users/amac/Downloads/pymc-devs-pymc-5b2d924/pymc/PyMCObjects.pyc in init(self, logp, doc, name, parents, random, trace, value, dtype, rseed, observed, cache_depth, plot, verbose, isdata, check_logp, logp_partial_gradients)
700 if check_logp:
701 # Check initial value

--> 702 if not isinstance(self.logp, float):
703 raise ValueError, "Stochastic " + self.name + "'s initial log-probability is %s, should be a float." %self.logp.repr()
704

/Users/amac/Downloads/pymc-devs-pymc-5b2d924/pymc/PyMCObjects.pyc in get_logp(self)
832 raise ZeroProbability, self.errmsg + "\nValue: %s\nParents' values:%s" % (self._value, self._parents.value)
833 else:
--> 834 raise ZeroProbability, self.errmsg
835
836 return logp

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

[mamacneil]

Just updated to the latest git (pymc-devs-pymc-055b94f) with a slightly modified result:

E.G. 2

In [14]: TruncatedNormal('Yi',mu=x1,tau=0.01,a=0,b=100,value=y1,observed=T)

FloatingPointError Traceback (most recent call last)
/Volumes/ma_macneil/Documents/Current_work/AIMS_survey_impacts/rates_pymc/ in ()
----> 1 TruncatedNormal('Yi',mu=x1,tau=0.01,a=0,b=100,value=y1,observed=T)

/Library/Python/2.7/site-packages/pymc/distributions.pyc in init(self, _args, *_kwds)
267 random = debug_wrapper(random)
268 else:
--> 269 Stochastic.init(self, logp=logp, random=random, logp_partial_gradients = logp_partial_gradients, dtype=dtype, **arg_dict_out)
270
271 new_class.name = name

/Library/Python/2.7/site-packages/pymc/PyMCObjects.pyc in init(self, logp, doc, name, parents, random, trace, value, dtype, rseed, observed, cache_depth, plot, verbose, isdata, check_logp, logp_partial_gradients)
700 dtype=dtype,
701 plot=plot,
--> 702 verbose=verbose)
703
704 # self._logp.force_compute()

/Library/Python/2.7/site-packages/pymc/Node.pyc in init(self, doc, name, parents, cache_depth, trace, dtype, plot, verbose)
191 self.extended_children = set()
192
--> 193 Node.init(self, doc, name, parents, cache_depth, verbose=verbose)
194
195 if self.dtype is None:

/Library/Python/2.7/site-packages/pymc/Node.pyc in init(self, doc, name, parents, cache_depth, verbose)
109
110 # Initialize

--> 111 self.parents = parents
112
113 def _get_parents(self):

/Library/Python/2.7/site-packages/pymc/Node.pyc in _set_parents(self, new_parents)
130
131 # Get new lazy function

--> 132 self.gen_lazy_function()
133
134 parents = property(_get_parents, _set_parents, doc="Self's parents: the variables referred to in self's declaration.")

/Library/Python/2.7/site-packages/pymc/PyMCObjects.pyc in gen_lazy_function(self)
740 ultimate_args = self.extended_parents | set([self]),
741 cache_depth = self._cache_depth)
--> 742 self._logp.force_compute()
743
744

/Library/Python/2.7/site-packages/pymc/LazyFunction.so in pymc.LazyFunction.LazyFunction.force_compute()

/Library/Python/2.7/site-packages/pymc/distributions.pyc in wrapper(**kwds)
2957 def wrapper(**kwds):
2958 value = kwds.pop('value')
-> 2959 return f(value, **kwds)
2960
2961 if arguments is None:

/Library/Python/2.7/site-packages/pymc/distributions.pyc in truncated_normal_like(x, mu, tau, a, b)
2486 lPhib = pymc.utils.normcdf((b-mu)/sigma, log=True)
2487 try:
-> 2488 d = utils.log_difference(lPhib, lPhia)
2489 except ValueError:
2490 return -np.inf

/Library/Python/2.7/site-packages/pymc/utils.pyc in log_difference(lx, ly)
732 # Negative log of double-precision infinity

733     li=-709.78271289338397

--> 734 diff = ly - lx
735 # Make sure log-difference can succeed

736     if np.any(diff>=0):

FloatingPointError: invalid value encountered in subtract

[mamacneil]

Ah I see this might be due to behaviour I didn't understand - the x1 values here are taking on the fill_value from ma and not being updated; is this correct?

[fonnesbeck]

Just getting back to this now. Yes, its only the value that should be masked array, not the parameters. So, the following works:

from pymc import *
from numpy import ma, array

x = array([1,2,3])

y = array([-999,2,4])
y1 = ma.masked_array(y,mask=y==-999)

# Normal('Yi',mu=x1,tau=0.01,value=y1,observed=T)
TruncatedNormal('Yi',mu=x,tau=0.01,a=0,b=100,value=y1,observed=T)