version 1.1

DESCRIPTION

@@ -0,0 +1,11 @@
+Package: polyapost
+Version: 1.1
+Date: 2006-03-27
+Imports: boot
+Title: Simulating from the Polya posterior
+Author:Glen Meeden <glen@stat.umn.edu> and  Radu Lazar <lazar@stat.umn.edu> 
+Maintainer: Glen Meeden <glen@stat.umn.edu>
+Description: Generate dependent samples from a non-full dimensional polytope
+          via a Markov Chain sampler
+License: GPL version 2 or newer
+Packaged: Wed Mar 29 11:12:01 2006; glen

NAMESPACE

@@ -0,0 +1,11 @@
+useDynLib(polyapost)
+
+importFrom(boot, simplex)
+
+export(polyap)
+export(wtpolyap)
+export(feasible)
+export(constrppmn)
+export(constrppprob)
+
+

R/checkconstr.R

@@ -0,0 +1,20 @@
+#checks if the matrices and vectors defining the
+#constaints are set up properly
+
+checkconstr<-function(A1,A2,A3,b1,b2,b3)
+  {
+    if( !is.matrix(A1)) stop("A1 is not a matrix")
+    if( !is.vector(b1)) stop("b1 is not a vector")
+    if( min(b1)< 0) stop("b1 has negative values")
+    if( nrow(A1) != length(b1)) stop("nrow(A1) != length(b1)")
+    if( !is.matrix(A2)) stop("A2 is not a matrix")
+    if( !is.vector(b2)) stop("b2 is not a vector")
+    if( min(b2)< 0) stop( "b2 has negative values")
+    if( nrow(A2) != length(b2)) stop("nrow(A2) != length(b2)")
+    if( !is.null(A3)){
+      if( !is.matrix(A3)) stop("A3 is not a matrix")
+      if( !is.vector(b3)) stop("b3 is not a vector")
+      if( min(b3)< 0) stop("b3 has negative values")
+      if( nrow(A3) != length(b3)) stop("nrow(A3) != length(b3)")
+    }
+  }

R/constrmat.R

@@ -0,0 +1,45 @@
+#this function transforms the constraints A1x=b1,A2x<=b2,A3%*x>=b3,x>=0 into #A1x=b1,A4x<=b4
+constrmat<-function(A1,A2=NULL,A3=NULL,b2=NULL,b3=NULL)
+{
+  if (is.matrix(A1)){
+      m1 <- nrow(A1)
+      q<-ncol(A1)
+  }
+  else{
+      m1 <- 1
+      m1<-length(A1)
+  }
+  u<-rep(0,q)
+  U<-diag(-1,q)
+  if (!is.null(A2))
+      if (is.matrix(A2))
+          m2 <- nrow(A2)
+      else m2 <- 1
+  else m2 <- 0
+  if (!is.null(A3))
+      if (is.matrix(A3))
+          m3 <- nrow(A3)
+      else m3 <- 1
+  else m3 <- 0
+  if(m2+m3==0){
+     A4<-U
+     b4<-u
+  }
+  else {
+     if (m2==0){
+       A4<-rbind(-A3,U)
+       b4<-c(-b3,u)
+     }
+     else {
+        if(m3==0){
+           A4<-rbind(A2,U)
+	   b4<-c(b2,u)
+        }
+	else{
+	  A4<-rbind(A2,-A3,U)
+	  b4<-c(b2,-b3,u)
+	}
+     }
+  }
+  list(A1,A4,b4)
+}

R/constrppmn.R

@@ -0,0 +1,21 @@
+#For the subset of the simplex defined by
+# A1 p = b1, A2 p <= b2 and A3 p >= b3
+# where the Ai's are matrices and the bi's
+#vectors of nonnegative real numbers this
+#function uses the Metroplis-Hastings algorithm
+
+constrppmn<-function(A1,A2,A3,b1,b2,b3,initsol,reps,ysamp,burnin)
+  {
+    checkconstr(A1,A2,A3,b1,b2,b3)
+    if(!is.null(A3)) {
+      A4<-rbind(A2,-A3)
+      b4<-c(b2,-b3)
+    }
+     else {
+       A4<-A2
+       b4<-b2
+     }
+     out<-polyaest(A1,A4,b4,initsol,reps,ysamp,burnin)
+     return(out)
+  }
+

