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
################################################################################
# FUNCTION: DESCRIPTION:
# dnig Returns density for inverse Gaussian DF
# pnig Returns probability for for inverse Gaussian DF
# qnig Returns quantiles for for inverse Gaussian DF
# rnig Returns random variates for inverse Gaussian DF
# FUNCTION: DESCRIPTION:
# .pnigC Fast C implementation
# .qnigC Fast C implementation
################################################################################
.dnig <- function(x, alpha, beta, delta, mu, log = FALSE)
{
x. <- x - mu
log.a <- delta*sqrt(alpha^2-beta^2) + log(delta*alpha/pi)
Sqrt <- sqrt(delta^2 + x.^2)
log.K1 <- log(besselK(alpha * Sqrt, 1, expon.scaled = TRUE)) - alpha*Sqrt
dnig <- log.a -log(Sqrt) + log.K1 + beta*x.
if(log) dnig else exp(dnig)
}
dnig <-
function(x, alpha = 1, beta = 0, delta = 1, mu = 0, log = FALSE)
{
# A function implemented by Diethelm Wuertz
# Description:
# Returns density for inverse Gaussian DF
# Example:
# x = rnorm(1000); u = dgh(x, 1.1, 0.2, 0.8, 0.4, -0.5)
# v = dnig(rnorm(x), 1.1, 0.2, 0.8, 0.4); u -v
# FUNCTION:
# Parameters:
if (length(alpha) == 4) {
mu = alpha[4]
delta = alpha[3]
beta = alpha[2]
alpha = alpha[1]
}
# Checks:
if (alpha <= 0) stop("alpha must be greater than zero")
if (delta <= 0) stop("delta must be greater than zero")
if (abs(beta) >= alpha) stop("abs value of beta must be less than alpha")
.dnig(x, alpha=alpha, beta=beta, delta=delta, mu=mu, log=log)
}
# ------------------------------------------------------------------------------
pnig <-
function(q, alpha = 1, beta = 0, delta = 1, mu = 0)
{
# A function implemented by Diethelm Wuertz
# Description:
# Returns probability for for inverse Gaussian DF
# Function:
# Probability:
# pgh(q = q, alpha = alpha, beta = beta, delta = delta, mu = mu,
# lambda = -0.5)
if (alpha <= 0)
stop("alpha must be greater than zero")
if (delta <= 0)
stop("delta must be greater than zero")
if (abs(beta) >= alpha)
stop("abs value of beta must be less than alpha")
ans <- q
for(i in seq_along(q)) {
ans[i] <- integrate(.dnig, -Inf, q[i],
stop.on.error = FALSE,
alpha = alpha, beta = beta,
delta = delta, mu = mu)$value
}
ans
}
# ------------------------------------------------------------------------------
qnig <-
function(p, alpha = 1, beta = 0, delta = 1, mu = 0)
{
# A function implemented by Diethelm Wuertz
# Description:
# Returns quantiles for for inverse Gaussian DF
# FUNCTION:
# Quantiles:
# qgh(p = p, alpha = alpha, beta = beta, delta = delta, mu = mu,
# lambda = -0.5)
# Checks:
if (alpha <= 0) stop("alpha must be greater than zero")
if (delta <= 0) stop("delta must be greater than zero")
if (abs(beta) >= alpha) stop("abs value of beta must be less than alpha")
# Internal Function:
.froot <- function(x, alpha, beta, delta, mu, p)
pnig(q = x, alpha = alpha, beta = beta, delta = delta, mu = mu) - p
# Quantiles:
ans <- p
for(i in seq_along(p)) {
lower = -1
upper = +1
counter = 0
v <- NA
while(is.na(v) && counter < 1000) {
v <- .unirootNA(f = .froot, interval = c(lower, upper),
alpha = alpha, beta = beta, delta = delta,
mu = mu, p = p[i])
counter = counter + 1
lower = lower - 2^counter
upper = upper + 2^counter
}
ans[i] <- v
}
# Return Value:
ans
}
# ------------------------------------------------------------------------------
rnig <-
function(n, alpha = 1, beta = 0, delta = 1, mu = 0)
{
# A function implemented by Diethelm Wuertz
# Description:
# Return normal inverse Gaussian distributed random variates
# Arguments:
# n - number of deviates to be generated
# alpha, beta - Shape Parameter, |beta| <= alpha
# delta - Scale Parameter, 0 <= delta
# mu - Location Parameter
# FUNCTION:
# Settings:
gamma = sqrt(alpha*alpha - beta*beta)
# GAMMA:
if (gamma == 0) {
# GAMMA = 0:
V = rnorm(n)^2
Z = delta*delta / V
X = sqrt(Z)*rnorm(n)
} else {
# GAMMA > 0:
U = runif(n)
V = rnorm(n)^2
# FIXED ...
z1 <- function(v, delta, gamma) {
delta/gamma + v/(2*gamma^2) - sqrt(v*delta/(gamma^3) +
(v/(2*gamma^2))^2 )
}
z2 <- function(v, delta, gamma) {
(delta/gamma)^2 / z1(v = v, delta = delta, gamma = gamma)
}
pz1 <- function(v, delta, gamma) {
delta / (delta + gamma * z1(v = v, delta = delta, gamma = gamma) )
}
s = (1-sign(U-pz1(v = V, delta = delta, gamma = gamma)))/2
Z = z1(v = V, delta = delta, gamma = gamma)*s + z2(v = V, delta =
delta, gamma = gamma)*(1-s)
X = mu + beta*Z + sqrt(Z)*rnorm(n)
}
# Return Value:
X
}
################################################################################
.qnigC <-
function(p, alpha = 1, beta = 0, delta = 1, mu = 0)
{
# Description:
# Returns quantiles for for inverse Gaussian DF
# Example:
# p = runif(10); .qnigC(p)
# FUNCTION:
# Checks:
if(alpha <= 0) stop("Invalid parameters: alpha <= 0.0\n")
if(alpha^2 <= beta^2) stop("Invalid parameters: alpha^2 <= beta^2.0\n")
if(delta <= 0) stop("Invalid parameters: delta <= 0.0\n")
if((sum(is.na(p)) > 0))
stop("Invalid probabilities:\n",p,"\n")
else
if(sum(p < 0)+sum(p > 1) > 0) stop("Invalid probabilities:\n",p,"\n")
# Evaluate NIG cdf by calling C function
n <- length(p)
q <- rep(0, n)
retValues <- .C(C_qNIG,
as.double(p),
as.double(mu),
as.double(delta),
as.double(alpha),
as.double(beta),
as.integer(n),
as.double(q))
quantiles <- retValues[[7]]
quantiles[quantiles <= -1.78e+308] <- -Inf
quantiles[quantiles >= 1.78e+308] <- Inf
# Return Value:
quantiles
}
# ------------------------------------------------------------------------------
.pnigC <-
function(q, alpha = 1, beta = 0, delta = 1, mu = 0)
{
# Description:
# Returns probabilities for for inverse Gaussian DF
# IMPORTANT NOTE:
# DW: C program fails, remains to check
# Example:
# q = sort(rnorm(10)); .pnigC(q) # FAILS
# FUNCTION:
# Checks:
if(alpha <= 0) stop("Invalid parameters: alpha <= 0.0\n")
if(alpha^2 <= beta^2) stop("Invalid parameters: alpha^2 <= beta^2.0\n")
if(delta <= 0) stop("Invalid parameters: delta <= 0.0\n")
if(sum(is.na(q)) > 0) stop("Invalid quantiles:\n", q)
# Evaluate NIG cdf by calling C function
n <- length(q)
p <- rep(0, n)
retValues <- .C(C_pNIG,
as.double(q),
as.double(mu),
as.double(delta),
as.double(alpha),
as.double(beta),
as.integer(n),
as.double(p))
probs <- retValues[[7]]
# Return Value:
probs
}
################################################################################
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.