sefm {mixbox}R Documentation

Approximating the asymptotic standard error for parameters of the finite mixture models based on the observed Fisher information matrix.

Description

The density function of each finite mixture model can be represented as

M(yΨ)=g=1GωgfY(y,Θg), {\cal{M}}(\bold{y}|\bold{\Psi})=\sum_{g=1}^{G} \omega_{g} f_{\bold{Y}}(\bold{y}, \bold{\Theta}_g),

where positive constants ω1,ω2,,ωG\omega_{1}, \omega_{2},\cdots,\omega_{G} are called weight (or mixing proportions) parameters with this properties that g=1Gωg=1\sum_{g=1}^{G}\omega_{g}=1 and Ψ=(Θ1,,ΘG)\bold{\Psi} = \bigl(\bold{\Theta}_{1},\cdots, \bold{\Theta}_{G}\bigr)^{\top} with Θg=(ωg,μg,Σg,λg)\bold{\Theta}_g=\bigl({\bold{\omega}}_g, {\bold{\mu}}_g, {{\Sigma}}_g, {\bold{\lambda}}_g\bigr)^{\top}. Herein, fY(y,Θg)f_{\bold{Y}}(\bold{y}, \bold{\Theta}_g) accounts for the density function of random vector Y\bold{Y} within gg-th component that admits the representation given by

Y=dμg+WλgZ0+WΣg12Z1, {\bold{Y}} \mathop=\limits^d {\boldsymbol{\mu}}_g+\sqrt{W}{\boldsymbol{\lambda}}_g\vert{Z}_0\vert + \sqrt{W}{\Sigma}_g^{\frac{1}{2}}{\bold{Z}}_1,

where μgRd{\bold{\mu}_g} \in {R}^{d} is location vector, λgRd{\bold{\lambda}_g} \in {R}^{d} is skewness vector, Σg\Sigma_g is a positive definite symmetric dispersion matrix for g=1,,Gg=1,\cdots,G. Further, WW is a positive random variable with mixing density function fW(wθg)f_W(w| \bold{\theta}_g), Z0N(0,1){Z}_{0}\sim N({0},1), and Z1Nd(0,Σg){\bold{Z}}_{1}\sim N_{d}\bigl({\bold{0}}, \Sigma_g\bigr). We note that WW, Z0Z_0, and Z1{\bold{Z}}_{1} are mutually independent. For approximating the asymptotic standard error for parameters of the finite mixture model based on observed Fisher information matrix, we use the method of Basford et al. (1997). In fact, the covariance matrix of maximum likelihood (ML) estimator Ψ^\hat{\bold{\Psi}}, can be approximated by the inverse of the observed information matrix as

i=1nh^ih^i, \sum_{i = 1}^{n} \hat{{\bold{h}}}_{i} \hat{{\bold{h}}}^{\top}_{i},

where

h^i=ΨlogLi(Ψ^), \hat{\bold{h}}_{i} = \frac{\partial}{\partial \bold{\Psi} } \log L_{i}(\hat{\bold{\Psi}}),

and logL(Ψ^)=i=1nlogLi(Ψ^)=i=1nlog{g=1Gωg^fY(yiΘg^)} \log L(\hat{\bold{\Psi}}) = \sum_{i =1}^{n} \log L_{i}(\hat{\bold{\Psi}}) = \sum_{i =1}^{n} \log \Bigl\{ \sum_{g=1}^{G} \widehat{\omega_{g}} f_{\bold{Y}}\bigl({\bold{y}}_{i}|\widehat{\bold{\Theta}_{g}}\bigr)\Bigr\}. Herein ω^g\widehat{\omega}_g and Θg^\widehat{\bold{\Theta}_g}, for g=1,,Gg=1,\cdots,G, denote the ML estimator of ωg\omega_g and Θg\bold{\Theta}_g, respectively.

Usage

sefm(Y, G, weight, model = "restricted", mu, sigma, lambda, family = "constant",
    skewness = "FALSE", param = NULL, theta = NULL, tick = NULL, h = 0.001, N = 3000,
  level = 0.05, PDF = NULL)

