# 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))
}