[SciPy-User] expectation function for stats.distributions

[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))