# Fast integration of a posterior distribution

I wish to infer the posterior distribution on the probability of success $\theta$ in some binomial process, the twist being that I know that $\theta$ lies in the interval [0.5, 1].

The trouble is that a Beta distribution that is supported on [0.5, 1] is no longer conjugate. In particular, the prior looks like this ($p$ is the pdf):

$p(\theta) \propto {(\theta-0.5)}^{\alpha-1} {(1-\theta)}^{\beta-1}$

The likelihood, given $x$ successes in $n$ trials, is the usual binomial likelihood:

$p(n,x | \theta) \propto \theta^x {(1-\theta)}^{n-x}$

The posterior, as it stands, looks like this:

$p( \theta | n, x) \propto \theta^x {(\theta-0.5)}^{\alpha-1} {(1-\theta)}^{n-x+\beta-1}$

I would like to be able to make CDF calculations (i.e. compute integrals) of the posterior efficiently. (numerical integration, e.g., may be too slow for my purposes). I am unable to figure out how to do this.

Much appreciate any help the forum can provide.

• It sounds like the only reason you use beta is because its conjugate? uniform over the interval is fine if you have a reasonable amount of data ($\alpha=\beta=1$) - then you just have an incomplete beta function to evaluate, most software packages can do this for you. Plus numerical integration is very quick because you have a short interval $[0.5,1]$ - uniform grid sampling would work well. – probabilityislogic Apr 18 '11 at 10:26
• The indefinite integral can be computed exactly in terms of hypergeometric functions and gamma functions. With integral values of $\alpha$, $\beta$, $n$, and $x$, this reduces to a polynomial in the limit of integration. If you don't need to perform general analysis of the expressions, though, then numerical quadrature as suggested by @probability will be fast and effective. It should be adaptive sampling, not uniform, though: the singularities at 1/2 and 1 for small values of $\alpha$ and $\beta$ might otherwise cause havoc. – whuber Apr 18 '11 at 15:43
• @probabilityislogic: if I can, I'd prefer to use something other than a uniform prior, since my prior is quite informative, and for speed, I'd like to observe as few actual binomial trials as possible. Beta seemed the natural alternative. But your answer induced me to consider a simple piece-wise uniform prior, which may have enough information for my purposes. – Venu Satuluri Apr 18 '11 at 22:04

How accurate does your posterior cdf need to be? You might consider replacing the continuous prior with a discrete approximation:

$p^*(\theta) \propto p(\theta) 1(\theta\in t_1, \dots, t_k)$

where $p(\theta)$ is your original continuous prior.

Then to compute the posterior you just calculate likelihood x prior

$p(\theta|x) \propto p^*(\theta)p(x|\theta)$

over the support of the prior $t_1, \dots, t_k$ and renormalize.

This is called "griddy Gibbs" by some. It can be quite effective if you have an informative prior in which case you can choose the grid points non-uniformly (and, of course, if you can live with a discrete approximation coarse enough to be computationally feasible).

• Do you have any suggestions for how to choose the grid points (non-uniformly)? If I need to choose k grid points, and say I restrict the grid points to only n possibilities, then the brute-force approach is to evaluate the likelihood of each (n choose k) combination, and choose the max. likelihood one. Are there alternatives? (Sorry for the late comment, I just saw your answer today.) – Venu Satuluri Sep 15 '11 at 18:24

There may be a simpler approach, simply by applying the usual Beta conjugate to the binomial, and then requiring $\theta \in [\frac12,1]$. You can do this with an indicator function, for example as in

$$p(\theta) \propto \theta^{\alpha-1} {(1-\theta)}^{\beta-1} \mathbb{1}[{\tfrac12 \le \theta \le 1}]$$

Now apply your $p(n,x | \theta) \propto \theta^x {(1-\theta)}^{n-x}$ to get the posterior density

$$p(\theta|x) \propto \theta^{\alpha+x-1} {(1-\theta)}^{\beta+n-x-1} \mathbb{1}[{\tfrac12 \le \theta \le 1}].$$

The posterior cumulative distribution function for $\theta \in [\frac12,1]$ is then $\dfrac{I_\theta(\alpha+x, \beta+n-x)-I_{\frac12}(\alpha+x, \beta+n-x)}{1-I_{\frac12}(\alpha+x, \beta+n-x)}$ with $I$ representing a regularised incomplete beta function or the cumulative distribution function of a Beta distribution, which any decent statistical program will calculate quickly, such as R's pbeta function.

You can always use Monte Carlo Integration or the Midpoint method. With Monte Carlo, you simply generate a bunch of points in your parameter space and see if they are in the area or volume or hyper-dimensional space you are trying to integrate.

From: http://farside.ph.utexas.edu/teaching/329/lectures/node109.html "Let us now consider the so-called Monte-Carlo method for evaluating multi-dimensional integrals. Consider, for example, the evaluation of the area, , enclosed by a curve, . Suppose that the curve lies wholly within some simple domain of area , as illustrated in Fig. 97. Let us generate points which are randomly distributed throughout . Suppose that of these points lie within curve . Our estimate for the area enclosed by the curve is simply" the ratio of random points in the space times the size of the space. The link has a nice picture and a description of the inferior midpoint method that I would suggest skipping.

• Because this is a 1D, not multidimensional, integral, numerical methods (as suggested by @probabilityislogic) tend to be much faster and more accurate than Monte Carlo methods. – whuber Apr 18 '11 at 15:37