Nothing
##
## Copyright (C) 1997-2023 Martin Maechler
##
## This program is free software; you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published by
## the Free Software Foundation; either version 2 of the License, or
## (at your option) any later version.
##
## A copy of the GNU General Public License is available at
## http://www.r-project.org/Licenses/
####---- Normal Mixtures "norMix" -------
####---- ~~~~~~~~~~~~~~~ ###### -------
#### Object-oriented S/R - functions for dealing with 1D normal mixtures.
#### -------------------------------------------------------------------------
#### Author: Martin Maechler, 20 Mar 1997
#### -------------------------------------------------------------------------
if(getRversion() < "2.15") paste0 <- function(...) paste(..., sep = '')
norMix <- function(mu, sig2 = rep(1, m), sigma = rep(1, m), w = NULL,
name = NULL, long.name = FALSE)
{
## Purpose: Constructor for 'norMix' (normal mixture) objects
## -------------------------------------------------------------------------
## Arguments: mu: vector of means; sig2: vector of variances sigma^2
## w : vector of weights (adding to 1) -- default: equal
## name : name attribute; constructed from (mu,sig2) by default
## long.name : logical used for default \code{name} construction
## -------------------------------------------------------------------------
## Author: Martin Maechler, Date: 20 Mar 97, 14:58
if(!is.numeric(mu)) stop("'mu' must be numeric!")
m <- length(mu)
if(!missing(sig2)) {
if(!missing(sigma))
stop("you must not specify both 'sig2' and 'sigma'; the latter is preferred now")
sigma <- sqrt(sig2)
.Deprecated(msg =
"The use of 'sig2' is deprecated; do specify 'sigma' (= sqrt(sig2)) instead")
}
if(length(sigma) == 1) sigma <- rep.int(sigma, m)
if(length(sigma) != m || !is.numeric(sigma)|| any(sigma <=0))
stop("'sigma' must be > 0 with same length as 'mu'")
if(is.null(w))
w <- rep.int(1/m, m)
else {
if(length(w) != m || !is.numeric(w) || any(w<0))
stop("'w' must be >= 0 with same length as 'mu'")
s <- sum(w)
if(abs(s-1) > 10*.Machine$double.eps) w <- w / s
}
if(is.null(name)) {
sformat <- function(v) sapply(v, format, digits=1)
pPar <- function(pp) {
pp <- if(m >= 10) c(sformat(pp[1:9]), "....") else sformat(pp)
if(long.name)
paste0("(", paste(pp, collapse=","), ")")
else
paste(pp, collapse="")
}
name <- paste0("NM",format(m),".", pPar(mu), "_", pPar(sigma))
}
structure(name = name, class = "norMix",
.Data = cbind(mu = mu, sigma = sigma, w = w))
}
`[.norMix` <- function (x, i, j, drop = TRUE) {
if(!missing(j) && "sig2" %in% j) { ## back-compatibility hack
.Deprecated(msg =
"The use of 'sig2' is deprecated; do specify 'sigma' (= sqrt(sig2)) instead")
if(length(j) == 1L) ## return sig2 = sigma^2:
if(missing(i)) x[,"sigma", drop=drop]^2 else x[i, "sigma", drop=drop]^2
else
stop("subsetting with 'sig2' together with other columns is no longer allowed; use 'sigma'")
}
else if(missing(j) && (nargs() == 3 || !drop)) { ## a ("subset mixture") norMix object:
r <- NextMethod("[")
if(!is.matrix(r)) ## e.g. for x[1, ]
r
else {
r[,"w"] <- r[,"w"] / sum(r[,"w"]) # renormalize
structure(r, class = class(x), name =
paste0(attr(x,"name"), "[", deparse(substitute(i), 20)[1L], ",]"))
}
}
else NextMethod("[")
}
is.norMix <- function(obj)
{
## Purpose: is 'obj' a "norMix", i.e. Normal Mixture object ?
## Author: Martin Maechler, Date: 20 Mar 97, 10:38
inherits(obj, "norMix") &&
(!is.null(w <- obj[,"w"])) &&
is.numeric(w) && all(w >= 0) && abs(sum(w)-1) < 1000*.Machine$double.eps
}
m.norMix <- function(obj) nrow(obj) ## Number of components of normal mixture
mean.norMix <- function(x, ...)
{
## Purpose: Return "true mean", i.e., the expectation of a normal mixture.
if(!is.norMix(x)) stop("'x' must be a 'Normal Mixture' object!")
x <- unclass(x)
drop(x[,"w"] %*% x[,"mu"])
}
var.norMix <- function(x, ...)
{
## Purpose: 'true' Variance, i.e. E[(X- E[X])^2] for X ~ normal mixture.
if(!is.norMix(x)) stop("'x' must be a 'Normal Mixture' object!")
x <- unclass(x)
w <- x[,"w"]
mj <- x[,"mu"]
mu <- drop(w %*% mj)
drop(w %*% (x[,"sigma"]^2 + (mj - mu)^2))
}
print.norMix <- function(x, ...)
{
## Purpose: print method for "norMix" objects (normal mixtures)
## -------------------------------------------------------------------------
## Author: Martin Maechler, Date: 20 Mar 97, 10:02
ox <- x
has.nam <- !is.null(nam <- attr(x,"name"))
cat("'Normal Mixture' object",
if(has.nam) paste("\t ``", nam, "''", sep=''), "\n")
if(has.nam) attr(x, "name") <- NULL
cl <- class(x); cl <- cl[ cl != "norMix"] #- the remaining classes
class(x) <- if(length(cl) > 0) cl ## else NULL
NextMethod("print", ...)
invisible(ox)
}
sort.norMix <- function(x, decreasing = FALSE, ...) {
## sort according to 'mu' (and ensure attributes as "name" are not changed):
x[] <- x[sort.list(x[,"mu"], decreasing = decreasing, ...) , ]
x
}
##' Return a call, e.g., something deriv() can use
norMix2call <- function(obj, oneArg = TRUE) {
w <- obj[,"w"]; mu <- obj[,"mu"]; sd <- obj[,"sigma"]
m <- length(w) #-- number of components
if(oneArg) {
ex <- substitute(W * dnorm((x - M)/S),
list(W = unname(w[1]), M = unname(mu[1]), S = unname(sd[1])))
for(j in seq_len(m)[-1L]) # j in 2:m, but empty if(m == 1)
ex <- substitute(EX + W * dnorm((x - M)/S),
list(EX = ex, W = w[j], M = mu[j], S = sd[j]))
} else {
## the first term
ex <- substitute(W * dnorm(x, mean = M, sd = S),
list(W = unname(w[1]), M = unname(mu[1]), S = unname(sd[1])))
for(j in seq_len(m)[-1L]) # j in 2:m, but empty if(m == 1)
ex <- substitute(EX + W * dnorm(x, mean = M, sd = S),
list(EX = ex, W = w[j], M = mu[j], S = sd[j]))
}
## return
ex
}
as.expression.norMix <- function(x, oneArg = TRUE, ...) as.expression(norMix2call(x, oneArg), ...)
as.function.norMix <- function(x, oneArg = TRUE, envir = parent.frame(), ...)
`body<-`(function(x) {}, envir=envir, value = norMix2call(x, oneArg))
dnorMix <- function(x, obj, log = FALSE)
{
## Purpose: density evaluation for "norMix" objects (normal mixtures)
## -------------------------------------------------------------------------
## Arguments: x: numeric; obj: Normal Mixture object
## -------------------------------------------------------------------------
## Author: Martin Maechler, Date: 20 Mar 97, 10:14
if(!is.norMix(obj)) {
## Old version had (obj, x, ..):
if(is.norMix(x)) { ## swap the first two arguments
tmp <- x ; x <- obj; obj <- tmp
.Deprecated(msg = "Deprecated use of dnorMix(obj, x, ..);
Either use dnorMixL(), or the new argument order (x, obj, ...) and
note that dnorMix() returns a numeric vector (not a list).")
}
else stop("'obj' must be a 'Normal Mixture' object!")
}
w <- obj[,"w"]; mu <- obj[,"mu"]; sd <- obj[,"sigma"]
m <- length(w) #-- number of components
if(m == 1)
return(dnorm(x, mean = mu[1], sd = sd[1], log = log))
## else
y <- numeric(length(x))
for(j in 1:m)
y <- y + w[j] * dnorm(x, mean = mu[j], sd = sd[j])
if(log) log(y) else y
}
dnorMixL <- function(obj, x = NULL, log = FALSE, xlim = NULL, n = 511)
{
## Purpose: density evaluation for "norMix" objects (normal mixtures)
## -------------------------------------------------------------------------
## Arguments: obj: Normal Mixture object; x:
## -------------------------------------------------------------------------
## Author: Martin Maechler, Date: 20 Mar 97, 10:14
if(!is.norMix(obj))
stop("'obj' must be a 'Normal Mixture' object!")
if(is.null(x)) {
if(is.null(xlim)) ##-- construct "reasonable" abscissa values
xlim <- mean.norMix(obj) + c(-3,3)*sqrt(var.norMix(obj))
x <- seq.int(xlim[1], xlim[2], length.out = n)
}
list(x = x, y = dnorMix(x, obj, log=log))
}
rnorMix <- function(n, obj)
{
## Purpose: Generate random numbers according to "norMix" object `obj'
## -------------------------------------------------------------------------
## Author: Martin Maechler, Date: 27 Jun 2002, 16:03
mu <- obj[,"mu"]
sd <- obj[,"sigma"]
if(n == 1) {
j <- sample(length(mu), size = 1, prob = obj[,"w"])
rnorm(1, mean = mu[j], sd = sd[j])
} else {
nj <- as.vector(rmultinom(n=1, size = n, prob = obj[,"w"]))
sample(unlist(lapply(seq_along(nj), function(j)
rnorm(nj[j], mean = mu[j], sd = sd[j]))))
}
}
## From: Erik Jørgensen <Erik.Jorgensen@agrsci.dk>
## Date: Thu, 13 Nov 2003 02:06:27 +0100
##
## ....... Please, feel free to use them.
##
## Erik Jørgensen
## Danish Institute of Agricultural Sciences
pnorMix <- function(q, obj, lower.tail = TRUE, log.p = FALSE)
{
if (!is.norMix(obj)) {
## Old version had (obj, q):
if(is.norMix(q)) { ## swap the first two arguments
tmp <- q ; q <- obj; obj <- tmp
.Deprecated(msg =
"Deprecated use of pnorMix(obj, q, ..); NEW argument order is (q, obj, ...)")
}
else stop("'obj' must be a 'Normal Mixture' object!")
}
sd <- obj[,"sigma"]
## if log.p just log(.) at the end [to be more accurate, need much more..]
cc <- if(log.p) function(m) log(c(m)) else c
## q can be a vector: -> outer
cc(pnorm(sweep(outer(q, obj[,"mu"], "-"), 2, sd, "/"),
lower.tail= lower.tail) %*% obj[, "w"])
}
dpnorMix <- function(x, obj, lower.tail = TRUE)
{
## Purpose: compute dnorMix() and pnorMix() simultaneously
## ----------------------------------------------------------------------
## Arguments: x: numeric; obj: 'norMix' object
## ----------------------------------------------------------------------
## Author: Martin Maechler, Date: 3 Jan 2008
stopifnot (is.norMix(obj))
obj <- unclass(obj)
w <- obj[,"w"]; mu <- obj[,"mu"]; sd <- obj[,"sigma"]
## This looks smarter, but really is slower :
## z <- sweep(outer(x, obj[,"mu"], "-"), 2, obj[,"sigma"], "/")
## list(d = c(dnorm(z) %*% w),
## p = c(pnorm(z, lower.tail= lower.tail) %*% w))
m <- length(w) #-- number of components
d <- p <- numeric(length(x))
for(j in 1:m) {
d <- d + w[j] * dnorm(x, mean= mu[j], sd= sd[j])
p <- p + w[j] * pnorm(x, mean= mu[j], sd= sd[j], lower.tail=lower.tail)
}
list(d = d, p = p)
}
##' <description>
##'
##' <details>
##' @title Sub Sequence, regularly from 1:m
##' @param m integer >= 0 (typically >= M)
##' @param M integer >= 2
##' @return an approximately arithmetic sub-sequence of 1:m of length <= M
##' @author Martin Maechler
sub_seq <- function(m, M) {
stopifnot((M1 <- as.integer(M-1L)) >= 1L, length(m <- as.integer(m)) == 1L)
if(m <= 2) seq_len(m)
else pmax(1L, pmin(m, unique(as.integer(
ceiling((m*floor(seq_len(m)*M1/m))/M1)))))
}
qnorMix <-
function(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 = pmax(1, pmin(20, 1000 / m)), n.mu.interp = 100)
## NOTE: keep defaults consistent with 'uniroot':
{
if (!is.norMix(obj)) {
## Old version had (obj, p):
if(is.norMix(p)) { ## swap the first two arguments
tmp <- p ; p <- obj; obj <- tmp
.Deprecated(msg =
"Deprecated use of qnorMix(obj, p, ..); NEW argument order is (p, obj, ...)")
}
else stop("'obj' must be a 'Normal Mixture' object!")
}
mu <- obj[, "mu"]
sd <- obj[, "sigma"]
m <- m.norMix(obj)
if(m == 1) # one component
return(qnorm(p, mu, sd, lower.tail=lower.tail, log.p=log.p))
## else
S <- if(lower.tail) 1 else -1
## vectorize in 'p'
r <- p
## Solve the left/right extremes p \in {0 , 1}
left <- if(log.p) rep(FALSE, length(p)) else p <= 0
right <- p >= if(log.p) 0 else 1
r[left] <- -Inf*S ; r[right] <- Inf*S
imid <- which(mid <- !left & !right) # 0 < p < 1
if(length(imid)) {
f.make <- function(p.i) {
if(traceRootsearch >= 3)
function(l) {
p <- pnorMix(l, obj, lower.tail=lower.tail, log.p=log.p)
cat(sprintf("p(%-19.16g) = %-19.16g\n", l, p))
p - p.i
}
else
function(l) pnorMix(l, obj, lower.tail=lower.tail,
log.p=log.p) - p.i
}
outRange <- function(p.i)
range(qnorm(p.i, mu, sd, lower.tail=lower.tail, log.p=log.p))
## sort p[] increasingly for easier root finding start:
p <- sort(p[mid], index.return = TRUE)
ip <- imid[p$ix]
pp <- p$x
hasDup <- any(iDup <- duplicated(pp))
if(hasDup) {
isUniq <- !iDup
## want *strictly* increasing sometimes; save CPU anyway
pp <- pp[isUniq]
i.pp <- cumsum(isUniq) ## pp[i.pp] |--> original pp[]
}
np <- length(pp)
rr <- pp # rr will contain = q..mix(pp, *)
missMeth <- missing(method)
method <- {
if(np <= 2) "eachRoot" ## in any case
else if(missMeth && m >= 100) "root2" else match.arg(method)
}
if(method == "eachRoot") { ## root finding from left to right ...
for(i in seq_along(pp)) {
ff <- f.make(pp[i])
rq <- outRange(pp[i])
## since pp[] is increasing, we can start from last 'root':
if(i > 1 && rq[1] < root)
rq[1] <- root
root <- safeUroot(ff, Sig = S, interval = rq, tol=tol, maxiter=maxiter,
trace = traceRootsearch)$root
rr[i] <- root
}
}
else { ## other 'method's => np > 2
rr[1] <- safeUroot(f.make(pp[1]), Sig = S, interval = outRange(pp[1]),
tol=tol, maxiter=maxiter,
trace = traceRootsearch)$root
rr[np] <- safeUroot(f.make(pp[np]), Sig = S, interval = outRange(pp[np]),
tol=tol, maxiter=maxiter,
trace = traceRootsearch)$root
ni <- length(iDone <- as.integer(c(1,np)))
if(any(method == c("interpQspline", "interpspline"))) {
## reverse interpolate, using relatively fast pnorMix()!
pp. <- pp[-iDone]
## those mu's that are inside our range:
rXtr <- rr[if(lower.tail) c(1L,np) else c(np,1L)]
mu. <- unique(sort(mu[rXtr[1] < mu & mu < rXtr[2]]))
## l.interp values between each mu
stopifnot(l.interp >= 1, n.mu.interp > 1)
## large m (== length(mu)) would give large k, and below,
## pnorMix() uses outer() --> a matrix of size m * k * l.interp
k <- length(qs <- c(rXtr[1L],
mu.[sub_seq(length(mu.), n.mu.interp)],
rXtr[2L]))
qs. <- qs[-k]
dq <- qs[-1] - qs. # == delta(qs)
qi <- c(t(dq %*% t(seq_len(l.interp)/l.interp) + qs.))
stopifnot(!is.unsorted(qi)) ## << FIXME remove if never triggering
ppi <- pnorMix(qi, obj, lower.tail=lower.tail, log.p=log.p)
## in an extreme case, pnorMix() is horizontal; hence
## qnorMix() has practically a discontinuity there.
## In that case, splinefun() completely "fails";
## we do need *monotone* (spline) interpolation:
mySfun <- function(x, y) # {specifying 'ties' also avoids warning}
splinefun(x,y, ties=mean.default, method="monoH.FC")
if(method == "interpspline")
qpp <- mySfun(ppi, qi)(pp.) ## is very fast
else { ## "interpQspline"
## logit() transform the P's --> interpolation is more linear
muT <- drop(obj[, "w"] %*% mu)
qp. <- qlogis(pp., muT, log.p=log.p)
qpp <- mySfun(qlogis(ppi, muT, log.p=log.p), qi)(qp.)
}
if(log.p)
warning("Newton steps for 'log.p = TRUE' not yet implemented") ## TODO!
else {
## now end with a few Newton steps
for(k in 1:maxiter) {
dp <- dpnorMix(qpp, obj, lower.tail=lower.tail)
del.p <- dp$p - pp.
## FIXME?: del.p may suffer from considerable cancellation
relE.f <- abs(del.p)
n0 <- relE.f > 0 & pp. > 0
relE.f[n0] <- (relE.f/pp.)[n0]
ii. <- dp$d > 0 ## & relE.f > tol
if(!any(ii.)) {
relErr <- mean(relE.f)
break # not converged though
}
## del.q := Delta(q) = F(q) / f(q) or 0 if f(q)=0
del.q <- numeric(length(pp.))
del.q[ii.] <- S*(del.p/dp$d)[ii.]
## only modify qpp[] where Newton step is ok:
## e.g. resulting qpp must remain increasing
while(length(iF <- which(S*diff(qNew <- qpp - del.q) <= 0))) {
iF <- c(iF,iF+1L)
del.q[iF] <- del.q[iF] / 2
if(traceRootsearch) cat(",")
}
qpp[ii.] <- qNew[ii.]
relErr <- sum(abs(del.q[ii.])) / sum(abs(qpp[ii.]))
if(traceRootsearch) {
cat(k,": relE =", formatC(relErr), sep='')
if(traceRootsearch >= 2) {
cat(" |\n")
if(traceRootsearch == 2 || length(qpp) <= 10)
print(summary(del.q / abs(qpp)))
else print(del.q / abs(qpp))
}
else cat("\n")
}
if(relErr < tol) break
}
if(relErr >= tol)
warning("Newton iterations have not converged")
}
rr[-iDone] <- qpp
}
else ## method == "root2"
while(ni < np) {
## not "done"; ni == length(iDone)
oi <- iDone
i.1 <- oi[-ni]
i.2 <- oi[-1]
l.new <- i.2 > i.1 + 1L # those "need new"
ii <- which(l.new)
iN <- (i.1 + i.2 + 1L) %/% 2
stopifnot( (i.1 < iN)[ii], (iN < i.2)[ii])
if(traceRootsearch) cat("ni new intervals, ni=",ni,"\n")
for(j in ii) {
## look in between i.1[j] .. i.2[j]
## NB: we can prove that i.1[j] < iN[j] < i.2[j]
rr[iN[j]] <- safeUroot(f.make(pp[iN [j]]), Sig = S,
lower= rr[i.1[j]],
upper= rr[i.2[j]],
tol=tol, maxiter=maxiter,
trace = traceRootsearch)$root
}
## update iDone[]:
seq_old <- seq_len(ni)
ni <- ni + length(ii)
iDone <- integer(ni)
iDone[seq_old + c(0L, cumsum(l.new))] <- oi
iDone[seq_along(ii) + ii ] <- iN[ii]
}
} ## else { method ..}
r[ip] <- if(hasDup) rr[i.pp] else rr
} ## end if(np)
r
}## end{qnorMix}
plot.norMix <-
function(x, type = "l", n = 511, xout = NULL, xlim = NULL, ylim,
xlab = "x", ylab = "f(x)", main = attr(x,"name"), lwd = 1.4,
p.norm = !p.comp, 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),
...)
{
## Purpose: plot method for "norMix" objects (normal mixtures)
## -------------------------------------------------------------------------
## Author: Martin Maechler, Date: 20 Mar 1997
if(!is.null(xlim) && is.null(xout)) ## determine xout
xout <- seq.int(xlim[1], xlim[2], length.out = n)
d.o <- dnorMixL(x, x = xout, n = n)
if(p.norm)
dn <- dnorm(d.o$x, mean = mean.norMix(x), sd = sqrt(var.norMix(x)))
if(!is.null(ll <- list(...)[["log"]]) && "y" %in% strsplit(ll,"")[[1]])
y0 <- max(1e-50, min(d.o$y, if(p.norm) dn))
else y0 <- 0
if(missing(ylim) || anyNA(ylim))
ylim <- c(y0, max(d.o$y, if(p.norm) dn))
plot(d.o, type = type, xlim = xlim, ylim = ylim,
main = main, xlab = xlab, ylab = ylab, lwd = lwd, ...)
if(p.norm) do.call(lines, c(list(x = d.o$x, y = dn), parNorm))
if(p.h0) do.call(abline, c(list(h = 0), parH0))
if(p.comp) {
m <- m.norMix(x) #-- number of components
w <- x[,"w"]; mu <- x[,"mu"]; sd <- x[,"sigma"]
for(j in 1:m)
do.call(lines,
c(list(x = d.o$x,
y = w[j] * dnorm(d.o$x, mean = mu[j], sd = sd[j])),
parComp))
}
invisible(x)
}
lines.norMix <-
function(x, type = "l", n = 511, xout = NULL, lwd = 1.4,
p.norm = FALSE, parNorm = list(col = 2, lty = 2, lwd = 0.4), ...)
{
## Purpose: lines method for "norMix" objects (normal mixtures)
## -------------------------------------------------------------
## Author: Martin Maechler, Date: 27 Jun 2002, 16:10
xlim <- if(is.null(xout)) par("usr")[1:2] # else NULL
d.o <- dnorMixL(x, x = xout, n = n, xlim = xlim)
lines(d.o, type = type, lwd = lwd, ...)
if(p.norm) {
dn <- dnorm(d.o$x, mean = mean.norMix(x), sd = sqrt(var.norMix(x)))
do.call(lines, c(list(x = d.o$x, y = dn), parNorm))
}
invisible()
}
r.norMix <- function(obj, x = NULL, xlim = NULL, n = 511, xy.return = TRUE)
{
## Purpose: Compute r := f / f0; f = normal mixture; f0 = "best" normal approx
## Author : Martin Maechler, Date: 20 Mar 97, 10:14
if(!is.norMix(obj)) stop("'obj' must be a 'Normal Mixture' object!")
d.o <- dnorMixL(obj, x, xlim = xlim, n = n)
dn <- dnorm(d.o$x, mean = mean.norMix(obj), sd = sqrt(var.norMix(obj)))
if(xy.return) list(x = d.o$x, y = d.o$y / dn, f0 = dn) else d.o$y / dn
}
## - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
## --> ../man/llnorMix.Rd for more documentation
## ~~~~~~~~~~~
nM2par <- function(obj, trafo = c("clr1", "logit"))
{
## Purpose: translate norMix object into our parametrization par.vector
## ----------------------------------------------------------------------
## Author: Martin Maechler, Date: 17 Dec 2007; "clr1": June 2019
stopifnot(is.norMix(obj))
trafo <- match.arg(trafo)
w <- switch(trafo,
"logit" = qlogis(obj[-1L, "w"]),
## clr := centered log ratio, by Aitchison (1986), see compositions :: clr
## clr1 : *omitting* the first entry
"clr1" = { ln <- log(obj[,"w"]); ln[-1L] - mean(ln) },
stop("invalid 'trafo': ", trafo))
## log(sqrt(.)) = log(.)/2
c(w, obj[,"mu"], log(obj[,"sigma"]))
}
.nM2par <- function(mu, sigma, trafo, w, check=TRUE)
{
## Purpose: Fast version of nM2par()
## -------------------------------------------------
## Author: Martin Maechler, Date: 18 Dec 2007
if(check) stopifnot(length(w) == (p <- length(mu)), length(sigma) == p)
c(switch(trafo,
"logit" = qlogis(w[-1L]),
"clr1" = { ln <- log(w); ln[-1L] - mean(ln) },
stop("invalid 'trafo': ", trafo)),
mu, log(sigma))
}
.par2nM <- function(p, trafo)
{
## Purpose: get (mu, sd, w) from our parametrization par.vector
## ----------------------------------------------------------------------
lp <- length(p)
stopifnot(is.numeric(p), lp %% 3 == 2)
m <- (lp + 1L) %/% 3
m1 <- m - 1L
names(p) <- NULL # so they are not transferred to mu,...
mu <- p[m:(m+m1)]
sd <- exp(p[(m+m):(m+m+m1)]) ## sigma = exp(tau)
if(m == 1)
list(mu=mu, sd=sd, w=1)
else { ## -- m >= 2
w <- switch(trafo,
"logit" = {
pi. <- plogis(p[1:m1]) ## \pi_j = inv_logit(\lambda_j)
if((sp <- sum(pi.)) > 1)
stop(sprintf("weights sum up to %.3g > 1 !", sp))
c(1 - sp, pi.)
},
"clr1" = {
## clr := centered log ratio, by Aitchison (1986), see compositions :: clr
## clr1 : *omitting* the first entry
## ln <- log(obj[,"w"]); ln[-1L] - mean(ln)
pM <- log(2) * .Machine$double.max.exp # = 709.7827 = log(.Machine$double.xmax)
p <- p[1:m1]
p <- c(-sum(p), p) # = (p_1,..., p_m)
if((mp <- max(p)) < pM) { ## normal case
sp <- sum(pi. <- exp(p)) # {TODO: exponential sum stable formula}
pi./sp
} else { ## extreme case, where exp(.) would overflow
kM <- which(p == mp)
p <- numeric(m) # all 0, apart from the max. value(s)
p[kM] <- 1/length(kM)
p
}
},
stop("invalid 'trafo': ", trafo))
list(mu=mu, sd=sd, w=w)
}
}
par2norMix <- function(p, trafo = c("clr1", "logit"),
name = sprintf("%s {from %s}",
trafo, deparse(substitute(p))[1]))
{
## Purpose: build norMix object from our parametrization par.vector
## ----------------------------------------------------------------------
## Author: Martin Maechler, Date: 17 Dec 2007
trafo <- match.arg(trafo)
force(name) # substitute(..)
nm <- .par2nM(p, trafo)
with(nm,
norMix(mu=mu, sigma = sd, w=w, name = name))
}
if(FALSE) ## this is not needed and mentioned in help(llnorMix) -- kept here as reminder
logLiknorMix <- function(obj, x) {
## Purpose: log-likelihood for 'norMix'
sum(dnorMix(x, obj, log=TRUE))
}
llnorMix <- function(p, x, m = (length(p)+1)/3, trafo = c("clr1", "logit"))
{
## Purpose: log-likelihood
## ----------------------------------------------------------------------
## Arguments: p : parameter vector, see below
## x : data vector
## m : number of mixture components
##
## 'p' is particularly parametrized,
## p = c( lambda_j, mu_j, tau_j) where
## depending on 'trafo' (before 2019: only "logit")
## trafo = "clr1" :
## \lambda_j = log(\pi_j) - mean_j'(log(\pi_j')), j=2,..,m
## trafo = "logit" [*the* only option before 2019-06]:
## \lambda_j = logit(\pi_j), j=2,..,m; and \pi_1 := 1- sum_j\pi_j
##
## and \tau_j = log(\sigma_j) such that parameters are unconstrained
## ----------------------------------------------------------------------
## Author: Martin Maechler, Date: 17 Dec 2007 (for hardwired "logit"); 2019-06
stopifnot(is.numeric(x), is.numeric(p), !is.matrix(p), 3*m == length(p)+1,
m == (m. <- as.integer(m)), (m <- m.) >= 1)
m1 <- m - 1L
mu <- p[m:(m+m1)]
sigma <- exp(p[(m+m):(m+m+m1)]) ## sigma = exp(tau)
if(m == 1)
return( sum(dnorm(x, mean = mu[1], sd = sigma[1], log = TRUE)) )
## else -- m >= 2
trafo <- match.arg(trafo)
w <- switch(trafo,
"logit" = {
pi. <- plogis(p[1:m1]) ## \pi_j = inv_logit(\lambda_j)
if((sp <- sum(pi.)) > 1) ## sum{1..K-1} pi[j] > 1
return(-Inf) # worst possible value
## as \pi_1 := 1 - sum_{j=2}^{m} \pi_j :
c(1 - sp, pi.)
},
"clr1" = {
## clr := centered log ratio, by Aitchison (1986), see compositions :: clr
## clr1 : *omitting* the first entry
## ln <- log(obj[,"w"]); ln[-1L] - mean(ln)
pM <- log(2) * .Machine$double.max.exp # = 709.7827, exp(pM) = double.xmax
p <- p[1:m1]
p <- c(-sum(p), p) # = (p_1,..., p_m)
pi. <- exp(if((mp <- max(p)) < pM) ## normal case
p
else ## extreme case, where exp(.) would overflow
p - mp)
pi. / sum(pi.)
},
stop("invalid 'trafo': ", trafo))
y <- 0
for(j in 1:m)
y <- y + w[j] * dnorm(x, mean = mu[j], sd = sigma[j])
## return
sum(log(y))
}
clus2norMix <- function(gr, x, name = deparse(sys.call()))
{
## Purpose: "Clustering to normal Mixture"
## ----------------------------------------------------------------------
## Arguments: gr: grouping/clustering vector (in {1,..,K}); possibly factor
## x : (original) data vector
## name : name for norMix() object; by default constructed
## ----------------------------------------------------------------------
## Author: Martin Maechler, Date: 31 Dec 2007, 10:24
if(length(gr) != (n <- length(x)))
stop("'gr' and 'x' are not of the same length")
r <- lapply(unname(split(x,gr)), ## << simple version of tapply(x, gr, *)
function(u){ nk <- length(u); m <- mean(u)
list(m, sum((u - m)^2)/(nk-1), nk) })
n. <- numeric(1)
norMix(mu = vapply(r, `[[`, n., 1L),
sigma= sqrt(vapply(r, `[[`, n., 2L)),
w = vapply(r, `[[`, n., 3L)/n,
name = name)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.