Nothing
# This library is free software; you can redistribute it and/or
# modify it under the terms of the GNU Library General Public
# License as published by the Free Software Foundation; either
# version 2 of the License, or (at your option) any later version.
#
# This library 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 Library General Public License for more details.
#
# You should have received a copy of the GNU Library General
# Public License along with this library; if not, write to the
# Free Foundation, Inc., 59 Temple Place, Suite 330, Boston,
# MA 02111-1307 USA
# Copyrights (C)
# for this R-port:
# 1999 - 2006, Diethelm Wuertz, GPL
# Diethelm Wuertz <wuertz@itp.phys.ethz.ch>
# info@rmetrics.org
# www.rmetrics.org
# for the code accessed (or partly included) from other R-ports:
# see R's copyright and license files
# for the code accessed (or partly included) from contributed R-ports
# and other sources
# see Rmetrics's copyright file
################################################################################
# FUNCTION: DESCRIPTION:
# 'fDISTFIT' S4 Class representation
# print.fDISTFIT Prints Results from a Fitted Distribution
# FUNCTION: NORMAL AND T DISTRIBUTION:
# nFit Fits parameters of a Normal density
# tFit Fits parameters of a Student-t density
# FUNCTION: STABLE DISTRIBUTION:
# stableFit Fits parameters of a stable density
# FUNCTION: GENERALIZED DISTRIBUTION:
# ghFit Fits parameters of a generalized hyperbolic density
# hypFit Fits parameters of a hyperbolic density
# nigFit Fits parameters of a normal inverse Gaussian density
################################################################################
test.nFit <- function()
{
# Graph Frame:
par(mfrow = c(1, 1))
# Simulate normal random variates:
## GNB: was
## RNGkind(kind = "Marsaglia-Multicarry", normal.kind = "Inversion")
## set.seed(4711, kind = "Marsaglia-Multicarry")
set.seed(4711)
s = rnorm(n = 2000, mean = 1, sd = 0.5)
# Fit:
cat("nFit():\n")
print(ans <- nFit(x = s))
# Precision within 10% ?
relErrorTest = c(
( (ans@fit$estimate[1] - 1.0)/1.0 < 0.10 ),
( (ans@fit$estimate[2] - 0.5)/0.5 < 0.10 ))
print(relErrorTest)
checkTrue(as.logical(mean(relErrorTest)))
# Return Value:
return()
}
# ------------------------------------------------------------------------------
test.tFit =
function()
{
# Graph Frame:
par(mfrow = c(1, 1))
# Simulate random variates t(4):
## GNB:
## RNGkind(kind = "Marsaglia-Multicarry", normal.kind = "Inversion")
## set.seed(4711, kind = "Marsaglia-Multicarry")
set.seed(4711)
s = rt(n = 2000, df = 4)
# Fit:
ans = tFit(x = s, df = 2*var(s)/(var(s)-1), trace = FALSE)
# Precision of df within 10% ?
relErrorTest = c(
( (ans@fit$estimate[1] - 4.0)/4.0 < 0.10 ))
print(ans)
print(relErrorTest)
checkTrue(as.logical(mean(relErrorTest)))
# Return Value:
return()
}
# ------------------------------------------------------------------------------
test.ghFit =
function()
{
# Graph Frame:
par(mfrow = c(1, 1))
# Simulate random variates:
## RNGkind(kind = "Marsaglia-Multicarry", normal.kind = "Inversion")
## set.seed(4711, kind = "Marsaglia-Multicarry")
set.seed(4711)
s = rgh(n = 2000, alpha = 0.8, beta = 0.2, delta = 2, mu = -0.4,
lambda = 1)
# Fit:
ans = ghFit(x = s, alpha = 1, beta = 0, delta = 1, mu = 0, lambda = 1,
trace = FALSE)
# Precision of parameters within 10% ?
relErrorTest = c(
( (ans@fit$estimate[1] - 0.8)/( 0.8) < 0.10 ),
( (ans@fit$estimate[2] - 0.2)/( 0.2) < 0.10 ),
( (ans@fit$estimate[3] - 2.0)/( 2.0) < 0.10 ),
( (ans@fit$estimate[4] + 0.4)/(-0.4) < 0.10 ),
( (ans@fit$estimate[5] - 2.0)/( 2.0) < 0.10 ))
print(ans)
print(relErrorTest)
checkTrue(as.logical(mean(relErrorTest)))
# Return Value:
return()
}
# ------------------------------------------------------------------------------
test.hypFit =
function()
{
# Graph Frame:
par(mfrow = c(1, 1))
# Simulate normal random variates:
## RNGkind(kind = "Marsaglia-Multicarry", normal.kind = "Inversion")
## set.seed(4711, kind = "Marsaglia-Multicarry")
set.seed(4711)
s = rhyp(n = 2000, alpha = 1.5, beta = 0.8, delta = 0.5, mu = -1)
# Fit:
ans = hypFit(s, alpha = 1, beta = 0, delta = 1, mu = mean(s),
trace = FALSE)
# Precision of parameters within 10% ?
relErrorTest = c(
( (ans@fit$estimate[1] - 1.5)/( 1.5) < 0.10 ),
( (ans@fit$estimate[2] - 0.8)/( 0.8) < 0.10 ),
( (ans@fit$estimate[3] - 0.5)/( 0.5) < 0.10 ),
( (ans@fit$estimate[4] + 1.0)/(-1.0) < 0.10 ))
print(ans)
print(relErrorTest)
checkTrue(as.logical(mean(relErrorTest)))
# Return Value:
return()
}
# ------------------------------------------------------------------------------
test.nigFit <- function()
{
# Graph Frame:
par(mfrow = c(1, 1))
# Simulate normal random variates:
## RNGkind(kind = "Marsaglia-Multicarry", normal.kind = "Inversion")
## set.seed(4711, kind = "Marsaglia-Multicarry")
set.seed(4712)
th.0 <- c(alpha = 1.5, beta = -0.7, delta = 0.5, mu = -1.0)
s <- do.call(rnig, c(list(n = 2000), th.0))
# Fit {includes plot !}
print(ans <- nigFit(s, alpha = 1, beta = 0, delta = 1, mu = mean(s),
trace = FALSE))
# Precision of parameters within 10% / 35% ?
est <- ans@fit$estimate
relErr <- (est - th.0) / th.0
checkTrue(all(abs(relErr) < 0.35))# 35%
## 10% / 35$ is silly, rather use SE scale:
print(SE <- ans@fit$se.coef)
print(cbind(th.0, est, SE))# NaN
SE[is.na(SE)] <- .1 # arbitrary
cat("t-ratios:\n"); print(t. <- (est - th.0)/ SE)
checkTrue(all(abs(t.) < 6))## 6 is too large, but have -5.476 for beta
# Return Value:
return()
}
# ------------------------------------------------------------------------------
test.stableFit <- function()
{
# Graph Frame:
par(mfrow = c(1, 1))
# Simulate stable random variates:
## RNGkind(kind = "Marsaglia-Multicarry", normal.kind = "Inversion")
## set.seed(4711, kind = "Marsaglia-Multicarry")
set.seed(4713)
th.0 <- c(alpha=1.8, beta=0.3, gamma = 1, delta = 0.1)
s <- do.call(stabledist::rstable, c(list(n = 500), th.0))
# Fit {includes plot !}
cat("stableFit():\n")
print(ans <- stableFit(x = s, alpha = 1.5))
# Precision of parameters within 5% / 10% / ...?
est <- ans@fit$estimate
arelErr <- abs(relErr <- (est - th.0) / th.0)
print(cbind(th.0, est, relErr))
checkTrue(all(arelErr[c("alpha","gamma")] < 0.05) &&
arelErr[["beta"]] < 0.20 && arelErr[["delta"]] < 0.50)
# MLE Fit:
if (FALSE) { # as this takes rather too long time ...
print(system.time(
ans <- stableFit(x = s, alpha = 1.5, type = "mle", trace = TRUE)
## .mleStableFit(s, 1.75, 0, 1, 0)
))
## user {lynne, 2012-09}
## 159.7
# The result would be:
print(ans)
## Precision of parameters within 5% / 10% / ...?
est <- ans@fit$estimate
arelErr <- abs(relErr <- (est - th.0) / th.0)
print(cbind(th.0, est, relErr))
checkTrue(all(arelErr[1:3] < 0.02) && arelErr[["delta"]] < 0.42)
}
# Return Value:
return()
}
################################################################################
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.