Normal (1-d) Mixture Models (S3 Classes and Methods)

Normal (1-d) Mixture Models (S3 Classes and Methods) [R logo]

Help pages for package ‘nor1mix’ version 1.1-1

clus2norMix | Transform Clustering / Grouping to Normal Mixture
dnorMix | Normal Mixture Density
dnorMixL | Normal Mixture Density
dpnorMix | Normal Mixture Density
is.norMix | Mixtures of Univariate Normal Distributions
lines.norMix | Plotting Methods for 'norMix' Objects
llnorMix | Likelihood and Parametrization of 1D Normal Mixtures
m.norMix | Mixtures of Univariate Normal Distributions
MarronWand | Marron-Wand Densities as 'norMix' Objects
mean.norMix | Mixtures of Univariate Normal Distributions
MW.nm1 | Marron-Wand Densities as 'norMix' Objects
MW.nm10 | Marron-Wand Densities as 'norMix' Objects
MW.nm11 | Marron-Wand Densities as 'norMix' Objects
MW.nm12 | Marron-Wand Densities as 'norMix' Objects
MW.nm13 | Marron-Wand Densities as 'norMix' Objects
MW.nm14 | Marron-Wand Densities as 'norMix' Objects
MW.nm15 | Marron-Wand Densities as 'norMix' Objects
MW.nm16 | Marron-Wand Densities as 'norMix' Objects
MW.nm2 | Marron-Wand Densities as 'norMix' Objects
MW.nm2.old | Marron-Wand Densities as 'norMix' Objects
MW.nm3 | Marron-Wand Densities as 'norMix' Objects
MW.nm4 | Marron-Wand Densities as 'norMix' Objects
MW.nm5 | Marron-Wand Densities as 'norMix' Objects
MW.nm6 | Marron-Wand Densities as 'norMix' Objects
MW.nm7 | Marron-Wand Densities as 'norMix' Objects
MW.nm8 | Marron-Wand Densities as 'norMix' Objects
MW.nm9 | Marron-Wand Densities as 'norMix' Objects
nM2par | Likelihood and Parametrization of 1D Normal Mixtures
norMix | Mixtures of Univariate Normal Distributions
par2norMix | Likelihood and Parametrization of 1D Normal Mixtures
plot.norMix | Plotting Methods for 'norMix' Objects
pnorMix | Normal Mixture Cumulative Distribution and Quantiles
print.norMix | Mixtures of Univariate Normal Distributions
qnorMix | Normal Mixture Cumulative Distribution and Quantiles
r.norMix | Ratio of Normal Mixture to Corresponding Normal
rnorMix | Generate 'Normal Mixture' Distributed Random Numbers
sort.norMix | Sort Method for "norMix" Objects
var.norMix | Mixtures of Univariate Normal Distributions


Transform Clustering / Grouping to Normal Mixture

Description

