# normgam.fit: Normal-gamma Maximum Likelihood Estimator In NormalGamma: Normal-gamma convolution model

## Description

Computes the Maximum Likelihood Estimator for the normal-gamma distribution, either from a normal-gamma distributed sample or from two samples respectively normal-gamma and normally distributed.

## Usage

 ```1 2``` ```normgam.fit(X, N = NULL, par.init = NULL, lower = NULL, upper = NULL, control = NULL, verbose = FALSE) ```

## Arguments

 `X` vector of normal-gamma distributed values. `N` vector of normal distributed values. `par.init` vector of initial values for parameters (optional). `par.init[1]` and `par.init[2]` are the mean and standard deviation of the normal distribution, and `par.init[3]` and `par.init[4]` are the shape and scale parameters of the gamma distribution. See details for default initial values. `lower, upper` Bounds on the variables for maximization (optional). `control` list of control parameters (see details). `verbose` logical; if `TRUE` initial values of the parameters are printed.

## Details

Likelihood maximization is run by the `R` function `optimx`.

By default, maximization is run with the following control parameters: the maximum number of iterations is equal to `1000` and the vector of scaling values for the parameters is `(par0[1],` ` par0[2],` `par0[3]*par0[4],` `sqrt(par0[3])*par0[4])/10` where `par0` is the vector of default initial parameters. In case of unsuccessful convergence, maximization is run with `optimx` default control parameters. A list of control parameters can also be chosen by the user (see `optimx`).

If `par.init == NULL`, the initial parameters are computed in two ways depending if `N` is `NULL` or not. If `N != NULL`, the initial parameters are computed following the method of the moments (see Plancade S., Rozenholc Y. and Lund E., BMC Bioinfo 2012). If `N == NULL`, the initial parameters `(par0[1], par0[2])` of the normal distribution are computed following the RMA procedure of Xie Y., Wang X. and Story M. (2009) for the normal-exponential convolution model, and the initial parameters of the gamma distribution, computed following the method of the moments, are `(par0[4]=``(sd(X)^2-par0[2])/``(mean(X)-par0[1]),` `par0[3]=``(mean(X)-``par0[1]/par0[4])`. Note that the RMA procedure for initial parameter computation when `N == NULL` stems from an heuristic adapted to microarray data. For parameters with different magnitude, user should specify initial parameters.

## Value

 `par ` vector of estimated parameters; `par[1]` and `par[2]` are the mean and standard deviation of the normal distribution and `par[3]` and `par[4]` are the shape and scale parameters of the gamma distribution `lik ` value of the normal-gamma log-likelihood corresponding to `par`. `conv ` integer code: `0` indicates successful convergence. This parameter has the value of the output parameter `conv` from the procedure `optimx` used for likelihood maximization (see `optimx` for details).

## References

Plancade S., Rozenholc Y. and Lund E. "Generalization of the normal-exponential model : exploration of a more accurate parametrisation for the signal distribution on Illumina BeadArrays", BMC Bioinfo 2012, 13(329).

Xie Y. Wang X. and Story M. (2009), Statistical methods of background correction for Illumina BeadArray data, Bioinformatics, 25(6), 751–757.

`dnormgam` computes the density of the normal-gamma distribution and `normgam.signal` implements the background correction using the normal-gamma model. `Intensities` provides an example of Illumina microarray data.
 ``` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54``` ```# Example 1: simulated data ## Not run: n = 1000 par = c(60,5,0.15,400) X = rnorm(n, mean=par[1], sd=par[2]) + rgamma(n, shape=par[3], scale=par[4]) N = rnorm(100, mean=par[1], sd=par[2]) par1 = normgam.fit(X, N)\$par par2 = normgam.fit(X)\$par F1 = dnormgam(par1, plot=FALSE) F2 = dnormgam(par2, plot=FALSE) par(mfrow=c(2,1)) H = histogram(X, type='irregular', verbose=FALSE, plot=FALSE) plot(H, xlim=c(0,500)) lines(F1\$xout, F1\$dout,col='red') plot(H, xlim=c(0,500)) lines(F2\$xout, F2\$dout,col='blue') ## End(Not run) # Example 2: Illumina data ## Not run: data(RegNegIntensities_Example) X = Intensities\$Regular N = Intensities\$Negative par1 = normgam.fit(X, N)\$par par2 = normgam.fit(X)\$par F1 = dnormgam(par1, plot=FALSE) F2 = dnormgam(par2, plot=FALSE) par(mfrow=c(2,1)) H = histogram(X, type='irregular', verbose=FALSE, plot=FALSE) plot(H, xlim=c(0,500)) lines(F1\$xout, F1\$dout, col='red') plot(H, xlim=c(0,500)) lines(F2\$xout, F2\$dout, col='blue') ## End(Not run) ```