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

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.

• Welcome to the site, @JoshHansen. Do you see this Q primarily as a programming question or a statistics question? It looks like the former to me. Note that CV isn't a Python or programming Q&A site (although, Stack Overflow is). If you are wondering about the statistical aspects of expected values, etc, would you edit this to clarify your question? If you mostly want to know how to implement something you understand in Python, flag your Q for moderator attention, & we can migrate it (please don't cross-post). Commented Nov 28, 2012 at 17:00
• It feels like whoever is able to answer this question would need expertise in both the stats and the python details, so I'm not sure where it would be best posted. Where are there more statistically-minded python programmers lurking? There are plenty of scipy questions on this site, and it seems useful for CV to deal not only with the theory but also with the practice of statistics, but I'll defer to others' judgment regarding where to post this question to best get it answered. Commented Nov 28, 2012 at 17:17
• If you have a data set, then you can approximate the expectation by simply using the arithmetic mean $\bar{x}=\dfrac{1}{n}\sum_{j=1}^n x_j$, which converges to the true mean (if it exists) by the LLN. This is a "distribution-free" method, but not necessarily optimal.
– user10525
Commented Nov 28, 2012 at 17:33
• You're right that this question (& others like it--this issue has come up before) lies on the border between the 2 sites. There are scipy questions here, but it seems, only 8, whereas there are 2,406 on SO. My guess, based on my experience, is that you would do better on SO. Most of the users there are programmers, but it's a very large community, & there are many statistically sophisticated users on SO. Nonetheless, you could also leave it here for a while & see if you get a sufficient answer, & flag to migrate it if you don't. (NB, CVers may migrate it anyway.) Commented Nov 28, 2012 at 17:38

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.

• I came up with something like this: def expectation(data, dist=scipy.stats.gamma): return dist(*dist.fit(data)).mean() Commented Nov 29, 2012 at 23:14