R/constrppprob.R

@@ -0,0 +1,19 @@
+#Generates k probability vectors in "steps" of
+#size step.
+
+constrppprob<-function(A1,A2,A3,b1,b2,b3,initsol,step,k)
+  {
+    checkconstr(A1,A2,A3,b1,b2,b3)
+    if(!is.null(A3)) {
+      A4<-rbind(A2,-A3)
+      b4<-c(b2,-b3)
+    }
+     else {
+       A4<-A2
+       b4<-b2
+     }
+    nrows<-nrow(A4)
+    out<-probvect(A1,A4,nrow(A4),b4,initsol,step,k)
+    return(out)
+  }
+

R/cwpolya.R

@@ -0,0 +1,18 @@
+
+cwpolya <- function(samp, wts, k) {
+	if (missing(wts)) wts <- rep(1, length(samp))
+	if (length(samp) != length(wts))
+		stop("length(samp) != length(wts)")
+	if (any(wts < 0) || all(wts <= 0))
+		stop("wts are foobar")
+	k <- as.integer(k)
+	if (k <= 0)
+		stop("k must be positive integer")
+       out<-.C("cwpolya",
+       x=as.double(c(samp, rep(0, k))),
+       w=as.double(cumsum(wts)),
+       n=as.integer(length(samp)),
+       k=as.integer(k),
+       PACKAGE="polyapost")
+       return(out$x)
+}

R/feasible.R

@@ -0,0 +1,54 @@
+#this function finds a feasible solution of A1x=b1,A2x<=b2,A3%*x>=b3,x>=eps or returns a negative vector when there is no feasible solution
+#A1 is m1xq,A2 is m2xq,A3 is m3xq, b1 is m1x1 , b2 is m2x1,b3 is m3xq, eps is a scalar between 0 and 1 and has to be bigger than min(b3)
+
+feasible<-function(A1,A2=NULL,A3=NULL,b1,b2=NULL,b3=NULL,eps)
+{
+  if (is.matrix(A1)){
+      m1 <- nrow(A1)
+      q<-ncol(A1)
+  }
+  else {
+      m1 <- 1
+      q<-length(A1)
+  }
+  if (!is.null(A2))
+      if (is.matrix(A2))
+          m2 <- nrow(A2)
+      else m2 <- 1
+  else m2 <- 0
+  if (!is.null(A3))
+      if (is.matrix(A3))
+          m3 <- nrow(A3)
+      else m3 <- 1
+  else m3 <- 0
+
+  if(m2+m3==0){
+     A<-A1
+     b<-b1
+     out<-feasible1(A,b,rep(eps,q))
+  }
+  else {
+     if (m2==0){
+       A<-cbind(rbind(A1,A3),rbind(matrix(0,m1,m3),diag(-1,m3)))
+       b<-c(b1,b3)
+       out<-feasible1(A,b,rep(eps,q+m3))[1:q]
+     }
+     else {
+        if(m3==0){
+	  A<-cbind(rbind(A1,A2),rbind(matrix(0,m1,m2),diag(1,m2)))
+	  b<-c(b1,b2)
+	  out<-feasible1(A,b,rep(eps,q+m2))[1:q]
+
+	}
+	else{
+	  A<-cbind(rbind(A1,A2,A3),rbind(matrix(0,m1,m2+m3),cbind(diag(1,m2),
+	     matrix(0,m2,m3)),cbind(matrix(0,m3,m2),diag(-1,m3))))
+          b<-c(b1,b2,b3)
+        out<-feasible1(A,b,rep(eps,q+m2+m3))[1:q]
+	}
+     }
+  }
+  return(out)
+}
+
+

R/feasible1.R

