| llnorMix {nor1mix} | R Documentation |
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.
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)
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 |
name |
(for |
obj |
a |
z |
a n * m |
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]),
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]).
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
sig2, corresponding to the arguments of norMix().
emstep.nm() returns an updated "norMix" object.
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.
(obj <- MW.nm10) # "the Claw" -- m = 6 components
length(pp <- nM2par(obj)) # 17 == (3*6) - 1
par2norMix(pp)
## really the same as the initial \code{obj} (see below)
## 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 :
all.EQ <- function(x,y, tol = 1e-15, ...) all.equal(x,y, tol=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)
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=""))))