Arguments

Y

an n×dn\times d matrix of observations gives the coordinates of the data points.

G

number of components.

weight

a vector of weight parameters (or mixing proportions).

model

it must be "canonical", "restricted", or "unrestricted". By default model = "restricted".

mu

a list of location vectors of G components.

sigma

a list of dispersion matrices of G components.

lambda

a list of skewness vectors of G components. If model is either "canonical" or "unrestricted", then skewness vactor must be given in matrix form of appropriate size.

family

name of mixing distribution. By default family = "constant" that corresponds to the finite mixture of multivariate normal (or skew normal) distribution. Other candidates for family name are: "bs" (for Birnbaum-Saunders), "burriii" (for Burr type iii), "chisq" (for chi-square), "exp" (for exponential), "f" (for Fisher), "gamma" (for gamma), "gig" (for generalized inverse-Gaussian), "igamma" (for inverse-gamma), "igaussian" (for inverse-Gaussian), "lindley" (for Lindley), "loglog" (for log-logistic), "lognorm" (for log-normal), "lomax" (for Lomax), "pstable" (for positive α\alpha-stable), "ptstable" (for polynomially tilted α\alpha-stable), "rayleigh" (for Rayleigh), and "weibull" (for Weibull).

skewness

a logical statement. By default skewness = "FALSE" which means that a symmetric model is fitted to each component (cluster). If skewness = "TRUE", then a skewed model is fitted to each component.

param

name of the elements of θ\bold{\theta} as the parameter vector of mixing distribution with density function fW(wθ)f_W(w| \bold{\theta}). By default it is NULL.

PDF

mathematical expression for mixing density function fW(wθ)f_W(w\vert \bold{\theta}). By default it is NULL.

theta

a list of maximum likelihood estimator for θ\bold{\theta} across G components. By default it is NULL.

tick

a binary vector whose length depends on type of family. The elements of tick are either 0 or 1. If element of tick is 0, then the corresponding element of θ\bold{\theta} is not considered in the formula of fW(wθ)f_W(w\vert{\bold{\theta)}} for computing the required posterior expectations. If element of tick is 1, then the corresponding element of θ\bold{\theta} is considered in the formula of fW(wθ)f_W(w|{\bold{\theta)}}. For instance, if family = "gamma" and either its shape or rate parameter is one, then tick = c(1). This is while, if family = "gamma" and both of the shape and rate parameters are in the formula of fW(wθ)f_W(w\vert{\bold{\theta)}}, then tick = c(1, 1). By default tick = NULL.

h

a positive small value for computing numerical derivative of fW(wθ)f_W(w\vert \bold{\theta}) with respect to θ\bold{\theta}, that is /θfW(wθ)\partial/ \partial \theta f_W(w\vert \bold{\theta}). By default h=0.001h = 0.001.

N

an integer number for approximating the posterior expected values within the E-step of the EM algorithm through the Monte Carlo method. By default N=3000N = 3000.

level

significance level α\alpha for constructing 100(1α)%100(1-\alpha)\% confidence interval. By default α=0.05\alpha = 0.05.

Details

Mathematical expressions for density function of mixing distributions fW(wθ)f_W(w\vert{\bold{\theta}}), are "bs" (for Birnbaum-Saunders), "burriii" (for Burr type iii), "chisq" (for chi-square), "exp" (for exponential), "f" (for Fisher), "gamma" (for gamma), "gig" (for generalized inverse-Gaussian), "igamma" (for inverse-gamma), "igaussian" (for inverse-Gaussian), "lindley" (for Lindley), "loglog" (for log-logistic), "lognorm" (for log-normal), "lomax" (for Lomax), "pstable" (for positive α\alpha-stable), "ptstable" (for polynomially tilted α\alpha-stable), "rayleigh" (for Rayleigh), and "weibull" (for Weibull). We note that the density functions of "pstable" and "ptstable" families have no closed form and so are not represented here. The pertinent and given by the following, respectively.

