constrppprob {polyapost}

Dependent Sampling from the Uniform Distribution on a Polytope.

Description

Let p=(p_1,\ldots,p_n) be a probability distribution which belongs to a lower dimensional polytope of the n-dimensional simplex. The polytope is defined by a collection of linear equality and inequality constraints. A dependent sequence of the p's are generated by a Markov chain using the Metropolis-Hastings algorithm whose stationary distribution is the uniform distribution over the polytope. This is done by generating k blocks of size step where the last member of each is returned.

Usage

constrppprob(A1,A2,A3,b1,b2,b3,initsol,step,k)

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.

step

The number of p's generated in a block before the last member of a block is returned.

k

The total number of blocks generated and hence the number of p's returned.

Value

The returned value is a k by n matrix of probability vectors.

Examples

A1<-rbind(rep(1,6),1:6)
A2<-rbind(c(2,5,7,1,10,8),diag(-1,6))
A3<-matrix(c(1,1,1,0,0,0),1,6)
b1<-c(1,3.5)
b2<-c(6,rep(0,6))
b3<-0.45
initsol<-rep(1/6,6)
constrppprob(A1,A2,A3,b1,b2,b3,initsol,2000,5)