3 The constrained Polya posterior

Jeremy Strief and Glen Meeden

Previous | Next

We begin by recalling a well known approximation to the Polya posterior. If $n/N$ is small then under the Polya posterior, $p=\left(p_1,\dots ,p_n\right)$ has approximately a Dirichlet distribution with a parameter vector of all ones, i.e., it is uniform on the $n-1$ dimensional simplex, where ${\displaystyle {\sum }_{j=1}^n{p}_j=1.}$ It is usually more efficient to generate complete copies of the population using this approximation than the urn model described in the previous section. In addition this approximation will be useful when we consider the constrained Polya posterior, a generalization of the Polya posterior which arises when prior information about auxiliary variables are available to the sampler.

In many problems, in addition to the variable of interest, $y,$ the sampler has in hand auxiliary variables for which prior information is available. A very common case is when the population mean of an auxiliary variable is known. More generally, we will assume that prior information about the population can be expressed by a set of linear equality and inequality constraints on a collection of auxiliary variables.

We assume that in addition to the characteristic of interest $y$ there is a set of auxiliary variables $x^1\mathrm{,}{x}^2\mathrm{,}\dots \mathrm{,}{x}^m.$ For unit $i$ let

$\left(y_i,x_i\right)=\left(y_i,x_i^1,x_i^2,\dots x_i^m\right)$

be the vector of values for $y$ and the auxiliary variables. We suppose that for any unit in the sample this vector of values is observed. We assume the prior information about the population can be expressed through a set of linear equality and inequality constraints on the population values of the auxiliary variables. For the set of possible values for a given auxiliary variable the coefficients defining a constraint will correspond to the proportions of units in the population taking on these values. We now illustrate this more precisely by explaining how we translate this prior information about the population to the observed sample values. Given a sample this will allow us to construct simulated copies of the population consistent with the prior information.

Given a sample $s,$ for $i=1,2,\dots ,n,$ let $\left(y_i,x_i\right)$ be the observed values which, for simplicity, we assume are distinct. Let $p_i$ be the proportion of units which are assigned the value $\left(y_i,x_i\right)$ in a simulated complete copy of the population. Any linear constraint on the population value of an auxiliary variable translates in an obvious way to a linear constraint on these observed values. For example, if the population mean of $x^1$ is known to be less than or equal to some value, say $b_1,$then for the simulated population this translates to the constraint

${\displaystyle \sum _{i=1}^n{p}_i{x}_i^1\le b_1.}$

If the population median of $x^2$ is known to be equal to $b_2$then for the simulated population this becomes the constraint

${\displaystyle \sum _{i=1}^n{p}_i{u}_i=0.5}$

where $u_i=1$ if $x_i^2\le b_2$ and it is zero otherwise. Hence, given a collection of population constraints based on prior information and a sample we will be able to represent the corresponding constraints on a simulated value of $p$ by two systems of equations

$A_{1,s}p=b_1$ (3.1)

$A_{\mathrm{2,}s}p\le b_2$ (3.2)

where $A_{\mathrm{1,}s}$ and $A_{\mathrm{2,}s}$ are $m_1\times n$ and $m_2\times n$ matrices and $b_1$ and $b_2$are vectors of the appropriate dimensions.

Let $P$ denote the subset of the $n$ dimensional simplex which is defined by equations (3.1) and (3.2). We assume the sample is such that $P$ is non-empty and hence it is a non-full dimensional polytope. In this case the appropriate approximate version of the Polya posterior should just be the uniform distribution over $P.$We call this distribution the constrained Polya posterior (CPP). If one could generate independent observations from the CPP then one could find approximately the posterior expectation of population parameters of interest and find approximate 0.95 stepwise Bayes credible intervals. Unfortunately we do not know how to do this. Instead, one can use Markov chain Monte Carlo (MCMC) methods to find such estimates approximately. This can done in R (R Development Core Team 2005) and using the R package polypost which is available in CRAN. More details on the CPP and simulating from it are available in Lazar et al. (2008).

Previous | Next