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(pi*yi) is found.

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

`A1` |
The matrix for the equality constraints.This must always contain the constraint that the sum of the pi's is one. |

`A2` |
The matrix for the <= inequality constraints. This must always contain the constraints -pi <= 0, i.e. that the pi's must be nonnegative. |

`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. |

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).

1 2 3 4 5 6 7 8 9 10 11 12 | ```
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]]
``` |

Questions? Problems? Suggestions? Tweet to @rdrrHQ or email at ian@mutexlabs.com.

All documentation is copyright its authors; we didn't write any of that.