fbs(wθ)=wβ+βw22παwexp{12α2[wβ+βw2]},f_{bs}(w\vert{\bold{\theta}}) = \frac {\sqrt{\frac{w}{\beta}}+\sqrt{\frac {\beta}{w}}}{2\sqrt{2\pi}\alpha w}\exp\Biggl\{-\frac {1}{2\alpha^2}\Bigl[\frac{w}{\beta}+\frac{\beta}{w}-2\Bigr]\Biggr\},

where θ=(α,β){\bold{\theta}}=(\alpha,\beta)^{\top}. Herein α>0\alpha> 0 and β>0\beta> 0 are the first and second parameters of this family, respectively.

fburrii(wθ)=αβwβ1(1+wβ)α1,f_{burrii}(w\vert {\bold{\theta}}) = \alpha \beta w^{-\beta-1} \bigl( 1+w^{-\beta} \bigr) ^{-\alpha-1},

where w>0w>0 and θ=(α,β){\bold{\theta}}=(\alpha, \beta)^{\top}. Herein α>0\alpha> 0 and β>0\beta> 0 are the first and second parameters of this family, respectively.

fchisq(wθ)=2α2Γ(α2)wα21exp{w2}, f_{chisq}(w\vert{{\theta}}) = \frac{2^{-\frac {\alpha}{2}}}{\Gamma\bigl(\frac{\alpha}{2}\bigr)} w^{\frac{\alpha}{2}-1}\exp\Bigl\{-\frac {w}{2} \Bigr\},

where w>0w>0 and θ=α{{\theta}}=\alpha. Herein α>0\alpha> 0 is the degrees of freedom parameter of this family.

fexp(wθ)=αexp{αw},f_{exp}(w\vert{{\theta}}) =\alpha \exp \bigl\{-\alpha w\bigr\},

where w>0w>0 and θ=α{{\theta}}=\alpha where α>0\alpha> 0 is the rate parameter of this family.

ff(wθ)=B1(α2,β2)(αβ)α2wα21(1+αwβ)(α+β2), f_{f}(w\vert{\bold{\theta}}) = B^{-1}\Bigl(\frac {\alpha}{2}, \frac {\beta}{2}\Bigr)\Bigl( \frac {\alpha}{\beta} \Bigr)^{\frac {\alpha}{2}} w^{\frac {\alpha}{2}-1}\Bigl(1 + \alpha\frac {w}{\beta} \Bigr)^{-\left(\frac {\alpha+\beta}{2} \right)},

where w>0w>0 and B(.,.)B(.,.) denotes the ordinary beta function. Herein θ=(α,β){\bold{\theta}}=(\alpha, \beta)^{\top} where α>0\alpha> 0 and β>0\beta> 0 are the first and second degrees of freedom parameters of this family, respectively.

fgamma(wθ)=βαΓ(α)(wβ)α1exp{βw}, f_{gamma}(w\vert{\bold{\theta}}) = \frac {\beta^{\alpha}}{\Gamma(\alpha)} \Bigl( \frac{w}{\beta}\Bigr)^{\alpha-1}\exp\bigl\{ - \beta w \bigr\},

where w>0w>0 and θ=(α,β){\bold{\theta}}=(\alpha, \beta)^{\top}. Herein α>0\alpha> 0 and β>0\beta> 0 are the shape and rate parameters of this family, respectively.

fgig(wθ)=12Kα(βδ)(βδ)α/2wα1exp{δ2wβw2}, f_{gig}(w\vert{\bold{\theta}}) =\frac{1}{2{\cal{K}}_{\alpha}( \sqrt{\beta \delta})}\Bigl(\frac{\beta}{\delta}\Bigr)^{\alpha/2}w^{\alpha-1} \exp\biggl\{-\frac{\delta}{2w}-\frac{\beta w}{2}\biggr\},

where Kα(.){\cal{K}}_{\alpha}(.) denotes the modified Bessel function of the third kind with order index α\alpha and θ=(α,δ,β){\bold{\theta}}=(\alpha, \delta, \beta)^{\top}. Herein <α<+-\infty <\alpha <+\infty, δ>0\delta> 0, and β>0\beta> 0 are the first, second, and third parameters of this family, respectively.

