[SciPy-User] expectation function for stats.distributions

josef.pktd@gmai... josef.pktd@gmai...
Fri Jul 31 11:01:02 CDT 2009

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

copy of file below


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

        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.

        expected value : float
    if fn is None:
        def fun(x, *args):
            return x*self.pdf(x, *args)
        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)
        invfac = 1.0
    return integrate.quad(fun, lb, ub,

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

More information about the SciPy-User mailing list