[SciPy-User] expectation function for stats.distributions

josef.pktd@gmai... josef.pktd@gmai...
Fri Jul 31 15:10:07 CDT 2009


On Fri, Jul 31, 2009 at 3:10 PM, nicky van foreest<vanforeest@gmail.com> wrote:
> 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?

Thanks for the interest.

No, integrate.quad only work for continuous random variables. But
there are the corresponding generic functions also for the discrete
distributions in stats.distribution.

An expectation function could be based on the generic calculation of
the uncentered moments of the discrete distributions. I fixed the
functions last year, but haven't looked at them in a while. I think
the main point was to decide which and how many points to use in the
calculation, otherwise it is just a vector inner product.


To the name: I'm up for ideas
I don't like integral so much because it is not obvious that it
calculates an expectation and not just the integral of the pdf, as in
stats.kde and my wrapper for the multivariate normal cdf.
Additionally, my function also calculates conditional expectations,
which are not just integrals.

I like expectation, but it is a noun and might refer to the
expectation of the random variable itself, i.e. the mean. That's why I
ended up with a variation on expectation_function.

Josef


>
> Nicky
>
> 2009/7/31  <josef.pktd@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@scipy.org
>> http://mail.scipy.org/mailman/listinfo/scipy-user
>>
> _______________________________________________
> SciPy-User mailing list
> SciPy-User@scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-user
>


More information about the SciPy-User mailing list