umontreal.iro.lecuyer.probdistmulti.BiNormalGenzDist 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: BiNormalGenzDist
* Description: bivariate normal distribution using Genz's algorithm
* Environment: Java
* Software: SSJ
* Organization: DIRO, Université de Montréal
* @author
* @since
*/
package umontreal.iro.lecuyer.probdistmulti;
import umontreal.iro.lecuyer.probdist.NormalDist;
/**
* Extends the class {@link BiNormalDist} for the bivariate
* normal distribution
* using Genz's algorithm as described in.
*
*/
public class BiNormalGenzDist extends BiNormalDist {
private static final double[][] W = {
// Gauss Legendre points and weights, n = 6
{ 0.1713244923791705, 0.3607615730481384, 0.4679139345726904},
// Gauss Legendre points and weights, n = 12
{ 0.4717533638651177e-1, 0.1069393259953183, 0.1600783285433464,
0.2031674267230659, 0.2334925365383547, 0.2491470458134029 },
// Gauss Legendre points and weights, n = 20
{ 0.1761400713915212e-1, 0.4060142980038694e-1, 0.6267204833410906e-1,
0.8327674157670475e-1, 0.1019301198172404, 0.1181945319615184,
0.1316886384491766, 0.1420961093183821, 0.1491729864726037,
0.1527533871307259 }
};
private static final double[][] X = {
// Gauss Legendre points and weights, n = 6
{ 0.9324695142031522, 0.6612093864662647, 0.2386191860831970},
// Gauss Legendre points and weights, n = 12
{ 0.9815606342467191, 0.9041172563704750, 0.7699026741943050,
0.5873179542866171, 0.3678314989981802, 0.1252334085114692 },
// Gauss Legendre points and weights, n = 20
{ 0.9931285991850949, 0.9639719272779138, 0.9122344282513259,
0.8391169718222188, 0.7463319064601508, 0.6360536807265150,
0.5108670019508271, 0.3737060887154196, 0.2277858511416451,
0.7652652113349733e-1 }
};
/**
* Constructs a BiNormalGenzDist object with default parameters
*
* μ1 = μ2 = 0,
* σ1 = σ2 = 1 and correlation
* ρ = rho.
*
*/
public BiNormalGenzDist (double rho) {
super (rho);
}
/**
* Constructs a BiNormalGenzDist object with parameters μ1 = mu1,
* μ2 = mu2, σ1 = sigma1, σ2 = sigma2
* and ρ = rho.
*
*/
public BiNormalGenzDist (double mu1, double sigma1,
double mu2, double sigma2, double rho) {
super (mu1, sigma1, mu2, sigma2, rho);
}
/**
* Computes the standard binormal distribution
* with the method described in. The code for the cdf
* was translated directly from the Matlab code written by Alan Genz
* and available from his web page at
* http://www.math.wsu.edu/faculty/genz/homepage (the code is copyrighted by Alan Genz
* and is included in this package with the kind permission of the author).
* The absolute error is expected to be smaller than
* 0.5⋅10-15.
*
*/
public static double cdf (double x, double y, double rho) {
double bvn = specialCDF (x, y, rho, 40.0);
if (bvn >= 0.0)
return bvn;
/*
// Copyright (C) 2005, Alan Genz, All rights reserved.
//
// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided the following conditions are met:
// 1. Redistributions of source code must retain the above copyright
// notice, this list of conditions and the following disclaimer.
// 2. Redistributions in binary form must reproduce the above copyright
// notice, this list of conditions and the following disclaimer in the
// documentation and/or other materials provided with the distribution.
// 3. The contributor name(s) may not be used to endorse or promote
// products derived from this software without specific prior written
// permission.
// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
// "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
// LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
// FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
// COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
// INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
// BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS
// OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
// ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR
// TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE
// USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
//
// function p = bvnl( dh, dk, r )
//
// A function for computing bivariate normal probabilities.
// bvnl calculates the probability that x < dh and y < dk.
// parameters
// dh 1st upper integration limit
// dk 2nd upper integration limit
// r correlation coefficient
//
// Author
// Alan Genz
// Department of Mathematics
// Washington State University
// Pullman, Wa 99164-3113
// Email : [email protected]
// This function is based on the method described by
// Drezner, Z and G.O. Wesolowsky, (1989),
// On the computation of the bivariate normal inegral,
// Journal of Statist. Comput. Simul. 35, pp. 101-107,
// with major modifications for double precision, for |r| close to 1,
// and for matlab by Alan Genz - last modifications 7/98.
//
// p = bvnu( -dh, -dk, r );
// return
//
// end bvnl
//
// function p = bvnu( dh, dk, r )
//
// A function for computing bivariate normal probabilities.
// bvnu calculates the probability that x > dh and y > dk.
// parameters
// dh 1st lower integration limit
// dk 2nd lower integration limit
// r correlation coefficient
//
// Author
// Alan Genz
// Department of Mathematics
// Washington State University
// Pullman, Wa 99164-3113
// Email : [email protected]
//
// This function is based on the method described by
// Drezner, Z and G.O. Wesolowsky, (1989),
// On the computation of the bivariate normal inegral,
// Journal of Statist. Comput. Simul. 35, pp. 101-107,
// with major modifications for double precision, for |r| close to 1,
// and for matlab by Alan Genz - last modifications 7/98.
// Note: to compute the probability that x < dh and y < dk, use
// bvnu( -dh, -dk, r ).
//
*/
final double TWOPI = 2.0 * Math.PI;
final double sqrt2pi = 2.50662827463100050241; // sqrt(TWOPI)
double h, k, hk, hs, asr, sn, as, a, b, c, d, sp, rs, ep, bs, xs;
int i, lg, ng, is;
if (Math.abs (rho) < 0.3) {
ng = 0;
lg = 3;
} else if (Math.abs (rho) < 0.75) {
ng = 1;
lg = 6;
} else {
ng = 2;
lg = 10;
}
h = -x;
k = -y;
hk = h * k;
bvn = 0;
if (Math.abs (rho) < 0.925) {
hs = (h * h + k * k) / 2.0;
asr = Math.asin (rho);
for (i = 0; i < lg; ++i) {
sn = Math.sin (asr * (1.0 - X[ng][i]) / 2.0);
bvn += W[ng][i] * Math.exp ((sn * hk - hs) / (1.0 - sn * sn));
sn = Math.sin (asr * (1.0 + X[ng][i]) / 2.0);
bvn += W[ng][i] * Math.exp ((sn * hk - hs) / (1.0 - sn * sn));
}
bvn = bvn * asr /(4.0*Math.PI) +
NormalDist.cdf01 (-h) * NormalDist.cdf01 (-k);
} else {
if (rho < 0.0) {
k = -k;
hk = -hk;
}
if (Math.abs (rho) < 1.0) {
as = (1.0 - rho) * (1.0 + rho);
a = Math.sqrt (as);
bs = (h - k) * (h - k);
c = (4.0 - hk) / 8.0;
d = (12.0 - hk) / 16.0;
asr = -(bs / as + hk) / 2.0;
if (asr > -100.0)
bvn = a * Math.exp (asr) * (1.0 - c * (bs - as) * (1.0 -
d * bs / 5.0) / 3.0 + c * d * as * as / 5.0);
if (-hk < 100.0) {
b = Math.sqrt (bs);
sp = sqrt2pi * NormalDist.cdf01 (-b / a);
bvn = bvn - Math.exp (-hk / 2.0) * sp * b * (1.0 - c * bs * (1.0 -
d * bs / 5.0) / 3.0);
}
a = a / 2.0;
for (i = 0; i < lg; ++i) {
for (is = -1; is <= 1; is += 2) {
xs = (a * (is * X[ng][i] + 1.0));
xs = xs * xs;
rs = Math.sqrt (1.0 - xs);
asr = -(bs / xs + hk) / 2.0;
if (asr > -100.0) {
sp = (1.0 + c * xs * (1.0 + d * xs));
ep = Math.exp (-hk * (1.0 - rs) / (2.0 * (1.0 + rs))) / rs;
bvn += a * W[ng][i] * Math.exp (asr) * (ep - sp);
}
}
}
bvn = -bvn / TWOPI;
}
if (rho > 0.0) {
if (k > h)
h = k;
bvn += NormalDist.cdf01 (-h);
}
if (rho < 0.0) {
xs = NormalDist.cdf01(-h) - NormalDist.cdf01(-k);
if (xs < 0.0)
xs = 0.0;
bvn = -bvn + xs;
}
}
if (bvn <= 0.0)
return 0.0;
if (bvn >= 1.0)
return 1.0;
return bvn;
}
public static double cdf (double mu1, double sigma1, double x,
double mu2, double sigma2, double y,
double rho) {
if (sigma1 <= 0)
throw new IllegalArgumentException ("sigma1 <= 0");
if (sigma2 <= 0)
throw new IllegalArgumentException ("sigma2 <= 0");
double X = (x - mu1)/sigma1;
double Y = (y - mu2)/sigma2;
return cdf (X, Y, rho);
}
public double cdf (double x, double y) {
return cdf ((x-mu1)/sigma1, (y-mu2)/sigma2, rho);
}
public double barF (double x, double y) {
return barF ((x-mu1)/sigma1, (y-mu2)/sigma2, rho);
}
public static double barF (double mu1, double sigma1, double x,
double mu2, double sigma2, double y,
double rho) {
if (sigma1 <= 0)
throw new IllegalArgumentException ("sigma1 <= 0");
if (sigma2 <= 0)
throw new IllegalArgumentException ("sigma2 <= 0");
double X = (x - mu1)/sigma1;
double Y = (y - mu2)/sigma2;
return barF (X, Y, rho);
}
public static double barF (double x, double y, double rho) {
return cdf (-x, -y, rho);
}
}
© 2015 - 2025 Weber Informatics LLC | Privacy Policy