[pymc] Development version breaks stochastic_from_dist() for non-numerical Stochastics

[MaxyB]

One reason I love pymc (over JAGS, BUGS, etc.) is that it allows me to
easily define "exotic" models with non-numerical parameters. In this
case, I have a model parameter that is a combinatorial object (a weak
order, or "stratified hierarchy"). With properly defined Stochastic
classes and step methods, posterior sampling of such objects works
beautifully. In pymc 2.0rc2, at least.

I recently upgraded to the development version, and was disappointed
to find that my custom Stochastic classes no longer work. Below is the
code for a simple example. This works fine in 2.0rc2, but not in svn.
In 2.0rc2, I can instantiate the custom Stochastic like so:
su = StrataUniform('test', ['a','b','c','d'])
Then su.random() generates uniformly random weak orders of the
elements ['a', 'b', 'c', 'd'], and su.logp returns the uniform
density. Beautiful.

But in svn, as soon as I try to instantiate the StrataUniform class, I
get an exception:
---------------------------------------------------------------------------
ValueError Traceback (most recent call
last)

/Users/max/projects/typfreq/src/<ipython console> in <module>()

/Library/Frameworks/Python.framework/Versions/2.5/lib/python2.5/site-
packages/pymc/PyMCObjects.pyc in random(***failed resolving
arguments***)
672 if self._random:
673 # Get current values of parents for use as
arguments for _random()
--> 674 r = self._random(**self.parents.value)
675 else:
676 raise AttributeError, 'Stochastic '+self.__name__
+' does not know how to draw its value, see documentation'

/Library/Frameworks/Python.framework/Versions/2.5/lib/python2.5/site-
packages/pymc/distributions.pyc in newfun(*args, **kwargs)
80 def bind_size(randfun, shape):
81 def newfun(*args, **kwargs):
---> 82 return np.reshape(randfun(size=shape, *args,
**kwargs),shape)
83 newfun.scalar_version = randfun
84 return newfun

/Library/Frameworks/Python.framework/Versions/2.5/lib/python2.5/site-
packages/numpy/core/fromnumeric.pyc in reshape(***failed resolving
arguments***)
149 reshape = a.reshape
150 except AttributeError:
--> 151 return _wrapit(a, 'reshape', newshape, order=order)
152 return reshape(newshape, order=order)
153

/Library/Frameworks/Python.framework/Versions/2.5/lib/python2.5/site-
packages/numpy/core/fromnumeric.pyc in _wrapit(***failed resolving
arguments***)
35 except AttributeError:
36 wrap = None
---> 37 result = getattr(asarray(obj),method)(*args, **kwds)
38 if wrap:
39 if not isinstance(result, mu.ndarray):

ValueError: total size of new array must be unchanged

It seems to me that the svn code has some sort of new fanciness that
assumes the stochastic's value is a numerical array. Is there a good
way for me to work around this? Would it be possible for the pymc
developers to relax this assumption; seems like it shouldn't be hard
to infer that the Stochastic generates a non-numerical, non-array
value, and refrain from trying to reshape it.

from pymc import *
import numpy as np
import random

class Strata(object):
"""
A stratified hierarchy (weak order) over the given constraint set
("con").
"""
def __init__(self, con, strataList=None, membership=None):
self.con = list(con)
self.strata = strataList or []
if membership is not None:
d = {}
for i,j in enumerate(membership):
if j not in d:
d[j] = set()
d[j].add(con[i])
for s in xrange(len(self.con)):
if s in d:
self.strata.append(d[s])
elif not strataList:
self.strata = [set(con)]

def __str__(self):
return ' >> '.join([str(s) for s in self.strata])

def __cmp__(self, other):
return cmp((self.con, self.strata), (other.con, other.strata))

def __hash__(self):
hashableStrata = tuple([frozenset(s) for s in self.strata])
return hash((tuple(self.con), hashableStrata))

def membership(self):
mem = np.zeros(len(self.con))
for stratNum, strat in enumerate(self.strata):
for c in strat:
cNum = self.con.index(c)
mem[cNum] = stratNum
return mem


def logp_StrataUniform(value=None, con=None):
"""
Uniform likelihood of, and generation of (via random_StrataUniform
below),
stratified hierarchies. NB uniform means most generated
stratified
hierarchies will be complex.
"""
if not value:
raise ValueError('Requested likelihood of null strata')
if not con:
raise ValueError('Requested likelihood of strata with null
con')
k = len(con)
return -np.log(k**k)

def random_StrataUniform(con, **kwargs):
if not con:
raise ValueError('Requested random strata with null con')
k = len(con)
membership = [random.choice(range(k)) for c in con]
return Strata(con, membership=membership)

StrataUniform = stochastic_from_dist('StrataUniform',
logp_StrataUniform,
random_StrataUniform, dtype=object)

--
Max Bane
max....@gmail.com

--
You received this message because you are subscribed to the Google Groups "PyMC" group.
To post to this group, send email to py...@googlegroups.com.
To unsubscribe from this group, send email to pymc+uns...@googlegroups.com.
For more options, visit this group at http://groups.google.com/group/pymc?hl=en.

[Chris Fonnesbeck]

Hi Max,

Sorry for the late reply. Have a look at the latest update (commit
96287af7 in GitHub); I think I have generalized the bind_size wrapper
that was causing the problem, without breaking anything else. Have a
look and see if it solves the issue.

Cheers,
Chris

> max.b...@gmail.com