pymc-devs / pymc2

THIS IS THE **OLD** PYMC PROJECT (VERSION 2). PLEASE USE PYMC INSTEAD:
http://pymc-devs.github.com/pymc/
Other
879 stars 229 forks source link

rv_frozen fails argcheck(..) because it thinks the shape argument is a list #117

Closed rogermosher closed 8 years ago

rogermosher commented 8 years ago

I am wrapping scipy's genextreme distribution, and setting priors for shape (normal), loc (normal), and scale (uniform). I then set up my observations, where gev is the PYMC wrapped version of genextreme:

# create observations
obs = gev("obs", c=c, loc=loc, scale=scale, value=data, observed=True)

It complains of a type error, because you can't compare a list with an int. Following the code with the debugger I see that at some point the "c" argument (the shape), has been separated from the other args and buried in a list within a list. In the rv_frozen init method it is passed to argcheck:

self.dist._argcheck(*shapes)

which removes one level of list, but the 'c' prior arrives inside argcheck buried in a list of length 1. So of course argcheck fails:

        self.b = where(c > 0, 1.0 / max(c, _XMIN), inf)

It thinks c is a list.

Not sure what is going wrong here. For more detail, I show my code at stackoverflow:

Hoping it was my error :)

rogermosher commented 8 years ago

Somehow the link to stackoverflow got lost, here is my code:

def genExtremeParams(data,oparams): """ genExtremeParams: use Bayesian techniques to estimate genExtremeParams from data :param data: an ndarray used to find the parameter estimates :param oparams: scipy fitted params (shape,loc,scale] :return: [(50%, 2.5%, and 97.5% quantiles) for shape, location, scale] """

# create stochastic variable
gev = pmscipy.stochastic_from_scipy_dist(ss.genextreme)

# create priors...
c     = pm.Normal("shape", oparams[0], 0.25)  #shape
loc   = pm.Normal("loc", oparams[1], 0.25)
scale = pm.Uniform("scale", 0,2.0)

# create observations
obs = gev("obs", c=c, loc=loc, scale=scale, value=data, observed=True)

# create model
mcmc = pm.MCMC(set([obs, c, loc, scale]))
mcmc.sample(40000,20000,10)

# generate estimates
shTrace = mcmc.trace("shape")       #shape
shQ = ssm.mquantiles(shTrace, [0.5, 0.025, 0.975], axis=0)

locTrace = mcmc.trace("loc")
locQ = ssm.mquantiles(locTrace, [0.5, 0.025, 0.975], axis=0)

sclTrace = mcmc.trace("scale")
sclQ = ssm.mquantiles(sclTrace, [0.5, 0.025, 0.975], axis=0)

# return quantiles (shape, loc, scale)
return shQ, locQ, sclQ
fonnesbeck commented 8 years ago

Those errors appear to be related to SciPy code, rather than PyMC. I'm not familiar enough with it to provide support, so you might post an issue with the SciPy project. However, jf this issue is with lists, you might consider using NumPy arrays throughout.

rogermosher commented 8 years ago

Thanks I think will update to PYMC3 and go on from there.