Nothing
## last modified August 29, 2002
mix <- function(mixdat, mixpar, dist = "norm", constr = list(conpi = "NONE",
conmu = "NONE", consigma = "NONE", fixpi = NULL, fixmu = NULL,
fixsigma = NULL, cov = NULL, size = NULL), emsteps = 1, usecondit = FALSE,
exptol = 5e-06, print.level = 0, ...)
{
k <- nrow(mixpar)
if (usecondit & ncol(mixdat) - 2 != k)
stop("Conditional data are not consistent with mixpar.")
if (usecondit & k == 1)
warning("Only 1 component, conditional data will not be used")
if (is.na(match(dist, c("norm", "lnorm", "gamma", "weibull",
"binom", "nbinom", "pois"))))
stop(paste("Unknown distribution ", dist, ".", sep = ""))
if (!testconstr(mixdat, mixpar, dist, constr))
stop("Invalid constraints.")
if (!testpar(mixpar, dist, constr))
stop("Invalid parameters.")
fitpar <- mixpar
if (print.level > 0)
print(fitpar)
estep <- function(emixdat, emixpar, edist, econstr) {
prj.i <- grpintprob(emixdat, emixpar, edist, econstr)
prji <- sweep(prj.i, 2, emixpar[, 1], "*")
pri.j <- sweep(prji, 1, apply(prji, 1, sum), "/")
pseudodat <- sweep(pri.j, 1, emixdat[, 2], "*")
nplusi <- apply(pseudodat, 2, sum)
pi <- nplusi/sum(nplusi)
if (econstr$conpi == "PFX") {
pi[econstr$fixpi] <- emixpar[econstr$fixpi, 1]
pi[!econstr$fixpi] <- pi[!econstr$fixpi] * (1 - sum(pi[econstr$fixpi]))/sum(pi[!econstr$fixpi])
}
list(pi = pi, pseudomixdat = data.frame(x = emixdat[,
1], pseudodat = pseudodat))
}
mstep <- function(mspseudomixdat, msmixpar, msdist, msconstr,
msprint.level, ...) {
multineg2llg <- function(mulpseudomixdat, mulmixpar,
muldist, mulconstr, multheta) {
muly <- as.matrix(mulpseudomixdat[, -1])
muln <- apply(muly, 2, sum)
mulpara <- mixtheta2par(multheta, mulmixpar, mulconstr,
mixprop = FALSE)
upmixpar <- data.frame(mulpi = mulmixpar[, 1], mulmu = mulpara[,
2], mulsigma = mulpara[, 3])
mulparvalid <- testpar(upmixpar, muldist, mulconstr)
if (mulparvalid) {
mulp <- grpintprob(mulpseudomixdat, upmixpar,
muldist, mulconstr)
mulnp <- sweep(mulp, 2, muln, "*")
mulyu <- as.vector(muly)
mulnpu <- as.vector(mulnp)[mulyu > 0]
mulyu <- mulyu[mulyu > 0]
min(-2 * sum(mulyu * log(mulnpu/mulyu)), 1e+16)
}
else 1e+16
}
theta0 <- mixpar2theta(msmixpar, msconstr, mixprop = FALSE)
msnlmobj <- nlm(multineg2llg, mulpseudomixdat = mspseudomixdat,
mulmixpar = msmixpar, muldist = msdist, mulconstr = msconstr,
p = theta0, print.level = msprint.level, ...)
msk <- nrow(msmixpar)
mspara <- mixtheta2par(msnlmobj$est, msmixpar, msconstr,
mixprop = FALSE)
msmu <- mspara[, 2]
mssigma <- mspara[, 3]
data.frame(pi = msmixpar[, 1], mu = msmu, sigma = mssigma)
}
if (!usecondit & emsteps > 0 & k > 1 & (constr$conpi == "NONE" |
(constr$conpi == "PFX" & sum(constr$fixpi) < k - 1))) {
for (i in 1:emsteps) {
fitpi <- estep(mixdat, fitpar, dist, constr)
fitpar[, 1] <- fitpi$pi
fitpar <- mstep(fitpi$pseudomixdat, fitpar, dist,
constr, print.level, ...)
if (print.level > 0)
print(fitpar)
}
}
mixtheta0 <- mixpar2theta(fitpar, constr)
mixlike <- function(lmixdat, lmixpar, ldist, lconstr, lmixtheta,
lusecondit, lexptol) {
lk <- nrow(lmixpar)
lm <- nrow(lmixdat)
uppar <- mixtheta2par(lmixtheta, lmixpar, lconstr)
parvalid <- testpar(uppar, ldist, lconstr)
if (!parvalid) {
like <- 1e+16
ldf <- lm - 1
}
else {
pmat <- grpintprob(lmixdat, uppar, ldist, lconstr)
ln <- sum(lmixdat[, 2])
joint <- sweep(pmat, 2, uppar[, 1], "*")
mixed <- apply(joint, 1, sum)
ly <- as.vector(as.matrix(lmixdat[, 2]))
if (lusecondit & lk > 1) {
conditional <- sweep(joint, 1, mixed, "/")
conditdat <- as.matrix(lmixdat[, -(1:2)])
subsam <- apply(conditdat, 1, sum)
nconditpr <- sweep(conditional, 1, subsam, "*")
conditdatu <- as.vector(conditdat)
nconditpru <- as.vector(nconditpr)[conditdatu >
0]
conditdatu <- conditdatu[conditdatu > 0]
}
pu <- as.vector(mixed)[ly > 0]
npu <- ln * pu
yu <- ly[ly > 0]
if (lusecondit) {
like <- min(-2 * (sum(yu * log(npu/yu)) + sum(conditdatu *
log(nconditpru/conditdatu))), 1e+16)
ldf <- sum(mixed > lexptol) - 1 + sum(nconditpr >
lexptol) - sum(subsam > 0) - length(lmixtheta)
}
else {
like <- min(-2 * sum(yu * log(npu/yu)), 1e+16)
ldf <- sum(mixed > lexptol) - 1 - length(lmixtheta)
}
}
attr(like, "df") <- ldf
like
}
nlmobj <- nlm(mixlike, lmixdat = mixdat, lmixpar = fitpar,
ldist = dist, lconstr = constr, lusecondit = usecondit,
lexptol = exptol, p = mixtheta0, hessian = TRUE, print.level = print.level,
...)
if (nlmobj$code == 3)
warning("The optimization process terminated because either the estimates are approximate local optimal solution or steptol is too small")
else if (nlmobj$code == 4)
warning("The optimization process terminated because iteration limit exceeded")
else if (nlmobj$code == 5)
warning("The optimization process terminated because either the function is unbounded below or stepmax is too small")
if (nlmobj$minimum == 1e+16)
warning("Invalid parameters")
param <- mixtheta2par(nlmobj$est, mixpar, constr)
covmat <- function(covmixpar, covconstr, covhessian) {
invmat <- try(solve(covhessian/2))
if (inherits(invmat, "try-error"))
invmat <- matrix(NA, nrow = nrow(covhessian), ncol = ncol(covhessian))
dr <- nrow(invmat)
covk <- nrow(covmixpar)
if (covconstr$conpi == "NONE")
lpi <- covk - 1
else if (covconstr$conpi == "PFX" & sum(covconstr$fixpi) <
covk - 1) {
pi.e <- covmixpar[!covconstr$fixpi, 1]
lpi <- length(pi.e) - 1
}
else lpi <- 0
if (covconstr$conmu == "NONE")
mu.e <- covmixpar[, 2]
else if (covconstr$conmu == "MFX")
mu.e <- covmixpar[!covconstr$fixmu, 2]
else if (covconstr$conmu == "MEQ")
mu.e <- covmixpar[1, 2]
else if (covconstr$conmu == "MES")
mu.e <- covmixpar[1:2, 2]
else if (covconstr$conmu == "MGC")
mu.e <- covmixpar[1:3, 2]
lmu <- length(mu.e)
if (covconstr$consigma == "NONE")
sigma.e <- covmixpar[, 3]
else if (covconstr$consigma == "SFX")
sigma.e <- covmixpar[!constr$fixsigma, 3]
else if (!is.na(match(covconstr$consigma, c("SEQ", "CCV"))))
sigma.e <- covmixpar[1, 3]
else if (!is.na(match(covconstr$consigma, c("FCV", "BINOM",
"NBINOM", "POIS"))))
sigma.e <- NULL
lsigma <- length(sigma.e)
sigmat <- diag(c(rep(1, dr - lsigma), sigma.e))
vmat <- sigmat %*% invmat %*% sigmat
se <- sqrt(diag(vmat))
pi.se <- rep(NA, covk)
if (lpi > 0) {
if (covconstr$conpi == "PFX")
pi.se[!covconstr$fixpi] <- c(se[1:lpi], sqrt(sum(vmat[1:lpi,
1:lpi])))
else pi.se[1:covk] <- c(se[1:lpi], sqrt(sum(vmat[1:lpi,
1:lpi])))
}
mu.se <- rep(NA, covk)
if (lmu > 0) {
if (covconstr$conmu == "MFX")
mu.se[!covconstr$fixmu] <- se[(lpi + 1):(lpi +
lmu)]
else mu.se[1:lmu] <- se[(lpi + 1):(lpi + lmu)]
}
sigma.se <- rep(NA, covk)
if (lsigma > 0) {
if (covconstr$consigma == "SFX")
sigma.se[!covconstr$fixsigma] <- se[(lpi + lmu +
1):(lpi + lmu + lsigma)]
else sigma.se[1:lsigma] <- se[(lpi + lmu + 1):(lpi +
lmu + lsigma)]
}
list(vmat = vmat, mixse = data.frame(pi.se = pi.se, mu.se = mu.se,
sigma.se = sigma.se))
}
vmat <- covmat(param, constr, nlmobj$hessian)
chisq <- nlmobj$minimum
df <- attr(mixlike(mixdat, fitpar, dist, constr, nlmobj$est,
usecondit, exptol), "df")
P <- 1 - pchisq(chisq, df)
mixfit <- list(parameters = param, se = vmat$mixse, distribution = dist,
constraint = constr, chisq = chisq, df = df, P = P, vmat = vmat$vmat,
mixdata = mixdat, usecondit = usecondit)
class(mixfit) <- "mix"
mixfit
}
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.