constrppmn: Estimating a Population Mean Using the Constrained Polya...

In polyapost: Simulating from the Polya Posterior

Description

Let p=(p1,...,pn) be a probability distribution defined on ysamp, the set of observed values, in a sample of size n from some population. p is assumed to belong to a polytope which is a lower dimensional subset of the n-dimensional simplex. The polytope is defined by a collection of linear equality and inequality constraints. A dependent sequence of values for p are generated by a Markov chain using the Metropolis-Hastings algorithm whose stationary distribution is the uniform distribution over the polytope. For each generated value of p the corresponding mean, sum(p*y) is found.

Usage

constrppmn(A1,A2,A3,b1,b2,b3,initsol,reps,ysamp,burnin)

Arguments

A1	The matrix for the equality constraints.This must always contain the constraint sum(p) == 1.
A2	The matrix for the <= inequality constraints. This must always contain the constraints -p <= 0.
A3	The matrix for the >= inequality constraints. If there are no such constraints A3 must be set equal to NULL.
b1	The rhs vector for A1, each component must be nonnegative.
b2	The rhs vector for A2, each component must be nonnegative.
b3	The rhs vector for A3, each component must be nonnegative. If A3 is NULL then b3 must be NULL.
initsol	A vector which lies in the interior of the polytope.
reps	The total length of the chain that is generated.
ysamp	The observed sample from the population of interest.
burnin	The point in the chain at which the set of computed means begins.

Value

The returned value is a list whose first component is the chain of the means of length reps - burnin -1, whose second component is the mean of the first component (i.e. the Polya estimate of the population mean) and whose third component is the 2.5th and 97.5th quantiles of the first component (i.e. an approximate 95 percent confidence interval of the population mean).

Examples

A1<-rbind(rep(1,6),1:6)
A2<-rbind(c(2,5,7,1,10,8),diag(-1,6))
b1<-c(1,3.5)
b2<-c(6,rep(0,6))
initsol<-rep(1/6,6)
rep<-1006
burnin<-1000
ysamp<-c(1,2.5,3.5,7,4.5,6)
out<-constrppmn(A1,A2,NULL,b1,b2,NULL,initsol,rep,ysamp,burnin)
out[[1]] # the Markov chain of the means.
out[[2]] # the average of out[[1]]
out[[3]] # the 2.5th and 97.5th quantiles of out[[1]]

polyapost documentation built on Oct. 7, 2021, 5:10 p.m.