How do you compute expected value of arbitrary distributions in scipy?

Question

Suppose I have a python function using scipy that returns the expectation $E\left[ X \right]$ for some data assuming it is gamma distributed:

def expectation(data):
    shape,loc,scale=scipy.stats.gamma.fit(data)
    expected_value = shape * scale
    return expected_value

(My understanding is that scipy's parameterization of the gamma leaves us with $E\left[ X \right] = shape \cdot scale$.) However, I would like to generalize my code so I can drop in different distributions in place of the gamma -- for example, the log-normal distribution. Is there a way to write that code in a general way? In other words, how do I finish this function:

def expectation(data, dist=scipy.stats.gamma):
    ???

I see a few possible approaches:

1. Use the scipy.stats.*.expect method. Thus far I haven't been able to figure out how to use it. How would I parameterize the method given the shape,loc,scale parameters above?

2. Use the mean method of a "frozen" random variable object. In scipy-speak, is "mean" equivalent to $E\left[ X \right]$?

3. Give up on writing general code and just compute $E\left[ X \right]$ directly for each distribution. I don't want to do this if I can avoid it.

Additionally, please address whether under your suggested method I would pay any performance penalty, i.e. because it uses a numerical rather than analytical approach to the integral in evaluating the expectation.

Answer 1

The docs say that .mean() of the frozen random variable is what your are looking for:

>>> help(scipy.stats.norm)
...
|  mean(self, *args, **kwds)
|      Mean of the distribution    
...

I have just tested your second idea and to me it seems as if it works:

>>> rv = scipy.stats.norm(3, 2)
>>> rv.mean()
>>> 3.0

My guess is that it works with arbitrary distributions of scipy.stats.