NormalGamma (version 1.1)

normgam.fit: Normal-gamma Maximum Likelihood Estimator

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

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.

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

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.

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.

See Also

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.

Examples

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