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, trafo = c("clr1", "logit")) par2norMix(p, trafo = c("clr1", "logit"), name = ) nM2par(obj, trafo = c("clr1", "logit")) 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 |
trafo |
|
name |
(for |
obj |
a |
z |
a n * m |
We use a parametrization of a (finite) m-component univariate normal mixture which is particularly apt for likelihood maximization, namely, one whose parameter space is almost a full R^M, M = 3m-1.
For an m-component mixture,
we map to and from a parameter vector θ (== p
as R-vector)
of length 3m-1. For mixture density
sum[j=1..m] pi[j] phi((t - mu[j])/s[j])/s[j],
we transform the pi[j] (for j =
1,…,m) via the transform specified by trafo
(see below),
and log-transform the s[j]. Consequently, θ is
partitioned into
p[ 1:(m-1)]
: For
trafo = "logit"
:p[j]
= logit(pi[j+1]) and
pi[1] is given implicitly as
pi[1] = 1 - sum[j=2..m] pi[j].
trafo = "clr1"
:(centered log
ratio, omitting 1st element):
Set l_j := ln(pi[j]) for j = 1, …, m, and
p[j]
= l_{j+1} - mean_j' l_j'
for j = 1, …, m-1.
p[ m:(2m-1)]
: p[m-1+ j]
= μ_j, for j=1:m.
p[2m:(3m-1)]
: p[2*m-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 3m-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.
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
.
(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=""))))