Simple transformation of a clustering or grouping to a normal mixture object (class "norMix", see, norMix.

Usage

clus2norMix(gr, x, name = deparse(sys.call()))

Arguments

gr | a grouping/clustering vector with values in {1,..,K}; possibly a factor.
x | numeric vector of (original) data (of the same length as gr).
name | name for norMix() object; constructed from the call by default.

Value

A call to norMix() with (mu, sig2, w) set to the empirical values of the groups (as defined by split(x,gr).

Note

Via this function, any simple clustering algorithm (such pam) can be used as simple mixture model fitting procedure.

Author(s)

Martin Maechler, Dec. 2007

See Also

norMix; further pam() (or clara()) from package cluster for sensible clusterings.

Examples

x9 <- rnorMix(500, MW.nm9)
library(cluster)
pxc <- pam(x9, k=3)
plot(pxc, which = 2)# silhouette

(nm.p9 <- clus2norMix(pxc$clustering, x9))
plot(nm.p9, p.norm=FALSE)
lines(MW.nm9, col="thistle")
[Package nor1mix version 1.1-1 Index]


Normal Mixture Density

Description

Evaluate the density function of the normal mixture specified as norMix object.

Usage

dnorMix (x, obj,        log = FALSE)

dnorMixL(obj, x = NULL, log = FALSE, xlim = NULL, n = 511)
dpnorMix(x, obj, lower.tail = TRUE)

Arguments

obj | an object of class norMix.
x | numeric vector with abscissa values where to evaluate the density (and probability, for dpnorMix()). For dnorMixL() by default, when NULL, it is constructed from n (and xlim if that is specified).
log | logical indicating log-density values should be returned.
xlim | range of abscissa values, used if x == NULL. By default, xlim is taken as mean plus/minus 3 standard deviations of the normal mixture.
n | number of abscissa values to generate if x is not specified.
lower.tail | logical; if TRUE (default), probabilities are P[X <= x], otherwise, P[X > x].

Value

dnorMix(x) returns the numeric vector of density values f(x), logged if log is TRUE.
dnorMixL() returns a list with components

x | the abscissa values.
y | the density values f(x) as for dnorMix().
d | the density values f(x) as for dnorMix().
p | the probability values F(x) as for pnorMix().

See Also

rnorMix for random number generation, and norMix for the construction and further methods, particularly plot.norMix which makes use dnorMix.

Examples

 ff <- dnorMixL(MW.nm7)
 str(ff)
 plot(ff, type = "h", ylim = c(0,1)) # rather use plot(ff, ...)

 x <- seq(-4,5, length=501)
 dp <- dpnorMix(x, MW.nm7)
 lines(x, dp$d, col = "tomato", lwd=3)
 lines(x, dp$p, col = 3, lwd=2)# does not fit y-wise
 stopifnot(all.equal(dp$d, dnorMix(x, MW.nm7), tol=1e-12),
           all.equal(dp$p, pnorMix(x, MW.nm7), tol=1e-12))


Likelihood and Parametrization of 1D Normal Mixtures

Description

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

llnorMix(p, *) computes the log likelihood, where as
obj <- par2norMix(p) and
p <- nM2par(obj)
map to and from norMix objects obj and parameter vector p in our parametrization.

Usage

llnorMix(p, x, m = (length(p) + 1)/3)

par2norMix(p, name = sprintf("{from %s}", deparse(substitute(p))[1]))
nM2par(obj)

Arguments

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

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 theta (== 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 theta 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 theta of length 3k-1.

Author(s)

Martin Maechler

See Also

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

Examples

(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 (see below)

## Log likelihood (of very artificial data):
llnorMix(pp, x = seq(-2, 2, length=1000))
## of more realistic data:
x <- rnorMix(1000, obj)
llnorMix(pp, x)

## Consistency check :
stopifnot(all.equal(pp, nM2par(par2norMix(pp)), tol= 1e-15),
          all.equal(obj, par2norMix(nM2par(obj)),
                    check.attributes=FALSE, tol=1e-15),
          ## Direct computation of log-likelihood:
          all.equal(sum(dnorMix(x, obj, log=TRUE)),
                    llnorMix(pp, x), tol= 1e-15)  )


Marron-Wand Densities as 'norMix' Objects

Description

The fifteen density examples used in Marron and Wand (1992)'s simulation study have been used in quite a few subsequent studies, can all be written as normal mixtures and are provided here for convenience and didactical examples of normal mixtures. Number 16 has been added by Jansen et al.

Usage

 MW.nm1  # Gaussian
 MW.nm2  # Skewed
 MW.nm2.old # Skewed(old)
 MW.nm3  # Str Skew
 MW.nm4  # Kurtotic
 MW.nm5  # Outlier
 MW.nm6  # Bimodal
 MW.nm7  # Separated (bimodal)
 MW.nm8  # Asymmetric Bimodal
 MW.nm9  # Trimodal
 MW.nm10 # Claw
 MW.nm11 # Double Claw
 MW.nm12 # Asymmetric Claw
 MW.nm13 # Asymm. Double Claw
 MW.nm14 # Smooth   Comb
 MW.nm15 # Discrete Comb
 MW.nm16 # Distant Bimodal

Author(s)

Martin Maechler

Source

They are translated from Steve Marron's Matlab code at http://www.stat.unc.edu/postscript/papers/marron/parameters/nmpar.m, however for number 2, the Matlab code had MW.nm2.old; and I've defined MW.nm2 as from the Annals paper; see also the last example below.

References

Marron, S. and Wand, M. (1992) Exact Mean Integrated Squared Error; Annals of Statistcs 20, 712–736.

For number 16,
Janssen, Marron, Verb..., Sarle (1995) ....

Examples

MW.nm10
plot(MW.nm14)

## These are defined as norMix() calls in  ../R/zMarrWand-dens.R
nms <- ls(pat="^MW.nm", "package:nor1mix")
nms <- nms[order(as.numeric(substring(nms,6)))]
for(n in nms) {
   cat("\n",n,":\n"); print(get(n, "package:nor1mix"))
}

## Plot all of them:
op <- par(mfrow=c(4,4), mgp = c(1.2, 0.5, 0), tcl = -0.2,
          mar = .1 + c(2,2,2,1), oma = c(0,0,3,0))
for(n in nms[-17]) plot(get(n, "package:nor1mix"))
mtext("The Marron-Wand Densities", outer= TRUE, font= 2, cex= 1.6)

## and their Q-Q-plots (not really fast):
prob <- ppoints(N <- 100)
for(n in nms[-17])
   qqnorm(qnorMix(prob, get(n, "package:nor1mix")), main = n)
mtext("QQ-plots of Marron-Wand Densities", outer = TRUE,
      font = 2, cex = 1.6)
par(op)

## "object" overview:
cbind(sapply(nms, function(n) { o <- get(n)
      sprintf("%-18s: K =%2d; rng = [%3.1f, %2.1f]",
              attr(o, "name"), nrow(o),
              min(o[,"mu"] - 3*sqrt(o[,"sig2"])),
              max(o[,"mu"] + 3*sqrt(o[,"sig2"])) )
             }))

## Note that Marron-Wand (1992), p.720  give #2 as
MW.nm2
## the parameters of which at first look quite different from
MW.nm2.old
## which has been the definition in the above "Source" Matlab code.
## It's easy to see that mu_{nm2} = -.3 + 1.2   * mu_{paper},
## and correspondigly,   s2_{nm2} =       1.2^2 * s2_{paper}
## such that they are "identical" apart from scale and location:
op. <- par(mfrow=2:1, mgp= c(1.2,0.5,0), tcl= -0.2, mar=.1+c(2,2,2,1))
plot(MW.nm2)
plot(MW.nm2.old)
par(op.)


Mixtures of Univariate Normal Distributions

Description

Objects of class norMix represent finite mixtures of (univariate) normal (aka Gaussian) distributions. Methods for construction, printing, plotting, and basic computations are provided.

Usage

norMix(mu, sig2 = rep(1,m), w = NULL, name = NULL, long.name = FALSE)

is.norMix(obj)
m.norMix(obj)
var.norMix(x, ...)
## S3 method for class 'norMix':
mean(x, ...)
## S3 method for class 'norMix':
print(x, ...)

Arguments

mu | numeric vector of length K, say, specifying the means μ of the K normal components.
sig2 | numeric vector of length K, specifying the variances σ^2 of the K normal components.
w | numeric vector of length K, specifying the mixture proportions p[j] of the normal components, j = 1,...,K. Defaults to equal proportions
name | optional name tag of the result (used for printing).
long.name | logical indicating if the name attribute should use punctuation and hence be slightly larger than by default.
obj,x | an object of class norMix.
... | further arguments passed to methods.

Details

The (one dimensional) normal mixtures, R objects of class "norMix", are constructed by norMix and tested for by is.norMix. m.norMix() returns the number of mixture components; the mean() method (for class "norMix" returns the mu vector of means and var.norMix() (not a method, call the function explicitly!) the sig2 vector of variances.

For further methods see below.

Value

norMix returns objects of class "norMix" which are currently implemented as 3-column matrix with column names mu, sig2, and w, and further attributes. The user should rarely need to access the underlying structure directly.

Note

For estimation of the parameters of a normal mixture distribution, I recommend using other R packages, such as flexmix.

Author(s)

Martin Maechler

See Also

dnorMix for the density, pnorMix for the cumulative distribution and the quantile function (qnorMix), and rnorMix for random numbers and plot.norMix, the plot method.

MarronWand has the Marron-Wand densities as normal mixtures.

Examples

ex <- norMix(mu = c(1,2,5))# s^2 = 1, equal proportions
ex
plot(ex)# looks like a mixture of only 2

plot(ex, log = "y")# maybe "revealing"


Plotting Methods for 'norMix' Objects

Description

The plot and lines methods for norMix objects draw the normal mixture density, optionally additonally with a fitted normal density.

Usage

## S3 method for class 'norMix':
plot(x, type = "l", n = 511, xout = NULL, xlim = NULL,
      xlab = "x", ylab = "f(x)", main = attr(x, "name"), lwd = 1.4,
      p.norm = TRUE, p.h0 = TRUE, p.comp = FALSE,
      parNorm = list(col = 2, lty = 2, lwd = 0.4),
      parH0 = list(col = 3, lty = 3, lwd = 0.4),
      parComp = list(col= "blue3", lty = 3, lwd = 0.4), ...)

## S3 method for class 'norMix':
lines(x, type = "l", n = 511, xout = NULL,
      lwd = 1.4, p.norm = FALSE, parNorm = list(col = 2, lty = 2, lwd = 0.4),
      ...)

Arguments

x | object of class norMix.
type | character denoting type of plot, see, e.g. lines.
n | number of points to generate if xout is unspecified.
xout | numeric or NULL giving the abscissae at which to draw the density.
xlim | range of x values to use; particularly important if xout is not specified where xlim is passed to dnorMix and gets a smart default if unspecified.
xlab,ylab | labels for the x and y axis with defaults.
main | main title of plot, defaulting to the norMix name.
lwd | line width for plotting with a non-standard default.
p.norm | logical indicating if the normal density with the same mean and variance should be drawn as well.
p.h0 | logical indicating if the line y = 0 should be drawn.
p.comp | logical indicating if the Gaussian components should also be drawn individually.
parNorm | graphical parameters for drawing the normal density if p.norm is true.
parH0 | graphical parameters for drawing the line y = 0 if p.h0 is true.
parComp | graphical parameters for drawing the single components if p.comp is true.
... | further arguments passed to and from methods.

Author(s)

Martin Maechler

See Also

norMix for the construction and further methods, particularly dnorMix which is used here.

Examples

plot(norMix(m=c(0,3),s=c(4,1)))

plot(MW.nm4, p.norm=FALSE, p.comp = TRUE)
## Further examples in  ?norMix and  ?rnorMix


Normal Mixture Cumulative Distribution and Quantiles

Description

Compute cumulative probabilities or quantiles (the inverse) for a normal mixture specified as norMix object.

Usage

pnorMix(q, obj, lower.tail = TRUE, log.p = FALSE)

qnorMix(p, obj, lower.tail = TRUE, log.p = FALSE,
        tol = .Machine$double.eps^0.25, maxiter = 1000, traceRootsearch = 0,
        method = c("interpQspline", "interpspline", "eachRoot", "root2"),
        l.interp = 20)

Arguments

obj | an object of class norMix.
p | numeric vector of probabilities. Note that for all methods but "eachRoot", qnorMix(p, *) works with the full vector p, typically using (inverse) interpolation approaches; consequently the result is very slightly dependent on p as a whole.
q | numeric vector of quantiles
lower.tail | logical; if TRUE (default), probabilities are P[X <= x], otherwise, P[X > x].
log.p | logical; if TRUE, probabilities p are given as log(p).
tol, maxiter | tolerance and maximal number of iterations for the root search algorithm, see method below and uniroot.
traceRootsearch | logical or integer in {0,1,2,3}, determining the amount of information printed during root search.
method | a string specifying which algorithm is used for the “root search”. .......
l.interp | positive integer for method = "interQpspline" or "interpspline", determining the number of values in each “mu-interval”.

Details

Whereas the distribution function pnorMix is the trivial sum of weighted normal probabilities (pnorm), its inverse, qnorMix is computed numerically: For each p we search for q such that pnorMix(obj, q) == p, i.e., f(q) = 0 for f(q) := pnorMix(obj, q) - p. This is a root finding problem which can be solved by uniroot(f, lower,upper,*). If length(p) <= 2 or method = "eachRoot", this happens one for one for the sorted p's. Otherwise, we start by doing this for the outermost non-trivial (0 < p < 1) values of p.

For method = "interQpspline" or "interpspline", we now compute p. <- pnorMix(q., obj) for values q. which are a grid of length l.interp in each interval [q_j,q_{j+1}] and then use montone inverse interpolation (splinefun(q., p., method="monoH.FC")) plus a few (maximally maxiter, typically one!) Newton steps. The default, "interQpspline", additionally logit-transforms the p. values to make the interpolation more linear. This method is faster, particularly for large length(p).

Value

a numeric vector of the same length as p or q, respectively.

Author(s)

Very first version (for length-1 p,q) by Erik Jørgensen Erik.Jorgensen@agrsci.dk.

See Also

dnorMix for the density function.

Examples

MW.nm3 # the "strange skew" one
plot(MW.nm3)
## now the cumlative :
x <- seq(-4,4, length=1001)
plot(x, pnorMix(x, MW.nm3), type="l", col=2)
## and some of its inverse :
pp <- seq(.1, .9, by=.1)
plot(qnorMix(pp, MW.nm3), pp)

## The "true" median of a normal mixture:
median.norMix <- function(x) qnorMix(1/2, x)
median.norMix(MW.nm3) ## -2.32


Ratio of Normal Mixture to Corresponding Normal

Description

Compute r(x) = f(x)/ f0(x) where f() is a normal mixture density and f0 the normal density with the same mean and variance as f.

Usage

r.norMix(obj, x = NULL, xlim = NULL, n = 511, xy.return = TRUE)

Arguments

obj | an object of class norMix.
x | numeric vector with abscissa values where to evaluate the density. Default is constructed from n (and xlim if specified).
xlim | range of abscissa values, used if x == NULL. By default, xlim taken as mean plus/minus 3 standard deviations of the normal mixture.
n | number of abscissa values to generate if x is not specified.
xy.return | logical indicating if the result should be a list or just a numeric vector, see below.

Value

It depends on xy.return. If it's false, a numeric vector of the same length as x, if true (as per default), a list that can be plotted, with components

x | abscissa values corresponding to argument x.
y | corresponding values r(x).
f0 | values of the moment matching normal density f0(x).

Note

The ratio function is used in certain semi-parametric density estimation methods (and theory).

Examples

  d3 <- norMix(m = 5*(0:2), w = c(0.6, 0.3, 0.1))
  plot(d3)
  rd3 <- r.norMix(d3)
  str(rd3)
  stopifnot(rd3 $ y  == r.norMix(d3, xy.ret = FALSE))
  par(new = TRUE)
  plot(rd3, type = "l", col = 3, axes = FALSE, xlab = "", ylab="")
  axis(4, col.axis=3)


Generate 'Normal Mixture' Distributed Random Numbers

Description

Generate n random numbers, distributed according to a normal mixture.

Usage

rnorMix(n, obj)

Arguments

n | the number of random numbers desired.
obj | an object of class norMix.

Details

For a mixture of m, i.e., m.norMix(obj), components, generate the number in each component as multinomial, and then use rnorm for each.

Value

numeric vector of length n.

See Also

dnorMix for the density, and norMix for the construction and further methods.

Examples

x <- rnorMix(5000, MW.nm10)
hist(x)# you don't see the claw
plot(density(x), ylim = c(0,0.6),
     main = "Estim. and true 'MW.nm10' density")
lines(MW.nm10, col = "orange")


Sort Method for "norMix" Objects

sort.norMix(nor1mix) R Documentation

Sort Method for 'norMix' Objects

Description

Sorting a "norMix" object (see norMix), sorts along the mu values; i.e., for the default decreasing = FALSE the resulting x[,"mu"] are sorted from left to right.

Usage

## S3 method for class 'norMix':
sort(x, decreasing = FALSE, ...)

Arguments

x | an object of class "norMix".
decreasing | logicial indicating if sorting should be up or down.
... | further arguments passed to sort(x[,"mu"],*).

Value

a "norMix" object like x.

Examples

sort(MW.nm9)
stopifnot(identical(MW.nm2, sort(MW.nm2)))