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