figamma(wθ)=1Γ(α)(wβ)α1exp{βw}, f_{igamma}(w\vert{\bold{\theta}}) = \frac{1}{\Gamma(\alpha)} \Bigl( \frac{w}{\beta}\Bigr)^{-\alpha-1}\exp\Bigl\{ - \frac{\beta}{w} \Bigr\},

where w>0w>0 and θ=(α,β){\bold{\theta}}=(\alpha, \beta)^{\top}. Herein α>0\alpha> 0 and β>0\beta> 0 are the shape and scale parameters of this family, respectively.

figaussian(wθ)=β2πw3exp{β(wα)22α2w}, f_{igaussian}(w\vert{\bold{\theta}}) =\sqrt{\frac{\beta}{2 \pi w^3}} \exp\biggl\{-\frac{\beta(w - \alpha)^2}{2\alpha^2 w}\biggr\},

where w>0w>0 and θ=(α,β){\bold{\theta}}=(\alpha, \beta)^{\top}. Herein α>0\alpha>0 and β>0\beta> 0 are the first (mean) and second (shape) parameters of this family, respectively.

flidley(wθ)=α2α+1(1+w)exp{αw},f_{lidley}(w\vert{{\theta}}) =\frac{\alpha^2}{\alpha+1} (1+w)\exp \bigl\{-\alpha w\bigr\},

where w>0w>0 and θ=α{{\theta}}=\alpha where α>0\alpha> 0 is the only parameter of this family.

floglog(wθ)=αβαwα1[(wβ)α+1]2, f_{loglog}(w\vert{\bold{\theta}}) =\frac{\alpha}{ \beta^{\alpha}} w^{\alpha-1}\left[ \Bigl( \frac {w}{\beta}\Bigr)^\alpha +1\right] ^{-2},

where w>0w>0 and θ=(α,β){\bold{\theta}}=(\alpha, \beta)^{\top}. Herein α>0\alpha> 0 and β>0\beta> 0 are the shape and scale (median) parameters of this family, respectively.

flognorm(wθ)=(2πσw)1exp{12(logwμσ)2}, f_{lognorm}(w\vert{\bold{\theta}}) = \bigl(\sqrt{2\pi} \sigma w \bigr)^{-1} \exp\biggl\{ -\frac{1}{2}\left( \frac {\log w - \mu}{\sigma}\right) ^2\biggr\},

where w>0w>0 and θ=(μ,σ){\bold{\theta}}=(\mu, \sigma)^{\top}. Herein <μ<+-\infty<\mu<+\infty and σ>0\sigma> 0 are the first and second parameters of this family, respectively.

flomax(wθ)=αβ(1+βw)(α+1), f_{lomax}(w\vert{\bold{\theta}}) = \alpha \beta \bigl( 1+\beta w\bigr)^{-(\alpha+1)},

where w>0w>0 and θ=(α,β){\bold{\theta}}=(\alpha, \beta)^{\top}. Herein α>0\alpha>0 and β>0\beta> 0 are the shape and rate parameters of this family, respectively.

frayleigh(wθ)=2wβ2exp{(wβ)2}, f_{rayleigh}(w\vert{{\theta}}) = 2\frac {w}{\beta^2}\exp\biggl\{ -\Bigl( \frac {w}{\beta}\Bigr)^2 \biggr\},

where w>0w>0 and θ=β{{\theta}}=\beta. Herein β>0\beta>0 is the scale parameter of this family.

fweibull(wθ)=αβ(wβ)α1exp{(wβ)α}, f_{weibull}(w\vert{\bold{\theta}}) = \frac {\alpha}{\beta}\Bigl( \frac {w}{\beta} \Bigr)^{\alpha - 1}\exp\biggl\{ -\Bigl( \frac{w}{\beta}\Bigr)^\alpha \biggr\},

where w>0w>0 and θ=(α,β){\bold{\theta}}=(\alpha, \beta)^{\top}. Herein α>0\alpha>0 and β>0\beta> 0 are the shape and scale parameters of this family, respectively.

