umontreal.iro.lecuyer.probdist.NormalInverseGaussianDist Maven / Gradle / Ivy
Go to download
Show more of this group Show more artifacts with this name
Show all versions of ssj Show documentation
Show all versions of ssj Show documentation
SSJ is a Java library for stochastic simulation, developed under the direction of Pierre L'Ecuyer,
in the Département d'Informatique et de Recherche Opérationnelle (DIRO), at the Université de Montréal.
It provides facilities for generating uniform and nonuniform random variates, computing different
measures related to probability distributions, performing goodness-of-fit tests, applying quasi-Monte
Carlo methods, collecting (elementary) statistics, and programming discrete-event simulations with both
events and processes.
The newest version!
/*
* Class: NormalInverseGaussianDist
* Description: normal inverse gaussian distribution
* Environment: Java
* Software: SSJ
* Copyright (C) 2001 Pierre L'Ecuyer and Université de Montréal
* Organization: DIRO, Université de Montréal
* @author
* @since
* SSJ is free software: you can redistribute it and/or modify it under
* the terms of the GNU General Public License (GPL) as published by the
* Free Software Foundation, either version 3 of the License, or
* any later version.
* SSJ 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.
* A copy of the GNU General Public License is available at
GPL licence site.
*/
package umontreal.iro.lecuyer.probdist;
import umontreal.iro.lecuyer.probdist.NormalDist;
import umontreal.iro.lecuyer.util.Num;
/**
* Extends the class {@link ContinuousDistribution} for
* the normal inverse gaussian distribution with location parameter
* μ, scale parameter
* δ > 0, tail heavyness
* α > 0, and
* asymmetry parameter β such that
* 0 <= | β| < α.
* Its density is
*
*
*
* f (x) = αδeδγ+β(x-μ)K1(α(δ^2 + (x - μ)^2)1/2)/π(δ^2 + (x - μ)^2)1/2,
*
* where K1 is the modified Bessel function of the second kind of order 1,
* and
* γ = (α^2 - β^2)1/2.
*
*
* The distribution function is given by
*
*
*
* F(x) = ∫-∞xdtf (t),
*
*
*/
public class NormalInverseGaussianDist extends ContinuousDistribution {
protected double alpha;
protected double beta;
protected double gamma;
protected double delta;
protected double mu;
/**
* Constructor for a normal inverse gaussian distribution with parameters α = alpha,
* β = beta, μ = mu and δ = delta.
*
*/
public NormalInverseGaussianDist (double alpha, double beta, double mu,
double delta) {
setParams (alpha, beta, mu, delta);
}
public double density (double x) {
return density (alpha, beta, mu, delta, x);
}
public double cdf (double x) {
return cdf (alpha, beta, mu, delta, x);
}
public double barF (double x) {
return barF (alpha, beta, mu, delta, x);
}
public double getMean() {
return getMean (alpha, beta, mu, delta);
}
public double getVariance() {
return getVariance (alpha, beta, mu, delta);
}
public double getStandardDeviation() {
return getStandardDeviation (alpha, beta, mu, delta);
}
/**
* Computes the density function
* for the normal inverse gaussian distribution with parameters α, β, μ
* and δ, evaluated at x.
*
*/
public static double density (double alpha, double beta, double mu,
double delta, double x) {
if (delta <= 0.0)
throw new IllegalArgumentException ("delta <= 0");
if (alpha <= 0.0)
throw new IllegalArgumentException ("alpha <= 0");
if (Math.abs(beta) >= alpha)
throw new IllegalArgumentException ("|beta| >= alpha");
double gamma = Math.sqrt(alpha*alpha - beta*beta);
double z = (x - mu)/delta;
double w;
if (Math.abs(z) <= 1.0e10)
w = Math.sqrt (1.0 + z*z);
else
w = Math.abs(z);
double y = alpha*delta*w;
double v = delta*(gamma + beta*z);
double R = Num.expBesselK1(v, y);
return alpha * R / (Math.PI*w);
}
/**
* NOT IMPLEMENTED.
* Computes the distribution function
* of the normal inverse gaussian distribution with parameters α,
* β, μ and δ, evaluated at x.
*
*/
public static double cdf (double alpha, double beta, double mu,
double delta, double x) {
if (delta <= 0.0)
throw new IllegalArgumentException ("delta <= 0");
if (alpha <= 0.0)
throw new IllegalArgumentException ("alpha <= 0");
if (Math.abs(beta) >= alpha)
throw new IllegalArgumentException ("|beta| >= alpha");
double gamma = Math.sqrt(alpha*alpha - beta*beta);
double z = (x - mu)/delta;
if (z > 0.0 && (gamma + (beta - alpha)*z >= XBIG))
return 1.0;
if (z < 0.0 && (gamma + (beta + alpha)*z <= -XBIGM))
return 0.0;
// double w = Math.sqrt (1.0 + z*z);
throw new UnsupportedOperationException
("NormalInverseGaussianDist: cdf NOT IMPLEMENTED");
}
/**
* NOT IMPLEMENTED.
* Computes the complementary distribution function of the normal inverse gaussian distribution
* with parameters α, β, μ and δ, evaluated at x.
*
*/
public static double barF (double alpha, double beta, double mu,
double delta, double x) {
return 1.0 - cdf (alpha, beta, mu, delta, x);
}
/**
* NOT IMPLEMENTED. Computes the inverse of the normal inverse gaussian distribution
* with parameters α, β, μ and δ.
*
*/
public static double inverseF (double alpha, double beta, double mu,
double delta, double u) {
throw new UnsupportedOperationException(" Inversion NOT IMPLEMENTED");
}
/**
* NOT IMPLEMENTED.
*
* @param x the list of observations used to evaluate parameters
*
* @param n the number of observations used to evaluate parameters
*
* @return returns the parameters
* [
* hat(α),
* hat(β), hat(μ),
* hat(δ)]
*
*/
public static double[] getMLE (double[] x, int n) {
if (n <= 0)
throw new IllegalArgumentException ("n <= 0");
/*
double[] parameters = new double[4];
double sum = 0;
for (int i = 0; i < n; i++) {
sum += x[i];
}
*/
throw new UnsupportedOperationException("getMLE is not implemented");
// return parameters;
}
/**
* NOT IMPLEMENTED.
*
* @param x the list of observations to use to evaluate parameters
*
* @param n the number of observations to use to evaluate parameters
*
*
*/
public static NormalInverseGaussianDist getInstanceFromMLE (double[] x,
int n) {
double parameters[] = getMLE (x, n);
return new NormalInverseGaussianDist (parameters[0], parameters[1],
parameters[2], parameters[3]);
}
/**
* Returns the mean
* E[X] = μ + δβ/γ of the
* normal inverse gaussian distribution with parameters α, β, μ and δ.
*
* @return the mean of the normal inverse gaussian distribution
*
* E[X] = μ + δβ/γ
*
*/
public static double getMean (double alpha, double beta, double mu,
double delta) {
if (delta <= 0.0)
throw new IllegalArgumentException ("delta <= 0");
if (alpha <= 0.0)
throw new IllegalArgumentException ("alpha <= 0");
if (Math.abs(beta) >= alpha)
throw new IllegalArgumentException ("|beta| >= alpha");
double gamma = Math.sqrt(alpha*alpha - beta*beta);
return mu + delta*beta/gamma;
}
/**
* Computes and returns the variance
* Var[X] = δα2/γ3 of the normal inverse gaussian distribution with parameters
* α, β, μ and δ.
*
* @return the variance of the normal inverse gaussian distribution
*
* Var[X] = δα2/γ3
*
*/
public static double getVariance (double alpha, double beta, double mu,
double delta) {
if (delta <= 0.0)
throw new IllegalArgumentException ("delta <= 0");
if (alpha <= 0.0)
throw new IllegalArgumentException ("alpha <= 0");
if (Math.abs(beta) >= alpha)
throw new IllegalArgumentException ("|beta| >= alpha");
double gamma = Math.sqrt(alpha*alpha - beta*beta);
return delta*alpha*alpha / (gamma*gamma*gamma);
}
/**
* Computes and returns the standard deviation of the normal inverse gaussian
* distribution with parameters α, β, μ and δ.
*
* @return the standard deviation of the normal inverse gaussian distribution
*
*/
public static double getStandardDeviation (double alpha, double beta,
double mu, double delta) {
return Math.sqrt (getVariance (alpha, beta, mu, delta));
}
/**
* Returns the parameter α of this object.
*
*/
public double getAlpha() {
return alpha;
}
/**
* Returns the parameter β of this object.
*
*/
public double getBeta() {
return beta;
}
/**
* Returns the parameter μ of this object.
*
*/
public double getMu() {
return mu;
}
/**
* Returns the parameter δ of this object.
*
*/
public double getDelta() {
return delta;
}
/**
* Sets the parameters α, β, μ and δ of this object.
*
*/
public void setParams (double alpha, double beta, double mu,
double delta) {
if (delta <= 0.0)
throw new IllegalArgumentException ("delta <= 0");
if (alpha <= 0.0)
throw new IllegalArgumentException ("alpha <= 0");
if (Math.abs(beta) >= alpha)
throw new IllegalArgumentException ("|beta| >= alpha");
gamma = Math.sqrt(alpha*alpha - beta*beta);
this.mu = mu;
this.delta = delta;
this.beta = beta;
this.alpha = alpha;
}
/**
* Returns a table containing the parameters of the current distribution.
* This table is put in regular order: [α, β, μ, δ].
*
*
*/
public double[] getParams () {
double[] retour = {alpha, beta, mu, delta};
return retour;
}
public String toString () {
return getClass().getSimpleName() + ": alpha = " + alpha + ", beta = " + beta +
", mu = " + mu + ", delta = " + delta;
}
}