[SciPy-dev] arithmetic on distributions

Robert Kern rkern at ucsd.edu
Wed Aug 10 17:40:59 CDT 2005

Perry, Alexander (GE Infrastructure) wrote:
> Has anyone implemented, or considered implementing, the arithmetic methods for rv_continuous?
> First suggestion.  When calling an object to freeze it, if the loc is already a distribution, we should merge its variance into the variance of the newly frozen distribution:
>>ten = norm ( 10, 3 )
>>print norm ( ten, 4 ).stats()
> (10.0, 25.0)
> Currently, an error is triggered.  Because of the difficulty of merging the standard deviation of the old value into the new desired standard deviation, I suspect that most people who encounter the error will take the lazy way out and simply discard the variance attached to the value:
>>try: value=value.stats()[0]
>>except: pass # the above is wrong!

I can't imagine why. I think most people would simply conclude that
passing in a frozen distribution doesn't work (especially since it is
certainly not documented to do so).

> Second suggestion.  When a binary operator such as __add__() is called on an object that derives from rv_continuous, currently an error is generated.  There should be a function that checks whether the other parameter is a simple number and, if it is, modifies its own loc and scale.  If both parameters of the binary operator are instances of rv_continuous, the function returns an instance of the new subclass rv_composite that simply stores those two parameters and the desired operator.  The new subclass passes through all operations that can be directly decomposed (such as rvs and stats).  This allows a real world distribution to be assembled as a single object, and that object directly supports both linear approximation and monte-carlo methods:
>>fifteen = norm ( 10, 3 ) + uniform( 3, 4 )
>># Linearized estimate
>>print "%7.5f" % norm ( fifteen, 0 ).cdf ( 20 )
> 0.94008
>># Monte carlo estimate
>>print sum ( [ x<20 for x in fifteen.rvs(100000) ])
> 94125
> The composite distribution instance implements the remainder of the standard "rv" features through integration and differentiation of the underlying distributions.  This is no different than what the users would have to do themselves, but saves the effort of explicitly calling those scipy features.
>># Accurate integrated value
>>print "%7.5f" % fifteen.cdf(20)
> 0.94000

I would suggest writing rv_composite first and implementing the
appropriate algorithms and shortcuts (like norm(0, 3) + norm(0, 4) =>
norm(0, 5)). The general case for most operations is going to entail
sampling which is going to require some more input from the user (e.g.
the requested accuracy of the estimated parameters). Once the algorithms
are implemented, and we know what information needs to flow where, then
is the right time to think about implementing a nicer syntax.

> Third suggestion.  The __cmp__ method should exist and can simply compare the mean since its python definition is that it returns a Boolean upon which flow of control can be decided.  The more accurate comparison functions (such as __lt__) should exist and return probabilities because their python definition is that they return a data result with no flow control interpretation.

If the rich comparisons (__lt__ and __eq__, at least) are defined, then
__cmp__ will never be called.

Robert Kern
rkern at ucsd.edu

"In the fields of hell where the grass grows high
  Are the graves of dreams allowed to die."
   -- Richard Harter

More information about the Scipy-dev mailing list