Online estimation of an unknown gamma distribution. I want to quickly create probability distributions for the two parameters to the gamma distribution (shape k and scale θ) after each sample taken.  How would you go about this?
My goal is to estimate the likelihood that a future sample will fall within a given range or, equivalently, what proportion of the population fall between two values. Once I have an estimate for my two parameters, I will need to determine what area of the parameter space meets my criterion, then double integrate over their joint distribution to produce my likelihood. How would you go about this?
If you see a more straightforward or efficient way to achieve this estimate please tell me as well. Keep in mind that two of the goals are minimization of samples taken and efficiency of calculation.
 A: Supposing that efficiency is a serious concern, you can use the method of moments estimators, which can be easily updated. (See the engineering statistics handbook here.
See also here.) Periodically, though, you can use the maximum likelihood estimates (when you have computational time available).
Suppose we are working with the PDF
$$
f(x\,|\,k,\theta) = \frac{x^{k - 1}}{\Gamma(k)\theta^k}\exp\left( \frac{-x}{\theta} \right)
$$
and have IID $x_i\sim \Gamma[k,\theta]$.
The method of moments estimators are then:
$$
\hat{\theta} = \frac{1}{\bar{x}n}\sum_i (x_i - \bar{x})^2
\;\;\;\;\&\;\;\;\;
\hat{k} = \frac{\bar{x}}{\hat{\theta}}
$$
Now to do online updates. At update $j$, let $\bar{x}_j$ be the sample mean and define $\gamma_j$ as the sum of squared differences:
$$
\gamma_j = \sum_i (x_i - \bar{x}_j)^2
$$
Suppose we have the estimates for data points $x_1,\ldots,x_{n-1}$ (so that we have $n$ points in total now).
Then, we can do the update as:
\begin{align}
\bar{x}_{n} &= \bar{x}_{n-1} + \frac{x_n - \bar{x}_{n-1}}{n} \\[2mm]
\gamma_{n} &= \gamma_{n-1} + (x_n - \bar{x}_{n-1})(x_n - \bar{x}_{n}) \\[2mm]
\hat{\theta}_n &= \frac{\gamma_n}{\bar{x}_n\,n} \\[2mm]
\hat{k}_n &= \frac{\bar{x}_n}{\hat{\theta}_n}
\end{align}
You can then integrate $f(x|\hat{k}_n,\hat{\theta}_n)$ to estimate the probability of a sample falling in a certain interval, for example.
