Create likelihood matrix from simulated data set

Simulate a data set, then compute the conditional likelihood matrix under a univariate normal likelihood and a mixture-of-normals prior. This models a simple nonparametric Empirical Bayes method applied to simulated data.

simulatemixdata(
  n,
  m,
  simtype = c("n", "nt"),
  log = FALSE,
  normalize.rows = !log
)

Arguments

n	Positive integer specifying the number of samples to generate and, consequently, the number of rows of the likelihood matrix L.
m	Integer 2 or greater specifying the number of mixture components.
simtype	The type of data to simulate. If simtype = "n", simulate n random numbers from a mixture of three univariate normals with mean zero and standard deviation 1, 3 and 6. If simtype = "nt", simulate from a mixture of three univariate normals (with zero mean and standard deviations 1, 3 and 5), and a t-distribution with 2 degrees of freedom.
log	If log = TRUE, return the
normalize.rows	If normalize.rows = TRUE, normalize the rows of the likelihood matrix so that the largest entry in each row is 1. The maximum-likelihood estimate of the mixture weights should be invariant to this normalization, and can improve the numerical stability of the optimization.

Value

simulatemixdata returns a list with three list elements:

x

The vector of simulated random numbers (it has length n).

s

The standard deviations of the mixture components in the mixture-of-normals prior. The rules for selecting the standard deviations are based on the autoselect.mixsd function from the ashr package.

L

The n x m conditional likelihood matrix, in which individual entries (i,j) of the likelihood matrix are given by the normal density function with mean zero and variance 1 + s[j]^2. If normalize.rows = TRUE, the entries in each row are normalized such that the larger entry in each row is 1. If log = TRUE, the matrix of log-likelihoods is returned.

Examples

# Generate the likelihood matrix for a data set with 1,000 samples
# and a nonparametric Empirical Bayes model with 20 mixture
# components.
dat <- simulatemixdata(1000,20)