Nothing
## Copyright (C) 2013 Marius Hofert, Bernhard Pfaff
##
## 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 3 of the License, or (at your option) any later
## version.
##
## This program is distributed in the hope that it will be useful, but WITHOUT
## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
## FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
## details.
##
## You should have received a copy of the GNU General Public License along with
## this program; if not, see <http://www.gnu.org/licenses/>.
##
## Student's t distribution
##
## Random variates
rmt <- function(n, df = 4, mu = 0, Sigma){
d <- dim(Sigma)[1.]
chi <- 2.0 * rgamma(n, shape = df/2.0)
m1 <- rmnorm(n, mu = mu, Sigma = Sigma)
m2 <- matrix(rep(sqrt(df) / sqrt(chi), d), ncol = d)
mu.matrix <- matrix(mu, nrow = n, ncol = d, byrow = TRUE)
return(m1 * m2 + mu.matrix)
}
## Density
dmt <- function(x, df, mu, Sigma, log = FALSE){
d <- dim(x)[2]
Q <- mahalanobis(x, mu, Sigma)
log.const.top <- lgamma((df + d)/2)
log.const.bottom <- lgamma(df/2) + (d * log(pi * df)) / 2 + 0.5 * log(det(Sigma))
log.top <- (- (df + d) * log(1 + Q / df)) / 2
out <- log.const.top + log.top - log.const.bottom
if(!(log)) out <- exp(out)
out
}
##
qst <- function(p, mu = 0, sd = 1, df, scale = FALSE){
quant <- qt(p, df)
if (scale) quant <- quant * sqrt((df - 2) / df)
mu + sd * quant
}
## Fitting
fit.mst <- function(data, nit = 2000, tol = 1e-10, ...){
if(is.matrix(data) == FALSE) data <- as.matrix(data)
mix.pars <- c(-4, 8, 0) # initial mixture parameters lambda, chi, psi; for t: lambda = -nu/2, chi = nu, psi = 0; see MFE (2015; Section 6.2.3)
optpars <- c(2)
optfunc <- function(thepars, mixparams, greeks){
MCECM.Qfunc(-thepars / 2, thepars, mixparams[3], greeks[[1]], greeks[[2]], greeks[[3]])
}
n <- dim(data)[1]
d <- dim(data)[2]
Xbar <- apply(data, 2, mean)
mu <- Xbar # initial mu
Sigma <- var(data) # initial Sigma; should (?) probably multiply with E(W) = E(GIG) = \sqrt{\chi/\psi} K_{\lambda+1}(\sqrt{\chi\psi}) / K_{\lambda}(\sqrt{\chi\psi}) based on the given \chi, \psi, \lambda
gamma <- rep(0, d) # initial gamma; for t: gamma = 0
i <- 0
ll <- sum(dmt(data, mix.pars[2], mu, Sigma, log = TRUE)) # initial log-likelihood
closeness <- 100
while ((closeness > tol) & (i < nit)){
i <- i + 1
EMresult <- EMupdate(data, mix.pars, mu, Sigma, gamma, symmetric = TRUE, scaling = FALSE) # EM update of mu, Sigma, gamma; for t: mu, Sigma
mu <- EMresult$mu
Sigma <- EMresult$Sigma
gamma <- EMresult$gamma
MCECMresult <- MCECMupdate(data, mix.pars, mu, Sigma, gamma, optpars, optfunc, xieval = TRUE, ...) # MCECM update of mixture parameters
mix.pars <- MCECMresult$mix.pars # updated mixture parameters
mix.pars[1] <- -mix.pars[2] / 2
conv <- MCECMresult$conv
conv.type <- MCECMresult$convtype
ll.old <- ll # old log-likelihood
ll <- sum(dmt(data, mix.pars[2], mu, Sigma, log = TRUE)) # updated log-likelihood
closeness <- abs((ll - ll.old) / ll.old) # tolerance (relative error between log-likelihoods)
}
lambda <- mix.pars[1]
chi <- mix.pars[2]
psi <- mix.pars[3]
EW <- EGIG(lambda, chi, psi)
EW2 <- EGIG(lambda, chi, psi, 2)
varW <- EW2 - EW^2 # var(W)
Sigma <- forceSymmetric(Sigma)
beta <- as.vector(solve(Sigma) %*% gamma)
mean <- as.numeric(mu + EW * gamma)
covariance <- EW * Sigma + varW * outer(gamma, gamma) # for t: cov(X) = EW * Sigma
if (d > 1){
cor <- cov2cor(Sigma)
} else {
cor <- 1
}
nu <- chi
list(df = nu, mu = mu, Sigma = Sigma, gamma = gamma, ll.max = ll, mean = mean,
covariance = covariance, correlation = cor)
}
##
fit.st <- function(data, ...){
if(is.timeSeries(data)) data <- series(data)
mu <- mean(data)
m2 <- mean((data - mu)^2)
m4 <- mean((data - mu)^4)
nu <- 4 + (6 * m2^2) / (m4 - 3 * m2^2)
sigma <- sqrt((nu - 2) * m2 / nu)
theta <- c(nu, mu, sigma)
negloglik <- function(theta, y){
- sum(log(dt((y - theta[2]) / abs(theta[3]), df = abs(theta[1]))) - log(abs(theta[3])))
}
optimfit <- optim(theta, fn = negloglik, y = data, ...)
par.ests <- optimfit$par
ifelse(optimfit$convergence == 0, converged <- TRUE, converged <- FALSE)
par.ests[1] <- abs(par.ests[1])
par.ests[2] <- abs(par.ests[2])
nItheta <- hessian(negloglik, par.ests, y = data)
asymp.cov <- solve(nItheta)
loglh.max <- -negloglik(par.ests, y = data)
par.ses <- sqrt(diag(asymp.cov))
names(par.ests) <- c("nu", "mu", "sigma")
names(par.ses) <- names(par.ests)
dimnames(asymp.cov) <- list(names(par.ests), names(par.ests))
list(converged = converged, par.ests = par.ests, par.ses = par.ses,
asymp.cov = asymp.cov, ll.max = loglh.max)
}
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.