u5_mixnorm: Mixed Normal Optimization
In cubfits: Codon Usage Bias Fits

Description Usage Arguments Details Value Note Author(s) References See Also Examples

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

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

`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.

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.

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().

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

Wei-Chen Chen wccsnow@gmail.com.

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

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

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