@@ -0,0 +1,26 @@
+#feasible1 returns a feasible solution for Ax=b, x>=eps, where eps is a small positive number and a negative vector if there is no feasible solution 
+
+feasible1<-function(A, b, eps)
+{
+  if(! is.matrix(A)) stop("A not a matrix")
+  n<-nrow(A)
+  p<-ncol(A)
+  if(n != length(b)) stop(" no. of rows does not match length of second parameter ")
+  if(p != length(eps)) stop("no. of columns  does not match length of third parameter")
+  rside<-b-A%*%eps
+  s<-0
+  for(i in 1:n){
+    if(rside[i] >= 0) s<-s+1
+  }
+  if(s == n){
+    d<-diag(rep(1,n))
+    linsol<-simplex(a=c(rep(0,p), rep(1,n)), A3=cbind(A,d), b3=rside, maxi=FALSE)
+    if(linsol$value<1e-7){
+      x0<-as.vector(linsol$soln[1:p])
+      x<-x0+eps
+    }
+    else  x<--1+eps
+  }
+  else x<--1+eps
+  return(x)
+}

R/means.R

@@ -0,0 +1,35 @@
+#description of the function
+#generates points from A1x=b1,A2x<=b2 and then gets the innerproduct between them
+#and ysamp
+# P=B%*t(B), B spans N(A1)
+#matsize=ncol(A)=ncol(p)=nrow(p)
+#nrow: no. of rows of A2
+#initsol: the initial feasible solution
+#rep: number of points we want to generate
+#ysamp the sample from the population y
+means<-function(P, matsize, A2, nrow, b2, initsol, rep, ysamp)
+{
+ if(! is.matrix(P)) stop("P not a matrix")
+ if(! is.matrix(A2)) stop("A2 not a matrix")
+ if(ncol(P)!=ncol(A2)) {stop("the no. of rows of the matrices are not equal")}
+ if(ncol(P)!=matsize) {stop("the dimension of the matrix that spans the null space is
+ not equal to the second parameter")}
+ if(nrow(A2)!=nrow) {stop("wrong no. of rows")}
+ if(nrow!=length(b2)) {stop("the no. of rows of the constraint matrix does not match
+the length of the rhs vector")}
+ if(matsize!=length(ysamp)){stop("the size of the sample does not match the no. of
+columns of the constraint matrix")}
+ foo<-.C("means",
+        P=as.double(P),
+        matsize=as.integer(matsize),
+	A2=as.double(A2),
+	nrow=as.integer(nrow),
+	b2=as.double(b2),
+        initsol=as.double(initsol),
+        rep=as.integer(rep),
+	ysamp=as.double(ysamp),
+        estimate=double(rep),
+	PACKAGE="polyapost")
+ return(foo$estimate)
+}
+

R/nullspace.R

@@ -0,0 +1,16 @@
+#finds the matrix that spans the null space of a known matrix with the no. of columns >= the no. of rows
+nullspace<-function(A)
+{
+  if(! is.matrix(A)) stop("A not matrix")
+  n<-nrow(A)
+  p<-ncol(A)
+  if(n >= p) stop("no. of rows greater or equal to the no. of columns")
+  s<-svd(A, nu=n, nv=p)
+  k<-0
+  for(i in 1:n) {
+     if (s$d[i] >1e-6) k<-k+1
+  }
+  v2<-NULL
+  for(i in (k+1):p) v2<-cbind(v2,s$v[,i])
+  return(v2)
+}

R/polyaest.R

@@ -0,0 +1,13 @@
+polyaest<-function(A1, A2, b2, initsol, rep, ysamp, burnin)
+{
+  if(! is.matrix(A1)) stop("A1 not matrix")
+  if(! is.matrix(A2)) stop("A2 not matrix")
+  if(nrow(A2)!=length(b2)) {stop(" the no. of rows of the constraint matrix does not match the length of the rhs vector ")}
+  P<-nullspace(A1)%*%t(nullspace(A1))
+  chain<-means(P, ncol(P), A2, nrow(A2), b2, initsol, rep, ysamp)
+  chain<-chain[burnin:rep]
+  estimate<-mean(chain)
+  quantiles<-quantile(chain, c(.025,.975))
+  return(list(chain, estimate, quantiles))
+}
+

R/polyap.R

@@ -0,0 +1,11 @@
+#Given a sample of y values, say ysamp,  of size n
+#from a population of size N this function uses the
+#Polya posterior to  generate one completed copy of
+#the entire population. With k=N-n it returns a vector
+#of length N with ysamp being the first n elements.
+
+polyap<-function(ysamp,k)
+  {
+    out<-cwpolya(ysamp,rep(1,length(ysamp)),k)
+    return(out)
+  }

R/probvect.R

@@ -0,0 +1,16 @@
+#description of the function probvect
+#It generates as many points as we want from the polytope A1x=b1, A2x<=b2. It returns #every step-th point.
+#length of chain=n*step
+probvect<-function(A1, A2, nrow, b2, initsol, step, n)
+{
+  if(! is.matrix(A1)) stop("A1 not a matrix")
+  if(! is.matrix(A2)) stop("A2 not a matrix")
+  if(nrow(A2)!=nrow) {stop("wrong no. of rows")}
+  P<-nullspace(A1)%*%t(nullspace(A1))
+  mat<-NULL
+  for(i in 1:n){
+    initsol<-probvect1(P, ncol(P), A2, nrow, b2, initsol, step)
+    mat<-rbind(mat, initsol)
+  }
+  return(mat)
+}

R/probvect1.R

@@ -0,0 +1,29 @@
+#description of the function probvect1
+#generates points from A1x=b1,A2x<=b2 and returns the last one
+# P=B%*t(B), B spans N(A1)
+#matsize=ncol(A)=ncol(p)=nrow(p)
+#nrow: no. of rows of A2
+#initsol the initial feasible solution
+#length the length of the Markov chain
+
+probvect1<-function(P, matsize, A2, nrow, b2, initsol, length)
+{
+ if(! is.matrix(P)) stop("P not a matrix")
+ if(! is.matrix(A2)) stop("A2 not a matrix")
+ if(ncol(P)!=ncol(A2)) {stop("the no. of rows of the matrices are not equal")}
+ if(ncol(P)!=matsize) {stop("the dimension of the matrix that spans the null space is not equal to the second parameter")}
+ if(nrow(A2)!=nrow) {stop("wrong no. of rows")}
+ if(nrow!=length(b2)) {stop(" the no. of rows of the constraint matrix does not match the length of the rhs vector ")}
+ foo<-.C("probvect1",
+        P=as.double(P),
+        matsize=as.integer(matsize),
+	A2=as.double(A2),
+	nrow=as.integer(nrow),
+	b2=as.double(b2),
+        initsol=as.double(initsol),
+        length=as.integer(length),
+        estimate=double(matsize),
+	PACKAGE="polyapost")
+ return(foo$estimate)
+}
+

R/wtpolyap.R

@@ -0,0 +1,14 @@
+#Given a vector of sampled  y values, ysamp, of size n
+#from a population of size N this function uses the
+#weighted Polya posterior to  generate one completed
+#copy of the entire population. Each member of the sample
+#is assigned a nonnegative weight. These weights are given
+#in the vector wts which is of length n. With k=N-n it returns
+#a vector of length N with ysamp being the first n elements.
+
+wtpolyap<-function(ysamp,wts,k)
+  {
+    out<-cwpolya(ysamp,wts,k)
+    return(out)
+  }
+

inst/doc/Makefile

@@ -0,0 +1,12 @@
+all: pp1.pdf
+pp1.pdf: pp1.tex pp1.bbl pp1.aux
+	pdflatex pp1
+pp1.tex: pp1.Rnw
+	echo 'Sweave("pp1.Rnw")' | R --no-save --no-restore 
+pp1.bbl: pp1.tex ref.bib pp1.aux
+	bibtex pp1
+	latex pp1
+pp1.aux: pp1.tex 
+	latex pp1
+clean: 
+	rm -f pp1.tex pp1.aux pp1.log pp1.dvi pp1-* Rplots.ps pp1.blg