eBayes.R

# eBayes.R

# functions

data.fun <- function(y, K=10, eps=1e-08) {

# y is a N x M matrix

# N is the number of observations

# M is the number of samples.

# takes a set of samples and creates

# a matrix of sample probabilities N x K

# a vector of outcomes K x 1

# a vector of samples sizes M x 1

min_y <- min(y, na.rm = TRUE)

max_y <- max(y, na.rm = TRUE)

diff_y <- (max_y - min_y)/K

p_hat <- matrix(NA,dim(y)[2],K)

y_hat <- matrix(NA,K,1)

min_t <- min_y - eps

max_t <- min_t + diff_y + eps

for (k in 1:K) {

p_hat[,k] <- colMeans(y > min_t & y <= max_t, na.rm = TRUE)

y_hat[k] <- (min_t + max_t)/2

min_t <- max_t - eps

max_t <- min_t + diff_y + eps

}

N_hat <- colSums(is.na(y)==0)

return(list(p_hat=p_hat, y_hat=y_hat, N_hat=N_hat))

}

mixmod.fun <- function(p_hat, N_hat, R=10, eps=1e-8, maxiter=500,

noise=1, verb=FALSE) {

# this estimator is based on the npEM algorithm in the

# mixtools library.

p_hat <- as.matrix(p_hat) # matrix of observed probabilities.

M <- dim(p_hat)[1] # number of samples

N_hat <- matrix(rep(N_hat,dim(p_hat)[2]),nrow=M)

# Number of observations in each sample.

index_p <- round(runif(R,min=1,max=M))

mu2 <- exp(log(t(p_hat[index_p,]) + eps) +

matrix(rnorm(R*dim(p_hat)[2],mean=0,sd=noise),

nrow=dim(p_hat)[2]))

mu <- t(t(mu2)/rowSums(t(mu2)))

# finds the set of multinomial distributions that are to be used.

# Initial values

log_stuff <- lfactorial(N_hat[,1]) - rowSums(lfactorial(p_hat*N_hat))

Lf_hat_fixed <- log_stuff + (p_hat*N_hat)%*%log(mu)

# calculates the likelihood function using the multinomial function.

R <- dim(mu)[2]

Lf_hat_fixed <- log_stuff + (p_hat*N_hat)%*%log(mu)

iter <- 1

finished <- FALSE

lambda <- matrix(1/R,1,R)

while (!finished) {

iter <- iter + 1

lambda_old <- lambda

exp1 <- exp(Lf_hat_fixed + log(t(matrix(rep(lambda,M),ncol=M))))

lambda <- colMeans(exp1/rowSums(exp1,na.rm = TRUE))

finished <- iter >= maxiter

maxchange <- sum((lambda - lambda_old)^2)

finished <- finished | (maxchange < eps)

if (verb) {

print(iter)

print(maxchange)

}

}

exp1 <- exp(Lf_hat_fixed + log(t(matrix(rep(lambda,M),ncol=M))))

z_hat <- exp1/rowSums(exp1,na.rm = TRUE)

return(list(posteriors = z_hat, iter= iter, mu=mu,

lambdahat = lambda))

}