In what follows, we give four examples. In the first, second, and third examples, we consider three mixture models including: two-component normal, two-component restricted skew tt, and two-component restricted skew sub-Gaussian α\alpha-stable (SSG) mixture models are fitted to iris, AIS, and bankruptcy data, respectively. In order to approximate the asymptotic standard error of the model parameters, the ML estimators for parameters of skew tt and SSG mixture models have been computed through the R packages EMMIXcskew (developed by Lee and McLachlan (2018) for skew tt) and mixSSG (developed by Teimouri (2022) for skew sub-Gaussian α\alpha-stable). To avoid running package mixSSG, we use the ML estimators correspond to bankruptcy data provided by Teimouri (2022). The package mixSSG is available at https://CRAN.R-project.org/package=mixSSG. In the fourth example, we apply a three-component generalized hyperbolic mixture model to Wheat data. The ML estimators of this mixture model have been obtained using the R package MixGHD available at https://cran.r-project.org/package=MixGHD. Finally, we note that if parameter h is very small (less than 0.001, say), then the approximated observed Fisher information matrix may not be invertible.

Value

A list consists of the maximum likelihood estimator, approximated asymptotic standard error, upper, and lower bounds of 100(1α)%100(1-\alpha)\% asymptotic confidence interval for parameters of the finite mixture model.

Author(s)

Mahdi Teimouri

References

K. E. Basford, D. R. Greenway, G. J. McLachlan, and D. Peel, (1997). Standard errors of fitted means under normal mixture, Computational Statistics, 12, 1-17.

S. X. Lee and G. J. McLachlan, (2018). EMMIXcskew: An R package for the fitting of a mixture of canonical fundamental skew t-distributions, Journal of Statistical Software, 83(3), 1-32, doi: 10.18637/jss.v083.i03.

M. Teimouri, (2022). Finite mixture of skewed sub-Gaussian stable distributions, https://arxiv.org/abs/2205.14067.

C. Tortora, R. P. Browne, A. ElSherbiny, B. C. Franczak, and P. D. McNicholas, (2021). Model-based clustering, classification, and discriminant analysis using the generalized hyperbolic distribution: MixGHD R package. Journal of Statistical Software, 98(3), 1-24, doi: 10.18637/jss.v098.i03.

Examples


# Example 1: Approximating the asymptotic standard error and 95 percent confidence interval
#            for the parameters of fitted three-component normal mixture model to iris data.
      Y <- as.matrix( iris[, 1:4] )
colnames(Y) <- NULL
rownames(Y) <- NULL
      G <- 3
 weight <- c( 0.334, 0.300, 0.366         )
    mu1 <- c( 5.0060, 3.428, 1.462, 0.246 )
    mu2 <- c( 5.9150, 2.777, 4.204, 1.298 )
    mu3 <- c( 6.5468, 2.949, 5.482, 1.985 )
 sigma1 <- matrix( c( 0.133, 0.109, 0.019, 0.011, 0.109, 0.154, 0.012, 0.010,
                      0.019, 0.012, 0.028, 0.005, 0.011, 0.010, 0.005, 0.010 ), nrow = 4 , ncol = 4)
 sigma2 <- matrix( c( 0.225, 0.076, 0.146, 0.043, 0.076, 0.080, 0.073, 0.034,
                      0.146, 0.073, 0.166, 0.049, 0.043, 0.034, 0.049, 0.033 ), nrow = 4 , ncol = 4)
 sigma3 <- matrix( c( 0.429, 0.107, 0.334, 0.065, 0.107, 0.115, 0.089, 0.061,
                      0.334, 0.089, 0.364, 0.087, 0.065, 0.061, 0.087, 0.086 ), nrow = 4 , ncol = 4)
     mu <- list(    mu1,    mu2,    mu3 )
  sigma <- list( sigma1, sigma2, sigma3 )
  sigma <- list( sigma1, sigma2, sigma3 )
 lambda <- list( rep(0, 4), rep(0, 4), rep(0, 4) )
   out1 <- sefm( Y, G, weight, model = "restricted", mu, sigma, lambda, family = "constant",
    skewness = "FALSE")
