llnorMix: Likelihood, Parametrization and EM-Steps For 1D Normal... In nor1mix: Normal (1-d) Mixture Models (S3 Classes and Methods)

Description

These functions work with an almost unconstrained parametrization of univariate normal mixtures.

`llnorMix(p, *)`

computes the log likelihood,

`obj <- par2norMix(p)`

maps parameter vector `p` to a `norMix` object `obj`,

`p <- nM2par(obj)`

maps from a `norMix` object `obj` to parameter vector `p`,

where `p` is always a parameter vector in our parametrization.

Partly for didactical reasons, the following functions provide the basic ingredients for the EM algorithm (see also `norMixEM`) to parameter estimation:

`estep.nm(x, obj, p)`

computes 1 E-step for the data `x`, given either a `"norMix"` object `obj` or parameter vector `p`.

`mstep.nm(x, z)`

computes 1 M-step for the data `x` and the probability matrix `z`.

`emstep.nm(x, obj)`

computes 1 E- and 1 M-step for the data `x` and the `"norMix"` object `obj`.

where again, `p` is a parameter vector in our parametrization, `x` is the (univariate) data, and `z` is a n * m `matrix` of (posterior) conditional probabilities, and θ is the full parameter vector of the mixture model.

Usage

 ```1 2 3 4 5 6 7 8``` ```llnorMix(p, x, m = (length(p) + 1)/3) par2norMix(p, name = sprintf("{from %s}", deparse(substitute(p))[1])) nM2par(obj) estep.nm(x, obj, par) mstep.nm(x, z) emstep.nm(x, obj) ```

Arguments

 `p, par` numeric vector: our parametrization of a univariate normal mixture, see details. `x` numeric: the data for which the likelihood is to be computed. `m` integer number of mixture components; this is not to be changed for a given `p`. `name` (for `par2norMix()`:) a name for the `"norMix"` object that is returned. `obj` a `"norMix"` object, see `norMix`. `z` a n * m `matrix` of (posterior) conditional probabilities, z[i,j]= P(x[i] in C_j | theta), where C_j denotes the j-th group (“cluster”).

Details

We use a parametrization of a (finite) univariate normal mixture which is particularly apt for likelihood maximization, namely, one whose parameter space is almost a full R^m, m = 3k-1.

For a k-component mixture, we map to and from a parameter vector θ (`== p` as R-vector) of length 3k-1. For mixture density

sum[j=1..k] pi[j] phi((t - mu[j])/s[j])/s[j],

we logit-transform the pi[j] (for j >= 2) and log-transform the s[j], such that θ is partitioned into

`p[ 1:(k-1)]`:

`p[j]`= logit(pi[j+1]) and pi[1] is given implicitly as pi[1] = 1 - sum[j=2..k] pi[j].

`p[ k:(2k-1)]`:

`p[k-1+ j]`= μ_j, for j=1:k.

`p[2k:(3k-1)]`:

`p[2*k-1+ j]` = log(s[j]), i.e., σ_j^2 = exp(2*p[.+j]).

Value

`llnorMix()` returns a number, namely the log-likelihood.

`par2norMix()` returns `"norMix"` object, see `norMix`.

`nM2par()` returns the parameter vector θ of length 3k-1.

`estep.nm()` returns `z`, the matrix of (conditional) probabilities.

`mstep.nm()` returns the model parameters as a `list` with components `w`, `mu`, and `sigma`, corresponding to the arguments of `norMix()`. (and see the 'Examples' on using `do.call(norMix, *)` with it.)

`emstep.nm()` returns an updated `"norMix"` object.

Author(s)

Martin Maechler

`norMix`, `logLik`. Note that the log likelihood of a `"norMix"` object is directly given by `sum(dnorMix(x, obj, log=TRUE))`.

To fit, using the EM algorithm, rather use `norMixEM()` than the `e.step`, `m.step`, or `em.step` functions.

Note that direct likelihood maximization, i.e., MLE, is typically considerably more efficient than the EM, and typically converges well with our parametrization, see `norMixMLE`.

Examples

 ``` 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``` ```(obj <- MW.nm10) # "the Claw" -- m = 6 components length(pp <- nM2par(obj)) # 17 == (3*6) - 1 par2norMix(pp) ## really the same as the initial 'obj' above ## Log likelihood (of very artificial data): llnorMix(pp, x = seq(-2, 2, length=1000)) set.seed(47)## of more realistic data: x <- rnorMix(1000, obj) llnorMix(pp, x) ## Consistency check : nM2par() and par2norMix() are inverses all.EQ <- function(x,y, tol = 1e-15, ...) all.equal(x,y, tolerance=tol, ...) stopifnot(all.EQ(pp, nM2par(par2norMix(pp))), all.EQ(obj, par2norMix(nM2par(obj)), check.attributes=FALSE), ## Direct computation of log-likelihood: all.EQ(sum(dnorMix(x, obj, log=TRUE)), llnorMix(pp, x)) ) ## E- and M- steps : ------------------------------ rE1 <- estep.nm(x, obj) rE2 <- estep.nm(x, par=pp) # the same as rE1 z <- rE1 str( rM <- mstep.nm(x, z)) (rEM <- emstep.nm(x, obj)) stopifnot(all.EQ(rE1, rE2), all.EQ(rEM, do.call(norMix, c(rM, name="")))) ```

nor1mix documentation built on May 31, 2017, 5:10 a.m.