Random() method can't be assigned by the decorator?

42 views
Skip to first unread message

Rich

unread,
Nov 27, 2009, 2:00:51 AM11/27/09
to PyMC
Hi all;
I've been playing around with implementing some custom stochastic
elements, and there seems to be an issue with providing a new
stochastic with a 'random' method via the '@stochastic' decorator. I
don't seem to have this problem when implementing 'Stochastic'
directly. I'm not sure if I'm doing something wrong or expecting the
decorator to do something it's not supposed to.

Specifically, if I enter code such as:
<code>
@stochastic(dtype=int)
def s(value=1900, t_l=1851, t_h=1962):
"""The switchpoint for the rate of disaster occurrence."""
def logp(value, t_l, t_h):
if value > t_h or value < t_l:
return -Inf
else:
return -log(t_h - t_l + 1)
def random(t_l, t_h):
return round( (t_l - t_h) * random() ) + t_l
</code>

(copied and pasted from the documentation), load it, and then try to
call s.rand() or s.random(), I get the error "TypeError: random()
takes exactly 2 arguments (0 given)". If I oblige and provide it with
endpoints -- s.random(1851,1962) -- then I get the error: "TypeError:
random() takes exactly 1 argument (3 given)"

I've tried this with my own functions, and I run into exactly the same
issue (actually I ran into it there first, and went back to the
tutorial code to try and figure it out).

If I enter the code directly implementing Stochastic (eg "s =
Stochastic( logp = s_logp, ...." and so on, from the tail end of
section 4.1 in the documentation) then the call s.random() works just
fine and returns a random number in the correct range.

This is all with the latest project source (I think I checked it out
Friday or Saturday) built on Windows 7 for Python 2.6.3 using MinGW.

On a completely different note, is there any interest in having
another example contributed? I put together a simple Bayesian network
example for practice (the horse genetics one from Jensen's 1996 book)
which -- while it is a slightly silly way to do that example -- might
be interesting to someone.

Regards,
Rich Harang

Anand Patil

unread,
Nov 27, 2009, 8:05:51 AM11/27/09
to py...@googlegroups.com
Hi Rich,

On Fri, Nov 27, 2009 at 7:00 AM, Rich <rich....@gmail.com> wrote:
> Specifically, if I enter code such as:
> <code>
> @stochastic(dtype=int)
> def s(value=1900, t_l=1851, t_h=1962):
>    """The switchpoint for the rate of disaster occurrence."""
>    def logp(value, t_l, t_h):
>        if value > t_h or value < t_l:
>            return -Inf
>        else:
>            return -log(t_h - t_l + 1)
>    def random(t_l, t_h):
>        return round( (t_l - t_h) * random() ) + t_l
> </code>
>
> (copied and pasted from the documentation), load it, and then try to
> call s.rand() or s.random(), I get the error "TypeError: random()
> takes exactly 2 arguments (0 given)".  If I oblige and provide it with
> endpoints -- s.random(1851,1962) -- then I get the error: "TypeError:
> random() takes exactly 1 argument (3 given)"

It's a namespace thing, try

@stochastic(dtype=int)
def s(value=1900, t_l=1851, t_h=1962):
"""The switchpoint for the rate of disaster occurrence."""
def logp(value, t_l, t_h):
if value > t_h or value < t_l:
return -Inf
else:
return -log(t_h - t_l + 1)
def random(t_l, t_h):
return round( (t_l - t_h) * numpy.random.random() ) + t_l

> On a completely different note, is there any interest in having
> another example contributed?  I put together a simple Bayesian network
> example for practice (the horse genetics one from Jensen's 1996 book)
> which -- while it is a slightly silly way to do that example -- might
> be interesting to someone.

Absolutely! There's a wiki page,
http://code.google.com/p/pymc/wiki/TutorialsAndRecipes, which hasn't
_exactly_ taken off yet... but your example would be a step in the
right direction!

Anand

Rich

unread,
Nov 28, 2009, 7:05:04 PM11/28/09
to PyMC
On Nov 27, 5:05 am, Anand Patil <anand.prabhakar.pa...@gmail.com>
wrote:
> It's a namespace thing, try
>
> @stochastic(dtype=int)
> def s(value=1900, t_l=1851, t_h=1962):
>    """The switchpoint for the rate of disaster occurrence."""
>    def logp(value, t_l, t_h):
>        if value > t_h or value < t_l:
>            return -Inf
>        else:
>            return -log(t_h - t_l + 1)
>    def random(t_l, t_h):
>        return round( (t_l - t_h) * numpy.random.random() ) + t_l

That did the trick, thank you.

Unfortunately I still have a (minor) issue. If I try to initialize a
custom stochastic without a value, I get one of two new error
messages, depending on whether I use the constructor or the
decorator. The decorator gives me the error: "ValueError: Stochastic
s's value initialized to None; no initial value or random method
provided." while the constructor gives me the error:
"ZeroProbability: Stochastic s's value is outside its support, or it
forbids its parents' current values." I'd assumed that it was related
to my previous issue, but it seems to have persisted.

This isn't a huge deal, since I can set it to some arbitrary value and
let it wash out in the burn-in, but it seems to not be performing as
the documentation suggests it might (although the problem, as before,
may well exist somewhere between the keyboard and the chair).

> Absolutely! There's a wiki page,http://code.google.com/p/pymc/wiki/TutorialsAndRecipes, which hasn't
> _exactly_ taken off yet... but your example would be a step in the
> right direction!

I'd be delighted to contribute it, but it doesn't seem like I can do
anything to the page you've linked but post a comment to it, and the
code is somewhat long for a single comment. I'm not familiar with
Google code, is there some trick I'm missing?

Thanks again;
Rich

Anand Patil

unread,
Dec 2, 2009, 6:30:54 PM12/2/09
to py...@googlegroups.com
On Sun, Nov 29, 2009 at 12:05 AM, Rich <rich....@gmail.com> wrote:
> I'd be delighted to contribute it, but it doesn't seem like I can do
> anything to the page you've linked but post a comment to it, and the
> code is somewhat long for a single comment.  I'm not familiar with
> Google code, is there some trick I'm missing?

Wow, I didn't realize that wasn't editable by anyone. This may explain
why the wiki page hasn't taken off. :) If you have a blog or github
page or somesuch where you'd normally post it, go ahead and post it
there, then post the link here... otherwise just email me the example
directly and I'll put it up.

Anand
Reply all
Reply to author
Forward
0 new messages