# Example 2: Approximating the asymptotic standard error and 95 percent confidence interval
#            for the parameters of fitted two-component restricted skew t mixture model to
#            AIS data.
      data( AIS )
      Y <- as.matrix( AIS[, 2:3] )
      G <- 2
 weight <- c(  0.5075,  0.4925 )
    mu1 <- c( 19.9827, 17.8882 )
    mu2 <- c( 21.7268,  5.7518 )
 sigma1 <- matrix( c(3.4915, 8.3941, 8.3941, 28.8113 ), nrow = 2, ncol = 2 )
 sigma2 <- matrix( c(2.2979, 0.0622, 0.0622,  0.0120 ), nrow = 2, ncol = 2 )
lambda1 <- ( c( 2.5186, -0.2898 ) )
lambda2 <- ( c( 2.1681,  3.5518 ) )
 theta1 <- c( 68.3088 )
 theta2 <- c(  3.8159 )
     mu <- list(     mu1,     mu2 )
  sigma <- list(  sigma1,  sigma2 )
 lambda <- list( lambda1, lambda2 )
  theta <- list(  theta1,  theta2 )
  param <- c( "nu" )
    PDF <- quote( (nu/2)^(nu/2)*w^(-nu/2 - 1)/gamma(nu/2)*exp( -nu/(w*2) ) )
  tick  <- c( 1, 1 )
   out2 <- sefm( Y, G, weight, model = "restricted", mu, sigma, lambda, family = "igamma",
            skewness = "TRUE", param, theta, tick, h = 0.001, N = 3000, level = 0.05, PDF )
# Example 3: Approximating the asymptotic standard error and 95 percent confidence interval
#            for the parameters of fitted two-component restricted skew sub-Gaussian
#            alpha-stable mixture model to bankruptcy data.
      data( bankruptcy )
      Y <- as.matrix( bankruptcy[, 2:3] );  colnames(Y) <- NULL; rownames(Y) <- NULL
      G <- 2
 weight <- c(  0.553,  0.447 )
    mu1 <- c( -3.649, -0.085 )
    mu2 <- c( 40.635, 19.042 )
 sigma1 <- matrix( c(1427.071, -155.356, -155.356, 180.991 ), nrow = 2, ncol = 2 )
 sigma2 <- matrix( c( 213.938,    9.256,    9.256,  74.639 ), nrow = 2, ncol = 2 )
lambda1 <- c( -41.437, -21.750 )
lambda2 <- c(  -3.666,  -1.964 )
 theta1 <- c( 1.506 )
 theta2 <- c( 1.879 )
     mu <- list(     mu1,     mu2 )
  sigma <- list(  sigma1,  sigma2 )
 lambda <- list( lambda1, lambda2 )
  theta <- list(  theta1,  theta2 )
  param <- c( "alpha" )
  tick  <- c( 1 )
   out3 <- sefm( Y, G, weight, model = "restricted", mu, sigma, lambda, family = "pstable",
            skewness = "TRUE", param, theta, tick, h = 0.01, N = 3000, level = 0.05 )
# Example 4: Approximating the asymptotic standard error and 95 percent confidence interval
#            for the parameters of fitted two-component restricted generalized inverse-Gaussian
#            mixture model to AIS data.
      data( wheat )
      Y <- as.matrix( wheat[, 1:7] ); colnames(Y) <- NULL; rownames(Y) <- NULL
      G <- 3
 weight <- c( 0.325, 0.341, 0.334 )
    mu1 <- c( 18.8329, 16.2235, 0.9001, 6.0826, 3.8170, 1.6604, 6.0260 )
    mu2 <- c( 11.5607, 13.1160, 0.8446, 5.1873, 2.7685, 4.9884, 5.2203 )
    mu3 <- c( 13.8071, 14.0720, 0.8782, 5.5016, 3.1513, 0.6575, 4.9111 )
