Mixed Normal Optimization {cubfits}

Mixed Normal Optimization

Description

Constrained optimization for mixed normal in 1D and typically for 2 components.

Usage

  mixnormerr.optim(X, K = 2, param = NULL)
  dmixnormerr(x, param)

Arguments

X

a gene expression data matrix of dimension N * R which has N genes and R replicates.

K

number of components to fit.

x

vector of quantiles.

param

parameters of mixnormerr, typically the element param of the mixnormerr.optim() returning object.

Details

The function mixnormerr.optim() maximizes likelihood using constrOptim() based on the gene expression data X (usually in log scale) for N genes and R replicates (NA is allowed). The likelihood of each gene expression is a K = 2 component mixed normal distribution (\sum_k p_k N(mu_k, \sigma_k^2 + \sigma_e^2)) with measurement errors of the replicates (N(0, \sigma_e^2)).

The sigma_k^2 is as the error of random component and the sigma_e^2 is as the error of fixed component. Both are within a mixture model of two normal distributions.

The function dmixnormerr() computes the density of the mixed normal distribution.

param is a parameter list and contains five elements: K for number of components, prop for proportions, mu for centers of components, sigma2 for variance of components, and sigma2.e for variance of measurement errors.

Value

mixnormerr.optim() returns a list containing three main elements param is the final results (MLEs), param.start is the starting parameters, and optim.ret is the original returns of constrOptim().

Note

This function is limited for small K. An equivalent EM algorithm should be done in a more stable way for large K.

Author(s)

Wei-Chen Chen [e-mail].

References

https://github.com/snoweye/cubfits/

See Also

print.mixnormerr(), simu.mixnormerr().

Examples

## Not run: 
suppressMessages(library(cubfits, quietly = TRUE))

### Get individual of phi.Obs.
GM <- apply(yassour[, -1], 1, function(x) exp(mean(log(x[x != 0]))))
phi.Obs.all <- yassour[, -1] / sum(GM) * 15000
phi.Obs.all[phi.Obs.all == 0] <- NA

### Run optimization.
X <- log(as.matrix(phi.Obs.all))
param.init <- list(K = 2, prop = c(0.95, 0.05), mu = c(-0.59, 3.11),
                   sigma2 = c(1.40, 0.59), sigma2.e = 0.03)
ret <- mixnormerr.optim(X, K = 2, param = param.init)
print(ret)

## End(Not run)