[SciPy-User] expectation function for stats.distributions

[nicky van foreest]

Hi Josef,

I would like such a function. As a matter of fact, just today I needed
the expection of min(X,k), where X is a Poisson RV and k a number. Of
course, this is not particularly difficult, but a function that does
the work for me is better yet. BTW: will integrate.quad also work for
discrete distributions?

Nicky

2009/7/31  <josef.pktd at gmail.com>:
> I just needed to verify an expectation for fixing of the stats.models code.
>
> A while ago a wrote a function that can be attached as method to the
> distributions classes. It's pretty simple, but comes in handy when an
> expectation or conditional expectation has to be checked.
>
> Is there an interest in adding this as a new method to the distributions?
>
> explanations are here
> http://jpktd.blogspot.com/2009/04/having-fun-with-expectations.html
>
> copy of file below
>
> Josef
>
> '''
> copy from webpage, where is the original?
>
> Author jpktd
>
> '''
>
>
> import numpy as np
> from scipy import stats, integrate
>
> def expectedfunc(self, fn=None, args=(), lb=None, ub=None, conditional=False):
>    '''calculate expected value of a function with respect to the distribution
>
>    only for standard version of distribution,
>    location and scale not tested
>
>    Parameters
>    ----------
>        all parameters are keyword parameters
>        fn : function (default: identity mapping)
>           Function for which integral is calculated. Takes only one argument.
>        args : tuple
>           argument (parameters) of the distribution
>        lb, ub : numbers
>           lower and upper bound for integration, default is set to the support
>           of the distribution
>        conditional : boolean (False)
>           If true then the integral is corrected by the conditional probability
>           of the integration interval. The return value is the expectation
>           of the function, conditional on being in the given interval.
>
>    Returns
>    -------
>        expected value : float
>    '''
>    if fn is None:
>        def fun(x, *args):
>            return x*self.pdf(x, *args)
>    else:
>        def fun(x, *args):
>            return fn(x)*self.pdf(x, *args)
>    if lb is None:
>        lb = self.a
>    if ub is None:
>        ub = self.b
>    if conditional:
>        invfac = self.sf(lb,*args) - self.sf(ub,*args)
>    else:
>        invfac = 1.0
>    return integrate.quad(fun, lb, ub,
>                                args=args)[0]/invfac
>
>
> print expectedfunc(stats.norm, lambda(x): (x)**4)
> print expectedfunc(stats.norm, lambda(x): min((x)**2,1.5**2))
>
> #Jonathans version
> from scipy.stats import norm as Gaussian
> c = 1.5 # our "cutoff" point
> c = 0.5 # try another value
> tmp = 2 * Gaussian.cdf(c) - 1
> gamma = tmp + c**2 * (1 - tmp) - 2 * c * Gaussian.pdf(c)
> print gamma
>
> print expectedfunc(stats.norm, lambda(x): min((x)**2,c**2))
> _______________________________________________
> SciPy-User mailing list
> SciPy-User at scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-user
>