lambda1 <- diag( c( 0.1308, 0.2566,-0.0243, 0.2625,-0.1259, 3.3111, 0.1057) )
lambda2 <- diag( c( 0.7745, 0.3084, 0.0142, 0.0774, 0.1989,-1.0591,-0.2792) )
lambda3 <- diag( c( 2.0956, 0.9718, 0.0042, 0.2137, 0.2957, 3.9484, 0.6209) )
 theta1 <- c( -3.3387, 4.2822 )
 theta2 <- c( -3.6299, 4.5249 )
 theta3 <- c( -3.9131, 5.8562 )
 sigma1 <- matrix( c(
 1.2936219, 0.5841467,-0.0027135, 0.2395983, 0.1271193, 0.2263583, 0.2105204,
 0.5841467, 0.2952009,-0.0045937, 0.1345133, 0.0392849, 0.0486487, 0.1222547,
-0.0027135,-0.0045937, 0.0003672,-0.0033093, 0.0016788, 0.0056345,-0.0033742,
 0.2395983, 0.1345133,-0.0033093, 0.0781141, 0.0069283,-0.0500718, 0.0747912,
 0.1271193, 0.0392849, 0.0016788, 0.0069283, 0.0266365, 0.0955757, 0.0002497,
 0.2263583, 0.0486487, 0.0056345,-0.0500718, 0.0955757, 1.9202036,-0.0455763,
 0.2105204, 0.1222547,-0.0033742, 0.0747912, 0.0002497,-0.0455763, 0.0893237 ), nrow = 7, ncol = 7 )
 sigma2 <- matrix( c(
 0.9969975, 0.4403820, 0.0144607, 0.1139573, 0.1639597,-0.2216050, 0.0499885,
 0.4403820, 0.2360065, 0.0010769, 0.0817149, 0.0525057,-0.0320012, 0.0606147,
 0.0144607, 0.0010769, 0.0008914,-0.0023864, 0.0049263,-0.0122188,-0.0042375,
 0.1139573, 0.0817149,-0.0023864, 0.0416206, 0.0030268, 0.0490919, 0.0407972,
 0.1639597, 0.0525057, 0.0049263, 0.0030268, 0.0379771,-0.0384626,-0.0095661,
-0.2216050,-0.0320012,-0.0122188, 0.0490919,-0.0384626, 4.0868766, 0.1459766,
 0.0499885, 0.0606147,-0.0042375, 0.0407972,-0.0095661, 0.1459766, 0.0661900 ), nrow = 7, ncol = 7 )
 sigma3 <- matrix( c(
 1.1245716, 0.5527725,-0.0005064, 0.2083688, 0.1190222,-0.4491047, 0.2494994,
 0.5527725, 0.3001219,-0.0036794, 0.1295874, 0.0419470,-0.1926131, 0.1586538,
-0.0005064,-0.0036794, 0.0004159,-0.0034247, 0.0019652,-0.0026687,-0.0044963,
 0.2083688, 0.1295874,-0.0034247, 0.0715283, 0.0055925,-0.0238820, 0.0867129,
 0.1190222, 0.0419470, 0.0019652, 0.0055925, 0.0243991,-0.0715797, 0.0026836,
-0.4491047,-0.1926131,-0.0026687,-0.0238820,-0.0715797, 1.5501246,-0.0048728,
 0.2494994, 0.1586538,-0.0044963, 0.0867129, 0.0026836,-0.0048728, 0.1509183 ), nrow = 7, ncol = 7 )
    mu <- list( mu1, mu2, mu3 )
 sigma <- list( sigma1 , sigma2, sigma3 )
lambda <- list( lambda1, lambda2, lambda3 )
 theta <- list( theta1 , theta2, theta3 )
  tick <- c( 1, 1, 0 )
 param <- c( "a", "b" )
   PDF <- quote( 1/( 2*besselK( b, a ) )*w^(a - 1)*exp( -b/2*(1/w + w) ) )
  out4 <- sefm( Y, G, weight, model = "unrestricted", mu, sigma, lambda, family = "gigaussian",
            skewness = "TRUE", param, theta, tick, h = 0.001, N = 3000, level = 0.05, PDF )


[Package mixbox version 1.2.3 Index]