[r-cran-eco] 02/30: Import Upstream version 2.2-1
Andreas Tille
tille at debian.org
Thu Sep 7 07:20:57 UTC 2017
This is an automated email from the git hooks/post-receive script.
tille pushed a commit to branch master
in repository r-cran-eco.
commit bb8a214ea9b9287fd1f4334017faa7ed83c26eb3
Author: Andreas Tille <tille at debian.org>
Date: Thu Sep 7 08:59:18 2017 +0200
Import Upstream version 2.2-1
---
DESCRIPTION | 18 +
NAMESPACE | 34 ++
R/checkdata.R | 51 ++
R/coef.eco.R | 9 +
R/coef.ecoNP.R | 14 +
R/cov.eco.R | 24 +
R/eco.R | 134 ++++++
R/ecoBD.R | 117 +++++
R/ecoNP.R | 153 ++++++
R/logit.R | 5 +
R/onLoad.R | 7 +
R/predict.eco.R | 34 ++
R/predict.ecoNP.R | 32 ++
R/predict.ecoNPX.R | 42 ++
R/predict.ecoX.R | 43 ++
R/print.eco.R | 19 +
R/print.ecoBD.R | 23 +
R/print.summary.eco.R | 27 ++
R/print.summary.ecoNP.R | 33 ++
R/print.summary.predict.eco.R | 10 +
R/summary.eco.R | 66 +++
R/summary.ecoNP.R | 86 ++++
R/summary.predict.eco.R | 20 +
R/varcov.R | 71 +++
data/census.txt | 1041 +++++++++++++++++++++++++++++++++++++++++
data/reg.txt | 276 +++++++++++
eco.pdf | Bin 0 -> 192714 bytes
man/census.Rd | 39 ++
man/eco.Rd | 212 +++++++++
man/ecoBD.Rd | 127 +++++
man/ecoNP.Rd | 222 +++++++++
man/predict.eco.Rd | 78 +++
man/predict.ecoNP.Rd | 81 ++++
man/reg.Rd | 33 ++
man/summary.eco.Rd | 72 +++
man/summary.ecoNP.Rd | 71 +++
src/Makevars | 2 +
src/bayes.c | 73 +++
src/bayes.h | 10 +
src/gibbsBase.c | 273 +++++++++++
src/gibbsDP.c | 417 +++++++++++++++++
src/gibbsXBase.c | 305 ++++++++++++
src/gibbsXDP.c | 477 +++++++++++++++++++
src/preBaseX.c | 84 ++++
src/preDP.c | 83 ++++
src/preDPX.c | 85 ++++
src/rand.c | 184 ++++++++
src/rand.h | 13 +
src/sample.c | 237 ++++++++++
src/sample.h | 16 +
src/subroutines.c | 122 +++++
src/subroutines.h | 11 +
src/vector.c | 103 ++++
src/vector.h | 22 +
54 files changed, 5841 insertions(+)
diff --git a/DESCRIPTION b/DESCRIPTION
new file mode 100644
index 0000000..6604407
--- /dev/null
+++ b/DESCRIPTION
@@ -0,0 +1,18 @@
+Package: eco
+Version: 2.2-1
+Date: 2005-09-28
+Title: R Package for Fitting Bayesian Models of Ecological Inference in 2x2 Tables
+Author: Kosuke Imai <kimai at princeton.edu>,
+ Ying Lu <ylu at iq.harvard.edu>.
+Maintainer: Kosuke Imai <kimai at princeton.edu>
+Depends: R (>= 2.0), MASS
+Description: eco is a publicly available R package that fits the
+ parametric and nonparametric Bayesian models for
+ ecological inference in 2x2 tables. The models are fit
+ using the Markov chain Monte Carlo algorithms that are
+ described in Imai and Lu (2004, 2005).
+LazyLoad: yes
+LazyData: yes
+License: GPL (version 2 or later)
+URL: http://imai.princeton.edu/research/eco.html
+Packaged: Thu Sep 29 20:44:14 2005; kimai
diff --git a/NAMESPACE b/NAMESPACE
new file mode 100644
index 0000000..63640e8
--- /dev/null
+++ b/NAMESPACE
@@ -0,0 +1,34 @@
+useDynLib(eco)
+
+importFrom(MASS, mvrnorm)
+
+export(
+ eco,
+ ecoBD,
+ ecoNP,
+ summary.eco,
+ summary.ecoNP,
+ print.summary.eco,
+ print.summary.ecoNP,
+ predict.eco,
+ predict.ecoX,
+ predict.ecoNP,
+ predict.ecoNPX)
+
+S3method(varcov, eco)
+S3method(varcov, ecoNP)
+S3method(coef, eco)
+S3method(coef, ecoNP)
+S3method(predict, eco)
+S3method(predict, ecoX)
+S3method(predict, ecoNP)
+S3method(predict, ecoNPX)
+S3method(summary, eco)
+S3method(summary, ecoNP)
+S3method(summary, predict.eco)
+S3method(print, eco)
+S3method(print, ecoBD)
+S3method(print, summary.eco)
+S3method(print, summary.ecoNP)
+S3method(print, summary.predict.eco)
+
diff --git a/R/checkdata.R b/R/checkdata.R
new file mode 100644
index 0000000..683f070
--- /dev/null
+++ b/R/checkdata.R
@@ -0,0 +1,51 @@
+checkdata <- function(X,Y, supplement, ndim) {
+ # check and reorganize inputs
+ if (any(X<0) || any(X>1) || any(Y<0) || any(Y>1))
+ stop("Values of X and Y have to be between 0 and 1.")
+ ind <- 1:length(X)
+ res <- list()
+ res$X1type <- res$samp.X1 <- res$X1.W1 <- 0
+ res$X0type <- res$samp.X0 <- res$X0.W2 <- 0
+
+ ## X = 1
+ X1.ind <- ind[along=(X==1)]
+ if (length(X[X!=1])<length(X)){
+ res$X1type <- 1
+ res$samp.X1 <- length(X1.ind)
+ res$X1.W1 <- Y[X1.ind]
+ }
+
+ ## X = 0
+ X0.ind <- ind[along=(X==0)]
+ if (length(X[X!=0])<length(X)){
+ res$X0type <- 1
+ res$samp.X0 <- length(X0.ind)
+ res$X0.W2 <- Y[X0.ind]
+ }
+
+ XX.ind <- setdiff(ind, union(X0.ind, X1.ind))
+ res$X.use <- X[XX.ind]
+ res$Y.use <- Y[XX.ind]
+
+ res$order.old<-order(c(XX.ind, X0.ind, X1.ind))
+ res$n.samp <- length(res$Y.use)
+ res$d <- cbind(res$X.use, res$Y.use)
+
+ ## check survey data
+
+ if (any(supplement <0) || any(supplement >1))
+ stop("survey data have to be between 0 and 1.")
+ if ((dim(supplement)[2] != ndim) && (length(supplement)>0))
+ stop("when context=TRUE, use n by 3 otherwise use n by 2 matrix\
+ for survey data, when context=TRUE")
+ if (is.null(supplement))
+ res$survey.samp <- res$survey.data <- res$survey.yes <- 0
+ else {
+ res$survey.samp <- length(supplement[,1])
+ res$survey.data <- as.matrix(supplement)
+ res$survey.yes <- 1
+ }
+
+ return(res)
+
+}
diff --git a/R/coef.eco.R b/R/coef.eco.R
new file mode 100644
index 0000000..cf62409
--- /dev/null
+++ b/R/coef.eco.R
@@ -0,0 +1,9 @@
+coef.eco <- function(object, subset = NULL, ...) {
+ mu <- object$mu
+ if (is.null(subset))
+ subset <- 1:nrow(mu)
+ else if (max(subset) > nrow(mu))
+ stop(paste("invalid input for `subset.' only", nrow(mu), "draws are stored."))
+
+ return(mu[subset,])
+}
diff --git a/R/coef.ecoNP.R b/R/coef.ecoNP.R
new file mode 100644
index 0000000..0957f6b
--- /dev/null
+++ b/R/coef.ecoNP.R
@@ -0,0 +1,14 @@
+coef.ecoNP <- function(object, subset = NULL, obs = NULL, ...) {
+ mu <- object$mu
+ if (is.null(subset))
+ subset <- 1:nrow(mu)
+ else if (max(subset) > nrow(mu))
+ stop(paste("invalid input for `subset.' only", nrow(mu), "draws are stored."))
+
+ if (is.null(obs))
+ obs <- 1:dim(object$mu)[3]
+ else if (max(subset) > dim(object$mu)[3])
+ stop(paste("invalid input for `obs.' only", dim(object$mu)[3], "draws are stored."))
+
+ return(mu[subset,,obs])
+}
diff --git a/R/cov.eco.R b/R/cov.eco.R
new file mode 100644
index 0000000..d59bab6
--- /dev/null
+++ b/R/cov.eco.R
@@ -0,0 +1,24 @@
+cov.eco <- function(object, subset = NULL, ...) {
+ if (is.null(subset))
+ subset <- 1:nrow(object$mu)
+ else if (max(subset) > nrow(object$mu))
+ stop(paste("invalid input for `subset.' only", nrow(mu), "draws are stored."))
+
+ p <- ncol(object$mu)
+ n <- length(subset)
+ Sigma <- array(0, c(p, p, n))
+ cov <- object$Sigma
+ for (i in 1:n) {
+ count <- 1
+ for (j in 1:p) {
+ Sigma[j,j:p,i] <- cov[subset[i],count:(count+p-j)]
+ count <- count + p - j + 1
+ }
+ diag(Sigma[,,i]) <- diag(Sigma[,,i]/2)
+ Sigma[,,i] <- Sigma[,,i] + t(Sigma[,,i])
+ }
+ if (n > 1)
+ return(Sigma)
+ else
+ return(Sigma[,,1])
+}
diff --git a/R/eco.R b/R/eco.R
new file mode 100644
index 0000000..851a99d
--- /dev/null
+++ b/R/eco.R
@@ -0,0 +1,134 @@
+eco <- function(formula, data = parent.frame(), N = NULL, supplement = NULL,
+ context = FALSE, mu0 = 0, tau0 = 2, nu0 = 4, S0 = 10,
+ mu.start = 0, Sigma.start = 10, parameter = TRUE,
+ grid = FALSE, n.draws = 5000, burnin = 0, thin = 0,
+ verbose = FALSE){
+
+ ## contextual effects
+ if (context)
+ ndim <- 3
+ else
+ ndim <- 2
+
+ ## checking inputs
+ if (burnin >= n.draws)
+ stop("n.draws should be larger than burnin")
+ if (length(mu0)==1)
+ mu0 <- rep(mu0, ndim)
+ else if (length(mu0)!=ndim)
+ stop("invalid inputs for mu0")
+ if (is.matrix(S0)) {
+ if (any(dim(S0)!=ndim))
+ stop("invalid inputs for S0")
+ }
+ else
+ S0 <- diag(S0, ndim)
+ if (length(mu.start)==1)
+ mu.start <- rep(mu.start, ndim)
+ else if (length(mu.start)!=ndim)
+ stop("invalid inputs for mu.start")
+ if (is.matrix(Sigma.start)) {
+ if (any(dim(Sigma.start)!=ndim))
+ stop("invalid inputs for Sigma.start")
+ }
+ else
+ Sigma.start <- diag(Sigma.start, ndim)
+
+ ## getting X, Y, and N
+ mf <- match.call()
+ tt <- terms(formula)
+ attr(tt, "intercept") <- 0
+ if (is.matrix(eval.parent(mf$data)))
+ data <- as.data.frame(data)
+ X <- model.matrix(tt, data)
+ Y <- model.response(model.frame(tt, data = data))
+ N <- eval(mf$N, data)
+
+ # check data and modify inputs
+ tmp <- checkdata(X,Y, supplement, ndim)
+ bdd <- ecoBD(formula=formula, data=data)
+
+ ## fitting the model
+ n.store <- floor((n.draws-burnin)/(thin+1))
+ unit.par <- 1
+ unit.w <- tmp$n.samp+tmp$samp.X1+tmp$samp.X0
+ n.w <- n.store * unit.w
+
+ if (context)
+ res <- .C("cBaseecoX", as.double(tmp$d), as.integer(tmp$n.samp),
+ as.integer(n.draws), as.integer(burnin), as.integer(thin+1),
+ as.integer(verbose), as.integer(nu0), as.double(tau0),
+ as.double(mu0), as.double(S0), as.double(mu.start),
+ as.double(Sigma.start), as.integer(tmp$survey.yes),
+ as.integer(tmp$survey.samp), as.double(tmp$survey.data),
+ as.integer(tmp$X1type), as.integer(tmp$samp.X1),
+ as.double(tmp$X1.W1), as.integer(tmp$X0type),
+ as.integer(tmp$samp.X0), as.double(tmp$X0.W2),
+ as.double(bdd$Wmin[,1,1]), as.double(bdd$Wmax[,1,1]),
+ as.integer(parameter), as.integer(grid),
+ pdSMu0 = double(n.store), pdSMu1 = double(n.store), pdSMu2 = double(n.store),
+ pdSSig00=double(n.store), pdSSig01=double(n.store), pdSSig02=double(n.store),
+ pdSSig11=double(n.store), pdSSig12=double(n.store), pdSSig22=double(n.store),
+ pdSW1=double(n.w), pdSW2=double(n.w), PACKAGE="eco")
+ else
+ res <- .C("cBaseeco", as.double(tmp$d), as.integer(tmp$n.samp),
+ as.integer(n.draws), as.integer(burnin), as.integer(thin+1),
+ as.integer(verbose), as.integer(nu0), as.double(tau0),
+ as.double(mu0), as.double(S0), as.double(mu.start),
+ as.double(Sigma.start), as.integer(tmp$survey.yes),
+ as.integer(tmp$survey.samp), as.double(tmp$survey.data),
+ as.integer(tmp$X1type), as.integer(tmp$samp.X1),
+ as.double(tmp$X1.W1), as.integer(tmp$X0type),
+ as.integer(tmp$samp.X0), as.double(tmp$X0.W2),
+ as.double(bdd$Wmin[,1,1]), as.double(bdd$Wmax[,1,1]),
+ as.integer(parameter), as.integer(grid),
+ pdSMu0=double(n.store), pdSMu1=double(n.store),
+ pdSSig00=double(n.store),
+ pdSSig01=double(n.store), pdSSig11=double(n.store),
+ pdSW1=double(n.w), pdSW2=double(n.w),
+ PACKAGE="eco")
+
+ W1.post <- matrix(res$pdSW1, n.store, unit.w, byrow=TRUE)[,tmp$order.old]
+ W2.post <- matrix(res$pdSW2, n.store, unit.w, byrow=TRUE)[,tmp$order.old]
+ W <- array(rbind(W1.post, W2.post), c(n.store, 2, unit.w))
+ colnames(W) <- c("W1", "W2")
+ res.out <- list(call = mf, X = X, Y = Y, N = N, W = W,
+ Wmin=bdd$Wmin[,1,], Wmax = bdd$Wmax[,1,],
+ burin = burnin, thin = thin, nu0 = nu0,
+ tau0 = tau0, mu0 = mu0, S0 = S0)
+
+ if (parameter)
+ if (context) {
+ res.out$mu <- cbind(matrix(res$pdSMu0, n.store, unit.par, byrow=TRUE),
+ matrix(res$pdSMu1, n.store, unit.par, byrow=TRUE),
+ matrix(res$pdSMu2, n.store, unit.par, byrow=TRUE))
+ colnames(res.out$mu) <- c("mu1", "mu2", "mu3")
+ res.out$Sigma <- cbind(matrix(res$pdSSig00, n.store, unit.par, byrow=TRUE),
+ matrix(res$pdSSig01, n.store, unit.par, byrow=TRUE),
+ matrix(res$pdSSig02, n.store, unit.par, byrow=TRUE),
+ matrix(res$pdSSig11, n.store, unit.par, byrow=TRUE),
+ matrix(res$pdSSig12, n.store, unit.par, byrow=TRUE),
+ matrix(res$pdSSig22, n.store, unit.par, byrow=TRUE))
+ colnames(res.out$Sigma) <- c("Sigma11", "Sigma12", "Sigma13",
+ "Sigma22", "Sigma23", "Sigma33")
+ }
+ else {
+ res.out$mu <- cbind(matrix(res$pdSMu0, n.store, unit.par, byrow=TRUE),
+ matrix(res$pdSMu1, n.store, unit.par, byrow=TRUE))
+ colnames(res.out$mu) <- c("mu1", "mu2")
+ res.out$Sigma <- cbind(matrix(res$pdSSig00, n.store, unit.par, byrow=TRUE),
+ matrix(res$pdSSig01, n.store, unit.par, byrow=TRUE),
+ matrix(res$pdSSig11, n.store, unit.par, byrow=TRUE))
+ colnames(res.out$Sigma) <- c("Sigma11", "Sigma12", "Sigma22")
+ }
+
+ if (context)
+ class(res.out) <- c("ecoX","eco")
+ else
+ class(res.out) <- c("eco")
+
+ return(res.out)
+
+}
+
+
diff --git a/R/ecoBD.R b/R/ecoBD.R
new file mode 100644
index 0000000..e9dbdcf
--- /dev/null
+++ b/R/ecoBD.R
@@ -0,0 +1,117 @@
+ecoBD <- function(formula, data = parent.frame(), N=NULL){
+ mf <- match.call()
+ tt <- terms(formula)
+ attr(tt, "intercept") <- 0
+
+ if (is.matrix(eval.parent(mf$data)))
+ data <- as.data.frame(data)
+ X <- model.matrix(tt, data)
+ Y <- as.matrix(model.response(model.frame(tt, data = data)))
+ N <- eval(mf$N, data)
+ n.obs <- nrow(X)
+
+ ## counts
+ if (all(X>1) & all(Y>1)) {
+ if (!is.null(N)) {
+ if (!all(apply(X, 1, sum) == N))
+ X <- cbind(X, N-apply(X, 1, sum))
+ if (!all(apply(Y, 1, sum) == N))
+ Y <- cbind(Y, N-apply(Y, 1, sum))
+ if(any(X<0) || any(Y<0))
+ stop("Invalid inputs for X, Y, or/and N")
+ }
+ else {
+ if (!all(apply(X, 1, sum) == apply(Y, 1, sum)))
+ stop("X and Y do not sum to the same number. Input N.")
+ N <- apply(X, 1, sum)
+ }
+ C <- ncol(X)
+ R <- ncol(Y)
+ Wmin <- Wmax <- Nmin <- Nmax <- array(NA, c(n.obs, R, C))
+ clab <- rlab <- NULL
+ for (j in 1:C)
+ clab <- c(clab, paste("c", j, sep=""))
+ for (i in 1:R) {
+ rlab <- c(rlab, paste("r", i, sep=""))
+ for (j in 1:C) {
+ Nmin[,i,j] <- apply(cbind(0, X[,j]+Y[,i]-N), 1, max)
+ Nmax[,i,j] <- apply(cbind(Y[,i], X[,j]), 1, min)
+ Wmin[,i,j] <- Nmin[,i,j]/X[,j]
+ Wmax[,i,j] <- Nmax[,i,j]/X[,j]
+ }
+ }
+ dimnames(Wmin) <- dimnames(Wmax) <- dimnames(Nmin) <-
+ dimnames(Nmax) <-
+ list(if (is.null(rownames(X))) 1:n.obs else rownames(X),
+ rlab, clab)
+ }
+ else { ## proportions
+ if (any(apply(X, 1, sum) > 1.000000001))
+ stop("invalid input for X")
+ if (any(apply(X, 1, sum) < 0.9999999999))
+ X <- cbind(X, 1-X)
+ if (any(apply(Y, 1, sum) > 1.0000000001))
+ stop("invalid input for Y")
+ if (any(apply(Y, 1, sum) < 0.9999999999))
+ Y <- cbind(Y, 1-Y)
+ C <- ncol(X)
+ R <- ncol(Y)
+ Wmin <- Wmax <- array(NA, c(n.obs, R, C))
+ clab <- rlab <- NULL
+ for (j in 1:C)
+ clab <- c(clab, paste("c", j, sep=""))
+ for (i in 1:R) {
+ rlab <- c(rlab, paste("r", i, sep=""))
+ for (j in 1:C) {
+ Wmin[,i,j] <- apply(cbind(0, (X[,j]+Y[,i]-1)/X[,j]), 1, max)
+ Wmax[,i,j] <- apply(cbind(1, Y[,i]/X[,j]), 1, min)
+ }
+ }
+ dimnames(Wmin) <- dimnames(Wmax) <-
+ list(if (is.null(rownames(X))) 1:n.obs else rownames(X),
+ rlab, clab)
+ if (!is.null(N)) {
+ Nmin <- Nmax <- array(NA, c(n.obs, R, C), dimnames =
+ dimnames(Wmin))
+ for (i in 1:R)
+ for (j in 1:C) {
+ Nmin[,i,j] <- Wmin[,i,j]*X[,j]*N
+ Nmax[,i,j] <- Wmax[,i,j]*X[,j]*N
+ }
+ }
+ else
+ Nmin <- Nmax <- NULL
+ }
+
+ ## aggregate bounds
+ aggWmin <- aggWmax <- matrix(NA, R, C, dimnames =
+ list(dimnames(Wmin)[[2]], dimnames(Wmin)[[3]]))
+ if (is.null(N))
+ for (j in 1:C) {
+ aggWmin[,j] <- apply(Wmin[,,j]*X[,j], 2, mean)
+ aggWmax[,j] <- apply(Wmax[,,j]*X[,j], 2, mean)
+ }
+ else
+ for (j in 1:C) {
+ aggWmin[,j] <- apply(Wmin[,,j]*X[,j], 2, weighted.mean, N)
+ aggWmax[,j] <- apply(Wmax[,,j]*X[,j], 2, weighted.mean, N)
+ }
+
+ if (!is.null(Nmin) & !is.null(Nmax)) {
+ aggNmin <- aggNmax <- matrix(NA, R, C, dimnames =
+ list(dimnames(Nmin)[[2]], dimnames(Nmin)[[3]]))
+ for (j in 1:C) {
+ aggNmin[,j] <- apply(Nmin[,,j], 2, sum)
+ aggNmax[,j] <- apply(Nmax[,,j], 2, sum)
+ }
+ }
+ else
+ aggNmin <- aggNmax <- NULL
+
+ ## output
+ res <- list(call = mf, X = X, Y = Y, N = N, aggWmin = aggWmin,
+ aggWmax = aggWmax, aggNmin = aggNmin, aggNmax = aggNmax,
+ Wmin = Wmin, Wmax = Wmax, Nmin = Nmin, Nmax = Nmax)
+ class(res) <- c("ecoBD", "eco")
+ return(res)
+}
diff --git a/R/ecoNP.R b/R/ecoNP.R
new file mode 100644
index 0000000..a359e89
--- /dev/null
+++ b/R/ecoNP.R
@@ -0,0 +1,153 @@
+ecoNP <- function(formula, data = parent.frame(), N = NULL, supplement = NULL,
+ context = FALSE, mu0 = 0, tau0 = 2, nu0 = 4, S0 = 10,
+ alpha = NULL, a0 = 1, b0 = 0.1, parameter = FALSE,
+ grid = FALSE, n.draws = 5000, burnin = 0, thin = 0,
+ verbose = FALSE){
+
+ ## contextual effects
+ if (context)
+ ndim <- 3
+ else
+ ndim <- 2
+
+ ## checking inputs
+ if (burnin >= n.draws)
+ stop("n.draws should be larger than burnin")
+
+ if (length(mu0)==1)
+ mu0 <- rep(mu0, ndim)
+ else if (length(mu0)!=ndim)
+ stop("invalid inputs for mu0")
+ if (is.matrix(S0)) {
+ if (any(dim(S0)!=ndim))
+ stop("invalid inputs for S0")
+ }
+ else
+ S0 <- diag(S0, ndim)
+
+ mf <- match.call()
+
+ ## getting X, Y and N
+ tt <- terms(formula)
+ attr(tt, "intercept") <- 0
+ if (is.matrix(eval.parent(mf$data)))
+ data <- as.data.frame(data)
+ X <- model.matrix(tt, data)
+ Y <- model.response(model.frame(tt, data = data))
+ N <- eval(mf$N, data)
+
+ ## alpha
+ if (is.null(alpha)) {
+ alpha.update <- TRUE
+ alpha <- 0
+ }
+ else
+ alpha.update <- FALSE
+
+ ## checking the data and calculating the bounds
+ tmp <- checkdata(X, Y, supplement, ndim)
+ bdd <- ecoBD(formula, data=data)
+
+ ## fitting the model
+ n.store <- floor((n.draws-burnin)/(thin+1))
+ unit.par <- unit.w <- tmp$n.samp+tmp$samp.X1+tmp$samp.X0
+ n.par <- n.store * unit.par
+ n.w <- n.store * unit.w
+ unit.a <- 1
+
+ if (context)
+ res <- .C("cDPecoX", as.double(tmp$d), as.integer(tmp$n.samp),
+ as.integer(n.draws), as.integer(burnin), as.integer(thin+1),
+ as.integer(verbose), as.integer(nu0), as.double(tau0),
+ as.double(mu0), as.double(S0), as.double(alpha),
+ as.integer(alpha.update), as.double(a0), as.double(b0),
+ as.integer(tmp$survey.yes), as.integer(tmp$survey.samp),
+ as.double(tmp$survey.data), as.integer(tmp$X1type),
+ as.integer(tmp$samp.X1), as.double(tmp$X1.W1),
+ as.integer(tmp$X0type), as.integer(tmp$samp.X0),
+ as.double(tmp$X0.W2),
+ as.double(bdd$Wmin[,1,1]), as.double(bdd$Wmax[,1,1]),
+ as.integer(parameter), as.integer(grid),
+ pdSMu0=double(n.par), pdSMu1=double(n.par),
+ pdSMu2=double(n.par),
+ pdSSig00=double(n.par), pdSSig01=double(n.par),
+ pdSSig02=double(n.par), pdSSig11=double(n.par),
+ pdSSig12=double(n.par), pdSSig22=double(n.par),
+ pdSW1=double(n.w), pdSW2=double(n.w),
+ pdSa=double(n.store), pdSn=integer(n.store), PACKAGE="eco")
+ else
+ res <- .C("cDPeco", as.double(tmp$d), as.integer(tmp$n.samp),
+ as.integer(n.draws), as.integer(burnin), as.integer(thin+1),
+ as.integer(verbose), as.integer(nu0), as.double(tau0),
+ as.double(mu0), as.double(S0), as.double(alpha),
+ as.integer(alpha.update), as.double(a0), as.double(b0),
+ as.integer(tmp$survey.yes), as.integer(tmp$survey.samp),
+ as.double(tmp$survey.data), as.integer(tmp$X1type),
+ as.integer(tmp$samp.X1), as.double(tmp$X1.W1),
+ as.integer(tmp$X0type), as.integer(tmp$samp.X0),
+ as.double(tmp$X0.W2),
+ as.double(bdd$Wmin[,1,1]), as.double(bdd$Wmax[,1,1]),
+ as.integer(parameter), as.integer(grid),
+ pdSMu0=double(n.par), pdSMu1=double(n.par),
+ pdSSig00=double(n.par), pdSSig01=double(n.par),
+ pdSSig11=double(n.par), pdSW1=double(n.w), pdSW2=double(n.w),
+ pdSa=double(n.store), pdSn=integer(n.store), PACKAGE="eco")
+
+ ## output
+ W1.post <- matrix(res$pdSW1, n.store, unit.w, byrow=TRUE)[,tmp$order.old]
+ W2.post <- matrix(res$pdSW2, n.store, unit.w, byrow=TRUE)[,tmp$order.old]
+ W <- array(rbind(W1.post, W2.post), c(n.store, 2, unit.w))
+ colnames(W) <- c("W1", "W2")
+ res.out <- list(call = mf, X = X, Y = Y, N = N, W = W,
+ Wmin = bdd$Wmin[,1,], Wmax = bdd$Wmax[,1,],
+ burin = burnin, thin = thin, nu0 = nu0, tau0 = tau0,
+ mu0 = mu0, a0 = a0, b0 = b0, S0 = S0)
+
+ ## optional outputs
+ if (parameter){
+ if (context) {
+ mu1.post <- matrix(res$pdSMu0, n.store, unit.par, byrow=TRUE)[,tmp$order.old]
+ mu2.post <- matrix(res$pdSMu1, n.store, unit.par, byrow=TRUE)[,tmp$order.old]
+ mu3.post <- matrix(res$pdSMu2, n.store, unit.par, byrow=TRUE)[,tmp$order.old]
+ Sigma11.post <- matrix(res$pdSSig00, n.store, unit.par, byrow=TRUE)[,tmp$order.old]
+ Sigma12.post <- matrix(res$pdSSig01, n.store, unit.par, byrow=TRUE)[,tmp$order.old]
+ Sigma13.post <- matrix(res$pdSSig02, n.store, unit.par, byrow=TRUE)[,tmp$order.old]
+ Sigma23.post <- matrix(res$pdSSig12, n.store, unit.par, byrow=TRUE)[,tmp$order.old]
+ Sigma22.post <- matrix(res$pdSSig11, n.store, unit.par, byrow=TRUE)[,tmp$order.old]
+ Sigma33.post <- matrix(res$pdSSig22, n.store, unit.par, byrow=TRUE)[,tmp$order.old]
+ res.out$mu <- array(rbind(mu1.post, mu2.post, mu3.post),
+ dim=c(n.store, 3, unit.par),
+ dimnames=list(1:n.store, c("mu1", "mu2", "mu3"), 1:unit.par))
+ res.out$Sigma <- array(rbind(Sigma11.post, Sigma12.post, Sigma13.post,
+ Sigma22.post, Sigma23.post, Sigma33.post),
+ dim=c(n.store, 6, unit.par),
+ dimnames=list(1:n.store, c("Sigma11",
+ "Sigma12", "Sigma13", "Sigma22", "Sigma23", "Sigma33"), 1:unit.par))
+ }
+ else {
+ mu1.post <- matrix(res$pdSMu0, n.store, unit.par, byrow=TRUE)[,tmp$order.old]
+ mu2.post <- matrix(res$pdSMu1, n.store, unit.par, byrow=TRUE)[,tmp$order.old]
+ Sigma11.post <- matrix(res$pdSSig00, n.store, unit.par, byrow=TRUE)[,tmp$order.old]
+ Sigma12.post <- matrix(res$pdSSig01, n.store, unit.par, byrow=TRUE)[,tmp$order.old]
+ Sigma22.post <- matrix(res$pdSSig11, n.store, unit.par, byrow=TRUE)[,tmp$order.old]
+ res.out$mu <- array(rbind(mu1.post, mu2.post), dim=c(n.store, 2, unit.par),
+ dimnames=list(1:n.store, c("mu1", "mu2"), 1:unit.par))
+ res.out$Sigma <- array(rbind(Sigma11.post, Sigma12.post, Sigma22.post),
+ dim=c(n.store, 3, unit.par),
+ dimnames=list(1:n.store, c("Sigma11", "Sigma12", "Sigma22"), 1:unit.par))
+ }
+ if (alpha.update)
+ res.out$alpha <- matrix(res$pdSa, n.store, unit.a, byrow=TRUE)
+ else
+ res.out$alpha <- alpha
+ res.out$nstar <- matrix(res$pdSn, n.store, unit.a, byrow=TRUE)
+ }
+
+ if (context)
+ class(res.out) <- c("ecoNPX", "ecoNP", "eco")
+ else
+ class(res.out) <- c("ecoNP", "eco")
+ return(res.out)
+}
+
+
diff --git a/R/logit.R b/R/logit.R
new file mode 100644
index 0000000..6cf6e90
--- /dev/null
+++ b/R/logit.R
@@ -0,0 +1,5 @@
+logit <- function(x)
+ return(log(x)-log(1-x))
+
+invlogit <- function(x)
+ return(exp(x)/(1+exp(x)))
diff --git a/R/onLoad.R b/R/onLoad.R
new file mode 100644
index 0000000..97b17c1
--- /dev/null
+++ b/R/onLoad.R
@@ -0,0 +1,7 @@
+".onLoad" <- function(lib, pkg) {
+ mylib <- dirname(system.file(package = pkg))
+ title <- packageDescription(pkg, lib = mylib)$Title
+ ver <- packageDescription(pkg, lib = mylib)$Version
+ cat(paste(pkg, ": ", title, "\nVersion ", ver, "\n", sep=""))
+}
+
diff --git a/R/predict.eco.R b/R/predict.eco.R
new file mode 100644
index 0000000..f0f9819
--- /dev/null
+++ b/R/predict.eco.R
@@ -0,0 +1,34 @@
+predict.eco <- function(object, newdraw = NULL, subset = NULL,
+ verbose = FALSE, ...){
+
+ if (is.null(newdraw) && is.null(object$mu))
+ stop("Posterior draws of mu and Sigma must be supplied")
+ else if (!is.null(newdraw)){
+ if (is.null(newdraw$mu) && is.null(newdraw$Sigma))
+ stop("Posterior draws of both mu and Sigma must be supplied.")
+ object <- newdraw
+ }
+
+ mu <- coef(object, subset = subset)
+ n.draws <- nrow(mu)
+ p <- ncol(mu)
+ Sigma <- varcov(object, subset = subset)
+
+ Wstar <- matrix(NA, nrow=n.draws, ncol=p)
+ tmp <- floor(n.draws/10)
+ inc <- 1
+ for (i in 1:n.draws) {
+ Wstar[i,] <- mvrnorm(1, mu = mu[i,], Sigma = Sigma[,,i])
+ if (i == inc*tmp & verbose) {
+ cat("", inc*10, "percent done.\n")
+ inc <- inc + 1
+ }
+ }
+ res <- apply(Wstar, 2, invlogit)
+ if (ncol(res) == 2)
+ colnames(res) <- c("W1", "W2")
+ else # this is called from predict.ecoX
+ colnames(res) <- c("W1", "W2", "X")
+ class(res) <- c("predict.eco", "matrix")
+ return(res)
+}
diff --git a/R/predict.ecoNP.R b/R/predict.ecoNP.R
new file mode 100644
index 0000000..32a77b6
--- /dev/null
+++ b/R/predict.ecoNP.R
@@ -0,0 +1,32 @@
+predict.ecoNP <- function(object, newdraw = NULL, subset = NULL,
+ obs = NULL, verbose = FALSE, ...){
+
+ if (is.null(newdraw) && is.null(object$mu))
+ stop("Posterior draws of mu and Sigma must be supplied")
+ else if (!is.null(newdraw)){
+ if (is.null(newdraw$mu) && is.null(newdraw$Sigma))
+ stop("Posterior draws of both mu and Sigma must be supplied.")
+ object <- newdraw
+ }
+
+ n.draws <- dim(object$mu)[1]
+ p <- dim(object$mu)[2]
+ n <- dim(object$mu)[3]
+ mu <- aperm(coef(object, subset = subset, obs = obs), c(2,3,1))
+
+ if (is.null(subset))
+ subset <- 1:n.draws
+ if (is.null(obs))
+ obs <- 1:n
+ Sigma <- aperm(object$Sigma[subset,,obs], c(2,3,1))
+
+ res <- .C("preDP", as.double(mu), as.double(Sigma), as.integer(n),
+ as.integer(n.draws), as.integer(p), as.integer(verbose),
+ pdStore = double(n.draws*p*n), PACKAGE="eco")$pdStore
+
+ res <- matrix(res, ncol=p, nrow=n.draws*n, byrow=TRUE)
+ colnames(res) <- c("W1", "W2")
+
+ class(res) <- c("predict.eco", "matrix")
+ return(res)
+}
diff --git a/R/predict.ecoNPX.R b/R/predict.ecoNPX.R
new file mode 100644
index 0000000..adf64ad
--- /dev/null
+++ b/R/predict.ecoNPX.R
@@ -0,0 +1,42 @@
+predict.ecoNPX <- function(object, newdraw = NULL, subset = NULL,
+ obs = NULL, cond = FALSE, verbose = FALSE, ...){
+
+ if (is.null(newdraw) && is.null(object$mu))
+ stop("Posterior draws of mu and Sigma must be supplied")
+ else if (!is.null(newdraw)){
+ if (is.null(newdraw$mu) && is.null(newdraw$Sigma))
+ stop("Posterior draws of both mu and Sigma must be supplied.")
+ object <- newdraw
+ }
+
+ n.draws <- dim(object$mu)[1]
+ n <- dim(object$mu)[3]
+ mu <- aperm(coef(object, subset = subset, obs = obs), c(2,3,1))
+
+ if (is.null(subset))
+ subset <- 1:n.draws
+ if (is.null(obs))
+ obs <- 1:n
+ Sigma <- aperm(object$Sigma[subset,,obs], c(2,3,1))
+
+ if (cond) { # conditional prediction
+ X <- object$X
+ res <- .C("preDPX", as.double(mu), as.double(Sigma), as.double(X),
+ as.integer(n), as.integer(n.draws), as.integer(2),
+ as.integer(verbose), pdStore = double(n.draws*2*n),
+ PACKAGE="eco")$pdStore
+ res <- matrix(res, ncol=2, nrow=n.draws*n, byrow=TRUE)
+ colnames(res) <- c("W1", "W2")
+ }
+ else { # unconditional prediction
+ res <- .C("preDP", as.double(mu), as.double(Sigma), as.integer(n),
+ as.integer(n.draws), as.integer(3), as.integer(verbose),
+ pdStore = double(n.draws*3*n), PACKAGE="eco")$pdStore
+
+ res <- matrix(res, ncol=3, nrow=n.draws*n, byrow=TRUE)
+ colnames(res) <- c("W1", "W2", "X")
+ }
+
+ class(res) <- c("predict.eco", "matrix")
+ return(res)
+}
diff --git a/R/predict.ecoX.R b/R/predict.ecoX.R
new file mode 100644
index 0000000..fbe2292
--- /dev/null
+++ b/R/predict.ecoX.R
@@ -0,0 +1,43 @@
+predict.ecoX <- function(object, newdraw = NULL, subset = NULL,
+ newdata = NULL, cond = FALSE, verbose = FALSE, ...){
+
+ if (is.null(newdraw) && is.null(object$mu))
+ stop("Posterior draws of mu and Sigma must be supplied")
+ else if (!is.null(newdraw)){
+ if (is.null(newdraw$mu) && is.null(newdraw$Sigma))
+ stop("Posterior draws of both mu and Sigma must be supplied.")
+ object <- newdraw
+ }
+
+ if (cond) { ## conditional prediction
+ mu <- coef(object, subset = subset)
+ n.draws <- nrow(mu)
+ if (is.null(subset))
+ subset <- 1:n.draws
+ Sigma <- object$Sigma[subset,]
+ if (is.null(newdata))
+ X <- object$X
+ else {
+ mf <- match.call()
+ if (is.matrix(eval.parent(mf$newdata)))
+ data <- as.data.frame(data)
+ tt <- terms(object)
+ attr(tt, "intercept") <- 0
+ X <- model.matrix(tt, newdata)
+ }
+ n <- nrow(X)
+ res <- .C("preBaseX", as.double(X), as.double(mu), as.double(t(Sigma)),
+ as.integer(length(c(X))), as.integer(nrow(mu)),
+ as.integer(verbose),
+ pdStore = double(n.draws*n*2), PACKAGE="eco")$pdStore
+ res <- array(res, c(2, n, n.draws), dimnames=list(c("W1", "W2"),
+ rownames(X), 1:n.draws))
+ class(res) <- c("predict.ecoX", "array")
+ }
+ else {
+ res <- predict.eco(object, newdraw = newdraw, subset = subset,
+ newdata = newdata, verbose = verbose, ...)
+ colnames(res) <- c("W1", "W2", "X")
+ }
+ return(res)
+}
diff --git a/R/print.eco.R b/R/print.eco.R
new file mode 100644
index 0000000..c2554a7
--- /dev/null
+++ b/R/print.eco.R
@@ -0,0 +1,19 @@
+print.eco <- function(x, digits = max(3, getOption("digits") -3),
+ ...){
+ cat("\nCall:\n", deparse(x$call), "\n\n", sep="")
+
+ if (is.null(x$N))
+ N <- rep(1, nrow(x$X))
+ else N <- x$N
+
+ W.mean <- cbind(mean(x$W[,1,] %*% (x$X*N/sum(x$X*N))),
+ mean(x$W[,2,] %*% ((1-x$X)*N/sum((1-x$X)*N))))
+ colnames(W.mean) <- c("W1", "W2")
+ rownames(W.mean) <- "posterior mean"
+
+ cat("Aggregate In-sample Estimates:\n\n")
+ print.default(format(W.mean, digits = digits), print.gap = 2, quote =
+ FALSE)
+ cat("\n")
+ invisible(x)
+}
diff --git a/R/print.ecoBD.R b/R/print.ecoBD.R
new file mode 100644
index 0000000..26370a7
--- /dev/null
+++ b/R/print.ecoBD.R
@@ -0,0 +1,23 @@
+print.ecoBD <- function(x, digits = max(3, getOption("digits") -3),
+ ...) {
+ cat("\nCall:\n", deparse(x$call), "\n\n", sep="")
+
+ cat("Aggregate Lower Bounds (Proportions):\n")
+ print.default(format(x$aggWmin, digits = digits), print.gap = 2, quote =
+ FALSE)
+ cat("\nAggregate Upper Bounds (Proportions):\n")
+ print.default(format(x$aggWmax, digits = digits), print.gap = 2, quote =
+ FALSE)
+
+ if (!is.null(x$aggNmin)) {
+ cat("\nAggregate Lower Bounds (Counts):\n")
+ print.default(format(x$aggNmin, digits = digits), print.gap = 2, quote =
+ FALSE)
+ cat("\nAggregate Upper Bounds (Counts):\n")
+ print.default(format(x$aggNmax, digits = digits), print.gap = 2, quote =
+ FALSE)
+ }
+
+ cat("\n")
+ invisible(x)
+}
diff --git a/R/print.summary.eco.R b/R/print.summary.eco.R
new file mode 100644
index 0000000..8c8df40
--- /dev/null
+++ b/R/print.summary.eco.R
@@ -0,0 +1,27 @@
+print.summary.eco <- function(x, digits=max(3, getOption("digits")-3), ...) {
+ cat("\nCall: ")
+ cat(paste(deparse(x$call), sep="\n", collapse="\n"))
+
+ cat("\n")
+ if (!is.null(x$param.table)) {
+ cat("\nParameter Estimates:\n")
+ printCoefmat(x$param.table, digits=digits, na.print="NA",...)
+ }
+
+ cat("\nAggregate Estimates:\n")
+ printCoefmat(x$agg.table, digits=digits, na.print="NA",...)
+
+ cat("\nNumber of Units:", x$n.obs)
+ cat("\nNumber of Monte Carlo Draws:", x$n.draws)
+
+
+ if (!is.null(x$W1.table)) {
+ cat("\n\nUnit-level Estimates of W1:\n")
+ printCoefmat(x$W1.table, digits=digits, na.print="NA",...)
+ cat("\n\nUnit-level Estimates of W2:\n")
+ printCoefmat(x$W2.table, digits=digits, na.print="NA",...)
+ }
+
+ cat("\n")
+ invisible(x)
+}
diff --git a/R/print.summary.ecoNP.R b/R/print.summary.ecoNP.R
new file mode 100644
index 0000000..c9c3b2b
--- /dev/null
+++ b/R/print.summary.ecoNP.R
@@ -0,0 +1,33 @@
+print.summary.ecoNP <- function(x, digits=max(3, getOption("digits")-3), ...)
+ {
+ cat("\nCall: ")
+ cat(paste(deparse(x$call), sep="\n", collapse="\n"))
+
+ cat("\n\nAggregate Estimates:\n")
+ printCoefmat(x$agg.table, digits=digits, na.print="NA",...)
+ cat("\nNumber of Units:", x$n.obs)
+ cat("\nNumber of Monte Carlo Draws:", x$n.draws)
+ if (!is.null(x$param.table)) {
+ tt <- x$param.table
+ cat("\nParameter Estimates of mu1:\n")
+ printCoefmat(tt$mu1.table, digits=digits, na.print="NA",...)
+ cat("\nParameter Estimates of mu2:\n")
+ printCoefmat(tt$mu2.table, digits=digits, na.print="NA",...)
+ cat("\nParameter Estimates of Sigma11:\n")
+ printCoefmat(tt$Sigma11.table, digits=digits, na.print="NA",...)
+ cat("\nParameter Estimates of Sigma12:\n")
+ printCoefmat(tt$Sigma12.table, digits=digits, na.print="NA",...)
+ cat("\nParameter Estimates of Sigma22:\n")
+ printCoefmat(tt$Sigma22.table, digits=digits, na.print="NA",...)
+ }
+
+ if (!is.null(x$W1.table)) {
+ cat("\n\nUnit-level Estimates of W1:\n")
+ printCoefmat(x$W1.table, digits=digits, na.print="NA",...)
+ cat("\n\nUnit-level Estimates of W2:\n")
+ printCoefmat(x$W2.table, digits=digits, na.print="NA",...)
+ }
+
+ cat("\n")
+ invisible(x)
+}
diff --git a/R/print.summary.predict.eco.R b/R/print.summary.predict.eco.R
new file mode 100644
index 0000000..5a546c5
--- /dev/null
+++ b/R/print.summary.predict.eco.R
@@ -0,0 +1,10 @@
+print.summary.predict.eco <- function(x, digits=max(3, getOption("digits")
+ -3), ...) {
+ cat("\nOut-of-sample Prediction:\n")
+ printCoefmat(x$W.table, digits=digits, na.print="NA",...)
+
+ cat("\nNumber of Monte Carlo Draws:", x$n.draws)
+
+ cat("\n")
+ invisible(x)
+}
diff --git a/R/summary.eco.R b/R/summary.eco.R
new file mode 100644
index 0000000..29785ff
--- /dev/null
+++ b/R/summary.eco.R
@@ -0,0 +1,66 @@
+summary.eco <- function(object, CI = c(2.5, 97.5), param = TRUE,
+ units = FALSE, subset = NULL,...) {
+
+ n.obs <- ncol(object$W[,1,])
+ n.draws <- nrow(object$W[,1,])
+
+ if (is.null(subset)) subset <- 1:n.obs
+ else if (!is.numeric(subset))
+ stop("Subset should be a numeric vector.")
+ else if (!all(subset %in% c(1:n.obs)))
+ stop("Subset should be any numbers in 1:obs.")
+
+ table.names<-c("mean", "std.dev", paste(min(CI), "%", sep=" "),
+ paste(max(CI), "%", sep=" "))
+
+ if (is.null(object$N))
+ N <- rep(1, nrow(object$X))
+ else N <- object$N
+
+ W1.agg.mean <- object$W[,1,] %*% (object$X*N/sum(object$X*N))
+ W2.agg.mean <- object$W[,2,] %*% ((1-object$X)*N/sum((1-object$X)*N))
+ agg.table <- rbind(cbind(mean(W1.agg.mean), sd(W1.agg.mean),
+ quantile(W1.agg.mean, min(CI)/100),
+ quantile(W1.agg.mean, max(CI)/100)),
+ cbind(mean(W2.agg.mean), sd(W2.agg.mean),
+ quantile(W2.agg.mean, min(CI)/100),
+ quantile(W2.agg.mean, max(CI)/100)))
+ colnames(agg.table) <- table.names
+ rownames(agg.table) <- c("W1", "W2")
+
+ if (units) {
+ W1.table <- cbind(apply(object$W[,1,subset], 2, mean),
+ apply(object$W[,1,subset], 2, sd),
+ apply(object$W[,1,subset], 2, quantile, min(CI)/100),
+ apply(object$W[,1,subset], 2, quantile, max(CI)/100))
+ W2.table <- cbind(apply(object$W[,2,subset], 2, mean),
+ apply(object$W[,2,subset], 2, sd),
+ apply(object$W[,2,subset], 2, quantile, min(CI)/100),
+ apply(object$W[,2,subset], 2, quantile, max(CI)/100))
+ colnames(W2.table) <- colnames(W1.table) <- table.names
+ rownames(W1.table) <- rownames(W2.table) <- row.names(object$X[subset])
+ }
+ else
+ W1.table <- W2.table <- NULL
+
+ if (param) {
+ if (is.null(object$mu) || is.null(object$Sigma))
+ stop("Parameters are missing values.")
+ else {
+ param <- cbind(object$mu, object$Sigma)
+ param.table <- cbind(apply(param, 2, mean), apply(param, 2, sd),
+ apply(param, 2, quantile, min(CI)/100),
+ apply(param, 2, quantile, max(CI)/100))
+ colnames(param.table) <- table.names
+ }
+ }
+ else
+ param.table <- NULL
+
+ ans <- list(call = object$call, W1.table = W1.table, W2.table = W2.table,
+ agg.table = agg.table, param.table = param.table,
+ n.draws = n.draws, n.obs = n.obs)
+
+ class(ans) <-"summary.eco"
+ return(ans)
+}
diff --git a/R/summary.ecoNP.R b/R/summary.ecoNP.R
new file mode 100644
index 0000000..2134a88
--- /dev/null
+++ b/R/summary.ecoNP.R
@@ -0,0 +1,86 @@
+summary.ecoNP <- function(object, CI=c(2.5, 97.5), param=FALSE, units=FALSE, subset=NULL,...) {
+
+
+ n.obs <- ncol(object$W[,1,])
+ n.draws <- nrow(object$W[,1,])
+
+ if (is.null(subset)) subset <- 1:n.obs
+ else if (!is.numeric(subset)) stop("Subset should be a numeric vector.")
+ else if (!all(subset %in% c(1:n.obs))) stop("Subset should be any numbers in 1:obs.")
+
+ table.names<-c("mean", "std.dev", paste(min(CI), "%", sep=" "), paste(max(CI), "%", sep=" "))
+
+
+ if (is.null(object$N))
+ N <- rep(1, nrow(object$X))
+ else N <- object$N
+
+ W1.agg.mean <- object$W[,1,] %*% (object$X*N/sum(object$X*N))
+ W2.agg.mean <- object$W[,2,] %*% ((1-object$X)*N/sum((1-object$X)*N))
+ agg.table <- rbind(cbind(mean(W1.agg.mean), sd(W1.agg.mean),
+ quantile(W1.agg.mean, min(CI)/100),
+ quantile(W1.agg.mean, max(CI)/100)),
+ cbind(mean(W2.agg.mean), sd(W2.agg.mean),
+ quantile(W2.agg.mean, min(CI)/100),
+ quantile(W2.agg.mean, max(CI)/100)))
+ colnames(agg.table) <- table.names
+ rownames(agg.table) <- c("W1", "W2")
+
+ if (units) {
+ W1.table <- cbind(apply(object$W[,1,subset], 2, mean),
+ apply(object$W[,1,subset], 2, sd),
+ apply(object$W[,1,subset], 2, quantile, min(CI)/100),
+ apply(object$W[,1,subset], 2, quantile, max(CI)/100))
+ W2.table <- cbind(apply(object$W[,2,subset], 2, mean),
+ apply(object$W[,2,subset], 2, sd),
+ apply(object$W[,2,subset], 2, quantile, min(CI)/100),
+ apply(object$W[,2,subset], 2, quantile, max(CI)/100))
+ colnames(W2.table) <- colnames(W1.table) <- table.names
+ rownames(W1.table) <- rownames(W2.table) <- row.names(object$X[subset])
+ }
+ else
+ W1.table <- W2.table <- NULL
+
+ if (is.null(param)) param <- FALSE
+ if (param) {
+ if (is.null(object$mu) || is.null(object$Sigma))
+ stop("Parameters are missing values.")
+ }
+
+
+ if (param) {
+ mu1.table <- cbind(apply(object$mu[,1,subset], 2, mean),
+ apply(object$mu[,1,subset], 2, sd),
+ apply(object$mu[,1,subset], 2, quantile, min(CI)/100),
+ apply(object$mu[,1,subset], 2, quantile, max(CI)/100))
+ mu2.table <- cbind(apply(object$mu[,2,subset], 2, mean),
+ apply(object$mu[,2,subset], 2, sd),
+ apply(object$mu[,2,subset], 2, quantile, min(CI)/100),
+ apply(object$mu[,2,subset], 2, quantile, max(CI)/100))
+ Sigma11.table <- cbind(apply(object$Sigma[,1,subset], 2, mean),
+ apply(object$Sigma[,1,subset], 2, sd),
+ apply(object$Sigma[,1,subset], 2, quantile, min(CI)/100),
+ apply(object$Sigma[,1,subset], 2, quantile, max(CI)/100))
+ Sigma12.table <- cbind(apply(object$Sigma[,2,subset], 2, mean),
+ apply(object$Sigma[,2,subset], 2, sd),
+ apply(object$Sigma[,2,subset], 2, quantile, min(CI)/100),
+ apply(object$Sigma[,2,subset], 2, quantile, max(CI)/100))
+ Sigma22.table <- cbind(apply(object$Sigma[,3,subset], 2, mean),
+ apply(object$Sigma[,3,subset], 2, sd),
+ apply(object$Sigma[,3,subset], 2, quantile, min(CI)/100),
+ apply(object$Sigma[,3,subset], 2, quantile, max(CI)/100))
+
+ colnames(mu1.table) <- colnames(mu2.table) <- table.names
+ colnames(Sigma11.table) <- colnames(Sigma12.table) <- colnames(Sigma22.table) <- table.names
+ param.table=list(mu1.table=mu1.table,mu2.table=mu2.table,Sigma11.table=Sigma11.table,Sigma12.table=Sigma12.table,Sigma22.table=Sigma22.table)
+ }
+ else
+ param.table <- NULL
+
+ ans <- list(call = object$call, W1.table = W1.table, W2.table = W2.table,
+ agg.table = agg.table, param.table = param.table,
+ n.draws = n.draws, n.obs = n.obs)
+
+ class(ans) <-c("summary.eco", "summary.ecoNP")
+ return(ans)
+}
diff --git a/R/summary.predict.eco.R b/R/summary.predict.eco.R
new file mode 100644
index 0000000..b9d19e3
--- /dev/null
+++ b/R/summary.predict.eco.R
@@ -0,0 +1,20 @@
+summary.predict.eco <- function(object, CI=c(2.5, 97.5), ...) {
+
+ if (any(CI < 0) || any(CI > 100))
+ stop("Invalid input for CI")
+ n.draws <- nrow(object)
+ n.var <- ncol(object)
+ table.names<-c("mean", "std.dev", paste(min(CI), "%", sep=" "),
+ paste(max(CI), "%", sep=" "))
+ W.table <- matrix(NA, ncol=length(table.names), nrow=n.var)
+ for (i in 1:n.var)
+ W.table[i,] <- cbind(mean(object[,i]), sd(object[,i]),
+ quantile(object[,i], min(CI)/100),
+ quantile(object[,i], max(CI)/100))
+ colnames(W.table) <- table.names
+ rownames(W.table) <- colnames(object)
+
+ res <- list(W.table = W.table, n.draws = n.draws)
+ class(res) <- "summary.predict.eco"
+ return(res)
+}
diff --git a/R/varcov.R b/R/varcov.R
new file mode 100644
index 0000000..46d36d9
--- /dev/null
+++ b/R/varcov.R
@@ -0,0 +1,71 @@
+varcov <- function(object, ...)
+ UseMethod("varcov")
+
+varcov.eco <- function(object, subset = NULL, ...) {
+ if (is.null(subset))
+ subset <- 1:nrow(object$Sigma)
+ else if (max(subset) > nrow(object$Sigma))
+ stop(paste("invalid input for `subset.' only", nrow(object$Sigma),
+ "draws are stored."))
+
+ p <- ncol(object$mu)
+ n <- length(subset)
+ Sigma <- array(0, c(p, p, n))
+ cov <- object$Sigma[subset,]
+
+ for (i in 1:n) {
+ count <- 1
+ for (j in 1:p) {
+ Sigma[j,j:p,i] <- cov[i,count:(count+p-j)]
+ count <- count + p - j + 1
+ }
+ diag(Sigma[,,i]) <- diag(Sigma[,,i]/2)
+ Sigma[,,i] <- Sigma[,,i] + t(Sigma[,,i])
+ }
+ if (n > 1)
+ return(Sigma)
+ else
+ return(Sigma[,,1])
+}
+
+varcov.ecoNP <- function(object, subset = NULL, obs = NULL, ...) {
+ if (is.null(subset))
+ subset <- 1:nrow(object$Sigma)
+ else if (max(subset) > nrow(object$Sigma))
+ stop(paste("invalid input for `subset.' only", nrow(object$Sigma),
+ "draws are stored."))
+
+ if (is.null(obs))
+ obs <- 1:dim(object$Sigma)[3]
+ else if (max(subset) > dim(object$Sigma)[3])
+ stop(paste("invalid input for `obs.' only", dim(object$Sigma)[3],
+ "draws are stored."))
+
+ p <- ncol(object$mu)
+ n <- length(subset)
+ m <- length(obs)
+ Sigma <- array(0, c(p, p, n, m))
+ cov <- object$Sigma[subset,,obs]
+
+ for (k in 1:m) {
+ for (i in 1:n) {
+ count <- 1
+ for (j in 1:p) {
+ Sigma[j,j:p,i,k] <- cov[i,count:(count+p-j),k]
+ count <- count + p - j + 1
+ }
+ diag(Sigma[,,i,k]) <- diag(Sigma[,,i,k]/2)
+ Sigma[,,i,k] <- Sigma[,,i,k] + t(Sigma[,,i,k])
+ }
+ }
+ if (n > 1)
+ if (m > 1)
+ return(Sigma)
+ else
+ return(Sigma[,,,1])
+ else
+ if (m > 1)
+ return(Sigma[,,1,])
+ else
+ return(Sigma[,,1,1])
+}
diff --git a/data/census.txt b/data/census.txt
new file mode 100644
index 0000000..1207217
--- /dev/null
+++ b/data/census.txt
@@ -0,0 +1,1041 @@
+Y X N W1 W2
+0.7372 0.5905 14497 0.6061 0.9262
+0.8518 0.2803 13263 0.646 0.932
+0.6833 0.6147 23193 0.5374 0.916
+0.8323 0.3533 16030 0.7078 0.9003
+0.8529 0.0585 14649 0.6161 0.8676
+0.6412 0.8279 21242 0.574 0.9647
+0.7407 0.5213 20618 0.5785 0.9172
+0.8366 0.2803 28537 0.6885 0.8942
+0.767 0.5082 25580 0.608 0.9313
+0.8011 0.1293 13924 0.6211 0.8279
+0.8212 0.2048 16309 0.653 0.8645
+0.6967 0.6109 12695 0.5598 0.9115
+0.7303 0.5513 22088 0.5787 0.9164
+0.8628 0.1275 14452 0.6308 0.8967
+0.8265 0.0542 9164 0.5936 0.8399
+0.7738 0.2223 18122 0.5492 0.8379
+0.7925 0.3867 17782 0.6247 0.8984
+0.7538 0.4603 15285 0.5975 0.8871
+0.7963 0.363 11529 0.5885 0.9148
+0.7913 0.2624 22830 0.6087 0.8563
+0.8267 0.3106 16489 0.6658 0.8993
+0.7934 0.263 15138 0.5901 0.866
+0.6144 0.8081 40389 0.5286 0.976
+0.7807 0.4734 20305 0.6074 0.9365
+0.8126 0.3024 13274 0.6522 0.8822
+0.8371 0.184 27995 0.6943 0.8693
+0.8345 0.1145 11133 0.6314 0.8607
+0.8104 0.1048 13220 0.6332 0.8312
+0.7971 0.1642 17927 0.5751 0.8407
+0.5676 0.8566 16268 0.5002 0.9704
+0.6269 0.7792 19983 0.5431 0.9227
+0.7477 0.4678 14306 0.5862 0.8895
+0.8149 0.291 22857 0.6554 0.8804
+0.8116 0.0952 22956 0.6482 0.8288
+0.8664 0.4205 174724 0.7458 0.9539
+0.8584 0.1776 11968 0.6226 0.9093
+0.8272 0.2357 21832 0.6413 0.8846
+0.7578 0.3262 15324 0.6279 0.8207
+0.7032 0.5794 23629 0.5257 0.9478
+0.7123 0.3861 18981 0.4527 0.8756
+0.5491 0.8739 23040 0.4894 0.9628
+0.6893 0.8353 18755 0.6368 0.9556
+0.7698 0.4024 34607 0.6003 0.8839
+0.562 0.7699 28845 0.4477 0.9442
+0.8555 0.4382 63959 0.7108 0.9684
+0.6888 0.5749 18840 0.5159 0.9227
+0.7518 0.6843 63652 0.6457 0.9817
+0.8468 0.2515 24621 0.6652 0.9078
+0.6464 0.7817 22403 0.5636 0.9427
+0.6518 0.5172 17805 0.4099 0.911
+0.7444 0.4599 22004 0.5545 0.906
+0.8281 0.2322 16792 0.6986 0.8672
+0.6315 0.7607 18481 0.5424 0.9145
+0.8552 0.1821 14523 0.7259 0.884
+0.7966 0.2922 19093 0.6078 0.8745
+0.5685 0.8022 20787 0.4666 0.982
+0.7626 0.4731 26777 0.624 0.887
+0.7641 0.3608 21781 0.5399 0.8906
+0.7606 0.4047 34799 0.5841 0.8805
+0.8377 0.1951 25780 0.6962 0.872
+0.7689 0.4204 10163 0.589 0.8993
+0.6345 0.8064 24815 0.5552 0.9644
+0.911 0.2591 11866 0.7625 0.9629
+0.8668 0.5366 18110 0.7839 0.9628
+0.8619 0.3161 10045 0.6844 0.944
+0.8228 0.3487 6785 0.6095 0.9371
+0.6932 0.8067 16691 0.6601 0.8314
+0.8846 0.3042 16784 0.7391 0.9482
+0.8431 0.3095 9215 0.6473 0.9309
+0.8262 0.4411 16580 0.6659 0.9527
+0.844 0.3692 15765 0.7856 0.8782
+0.9075 0.0531 19911 0.7483 0.9164
+0.8959 0.0927 17007 0.7502 0.9108
+0.8133 0.8471 17194 0.7887 0.9498
+0.8511 0.4415 10208 0.7464 0.9339
+0.8942 0.3574 8830 0.743 0.9783
+0.6898 0.7962 11585 0.6212 0.9576
+0.8068 0.5353 15609 0.6723 0.9617
+0.9039 0.1903 16328 0.7909 0.9305
+0.9259 0.187 21353 0.8428 0.945
+0.9151 0.1073 6340 0.7809 0.9313
+0.8581 0.4869 20002 0.7361 0.9738
+0.9141 0.1373 10633 0.8103 0.9307
+0.8984 0.2031 11792 0.7445 0.9376
+0.9189 0.056 17568 0.7713 0.9277
+0.8625 0.2736 16954 0.776 0.8951
+0.8076 0.7142 39918 0.7413 0.9734
+0.8552 0.5264 9985 0.7574 0.9638
+0.8562 0.7843 18208 0.8291 0.9547
+0.8138 0.674 11030 0.7382 0.97
+0.8238 0.4197 9735 0.662 0.9409
+0.8448 0.4056 19982 0.6896 0.9507
+0.8348 0.3763 14424 0.6588 0.941
+0.7861 0.4602 22686 0.6669 0.8877
+0.8006 0.6327 14707 0.7234 0.9335
+0.843 0.341 13251 0.631 0.9527
+0.8652 0.5541 15419 0.7856 0.9642
+0.8839 0.1016 6406 0.745 0.8996
+0.8105 0.7838 25841 0.7659 0.9723
+0.9037 0.0767 8644 0.7541 0.9162
+0.938 0.1765 9305 0.8234 0.9625
+0.8906 0.0787 16909 0.7911 0.8991
+0.8888 0.3233 9917 0.7723 0.9444
+0.8953 0.412 68630 0.7876 0.9707
+0.7962 0.6837 16656 0.7335 0.9319
+0.8824 0.1197 11699 0.7407 0.9016
+0.9426 0.1122 38881 0.8221 0.9579
+0.914 0.1355 11422 0.7319 0.9426
+0.7948 0.4321 21282 0.5936 0.9478
+0.9266 0.0769 20369 0.7518 0.9411
+0.8267 0.5829 14471 0.7419 0.9452
+0.8958 0.0704 18151 0.7212 0.909
+0.9042 0.2204 26199 0.7164 0.9573
+0.9315 0.1303 100440 0.7765 0.9547
+0.8944 0.1623 36441 0.6987 0.9323
+0.8259 0.557 25308 0.7302 0.9461
+0.779 0.2536 3222 0.6463 0.8241
+0.8698 0.2847 9971 0.752 0.9167
+0.9373 0.2999 3715 0.8285 0.9839
+0.8303 0.2903 5098 0.7155 0.8773
+0.7789 0.5583 5157 0.654 0.9368
+0.8622 0.4062 4476 0.7618 0.9308
+0.8472 0.472 12667 0.726 0.9556
+0.9585 0.3567 9605 0.8943 0.994
+0.9427 0.1757 10399 0.8232 0.9682
+0.9154 0.5007 61589 0.8437 0.9873
+0.8909 0.4094 29170 0.7988 0.9548
+0.8049 0.4936 4151 0.673 0.9334
+0.7078 0.6601 16178 0.5967 0.9234
+0.7867 0.4645 8316 0.6259 0.9261
+0.8496 0.5545 3690 0.74 0.986
+0.9188 0.2208 60426 0.8528 0.9376
+0.7722 0.1139 7896 0.5929 0.7952
+0.7426 0.4821 20782 0.5896 0.8851
+0.7822 0.7521 12377 0.7188 0.9746
+0.8612 0.2118 4684 0.6825 0.9093
+0.8699 0.3794 7440 0.6989 0.9744
+0.956 0.1608 4752 0.9359 0.9599
+0.7522 0.7545 14723 0.681 0.9712
+0.8458 0.4617 7550 0.7372 0.939
+0.8165 0.4665 3220 0.7024 0.9162
+0.7513 0.5451 11722 0.6111 0.9194
+0.9231 0.2628 7245 0.782 0.9734
+0.8559 0.6004 20230 0.7714 0.9828
+0.9113 0.2695 16472 0.855 0.9321
+0.824 0.5412 7733 0.7317 0.9329
+0.9153 0.3841 14989 0.8043 0.9846
+0.9205 0.1701 4279 0.7088 0.964
+0.9443 0.4087 4434 0.8808 0.9882
+0.8469 0.3428 5717 0.6337 0.9582
+0.9333 0.3176 18370 0.8315 0.9806
+0.892 0.5149 10182 0.8163 0.9723
+0.8913 0.4235 10444 0.7834 0.9706
+0.9336 0.2237 3058 0.8114 0.9688
+0.8358 0.3005 10576 0.6611 0.9108
+0.89 0.3348 4883 0.7309 0.9701
+0.8108 0.4197 13110 0.6487 0.928
+0.7913 0.4058 5170 0.6192 0.9089
+0.8954 0.3928 13146 0.7711 0.9758
+0.8147 0.5034 3335 0.701 0.93
+0.8083 0.3048 11620 0.6242 0.889
+0.8502 0.3063 11445 0.7233 0.9062
+0.8341 0.241 8332 0.6713 0.8858
+0.566 0.7216 5647 0.4402 0.8919
+0.7222 0.5683 14236 0.5765 0.9139
+0.7975 0.2004 7674 0.6664 0.8303
+0.8055 0.2523 18106 0.6718 0.8506
+0.8474 0.41 8807 0.7081 0.9442
+0.825 0.2809 15620 0.6433 0.8959
+0.8774 0.4885 44705 0.7821 0.9684
+0.747 0.5787 16612 0.6164 0.9264
+0.7759 0.4991 4664 0.6207 0.9307
+0.8509 0.3937 17895 0.7075 0.9441
+0.6805 0.8168 19701 0.6183 0.9579
+0.8042 0.5142 9705 0.6525 0.9648
+0.5961 0.7272 8031 0.4574 0.9658
+0.764 0.6677 5441 0.6802 0.9325
+0.8285 0.3212 7831 0.6394 0.918
+0.8555 0.2025 21574 0.6635 0.9042
+0.9239 0.0641 5400 0.763 0.9349
+0.7759 0.2706 3355 0.5264 0.8684
+0.8274 0.557 64585 0.7059 0.9802
+0.6943 0.677 3876 0.5987 0.8946
+0.8557 0.1805 9675 0.7314 0.8831
+0.8752 0.0684 11733 0.6214 0.8939
+0.8001 0.4965 17574 0.6669 0.9314
+0.6916 0.7155 6320 0.5927 0.9405
+0.821 0.4378 7535 0.6372 0.9641
+0.7934 0.413 5974 0.6737 0.8777
+0.8756 0.26 20841 0.7223 0.9295
+0.7816 0.3666 15111 0.6289 0.87
+0.84 0.243 13780 0.6581 0.8984
+0.6916 0.7371 8734 0.6028 0.9403
+0.7654 0.5532 20643 0.6258 0.9382
+0.7272 0.5755 5715 0.5908 0.9122
+0.8042 0.5319 11877 0.6676 0.9595
+0.8553 0.0682 2917 0.6432 0.8709
+0.7547 0.575 20586 0.6217 0.9345
+0.8526 0.3003 20478 0.6249 0.9504
+0.7588 0.422 13806 0.582 0.8878
+0.6872 0.6128 14344 0.5371 0.9247
+0.5461 0.7445 12429 0.3963 0.983
+0.8832 0.2375 6421 0.7416 0.9273
+0.734 0.6144 12549 0.6156 0.9227
+0.8214 0.311 2363 0.6803 0.8851
+0.7886 0.4236 7132 0.5647 0.9531
+0.7817 0.4941 16765 0.6173 0.9422
+0.7615 0.3987 17188 0.5621 0.8937
+0.8119 0.34 7645 0.6252 0.908
+0.8547 0.2878 26924 0.7095 0.9134
+0.832 0.0908 8248 0.6088 0.8542
+0.817 0.2233 12282 0.6276 0.8715
+0.9084 0.334 142942 0.78 0.9728
+0.7596 0.3175 3219 0.5587 0.853
+0.8106 0.6173 12004 0.707 0.9778
+0.8812 0.0885 11242 0.6543 0.9033
+0.7792 0.3925 12665 0.6019 0.8937
+0.6964 0.6178 13529 0.5257 0.9723
+0.8364 0.1522 20377 0.6125 0.8766
+0.8613 0.0743 7050 0.7424 0.8708
+0.849 0.1566 18465 0.6531 0.8854
+0.7386 0.7277 13459 0.6521 0.9697
+0.8843 0.1521 9449 0.7397 0.9103
+0.7109 0.6983 12700 0.603 0.9608
+0.7711 0.318 11118 0.6168 0.843
+0.7927 0.3391 7625 0.6292 0.8766
+0.7819 0.5048 14327 0.6286 0.9381
+0.7001 0.7221 16765 0.5943 0.9749
+0.7649 0.4765 7202 0.6233 0.8939
+0.8238 0.2798 21107 0.6475 0.8922
+0.6857 0.6811 11742 0.5544 0.9663
+0.8197 0.2683 4204 0.641 0.8852
+0.7999 0.5946 14949 0.7159 0.9229
+0.688 0.6382 8257 0.5565 0.92
+0.7879 0.4282 8858 0.6494 0.8916
+0.7109 0.694 9150 0.6044 0.9525
+0.7667 0.491 24493 0.6166 0.9116
+0.5747 0.8537 8529 0.5097 0.9543
+0.8266 0.6461 9252 0.7405 0.9838
+0.742 0.5904 5999 0.5807 0.9744
+0.8343 0.5209 17687 0.7296 0.9482
+0.7838 0.0574 3886 0.6816 0.7901
+0.804 0.5774 7613 0.7138 0.9273
+0.8146 0.7633 4665 0.7616 0.9855
+0.6883 0.6868 10733 0.5678 0.9527
+0.7616 0.3033 11399 0.5618 0.8485
+0.6567 0.5686 6367 0.4691 0.9039
+0.7263 0.5716 17651 0.566 0.9402
+0.6969 0.4023 5362 0.4316 0.8755
+0.8675 0.0912 5079 0.7495 0.8793
+0.7086 0.5334 15533 0.5412 0.9001
+0.7223 0.6479 14315 0.5931 0.9601
+0.8391 0.3703 13378 0.7148 0.9123
+0.6791 0.677 14286 0.5464 0.9573
+0.8276 0.4579 28077 0.6796 0.9526
+0.797 0.5052 13188 0.6566 0.9404
+0.7449 0.4598 7707 0.566 0.8972
+0.7085 0.601 13043 0.5677 0.9206
+0.824 0.1105 9623 0.7187 0.837
+0.8239 0.0515 6214 0.7594 0.8275
+0.8682 0.2726 7366 0.7037 0.9298
+0.7589 0.5062 14044 0.5665 0.9562
+0.827 0.284 14232 0.6964 0.8788
+0.6933 0.588 16179 0.5321 0.9235
+0.6306 0.7207 9904 0.5013 0.9642
+0.6408 0.7591 3146 0.5595 0.8971
+0.6733 0.6679 13356 0.5364 0.9486
+0.8703 0.4911 47837 0.7737 0.9635
+0.8353 0.3907 6337 0.6874 0.9301
+0.7965 0.6158 3675 0.6911 0.9653
+0.7066 0.5884 13833 0.5409 0.9436
+0.7735 0.5007 14511 0.6245 0.9229
+0.8716 0.229 6729 0.8034 0.8919
+0.7663 0.7493 9446 0.6955 0.9776
+0.7516 0.7185 21263 0.6643 0.9744
+0.6935 0.6802 8249 0.57 0.9564
+0.7275 0.7232 6264 0.6296 0.9833
+0.8544 0.3148 12603 0.684 0.9327
+0.7476 0.4775 7571 0.5544 0.9242
+0.8414 0.3564 9292 0.6739 0.9341
+0.6202 0.7396 15841 0.4974 0.969
+0.7719 0.5836 21134 0.6322 0.9677
+0.8237 0.3362 8060 0.6336 0.92
+0.7842 0.3132 7762 0.5117 0.9085
+0.7666 0.578 19145 0.6234 0.9627
+0.7758 0.4023 6913 0.5592 0.9216
+0.7084 0.6783 7566 0.6165 0.9022
+0.7878 0.5372 9069 0.6619 0.934
+0.8939 0.1328 13446 0.7766 0.9118
+0.7977 0.3897 17844 0.5877 0.9318
+0.8504 0.4015 16823 0.7296 0.9314
+0.7713 0.6664 8471 0.6771 0.9597
+0.7567 0.6104 19938 0.6391 0.9412
+0.841 0.2648 8946 0.6885 0.896
+0.7245 0.6655 4257 0.6149 0.9424
+0.8481 0.0736 3562 0.7481 0.8561
+0.874 0.1069 11472 0.7423 0.8898
+0.7741 0.4134 9370 0.5844 0.9078
+0.6253 0.6992 16796 0.4766 0.9709
+0.7497 0.4966 7116 0.5716 0.9255
+0.7861 0.4963 13386 0.6255 0.9444
+0.9076 0.3555 18260 0.8098 0.9615
+0.9383 0.0827 26550 0.8246 0.9486
+0.9313 0.1869 10880 0.7875 0.9643
+0.8815 0.3815 11753 0.7732 0.9483
+0.9568 0.0569 22638 0.7919 0.9667
+0.9595 0.071 95278 0.8586 0.9672
+0.9624 0.0915 24256 0.8622 0.9725
+0.9705 0.0582 24765 0.7989 0.9811
+0.9703 0.0899 221162 0.8734 0.9799
+0.967 0.0882 63809 0.8116 0.9821
+0.966 0.0966 19152 0.9665 0.966
+0.9773 0.1106 22739 0.8978 0.9872
+0.9822 0.052 19032 0.8556 0.9891
+0.9761 0.0943 20496 0.8473 0.9895
+0.9706 0.0736 6265 0.9284 0.974
+0.9772 0.0566 25036 0.8326 0.9859
+0.9697 0.1032 34389 0.8611 0.9822
+0.9908 0.0657 3165 0.9615 0.9929
+0.9806 0.0605 39158 0.9076 0.9853
+0.9757 0.0934 50947 0.8997 0.9835
+0.9807 0.0547 9923 0.8785 0.9866
+0.9601 0.1159 80378 0.884 0.97
+0.8323 0.0906 11866 0.6074 0.8547
+0.8161 0.0599 11152 0.6512 0.8266
+0.9191 0.0797 7650 0.7656 0.9324
+0.901 0.1251 9446 0.7073 0.9287
+0.8574 0.1418 18797 0.6597 0.89
+0.8459 0.0989 10431 0.6686 0.8654
+0.7912 0.1167 19870 0.7558 0.7959
+0.8736 0.3304 14028 0.699 0.9598
+0.8685 0.2894 11600 0.7066 0.9345
+0.8877 0.0792 15497 0.7663 0.8982
+0.882 0.0713 7078 0.6812 0.8975
+0.8887 0.1824 10639 0.7311 0.9239
+0.9101 0.0558 14432 0.7043 0.9223
+0.9325 0.0636 6388 0.7217 0.9468
+0.8136 0.4066 30099 0.6792 0.9058
+0.9041 0.257 14154 0.7295 0.9645
+0.8032 0.11 6920 0.7135 0.8143
+0.9213 0.1339 31644 0.7572 0.9467
+0.8762 0.3176 40077 0.7146 0.9514
+0.8912 0.0649 12333 0.6737 0.9063
+0.8667 0.1921 17125 0.6812 0.9109
+0.8724 0.2383 10745 0.6595 0.939
+0.8845 0.0594 3635 0.6759 0.8976
+0.9057 0.2013 8847 0.7962 0.9333
+0.9084 0.0861 24610 0.74 0.9243
+0.8444 0.1137 8645 0.6612 0.8679
+0.8904 0.0693 6439 0.6816 0.9059
+0.9009 0.0818 17076 0.7156 0.9173
+0.6868 0.0729 6966 0.4626 0.7044
+0.9214 0.1077 13264 0.722 0.9454
+0.878 0.1083 13164 0.7363 0.8952
+0.8955 0.2359 22707 0.722 0.9491
+0.878 0.1338 10839 0.6076 0.9198
+0.8792 0.1493 8741 0.6582 0.918
+0.8684 0.1976 26234 0.7454 0.8987
+0.9435 0.1831 217506 0.8033 0.9749
+0.8635 0.2429 9957 0.704 0.9147
+0.7924 0.052 15089 0.7513 0.7947
+0.8739 0.0697 7932 0.6582 0.8901
+0.865 0.1703 13338 0.7252 0.8936
+0.8888 0.0618 7831 0.719 0.9
+0.8537 0.2163 19083 0.6804 0.9016
+0.8597 0.206 7047 0.6894 0.9038
+0.9173 0.2355 27783 0.7707 0.9624
+0.8996 0.0578 9741 0.6643 0.914
+0.8711 0.2159 20168 0.7542 0.9033
+0.8468 0.1411 12020 0.5825 0.8902
+0.9247 0.1566 15080 0.7039 0.9656
+0.9143 0.0707 7268 0.7335 0.928
+0.9194 0.1601 10904 0.6999 0.9612
+0.8356 0.0748 7499 0.6257 0.8526
+0.8626 0.2581 10137 0.698 0.9198
+0.8517 0.106 20921 0.7384 0.8651
+0.8807 0.1776 12664 0.703 0.9191
+0.8815 0.0894 8275 0.6365 0.9055
+0.9147 0.1515 5922 0.6678 0.9588
+0.8925 0.0659 10985 0.6395 0.9103
+0.7976 0.0614 4427 0.6838 0.8051
+0.8846 0.2496 13482 0.7082 0.9433
+0.8791 0.2265 14144 0.6133 0.9569
+0.8608 0.1829 8922 0.6446 0.9092
+0.8653 0.1038 5577 0.6356 0.892
+0.8594 0.1228 8842 0.7099 0.8804
+0.8162 0.3233 12557 0.6746 0.8838
+0.8733 0.2272 10483 0.7632 0.9057
+0.9155 0.1288 14868 0.7347 0.9423
+0.8507 0.2049 23798 0.6336 0.9066
+0.8845 0.1286 10278 0.6853 0.9139
+0.9124 0.1342 15573 0.7923 0.931
+0.8405 0.3077 9806 0.6606 0.9205
+0.4982 0.2086 21422 0.278 0.5563
+0.6643 0.4912 17375 0.4882 0.8343
+0.5638 0.4379 17347 0.4922 0.6195
+0.617 0.3602 23248 0.4338 0.7201
+0.7556 0.4213 14877 0.4645 0.9675
+0.5574 0.7682 15754 0.4347 0.9641
+0.7225 0.614 44618 0.569 0.9667
+0.7852 0.2738 45178 0.6018 0.8543
+0.8096 0.3973 6044 0.6135 0.9388
+0.6104 0.1217 2934 0.3894 0.6411
+0.6143 0.5005 7125 0.3029 0.9264
+0.6798 0.5723 17037 0.4648 0.9673
+0.452 0.8465 10858 0.3666 0.9232
+0.6541 0.6366 19481 0.4795 0.9599
+0.7042 0.6144 26028 0.5544 0.9427
+0.536 0.892 9066 0.4847 0.9591
+0.6035 0.6886 14157 0.4623 0.9158
+0.6418 0.4533 8356 0.3205 0.9083
+0.7472 0.316 11157 0.4101 0.9029
+0.5671 0.4682 22383 0.4488 0.6712
+0.5472 0.6403 23150 0.4222 0.7696
+0.8233 0.2817 9367 0.562 0.9258
+0.6922 0.3918 13658 0.5313 0.7958
+0.8708 0.2163 6431 0.6751 0.9248
+0.4644 0.3716 19620 0.2012 0.62
+0.5878 0.2663 23164 0.4863 0.6247
+0.7835 0.3796 12871 0.5059 0.9534
+0.8015 0.1319 7273 0.5308 0.8426
+0.4931 0.8873 8302 0.4411 0.9028
+0.5461 0.7423 13563 0.404 0.9557
+0.6026 0.5561 25443 0.3981 0.8588
+0.931 0.2684 274998 0.8166 0.9729
+0.7653 0.5441 19388 0.6064 0.9549
+0.5715 0.5311 9041 0.3659 0.8044
+0.4757 0.6832 18318 0.3417 0.7648
+0.7331 0.4871 32222 0.5332 0.923
+0.6135 0.5483 7747 0.3543 0.9283
+0.5465 0.6667 11175 0.3607 0.9181
+0.8218 0.2146 13739 0.6072 0.8805
+0.6162 0.3788 3788 0.4599 0.7114
+0.5879 0.6133 7960 0.4898 0.7437
+0.7599 0.4841 6180 0.5284 0.9771
+0.6233 0.5924 16717 0.5117 0.7854
+0.6103 0.5706 10463 0.4496 0.8239
+0.4391 0.4694 44930 0.2865 0.5742
+0.4898 0.4282 15632 0.3701 0.5794
+0.6181 0.5561 28751 0.5135 0.749
+0.8286 0.3552 13754 0.6517 0.926
+0.7835 0.3228 21137 0.5752 0.8828
+0.4055 0.9132 13106 0.3522 0.9675
+0.6127 0.4124 20165 0.5492 0.6572
+0.8663 0.3464 13894 0.7039 0.9523
+0.5494 0.1725 17781 0.3828 0.5841
+0.8648 0.2345 12341 0.7357 0.9044
+0.835 0.2957 13232 0.6463 0.9143
+0.7635 0.4963 13436 0.5576 0.9665
+0.5148 0.7378 9344 0.3938 0.8551
+0.8662 0.4437 4424 0.7896 0.9273
+0.5106 0.8011 9812 0.4032 0.9431
+0.8403 0.2279 12814 0.6065 0.9093
+0.8617 0.3463 29797 0.7029 0.9458
+0.9488 0.1064 97157 0.7827 0.9686
+0.9555 0.1569 456882 0.8684 0.9717
+0.8105 0.4735 7463 0.674 0.9333
+0.8747 0.2368 14570 0.7116 0.9253
+0.96 0.0553 26840 0.7916 0.9698
+0.9371 0.1346 18910 0.7129 0.972
+0.7653 0.5008 11577 0.5899 0.9413
+0.8585 0.318 21887 0.726 0.9202
+0.943 0.0984 40824 0.7518 0.9639
+0.937 0.1784 22296 0.7549 0.9765
+0.9249 0.2217 12364 0.7618 0.9714
+0.867 0.3524 13390 0.6886 0.964
+0.9059 0.2739 24652 0.7247 0.9742
+0.8837 0.3001 26933 0.6949 0.9647
+0.839 0.337 12860 0.6357 0.9424
+0.7645 0.4171 11936 0.58 0.8965
+0.833 0.3442 20112 0.6539 0.9271
+0.8648 0.3333 15462 0.6666 0.9639
+0.8624 0.222 20672 0.7144 0.9047
+0.8012 0.3024 16691 0.5941 0.8909
+0.9712 0.0679 17004 0.8787 0.9779
+0.6973 0.7386 19421 0.597 0.9807
+0.8809 0.2388 13271 0.7706 0.9156
+0.7571 0.5369 15637 0.5906 0.9502
+0.8153 0.4477 19991 0.669 0.9339
+0.8427 0.4826 7221 0.7337 0.9443
+0.6758 0.8797 37924 0.6398 0.9389
+0.8794 0.2083 12136 0.6511 0.9395
+0.7867 0.5727 16009 0.6591 0.9577
+0.8365 0.5453 15778 0.725 0.9703
+0.867 0.2794 9919 0.6828 0.9384
+0.7172 0.7737 12645 0.6422 0.9738
+0.8472 0.4658 15140 0.7306 0.9488
+0.7775 0.6848 14776 0.681 0.9873
+0.7124 0.8834 26268 0.6773 0.9781
+0.7787 0.5377 25355 0.6208 0.9624
+0.8696 0.3166 11622 0.7201 0.9389
+0.7151 0.7538 16775 0.633 0.9666
+0.8701 0.3792 15205 0.7085 0.9689
+0.7283 0.4404 10271 0.4851 0.9196
+0.8504 0.2911 4559 0.7212 0.9035
+0.898 0.2387 4177 0.7422 0.9469
+0.8 0.7001 11224 0.7214 0.9834
+0.8116 0.3968 8215 0.7279 0.8666
+0.867 0.3189 26089 0.7184 0.9365
+0.7717 0.6927 47397 0.6827 0.9725
+0.7405 0.7881 28340 0.6759 0.9807
+0.5989 0.9393 8017 0.5756 0.9589
+0.8895 0.0835 10023 0.724 0.9045
+0.8517 0.3629 11599 0.7059 0.9348
+0.8156 0.4741 12666 0.6839 0.9342
+0.6824 0.7755 12901 0.5984 0.9727
+0.7958 0.5156 8537 0.6592 0.9412
+0.8619 0.29 20937 0.6983 0.9287
+0.7817 0.5623 14143 0.6642 0.9326
+0.8549 0.4383 15500 0.732 0.9508
+0.8247 0.3287 8258 0.6164 0.9266
+0.8206 0.4564 34958 0.6393 0.9727
+0.8325 0.3975 8818 0.6636 0.9439
+0.8569 0.3258 12274 0.6839 0.9404
+0.8407 0.3709 20574 0.6902 0.9295
+0.6789 0.8389 27142 0.6199 0.9863
+0.7997 0.4206 19827 0.6218 0.9289
+0.7244 0.6973 22904 0.6092 0.9899
+0.7253 0.8045 24257 0.6616 0.9871
+0.8535 0.3935 10774 0.7179 0.9415
+0.738 0.7023 19129 0.6436 0.9607
+0.7991 0.5512 25458 0.6767 0.9495
+0.8112 0.4899 12434 0.6552 0.9611
+0.8625 0.1631 11996 0.6111 0.9114
+0.8479 0.3798 15906 0.663 0.9611
+0.6632 0.8307 20557 0.5987 0.9796
+0.7367 0.6316 13908 0.6061 0.9606
+0.6963 0.6685 22308 0.5663 0.9585
+0.8678 0.2401 7214 0.6992 0.921
+0.8242 0.358 5394 0.6991 0.894
+0.8119 0.4635 26344 0.6328 0.9667
+0.8763 0.2418 13792 0.7238 0.9249
+0.8551 0.1715 11887 0.6675 0.894
+0.7066 0.7734 8715 0.6341 0.9539
+0.7587 0.5829 16537 0.6096 0.9671
+0.8172 0.3978 11547 0.6377 0.9357
+0.691 0.8951 11856 0.6585 0.9686
+0.8117 0.3452 11632 0.61 0.9181
+0.8608 0.1748 11098 0.6392 0.9077
+0.6799 0.817 21531 0.6165 0.963
+0.674 0.7015 21145 0.5529 0.9585
+0.766 0.5776 14124 0.6204 0.9651
+0.8949 0.1799 10158 0.6924 0.9393
+0.8538 0.0889 9093 0.6386 0.8748
+0.6088 0.9093 14451 0.573 0.968
+0.8586 0.2187 13358 0.6571 0.915
+0.7765 0.6956 29630 0.6821 0.9921
+0.6396 0.8488 38287 0.5867 0.9366
+0.8141 0.4006 10209 0.6447 0.9273
+0.9002 0.2135 10375 0.6745 0.9615
+0.5899 0.7589 12461 0.4682 0.973
+0.8246 0.3803 11620 0.6135 0.9542
+0.8358 0.5057 15370 0.6982 0.9767
+0.7392 0.756 34182 0.665 0.9691
+0.9696 0.0739 17559 0.7404 0.9879
+0.9336 0.1397 24115 0.7201 0.9682
+0.973 0.0508 76767 0.8539 0.9794
+0.9198 0.069 15032 0.8862 0.9223
+0.9511 0.1413 19481 0.7751 0.98
+0.9415 0.0738 20912 0.7738 0.9549
+0.9561 0.096 18387 0.762 0.9767
+0.9611 0.053 16512 0.7554 0.9726
+0.9564 0.109 17659 0.8265 0.9723
+0.9496 0.1428 16152 0.8005 0.9744
+0.9099 0.2037 12555 0.6535 0.9755
+0.9758 0.0928 240901 0.8971 0.9839
+0.9421 0.0579 21466 0.7804 0.952
+0.9609 0.0994 23829 0.7965 0.9791
+0.9525 0.0785 13356 0.7235 0.972
+0.9601 0.1 24794 0.7701 0.9812
+0.9284 0.1472 10623 0.7807 0.9539
+0.9696 0.0673 14689 0.729 0.987
+0.9585 0.0734 12401 0.7264 0.9769
+0.8658 0.1158 14159 0.6543 0.8935
+0.8828 0.0878 14190 0.7697 0.8937
+0.964 0.0844 27330 0.7977 0.9793
+0.9538 0.147 18329 0.7769 0.9843
+0.9748 0.0513 11409 0.8068 0.9838
+0.9656 0.0937 21135 0.8141 0.9813
+0.9515 0.0564 16608 0.794 0.9609
+0.9645 0.0681 19432 0.7846 0.9777
+0.9713 0.0527 64617 0.7946 0.9811
+0.9631 0.0676 572262 0.876 0.9694
+0.9511 0.1263 23521 0.7993 0.973
+0.9537 0.0514 7211 0.6792 0.9686
+0.941 0.1622 59044 0.9196 0.9451
+0.943 0.0523 54084 0.8614 0.9475
+0.9592 0.0674 113938 0.8557 0.9667
+0.9467 0.0733 16015 0.8995 0.9504
+0.9601 0.0622 30474 0.914 0.9631
+0.9589 0.088 78181 0.9009 0.9645
+0.9589 0.1222 21771 0.844 0.9749
+0.8493 0.2394 21042 0.6665 0.9068
+0.8378 0.0757 8161 0.6699 0.8515
+0.7783 0.5093 17933 0.6581 0.9032
+0.7951 0.4154 22128 0.6615 0.8901
+0.8086 0.568 16236 0.7238 0.92
+0.747 0.4532 12761 0.591 0.8762
+0.7664 0.365 10125 0.6388 0.8398
+0.9041 0.1669 36853 0.8262 0.9197
+0.8208 0.1193 15407 0.6703 0.8412
+0.8431 0.2304 18730 0.7022 0.8853
+0.8135 0.1149 14422 0.6626 0.8331
+0.7924 0.3778 4052 0.659 0.8735
+0.8334 0.1654 10177 0.6744 0.865
+0.739 0.4938 10604 0.5903 0.8839
+0.8825 0.1206 19608 0.7961 0.8944
+0.821 0.3249 16315 0.6796 0.8891
+0.8137 0.5401 8069 0.7451 0.8944
+0.819 0.1893 20897 0.6595 0.8562
+0.7698 0.3104 19846 0.6886 0.8063
+0.7903 0.5582 18895 0.6865 0.9213
+0.8058 0.4197 25521 0.6919 0.8883
+0.8548 0.3228 5523 0.7173 0.9203
+0.9206 0.105 3515 0.8889 0.9243
+0.8436 0.1301 21258 0.6753 0.8688
+0.8121 0.166 9720 0.6561 0.8432
+0.8008 0.3509 17981 0.696 0.8575
+0.8371 0.3518 26828 0.6883 0.9179
+0.7253 0.5916 22790 0.6145 0.8858
+0.8305 0.3049 35847 0.693 0.8908
+0.7873 0.4595 17444 0.6712 0.886
+0.8133 0.2242 26250 0.6669 0.8556
+0.8624 0.4336 7560 0.7916 0.9166
+0.7863 0.4719 17787 0.6597 0.8994
+0.7669 0.4571 8933 0.6699 0.8487
+0.8808 0.2516 45016 0.758 0.922
+0.7571 0.6231 26846 0.6621 0.914
+0.8198 0.2761 15442 0.6751 0.8749
+0.8727 0.1156 11798 0.7141 0.8934
+0.7785 0.5718 10955 0.6908 0.8958
+0.7478 0.3925 6328 0.5419 0.8809
+0.8815 0.2126 24503 0.7207 0.9249
+0.801 0.2438 28777 0.6796 0.8401
+0.7891 0.4501 5928 0.6848 0.8745
+0.8537 0.3015 8264 0.7191 0.9118
+0.7959 0.4365 16170 0.6852 0.8817
+0.859 0.1604 12092 0.7758 0.8749
+0.8756 0.156 9603 0.7323 0.902
+0.7325 0.4837 12514 0.6027 0.854
+0.8624 0.3785 50047 0.7116 0.9542
+0.8046 0.2388 10346 0.6402 0.8561
+0.8358 0.3191 12196 0.7027 0.8982
+0.7632 0.4055 23739 0.6393 0.8477
+0.8479 0.4778 25049 0.7192 0.9656
+0.7475 0.5642 15903 0.6452 0.88
+0.8169 0.2885 9776 0.706 0.8618
+0.8689 0.3105 11049 0.7604 0.9178
+0.8291 0.361 7003 0.7215 0.8898
+0.836 0.4908 12127 0.729 0.9391
+0.7981 0.4765 10981 0.6661 0.9182
+0.8573 0.4908 7767 0.7825 0.9295
+0.7624 0.4166 12265 0.6 0.8784
+0.7408 0.4866 24956 0.5844 0.8891
+0.8331 0.1467 5380 0.7364 0.8497
+0.868 0.1099 21199 0.7382 0.884
+0.8146 0.4539 14114 0.7126 0.8993
+0.7288 0.4176 37004 0.6342 0.7967
+0.8145 0.2823 26139 0.6547 0.8773
+0.878 0.2491 27223 0.7311 0.9268
+0.8238 0.1484 19903 0.6895 0.8472
+0.8095 0.3224 21101 0.7114 0.8562
+0.7001 0.5275 10953 0.6007 0.811
+0.815 0.1122 13910 0.5964 0.8426
+0.7691 0.1285 13731 0.5896 0.7955
+0.791 0.0917 20814 0.5928 0.811
+0.8685 0.089 5132 0.7133 0.8836
+0.7958 0.2975 3775 0.6643 0.8514
+0.815 0.2791 22896 0.6296 0.8868
+0.8094 0.5037 13872 0.7017 0.9187
+0.8339 0.3989 47083 0.731 0.9022
+0.8057 0.6385 14256 0.7414 0.9193
+0.859 0.4919 7868 0.7747 0.9407
+0.8063 0.4283 25701 0.6753 0.9044
+0.7828 0.0856 21226 0.7822 0.7829
+0.7502 0.4255 20363 0.5997 0.8617
+0.8054 0.0716 10893 0.6179 0.8198
+0.9676 0.0527 20129 0.82 0.9758
+0.9785 0.0539 21689 0.9333 0.9811
+0.9792 0.0848 54885 0.9102 0.9855
+0.9582 0.0566 17402 0.8569 0.9643
+0.9702 0.0655 185700 0.9087 0.9745
+0.9333 0.0721 20353 0.8515 0.9396
+0.9682 0.1336 24459 0.8721 0.983
+0.9711 0.0548 387373 0.8535 0.9779
+0.9458 0.0585 31807 0.8806 0.9499
+0.8649 0.1604 9736 0.7382 0.8892
+0.9379 0.0819 13150 0.8124 0.9491
+0.917 0.0746 20949 0.6993 0.9345
+0.9152 0.1717 17888 0.7695 0.9454
+0.8477 0.0635 11655 0.7986 0.851
+0.8835 0.2004 15362 0.7212 0.9242
+0.9054 0.0622 11135 0.7633 0.9149
+0.96 0.07 12602 0.8673 0.967
+0.9562 0.1026 19354 0.8413 0.9693
+0.9312 0.0844 18324 0.7731 0.9458
+0.9715 0.0614 21916 0.8648 0.9785
+0.8891 0.0705 16607 0.7541 0.8993
+0.905 0.0576 11503 0.6767 0.9189
+0.97 0.1281 14149 0.8533 0.9872
+0.9234 0.0576 7731 0.8112 0.9303
+0.8812 0.0658 19919 0.7879 0.8878
+0.9627 0.107 24962 0.8558 0.9755
+0.944 0.2459 24223 0.8315 0.9807
+0.9074 0.0946 7007 0.7391 0.925
+0.9519 0.0719 10843 0.809 0.9629
+0.8063 0.2305 14168 0.7183 0.8327
+0.8623 0.2541 14492 0.8232 0.8756
+0.9081 0.0595 9442 0.7509 0.918
+0.9334 0.3127 40189 0.8793 0.9581
+0.9455 0.1461 10519 0.7879 0.9725
+0.8926 0.4097 13857 0.8358 0.932
+0.979 0.1061 69249 0.9065 0.9876
+0.9144 0.28 15215 0.8138 0.9535
+0.9751 0.0596 17534 0.8632 0.9822
+0.9073 0.1195 34279 0.8057 0.9211
+0.8567 0.2051 13650 0.7504 0.8841
+0.8426 0.1344 17020 0.7644 0.8547
+0.9707 0.082 27034 0.9175 0.9755
+0.9042 0.3976 15864 0.8343 0.9504
+0.9514 0.0949 87999 0.8944 0.9574
+0.9508 0.1016 96446 0.8843 0.9583
+0.9542 0.0571 1261132 0.9223 0.9562
+0.766 0.6339 24797 0.6601 0.9493
+0.7587 0.5347 30233 0.6416 0.8932
+0.7899 0.3763 48963 0.6315 0.8855
+0.7636 0.6753 12917 0.6701 0.958
+0.6956 0.7063 24417 0.5921 0.9445
+0.6228 0.8637 22398 0.5702 0.9561
+0.6362 0.7813 16678 0.5635 0.8956
+0.7294 0.7482 11496 0.6513 0.9613
+0.7884 0.6219 69560 0.6731 0.978
+0.7605 0.3172 18274 0.5943 0.8378
+0.6964 0.6368 21071 0.5588 0.9378
+0.7106 0.3946 18266 0.5646 0.8057
+0.7209 0.7149 21886 0.6405 0.9224
+0.5986 0.6269 25052 0.4178 0.9025
+0.7219 0.5756 25283 0.6364 0.8378
+0.7504 0.4965 15943 0.6728 0.827
+0.7379 0.6072 12635 0.6112 0.9337
+0.7062 0.6917 19764 0.5903 0.9664
+0.6501 0.7412 20473 0.5529 0.9285
+0.7592 0.5571 25032 0.641 0.908
+0.6021 0.7187 16142 0.4968 0.8711
+0.8215 0.3021 49157 0.6666 0.8885
+0.7725 0.61 24508 0.6649 0.9408
+0.6815 0.6346 17421 0.5455 0.9177
+0.7591 0.2398 18694 0.5889 0.8128
+0.6741 0.595 18783 0.5075 0.9188
+0.7141 0.4839 18401 0.5325 0.8844
+0.7399 0.5372 29578 0.5845 0.9203
+0.6221 0.6601 17236 0.4808 0.8967
+0.8152 0.3524 22800 0.6082 0.9278
+0.8083 0.526 14598 0.7374 0.8868
+0.6768 0.5879 21831 0.5775 0.8185
+0.7973 0.6158 24564 0.6939 0.963
+0.8073 0.2495 19168 0.691 0.846
+0.7761 0.6447 39223 0.6733 0.9626
+0.8024 0.2143 17587 0.6148 0.8535
+0.8056 0.5265 41900 0.6892 0.9351
+0.7793 0.5208 14394 0.625 0.9469
+0.7872 0.3134 59419 0.6331 0.8575
+0.758 0.7167 27138 0.6794 0.9568
+0.7344 0.5001 20758 0.6025 0.8663
+0.6706 0.6019 25660 0.546 0.8589
+0.727 0.5164 34063 0.5624 0.9027
+0.8524 0.0565 12584 0.782 0.8566
+0.899 0.2431 17245 0.7412 0.9497
+0.8527 0.0633 4421 0.6929 0.8636
+0.8911 0.0619 14890 0.7959 0.8973
+0.8836 0.106 12074 0.7453 0.9
+0.8402 0.079 18929 0.7298 0.8497
+0.8297 0.0552 7734 0.548 0.8462
+0.8887 0.2062 17791 0.7081 0.9356
+0.879 0.1477 7660 0.7515 0.9011
+0.8674 0.1692 6620 0.6384 0.914
+0.7938 0.0562 13529 0.7224 0.798
+0.8814 0.1029 11560 0.6773 0.9048
+0.8663 0.221 11601 0.6568 0.9257
+0.9019 0.3227 120823 0.7602 0.9694
+0.8613 0.101 7100 0.7406 0.8748
+0.8357 0.0566 11025 0.6795 0.8451
+0.877 0.1533 14171 0.744 0.9011
+0.8514 0.208 20103 0.6948 0.8925
+0.6996 0.7291 21159 0.6097 0.9414
+0.8719 0.1563 14857 0.6968 0.9043
+0.8946 0.2275 31030 0.7363 0.9412
+0.8265 0.3255 23835 0.6522 0.9106
+0.8929 0.1193 10207 0.7545 0.9117
+0.9091 0.3026 70051 0.802 0.9555
+0.7701 0.4271 16765 0.6032 0.8944
+0.8615 0.124 12422 0.763 0.8754
+0.8206 0.0838 16960 0.6749 0.834
+0.7973 0.669 18697 0.7194 0.9548
+0.8702 0.1112 12315 0.7516 0.8851
+0.8861 0.2295 18773 0.6866 0.9455
+0.8619 0.1479 11558 0.6953 0.8908
+0.9126 0.1472 4579 0.8472 0.9239
+0.8507 0.0883 9780 0.6343 0.8717
+0.8415 0.0956 3785 0.5608 0.8712
+0.8556 0.0937 12950 0.6993 0.8718
+0.9117 0.1432 73016 0.8322 0.925
+0.7483 0.3915 6416 0.6246 0.8279
+0.8293 0.4518 15300 0.6946 0.9403
+0.8632 0.0566 12399 0.6752 0.8745
+0.8173 0.1512 4193 0.7192 0.8348
+0.8526 0.2097 19056 0.6927 0.895
+0.8992 0.0728 9889 0.8014 0.9069
+0.8802 0.0923 15161 0.7186 0.8967
+0.878 0.0939 11679 0.7311 0.8932
+0.8753 0.4007 29777 0.7473 0.9608
+0.8663 0.1329 13406 0.7872 0.8785
+0.9006 0.1935 12873 0.7077 0.9468
+0.8185 0.3982 30535 0.6601 0.9233
+0.851 0.0864 4236 0.7268 0.8628
+0.8373 0.0566 14748 0.7542 0.8423
+0.8291 0.399 25279 0.7122 0.9067
+0.8826 0.0704 3492 0.6992 0.8965
+0.8639 0.0822 8131 0.6527 0.8828
+0.9154 0.1809 22405 0.7705 0.9474
+0.8326 0.0691 6118 0.747 0.839
+0.847 0.088 11043 0.7459 0.8567
+0.8769 0.1054 16363 0.7587 0.8908
+0.8505 0.2536 19196 0.6728 0.9109
+0.8395 0.3408 24906 0.6704 0.9269
+0.8879 0.4762 155103 0.7837 0.9825
+0.8607 0.1219 13426 0.6644 0.888
+0.8316 0.1197 10490 0.6457 0.8568
+0.8693 0.0585 20658 0.7444 0.877
+0.8484 0.2097 19310 0.6621 0.8978
+0.8037 0.4509 21158 0.6589 0.9226
+0.8167 0.2976 4358 0.6816 0.8739
+0.8915 0.1189 12300 0.6904 0.9186
+0.8999 0.0808 21992 0.7373 0.9142
+0.8496 0.0757 8331 0.71 0.861
+0.9118 0.1082 23865 0.7456 0.932
+0.831 0.0615 10955 0.6261 0.8445
+0.7976 0.3237 18017 0.6116 0.8866
+0.8806 0.2414 19503 0.6782 0.945
+0.8852 0.3815 21599 0.735 0.9778
+0.9019 0.1466 12350 0.7189 0.9334
+0.9772 0.071 1620 0.8174 0.9894
+0.8903 0.2708 13131 0.7213 0.9531
+0.8488 0.3647 18445 0.7584 0.9007
+0.9501 0.1318 36070 0.8111 0.9712
+0.8771 0.1018 94021 0.8644 0.8785
+0.9229 0.0797 3138 0.652 0.9463
+0.871 0.3646 25262 0.7184 0.9586
+0.8336 0.4579 10309 0.6609 0.9796
+0.8032 0.4684 13929 0.7102 0.8852
+0.84 0.4641 13517 0.7373 0.9289
+0.7444 0.2183 17126 0.7614 0.7396
+0.9507 0.1319 2737 0.8698 0.963
+0.8556 0.4496 6793 0.7109 0.9738
+0.8532 0.3507 19013 0.6855 0.9438
+0.9254 0.2443 3135 0.7963 0.9671
+0.9115 0.2568 20751 0.7429 0.9698
+0.8545 0.3698 14167 0.6982 0.9462
+0.9547 0.0701 19367 0.8143 0.9652
+0.9502 0.1863 107583 0.8235 0.9792
+0.8346 0.2059 17214 0.7774 0.8494
+0.9481 0.0567 10031 0.7909 0.9575
+0.9436 0.0717 22696 0.7529 0.9583
+0.9352 0.1825 38652 0.7907 0.9675
+0.8704 0.3573 25403 0.7598 0.9319
+0.9372 0.122 32488 0.7427 0.9643
+0.9108 0.2366 21810 0.7702 0.9544
+0.7071 0.6336 13806 0.5807 0.9257
+0.9201 0.0833 6494 0.6118 0.9481
+0.8897 0.4191 14527 0.7608 0.9827
+0.9443 0.2052 36259 0.8694 0.9636
+0.7355 0.151 7086 0.7832 0.7271
+0.815 0.2858 19966 0.7365 0.8465
+0.9491 0.1225 49570 0.7738 0.9735
+0.8424 0.5407 10352 0.7263 0.9792
+0.8122 0.461 15289 0.6862 0.92
+0.779 0.2219 17847 0.7677 0.7822
+0.9127 0.2123 9528 0.8028 0.9423
+0.9316 0.274 93247 0.8196 0.9739
+0.8158 0.6205 26929 0.7219 0.9694
+0.8275 0.142 11395 0.8084 0.8306
+0.9063 0.2072 14121 0.7344 0.9512
+0.9439 0.1075 33318 0.7508 0.9671
+0.9248 0.1066 21801 0.6902 0.9528
+0.8874 0.4158 20513 0.7844 0.9608
+0.9479 0.0973 34549 0.7503 0.9692
+0.8824 0.3175 4737 0.7547 0.9418
+0.8773 0.3388 10089 0.7332 0.9511
+0.9231 0.2879 30220 0.8252 0.9627
+0.7481 0.0557 10444 0.7577 0.7475
+0.8988 0.2348 26157 0.7257 0.9519
+0.9499 0.0513 3451 0.8927 0.953
+0.8978 0.2324 33609 0.7186 0.952
+0.922 0.1598 19161 0.7593 0.953
+0.9093 0.3038 9195 0.7834 0.9642
+0.8472 0.4008 11673 0.6858 0.9551
+0.8127 0.3208 7945 0.601 0.9127
+0.9117 0.2734 24452 0.7867 0.9587
+0.9319 0.2403 54615 0.8047 0.9721
+0.7381 0.0526 798 0.6429 0.7434
+0.897 0.2586 7242 0.6946 0.9676
+0.8349 0.6267 7474 0.7541 0.9706
+0.9077 0.3239 10234 0.7768 0.9704
+0.8951 0.2557 26153 0.7585 0.942
+0.8392 0.4472 11231 0.6952 0.9557
+0.8986 0.3411 7308 0.7485 0.9763
+0.8867 0.2497 19101 0.701 0.9485
+0.9237 0.2358 34322 0.756 0.9754
+0.9001 0.3453 7546 0.7985 0.9536
+0.8939 0.2116 7088 0.7533 0.9316
+0.8596 0.42 13970 0.7081 0.9693
+0.834 0.3797 12470 0.6969 0.9179
+0.9041 0.0935 4653 0.7195 0.9232
+0.8772 0.2981 20152 0.7172 0.9451
+0.7944 0.1697 2009 0.651 0.8237
+0.8098 0.5361 20069 0.6961 0.941
+0.9428 0.0966 5787 0.6351 0.9757
+0.8686 0.4081 18661 0.7218 0.9699
+0.914 0.1864 5869 0.745 0.9527
+0.8444 0.3012 7685 0.6592 0.9242
+0.8964 0.5389 6702 0.8378 0.965
+0.9061 0.196 18308 0.6898 0.9588
+0.8773 0.4019 29985 0.7341 0.9736
+0.9629 0.1481 85409 0.8667 0.9796
+0.8953 0.1921 11521 0.653 0.9529
+0.8912 0.2715 43159 0.7894 0.9291
+0.9076 0.2533 8989 0.8103 0.9406
+0.8846 0.2085 7294 0.7055 0.9318
+0.8833 0.2772 13836 0.6787 0.9618
+0.9447 0.0595 18008 0.7817 0.9551
+0.8254 0.2395 11171 0.7054 0.8632
+0.8379 0.5079 11578 0.7307 0.9486
+0.8118 0.5423 9001 0.6779 0.9704
+0.8547 0.4628 18716 0.765 0.932
+0.8443 0.4263 15519 0.72 0.9365
+0.9009 0.1759 30546 0.7205 0.9395
+0.818 0.056 11972 0.7791 0.8203
+0.9107 0.1702 16420 0.7206 0.9497
+0.8123 0.3369 27207 0.6578 0.8909
+0.8588 0.3315 27467 0.6932 0.9409
+0.9354 0.2737 20479 0.8139 0.9812
+0.8991 0.2147 14684 0.7523 0.9393
+0.7574 0.6146 6472 0.6247 0.9691
+0.7956 0.3832 13506 0.5777 0.931
+0.8962 0.3423 6433 0.7997 0.9463
+0.9175 0.1635 33337 0.7545 0.9494
+0.903 0.1925 4905 0.7818 0.9318
+0.8412 0.2809 21533 0.6317 0.923
+0.9013 0.1966 12968 0.762 0.9354
+0.7938 0.5774 13551 0.7096 0.9089
+0.7823 0.4914 11082 0.6702 0.8907
+0.8706 0.3494 40176 0.7003 0.9621
+0.7626 0.5143 11924 0.6443 0.8878
+0.748 0.6914 3781 0.7043 0.8458
+0.7882 0.5194 11359 0.6492 0.9385
+0.8524 0.3424 15828 0.6651 0.9499
+0.8638 0.2516 5625 0.6947 0.9207
+0.8579 0.3757 10075 0.6867 0.9609
+0.7551 0.6347 6528 0.6259 0.9795
+0.8413 0.5045 30508 0.7224 0.9624
+0.9059 0.3649 17670 0.7766 0.9801
+0.8101 0.5724 6626 0.7269 0.9216
+0.9195 0.2303 15755 0.7696 0.9644
+0.8727 0.3225 16870 0.7076 0.9513
+0.8807 0.0595 10056 0.6923 0.8926
+0.819 0.3879 6254 0.5812 0.9697
+0.7884 0.2057 18621 0.5717 0.8445
+0.9305 0.094 14624 0.7125 0.9531
+0.8616 0.0647 8251 0.6142 0.8787
+0.7961 0.462 9068 0.7262 0.8561
+0.7701 0.5523 6763 0.6337 0.9382
+0.6921 0.1902 4842 0.6037 0.7128
+0.7495 0.6107 8423 0.6297 0.9375
+0.7783 0.4962 28411 0.6598 0.8949
+0.8397 0.3962 12848 0.6704 0.9508
+0.9154 0.3585 122886 0.7946 0.9828
+0.782 0.3953 13019 0.607 0.8965
+0.741 0.4828 10780 0.5739 0.897
+0.7973 0.4494 5047 0.597 0.9608
+0.7802 0.5481 6948 0.6686 0.9156
+0.8003 0.4465 4762 0.6463 0.9245
+0.8457 0.5501 6273 0.7476 0.9656
+0.7583 0.5108 7220 0.6003 0.9233
+0.897 0.2361 16406 0.6783 0.9646
+0.8383 0.4624 12115 0.7051 0.9529
+0.8106 0.5166 9259 0.6787 0.9515
+0.7964 0.3152 7274 0.642 0.8675
+0.9108 0.2699 6870 0.7546 0.9685
+0.7977 0.5575 20755 0.695 0.9272
+0.8412 0.5122 6587 0.7389 0.9486
+0.882 0.1454 15796 0.714 0.9106
+0.7756 0.5647 20090 0.6355 0.9572
+0.7847 0.3051 11868 0.5852 0.8723
+0.7672 0.5781 3484 0.64 0.9415
+0.8796 0.4449 121680 0.7645 0.9718
+0.8343 0.5338 12245 0.7272 0.9569
+0.8567 0.3778 7915 0.7258 0.9362
+0.8658 0.5285 10127 0.7599 0.9845
+0.8574 0.3968 9990 0.698 0.9622
+0.8919 0.0832 10366 0.6833 0.9108
+0.8191 0.0955 11617 0.7286 0.8286
+0.7823 0.3735 51029 0.5881 0.8981
+0.7972 0.5716 4507 0.6506 0.9927
+0.8422 0.5781 10629 0.7409 0.981
+0.7793 0.5724 5606 0.6482 0.9549
+0.8703 0.2227 8958 0.7088 0.9166
+0.7887 0.4841 8568 0.6661 0.9038
+0.8455 0.1792 12508 0.6983 0.8776
+0.7753 0.262 5927 0.6948 0.8038
+0.6914 0.4045 5422 0.4355 0.8653
+0.9216 0.2162 42155 0.7636 0.9651
+0.8763 0.1678 17999 0.7192 0.9079
+0.9015 0.067 26205 0.6838 0.9171
+0.7191 0.592 18625 0.6058 0.8835
+0.8732 0.3123 11932 0.7013 0.9513
+0.8438 0.2081 6064 0.687 0.885
+0.7791 0.601 7048 0.6667 0.9484
+0.7503 0.6343 9866 0.6206 0.9751
+0.8555 0.1206 17711 0.7551 0.8692
+0.8973 0.1254 6484 0.6716 0.9296
+0.9038 0.4462 20420 0.8054 0.983
+0.8565 0.0926 28244 0.7683 0.8655
+0.7747 0.4804 6782 0.6025 0.9339
+0.8251 0.0958 23548 0.6906 0.8393
+0.8814 0.1129 14816 0.7327 0.9003
+0.8735 0.4814 5733 0.8207 0.9226
+0.9305 0.0556 11731 0.7531 0.941
+0.9453 0.0841 17159 0.7782 0.9607
+0.9346 0.0574 35863 0.8601 0.9391
+0.8827 0.193 37352 0.804 0.9015
+0.9227 0.0735 18247 0.7547 0.9361
+0.9233 0.2153 12311 0.8045 0.9559
+0.9026 0.0868 59693 0.8529 0.9073
+0.7982 0.3301 35604 0.7592 0.8174
+0.9 0.1664 27932 0.7732 0.9253
+0.8042 0.0777 13986 0.6532 0.817
+0.9463 0.0511 9718 0.7626 0.9562
+0.8955 0.0897 18108 0.7968 0.9052
+0.9294 0.0669 13263 0.8343 0.9362
+0.9818 0.0683 21722 0.9501 0.9841
diff --git a/data/reg.txt b/data/reg.txt
new file mode 100644
index 0000000..603cbcc
--- /dev/null
+++ b/data/reg.txt
@@ -0,0 +1,276 @@
+"X" "Y" "N"
+0.2475 0.8875 40000
+0.18181818 0.95454545 4400
+0.13550136 0.80758808 36900
+0.1884058 0.85507246 6900
+0.072625698 0.97206704 89500
+0.11002445 0.96821516 245400
+0.13636364 0.95454545 4400
+0.028571429 1 14000
+0.086021505 0.97849462 9300
+0.10569106 0.97560976 12300
+0.10769231 0.95384615 13000
+0.2421875 0.953125 12800
+0.1233279 0.7548124 613000
+0.17105263 0.72368421 7600
+0.12121212 1 3300
+0.22273074 0.78032037 262200
+0.19007392 0.84899683 94700
+0.28571429 0.85714286 2800
+0.18181818 0.93181818 4400
+0.5125 0.48333333 24000
+0.095238095 0.95238095 2100
+0.38888889 0.72222222 1800
+0.23636364 0.87272727 5500
+0.36363636 0.88636364 4400
+0.069444444 0.77777778 7200
+0.24489796 0.91836735 4900
+0.16 0.96 7500
+0.16197183 0.94366197 14200
+0.12688698 0.84210526 245100
+0.03030303 1 6600
+0.16352201 0.94339623 15900
+0.27692308 0.78974359 19500
+0.50980392 0.78431373 5100
+0.1 0.95 2000
+0.17344173 0.72086721 36900
+0.10375276 0.9602649 45300
+0.29146919 0.87203791 42200
+0.25806452 0.83870968 6200
+0.090909091 1 2200
+0.4025974 0.87012987 7700
+0.11134454 0.8592437 47600
+0.2971246 0.87859425 31300
+0.13953488 0.96899225 12900
+0.10211268 0.67957746 28400
+0.22826087 0.89130435 9200
+0.063829787 0.86322188 32900
+0.125 0.975 4000
+0.13659148 0.77318296 159600
+0.085271318 0.96899225 12900
+0.1887396 0.93090211 156300
+0.068181818 0.96306818 35200
+0.066179159 0.92431444 273500
+0.16480687 0.79828326 116500
+0.28021978 0.72527473 18200
+0.23756906 0.7679558 18100
+0.2742616 0.89873418 23700
+0.06875 0.98125 16000
+0.073476703 0.96057348 55800
+0.73899371 0.29874214 95400
+0.21126761 0.94366197 7100
+0.24418605 0.89534884 8600
+0.22972973 0.91891892 7400
+0.275 0.575 4000
+0.13519814 0.93822844 85800
+0.22222222 0.97222222 3600
+0.11702128 0.96808511 9400
+0.15625 0.953125 6400
+0.17037037 0.91851852 27000
+0.21495327 0.95327103 10700
+0.2745098 0.94117647 15300
+0.35164835 0.85714286 9100
+0.22926829 0.83902439 20500
+0.19444444 0.92592593 10800
+0.42268041 0.75257732 9700
+0.22295082 0.61311475 30500
+0.32200772 0.58378378 129500
+0.19127086 0.75224647 77900
+0.21818182 0.90909091 5500
+0.076923077 0.97435897 3900
+0.28787879 0.86363636 6600
+0.43859649 0.71052632 11400
+0.39316239 0.87179487 11700
+0.5112782 0.66165414 13300
+0.29543635 0.82786229 124900
+0.58333333 0.45833333 7200
+0.46564885 0.29770992 13100
+0.2173913 1 18400
+0.32835821 0.70149254 13400
+0.19277108 0.92771084 8300
+0.26277372 0.92335766 27400
+0.44936709 0.90506329 15800
+0.27472527 0.87912088 9100
+0.12541806 0.95986622 119600
+0.1835443 0.81012658 15800
+0.21111111 0.93111111 45000
+0.10652921 0.96907216 29100
+0.10126582 0.98734177 7900
+0.37254902 0.74509804 15300
+0.10169492 0.98870056 17700
+0.58426966 0.7752809 8900
+0.41142857 0.65714286 17500
+0.39572193 0.85026738 18700
+0.32820245 0.58361597 383300
+0.28975265 0.58833922 56600
+0.25217391 0.73043478 11500
+0.46491228 0.90350877 11400
+0.28775835 0.71701113 62900
+0.37288136 0.88135593 5900
+0.37704918 0.68852459 12200
+0.20192308 0.93269231 10400
+0.046218487 0.98739496 23800
+0.2231405 1 12100
+0.4375 0.97916667 4800
+0.43956044 0.92307692 9100
+0.43434343 0.96969697 9900
+0.36945813 0.89408867 40600
+0.30718954 0.9869281 15300
+0.27906977 0.9379845 25800
+0.17605634 0.95422535 28400
+0.29652997 0.87697161 31700
+0.59322034 0.74576271 5900
+0.18394649 0.80602007 29900
+0.2970297 0.88118812 10100
+0.11353712 1 22900
+0.098360656 0.96721311 12200
+0.28813559 0.8559322 23600
+0.30837004 0.77092511 22700
+0.43209877 0.87654321 8100
+0.18421053 0.76315789 7600
+0.62162162 0.2972973 7400
+0.27659574 0.87234043 9400
+0.14741036 0.74302789 50200
+0.054945055 0.98901099 9100
+0.021276596 0.55319149 4700
+0.4 0.66923077 13000
+0.008264463 0.99173554 12100
+0.015151515 0.48484848 6600
+0.31155779 0.73869347 19900
+0.48461538 0.9 13000
+0.35664336 0.68531469 14300
+0.2962963 0.87037037 10800
+0.10532838 0.70384139 80700
+0.060509554 0.77707006 31400
+0.13300493 0.66995074 40600
+0.0625 0.79779412 27200
+0.35483871 0.41935484 3100
+0.10614525 0.37988827 17900
+0.40594059 0.45544554 10100
+0.078947368 0.89473684 41800
+0.26315789 0.55263158 15200
+0.021505376 0.47311828 9300
+0.3968254 0.66666667 6300
+0.17615176 0.77777778 36900
+0.29365079 0.42857143 25200
+0.26282051 0.52564103 31200
+0.24383917 0.44487678 77100
+0.28205128 0.82051282 3900
+0.054054054 0.54054054 3700
+0.097826087 0.95217391 46000
+0.11 0.86 10000
+0.31531532 0.82882883 22200
+0.29279279 0.83033033 66600
+0.44244604 0.6618705 27800
+0.2228164 0.84402852 112200
+0.36363636 0.68181818 15400
+0.11570248 0.73002755 72600
+0.46 0.64 5000
+0.029411765 0.97058824 3400
+0.37634409 0.75268817 18600
+0.40740741 0.72839506 8100
+0.18958333 0.71319444 144000
+0.45544554 0.7359736 30300
+0.23574144 0.72623574 26300
+0.021186441 0.59322034 23600
+0.053811659 0.88340807 22300
+0.52136752 0.76923077 11700
+0.48051948 0.66233766 7700
+0.33333333 0.66666667 3300
+0.15027322 0.83879781 36600
+0.08 0.91 10000
+0.18443804 0.86455331 34700
+0.41818182 0.70909091 5500
+0.18918919 0.87162162 14800
+0.34797297 0.64527027 29600
+0.091463415 0.89634146 16400
+0.051612903 0.76774194 15500
+0.022727273 0.54545455 8800
+0.010309278 0.77319588 9700
+0.41304348 0.78985507 13800
+0.2164557 0.7943038 158000
+0.19811321 0.94339623 10600
+0.23414634 0.64390244 20500
+0.32716049 0.7191358 32400
+0.25118483 0.73459716 42200
+0.52898551 0.88405797 13800
+0.12820513 0.37435897 39000
+0.20491803 0.88934426 24400
+0.30188679 0.69811321 5300
+0.34265734 0.61538462 14300
+0.42268041 0.65979381 9700
+0.39215686 0.39215686 5100
+0.29577465 0.86619718 14200
+0.37569061 0.72375691 36200
+0.11594203 0.56521739 6900
+0.072222222 0.91388889 36000
+0.25581395 0.85581395 21500
+0.50591017 0.53900709 42300
+0.18137255 0.68872549 40800
+0.14371257 0.7005988 50100
+0.097744361 0.52255639 26600
+0.3203125 0.421875 25600
+0.376 0.704 12500
+0.090534979 0.81481481 24300
+0.078125 0.6484375 12800
+0.04964539 0.74822695 28200
+0.17021277 0.72340426 4700
+0.043956044 0.97802198 9100
+0.33333333 0.625 2400
+0.18032787 0.73770492 24400
+0.37142857 0.86857143 17500
+0.22968907 0.78335005 99700
+0.55555556 0.57575758 9900
+0.37142857 0.84285714 7000
+0.010309278 0.72164948 9700
+0.35033259 0.51884701 45100
+0.07751938 0.82170543 25800
+0.34394904 0.63694268 31400
+0.044117647 0.91176471 13600
+0.0125 0.4375 8000
+0.26890756 0.60504202 11900
+0.2293578 0.72706422 43600
+0.56140351 0.73684211 5700
+0.16812609 0.5061296 57100
+0.46341463 0.70731707 8200
+0.35955056 0.86516854 8900
+0.37305699 0.51295337 19300
+0.40641711 0.9197861 18700
+0.55932203 0.61016949 5900
+0.31305115 0.61463845 113400
+0.17525773 0.71134021 19400
+0.33727811 0.64497041 16900
+0.30057803 0.82080925 17300
+0.59689922 0.75968992 12900
+0.43055556 0.72916667 14400
+0.37218045 0.67669173 26600
+0.38732394 0.66901408 14200
+0.37762238 0.93706294 14300
+0.48101266 0.59493671 7900
+0.52380952 0.73333333 10500
+0.37209302 0.74418605 43000
+0.44720497 0.7515528 16100
+0.15371901 0.60165289 121000
+0.26153846 0.66153846 26000
+0.46590909 0.75 8800
+0.21203438 0.64756447 34900
+0.55 0.66666667 6000
+0.34302326 0.72674419 17200
+0.22857143 0.78571429 21000
+0.2556391 0.59022556 26600
+0.55102041 0.68367347 9800
+0.14285714 0.875 33600
+0.53658537 0.68292683 4100
+0.48734177 0.59493671 15800
+0.41843972 0.73049645 14100
+0.29069767 0.70348837 17200
+0.1 0.60454545 22000
+0.5147929 0.81065089 33800
+0.090909091 0.625 26400
+0.29274843 0.61862131 111700
+0.29113924 0.7721519 7900
+0.18826135 0.65891473 90300
+0.41176471 0.52673797 37400
+0.24260355 0.76331361 16900
+0.5801105 0.72375691 18100
+0.24285714 0.57142857 42000
diff --git a/eco.pdf b/eco.pdf
new file mode 100644
index 0000000..56b09ce
Binary files /dev/null and b/eco.pdf differ
diff --git a/man/census.Rd b/man/census.Rd
new file mode 100644
index 0000000..2394606
--- /dev/null
+++ b/man/census.Rd
@@ -0,0 +1,39 @@
+\name{census}
+
+\docType{data}
+
+\alias{census}
+
+\title{Black Illiteracy Rates in 1910 US Census}
+
+\description{
+ This data set contains the proportion of the residents who are black,
+ the proportion of those who can read, the total population as well as
+ the actual black literacy rate and white literacy rate for 1040
+ counties in the US. The dataset was originally analyzed by Robison
+ (1950) at the state level. King (1997) recoded the 1910 census at
+ county level. The data set only includes those who are older than 10
+ years of age.
+}
+
+\usage{data(census)}
+
+\format{A data frame containing 5 variables and 1040 observations
+ \tabular{lll}{
+ X \tab numeric \tab the proportion of Black residents in each county\cr
+ Y \tab numeric \tab the overall literacy rates in each county\cr
+ N \tab numeric \tab the total number of residents in each county \cr
+ W1 \tab numeric \tab the actual Black literacy rate \cr
+ W2 \tab numeric \tab the actual White literacy rate
+ }
+}
+\references{
+ Robinson, W.S. (1950). ``Ecological Correlations and the Behavior
+ of Individuals.'' \emph{American Sociological Review}, vol. 15,
+ pp.351-357. \cr \cr
+ King, G. (1997). \dQuote{A Solution to the Ecological
+ Inference Problem: Reconstructing Individual Behavior from Aggregate
+ Data}. Princeton University Press, Princeton, NJ.
+}
+
+\keyword{datasets}
diff --git a/man/eco.Rd b/man/eco.Rd
new file mode 100644
index 0000000..c47c559
--- /dev/null
+++ b/man/eco.Rd
@@ -0,0 +1,212 @@
+\name{eco}
+
+\alias{eco}
+
+\title{Fitting the Parametric Bayesian Model of Ecological Inference in
+ 2x2 Tables}
+
+\description{
+ \code{eco} is used to fit the parametric Bayesian model (based on a
+ Normal/Inverse-Wishart prior) for ecological inference in \eqn{2
+ \times 2} tables via Markov chain Monte Carlo. It gives the in-sample
+ predictions as well as the estimates of the model parameters. The
+ model and algorithm are described in Imai and Lu (2004). The
+ contextual effect can also be modeled by following the strategy
+ described in Imai and Lu (2005).
+}
+
+\usage{
+eco(formula, data = parent.frame(), N = NULL, supplement = NULL,
+ context = FALSE, mu0 = 0, tau0 = 2, nu0 = 4, S0 = 10,
+ mu.start = 0, Sigma.start = 10, parameter = TRUE,
+ grid = FALSE, n.draws = 5000, burnin = 0, thin = 0,
+ verbose = FALSE)
+}
+
+\arguments{
+ \item{formula}{A symbolic description of the model to be fit,
+ specifying the column and row margins of \eqn{2 \times
+ 2} ecological tables. \code{Y ~ X} specifies \code{Y} as the
+ column margin and \code{X} as the row margin. Details and specific
+ examples are given below.
+ }
+ \item{data}{An optional data frame in which to interpret the variables
+ in \code{formula}. The default is the environment in which
+ \code{eco} is called.
+ }
+ \item{N}{An optional variable representing the size of the unit; e.g.,
+ the total number of voters.}
+ \item{supplement}{An optional matrix of supplemental data. The matrix
+ has two columns, which contain additional individual-level data such
+ as survey data for \eqn{W_1} and \eqn{W_2}, respectively. If
+ \code{NULL}, no additional individual-level data are included in the
+ model. The default is \code{NULL}.
+ }
+ \item{context}{Logical. If \code{TRUE}, the contextual effect is also
+ modeled. See Imai and Lu (2005) for details. The default is
+ \code{FALSE}.
+ }
+ \item{mu0}{A scalar or a numeric vector that specifies the prior mean
+ for the mean parameter \eqn{\mu}. If it is a scalar, then its value
+ will be repeated to yield a vector of the length of \eqn{\mu}, otherwise,
+ it needs to be a vector of same length as \eqn{\mu}.
+ When \code{context=TRUE}, the length of \eqn{\mu} is 3,
+ otherwise it is 2. The default is \code{0}.
+ }
+ \item{tau0}{A positive integer representing the prior scale
+ for the mean parameter \eqn{\mu}. The default is \code{2}.
+ }
+ \item{nu0}{A positive integer representing the prior degrees of
+ freedom of the variance matrix \eqn{\Sigma}. the default is \code{4}.
+ }
+ \item{S0}{A postive scalar or a positive definite matrix that specifies
+ the prior scale matrix for the variance matrix \eqn{\Sigma}. If it is
+ a scalar, then the prior scale matrix will be a digonal matrix with
+ the same dimensions as \eqn{\Sigma} and the diagonal elements all take value
+ of \code{S0}, otherwise \code{S0} needs to have same dimensions as
+ \eqn{\Sigma}. When \code{context=TRUE}, \eqn{\Sigma} is a
+ \eqn{3 \times 3} matrix, otherwise, it is \eqn{2 \times 2}.
+ The default is \code{10}.
+ }
+ \item{mu.start}{A scalar or a numeric vector that specifies the
+ starting values of the mean parameter \eqn{\mu}.
+ If it is a scalar, then its value will be repeated to
+ yield a vector of the length of \eqn{\mu}, otherwise,
+ it needs to be a vector of same length as \eqn{\mu}.
+ When \code{context=FALSE}, the length of \eqn{\mu} is 2,
+ otherwise it is 3. The default is \code{0}.
+ }
+ \item{Sigma.start}{A scalar or a positive definite matrix
+ that specified the starting value of the variance matrix
+ \eqn{\Sigma}. If it is a scalar, then the prior scale
+ matrix will be a digonal matrix with the same dimensions
+ as \eqn{\Sigma} and the diagonal elements all take value
+ of \code{S0}, otherwise \code{S0} needs to have same dimensions as
+ \eqn{\Sigma}. When \code{context=TRUE}, \eqn{\Sigma} is a
+ \eqn{3 \times 3} matrix, otherwise, it is \eqn{2 \times 2}.
+ The default is \code{10}.
+ }
+ \item{parameter}{Logical. If \code{TRUE}, the Gibbs draws of the population
+ parameters, \eqn{\mu} and \eqn{\Sigma}, are returned in addition to
+ the in-sample predictions of the missing internal cells,
+ \eqn{W}. The default is \code{TRUE}.
+ }
+ \item{grid}{Logical. If \code{TRUE}, the grid method is used to sample
+ \eqn{W} in the Gibbs sampler. If \code{FALSE}, the Metropolis
+ algorithm is used where candidate draws are sampled from the uniform
+ distribution on the tomography line for each unit. Note that the
+ grid method is significantly slower than the Metropolis algorithm.
+ }
+ \item{n.draws}{A positive integer. The number of MCMC draws.
+ The default is \code{5000}.
+ }
+ \item{burnin}{A positive integer. The burnin interval for the Markov
+ chain; i.e. the number of initial draws that should not be stored. The
+ default is \code{0}.
+ }
+ \item{thin}{A positive integer. The thinning interval for the
+ Markov chain; i.e. the number of Gibbs draws between the recorded
+ values that are skipped. The default is \code{0}.
+ }
+ \item{verbose}{Logical. If \code{TRUE}, the progress of the Gibbs
+ sampler is printed to the screen. The default is \code{FALSE}.
+ }
+}
+
+\details{
+ An example of \eqn{2 \times 2} ecological table for racial voting is
+ given below:
+ \tabular{lccc}{
+ \tab black voters \tab white voters \tab \cr
+ Voted \tab \eqn{W_{1i}} \tab \eqn{W_{2i}} \tab \eqn{Y_i} \cr
+ Not voted \tab \eqn{1-W_{1i}} \tab \eqn{1-W_{2i}} \tab \eqn{1-Y_i} \cr
+ \tab \eqn{X_i} \tab \eqn{1-X_i} \tab
+ }
+ where \eqn{Y_i} and \eqn{X_i} represent the observed margins, and
+ \eqn{W_1} and \eqn{W_2} are unknown variables. All variables are
+ proportions and hence bounded between 0 and 1. For each \eqn{i}, the
+ following deterministic relationship holds,
+ \eqn{Y_i=X_i W_{1i}+(1-X_i)W_{2i}}.
+}
+
+\examples{
+
+## load the registration data
+data(reg)
+
+## NOTE: convergence has not been properly assessed for the following
+## examples. See Imai and Lu (2004, 2005) for more complete analyses.
+
+## fit the parametric model with the default prior specification
+res <- eco(Y ~ X, data = reg, verbose = TRUE)
+## summarize the results
+summary(res)
+
+## obtain out-of-sample prediction
+out <- predict(res, verbose = TRUE)
+## summarize the results
+summary(out)
+
+## load the Robinson's census data
+data(census)
+
+## fit the parametric model with contextual effects and N
+## using the default prior specification
+res1 <- eco(Y ~ X, N = N, context = TRUE, data = census, verbose = TRUE)
+## summarize the results
+summary(res1)
+
+## obtain out-of-sample prediction
+out1 <- predict(res1, verbose = TRUE)
+## summarize the results
+summary(out1)
+
+}
+
+\value{
+ An object of class \code{eco} containing the following elements:
+ \item{call}{The matched call.}
+ \item{X}{The row margin, \eqn{X}.}
+ \item{Y}{The column margin, \eqn{Y}.}
+ \item{N}{The size of each table, \eqn{N}.}
+ \item{burnin}{The number of initial burnin draws.}
+ \item{thin}{The thinning interval.}
+ \item{nu0}{The prior degrees of freedom.}
+ \item{tau0}{The prior scale parameter.}
+ \item{mu0}{The prior mean.}
+ \item{S0}{The prior scale matrix.}
+ \item{W}{A three dimensional array storing the posterior in-sample
+ predictions of \eqn{W}. The first dimension indexes the Monte Carlo
+ draws, the second dimension indexes the columns of the table, and the
+ third dimension represents the observations.}
+ \item{Wmin}{A numeric matrix storing the lower bounds of \eqn{W}.}
+ \item{Wmax}{A numeric matrix storing the upper bounds of \eqn{W}.}
+ The following additional elements are included in the output when
+ \code{parameter = TRUE}.
+ \item{mu}{The posterior draws of the population mean parameter,
+ \eqn{\mu}.}
+ \item{Sigma}{The posterior draws of the population variance matrix,
+ \eqn{\Sigma}.}
+}
+
+\author{
+ Kosuke Imai, Department of Politics, Princeton University
+ \email{kimai at Princeton.Edu}, \url{http://imai.princeton.edu};
+ Ying Lu, Institute for Quantitative Social Sciences,
+ Harvard University \email{ylu at Latte.Harvard.Edu}}
+
+\references{
+ Imai, Kosuke and Ying Lu. (2004) \dQuote{ Parametric and Nonparametric
+ Bayesian Models for Ecological Inference in \eqn{2 \times 2}
+ Tables.} Working Paper, Princeton University,
+ available at
+ \url{http://imai.princeton.edu/research/einonpar.html}
+
+ Imai, Kosuke and Ying Lu. (2005) \dQuote{An Incomplete Data Approach
+ to Ecological Inference.} Working Paper, Princeton University,
+ available at
+ \url{http://imai.princeton.edu/research/coarse.html} }
+}
+
+\seealso{\code{ecoNP}, \code{predict.eco}, \code{summary.eco}}
+\keyword{models}
diff --git a/man/ecoBD.Rd b/man/ecoBD.Rd
new file mode 100644
index 0000000..67ee96a
--- /dev/null
+++ b/man/ecoBD.Rd
@@ -0,0 +1,127 @@
+\name{ecoBD}
+
+\alias{ecoBD}
+
+\title{Calculating the Bounds for Ecological Inference in RxC Tables}
+
+\description{
+ \code{ecoBD} is used to calculate the bounds for missing internal cells
+ of \eqn{R \times C} ecological table. The data can be entered either
+ in the form of counts or proportions.
+}
+
+\usage{
+ecoBD(formula, data = parent.frame(), N = NULL)
+}
+
+\arguments{
+ \item{formula}{A symbolic description of ecological table to be used,
+ specifying the column and row margins of \eqn{R \times
+ C} ecological tables.
+ Details and specific
+ examples are given below.
+ }
+ \item{data}{An optional data frame in which to interpret the variables
+ in \code{formula}. The default is the environment in which
+ \code{ecoBD} is called.
+ }
+ \item{N}{An optional variable representing the size of the unit; e.g.,
+ the total number of voters. If \code{formula} is entered as counts
+ and the last row and/or column is omitted, this input is necessary.}
+}
+
+\details{
+ The data may be entered either in the form of counts or proportions.
+ If proportions are used, \code{formula} may omit the last row and/or
+ column of tables, which can be calculated from the remaining margins.
+ For example, \code{Y ~ X} specifies \code{Y} as the first column
+ margin and \code{X} as the first row margin in \eqn{2 \times 2} tables.
+ If counts are used, \code{formula} may omit the last row and/or column
+ margin of the table only if \code{N} is supplied. For larger tables,
+ one can use \code{cbind()} and \code{+}. For example, \code{cbind(Y1,
+ Y2, Y3) ~ X1 + X2 + X3 + X4)} specifies \eqn{3 \times 4} tables.
+
+ An \eqn{R \times C} ecological table in the form of counts:
+ \tabular{lcccc}{
+ \eqn{n_{i11}} \tab \eqn{n_{i12}} \tab \dots \tab \eqn{n_{i1C}} \tab \eqn{n_{i1.}} \cr
+ \eqn{n_{i21}} \tab \eqn{n_{i22}} \tab \dots \tab \eqn{n_{i2C}} \tab \eqn{n_{i2.}} \cr
+ \dots \tab \dots \tab \dots \tab \dots \tab \dots\cr
+ \eqn{n_{iR1}} \tab \eqn{n_{iR2}} \tab \dots \tab \eqn{n_{iRC}} \tab \eqn{n_{iR.}} \cr
+ \eqn{n_{i.1}} \tab \eqn{n_{i.2}} \tab \dots \tab \eqn{n_{i.C}} \tab \eqn{N_i}
+ }
+ where \eqn{n_{nr.}} and \eqn{n_{i.c}} represent the observed margins,
+ \eqn{N_i} represents the size of the table, and \eqn{n_{irc}} are
+ unknown variables. Note that for each \eqn{i}, the following
+ deterministic relationships hold; \eqn{n_{ir.} = \sum_{c=1}^C
+ n_{irc}} for \eqn{r=1,\dots,R}, and \eqn{n_{i.c}=\sum_{r=1}^R n_{irc}}
+ for \eqn{c=1,\dots,C}. Then, each of the unknown inner cells can be
+ bounded in the following manner,
+ \deqn{\max(0, n_{ir.}+n_{i.c}-N_i) \le n_{irc} \le \min(n_{ir.}, n_{i.c}).}
+ If the size of tables, \code{N}, is provided,
+
+ An \eqn{R \times C} ecological table in the form of proportions:
+ \tabular{lcccc}{
+ \eqn{W_{i11}} \tab \eqn{W_{i12}} \tab \dots \tab \eqn{W_{i1C}} \tab \eqn{Y_{i1}} \cr
+ \eqn{W_{i21}} \tab \eqn{W_{i22}} \tab \dots \tab \eqn{W_{i2C}} \tab \eqn{Y_{i2}} \cr
+ \dots \tab \dots \tab \dots \tab \dots \tab \dots \cr
+ \eqn{W_{iR1}} \tab \eqn{W_{iR2}} \tab \dots \tab \eqn{W_{iRC}} \tab \eqn{Y_{iR}} \cr
+ \eqn{X_{i1}} \tab \eqn{X_{i2}} \tab \dots \tab \eqn{X_{iC}} \tab
+ }
+ where \eqn{Y_{ir}} and \eqn{X_{ic}} represent the observed margins,
+ and \eqn{W_{irc}} are unknown variables. Note that for each \eqn{i},
+ the following deterministic relationships hold; \eqn{Y_{ir} =
+ \sum_{c=1}^C X_{ic} W_{irc}} for \eqn{r=1,\dots,R}, and \eqn{\sum_{r=1}^R
+ W_{irc}=1} for \eqn{c=1,\dots,C}. Then, each of
+ the inner cells of the table can be bounded in the following manner,
+ \deqn{\max(0, (X_{ic} + Y_{ir}-1)/X_{ic}) \le W_{irc}
+ \le \min(1, Y_{ir}/X_{ir}).}
+}
+
+\examples{
+
+## load the registration data
+data(reg)
+
+## calculate the bounds
+res <- ecoBD(Y ~ X, N = N, data = reg)
+## print the results
+print(res)
+}
+
+\value{
+ An object of class \code{ecoBD} containing the following elements
+ (When three dimensional arrays are used, the first dimension indexes
+ the observations, the second dimension indexes the row numbers, and
+ the third dimension indexes the column numbers):
+ \item{call}{The matched call.}
+ \item{X}{A matrix of the observed row margin, \eqn{X}.}
+ \item{Y}{A matrix of the observed column margin, \eqn{Y}.}
+ \item{N}{A vector of the size of ecological tables, \eqn{N}.}
+ \item{aggWmin}{A three dimensional array of aggregate lower bounds for
+ proportions.}
+ \item{aggWmax}{A three dimensional array of aggregate upper bounds for
+ proportions.}
+ \item{Wmin}{A three dimensional array of lower bounds for
+ proportions.}
+ \item{Wmax}{A three dimensional array of upper bounds for
+ proportions.}
+ \item{Nmin}{A three dimensional array of lower bounds for
+ counts.}
+ \item{Nmax}{A three dimensional array of upper bounds for
+ counts.}
+ The object can be printed through \code{print.ecoBD}.
+}
+
+\author{
+ Kosuke Imai, Department of Politics, Princeton University
+ \email{kimai at Princeton.Edu}, \url{http://www.princeton.edu/~kimai};
+ Ying Lu, Institute for Quantitative Social Sciences,
+ Harvard University \email{ylu at Latte.Harvard.Edu}}
+
+\references{
+ Imai, Kosuke. (2005) \dQuote{Ecological Inference in \eqn{R \times C}
+ Tables} Working Paper, Princeton University.
+}
+
+\seealso{\code{eco}, \code{ecoNP}}
+\keyword{models}
diff --git a/man/ecoNP.Rd b/man/ecoNP.Rd
new file mode 100644
index 0000000..b730821
--- /dev/null
+++ b/man/ecoNP.Rd
@@ -0,0 +1,222 @@
+\name{ecoNP}
+
+\alias{ecoNP}
+
+\title{Fitting the Nonparametric Bayesian Model of Ecological Inference
+ in 2x2 Tables}
+
+\description{
+ \code{ecoNP} is used to fit the nonparametric Bayesian model (based on
+ a Dirichlet process prior) for ecological inference in \eqn{2 \times
+ 2} tables via Markov chain Monte Carlo. It gives the in-sample
+ predictions as well as out-of-sample predictions for population
+ inference. The model and algorithm are described in Imai and Lu
+ (2004). The contextual effect can also be modeled by following the
+ strategy described in Imai and Lu (2005).
+}
+
+\usage{
+ecoNP(formula, data = parent.frame(), N = NULL, supplement = NULL,
+ context = FALSE, mu0 = 0, tau0 = 2, nu0 = 4, S0 = 10,
+ alpha = NULL, a0 = 1, b0 = 0.1, parameter = FALSE,
+ grid = FALSE, n.draws = 5000, burnin = 0, thin = 0,
+ verbose = FALSE)
+}
+
+\arguments{
+ \item{formula}{A symbolic description of the model to be fit,
+ specifying the column and row margins of \eqn{2 \times
+ 2} ecological tables. \code{Y ~ X} specifies \code{Y} as the
+ column margin and \code{X} as the row margin. Details and specific
+ examples are given below.
+ }
+ \item{data}{An optional data frame in which to interpret the variables
+ in \code{formula}. The default is the environment in which
+ \code{ecoNP} is called.
+ }
+ \item{N}{An optional variable representing the size of the unit; e.g.,
+ the total number of voters.}
+ \item{supplement}{An optional matrix of supplemental data. The matrix
+ has two columns, which contain additional individual-level data such
+ as survey data for \eqn{W_1} and \eqn{W_2}, respectively. If
+ \code{NULL}, no additional individual-level data are included in the
+ model. The default is \code{NULL}.
+ }
+ \item{context}{Logical. If \code{TRUE}, the contextual effect is also
+ modeled. See Imai and Lu (2005) for details. The default is
+ \code{FALSE}.
+ }
+ \item{mu0}{A scalar or a numeric vector that specifies the prior mean
+ for the mean parameter \eqn{\mu}. If it is a scalar, then its value
+ will be repeated to yield a vector of the length of \eqn{\mu}, otherwise,
+ it needs to be a vector of same length as \eqn{\mu}.
+ When \code{context=TRUE }, the length of \eqn{\mu} is 3,
+ otherwise it is 2. The default is \code{0}.
+ }
+ \item{tau0}{A positive integer representing the prior scale
+ for the mean parameter \eqn{\mu}. The default is \code{2}.
+ }
+ \item{nu0}{A positive integer representing the prior degrees of
+ freedom of the variance matrix \eqn{\Sigma}. the default is \code{4}.
+ }
+ \item{S0}{A postive scalar or a positive definite matrix that specifies
+ the prior scale matrix for the variance matrix \eqn{\Sigma}. If it is
+ a scalar, then the prior scale matrix will be a digonal matrix with
+ the same dimensions as \eqn{\Sigma} and the diagonal elements all take value
+ of \code{S0}, otherwise \code{S0} needs to have same dimensions as
+ \eqn{\Sigma}. When \code{context=TRUE}, \eqn{\Sigma} is a
+ \eqn{3 \times 3} matrix, otherwise, it is \eqn{2 \times 2}.
+ The default is \code{10}.
+ }
+ \item{alpha}{A positive scalar representing a user-specified fixed
+ value of the concentration parameter, \eqn{\alpha}. If \code{NULL},
+ \eqn{\alpha} will be updated at each Gibbs draw, and its prior
+ parameters \code{a0} and \code{b0} need to be specified. The default
+ is \code{NULL}.
+ }
+ \item{a0}{A positive integer representing the value of shape parameter
+ of the gamma prior distribution for \eqn{\alpha}. The default is \code{1}.
+ }
+ \item{b0}{A positive integer representing the value of the scale
+ parameter of the gamma prior distribution for \eqn{\alpha}. The
+ default is \code{0.1}.
+ }
+ \item{parameter}{Logical. If \code{TRUE}, the Gibbs draws of the population
+ parameters, \eqn{\mu} and \eqn{\Sigma}, are returned in addition to
+ the in-sample predictions of the missing internal cells,
+ \eqn{W}. The default is \code{FALSE}. This needs to be set to
+ \code{TRUE} if one wishes to make population inferences through
+ \code{predict.eco}. See an example below.
+ }
+ \item{grid}{Logical. If \code{TRUE}, the grid method is used to sample
+ \eqn{W} in the Gibbs sampler. If \code{FALSE}, the Metropolis
+ algorithm is used where candidate draws are sampled from the uniform
+ distribution on the tomography line for each unit. Note that the
+ grid method is significantly slower than the Metropolis algorithm.
+ }
+ \item{n.draws}{A positive integer. The number of MCMC draws.
+ The default is \code{5000}.
+ }
+ \item{burnin}{A positive integer. The burnin interval for the Markov
+ chain; i.e. the number of initial draws that should not be stored. The
+ default is \code{0}.
+ }
+ \item{thin}{A positive integer. The thinning interval for the
+ Markov chain; i.e. the number of Gibbs draws between the recorded
+ values that are skipped. The default is \code{0}.
+ }
+ \item{verbose}{Logical. If \code{TRUE}, the progress of the gibbs
+ sampler is printed to the screen. The default is \code{FALSE}.
+ }
+}
+
+\details{
+ An example of \eqn{2 \times 2} ecological table for racial voting is
+ given below:
+ \tabular{lccc}{
+ \tab black voters \tab white voters \tab \cr
+ Voted \tab \eqn{W_{1i}} \tab \eqn{W_{2i}} \tab \eqn{Y_i} \cr
+ Not voted \tab \eqn{1-W_{1i}} \tab \eqn{1-W_{2i}} \tab \eqn{1-Y_i} \cr
+ \tab \eqn{X_i} \tab \eqn{1-X_i} \tab
+ }
+ where \eqn{Y_i} and \eqn{X_i} represent the observed margins, and
+ \eqn{W_1} and \eqn{W_2} are unknown variables. All variables are
+ proportions and hence bounded between 0 and 1. For each \eqn{i}, the
+ following deterministic relationship holds,
+ \eqn{Y_i=X W_{1i}+(1-X_i)W_{2i}}.
+}
+
+\examples{
+
+## load the registration data
+data(reg)
+
+## NOTE: We set the number of MCMC draws to be a very small number in
+## the following examples; i.e., convergence has not been properly
+## assessed. See Imai and Lu (2004, 2005) for more complete examples.
+
+## fit the nonparametric model to give in-sample predictions
+## store the parameters to make population inference later
+res <- ecoNP(Y ~ X, data = reg, n.draws = 50, param = TRUE, verbose = TRUE)
+##summarize the results
+summary(res)
+
+## obtain out-of-sample prediction
+out <- predict(res, verbose = TRUE)
+## summarize the results
+summary(out)
+
+## density plots of the out-of-sample predictions
+par(mfrow=c(2,1))
+plot(density(out[,1]), main = "W1")
+plot(density(out[,2]), main = "W2")
+
+
+## load the Robinson's census data
+data(census)
+
+## fit the parametric model with contextual effects and N
+## using the default prior specification
+res1 <- ecoNP(Y ~ X, N = N, context = TRUE, param = TRUE, data = census,
+ n.draws = 25, verbose = TRUE)
+## summarize the results
+summary(res1)
+
+## out-of sample prediction
+pres1 <- predict(res1)
+summary(pres1)
+}
+
+\value{
+ An object of class \code{ecoNP} containing the following elements:
+ \item{call}{The matched call.}
+ \item{X}{The row margin, \eqn{X}.}
+ \item{Y}{The column margin, \eqn{Y}.}
+ \item{burnin}{The number of initial burnin draws.}
+ \item{thin}{The thinning interval.}
+ \item{nu0}{The prior degrees of freedom.}
+ \item{tau0}{The prior scale parameter.}
+ \item{mu0}{The prior mean.}
+ \item{S0}{The prior scale matrix.}
+ \item{a0}{The prior shape parameter.}
+ \item{b0}{The prior scale parameter.}
+ \item{W}{A three dimensional array storing the posterior in-sample
+ predictions of \eqn{W}. The first dimension indexes the Monte Carlo
+ draws, the second dimension indexes the columns of the table, and the
+ third dimension represents the observations.}
+ \item{Wmin}{A numeric matrix storing the lower bounds of \eqn{W}.}
+ \item{Wmax}{A numeric matrix storing the upper bounds of \eqn{W}.}
+ The following additional elements are included in the output when
+ \code{parameter = TRUE}.
+ \item{mu}{A three dimensional array storing the posterior draws of the
+ population mean parameter, \eqn{\mu}. The first dimension indexes the
+ Monte Carlo draws, the second dimension indexes the columns of the
+ table, and the third dimension represents the observations.}
+ \item{Sigma}{A three dimensional array storing the posterior draws of
+ the population variance matrix, \eqn{\Sigma}. The first dimension
+ indexes the Monte Carlo draws, the second dimension indexes the
+ parameters, and the third dimension represents the observations. }
+ \item{alpha}{The posterior draws of \eqn{\alpha}.}
+ \item{nstar}{The number of clusters at each Gibbs draw.}
+}
+
+\author{
+ Kosuke Imai, Department of Politics, Princeton University
+ \email{kimai at Princeton.Edu}, \url{http://www.princeton.edu/~kimai};
+ Ying Lu, Institute for Quantitative Social Sciences,
+ Harvard University \email{ylu at Latte.Harvard.Edu}}
+
+\references{
+ Imai, Kosuke and Ying Lu. (2004) \dQuote{ Parametric and Nonparametric
+ Bayesian Models for Ecological Inference in \eqn{2 \times 2}
+ Tables.} Proceedings of the American Statistical Association.
+\url{http://www.princeton.edu/~kimai/research/einonpar.html}
+
+Imai, Kosuke and Ying Lu. (2005) \dQuote{An Incomplete Data Approach
+ to Ecological Inference.} Working Paper, Princeton University,
+available at
+\url{http://www.princeton.edu/~kimai/research/einonpar.html}
+}
+
+\seealso{\code{eco}, \code{predict.eco}, \code{summary.ecoNP}}
+\keyword{models}
diff --git a/man/predict.eco.Rd b/man/predict.eco.Rd
new file mode 100644
index 0000000..ed8e53d
--- /dev/null
+++ b/man/predict.eco.Rd
@@ -0,0 +1,78 @@
+\name{predict.eco}
+
+\alias{predict.eco}
+\alias{predict.ecoX}
+
+\title{Out-of-Sample Posterior Prediction under the Parametric Bayesian
+Model for Ecological Inference in 2x2 Tables}
+
+\description{
+ Obtains out-of-sample posterior predictions under the fitted
+ parametric Bayesian model for ecological
+ inference. \code{predict} method for class \code{eco} and \code{ecoX}.
+}
+
+\usage{
+ \method{predict}{eco}(object, newdraw = NULL, subset = NULL, verbose = FALSE, ...)
+ \method{predict}{ecoX}(object, newdraw = NULL, subset = NULL, newdata = NULL, cond = FALSE, verbose = FALSE, ...)
+}
+
+\arguments{
+ \item{object}{An output object from \code{eco} or \code{ecoNP}.}
+ \item{newdraw}{An optional list containing two matrices (or three
+ dimensional arrays for the nonparametric model) of MCMC draws
+ of \eqn{\mu} and \eqn{\Sigma}. Those elements should be named as
+ \code{mu} and \code{Sigma}, respectively. The default is the
+ original MCMC draws stored in \code{object}.
+ }
+ \item{newdata}{An optional data frame containing a new data set for
+ which posterior predictions will be made. The new data set must have
+ the same variable names as those in the original data.
+ }
+ \item{subset}{A scalar or numerical vector specifying the row
+ number(s) of \code{mu} and \code{Sigma} in the output object from
+ \code{eco}. If specified, the posterior draws of parameters for
+ those rows are used for posterior prediction. The default is
+ \code{NULL} where all the posterior draws are used.
+ }
+ \item{cond}{logical. If \code{TRUE}, then the conditional prediction
+ will made for the parametric model with contextual effects. The
+ default is \code{FALSE}.
+ }
+ \item{verbose}{logical. If \code{TRUE}, helpful messages along with a
+ progress report on the Monte Carlo sampling from the posterior
+ predictive distributions are printed on the screen. The default is
+ \code{FALSE}.
+ }
+ \item{...}{further arguments passed to or from other methods.}
+}
+
+\details{The posterior predictive values are computed using the
+ Monte Carlo sample stored in the \code{eco} output (or other sample if
+ \code{newdraw} is specified). Given each Monte Carlo sample of the
+ parameters, we sample the vector-valued latent variable from the
+ appropriate multivariate Normal distribution. Then, we apply the
+ inverse logit transformation to obtain the predictive values of
+ proportions, \eqn{W}. The computation may be slow (especially for the
+ nonparametric model) if a large Monte Carlo sample of the model
+ parameters is used. In either case, setting \code{verbose = TRUE} may
+ be helpful in monitoring the progress of the code.
+}
+
+\value{
+ \code{predict.eco} yields a matrix of class \code{predict.eco}
+ containing the Monte Carlo sample from the posterior predictive
+ distribution of inner cells of ecological
+ tables. \code{summary.predict.eco} will summarize the output, and
+ \code{print.summary.predict.eco} will print the summary.
+}
+
+\seealso{\code{eco}, \code{predict.ecoNP}}
+
+\author{
+ Kosuke Imai, Department of Politics, Princeton University
+ \email{kimai at Princeton.Edu}; Ying Lu, Institute for Quantitative
+ Social Sciences, Harvard University \email{ylu at Latte.Harvard.Edu}
+}
+
+\keyword{methods}
diff --git a/man/predict.ecoNP.Rd b/man/predict.ecoNP.Rd
new file mode 100644
index 0000000..9e2e4dd
--- /dev/null
+++ b/man/predict.ecoNP.Rd
@@ -0,0 +1,81 @@
+\name{predict.eco}
+
+\alias{predict.ecoNP}
+\alias{predict.ecoNPX}
+
+\title{Out-of-Sample Posterior Prediction under the
+ Nonparametric Bayesian Model for Ecological Inference in 2x2 Tables}
+
+\description{
+ Obtains out-of-sample posterior predictions under the fitted
+ nonparametric Bayesian model for ecological
+ inference. \code{predict} method for class \code{ecoNP} and \code{ecoNPX}.
+}
+
+\usage{
+ \method{predict}{ecoNP}(object, newdraw = NULL, subset = NULL, obs = NULL, verbose = FALSE, ...)
+ \method{predict}{ecoNPX}(object, newdraw = NULL, subset = NULL, obs = NULL, cond = FALSE, verbose = FALSE, ...)
+}
+
+\arguments{
+ \item{object}{An output object from \code{ecoNP}.}
+ \item{newdraw}{An optional list containing two matrices (or three
+ dimensional arrays for the nonparametric model) of MCMC draws
+ of \eqn{\mu} and \eqn{\Sigma}. Those elements should be named as
+ \code{mu} and \code{Sigma}, respectively. The default is the
+ original MCMC draws stored in \code{object}.
+ }
+ \item{subset}{A scalar or numerical vector specifying the row
+ number(s) of \code{mu} and \code{Sigma} in the output object from
+ \code{eco}. If specified, the posterior draws of parameters for
+ those rows are used for posterior prediction. The default is
+ \code{NULL} where all the posterior draws are used.
+ }
+ \item{obs}{An integer or vector of integers specifying the observation
+ number(s) whose posterior draws will be used for predictions. The
+ default is \code{NULL} where all the observations in the data set
+ are selected.
+ }
+ \item{cond}{logical. If \code{TRUE}, then the conditional prediction
+ will made for the parametric model with contextual effects. The
+ default is \code{FALSE}.
+ }
+ \item{verbose}{logical. If \code{TRUE}, helpful messages along with a
+ progress report on the Monte Carlo sampling from the posterior
+ predictive distributions are printed on the screen. The default is
+ \code{FALSE}.
+ }
+ \item{...}{further arguments passed to or from other methods.}
+}
+
+\details{The posterior predictive values are computed using the
+ Monte Carlo sample stored in the \code{eco} or \code{ecoNP} output
+ (or other sample if
+ \code{newdraw} is specified). Given each Monte Carlo sample of the
+ parameters, we sample the vector-valued latent variable from the
+ appropriate multivariate Normal distribution. Then, we apply the
+ inverse logit transformation to obtain the predictive values of
+ proportions, \eqn{W}. The computation may be slow (especially for the
+ nonparametric model) if a large Monte Carlo sample of the model
+ parameters is used. In either case, setting \code{verbose = TRUE} may
+ be helpful in monitoring the progress of the code.
+}
+
+\value{
+ \code{predict.eco} yields a matrix of class \code{predict.eco}
+ containing the Monte Carlo sample from the posterior predictive
+ distribution of inner cells of ecological
+ tables. \code{summary.predict.eco} will summarize the output, and
+ \code{print.summary.predict.eco} will print the summary.
+}
+
+\seealso{\code{eco}, \code{ecoNP}, \code{summary.eco},
+ \code{summary.ecoNP}}
+
+\author{
+ Kosuke Imai, Department of Politics, Princeton University
+ \email{kimai at Princeton.Edu}; Ying Lu, Institute for Quantitative
+ Social Sciences, Harvard University \email{ylu at Latte.Harvard.Edu}
+}
+
+\keyword{methods}
diff --git a/man/reg.Rd b/man/reg.Rd
new file mode 100644
index 0000000..9a86276
--- /dev/null
+++ b/man/reg.Rd
@@ -0,0 +1,33 @@
+\name{reg}
+
+\docType{data}
+
+\alias{reg}
+
+\title{Voter Registration in US Southern States}
+
+\description{
+ This data set contains the racial composition, the registration rate,
+ the number of eligible voters as well as the actual observed racial
+ registration rates for every county in four US southern states:
+ Florida, Louisiana, North Carolina, and South Carolina.
+}
+
+\usage{data(reg)}
+
+\format{A data frame containing 5 variables and 275 observations
+ \tabular{lll}{
+ X \tab numeric \tab the fraction of Black voters \cr
+ Y \tab numeric \tab the fraction of voters who registered themselves\cr
+ N \tab numeric \tab the total number of voters in each county \cr
+ W1 \tab numeric \tab the actual fraction of Black voters who
+ registered themselves \cr
+ W2 \tab numeric \tab the actual fraction of White voters who registered themselves
+ }
+}
+\references{King, G. (1997). \dQuote{A Solution to the Ecological
+ Inference Problem: Reconstructing Individual Behavior from Aggregate
+ Data}. Princeton University Press, Princeton, NJ.
+}
+
+\keyword{datasets}
diff --git a/man/summary.eco.Rd b/man/summary.eco.Rd
new file mode 100644
index 0000000..60ddacf
--- /dev/null
+++ b/man/summary.eco.Rd
@@ -0,0 +1,72 @@
+\name{summary.eco}
+
+\alias{summary.eco}
+\alias{print.summary.eco}
+
+\title{Summarizing the Results for the Bayesian Parametric Model for
+ Ecological Inference in 2x2 Tables}
+
+\description{
+ \code{summary} method for class \code{eco}.
+}
+
+\usage{
+ \method{summary}{eco}(object, CI = c(2.5, 97.5), param = TRUE,
+ units = FALSE, subset = NULL, ...)
+
+ \method{print}{summary.eco}(x, digits = max(3, getOption("digits") - 3), ...)
+}
+
+\arguments{
+ \item{object}{An output object from \code{eco}.}
+ \item{CI}{A vector of lower and upper bounds for the Bayesian credible
+ intervals used to summarize the results. The default is the
+ equal tail 95 percent credible interval.
+ }
+ \item{x}{An object of class \code{summary.eco}.}
+ \item{digits}{the number of significant digits to use when printing.}
+ \item{param}{Logical. If \code{TRUE}, the posterior estimates of the
+ population parameters will be provided. The default value is
+ \code{TRUE}.
+ }
+ \item{units}{Logical. If \code{TRUE}, the in-sample predictions for
+ each unit or for a subset of units will be provided. The default
+ value is \code{FALSE}.
+ }
+ \item{subset}{A numeric vector indicating the subset of the units whose
+ in-sample predications to be provided when \code{units} is
+ \code{TRUE}. The default value is \code{NULL} where the in-sample
+ predictions for each unit will be provided.
+ }
+ \item{...}{further arguments passed to or from other methods.}
+}
+
+\value{
+ \code{summary.eco} yields an object of class \code{summary.eco}
+ containing the following elements:
+ \item{call}{The call from \code{eco}.}
+ \item{n.obs}{The number of units.}
+ \item{n.draws}{The number of Monte Carlo samples.}
+ \item{agg.table}{Aggregate posterior estimates of the marginal means
+ of \eqn{W_1} and \eqn{W_2} using \eqn{X} and \eqn{N} as weights.}
+ If \code{param = TRUE}, the following elements are also included:
+ \item{param.table}{Posterior estimates of model parameters: population
+ mean estimates of \eqn{W_1} and \eqn{W_2} and their logit
+ transformations.}
+ If \code{units = TRUE}, the following elements are also included:
+ \item{W1.table}{Unit-level posterior estimates for \eqn{W_1}.}
+ \item{W2.table}{Unit-level posterior estimates for \eqn{W_2}.}
+
+ This object can be printed by \code{print.summary.eco}
+}
+
+\seealso{\code{eco}, \code{predict.eco}}
+
+\author{
+ Kosuke Imai, Department of Politics, Princeton University
+ \email{kimai at Princeton.Edu}, \url{http://www.princeton.edu/~kimai};
+ Ying Lu, Institute for Quantitative Social Sciences,
+ Harvard University \email{ylu at Latte.Harvard.Edu}
+}
+
+\keyword{methods}
diff --git a/man/summary.ecoNP.Rd b/man/summary.ecoNP.Rd
new file mode 100644
index 0000000..8c5d759
--- /dev/null
+++ b/man/summary.ecoNP.Rd
@@ -0,0 +1,71 @@
+\name{summary.ecoNP}
+
+\alias{summary.ecoNP}
+\alias{print.summary.ecoNP}
+
+\title{Summarizing the Results for the Bayesian Nonparametric Model for
+Ecological Inference in 2x2 Tables }
+
+\description{
+ \code{summary} method for class \code{ecoNP}.
+}
+
+\usage{
+ \method{summary}{ecoNP}(object, CI = c(2.5, 97.5), param = FALSE,
+ units = FALSE, subset = NULL, ...)
+
+ \method{print}{summary.ecoNP}(x, digits = max(3, getOption("digits") - 3), ...)
+}
+
+\arguments{
+ \item{object}{An output object from \code{ecoNP}.}
+ \item{CI}{A vector of lower and upper bounds for the Bayesian credible
+ intervals used to summarize the results. The default is the equal
+ tail 95 percent credible interval.
+ }
+ \item{x}{An object of class \code{summary.ecoNP}.}
+ \item{digits}{the number of significant digits to use when printing.}
+ \item{param}{Logical. If \code{TRUE}, the posterior estimates of the
+ population parameters will be provided. The default value is
+ \code{FALSE}.
+ }
+ \item{units}{Logical. If \code{TRUE}, the in-sample predictions for
+ each unit or for a subset of units will be provided. The default
+ value is \code{FALSE}.
+ }
+ \item{subset}{A numeric vector indicating the subset of the units whose
+ in-sample predications to be provided when \code{units} is
+ \code{TRUE}. The default value is \code{NULL} where the in-sample
+ predictions for each unit will be provided.}
+ \item{...}{further arguments passed to or from other methods.}
+}
+
+\value{
+ \code{summary.ecoNP} yields an object of class \code{summary.ecoNP}
+ containing the following elements:
+ \item{call}{The call from \code{ecoNP}.}
+ \item{n.obs}{The number of units.}
+ \item{n.draws}{The number of Monte Carlo samples.}
+ \item{agg.table}{Aggregate posterior estimates of the marginal means
+ of \eqn{W_1} and \eqn{W_2} using \eqn{X} and \eqn{N} as weights.}
+ If \code{param = TRUE}, the following elements are also included:
+ \item{param.table}{Posterior estimates of model parameters: population
+ mean estimates of \eqn{W_1} and \eqn{W_2}. If \code{subset} is
+ specified, only a subset of the population parameters are included.}
+ If \code{unit = TRUE}, the following elements are also included:
+ \item{W1.table}{Unit-level posterior estimates for \eqn{W_1}.}
+ \item{W2.table}{Unit-level posterior estimates for \eqn{W_2}.}
+
+ This object can be printed by \code{print.summary.ecoNP}
+}
+
+\seealso{\code{ecoNP}, \code{predict.eco}}
+
+\author{
+ Kosuke Imai, Department of Politics, Princeton University
+ \email{kimai at Princeton.Edu}, \url{http://www.princeton.edu/~kimai};
+ Ying Lu, Institute for Quantitative Social Sciences,
+ Harvard University \email{ylu at Latte.Harvard.Edu}
+}
+
+\keyword{methods}
diff --git a/src/Makevars b/src/Makevars
new file mode 100644
index 0000000..c0c23fa
--- /dev/null
+++ b/src/Makevars
@@ -0,0 +1,2 @@
+ PKG_LIBS = $(LAPACK_LIBS) $(BLAS_LIBS) $(FLIBS)
+
diff --git a/src/bayes.c b/src/bayes.c
new file mode 100644
index 0000000..9093a71
--- /dev/null
+++ b/src/bayes.c
@@ -0,0 +1,73 @@
+/******************************************************************
+ This file is a part of eco: R Package for Fitting Bayesian Models
+ of Ecological Inference for 2x2 Tables
+ by Kosuke Imai and Ying Lu
+ Copyright: GPL version 2 or later.
+*******************************************************************/
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include <Rmath.h>
+#include <R_ext/Utils.h>
+#include <R.h>
+#include "vector.h"
+#include "subroutines.h"
+#include "rand.h"
+
+/** Normal-InvWishart updating
+ Y|mu, Sigma ~ N(mu, Sigma)
+ mu|Sigma ~ N(mu0, Sigma/tau0)
+ Sigma ~ InvWish(nu0, S0^{-1}) **/
+void NIWupdate(
+ double **Y, /* data */
+ double *mu, /* mean */
+ double **Sigma, /* variance */
+ double **InvSigma, /* precision */
+ double *mu0, /* prior mean */
+ double tau0, /* prior scale */
+ int nu0, /* prior df */
+ double **S0, /* prior scale */
+ int n_samp, /* sample size */
+ int n_dim) /* dimension */
+{
+ int i,j,k;
+ double *Ybar = doubleArray(n_dim);
+ double *mun = doubleArray(n_dim);
+ double **Sn = doubleMatrix(n_dim, n_dim);
+ double **mtemp = doubleMatrix(n_dim, n_dim);
+
+ for (j=0; j<n_dim; j++) {
+ Ybar[j] = 0;
+ for (i=0; i<n_samp; i++)
+ Ybar[j] += Y[i][j];
+ Ybar[j] /= n_samp;
+ for (k=0; k<n_dim; k++)
+ Sn[j][k] = S0[j][k];
+ }
+ for (j=0; j<n_dim; j++) {
+ mun[j] = (tau0*mu0[j]+n_samp*Ybar[j])/(tau0+n_samp);
+ for (k=0; k<n_dim; k++) {
+ Sn[j][k] += (tau0*n_samp)*(Ybar[j]-mu0[j])*(Ybar[k]-mu0[k])/(tau0+n_samp);
+ for (i=0; i<n_samp; i++)
+ Sn[j][k] += (Y[i][j]-Ybar[j])*(Y[i][k]-Ybar[k]);
+ /* conditioning on mu:
+ Sn[j][k]+=tau0*(mu[j]-mu0[j])*(mu[k]-mu0[k]);
+ Sn[j][k]+=(Y[i][j]-mu[j])*(Y[i][k]-mu[k]); */
+ }
+ }
+ dinv(Sn, n_dim, mtemp);
+ rWish(InvSigma, mtemp, nu0+n_samp, n_dim);
+ dinv(InvSigma, n_dim, Sigma);
+
+ for (j=0; j<n_dim; j++)
+ for (k=0; k<n_dim; k++)
+ mtemp[j][k] = Sigma[j][k]/(tau0+n_samp);
+
+ rMVN(mu, mun, mtemp, n_dim);
+
+ free(Ybar);
+ free(mun);
+ FreeMatrix(Sn, n_dim);
+ FreeMatrix(mtemp, n_dim);
+}
diff --git a/src/bayes.h b/src/bayes.h
new file mode 100644
index 0000000..8325f20
--- /dev/null
+++ b/src/bayes.h
@@ -0,0 +1,10 @@
+/******************************************************************
+ This file is a part of eco: R Package for Fitting Bayesian Models
+ of Ecological Inference for 2x2 Tables
+ by Kosuke Imai and Ying Lu
+ Copyright: GPL version 2 or later.
+*******************************************************************/
+
+void NIWupdate(double **Y, double *mu, double **Sigma, double **InvSigma,
+ double *mu0, double tau0, int nu0, double **S0,
+ int n_samp, int n_dim);
diff --git a/src/gibbsBase.c b/src/gibbsBase.c
new file mode 100644
index 0000000..b718510
--- /dev/null
+++ b/src/gibbsBase.c
@@ -0,0 +1,273 @@
+/******************************************************************
+ This file is a part of eco: R Package for Fitting Bayesian Models
+ of Ecological Inference for 2x2 Tables
+ by Kosuke Imai and Ying Lu
+ Copyright: GPL version 2 or later.
+*******************************************************************/
+
+#include <stddef.h>
+#include <stdio.h>
+#include <math.h>
+#include <Rmath.h>
+#include <R_ext/Utils.h>
+#include "vector.h"
+#include "subroutines.h"
+#include "rand.h"
+#include "bayes.h"
+#include "sample.h"
+
+/* Normal Parametric Model for 2x2 Tables */
+void cBaseeco(
+ /*data input */
+ double *pdX, /* data (X, Y) */
+ int *pin_samp, /* sample size */
+
+ /*MCMC draws */
+ int *n_gen, /* number of gibbs draws */
+ int *burn_in, /* number of draws to be burned in */
+ int *pinth, /* keep every nth draw */
+ int *verbose, /* 1 for output monitoring */
+
+ /* prior specification*/
+ int *pinu0, /* prior df parameter for InvWish */
+ double *pdtau0, /* prior scale parameter for Sigma */
+ double *mu0, /* prior mean for mu */
+ double *pdS0, /* prior scale for Sigma */
+ double *mustart, /* starting values for mu */
+ double *Sigmastart, /* starting values for Sigma */
+
+ /* incorporating survey data */
+ int *survey, /*1 if survey data available (set of W_1, W_2)
+ 0 not*/
+ int *sur_samp, /*sample size of survey data*/
+ double *sur_W, /*set of known W_1, W_2 */
+
+ /* incorporating homeogenous areas */
+ int *x1, /* 1 if X=1 type areas available
+ W_1 known, W_2 unknown */
+ int *sampx1, /* number X=1 type areas */
+ double *x1_W1, /* values of W_1 for X1 type areas */
+ int *x0, /* 1 if X=0 type areas available
+ W_2 known, W_1 unknown */
+ int *sampx0, /* number X=0 type areas */
+ double *x0_W2, /* values of W_2 for X0 type areas */
+
+ /* bounds of W1 */
+ double *minW1, double *maxW1,
+
+ /* flags */
+ int *parameter, /* 1 if save population parameter */
+ int *Grid, /* 1 if Grid algorithm is used; 0 for
+ Metropolis */
+
+ /* storage for Gibbs draws of mu/sigmat*/
+ double *pdSMu0, double *pdSMu1,
+ double *pdSSig00, double *pdSSig01, double *pdSSig11,
+
+ /* storage for Gibbs draws of W*/
+ double *pdSW1, double *pdSW2
+ ){
+
+ /* some integers */
+ int n_samp = *pin_samp; /* sample size */
+ int s_samp = *sur_samp; /* sample size of survey data */
+ int x1_samp = *sampx1; /* sample size for X=1 */
+ int x0_samp = *sampx0; /* sample size for X=0 */
+ int t_samp = n_samp+s_samp+x1_samp+x0_samp; /* total sample size */
+ int nth = *pinth;
+ int n_dim = 2; /* dimension */
+ int n_step = 1000; /* 1/The default size of grid step */
+
+ /* prior parameters */
+ double tau0 = *pdtau0; /* prior scale */
+ int nu0 = *pinu0; /* prior degrees of freedom */
+ double **S0 = doubleMatrix(n_dim, n_dim); /* The prior S parameter for InvWish */
+
+ /* data */
+ double **X = doubleMatrix(n_samp, n_dim); /* The Y and covariates */
+ double **W = doubleMatrix(t_samp, n_dim); /* The W1 and W2 matrix */
+ double **Wstar = doubleMatrix(t_samp, n_dim); /* logit tranformed W */
+ double **S_W = doubleMatrix(s_samp, n_dim); /* The known W1 and W2 matrix*/
+ double **S_Wstar = doubleMatrix(s_samp, n_dim); /* logit transformed S_W*/
+
+ /* grids */
+ double **W1g = doubleMatrix(n_samp, n_step); /* grids for W1 */
+ double **W2g = doubleMatrix(n_samp, n_step); /* grids for W2 */
+ int *n_grid = intArray(n_samp); /* grid size */
+
+ /* model parameters */
+ double *mu = doubleArray(n_dim); /* The mean */
+ double **Sigma = doubleMatrix(n_dim, n_dim); /* The covariance matrix */
+ double **InvSigma = doubleMatrix(n_dim, n_dim); /* The inverse covariance matrix */
+
+ /* misc variables */
+ int i, j, k, main_loop; /* used for various loops */
+ int itemp, itempS, itempC, itempA;
+ int progress = 1, itempP = ftrunc((double) *n_gen/10);
+ double dtemp, dtemp1;
+
+ /* get random seed */
+ GetRNGstate();
+
+ /* read the priors */
+ itemp=0;
+ for(k=0;k<n_dim;k++)
+ for(j=0;j<n_dim;j++) S0[j][k]=pdS0[itemp++];
+
+ /* read the data */
+ itemp = 0;
+ for (j = 0; j < n_dim; j++)
+ for (i = 0; i < n_samp; i++)
+ X[i][j] = pdX[itemp++];
+
+ /* Initialize W, Wstar for n_samp */
+ for (i=0; i< n_samp; i++) {
+ if (X[i][1]!=0 && X[i][1]!=1) {
+ W[i][0]=runif(minW1[i], maxW1[i]);
+ W[i][1]=(X[i][1]-X[i][0]*W[i][0])/(1-X[i][0]);
+ }
+ if (X[i][1]==0)
+ for (j=0; j<n_dim; j++) W[i][j]=0.0001;
+ if (X[i][1]==1)
+ for (j=0; j<n_dim; j++) W[i][j]=0.9999;
+ for (j=0; j<n_dim; j++)
+ Wstar[i][j]=log(W[i][j])-log(1-W[i][j]);
+ }
+
+ /* read homeogenous areas information */
+ if (*x1==1)
+ for (i=0; i<x1_samp; i++) {
+ W[(n_samp+i)][0]=x1_W1[i];
+ if (W[(n_samp+i)][0]==0) W[(n_samp+i)][0]=0.0001;
+ if (W[(n_samp+i)][0]==1) W[(n_samp+i)][0]=0.9999;
+ Wstar[(n_samp+i)][0]=log(W[(n_samp+i)][0])-log(1-W[(n_samp+i)][0]);
+ }
+ if (*x0==1)
+ for (i=0; i<x0_samp; i++) {
+ W[(n_samp+x1_samp+i)][1]=x0_W2[i];
+ if (W[(n_samp+x1_samp+i)][1]==0) W[(n_samp+x1_samp+i)][1]=0.0001;
+ if (W[(n_samp+x1_samp+i)][1]==1) W[(n_samp+x1_samp+i)][1]=0.9999;
+ Wstar[(n_samp+x1_samp+i)][1]=log(W[(n_samp+x1_samp+i)][1])-log(1-W[(n_samp+x1_samp+i)][1]);
+ }
+
+ /* read the survey data */
+ if (*survey==1) {
+ itemp = 0;
+ for (j=0; j<n_dim; j++)
+ for (i=0; i<s_samp; i++) {
+ S_W[i][j]=sur_W[itemp++];
+ if (S_W[i][j]==0) S_W[i][j]=0.0001;
+ if (S_W[i][j]==1) S_W[i][j]=0.9999;
+ S_Wstar[i][j]=log(S_W[i][j])-log(1-S_W[i][j]);
+ W[(n_samp+x1_samp+x0_samp+i)][j]=S_W[i][j];
+ Wstar[(n_samp+x1_samp+x0_samp+i)][j]=S_Wstar[i][j];
+ }
+ }
+
+ /* counters */
+ itempA=0; /* for alpha */
+ itempS=0; /* for storage */
+ itempC=0; /* control nth draw */
+
+ /*** calculate grids ***/
+ if (*Grid)
+ GridPrep(W1g, W2g, X, maxW1, minW1, n_grid, n_samp, n_step);
+
+ /* starting vales of mu and Sigma */
+ itemp = 0;
+ for(j=0;j<n_dim;j++){
+ mu[j] = mustart[j];
+ for(k=0;k<n_dim;k++)
+ Sigma[j][k]=Sigmastart[itemp++];
+ }
+ dinv(Sigma, n_dim, InvSigma);
+
+ /*** Gibbs sampler! ***/
+ if (*verbose)
+ Rprintf("Starting Gibbs Sampler...\n");
+ for(main_loop=0; main_loop<*n_gen; main_loop++){
+ /** update W, Wstar given mu, Sigma in regular areas **/
+ for (i=0;i<n_samp;i++){
+ if ( X[i][1]!=0 && X[i][1]!=1 ) {
+ if (*Grid)
+ rGrid(W[i], W1g[i], W2g[i], n_grid[i], mu, InvSigma, n_dim);
+ else
+ rMH(W[i], X[i], minW1[i], maxW1[i], mu, InvSigma, n_dim);
+ }
+ /*3 compute Wsta_i from W_i*/
+ Wstar[i][0]=log(W[i][0])-log(1-W[i][0]);
+ Wstar[i][1]=log(W[i][1])-log(1-W[i][1]);
+ }
+
+ /* update W2 given W1, mu and Sigma in x1 homeogeneous areas */
+ if (*x1==1)
+ for (i=0; i<x1_samp; i++) {
+ dtemp=mu[1]+Sigma[0][1]/Sigma[0][0]*(Wstar[n_samp+i][0]-mu[0]);
+ dtemp1=Sigma[1][1]*(1-Sigma[0][1]*Sigma[0][1]/(Sigma[0][0]*Sigma[1][1]));
+ dtemp1=sqrt(dtemp1);
+ Wstar[n_samp+i][1]=rnorm(dtemp, dtemp1);
+ W[n_samp+i][1]=exp(Wstar[n_samp+i][1])/(1+exp(Wstar[n_samp+i][1]));
+ }
+
+ /* update W1 given W2, mu and Sigma in x0 homeogeneous areas */
+ if (*x0==1)
+ for (i=0; i<x0_samp; i++) {
+ dtemp=mu[0]+Sigma[0][1]/Sigma[1][1]*(Wstar[n_samp+x1_samp+i][1]-mu[1]);
+ dtemp1=Sigma[0][0]*(1-Sigma[0][1]*Sigma[0][1]/(Sigma[0][0]*Sigma[1][1]));
+ dtemp1=sqrt(dtemp1);
+ Wstar[n_samp+x1_samp+i][0]=rnorm(dtemp, dtemp1);
+ W[n_samp+x1_samp+i][0]=exp(Wstar[n_samp+x1_samp+i][0])/(1+exp(Wstar[n_samp+x1_samp+i][0]));
+ }
+
+ /* update mu, Sigma given wstar using effective sample of Wstar */
+ NIWupdate(Wstar, mu, Sigma, InvSigma, mu0, tau0, nu0, S0, t_samp, n_dim);
+
+ /*store Gibbs draw after burn-in and every nth draws */
+ if (main_loop>=*burn_in){
+ itempC++;
+ if (itempC==nth){
+ pdSMu0[itempA]=mu[0];
+ pdSMu1[itempA]=mu[1];
+ pdSSig00[itempA]=Sigma[0][0];
+ pdSSig01[itempA]=Sigma[0][1];
+ pdSSig11[itempA]=Sigma[1][1];
+ itempA++;
+ for(i=0; i<(n_samp+x1_samp+x0_samp); i++){
+ pdSW1[itempS]=W[i][0];
+ pdSW2[itempS]=W[i][1];
+ itempS++;
+ }
+ itempC=0;
+ }
+ }
+ if (*verbose)
+ if (itempP == main_loop) {
+ Rprintf("%3d percent done.\n", progress*10);
+ itempP+=ftrunc((double) *n_gen/10); progress++;
+ R_FlushConsole();
+ }
+ R_CheckUserInterrupt();
+ } /* end of Gibbs sampler */
+
+ if(*verbose)
+ Rprintf("100 percent done.\n");
+
+ /** write out the random seed **/
+ PutRNGstate();
+
+ /* Freeing the memory */
+ FreeMatrix(X, n_samp);
+ FreeMatrix(W, t_samp);
+ FreeMatrix(Wstar, t_samp);
+ FreeMatrix(S_W, s_samp);
+ FreeMatrix(S_Wstar, s_samp);
+ FreeMatrix(S0, n_dim);
+ FreeMatrix(W1g, n_samp);
+ FreeMatrix(W2g, n_samp);
+ free(n_grid);
+ free(mu);
+ FreeMatrix(Sigma,n_dim);
+ FreeMatrix(InvSigma, n_dim);
+
+} /* main */
+
diff --git a/src/gibbsDP.c b/src/gibbsDP.c
new file mode 100644
index 0000000..eb75c88
--- /dev/null
+++ b/src/gibbsDP.c
@@ -0,0 +1,417 @@
+/******************************************************************
+ This file is a part of eco: R Package for Fitting Bayesian Models
+ of Ecological Inference for 2x2 Tables
+ by Kosuke Imai and Ying Lu
+ Copyright: GPL version 2 or later.
+*******************************************************************/
+
+#include <stddef.h>
+#include <stdio.h>
+#include <math.h>
+#include <Rmath.h>
+#include <R_ext/Utils.h>
+#include "vector.h"
+#include "subroutines.h"
+#include "rand.h"
+#include "bayes.h"
+#include "sample.h"
+
+void cDPeco(
+ /*data input */
+ double *pdX, /* data (X, Y) */
+ int *pin_samp, /* sample size */
+
+ /*MCMC draws */
+ int *n_gen, /* number of gibbs draws */
+ int *burn_in, /* number of draws to be burned in */
+ int *pinth, /* keep every nth draw */
+ int *verbose, /* 1 for output monitoring */
+
+ /* prior specification*/
+ int *pinu0, /* prior df parameter for InvWish */
+ double *pdtau0, /* prior scale parameter for Sigma under G0*/
+ double *mu0, /* prior mean for mu under G0 */
+ double *pdS0, /* prior scale for Sigma */
+
+ /* DP prior specification */
+ double *alpha0, /* precision parameter, can be fixed or updated*/
+ int *pinUpdate, /* 1 if alpha gets updated */
+ double *pda0, double *pdb0, /* prior for alpha if alpha updated*/
+
+ /*incorporating survey data */
+ int *survey, /* 1 if survey data available (set of W_1, W_2) */
+ /* 0 otherwise*/
+ int *sur_samp, /* sample size of survey data*/
+ double *sur_W, /* set of known W_1, W_2 */
+
+ /*incorporating homeogenous areas */
+ int *x1, /* 1 if X=1 type areas available
+ W_1 known, W_2 unknown */
+ int *sampx1, /* number X=1 type areas */
+ double *x1_W1, /* values of W_1 for X1 type areas */
+
+ int *x0, /* 1 if X=0 type areas available
+ W_2 known, W_1 unknown */
+ int *sampx0, /* number X=0 type areas */
+ double *x0_W2, /* values of W_2 for X0 type areas */
+
+ /* bounds of W1 */
+ double *minW1, double *maxW1,
+
+ /* storage */
+ int *parameter, /* 1 if save population parameter */
+ int *Grid, /* 1 if Grid algorithm used; \
+ 0 if Metropolis algorithm used*/
+
+ /* storage for Gibbs draws of mu/sigmat*/
+ double *pdSMu0, double *pdSMu1,
+ double *pdSSig00, double *pdSSig01, double *pdSSig11,
+ /* storage for Gibbs draws of W*/
+ double *pdSW1, double *pdSW2,
+ /* storage for Gibbs draws of alpha */
+ double *pdSa,
+ /* storage for nstar at each Gibbs draw*/
+ int *pdSn
+ ){
+ /*some integers */
+ int n_samp = *pin_samp; /* sample size */
+ int s_samp = *sur_samp; /* sample size of survey data */
+ int x1_samp = *sampx1; /* sample size for X=1 */
+ int x0_samp = *sampx0; /* sample size for X=0 */
+ int t_samp = n_samp+x1_samp+x0_samp+s_samp; /* total sample size */
+ int nth = *pinth; /* keep every nth draw */
+ int n_dim = 2; /* dimension */
+ int n_step=1000; /* The default size of grid step */
+
+ /*prior parameters */
+ double tau0 = *pdtau0; /* prior scale */
+ int nu0 = *pinu0; /* prior degree of freedom*/
+ double **S0 = doubleMatrix(n_dim,n_dim);/*The prior S parameter for InvWish*/
+ double alpha = *alpha0; /* precision parameter*/
+ double a0 = *pda0, b0 = *pdb0; /* hyperprior for alpha */
+
+ /* data */
+ double **X = doubleMatrix(n_samp,n_dim); /* The Y and covariates */
+ double **W = doubleMatrix(t_samp,n_dim); /* The W1 and W2 matrix */
+ double **Wstar = doubleMatrix(t_samp,n_dim); /* The pseudo data */
+ double **S_W = doubleMatrix(s_samp,n_dim); /* The known W1 and W2 matrix*/
+ double **S_Wstar = doubleMatrix(s_samp,n_dim); /* The logit transformed S_W*/
+
+ /* grids */
+ double **W1g = doubleMatrix(n_samp, n_step); /* grids for W1 */
+ double **W2g = doubleMatrix(n_samp, n_step); /* grids for W2 */
+ int *n_grid = intArray(n_samp); /* grids size */
+
+ /* Model parameters */
+ /* Dirichlet variables */
+ double **mu = doubleMatrix(t_samp,n_dim); /* mean matrix */
+ double ***Sigma = doubleMatrix3D(t_samp,n_dim,n_dim); /*covarince matrix*/
+ double ***InvSigma = doubleMatrix3D(t_samp,n_dim,n_dim); /* inv of Sigma*/
+
+ int nstar; /* # clusters with distict theta values */
+ int *C = intArray(t_samp); /* vector of cluster membership */
+ double *q = doubleArray(t_samp); /* Weights of posterior of Dirichlet */
+ double *qq = doubleArray(t_samp); /* cumulative weight vector of q */
+ double **S_bvt = doubleMatrix(n_dim,n_dim); /* S paramter for BVT in q0 */
+
+ /* variables defined in remixing step: cycle through all clusters */
+ double **Wstarmix = doubleMatrix(t_samp,n_dim); /*data matrix used */
+ double *mu_mix = doubleArray(n_dim); /*updated MEAN parameter */
+ double **Sigma_mix = doubleMatrix(n_dim,n_dim); /*updated VAR parameter */
+ double **InvSigma_mix = doubleMatrix(n_dim,n_dim); /* Inv of Sigma_mix */
+ int nj; /* record # of obs in each cluster */
+ int *sortC = intArray(t_samp); /* record (sorted)original obs id */
+ int *indexC = intArray(t_samp); /* record original obs id */
+ int *label = intArray(t_samp); /* store index values */
+
+ /* misc variables */
+ int i, j, k, l, main_loop; /* used for various loops */
+ int itemp;
+ int itempA=0; /* counter for alpha */
+ int itempS=0; /* counter for storage */
+ int itempC=0; /* counter to control nth draw */
+ int progress = 1, itempP = ftrunc((double) *n_gen/10);
+ double dtemp, dtemp1;
+ double *vtemp = doubleArray(n_dim);
+ double **mtemp = doubleMatrix(n_dim,n_dim);
+ double **mtemp1 = doubleMatrix(n_dim,n_dim);
+ double **onedata = doubleMatrix(1, n_dim);
+
+ /* get random seed */
+ GetRNGstate();
+
+ /* read priors under G0*/
+ itemp=0;
+ for(k=0;k<n_dim;k++)
+ for(j=0;j<n_dim;j++) S0[j][k]=pdS0[itemp++];
+
+ /* read the data set */
+ itemp = 0;
+ for (j = 0; j < n_dim; j++)
+ for (i = 0; i < n_samp; i++) X[i][j] = pdX[itemp++];
+
+ /*Intialize W, Wsatr for n_samp */
+ for (i=0; i< n_samp; i++) {
+ if (X[i][1]!=0 && X[i][1]!=1) {
+ W[i][0]=runif(minW1[i], maxW1[i]);
+ W[i][1]=(X[i][1]-X[i][0]*W[i][0])/(1-X[i][0]);
+ }
+ if (X[i][1]==0)
+ for (j=0; j<n_dim; j++) W[i][j]=0.0001;
+ if (X[i][1]==1)
+ for (j=0; j<n_dim; j++) W[i][j]=0.9999;
+ for (j=0; j<n_dim; j++)
+ Wstar[i][j]=log(W[i][j])-log(1-W[i][j]);
+ }
+
+ /*read homeogenous areas information */
+ if (*x1==1)
+ for (i=0; i<x1_samp; i++) {
+ W[(n_samp+i)][0]=x1_W1[i];
+ if (W[(n_samp+i)][0]==0) W[(n_samp+i)][0]=0.0001;
+ if (W[(n_samp+i)][0]==1) W[(n_samp+i)][0]=0.9999;
+ Wstar[(n_samp+i)][0]=log(W[(n_samp+i)][0])-log(1-W[(n_samp+i)][0]);
+ }
+
+ if (*x0==1)
+ for (i=0; i<x0_samp; i++) {
+ W[(n_samp+x1_samp+i)][1]=x0_W2[i];
+ if (W[(n_samp+x1_samp+i)][1]==0) W[(n_samp+x1_samp+i)][1]=0.0001;
+ if (W[(n_samp+x1_samp+i)][1]==1) W[(n_samp+x1_samp+i)][1]=0.9999;
+ Wstar[(n_samp+x1_samp+i)][1]=log(W[(n_samp+x1_samp+i)][1])-log(1-W[(n_samp+x1_samp+i)][1]);
+ }
+
+ /*read the survey data */
+ if (*survey==1) {
+ itemp = 0;
+ for (j=0; j<n_dim; j++)
+ for (i=0; i<s_samp; i++) {
+ S_W[i][j]=sur_W[itemp++];
+ if (S_W[i][j]==0) S_W[i][j]=0.0001;
+ if (S_W[i][j]==1) S_W[i][j]=0.9999;
+ S_Wstar[i][j]=log(S_W[i][j])-log(1-S_W[i][j]);
+ W[n_samp+x1_samp+x0_samp+i][j]=S_W[i][j];
+ Wstar[n_samp+x1_samp+x0_samp+i][j]=S_Wstar[i][j];
+ }
+ }
+
+ /* Calcualte grids */
+ if (*Grid)
+ GridPrep(W1g,W2g, X, maxW1, minW1, n_grid, n_samp, n_step);
+
+ /* parmeters for Bivaraite t-distribution-unchanged in MCMC */
+ for (j=0;j<n_dim;j++)
+ for(k=0;k<n_dim;k++)
+ mtemp[j][k]=S0[j][k]*(1+tau0)/(tau0*(nu0-n_dim+1));
+ dinv(mtemp, n_dim, S_bvt);
+
+ /**draw initial values of mu_i, Sigma_i under G0 for all effective sample**/
+ /*1. Sigma_i under InvWish(nu0, S0^-1) with E(Sigma)=S0/(nu0-3)*/
+ /* InvSigma_i under Wish(nu0, S0^-1 */
+ /*2. mu_i|Sigma_i under N(mu0, Sigma_i/tau0) */
+ dinv(S0, n_dim, mtemp);
+ for(i=0;i<t_samp;i++){
+ /*draw from wish(nu0, S0^-1) */
+ rWish(InvSigma[i], mtemp, nu0, n_dim);
+ dinv(InvSigma[i], n_dim, Sigma[i]);
+ for (j=0;j<n_dim;j++)
+ for(k=0;k<n_dim;k++) mtemp1[j][k]=Sigma[i][j][k]/tau0;
+ rMVN(mu[i], mu0, mtemp1, n_dim);
+ }
+
+ /* initialize the cluster membership */
+ nstar=t_samp; /* the # of disticnt values */
+ for(i=0;i<t_samp;i++)
+ C[i]=i; /*cluster is from 0...n_samp-1 */
+
+ if (*verbose)
+ Rprintf("Starting Gibbs Sampler...\n");
+ for(main_loop=0; main_loop<*n_gen; main_loop++){
+ /**update W, Wstar given mu, Sigma only for the unknown W/Wstar**/
+ for (i=0;i<n_samp;i++){
+ if (X[i][1]!=0 && X[i][1]!=1) {
+ if (*Grid)
+ rGrid(W[i], W1g[i], W2g[i], n_grid[i], mu[i], InvSigma[i], n_dim);
+ else
+ rMH(W[i], X[i], minW1[i], maxW1[i], mu[i], InvSigma[i], n_dim);
+ }
+ /*3 compute Wsta_i from W_i*/
+ Wstar[i][0]=log(W[i][0])-log(1-W[i][0]);
+ Wstar[i][1]=log(W[i][1])-log(1-W[i][1]);
+ }
+
+ if (*x1==1)
+ for (i=0; i<x1_samp; i++) {
+ dtemp=mu[n_samp+i][1]+Sigma[n_samp+i][0][1]/Sigma[n_samp+i][0][0]*(Wstar[n_samp+i][0]-mu[n_samp+i][0]);
+ dtemp1=Sigma[n_samp+i][1][1]*(1-Sigma[n_samp+i][0][1]*Sigma[n_samp+i][0][1]/(Sigma[n_samp+i][0][0]*Sigma[n_samp+i][1][1]));
+ Wstar[n_samp+i][1]=norm_rand()*sqrt(dtemp1)+dtemp;
+ W[n_samp+i][1]=exp(Wstar[n_samp+i][1])/(1+exp(Wstar[n_samp+i][1]));
+ }
+
+ /*update W1 given W2, mu_ord and Sigma_ord in x0 homeogeneous areas */
+ if (*x0==1)
+ for (i=0; i<x0_samp; i++) {
+ dtemp=mu[n_samp+x1_samp+i][0]+Sigma[n_samp+x1_samp+i][0][1]/Sigma[n_samp+x1_samp+i][1][1]*(Wstar[n_samp+x1_samp+i][1]-mu[n_samp+x1_samp+i][1]);
+ dtemp1=Sigma[n_samp+x1_samp+i][0][0]*(1-Sigma[n_samp+x1_samp+i][0][1]*Sigma[n_samp+x1_samp+i][0][1]/(Sigma[n_samp+x1_samp+i][0][0]*Sigma[n_samp+x1_samp+i][1][1]));
+ Wstar[n_samp+i][0]=norm_rand()*sqrt(dtemp1)+dtemp;
+ W[n_samp+x1_samp+i][0]=exp(Wstar[n_samp+x1_samp+i][0])/(1+exp(Wstar[n_samp+x1_samp+i][0]));
+ }
+
+ /**updating mu, Sigma given Wstar uisng effective sample size W_star**/
+ for (i=0; i<t_samp; i++){
+ /* generate weight vector q */
+ dtemp=0;
+ for (j=0; j<t_samp; j++){
+ if (j!=i)
+ q[j]=dMVN(Wstar[i], mu[j], InvSigma[j], n_dim, 0);
+ else
+ q[j]=alpha*dMVT(Wstar[i], mu0, S_bvt, nu0-n_dim+1, 2, 0);
+ dtemp+=q[j];
+ qq[j]=dtemp; /*compute qq, the cumlative of q*/
+ }
+ /*standardize q and qq */
+ for (j=0; j<t_samp; j++) qq[j]/=dtemp;
+
+ /** draw the configuration parameter **/
+ j=0; dtemp=unif_rand();
+ while (dtemp > qq[j]) j++;
+ /** Dirichlet update Sigma_i, mu_i|Sigma_i **/
+ /* j=i: posterior update given Wstar[i] */
+ if (j==i){
+ onedata[0][0] = Wstar[i][0];
+ onedata[0][1] = Wstar[i][1];
+ NIWupdate(onedata, mu[i], Sigma[i], InvSigma[i], mu0, tau0,nu0, S0, 1, n_dim);
+ C[i]=nstar;
+ nstar++;
+ }
+ /* j=i': replace with i' obs */
+ else {
+ /*1. mu_i=mu_j, Sigma_i=Sigma_j*/
+ /*2. update C[i]=C[j] */
+ for(k=0;k<n_dim;k++) {
+ mu[i][k]=mu[j][k];
+ for(l=0;l<n_dim;l++) {
+ Sigma[i][k][l]=Sigma[j][k][l];
+ InvSigma[i][k][l]=InvSigma[j][k][l];
+ }
+ }
+ C[i]=C[j];
+ }
+ sortC[i]=C[i];
+ } /* end of i loop*/
+
+ /** remixing step using effective sample of Wstar**/
+ for(i=0;i<t_samp;i++) indexC[i]=i;
+ R_qsort_int_I(sortC, indexC, 1, t_samp);
+
+ nstar=0; i=0;
+ while (i<t_samp){
+ j=sortC[i]; /*saves the first element in a block of same values */
+ nj=0; /* counter for a block of same values */
+
+ /* get data for remixing */
+ while ((sortC[i]==j) && (i<t_samp)) {
+ label[nj]=indexC[i];
+ for (k=0; k<n_dim; k++)
+ Wstarmix[nj][k]=Wstar[label[nj]][k];
+ nj++; i++;
+ } /* i records the current position in IndexC */
+ /* nj records the # of obs in Psimix */
+
+ /** posterior update for mu_mix, Sigma_mix based on Psimix **/
+ NIWupdate(Wstarmix, mu_mix,Sigma_mix, InvSigma_mix, mu0, tau0, nu0, S0, nj, n_dim);
+
+ /**update mu, Simgat with mu_mix, Sigmat_mix via label**/
+ for (j=0;j<nj;j++){
+ C[label[j]]=nstar; /* updating C vector with no gap */
+ for (k=0; k<n_dim; k++){
+ mu[label[j]][k]=mu_mix[k];
+ for (l=0;l<n_dim;l++){
+ Sigma[label[j]][k][l]=Sigma_mix[k][l];
+ InvSigma[label[j]][k][l]=InvSigma_mix[k][l];
+ }
+ }
+ }
+ nstar++; /*finish update one distinct value*/
+ } /* nstar is the number of distinct values */
+
+ /** updating alpha **/
+ if(*pinUpdate) {
+ dtemp=b0-log(rbeta(alpha+1, (double) t_samp));
+ dtemp1=(double)(a0+nstar-1)/(t_samp*dtemp);
+ if(unif_rand() < dtemp1)
+ alpha=rgamma(a0+nstar, 1/dtemp);
+ else
+ alpha=rgamma(a0+nstar-1, 1/dtemp);
+ }
+
+ /*store Gibbs draws after burn_in */
+ R_CheckUserInterrupt();
+ if (main_loop>=*burn_in) {
+ itempC++;
+ if (itempC==nth){
+ if(*pinUpdate) {
+ pdSa[itempA]=alpha;
+ }
+ pdSn[itempA]=nstar;
+ itempA++;
+
+ for(i=0; i<(n_samp+x1_samp+x0_samp); i++) {
+ pdSMu0[itempS]=mu[i][0];
+ pdSMu1[itempS]=mu[i][1];
+ pdSSig00[itempS]=Sigma[i][0][0];
+ pdSSig01[itempS]=Sigma[i][0][1];
+ pdSSig11[itempS]=Sigma[i][1][1];
+ pdSW1[itempS]=W[i][0];
+ pdSW2[itempS]=W[i][1];
+ itempS++;
+ }
+ itempC=0;
+ }
+ }
+
+ if (*verbose)
+ if (itempP == main_loop) {
+ Rprintf("%3d percent done.\n", progress*10);
+ itempP+=ftrunc((double) *n_gen/10); progress++;
+ R_FlushConsole();
+ }
+ } /*end of MCMC for DP*/
+
+ if (*verbose)
+ Rprintf("100 percent done.\n");
+
+ /** write out the random seed **/
+ PutRNGstate();
+
+ /* Freeing the memory */
+ FreeMatrix(S0, n_dim);
+ FreeMatrix(X, n_samp);
+ FreeMatrix(W, t_samp);
+ FreeMatrix(Wstar, t_samp);
+ FreeMatrix(S_W, s_samp);
+ FreeMatrix(S_Wstar, s_samp);
+ FreeMatrix(W1g, n_samp);
+ FreeMatrix(W2g, n_samp);
+ free(n_grid);
+ FreeMatrix(mu, t_samp);
+ Free3DMatrix(Sigma, t_samp,n_dim);
+ Free3DMatrix(InvSigma, t_samp, n_dim);
+ free(C);
+ free(q);
+ free(qq);
+ FreeMatrix(S_bvt, n_dim);
+ FreeMatrix(Wstarmix, t_samp);
+ free(mu_mix);
+ FreeMatrix(Sigma_mix, n_dim);
+ FreeMatrix(InvSigma_mix, n_dim);
+ free(sortC);
+ free(indexC);
+ free(label);
+ free(vtemp);
+ FreeMatrix(mtemp, n_dim);
+ FreeMatrix(mtemp1, n_dim);
+ free(onedata);
+} /* main */
+
+
diff --git a/src/gibbsXBase.c b/src/gibbsXBase.c
new file mode 100644
index 0000000..6acbf47
--- /dev/null
+++ b/src/gibbsXBase.c
@@ -0,0 +1,305 @@
+/******************************************************************
+ This file is a part of eco: R Package for Fitting Bayesian Models
+ of Ecological Inference for 2x2 Tables
+ by Kosuke Imai and Ying Lu
+ Copyright: GPL version 2 or later.
+*******************************************************************/
+
+#include <stddef.h>
+#include <stdio.h>
+#include <math.h>
+#include <Rmath.h>
+#include <R_ext/Utils.h>
+#include "vector.h"
+#include "subroutines.h"
+#include "rand.h"
+#include "bayes.h"
+#include "sample.h"
+
+/* Normal Parametric Model for 2x2 Tables with Contextual Effects */
+void cBaseecoX(
+ /*data input */
+ double *pdX, /* data (X, Y) */
+ int *pin_samp, /* sample size */
+
+ /*MCMC draws */
+ int *n_gen, /* number of gibbs draws */
+ int *burn_in, /* number of draws to be burned in */
+ int *pinth, /* keep every nth draw */
+ int *verbose, /* 1 for output monitoring */
+
+ /* prior specification*/
+ int *pinu0, /* prior df parameter for InvWish */
+ double *pdtau0, /* prior scale parameter for Sigma under G0*/
+ double *mu0, /* prior mean for mu under G0 */
+ double *pdS0, /* prior scale for Sigma */
+ double *mustart, /* starting values for mu */
+ double *Sigmastart, /* starting values for Sigma */
+
+ /*incorporating survey data */
+ int *survey, /*1 if survey data available (set of W_1, W_2)
+ 0 not*/
+ int *sur_samp, /*sample size of survey data*/
+ double *sur_W, /*set of known W_1, W_2 */
+
+ /* incorporating homeogenous areas */
+ int *x1, /* 1 if X=1 type areas available W_1 known, W_2 unknown */
+ int *sampx1, /* number X=1 type areas */
+ double *x1_W1, /* values of W_1 for X1 type areas */
+
+ int *x0, /* 1 if X=0 type areas available W_2
+ known, W_1 unknown */
+ int *sampx0, /* number X=0 type areas */
+ double *x0_W2, /* values of W_2 for X0 type areas */
+
+ /* bounds fo W1 */
+ double *minW1, double *maxW1,
+
+ /* flags */
+ int *parameter, /* 1 if save population parameter */
+ int *Grid, /* 1 if Grid algorithm is used; 0 for
+ Metropolis */
+
+ /* storage for Gibbs draws of mu/sigmat*/
+ double *pdSMu0, double *pdSMu1, double *pdSMu2,
+ double *pdSSig00, double *pdSSig01, double *pdSSig02,
+ double *pdSSig11, double *pdSSig12, double *pdSSig22,
+
+ /* storage for Gibbs draws of W*/
+ double *pdSW1, double *pdSW2
+ ){
+
+ /* some integers */
+ int n_samp = *pin_samp; /* sample size */
+ int s_samp = *sur_samp; /* sample size of survey data */
+ int x1_samp = *sampx1; /* sample size for X=1 */
+ int x0_samp = *sampx0; /* sample size for X=0 */
+ int t_samp = n_samp+s_samp+x1_samp+x0_samp; /* total sample size */
+ int nth = *pinth;
+ int n_dim = 2; /* dimension */
+ int n_step = 1000; /* 1/The default size of grid step */
+
+ /* prior parameters */
+ double tau0 = *pdtau0;
+ int nu0 = *pinu0;
+ double **S0 = doubleMatrix(n_dim+1,n_dim+1); /* The prior S parameter for InvWish */
+
+ /* data */
+ double **X = doubleMatrix(n_samp,n_dim); /* The Y and covariates */
+ double **W = doubleMatrix(t_samp,n_dim); /* The W1 and W2 matrix */
+ double **Wstar = doubleMatrix(t_samp,n_dim+1); /* logit transformed
+ W and X */
+ double **S_W = doubleMatrix(s_samp, n_dim+1); /* known W1, W2, X */
+ double **S_Wstar = doubleMatrix(s_samp, n_dim+1); /* logit
+ transformed S_W */
+ /* grids */
+ double **W1g = doubleMatrix(n_samp, n_step);
+ double **W2g = doubleMatrix(n_samp, n_step);
+ int *n_grid = intArray(n_samp); /* grid size */
+
+ /* ordinary model variables */
+ double *mu = doubleArray(n_dim+1);
+ double **Sigma = doubleMatrix(n_dim+1,n_dim+1);
+ double **InvSigma = doubleMatrix(n_dim+1,n_dim+1);
+
+ /* conditional mean & variance for (W1, W2) given X */
+ double *mu_w = doubleArray(n_dim);
+ double **Sigma_w = doubleMatrix(n_dim,n_dim);
+ double **InvSigma_w = doubleMatrix(n_dim,n_dim);
+
+ /* misc variables */
+ int i, j, k, t, main_loop; /* used for various loops */
+ int itemp, itempS, itempC, itempA;
+ int progress = 1, itempP = ftrunc((double) *n_gen/10);
+ double dtemp, dtemp1;
+
+ /* get random seed */
+ GetRNGstate();
+
+ /* priors */
+ itemp = 0;
+ for(k=0; k<(n_dim+1); k++)
+ for(j=0; j<(n_dim+1); j++)
+ S0[j][k] = pdS0[itemp++];
+
+ /* read the data set */
+ itemp = 0;
+ for (j = 0; j < n_dim; j++)
+ for (i = 0; i < n_samp; i++)
+ X[i][j] = pdX[itemp++];
+
+ /* Initialize W, Wstar for n_samp */
+ for (i=0; i< n_samp; i++) {
+ if (X[i][1]!=0 && X[i][1]!=1) {
+ W[i][0]=runif(minW1[i], maxW1[i]);
+ W[i][1]=(X[i][1]-X[i][0]*W[i][0])/(1-X[i][0]);
+ }
+ if (X[i][1]==0)
+ for (j=0; j<n_dim; j++) W[i][j]=0.0001;
+ if (X[i][1]==1)
+ for (j=0; j<n_dim; j++) W[i][j]=0.9999;
+ for (j=0; j<n_dim; j++)
+ Wstar[i][j]=log(W[i][j])-log(1-W[i][j]);
+ Wstar[i][n_dim] = log(X[i][0])-log(1-X[i][0]);
+ }
+
+ /*read homeogenous areas information */
+ if (*x1==1)
+ for (i=0; i<x1_samp; i++) {
+ W[(n_samp+i)][0]=x1_W1[i];
+ if (W[(n_samp+i)][0]==0) W[(n_samp+i)][0]=0.0001;
+ if (W[(n_samp+i)][0]==1) W[(n_samp+i)][0]=0.9999;
+ Wstar[(n_samp+i)][0]=log(W[(n_samp+i)][0])-log(1-W[(n_samp+i)][0]);
+ Wstar[n_samp+i][n_dim]=log(0.9999)-log(0.0001);
+ }
+
+ if (*x0==1)
+ for (i=0; i<x0_samp; i++) {
+ W[(n_samp+x1_samp+i)][1]=x0_W2[i];
+ if (W[(n_samp+x1_samp+i)][1]==0) W[(n_samp+x1_samp+i)][1]=0.0001;
+ if (W[(n_samp+x1_samp+i)][1]==1) W[(n_samp+x1_samp+i)][1]=0.9999;
+ Wstar[(n_samp+x1_samp+i)][1]=log(W[(n_samp+x1_samp+i)][1])-log(1-W[(n_samp+x1_samp+i)][1]);
+ Wstar[n_samp+x1_samp+i][n_dim]=log(0.0001)-log(0.9999);
+ }
+
+ /*read the survey data */
+ if (*survey==1) {
+ itemp = 0;
+ for (j=0; j<=n_dim; j++)
+ for (i=0; i<s_samp; i++) {
+ S_W[i][j]=sur_W[itemp++];
+ if (S_W[i][j]==0) S_W[i][j]=0.0001;
+ if (S_W[i][j]==1) S_W[i][j]=0.9999;
+ S_Wstar[i][j]=log(S_W[i][j])-log(1-S_W[i][j]);
+ if (j<n_dim) {
+ W[(n_samp+x1_samp+x0_samp+i)][j]=S_W[i][j];
+ Wstar[(n_samp+x1_samp+x0_samp+i)][j]=S_Wstar[i][j];
+ }
+ else
+ Wstar[(n_samp+x1_samp+x0_samp+i)][j]=S_Wstar[i][j];
+ }
+ }
+
+ /* counters */
+ itempA=0; /* for alpha */
+ itempS=0; /* for storage */
+ itempC=0; /* control nth draw */
+
+ /*** calculate grids ***/
+ if (*Grid)
+ GridPrep(W1g, W2g, X, maxW1, minW1, n_grid, n_samp, n_step);
+
+ /* starting values of mu and Sigma */
+ itemp = 0;
+ for(j=0;j<(n_dim+1);j++){
+ mu[j] = mustart[j];
+ for(k=0;k<(n_dim+1);k++)
+ Sigma[j][k]=Sigmastart[itemp++];
+ }
+ dinv(Sigma, n_dim+1, InvSigma);
+
+ /***Gibbs Sampler ***/
+ if (*verbose)
+ Rprintf("Starting Gibbs Sampler...\n");
+ for(main_loop=0; main_loop<*n_gen; main_loop++){
+ /* conditional variance */
+ for (j=0; j<n_dim; j++)
+ for (k=0; k<n_dim; k++)
+ Sigma_w[j][k]=Sigma[j][k]-Sigma[n_dim][j]/Sigma[n_dim][n_dim]*Sigma[n_dim][k];
+ dinv(Sigma_w, n_dim, InvSigma_w);
+
+ /**update W, Wstar given mu, Sigma in regular areas**/
+ for (i=0; i<n_samp; i++){
+ for (j=0; j<n_dim; j++)
+ mu_w[j]=mu[j]+Sigma[n_dim][j]/Sigma[n_dim][n_dim]*(Wstar[i][2]-mu[n_dim]);
+ if ( X[i][1]!=0 && X[i][1]!=1 ) {
+ if (*Grid)
+ rGrid(W[i], W1g[i],W2g[i], n_grid[i], mu_w, InvSigma_w,
+ n_dim);
+ else
+ rMH(W[i], X[i], minW1[i], maxW1[i], mu_w, InvSigma_w, n_dim);
+ }
+ /*3 compute Wsta_i from W_i*/
+ Wstar[i][0]=log(W[i][0])-log(1-W[i][0]);
+ Wstar[i][1]=log(W[i][1])-log(1-W[i][1]);
+ }
+
+ /*update W2 given W1, mu and Sigma in x1 homeogeneous areas */
+ if (*x1==1)
+ for (i=0; i<x1_samp; i++) {
+ dtemp=mu_w[1]+Sigma_w[0][1]/Sigma_w[0][0]*(Wstar[n_samp+i][0]-mu_w[0]);
+ dtemp1=Sigma_w[1][1]*(1-Sigma_w[0][1]*Sigma_w[0][1]/(Sigma_w[0][0]*Sigma_w[1][1]));
+ dtemp1=sqrt(dtemp1);
+ Wstar[n_samp+i][1]=rnorm(dtemp, dtemp1);
+ W[n_samp+i][1]=exp(Wstar[n_samp+i][1])/(1+exp(Wstar[n_samp+i][1]));
+ }
+
+ /*update W1 given W2, mu and Sigma in x0 homeogeneous areas */
+ if (*x0==1)
+ for (i=0; i<x0_samp; i++) {
+ dtemp=mu_w[0]+Sigma_w[0][1]/Sigma_w[1][1]*(Wstar[n_samp+x1_samp+i][1]-mu_w[1]);
+ dtemp1=Sigma_w[0][0]*(1-Sigma_w[0][1]*Sigma_w[0][1]/(Sigma_w[0][0]*Sigma_w[1][1]));
+ dtemp1=sqrt(dtemp1);
+ Wstar[n_samp+x1_samp+i][0]=rnorm(dtemp, dtemp1);
+ W[n_samp+x1_samp+i][0]=exp(Wstar[n_samp+x1_samp+i][0])/(1+exp(Wstar[n_samp+x1_samp+i][0]));
+ }
+
+ /* update mu, Sigma given wstar using effective sample of Wstar */
+ NIWupdate(Wstar, mu, Sigma, InvSigma, mu0, tau0, nu0, S0, t_samp, n_dim+1);
+
+ /*store Gibbs draw after burn-in and every nth draws */
+ R_CheckUserInterrupt();
+ if (main_loop>=*burn_in){
+ itempC++;
+ if (itempC==nth){
+ pdSMu0[itempA]=mu[0];
+ pdSMu1[itempA]=mu[1];
+ pdSMu2[itempA]=mu[2];
+ pdSSig00[itempA]=Sigma[0][0];
+ pdSSig01[itempA]=Sigma[0][1];
+ pdSSig02[itempA]=Sigma[0][2];
+ pdSSig11[itempA]=Sigma[1][1];
+ pdSSig12[itempA]=Sigma[1][2];
+ pdSSig22[itempA]=Sigma[2][2];
+ itempA++;
+ for(i=0; i<(n_samp+x1_samp+x0_samp); i++){
+ pdSW1[itempS]=W[i][0];
+ pdSW2[itempS]=W[i][1];
+ itempS++;
+ }
+ itempC=0;
+ }
+ } /*end of stroage *burn_in*/
+ if (*verbose)
+ if (itempP == main_loop) {
+ Rprintf("%3d percent done.\n", progress*10);
+ itempP+=ftrunc((double) *n_gen/10); progress++;
+ R_FlushConsole();
+ }
+ } /*end of MCMC for normal */
+
+ if(*verbose)
+ Rprintf("100 percent done.\n");
+
+
+ /** write out the random seed **/
+ PutRNGstate();
+
+ /* Freeing the memory */
+ FreeMatrix(X, n_samp);
+ FreeMatrix(W, t_samp);
+ FreeMatrix(Wstar, t_samp);
+ free(n_grid);
+ FreeMatrix(S0, n_dim+1);
+ FreeMatrix(W1g, n_samp);
+ FreeMatrix(W2g, n_samp);
+ FreeMatrix(S_W, s_samp);
+ FreeMatrix(S_Wstar, s_samp);
+ free(mu);
+ FreeMatrix(Sigma, n_dim+1);
+ FreeMatrix(InvSigma, n_dim+1);
+ free(mu_w);
+ FreeMatrix(Sigma_w, n_dim);
+ FreeMatrix(InvSigma_w, n_dim);
+} /* main */
+
diff --git a/src/gibbsXDP.c b/src/gibbsXDP.c
new file mode 100644
index 0000000..bbc8e50
--- /dev/null
+++ b/src/gibbsXDP.c
@@ -0,0 +1,477 @@
+/******************************************************************
+ This file is a part of eco: R Package for Estimating Fitting
+ Bayesian Models of Ecological Inference for 2X2 tables
+ by Ying Lu and Kosuke Imai
+ Copyright: GPL version 2 or later.
+*******************************************************************/
+
+#include <stddef.h>
+#include <stdio.h>
+#include <math.h>
+#include <Rmath.h>
+#include <R_ext/Utils.h>
+#include "vector.h"
+#include "subroutines.h"
+#include "rand.h"
+#include "bayes.h"
+#include "sample.h"
+
+void cDPecoX(
+ /*data input */
+ double *pdX, /* data (X, Y) */
+ int *pin_samp, /* sample size */
+
+ /*MCMC draws */
+ int *n_gen, /* number of gibbs draws */
+ int *burn_in, /* number of draws to be burned in */
+ int *pinth, /* keep every nth draw */
+ int *verbose, /* 1 for output monitoring */
+ /* prior specification*/
+ int *pinu0, /* prior df parameter for InvWish */
+ double *pdtau0, /* prior scale parameter for Sigma under G0*/
+ double *mu0, /* prior mean for mu under G0 (3x1) */
+ double *pdS0, /* prior scale for Sigma (3x3) */
+
+ double *alpha0, /* precision parameter, can be fixed or updated*/
+ int *pinUpdate, /* 1 if alpha gets updated */
+ double *pda0, double *pdb0, /* prior for alpha if alpha updated*/
+
+ /*incorporating survey data */
+ int *survey, /*1 if survey data available(set of W_1, W_2,X)*/
+ /*0 otherwise*/
+ int *sur_samp, /*sample size of survey data*/
+ double *sur_W, /*set of known W_1, W_2 */
+
+ /*incorporating homeogenous areas */
+ int *x1, /* 1 if X=1 type areas available W_1 known, W_2 unknown */
+ int *sampx1, /* number X=1 type areas */
+ double *x1_W1, /* values of W_1 for X1 type areas */
+
+ int *x0, /* 1 if X=0 type areas available W_2 known, W_1 unknown */
+ int *sampx0, /* number X=0 type areas */
+ double *x0_W2, /* values of W_2 for X0 type areas */
+
+ /* bounds fo W1 */
+ double *minW1, double *maxW1,
+
+ /* flags */
+ int *parameter, /* 1 if save population parameter */
+ int *Grid, /* 1 if Grid algorithm is used; 0 for
+ Metropolis */
+
+ /* storage for Gibbs draws of mu/sigmat*/
+ double *pdSMu0, double *pdSMu1, double *pdSMu2,
+ double *pdSSig00, double *pdSSig01, double *pdSSig02,
+ double *pdSSig11, double *pdSSig12, double *pdSSig22,
+ /* storage for Gibbs draws of W*/
+ double *pdSW1, double *pdSW2,
+ /* storage for Gibbs draws of alpha */
+ double *pdSa,
+ /* storage for nstar at each Gibbs draw*/
+ int *pdSn
+ ){
+ /*some integers */
+ int n_samp = *pin_samp; /* sample size */
+ int s_samp = *sur_samp; /* sample size of survey data */
+ int x1_samp = *sampx1; /* sample size for X=1 */
+ int x0_samp = *sampx0; /* sample size for X=0 */
+ int t_samp = n_samp+x1_samp+x0_samp+s_samp; /* total sample size */
+ int nth = *pinth; /* keep every nth draw */
+ int n_dim = 2; /* dimension */
+ int n_step=1000; /* The default size of grid step */
+
+ /*prior parameters */
+ double tau0 = *pdtau0; /* prior scale */
+ int nu0 = *pinu0; /* prior degree of freedom*/
+ double **S0 = doubleMatrix(n_dim+1,n_dim+1);/*The prior S parameter for InvWish*/
+ double alpha = *alpha0; /* precision parameter*/
+ double a0 = *pda0, b0 = *pdb0; /* hyperprior for alpha */
+
+
+ /* data */
+ double **X = doubleMatrix(n_samp,n_dim); /* The Y and covariates */
+ double **W = doubleMatrix(t_samp,n_dim); /* The W1 and W2 matrix */
+ double **Wstar = doubleMatrix(t_samp,(n_dim+1)); /* The pseudo data */
+ double **S_W = doubleMatrix(s_samp,n_dim+1); /* The known W1 and W2,X */
+ double **S_Wstar = doubleMatrix(s_samp,n_dim+1);/* The logit
+ transformed S_W*/
+
+ /* grids */
+ double **W1g = doubleMatrix(n_samp, n_step); /* grids for W1 */
+ double **W2g = doubleMatrix(n_samp, n_step); /* grids for W2 */
+ int *n_grid = intArray(n_samp); /* grids size */
+
+ /* Model parameters */
+ /* Dirichlet variables */
+ double **mu = doubleMatrix(t_samp,(n_dim+1)); /* mean matrix */
+ double ***Sigma = doubleMatrix3D(t_samp,(n_dim+1),(n_dim+1)); /*covarince matrix*/
+ double ***InvSigma = doubleMatrix3D(t_samp,(n_dim+1),(n_dim+1)); /* inv of Sigma*/
+
+ /*conditional distribution parameter */
+ double **Sigma_w=doubleMatrix(n_dim,n_dim);
+ double **InvSigma_w=doubleMatrix(n_dim,n_dim);
+ double *mu_w=doubleArray(n_dim);
+
+ int nstar; /* # clusters with distict theta values */
+ int *C = intArray(t_samp); /* vector of cluster membership */
+ double *q = doubleArray(t_samp); /* Weights of posterior of Dirichlet */
+ double *qq = doubleArray(t_samp); /* cumulative weight vector of q */
+ double **S_tvt = doubleMatrix((n_dim+1),(n_dim+1)); /* S paramter for BVT in q0 */
+
+ /* variables defined in remixing step: cycle through all clusters */
+ double **Wstarmix = doubleMatrix(t_samp,(n_dim+1)); /*data matrix used */
+ double *mu_mix = doubleArray((n_dim+1)); /*updated MEAN parameter */
+ double **Sigma_mix = doubleMatrix((n_dim+1),(n_dim+1)); /*updated VAR parameter */
+ double **InvSigma_mix = doubleMatrix((n_dim+1),(n_dim+1)); /* Inv of Sigma_mix */
+ int nj; /* record # of obs in each cluster */
+ int *sortC = intArray(t_samp); /* record (sorted)original obs id */
+ int *indexC = intArray(t_samp); /* record original obs id */
+ int *label = intArray(t_samp); /* store index values */
+
+ /* misc variables */
+ int i, j, k, l, main_loop; /* used for various loops */
+ int itemp;
+ int itempA=0; /* counter for alpha */
+ int itempS=0; /* counter for storage */
+ int itempC=0; /* counter to control nth draw */
+ int progress = 1, itempP = ftrunc((double) *n_gen/10);
+ double dtemp, dtemp1, dtemp2;
+ double *vtemp = doubleArray((n_dim+1));
+ double **mtemp = doubleMatrix((n_dim+1),(n_dim+1));
+ double **mtemp1 = doubleMatrix((n_dim+1),(n_dim+1));
+ double **onedata = doubleMatrix(1, (n_dim+1));
+
+ /* get random seed */
+ GetRNGstate();
+
+ /* read priors under G0*/
+ itemp=0;
+ for(k=0;k<(n_dim+1);k++)
+ for(j=0;j<(n_dim+1);j++) S0[j][k]=pdS0[itemp++];
+
+ /* read the data set */
+ itemp = 0;
+ for (j = 0; j < n_dim; j++)
+ for (i = 0; i < n_samp; i++) X[i][j] = pdX[itemp++];
+
+ /*Intialize W, Wsatr for n_samp */
+ for (i=0; i< n_samp; i++) {
+ if (X[i][1]!=0 && X[i][1]!=1) {
+ W[i][0]=runif(minW1[i], maxW1[i]);
+ W[i][1]=(X[i][1]-X[i][0]*W[i][0])/(1-X[i][0]);
+ }
+ if (X[i][1]==0)
+ for (j=0; j<n_dim; j++) W[i][j]=0.0001;
+ if (X[i][1]==1)
+ for (j=0; j<n_dim; j++) W[i][j]=0.9999;
+ for (j=0; j<n_dim; j++)
+ Wstar[i][j]=log(W[i][j])-log(1-W[i][j]);
+ if (X[i][0]==0)
+ { Wstar[i][2]=log(0.0001/0.9999); }
+ else {
+ if (X[i][0]==1)
+ { Wstar[i][2]=log(0.9999/0.0001); }
+ else
+ { Wstar[i][2]=log(X[i][0]/(1-X[i][0])); }
+ }
+ }
+ /*read homeogenous areas information */
+ if (*x1==1)
+ for (i=0; i<x1_samp; i++) {
+ W[(n_samp+i)][0]=x1_W1[i];
+ if (W[(n_samp+i)][0]==0) W[(n_samp+i)][0]=0.0001;
+ if (W[(n_samp+i)][0]==1) W[(n_samp+i)][0]=0.9999;
+ Wstar[(n_samp+i)][0]=log(W[(n_samp+i)][0])-log(1-W[(n_samp+i)][0]);
+ }
+
+ if (*x0==1)
+ for (i=0; i<x0_samp; i++) {
+ W[(n_samp+x1_samp+i)][1]=x0_W2[i];
+ if (W[(n_samp+x1_samp+i)][1]==0) W[(n_samp+x1_samp+i)][1]=0.0001;
+ if (W[(n_samp+x1_samp+i)][1]==1) W[(n_samp+x1_samp+i)][1]=0.9999;
+ Wstar[(n_samp+x1_samp+i)][1]=log(W[(n_samp+x1_samp+i)][1])-log(1-W[(n_samp+x1_samp+i)][1]);
+ }
+
+ /*read the survey data */
+ if (*survey==1) {
+ itemp = 0;
+ for (j=0; j<=n_dim; j++)
+ for (i=0; i<s_samp; i++) {
+ S_W[i][j]=sur_W[itemp++];
+ if (S_W[i][j]==0) S_W[i][j]=0.0001;
+ if (S_W[i][j]==1) S_W[i][j]=0.9999;
+ S_Wstar[i][j]=log(S_W[i][j])-log(1-S_W[i][j]);
+ if (j<n_dim) {
+ W[(n_samp+x1_samp+x0_samp+i)][j]=S_W[i][j];
+ Wstar[(n_samp+x1_samp+x0_samp+i)][j]=S_Wstar[i][j];
+ }
+ else
+ Wstar[(n_samp+x1_samp+x0_samp+i)][j]=S_Wstar[i][j];
+ }
+ }
+
+
+ /* Calcualte grids */
+ if (*Grid)
+ GridPrep(W1g,W2g, X, maxW1, minW1, n_grid, n_samp, n_step);
+
+ /* parmeters for Trivaraite t-distribution-unchanged in MCMC */
+ for (j=0;j<=n_dim;j++)
+ for(k=0;k<=n_dim;k++)
+ mtemp[j][k]=S0[j][k]*(1+tau0)/(tau0*(nu0-n_dim+1));
+ dinv(mtemp, (n_dim+1), S_tvt);
+
+ /**draw initial values of mu_i, Sigma_i under G0 for all effective sample**/
+ /*1. Sigma_i under InvWish(nu0, S0^-1) with E(Sigma)=S0/(nu0-3)*/
+ /* InvSigma_i under Wish(nu0, S0^-1 */
+ /*2. mu_i|Sigma_i under N(mu0, Sigma_i/tau0) */
+ dinv(S0, (n_dim+1), mtemp);
+
+ for(i=0;i<t_samp;i++){
+ /*draw from wish(nu0, S0^-1) */
+ rWish(InvSigma[i], mtemp, nu0, (n_dim+1));
+ dinv(InvSigma[i], (n_dim+1), Sigma[i]);
+ for (j=0;j<=n_dim;j++)
+ for(k=0;k<=n_dim;k++) mtemp1[j][k]=Sigma[i][j][k]/tau0;
+ rMVN(mu[i], mu0, mtemp1, (n_dim+1));
+ }
+
+
+ /* initialize the cluster membership */
+ nstar=t_samp; /* the # of disticnt values */
+ for(i=0;i<t_samp;i++)
+ C[i]=i; /*cluster is from 0...n_samp-1 */
+
+ if (*verbose)
+ Rprintf("Starting Gibbs Sampler...\n");
+
+
+ for(main_loop=0; main_loop<*n_gen; main_loop++){
+ /**update W, Wstar given mu, Sigma only for the unknown W/Wstar**/
+ for (i=0; i<t_samp; i++){
+ for (j=0; j<n_dim; j++) {
+ mu_w[j]=mu[i][j]+Sigma[i][n_dim][j]/Sigma[i][n_dim][n_dim]*(Wstar[i][n_dim]-mu[i][n_dim]);
+ }
+ for (j=0; j<n_dim; j++)
+ for (k=0; k<n_dim; k++) {
+ Sigma_w[j][k]=Sigma[i][j][k]-Sigma[i][n_dim][j]/Sigma[i][n_dim][n_dim]*Sigma[i][n_dim][k];
+ }
+
+ dinv(Sigma_w, n_dim, InvSigma_w);
+
+
+ if (i<n_samp)
+ if (X[i][1]!=0 && X[i][1]!=1) {
+
+ /*1 project BVN(mu_i, Sigma_i) on the inth tomo line */
+ /*2 sample W_i on the ith tomo line */
+
+ if (*Grid)
+ rGrid(W[i], W1g[i], W2g[i], n_grid[i], mu_w, InvSigma_w, n_dim);
+ else {
+
+ rMH(W[i], X[i], minW1[i], maxW1[i], mu_w, InvSigma_w, n_dim);
+
+ }
+ }
+
+
+ Wstar[i][0]=log(W[i][0])-log(1-W[i][0]);
+ Wstar[i][1]=log(W[i][1])-log(1-W[i][1]);
+
+ if (*x1==1 && i>=n_samp && i<(n_samp+x1_samp)) {
+ dtemp=mu_w[1]+Sigma_w[0][1]/Sigma_w[0][0]*(Wstar[i][0]-mu_w[0]);
+ dtemp1=Sigma_w[1][1]*(1-Sigma_w[0][1]*Sigma_w[0][1]/(Sigma_w[0][0]*Sigma_w[1][1]));
+ Wstar[i][1]=norm_rand()*sqrt(dtemp1)+dtemp;
+ W[i][1]=exp(Wstar[i][1])/(1+exp(Wstar[i][1]));
+ }
+
+ /*update W1 given W2, mu_ord and Sigma_ord in x0 homeogeneous areas */
+
+ if (*x0==1 && i>=(n_samp+x1_samp) && i<(n_samp+x1_samp+x0_samp)) {
+ dtemp=mu_w[0]+Sigma_w[0][1]/Sigma_w[1][1]*(Wstar[i][1]-mu_w[1]);
+ dtemp1=Sigma_w[0][0]*(1-Sigma_w[0][1]*Sigma_w[0][1]/(Sigma_w[0][0]*Sigma_w[1][1]));
+ Wstar[i][0]=norm_rand()*sqrt(dtemp1)+dtemp;
+ W[i][0]=exp(Wstar[i][0])/(1+exp(Wstar[i][0]));
+ }
+ }
+
+ /**updating mu, Sigma given Wstar uisng effective sample size t_samp**/
+ for (i=0; i<t_samp; i++){
+ /* generate weight vector q */
+ dtemp=0;
+ for (j=0; j<t_samp; j++){
+ if (j!=i)
+ q[j]=dMVN(Wstar[i], mu[j], InvSigma[j], (n_dim+1), 0);
+ else
+ q[j]=alpha*dMVT(Wstar[i], mu0, S_tvt, (nu0-(n_dim+1)+1), (n_dim+1), 0);
+ dtemp+=q[j];
+ qq[j]=dtemp; /*compute qq, the cumlative of q*/
+ }
+ /*standardize q and qq */
+ for (j=0; j<t_samp; j++) {
+ qq[j]/=dtemp;
+ }
+
+ /** draw the configuration parameter **/
+
+ /* j=i means to draw from posterior baseline */
+ /* j=i' means to replace with i' obs */
+ j=0; dtemp=unif_rand();
+ while (dtemp > qq[j]) j++;
+ /** Dirichlet update Sigma_i, mu_i|Sigma_i **/
+ if (j==i){
+ onedata[0][0] = Wstar[i][0];
+ onedata[0][1] = Wstar[i][1];
+ onedata[0][2] = Wstar[i][2];
+ NIWupdate(onedata, mu[i], Sigma[i], InvSigma[i], mu0, tau0,nu0, S0, 1, n_dim+1);
+ C[i]=nstar;
+ nstar++;
+ }
+ else {
+ /*1. mu_i=mu_j, Sigma_i=Sigma_j*/
+ /*2. update C[i]=C[j] */
+ for(k=0;k<=n_dim;k++) {
+ mu[i][k]=mu[j][k];
+ for(l=0;l<=n_dim;l++) {
+ Sigma[i][k][l]=Sigma[j][k][l];
+ InvSigma[i][k][l]=InvSigma[j][k][l];
+ }
+ }
+ C[i]=C[j];
+ }
+ sortC[i]=C[i];
+ } /* end of i loop*/
+ /** remixing step using effective sample**/
+ for(i=0;i<t_samp;i++)
+ indexC[i]=i;
+ R_qsort_int_I(sortC, indexC, 1, t_samp);
+
+ nstar=0;
+ i=0;
+ while (i<t_samp){
+ j=sortC[i]; /*saves the first element in a block of same values */
+ nj=0; /* counter for a block of same values */
+
+ /* get data for remixing */
+ while ((sortC[i]==j) && (i<t_samp)) {
+ label[nj]=indexC[i];
+ for (k=0; k<=n_dim; k++) {
+ Wstarmix[nj][k]=Wstar[label[nj]][k];
+ }
+ nj++;
+ i++;
+ } /* i records the current position in IndexC */
+ /* nj records the # of obs in Psimix */
+
+ /** posterior update for mu_mix, Sigma_mix based on Psimix **/
+ NIWupdate(Wstarmix, mu_mix,Sigma_mix, InvSigma_mix, mu0, tau0, nu0, S0, nj, (n_dim+1));
+
+ /**update mu, Simgat with mu_mix, Sigmat_mix via label**/
+ for (j=0;j<nj;j++){
+ C[label[j]]=nstar; /*updating C vector with no gap */
+ for (k=0; k<=n_dim; k++){
+ mu[label[j]][k]=mu_mix[k];
+ for (l=0;l<=n_dim;l++){
+ Sigma[label[j]][k][l]=Sigma_mix[k][l];
+ InvSigma[label[j]][k][l]=InvSigma_mix[k][l];
+ }
+ }
+ }
+ nstar++; /*finish update one distinct value*/
+ } /* nstar is the number of distinct values */
+
+ /** updating alpha **/
+ if(*pinUpdate) {
+ dtemp1=(double)(alpha+1);
+ dtemp2=(double)t_samp;
+ dtemp=b0-log(rbeta(dtemp1, dtemp2));
+
+ dtemp1=(double)(a0+nstar-1)/(t_samp*dtemp);
+ if(unif_rand() < dtemp1) {
+ dtemp2=(double)(a0+nstar);
+ alpha=rgamma(dtemp2, 1/dtemp);
+ }
+ else {
+ dtemp2=(double)(a0+nstar-1);
+ alpha=rgamma(dtemp2, 1/dtemp);
+ }
+ }
+
+ /*store Gibbs draws after burn_in */
+ R_CheckUserInterrupt();
+ if (main_loop>=*burn_in) {
+ itempC++;
+ if (itempC==nth){
+ if(*pinUpdate) {
+ pdSa[itempA]=alpha;
+ pdSn[itempA]=nstar;
+ itempA++;
+ }
+
+ for(i=0; i<(n_samp+x1_samp+x0_samp); i++) {
+ pdSMu0[itempS]=mu[i][0];
+ pdSMu1[itempS]=mu[i][1];
+ pdSMu2[itempS]=mu[i][2];
+ pdSSig00[itempS]=Sigma[i][0][0];
+ pdSSig01[itempS]=Sigma[i][0][1];
+ pdSSig02[itempS]=Sigma[i][0][2];
+ pdSSig11[itempS]=Sigma[i][1][1];
+ pdSSig12[itempS]=Sigma[i][1][2];
+ pdSSig22[itempS]=Sigma[i][2][2];
+ pdSW1[itempS]=W[i][0];
+ pdSW2[itempS]=W[i][1];
+ itempS++;
+ }
+ itempC=0;
+ }
+ }
+ if (*verbose)
+ if (itempP == main_loop) {
+ Rprintf("%3d percent done.\n", progress*10);
+ itempP+=ftrunc((double) *n_gen/10); progress++;
+ R_FlushConsole();
+ }
+ } /*end of MCMC for DP*/
+
+if (*verbose)
+ Rprintf("100 percent done.\n");
+
+ /** write out the random seed **/
+ PutRNGstate();
+
+ /* Freeing the memory */
+ FreeMatrix(S0, n_dim+1);
+ FreeMatrix(X, n_samp);
+ FreeMatrix(W, t_samp);
+ FreeMatrix(Wstar, t_samp);
+ FreeMatrix(S_W, s_samp);
+ FreeMatrix(S_Wstar, s_samp);
+ FreeMatrix(W1g, n_samp);
+ FreeMatrix(W2g, n_samp);
+ free(n_grid);
+ FreeMatrix(mu, t_samp);
+ Free3DMatrix(Sigma, t_samp,n_dim+1);
+ Free3DMatrix(InvSigma, t_samp, n_dim+1);
+ free(mu_w);
+ FreeMatrix(Sigma_w, n_dim);
+ FreeMatrix(InvSigma_w, n_dim);
+ free(C);
+ free(q);
+ free(qq);
+ FreeMatrix(S_tvt, n_dim+1);
+ FreeMatrix(Wstarmix, t_samp);
+ free(mu_mix);
+ FreeMatrix(Sigma_mix, n_dim+1);
+ FreeMatrix(InvSigma_mix, n_dim+1);
+ free(sortC);
+ free(indexC);
+
+ free(label);
+ free(vtemp);
+ FreeMatrix(mtemp, n_dim+1);
+ FreeMatrix(mtemp1, n_dim+1);
+ free(onedata);
+} /* main */
+
+
diff --git a/src/preBaseX.c b/src/preBaseX.c
new file mode 100644
index 0000000..75918e3
--- /dev/null
+++ b/src/preBaseX.c
@@ -0,0 +1,84 @@
+/******************************************************************
+ This file is a part of eco: R Package for Fitting Bayesian Models
+ of Ecological Inference for 2x2 Tables
+ by Kosuke Imai and Ying Lu
+ Copyright: GPL version 2 or later.
+*******************************************************************/
+
+#include <stddef.h>
+#include <stdio.h>
+#include <math.h>
+#include <Rmath.h>
+#include <R_ext/Utils.h>
+#include "vector.h"
+#include "subroutines.h"
+#include "rand.h"
+#include "bayes.h"
+#include "sample.h"
+
+/* Conditional Prediction for Normal Parametric Model for 2x2 Tables */
+void preBaseX(
+ double *X, /* data X */
+ double *pdmu,
+ double *pdSigma,
+ int *pin_samp,
+ int *pin_draw,
+ int *verbose, /* 1 for output monitoring */
+ double *pdStore
+ ){
+
+ /* some integers */
+ int n_samp = *pin_samp; /* sample size */
+ int n_draw = *pin_draw; /* sample size of survey data */
+ int n_dim = 2;
+
+ double *mu = doubleArray(n_dim); /* The mean */
+ double *Wstar = doubleArray(n_dim);
+ double **Sigma = doubleMatrix(n_dim, n_dim); /* The covariance matrix */
+
+ /* misc variables */
+ int i, j, k, main_loop; /* used for various loops */
+ int itemp=0;
+ int itempM=0;
+ int itempS=0;
+ int progress = 1, itempP = ftrunc((double) n_draw/10);
+
+ /* get random seed */
+ GetRNGstate();
+
+ for(main_loop=0; main_loop<n_draw; main_loop++){
+ Sigma[0][0] = pdSigma[itempS]-pdSigma[itempS+2]*pdSigma[itempS+2]/pdSigma[itempS+5];
+ Sigma[1][1] = pdSigma[itempS+3]-pdSigma[itempS+4]*pdSigma[itempS+4]/pdSigma[itempS+5];
+ Sigma[0][1] = pdSigma[itempS+1]-pdSigma[itempS+2]*pdSigma[itempS+4]/pdSigma[itempS+5];
+ Sigma[1][0] = Sigma[0][1];
+ for(i=0; i<n_samp; i++) {
+ mu[0] = pdmu[itempM]+pdSigma[itempS+2]/pdSigma[itempS+5]*(X[i]-pdmu[itempM+2]);
+ mu[1] = pdmu[itempM+1]+pdSigma[itempS+4]/pdSigma[itempS+5]*(X[i]-pdmu[itempM+2]);
+ rMVN(Wstar, mu, Sigma, n_dim);
+ for (j=0; j<n_dim; j++)
+ pdStore[itemp++] = exp(Wstar[j])/(1+exp(Wstar[j]));
+ }
+ itempS += 6;
+ itempM += 3;
+ if (*verbose)
+ if (itempP == main_loop) {
+ Rprintf("%3d percent done.\n", progress*10);
+ itempP+=ftrunc((double) n_draw/10); progress++;
+ R_FlushConsole();
+ }
+ R_CheckUserInterrupt();
+ }
+
+ if(*verbose)
+ Rprintf("100 percent done.\n");
+
+ /** write out the random seed **/
+ PutRNGstate();
+
+ /* Freeing the memory */
+ free(mu);
+ free(Wstar);
+ FreeMatrix(Sigma,n_dim);
+
+} /* main */
+
diff --git a/src/preDP.c b/src/preDP.c
new file mode 100644
index 0000000..b56aff8
--- /dev/null
+++ b/src/preDP.c
@@ -0,0 +1,83 @@
+/******************************************************************
+ This file is a part of eco: R Package for Fitting Bayesian Models
+ of Ecological Inference for 2x2 Tables
+ by Kosuke Imai and Ying Lu
+ Copyright: GPL version 2 or later.
+*******************************************************************/
+
+#include <stddef.h>
+#include <stdio.h>
+#include <math.h>
+#include <Rmath.h>
+#include <R_ext/Utils.h>
+#include "vector.h"
+#include "subroutines.h"
+#include "rand.h"
+#include "bayes.h"
+#include "sample.h"
+
+/* Prediction for Nonparametric Model for 2x2 Tables */
+void preDP(
+ double *pdmu,
+ double *pdSigma,
+ int *pin_samp,
+ int *pin_draw,
+ int *pin_dim,
+ int *verbose, /* 1 for output monitoring */
+ double *pdStore
+ ){
+
+ /* some integers */
+ int n_samp = *pin_samp; /* sample size */
+ int n_draw = *pin_draw; /* sample size of survey data */
+ int n_dim = *pin_dim; /* dimension */
+
+ double *mu = doubleArray(n_dim); /* The mean */
+ double *Wstar = doubleArray(n_dim);
+ double **Sigma = doubleMatrix(n_dim, n_dim); /* The covariance matrix */
+
+ /* misc variables */
+ int i, j, k, main_loop; /* used for various loops */
+ int itemp = 0;
+ int itempM = 0;
+ int itempS = 0;
+ int progress = 1, itempP = ftrunc((double) n_draw/10);
+
+ /* get random seed */
+ GetRNGstate();
+
+ for(main_loop=0; main_loop<n_draw; main_loop++){
+ for(i=0; i<n_samp; i++) {
+ for (j=0;j<n_dim;j++) {
+ mu[j] = pdmu[itempM++];
+ for (k=j;k<n_dim;k++) {
+ Sigma[j][k] = pdSigma[itempS++];
+ Sigma[k][j] = Sigma[j][k];
+ }
+ }
+ rMVN(Wstar, mu, Sigma, n_dim);
+ for (j=0; j<n_dim; j++)
+ pdStore[itemp++] = exp(Wstar[j])/(1+exp(Wstar[j]));
+ }
+ if (*verbose)
+ if (itempP == main_loop) {
+ Rprintf("%3d percent done.\n", progress*10);
+ itempP+=ftrunc((double) n_draw/10); progress++;
+ R_FlushConsole();
+ }
+ R_CheckUserInterrupt();
+ }
+
+ if(*verbose)
+ Rprintf("100 percent done.\n");
+
+ /** write out the random seed **/
+ PutRNGstate();
+
+ /* Freeing the memory */
+ free(mu);
+ free(Wstar);
+ FreeMatrix(Sigma,n_dim);
+
+} /* main */
+
diff --git a/src/preDPX.c b/src/preDPX.c
new file mode 100644
index 0000000..f913bc5
--- /dev/null
+++ b/src/preDPX.c
@@ -0,0 +1,85 @@
+/******************************************************************
+ This file is a part of eco: R Package for Fitting Bayesian Models
+ of Ecological Inference for 2x2 Tables
+ by Kosuke Imai and Ying Lu
+ Copyright: GPL version 2 or later.
+*******************************************************************/
+
+#include <stddef.h>
+#include <stdio.h>
+#include <math.h>
+#include <Rmath.h>
+#include <R_ext/Utils.h>
+#include "vector.h"
+#include "subroutines.h"
+#include "rand.h"
+#include "bayes.h"
+#include "sample.h"
+
+/* Conditional Prediction for Nonparametric Model for 2x2 Tables */
+void preDPX(
+ double *pdmu,
+ double *pdSigma,
+ double *X,
+ int *pin_samp,
+ int *pin_draw,
+ int *pin_dim,
+ int *verbose, /* 1 for output monitoring */
+ double *pdStore
+ ){
+
+ /* some integers */
+ int n_samp = *pin_samp; /* sample size */
+ int n_draw = *pin_draw; /* sample size of survey data */
+ int n_dim = *pin_dim; /* dimension */
+
+ double *mu = doubleArray(n_dim); /* The mean */
+ double *Wstar = doubleArray(n_dim);
+ double **Sigma = doubleMatrix(n_dim, n_dim); /* The covariance matrix */
+
+ /* misc variables */
+ int i, j, k, main_loop; /* used for various loops */
+ int itemp = 0;
+ int itempM = 0;
+ int itempS = 0;
+ int progress = 1, itempP = ftrunc((double) n_draw/10);
+
+ /* get random seed */
+ GetRNGstate();
+
+ for(main_loop=0; main_loop<n_draw; main_loop++){
+ for(i=0; i<n_samp; i++) {
+ mu[0] = pdmu[itempM]+pdSigma[itempS+2]/pdSigma[itempS+5]*(X[i]-pdmu[itempM+2]);
+ mu[1] = pdmu[itempM+1]+pdSigma[itempS+4]/pdSigma[itempS+5]*(X[i]-pdmu[itempM+2]);
+ Sigma[0][0] = pdSigma[itempS]-pdSigma[itempS+2]*pdSigma[itempS+2]/pdSigma[itempS+5];
+ Sigma[1][1] = pdSigma[itempS+3]-pdSigma[itempS+4]*pdSigma[itempS+4]/pdSigma[itempS+5];
+ Sigma[0][1] = pdSigma[itempS+1]-pdSigma[itempS+2]*pdSigma[itempS+4]/pdSigma[itempS+5];
+ Sigma[1][0] = Sigma[0][1];
+ rMVN(Wstar, mu, Sigma, n_dim);
+ for (j=0; j<n_dim; j++)
+ pdStore[itemp++] = exp(Wstar[j])/(1+exp(Wstar[j]));
+ itempS += 6;
+ itempM += 3;
+ }
+ if (*verbose)
+ if (itempP == main_loop) {
+ Rprintf("%3d percent done.\n", progress*10);
+ itempP+=ftrunc((double) n_draw/10); progress++;
+ R_FlushConsole();
+ }
+ R_CheckUserInterrupt();
+ }
+
+ if(*verbose)
+ Rprintf("100 percent done.\n");
+
+ /** write out the random seed **/
+ PutRNGstate();
+
+ /* Freeing the memory */
+ free(mu);
+ free(Wstar);
+ FreeMatrix(Sigma,n_dim);
+
+} /* main */
+
diff --git a/src/rand.c b/src/rand.c
new file mode 100644
index 0000000..00aea1e
--- /dev/null
+++ b/src/rand.c
@@ -0,0 +1,184 @@
+/******************************************************************
+ This file is a part of eco: R Package for Fitting Bayesian Models
+ of Ecological Inference for 2x2 Tables
+ by Kosuke Imai and Ying Lu
+ Copyright: GPL version 2 or later.
+*******************************************************************/
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include <Rmath.h>
+#include <R_ext/Utils.h>
+#include <R.h>
+#include "vector.h"
+#include "subroutines.h"
+#include "rand.h"
+#include "sample.h"
+
+/* Multivariate Normal density */
+double dMVN(
+ double *Y, /* The data */
+ double *MEAN, /* The parameters */
+ double **SIG_INV, /* inverse of the covariance matrix */
+ int dim, /* dimension */
+ int give_log){ /* 1 if log_scale 0 otherwise */
+
+ int j,k;
+ double value=0.0;
+
+ for(j=0;j<dim;j++){
+ for(k=0;k<j;k++)
+ value+=2*(Y[k]-MEAN[k])*(Y[j]-MEAN[j])*SIG_INV[j][k];
+ value+=(Y[j]-MEAN[j])*(Y[j]-MEAN[j])*SIG_INV[j][j];
+ }
+
+ value=-0.5*value-0.5*dim*log(2*M_PI)+0.5*ddet(SIG_INV, dim, 1);
+
+
+ if(give_log)
+ return(value);
+ else
+ return(exp(value));
+
+}
+
+/* the density of Multivariate T-distribution */
+double dMVT(
+ double *Y, /* The data */
+ double *MEAN, /* mean */
+ double **SIG_INV, /* inverse of scale matrix */
+ int nu, /* Degrees of freedom */
+ int dim, /* dimension */
+ int give_log) /* 1 if log_scale 0 otherwise */
+{
+ int j,k;
+ double value=0;
+
+ for(j=0;j<dim;j++){
+ for(k=0;k<j;k++)
+ value+=2*(Y[k]-MEAN[k])*(Y[j]-MEAN[j])*SIG_INV[j][k];
+ value+=(Y[j]-MEAN[j])*(Y[j]-MEAN[j])*SIG_INV[j][j];
+ }
+
+ value=0.5*ddet(SIG_INV, dim,1) - 0.5*dim*(log((double)nu)+log(M_PI)) -
+ 0.5*((double)dim+nu)*log(1+value/(double)nu) +
+ lgammafn(0.5*(double)(nu+dim)) - lgammafn(0.5*(double)nu);
+
+ if(give_log)
+ return(value);
+ else
+ return(exp(value));
+}
+
+
+/* Sample from the MVN dist */
+void rMVN(
+ double *Sample, /* Vector for the sample */
+ double *mean, /* The vector of means */
+ double **Var, /* The matrix Variance */
+ int size) /* The dimension */
+{
+ int j,k;
+ double **Model = doubleMatrix(size+1, size+1);
+ double cond_mean;
+
+ /* draw from mult. normal using SWP */
+ for(j=1;j<=size;j++){
+ for(k=1;k<=size;k++)
+ Model[j][k]=Var[j-1][k-1];
+ Model[0][j]=mean[j-1];
+ Model[j][0]=mean[j-1];
+ }
+ Model[0][0]=-1;
+ Sample[0]=(double)norm_rand()*sqrt(Model[1][1])+Model[0][1];
+ for(j=2;j<=size;j++){
+ SWP(Model,j-1,size+1);
+ cond_mean=Model[j][0];
+ for(k=1;k<j;k++) cond_mean+=Sample[k-1]*Model[j][k];
+ Sample[j-1]=(double)norm_rand()*sqrt(Model[j][j])+cond_mean;
+ }
+
+ FreeMatrix(Model,size+1);
+}
+
+
+/* Sample from a wish dist */
+/* Odell, P. L. and Feiveson, A. H. ``A Numerical Procedure to Generate
+ a Sample Covariance Matrix'' Journal of the American Statistical
+ Association, Vol. 61, No. 313. (Mar., 1966), pp. 199-203. */
+
+void rWish(
+ double **Sample, /* The matrix with to hold the sample */
+ double **S, /* The parameter */
+ int df, /* the degrees of freedom */
+ int size) /* The dimension */
+{
+ int i,j,k;
+ double *V = doubleArray(size);
+ double **B = doubleMatrix(size, size);
+ double **C = doubleMatrix(size, size);
+ double **N = doubleMatrix(size, size);
+ double **mtemp = doubleMatrix(size, size);
+
+ for(i=0;i<size;i++) {
+ V[i]=rchisq((double) df-i-1);
+ B[i][i]=V[i];
+ for(j=(i+1);j<size;j++)
+ N[i][j]=norm_rand();
+ }
+
+ for(i=0;i<size;i++) {
+ for(j=i;j<size;j++) {
+ Sample[i][j]=0;
+ Sample[j][i]=0;
+ mtemp[i][j]=0;
+ mtemp[j][i]=0;
+ if(i==j) {
+ if(i>0)
+ for(k=0;k<j;k++)
+ B[j][j]+=N[k][j]*N[k][j];
+ }
+ else {
+ B[i][j]=N[i][j]*sqrt(V[i]);
+ if(i>0)
+ for(k=0;k<i;k++)
+ B[i][j]+=N[k][i]*N[k][j];
+ }
+ B[j][i]=B[i][j];
+ }
+ }
+
+ dcholdc(S, size, C);
+ for(i=0;i<size;i++)
+ for(j=0;j<size;j++)
+ for(k=0;k<size;k++)
+ mtemp[i][j]+=C[i][k]*B[k][j];
+ for(i=0;i<size;i++)
+ for(j=0;j<size;j++)
+ for(k=0;k<size;k++)
+ Sample[i][j]+=mtemp[i][k]*C[j][k];
+
+ free(V);
+ FreeMatrix(B, size);
+ FreeMatrix(C, size);
+ FreeMatrix(N, size);
+ FreeMatrix(mtemp, size);
+}
+
+/* Sample from a Dirichlet distribution */
+void rDirich(
+ double *Sample, /* Vector for the sample */
+ double *theta, /* parameters */
+ int size) /* The dimension */
+{
+ int j;
+ double dtemp=0;
+
+ for (j=0; j<size; j++) {
+ Sample[j] = rgamma(theta[j], 1.0);
+ dtemp += Sample[j];
+ }
+ for (j=0 ; j<size; j++)
+ Sample[j] /= dtemp;
+}
diff --git a/src/rand.h b/src/rand.h
new file mode 100644
index 0000000..ed5b569
--- /dev/null
+++ b/src/rand.h
@@ -0,0 +1,13 @@
+/******************************************************************
+ This file is a part of eco: R Package for Fitting Bayesian Models
+ of Ecological Inference for 2x2 Tables
+ by Kosuke Imai and Ying Lu
+ Copyright: GPL version 2 or later.
+*******************************************************************/
+
+double dMVN(double *Y, double *MEAN, double **SIGMA, int dim, int give_log);
+double dMVT(double *Y, double *MEAN, double **SIG_INV, int nu, int dim, int give_log);
+void rMVN(double *Sample, double *mean, double **inv_Var, int size);
+void rWish(double **Sample, double **S, int df, int size);
+void rDirich(double *Sample, double *theta, int size);
+
diff --git a/src/sample.c b/src/sample.c
new file mode 100644
index 0000000..a264211
--- /dev/null
+++ b/src/sample.c
@@ -0,0 +1,237 @@
+/******************************************************************
+ This file is a part of eco: R Package for Fitting Bayesian Models
+ of Ecological Inference for 2x2 Tables
+ by Kosuke Imai and Ying Lu
+ Copyright: GPL version 2 or later.
+*******************************************************************/
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include <Rmath.h>
+#include <R_ext/Utils.h>
+#include <R.h>
+#include "vector.h"
+#include "subroutines.h"
+#include "rand.h"
+
+
+/* Grid method samping from tomography line*/
+void rGrid(
+ double *Sample, /* W_i sampled from each tomography line */
+ double *W1gi, /* The grid lines of W1[i] */
+ double *W2gi, /* The grid lines of W2[i] */
+ int ni_grid, /* number of grids for observation i*/
+ double *mu, /* mean vector for normal */
+ double **InvSigma, /* Inverse covariance matrix for normal */
+ int n_dim) /* dimension of parameters */
+{
+ int j;
+ double dtemp;
+ double *vtemp=doubleArray(n_dim);
+ double *prob_grid=doubleArray(ni_grid); /* density by grid */
+ double *prob_grid_cum=doubleArray(ni_grid); /* cumulative density by grid */
+
+ dtemp=0;
+ for (j=0;j<ni_grid;j++){
+ vtemp[0]=log(W1gi[j])-log(1-W1gi[j]);
+ vtemp[1]=log(W2gi[j])-log(1-W2gi[j]);
+ prob_grid[j]=dMVN(vtemp, mu, InvSigma, n_dim, 1) -
+ log(W1gi[j])-log(W2gi[j])-log(1-W1gi[j])-log(1-W2gi[j]);
+ prob_grid[j]=exp(prob_grid[j]);
+ dtemp+=prob_grid[j];
+ prob_grid_cum[j]=dtemp;
+ }
+ for (j=0;j<ni_grid;j++)
+ prob_grid_cum[j]/=dtemp; /*standardize prob.grid */
+
+ /*2 sample W_i on the ith tomo line */
+ j=0;
+ dtemp=unif_rand();
+ while (dtemp > prob_grid_cum[j]) j++;
+ Sample[0]=W1gi[j];
+ Sample[1]=W2gi[j];
+
+ free(vtemp);
+ free(prob_grid);
+ free(prob_grid_cum);
+
+}
+
+/* preparation for Grid */
+void GridPrep(
+ double **W1g, /* grids holder for W1 */
+ double **W2g, /* grids holder for W2 */
+ double **X, /* data: [X Y] */
+ double *maxW1, /* upper bound for W1 */
+ double *minW1, /* lower bound for W1 */
+ int *n_grid, /* number of grids */
+ int n_samp, /* sample size */
+ int n_step /* step size */
+)
+{
+ int i, j;
+ double dtemp;
+ double *resid = doubleArray(n_samp);
+
+ for(i=0; i<n_samp; i++)
+ for (j=0; j<n_step; j++){
+ W1g[i][j]=0;
+ W2g[i][j]=0;
+ }
+ for(i=0;i<n_samp;i++) {
+ if (X[i][1]!=0 && X[i][1]!=1) {
+ /* 1/n_step is the length of the grid */
+ dtemp=(double)1/n_step;
+ if ((maxW1[i]-minW1[i]) > (2*dtemp)) {
+ n_grid[i]=ftrunc((maxW1[i]-minW1[i])*n_step);
+ resid[i]=(maxW1[i]-minW1[i])-n_grid[i]*dtemp;
+ /*if (maxW1[i]-minW1[i]==1) resid[i]=dtemp/4; */
+ j=0;
+ while (j<n_grid[i]) {
+ W1g[i][j]=minW1[i]+(j+1)*dtemp-(dtemp+resid[i])/2;
+ if ((W1g[i][j]-minW1[i])<resid[i]/2) W1g[i][j]+=resid[i]/2;
+ if ((maxW1[i]-W1g[i][j])<resid[i]/2) W1g[i][j]-=resid[i]/2;
+ W2g[i][j]=(X[i][1]-X[i][0]*W1g[i][j])/(1-X[i][0]);
+ j++;
+ }
+ }
+ else {
+ W1g[i][0]=minW1[i]+(maxW1[i]-minW1[i])/3;
+ W2g[i][0]=(X[i][1]-X[i][0]*W1g[i][0])/(1-X[i][0]);
+ W1g[i][1]=minW1[i]+2*(maxW1[i]-minW1[i])/3;
+ W2g[i][1]=(X[i][1]-X[i][0]*W1g[i][1])/(1-X[i][0]);
+ n_grid[i]=2;
+ }
+ }
+ }
+
+ free(resid);
+}
+
+/* sample W via MH for 2x2 table */
+void rMH(
+ double *W, /* previous draws */
+ double *XY, /* X_i and Y_i */
+ double W1min, /* lower bound for W1 */
+ double W1max, /* upper bound for W1 */
+ double *mu, /* mean vector for normal */
+ double **InvSigma, /* Inverse covariance matrix for normal */
+ int n_dim) /* dimension of parameters */
+{
+ int j;
+ double dens1, dens2, ratio;
+ double *Sample = doubleArray(n_dim);
+ double *vtemp = doubleArray(n_dim);
+ double *vtemp1 = doubleArray(n_dim);
+
+ /* sample W_1 from unif(W1min, W1max) */
+ Sample[0] = runif(W1min, W1max);
+ Sample[1] = XY[1]/(1-XY[0])-Sample[0]*XY[0]/(1-XY[0]);
+ for (j = 0; j < n_dim; j++) {
+ vtemp[j] = log(Sample[j])-log(1-Sample[j]);
+ vtemp1[j] = log(W[j])-log(1-W[j]);
+ }
+ /* acceptance ratio */
+ dens1 = dMVN(vtemp, mu, InvSigma, n_dim, 1) -
+ log(Sample[0])-log(Sample[1])-log(1-Sample[0])-log(1-Sample[1]);
+ dens2 = dMVN(vtemp1, mu, InvSigma, n_dim, 1) -
+ log(W[0])-log(W[1])-log(1-W[0])-log(1-W[1]);
+ ratio = fmin2(1, exp(dens1-dens2));
+
+ /* accept */
+ if (unif_rand() < ratio)
+ for (j=0; j<n_dim; j++)
+ W[j]=Sample[j];
+
+ free(Sample);
+ free(vtemp);
+ free(vtemp1);
+}
+
+
+/* sample W via MH for 2xC table */
+void rMH2c(
+ double *W, /* W */
+ double *X, /* X_i */
+ double Y, /* Y_i */
+ double *minU, /* lower bound for U */
+ double *maxU, /* upper bound for U */
+ double *mu, /* mean vector for normal */
+ double **InvSigma, /* Inverse covariance matrix for normal */
+ int n_dim, /* dimension of parameters */
+ int maxit, /* max number of iterations for
+ rejection sampling */
+ int reject) /* if 1, use rejection sampling to
+ draw from the truncated Dirichlet
+ if 0, use Gibbs sampling
+ */
+{
+ int iter = 100; /* number of Gibbs iterations */
+ int i, j, exceed;
+ double dens1, dens2, ratio, dtemp;
+ double *Sample = doubleArray(n_dim);
+ double *param = doubleArray(n_dim);
+ double *vtemp = doubleArray(n_dim);
+ double *vtemp1 = doubleArray(n_dim);
+
+ /* set parent Dirichlet parameter to 1 */
+ for (j = 0; j < n_dim; j++)
+ param[j] = 1.0;
+
+ /* Sample a candidate draw of W from truncated Dirichlet */
+ if (reject) { /* rejection sampling */
+ i = 0; exceed = 1;
+ while (exceed > 0) {
+ rDirich(vtemp, param, n_dim);
+ exceed = 0;
+ for (j = 0; j < n_dim; j++)
+ if (vtemp[j] > maxU[j] || vtemp[j] < minU[j])
+ exceed++;
+ i++;
+ if (i > maxit)
+ error("rMH2c: rejection algorithm failed because bounds are too tight.\n increase maxit or use gibbs sampler instead.");
+ }
+ }
+ else { /* gibbs sampler */
+ for (j = 0; j < n_dim; j++)
+ vtemp[j] = W[j]*X[j]/Y;
+ for (i = 0; i < iter; i++) {
+ dtemp = vtemp[n_dim-1];
+ for (j = 0; j < n_dim-1; j++) {
+ dtemp += vtemp[j];
+ vtemp[j] = runif(fmax2(minU[j], dtemp-maxU[n_dim-1]),
+ fmin2(maxU[j], dtemp-minU[n_dim-1]));
+ dtemp -= vtemp[j];
+ }
+ vtemp[n_dim-1] = dtemp;
+ }
+ }
+ /* calcualte W and its logit transformation */
+ for (j = 0; j < n_dim; j++) {
+ Sample[j] = vtemp[j]*Y/X[j];
+ vtemp[j] = log(Sample[j])-log(1-Sample[j]);
+ vtemp1[j] = log(W[j])-log(1-W[j]);
+ }
+
+ /* acceptance ratio */
+ dens1 = dMVN(vtemp, mu, InvSigma, n_dim, 1);
+ dens2 = dMVN(vtemp1, mu, InvSigma, n_dim, 1);
+ for (j=0; j<n_dim; j++) {
+ dens1 -= (log(Sample[j])+log(1-Sample[j]));
+ dens2 -= (log(W[j])+log(1-W[j]));
+ }
+ ratio=fmin2(1, exp(dens1-dens2));
+
+ /* accept */
+ if (unif_rand() < ratio)
+ for (j = 0; j < n_dim; j++)
+ W[j] = Sample[j];
+
+ free(Sample);
+ free(param);
+ free(vtemp);
+ free(vtemp1);
+}
+
+
diff --git a/src/sample.h b/src/sample.h
new file mode 100644
index 0000000..37fd4e6
--- /dev/null
+++ b/src/sample.h
@@ -0,0 +1,16 @@
+/******************************************************************
+ This file is a part of eco: R Package for Fitting Bayesian Models
+ of Ecological Inference for 2x2 Tables
+ by Kosuke Imai and Ying Lu
+ Copyright: GPL version 2 or later.
+*******************************************************************/
+
+void rGrid(double *Sample, double *W1gi, double *W2gi, int ni_grid,
+ double *mu, double **InvSigma, int n_dim);
+void GridPrep(double **W1g, double **W2g, double **X, double *maxW1,
+ double *minW1, int *n_grid, int n_samp, int n_step);
+void rMH(double *W, double *XY, double W1min, double W1max,
+ double *mu, double **InvSigma, int n_dim);
+void rMH2c(double *W, double *X, double Y, double *minU,
+ double *maxU, double *mu, double **InvSigma, int n_dim,
+ int maxit, int reject);
diff --git a/src/subroutines.c b/src/subroutines.c
new file mode 100644
index 0000000..a3e0989
--- /dev/null
+++ b/src/subroutines.c
@@ -0,0 +1,122 @@
+/******************************************************************
+ This file is a part of eco: R Package for Fitting Bayesian Models
+ of Ecological Inference for 2x2 Tables
+ by Kosuke Imai and Ying Lu
+ Copyright: GPL version 2 or later.
+*******************************************************************/
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+#include <R_ext/Utils.h>
+#include <R.h>
+#include "vector.h"
+#include "rand.h"
+
+/* The Sweep operator */
+void SWP(
+ double **X, /* The Matrix to work on */
+ int k, /* The row to sweep */
+ int size) /* The dim. of X */
+{
+ int i,j;
+
+ if (X[k][k] < 10e-20)
+ error("SWP: singular matrix.\n");
+ else
+ X[k][k]=-1/X[k][k];
+ for(i=0;i<size;i++)
+ if(i!=k){
+ X[i][k]=-X[i][k]*X[k][k];
+ X[k][i]=X[i][k];
+ }
+ for(i=0;i<size;i++)
+ for(j=0;j<size;j++)
+ if(i!=k && j!=k)
+ X[i][j]=X[i][j]+X[i][k]*X[k][j]/X[k][k];
+
+}
+
+
+/* inverting a matrix */
+void dinv(double **X,
+ int size,
+ double **X_inv)
+{
+ int i,j, k, errorM;
+ double *pdInv = doubleArray(size*size);
+
+ for (i = 0, j = 0; j < size; j++)
+ for (k = 0; k <= j; k++)
+ pdInv[i++] = X[k][j];
+ F77_CALL(dpptrf)("U", &size, pdInv, &errorM);
+ if (!errorM) {
+ F77_CALL(dpptri)("U", &size, pdInv, &errorM);
+ if (errorM) {
+ Rprintf("LAPACK dpptri failed, %d\n", errorM);
+ error("Exiting from dinv().\n");
+ }
+ }
+ else {
+ Rprintf("LAPACK dpptrf failed, %d\n", errorM);
+ error("Exiting from dinv().\n");
+ }
+ for (i = 0, j = 0; j < size; j++) {
+ for (k = 0; k <= j; k++) {
+ X_inv[j][k] = pdInv[i];
+ X_inv[k][j] = pdInv[i++];
+ }
+ }
+
+ free(pdInv);
+}
+
+
+/* Cholesky decomposition */
+/* returns lower triangular matrix */
+void dcholdc(double **X, int size, double **L)
+{
+ int i, j, k, errorM;
+ double *pdTemp = doubleArray(size*size);
+
+ for (j = 0, i = 0; j < size; j++)
+ for (k = 0; k <= j; k++)
+ pdTemp[i++] = X[k][j];
+ F77_CALL(dpptrf)("U", &size, pdTemp, &errorM);
+ if (errorM) {
+ Rprintf("LAPACK dpptrf failed, %d\n", errorM);
+ error("Exiting from dcholdc().\n");
+ }
+ for (j = 0, i = 0; j < size; j++) {
+ for (k = 0; k < size; k++) {
+ if(j<k)
+ L[j][k] = 0.0;
+ else
+ L[j][k] = pdTemp[i++];
+ }
+ }
+
+ free(pdTemp);
+}
+
+/* calculate the determinant of the positive definite symmetric matrix
+ using the Cholesky decomposition */
+double ddet(double **X, int size, int give_log)
+{
+ int i;
+ double logdet=0.0;
+ double **pdTemp = doubleMatrix(size, size);
+
+ dcholdc(X, size, pdTemp);
+ for(i = 0; i < size; i++)
+ logdet += log(pdTemp[i][i]);
+
+ FreeMatrix(pdTemp, size);
+ if(give_log)
+ return(2.0*logdet);
+ else
+ return(exp(2.0*logdet));
+}
+
+
diff --git a/src/subroutines.h b/src/subroutines.h
new file mode 100644
index 0000000..528a266
--- /dev/null
+++ b/src/subroutines.h
@@ -0,0 +1,11 @@
+/******************************************************************
+ This file is a part of eco: R Package for Fitting Bayesian Models
+ of Ecological Inference for 2x2 Tables
+ by Kosuke Imai and Ying Lu
+ Copyright: GPL version 2 or later.
+*******************************************************************/
+
+void SWP( double **X, int k, int size);
+void dinv(double **X, int size, double **X_inv);
+void dcholdc(double **X, int size, double **L);
+double ddet(double **X, int size, int give_log);
diff --git a/src/vector.c b/src/vector.c
new file mode 100644
index 0000000..8c62425
--- /dev/null
+++ b/src/vector.c
@@ -0,0 +1,103 @@
+/******************************************************************
+ This file is a part of eco: R Package for Fitting Bayesian Models
+ of Ecological Inference for 2x2 Tables
+ by Kosuke Imai and Ying Lu
+ Copyright: GPL version 2 or later.
+*******************************************************************/
+
+#include <stdlib.h>
+#include <assert.h>
+#include <stdio.h>
+#include <R_ext/Utils.h>
+#include <R.h>
+
+int* intArray(int num) {
+ int *iArray = (int *)malloc(num * sizeof(int));
+ if (iArray)
+ return iArray;
+ else
+ error("Out of memory error in intArray\n");
+}
+
+int** intMatrix(int row, int col) {
+ int i;
+ int **iMatrix = (int **)malloc(row * sizeof(int *));
+ if (iMatrix) {
+ for (i = 0; i < row; i++) {
+ iMatrix[i] = (int *)malloc(col * sizeof(int));
+ if (!iMatrix[i])
+ error("Out of memory error in intMatrix\n");
+ }
+ return iMatrix;
+ }
+ else
+ error("Out of memory error in intMatrix\n");
+}
+
+double* doubleArray(int num) {
+ double *dArray = (double *)malloc(num * sizeof(double));
+ if (dArray)
+ return dArray;
+ else
+ error("Out of memory error in doubleArray\n");
+}
+
+double** doubleMatrix(int row, int col) {
+ int i;
+ double **dMatrix = (double **)malloc((size_t)(row * sizeof(double *)));
+ if (dMatrix) {
+ for (i = 0; i < row; i++) {
+ dMatrix[i] = (double *)malloc((size_t)(col * sizeof(double)));
+ if (!dMatrix[i])
+ error("Out of memory error in doubleMatrix\n");
+ }
+ return dMatrix;
+ }
+ else
+ error("Out of memory error in doubleMatrix\n");
+}
+
+double*** doubleMatrix3D(int x, int y, int z) {
+ int i;
+ double ***dM3 = (double ***)malloc(x * sizeof(double **));
+ if (dM3) {
+ for (i = 0; i < x; i++)
+ dM3[i] = doubleMatrix(y, z);
+ return dM3;
+ }
+ else
+ error("Out of memory error in doubleMatrix3D\n");
+}
+
+long* longArray(int num) {
+ long *lArray = (long *)malloc(num * sizeof(long));
+ if (lArray)
+ return lArray;
+ else
+ error("Out of memory error in longArray\n");
+}
+
+void FreeMatrix(double **Matrix, int row) {
+ int i;
+ for (i = 0; i < row; i++)
+ free(Matrix[i]);
+ free(Matrix);
+}
+
+void FreeintMatrix(int **Matrix, int row) {
+ int i;
+ for (i = 0; i < row; i++)
+ free(Matrix[i]);
+ free(Matrix);
+}
+
+void Free3DMatrix(double ***Matrix, int index, int row) {
+ int i;
+ for (i = 0; i < index; i++)
+ FreeMatrix(Matrix[i], row);
+ free(Matrix);
+}
+
+
+
+
diff --git a/src/vector.h b/src/vector.h
new file mode 100644
index 0000000..ba6fe0c
--- /dev/null
+++ b/src/vector.h
@@ -0,0 +1,22 @@
+/******************************************************************
+ This file is a part of eco: R Package for Fitting Bayesian Models
+ of Ecological Inference for 2x2 Tables
+ by Kosuke Imai and Ying Lu
+ Copyright: GPL version 2 or later.
+*******************************************************************/
+
+#include <stdlib.h>
+#include <assert.h>
+
+int *intArray(int num);
+int **intMatrix(int row, int col);
+
+double *doubleArray(int num);
+double **doubleMatrix(int row, int col);
+double ***doubleMatrix3D(int x, int y, int z);
+
+long *longArray(int num);
+
+void FreeMatrix(double **Matrix, int row);
+void FreeintMatrix(int **Matrix, int row);
+void Free3DMatrix(double ***Matrix, int index, int row);
--
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/debian-science/packages/r-cran-eco.git
More information about the debian-science-commits
mailing list