Module_Prob_Procedures

Copyright 2021 IRD

This file is part of statpack.

statpack is free software: you can redistribute it and/or modify it under the terms of the GNU Lesser General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version.

statpack 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 Lesser General Public License for more details.

You can find a copy of the GNU Lesser General Public License in the statpack/doc directory.

Authors: Pascal Terray (LOCEAN/IPSL, Paris, France)


MODULE EXPORTING SUBROUTINES AND FUNCTIONS FOR PROBABILITY DISTRIBUTION FUNCTIONS, INVERSES, AND OTHER PARAMETERS.

LATEST REVISION : 23/08/2021


function lngamma ( x )

Purpose

Evaluates the logarithm of the gamma function: ln( gamma(X) ) for a strictly positive real argument X.

Arguments

X (INPUT) real(stnd)
On entry, a strictly positive real argument X.

Further Details

This function uses a Lanczos-type approximation to ln(gamma(X)) for X > 0. Its accuracy is about 14 significant digits except for small regions in the vicinity of 1 and 2.

This function is adapted from:

  1. Lanczos, C., 1964:
    A precision approximation of the gamma function. J. SIAM Numer. Anal., B, 1, 86-96.

function lngamma ( x )

Purpose

Evaluates the logarithm of the gamma function: ln( gamma(X(:)) ) for a strictly positive real vector argument X(:).

Arguments

X (INPUT) real(stnd), dimension(:)
On entry, a strictly positive real vector argument X.

Further Details

This function uses a Lanczos-type approximation to ln(gamma(X)) for X > 0. Its accuracy is about 14 significant digits except for small regions in the vicinity of 1 and 2.

The function is parallelized if OPENMP is used.

This function is adapted from:

  1. Lanczos, C., 1964:
    A precision approximation of the gamma function. J. SIAM Numer. Anal., B, 1, 86-96.

function lngamma ( x )

Purpose

Evaluates the logarithm of the gamma function: ln( gamma(X(:,:)) ) for a strictly positive real matrix argument X(:,:).

Arguments

X (INPUT) real(stnd), dimension(:,:)
On entry, a strictly positive real matrix argument X.

Further Details

This function uses a Lanczos-type approximation to ln(gamma(X)) for X > 0. Its accuracy is about 14 significant digits except for small regions in the vicinity of 1 and 2.

The function is parallelized if OPENMP is used.

This function is adapted from

  1. Lanczos, C., 1964:
    A precision approximation of the gamma function. J. SIAM Numer. Anal., B, 1, 86-96.

function probgamma ( x, gamp, acu, maxiter, failure )

Purpose

Evaluates the gamma probability distribution function (e.g. the Incomplete Gamma Integral) for a positive real argument X and a strictly positive value GAMP of the parameter p of the Gamma distribution.

PROBGAMMA computes the probability that a random variable having a Gamma distribution with parameter GAMP will be less than or equal to X.

Arguments

X (INPUT) real(stnd)
On entry, a positive real argument X which is the value of the upper integral limit. X must be greater or equal to zero.
GAMP (INPUT) real(stnd)
On entry, a strictly positive real argument which is the value of the parameter p of the Gamma distribution. GAMP must be greater than zero.
ACU (INPUT, OPTIONAL) real(stnd)

On entry, the desired accuracy of the result. If l decimal places of accuracy are required then ACU should be set to 10**(-(l+1)). ACU is a small strictly positive integer. ACU should not be set smaller than the machine precision since the stated accuracy cannot be attained. In that case the machine precision is used instead.

The default value for ACU is epsilon( ACU ).

MAXITER (INPUT, OPTIONAL) integer(i4b)

On entry, MAXITER controls the maximum number of iterations when evaluating the Pearson’s series expansion of the incomplete Gamma integral.

The default value is 1000.

FAILURE (INPUT, OPTIONAL) logical(lgl)

On entry, if FAILURE is set to true, the values for which the “integrating” process did not converge to the desired accuracy are set to -1. The algorithm fails to converge if the number of iterations exceeds MAXITER.

The default value is false.

Further Details

For large GAMP (e.g. GAMP>1000), this function used a normal approximation, based on the Wilson-Hilferty transformation, see reference (3), Formula 26.4.14, p.941.

Otherwise, a Pearson’s series expansion is used for evaluating the integral, see reference (3), Formula 6.5.29, p.262. The “integrating” process is terminated when both the absolute and relative contributions to the integral is not greater than the value of ACU. The default value for ACU gives the maximum precision of this function.

The time taken by this function thus depends in the precision requested through ACU, and also varies slightly with the input arguments X and GAMP.

This function is more accurate than PROBGAMMA3, but it may be slower.

This function is adapted from:

  1. Lau, C.L., 1980:
    Algorithm AS 147: A simple series for the Incomplete Gamma Integral. Appl. Statist., Vol. 29, No. 1, pp. 113-114
  2. Shea, B.L., 1988:
    Algorithm AS 239: Chi-squared and Incomplete Gamma Integral. Appl. Statist., Vol. 37, No. 3, pp. 466-473
  3. Abramowitz, M., and Stegun, I.A., 1970:
    Handbook of Mathematical Functions, (formulae 6.5.29 and 26.4.14). New York, Dover Publications

function probgamma ( x, gamp, acu, maxiter, failure )

Purpose

Evaluates the gamma probability distribution function (e.g. the Incomplete Gamma Integral) for a positive real vector argument X and a strictly positive value GAMP of the parameter P of the Gamma distribution.

PROBGAMMA computes the probability that a random variable having a Gamma distribution with parameter GAMP will be less than or equal to X(i) for i=1 to size(X).

Arguments

X (INPUT) real(stnd), dimension(:)
On entry, a positive real vector argument X which gives the values of the upper integral limit. Elements in X(:) must be greater or equal to zero.
GAMP (INPUT) real(stnd)
On entry, a strictly positive real argument which is the value of the parameter p of the Gamma distribution. GAMP must be greater than zero.
ACU (INPUT, OPTIONAL) real(stnd)

On entry, the desired accuracy of the result. If l decimal places of accuracy are required then ACU should be set to 10**(-(l+1)). ACU is a small strictly positive integer. ACU should not be set smaller than the machine precision since the stated accuracy cannot be attained. In that case the machine precision is used instead.

The default value for ACU is epsilon( ACU ).

MAXITER (INPUT, OPTIONAL) integer(i4b)

On entry, MAXITER controls the maximum number of iterations when evaluating the Pearson’s series expansion of the incomplete Gamma integral.

The default value is 1000.

FAILURE (INPUT, OPTIONAL) logical(lgl)

On entry, if FAILURE is set to true, the values for which the “integrating” process did not converge to the desired accuracy are set to -1. The algorithm fails to converge if the number of iterations exceeds MAXITER.

The default value is false.

Further Details

For large GAMP (e.g. GAMP>1000), this function used a normal approximation, based on the Wilson-Hilferty transformation, see reference (3), Formula 26.4.14, p.941.

Otherwise, a Pearson’s series expansion is used for evaluating the integral, see reference (3), Formula 6.5.29, p.262. The “integrating” process is terminated when both the absolute and relative contributions to the integral is not greater than the value of ACU. The default value for ACU gives the maximum precision of this function.

The time taken by this function thus depends in the precision requested through ACU, and also varies slightly with the input arguments X and GAMP.

This function is more accurate than PROBGAMMA3, but it may be slower.

The function is parallelized if OPENMP is used.

This function is adapted from:

  1. Lau, C.L., 1980:
    Algorithm AS 147: A simple series for the Incomplete Gamma Integral. Appl. Statist., Vol. 29, No. 1, pp. 113-114
  2. Shea, B.L., 1988:
    Algorithm AS 239: Chi-squared and Incomplete Gamma Integral. Appl. Statist., Vol. 37, No. 3, pp. 466-473
  3. Abramowitz, M., and Stegun, I.A., 1970:
    Handbook of Mathematical Functions, (formulae 6.5.29 and 26.4.14). New York, Dover Publications

function probgamma ( x, gamp, acu, maxiter, failure )

Purpose

Evaluates the gamma probability distribution function (e.g. the Incomplete Gamma Integral) for a positive real vector argument X and a strictly positive vector argument GAMP of the parameter P of the Gamma distribution.

PROBGAMMA computes the probability that a random variable having a Gamma distribution with parameter GAMP(i) will be less than or equal to X(i) for i=1 to size(X).

Arguments

X (INPUT) real(stnd), dimension(:)
On entry, a positive real vector argument X which gives the values of the upper integral limit. Elements in X(:) must be greater or equal to zero.
GAMP (INPUT) real(stnd), dimension(:)

On entry, a strictly positive real vector argument which are the values of the parameter p of the Gamma distribution. Elements in GAMP(:) must be greater than zero.

The size of GAMP must verify size(GAMP) = size(X) .

ACU (INPUT, OPTIONAL) real(stnd)

On entry, the desired accuracy of the result. If l decimal places of accuracy are required then ACU should be set to 10**(-(l+1)). ACU is a small strictly positive integer. ACU should not be set smaller than the machine precision since the stated accuracy cannot be attained. In that case the machine precision is used instead.

The default value for ACU is epsilon( ACU ).

MAXITER (INPUT, OPTIONAL) integer(i4b)

On entry, MAXITER controls the maximum number of iterations when evaluating the Pearson’s series expansion of the incomplete Gamma integral.

The default value is 1000.

FAILURE (INPUT, OPTIONAL) logical(lgl)

On entry, if FAILURE is set to true, the values for which the “integrating” process did not converge to the desired accuracy are set to -1. The algorithm fails to converge if the number of iterations exceeds MAXITER.

The default value is false.

Further Details

For large GAMP(i) (e.g. GAMP(i)>1000), this function used a normal approximation, based on the Wilson-Hilferty transformation, see reference (3), Formula 26.4.14, p.941.

Otherwise, a Pearson’s series expansion is used for evaluating the integral, see reference (3), Formula 6.5.29, p.262. The “integrating” process is terminated when both the absolute and relative contributions to the integral is not greater than the value of ACU. The default value for ACU gives the maximum precision of this function.

The time taken by this function thus depends in the precision requested through ACU, and also varies slightly with the input arguments X and GAMP.

This function is more accurate than PROBGAMMA3, but it may be slower.

The function is parallelized if OPENMP is used.

This function is adapted from:

  1. Lau, C.L., 1980:
    Algorithm AS 147: A simple series for the Incomplete Gamma Integral. Appl. Statist., Vol. 29, No. 1, pp. 113-114
  2. Shea, B.L., 1988:
    Algorithm AS 239: Chi-squared and Incomplete Gamma Integral. Appl. Statist., Vol. 37, No. 3, pp. 466-473
  3. Abramowitz, M., and Stegun, I.A., 1970:
    Handbook of Mathematical Functions, (formulae 6.5.29 and 26.4.14). New York, Dover Publications

function probgamma ( x, gamp, acu, maxiter, failure )

Purpose

Evaluates the gamma probability distribution function (e.g. the Incomplete Gamma Integral) for a positive real matrix argument X and a strictly positive value GAMP of the parameter P of the Gamma distribution.

PROBGAMMA computes the probability that a random variable having a Gamma distribution with parameter GAMP will be less than or equal to X(i,j) for i=1 to size(X,1) and j=1 to size(X,2).

Arguments

X (INPUT) real(stnd), dimension(:,:)
On entry, a positive real matrix argument X which gives the values of the upper integral limit. Elements in X(:,:) must be greater or equal to zero.
GAMP (INPUT) real(stnd)
On entry, a strictly positive real argument which is the value of the parameter p of the Gamma distribution. GAMP must be greater than zero.
ACU (INPUT, OPTIONAL) real(stnd)

On entry, the desired accuracy of the result. If l decimal places of accuracy are required then ACU should be set to 10**(-(l+1)). ACU is a small strictly positive integer. ACU should not be set smaller than the machine precision since the stated accuracy cannot be attained. In that case the machine precision is used instead.

The default value for ACU is epsilon( ACU ).

MAXITER (INPUT, OPTIONAL) integer(i4b)

On entry, MAXITER controls the maximum number of iterations when evaluating the Pearson’s series expansion of the incomplete Gamma integral.

The default value is 1000.

FAILURE (INPUT, OPTIONAL) logical(lgl)

On entry, if FAILURE is set to true, the values for which the “integrating” process did not converge to the desired accuracy are set to -1. The algorithm fails to converge if the number of iterations exceeds MAXITER.

The default value is false.

Further Details

For large GAMP (e.g. GAMP>1000), this function used a normal approximation, based on the Wilson-Hilferty transformation, see reference (3), Formula 26.4.14, p.941.

Otherwise, a Pearson’s series expansion is used for evaluating the integral, see reference (3), Formula 6.5.29, p.262. The “integrating” process is terminated when both the absolute and relative contributions to the integral is not greater than the value of ACU. The default value for ACU gives the maximum precision of this function.

The time taken by this function thus depends in the precision requested through ACU, and also varies slightly with the input arguments X and GAMP.

This function is more accurate than PROBGAMMA3, but it may be slower.

The function is parallelized if OPENMP is used.

This function is adapted from:

  1. Lau, C.L., 1980:
    Algorithm AS 147: A simple series for the Incomplete Gamma Integral. Appl. Statist., Vol. 29, No. 1, pp. 113-114
  2. Shea, B.L., 1988:
    Algorithm AS 239: Chi-squared and Incomplete Gamma Integral. Appl. Statist., Vol. 37, No. 3, pp. 466-473
  3. Abramowitz, M., and Stegun, I.A., 1970:
    Handbook of Mathematical Functions, (formulae 6.5.29 and 26.4.14). New York, Dover Publications

function probgamma ( x, gamp, acu, maxiter, failure )

Purpose

Evaluates the gamma probability distribution function (e.g. the Incomplete Gamma Integral) for a positive real matrix argument X and a strictly positive matrix argument GAMP of the parameter P of the Gamma distribution.

PROBGAMMA computes the probability that a random variable having a Gamma distribution with parameter GAMP(i,j) will be less than or equal to X(i,j) for i=1 to size(X,1) and j=1 to size(X,2).

Arguments

X (INPUT) real(stnd), dimension(:,:)
On entry, a positive real matrix argument X which gives the values of the upper integral limit. Elements in X(:) must be greater or equal to zero.
GAMP (INPUT) real(stnd), dimension(:,:)

On entry, a strictly positive real matrix argument which are the values of the parameter p of the Gamma distribution. Elements in GAMP(:,:) must be greater than zero.

The shape of GAMP must verify:

  • size(GAMP,1) = size(X,1)
  • size(GAMP,2) = size(X,2).
ACU (INPUT, OPTIONAL) real(stnd)

On entry, the desired accuracy of the result. If l decimal places of accuracy are required then ACU should be set to 10**(-(l+1)). ACU is a small strictly positive integer. ACU should not be set smaller than the machine precision since the stated accuracy cannot be attained. In that case the machine precision is used instead.

The default value for ACU is epsilon( ACU ).

MAXITER (INPUT, OPTIONAL) integer(i4b)

On entry, MAXITER controls the maximum number of iterations when evaluating the Pearson’s series expansion of the incomplete Gamma integral.

The default value is 1000.

FAILURE (INPUT, OPTIONAL) logical(lgl)

On entry, if FAILURE is set to true, the values for which the “integrating” process did not converge to the desired accuracy are set to -1. The algorithm fails to converge if the number of iterations exceeds MAXITER.

The default value is false.

Further Details

For large GAMP(i,j) (e.g. GAMP(i,j)>1000), this function used a normal approximation, based on the Wilson-Hilferty transformation, see reference (3), Formula 26.4.14, p.941.

Otherwise, a Pearson’s series expansion is used for evaluating the integral, see reference (3), Formula 6.5.29, p.262. The “integrating” process is terminated when both the absolute and relative contributions to the integral is not greater than the value of ACU. The default value for ACU gives the maximum precision of this function.

The time taken by this function thus depends in the precision requested through ACU, and also varies slightly with the input arguments X and GAMP.

This function is more accurate than PROBGAMMA3, but it may be slower.

The function is parallelized if OPENMP is used.

This function is adapted from:

  1. Lau, C.L., 1980:
    Algorithm AS 147: A simple series for the Incomplete Gamma Integral. Appl. Statist., Vol. 29, No. 1, pp. 113-114
  2. Shea, B.L., 1988:
    Algorithm AS 239: Chi-squared and Incomplete Gamma Integral. Appl. Statist., Vol. 37, No. 3, pp. 466-473
  3. Abramowitz, M., and Stegun, I.A., 1970:
    Handbook of Mathematical Functions, (formulae 6.5.29 and 26.4.14). New York, Dover Publications

function probgamma2 ( x, gamp, acu, maxiter, failure )

Purpose

Evaluates the gamma probability distribution function (e.g. the Incomplete Gamma Integral) for a positive real argument X and a strictly positive value GAMP of the parameter p of the Gamma distribution.

PROBGAMMA2 computes the probability that a random variable having a Gamma distribution with parameter GAMP will be less than or equal to X.

Arguments

X (INPUT) real(stnd)
On entry, a positive real argument X which is the value of the upper integral limit. X must be greater or equal to zero.
GAMP (INPUT) real(stnd)
On entry, a strictly positive real argument which is the value of the parameter p of the Gamma distribution. GAMP must be greater than zero.
ACU (INPUT, OPTIONAL) real(stnd)

On entry, the desired accuracy of the result. If l decimal places of accuracy are required then ACU should be set to 10**(-(l+1)). ACU is a small strictly positive integer. ACU should not be set smaller than the machine precision since the stated accuracy cannot be attained. In that case the machine precision is used instead.

The default value for ACU is epsilon( ACU ).

MAXITER (INPUT, OPTIONAL) integer(i4b)

On entry, MAXITER controls the maximum number of iterations when evaluating the Pearson’s series or continued fraction expansion of the incomplete Gamma integral.

The default value is 1000.

FAILURE (INPUT, OPTIONAL) logical(lgl)

On entry, if FAILURE is set to true, the values for which the “integrating” process did not converge to the desired accuracy are set to -1. The algorithm fails to converge if the number of iterations exceeds MAXITER.

The default value is false.

Further Details

For large GAMP (e.g. GAMP>1000), this function used a normal approximation, based on the Wilson-Hilferty transformation, see reference (1), Formula 26.4.14, p.941.

For X<=1 ou X<GAMP, a Pearson’s series expansion is used, see reference (1), Formula 6.5.29, p.262. For other values of X, a continued fraction expansion is used since this expansion tends to converge more quickly than Pearson’s series expansion, see reference (1), Formula 6.5.31, p.263. In both cases, the “integrating” process is terminated when both the absolute and relative contributions to the integral is not greater than the value of ACU. The default value for ACU gives the maximum precision of this function.

The time taken by this function thus depends in the precision requested through ACU, and also varies slightly with the input arguments X and GAMP.

This function is more accurate than PROBGAMMA3, but is slower.

This function is adapted from:

  1. Abramowitz, M., and Stegun, I.A., 1970:
    Handbook of Mathematical Functions, (formulae 6.5.29, 6.5.31 and 26.4.14). New York, Dover Publications
  2. Shea, B.L., 1988:
    Algorithm AS 239: Chi-squared and Incomplete Gamma Integral. Appl. Statist., Vol. 37, No. 3, pp. 466-473

function probgamma2 ( x, gamp, acu, maxiter, failure )

Purpose

Evaluates the gamma probability distribution function (e.g. the Incomplete Gamma Integral) for a positive real vector argument X and a strictly positive value GAMP of the parameter p of the Gamma distribution.

PROBGAMMA2 computes the probability that a random variable having a Gamma distribution with parameter GAMP will be less than or equal to X(i) for i=1 to size(X).

Arguments

X (INPUT) real(stnd), dimension(:)
On entry, a positive real vector argument X which gives the values of the upper integral limit. Elements in X(:) must be greater or equal to zero.
GAMP (INPUT) real(stnd)
On entry, a strictly positive real argument which is the value of the parameter p of the Gamma distribution. GAMP must be greater than zero.
ACU (INPUT, OPTIONAL) real(stnd)

On entry, the desired accuracy of the result. If l decimal places of accuracy are required then ACU should be set to 10**(-(l+1)). ACU is a small strictly positive integer. ACU should not be set smaller than the machine precision since the stated accuracy cannot be attained. In that case the machine precision is used instead.

The default value for ACU is epsilon( ACU ).

MAXITER (INPUT, OPTIONAL) integer(i4b)

On entry, MAXITER controls the maximum number of iterations when evaluating the Pearson’s series or continued fraction expansion of the incomplete Gamma integral.

The default value is 1000.

FAILURE (INPUT, OPTIONAL) logical(lgl)

On entry, if FAILURE is set to true, the values for which the “integrating” process did not converge to the desired accuracy are set to -1. The algorithm fails to converge if the number of iterations exceeds MAXITER.

The default value is false.

Further Details

For large GAMP (e.g. GAMP>1000), this function used a normal approximation, based on the Wilson-Hilferty transformation, see reference (1), Formula 26.4.14, p.941.

For X<=1 ou X<GAMP, a Pearson’s series expansion is used, see reference (1), Formula 6.5.29, p.262. For other values of X, a continued fraction expansion is used since this expansion tends to converge more quickly than Pearson’s series expansion, see reference (1), Formula 6.5.31, p.263. In both cases, the “integrating” process is terminated when both the absolute and relative contributions to the integral is not greater than the value of ACU. The default value for ACU gives the maximum precision of this function.

The time taken by this function thus depends in the precision requested through ACU, and also varies slightly with the input arguments X and GAMP.

This function is more accurate than PROBGAMMA3, but is slower.

This function is adapted from:

  1. Abramowitz, M., and Stegun, I.A., 1970:
    Handbook of Mathematical Functions, (formulae 6.5.29, 6.5.31 and 26.4.14). New York, Dover Publications
  2. Shea, B.L., 1988:
    Algorithm AS 239: Chi-squared and Incomplete Gamma Integral. Appl. Statist., Vol. 37, No. 3, pp. 466-473

function probgamma2 ( x, gamp, acu, maxiter, failure )

Purpose

Evaluates the gamma probability distribution function (e.g. the Incomplete Gamma Integral) for a positive real vector argument X and a strictly positive vector argument GAMP of the parameter P of the Gamma distribution.

PROBGAMMA2 computes the probability that a random variable having a Gamma distribution with parameter GAMP(i) will be less than or equal to X(i) for i=1 to size(X).

Arguments

X (INPUT) real(stnd), dimension(:)
On entry, a positive real vector argument X which gives the values of the upper integral limit. Elements in X(:) must be greater or equal to zero.
GAMP (INPUT) real(stnd), dimension(:)

On entry, a strictly positive real vector argument which are the values of the parameter p of the Gamma distribution. Elements in GAMP(:) must be greater than zero.

The size of GAMP must verify size(GAMP) = size(X) .

ACU (INPUT, OPTIONAL) real(stnd)

On entry, the desired accuracy of the result. If l decimal places of accuracy are required then ACU should be set to 10**(-(l+1)). ACU is a small strictly positive integer. ACU should not be set smaller than the machine precision since the stated accuracy cannot be attained. In that case the machine precision is used instead.

The default value for ACU is epsilon( ACU ).

MAXITER (INPUT, OPTIONAL) integer(i4b)

On entry, MAXITER controls the maximum number of iterations when evaluating the Pearson’s series or continued fraction expansion of the incomplete Gamma integral.

The default value is 1000.

FAILURE (INPUT, OPTIONAL) logical(lgl)

On entry, if FAILURE is set to true, the values for which the “integrating” process did not converge to the desired accuracy are set to -1. The algorithm fails to converge if the number of iterations exceeds MAXITER.

The default value is false.

Further Details

For large GAMP(i) (e.g. GAMP(i)>1000), this function used a normal approximation, based on the Wilson-Hilferty transformation, see reference (1), Formula 26.4.14, p.941.

For X(i)<=1 ou X(i)<GAMP(i), a Pearson’s series expansion is used, see reference (1), Formula 6.5.29, p.262. For other values of X(i), a continued fraction expansion is used since this expansion tends to converge more quickly than Pearson’s series expansion, see reference (1), Formula 6.5.31, p.263. In both cases, the “integrating” process is terminated when both the absolute and relative contributions to the integral is not greater than the value of ACU. The default value for ACU gives the maximum precision of this function.

The time taken by this function thus depends in the precision requested through ACU, and also varies slightly with the input arguments X and GAMP.

This function is more accurate than PROBGAMMA3, but is slower.

This function is adapted from:

  1. Abramowitz, M., and Stegun, I.A., 1970:
    Handbook of Mathematical Functions, (formulae 6.5.29, 6.5.31 and 26.4.14). New York, Dover Publications
  2. Shea, B.L., 1988:
    Algorithm AS 239: Chi-squared and Incomplete Gamma Integral. Appl. Statist., Vol. 37, No. 3, pp. 466-473

function probgamma2 ( x, gamp, acu, maxiter, failure )

Purpose

Evaluates the gamma probability distribution function (e.g. the Incomplete Gamma Integral) for a positive real matrix argument X and a strictly positive value GAMP of the parameter p of the Gamma distribution.

PROBGAMMA2 computes the probability that a random variable having a Gamma distribution with parameter GAMP will be less than or equal to X(i,j) for i=1 to size(X,1) and j=1 to size(X,2).

Arguments

X (INPUT) real(stnd), dimension(:,:)
On entry, a positive real matrix argument X which gives the values of the upper integral limit. Elements in X(:,:) must be greater or equal to zero.
GAMP (INPUT) real(stnd)
On entry, a strictly positive real argument which is the value of the parameter p of the Gamma distribution. GAMP must be greater than zero.
ACU (INPUT, OPTIONAL) real(stnd)

On entry, the desired accuracy of the result. If l decimal places of accuracy are required then ACU should be set to 10**(-(l+1)). ACU is a small strictly positive integer. ACU should not be set smaller than the machine precision since the stated accuracy cannot be attained. In that case the machine precision is used instead.

The default value for ACU is epsilon( ACU ).

MAXITER (INPUT, OPTIONAL) integer(i4b)

On entry, MAXITER controls the maximum number of iterations when evaluating the Pearson’s series or continued fraction expansion of the incomplete Gamma integral.

The default value is 1000.

FAILURE (INPUT, OPTIONAL) logical(lgl)

On entry, if FAILURE is set to true, the values for which the “integrating” process did not converge to the desired accuracy are set to -1. The algorithm fails to converge if the number of iterations exceeds MAXITER.

The default value is false.

Further Details

For large GAMP (e.g. GAMP>1000), this function used a normal approximation, based on the Wilson-Hilferty transformation, see reference (1), Formula 26.4.14, p.941.

For X<=1 ou X<GAMP, a Pearson’s series expansion is used, see reference (1), Formula 6.5.29, p.262. For other values of X, a continued fraction expansion is used since this expansion tends to converge more quickly than Pearson’s series expansion, see reference (1), Formula 6.5.31, p.263. In both cases, the “integrating” process is terminated when both the absolute and relative contributions to the integral is not greater than the value of ACU. The default value for ACU gives the maximum precision of this function.

The time taken by this function thus depends in the precision requested through ACU, and also varies slightly with the input arguments X and GAMP.

This function is more accurate than PROBGAMMA3, but is slower.

The function is parallelized if OPENMP is used.

This function is adapted from:

  1. Abramowitz, M., and Stegun, I.A., 1970:
    Handbook of Mathematical Functions, (formulae 6.5.29, 6.5.31 and 26.4.14). New York, Dover Publications
  2. Shea, B.L., 1988:
    Algorithm AS 239: Chi-squared and Incomplete Gamma Integral. Appl. Statist., Vol. 37, No. 3, pp. 466-473

function probgamma2 ( x, gamp, acu, maxiter, failure )

Purpose

Evaluates the gamma probability distribution function (e.g. the Incomplete Gamma Integral) for a positive real matrix argument X and a strictly positive matrix argument GAMP of the parameter P of the Gamma distribution.

PROBGAMMA2 computes the probability that a random variable having a Gamma distribution with parameter GAMP(i,j) will be less than or equal to X(i,j) for i=1 to size(X,1) and j=1 to size(X,2).

Arguments

X (INPUT) real(stnd), dimension(:,:)
On entry, a positive real matrix argument X which gives the values of the upper integral limit. Elements in X(:,:) must be greater or equal to zero.
GAMP (INPUT) real(stnd), dimension(:,:)

On entry, a strictly positive real matrix argument which are the values of the parameter p of the Gamma distribution. Elements in GAMP(:,:) must be greater than zero.

The shape of GAMP must verify:

  • size(GAMP,1) = size(X,1)
  • size(GAMP,2) = size(X,2).
ACU (INPUT, OPTIONAL) real(stnd)

On entry, the desired accuracy of the result. If l decimal places of accuracy are required then ACU should be set to 10**(-(l+1)). ACU is a small strictly positive integer. ACU should not be set smaller than the machine precision since the stated accuracy cannot be attained. In that case the machine precision is used instead.

The default value for ACU is epsilon( ACU ).

MAXITER (INPUT, OPTIONAL) integer(i4b)

On entry, MAXITER controls the maximum number of iterations when evaluating the Pearson’s series or continued fraction expansion of the incomplete Gamma integral.

The default value is 1000.

FAILURE (INPUT, OPTIONAL) logical(lgl)

On entry, if FAILURE is set to true, the values for which the “integrating” process did not converge to the desired accuracy are set to -1. The algorithm fails to converge if the number of iterations exceeds MAXITER.

The default value is false.

Further Details

For large GAMP(i,j) (e.g. GAMP(i,j)>1000), this function used a normal approximation, based on the Wilson-Hilferty transformation, see reference (1), Formula 26.4.14, p.941.

For X(i,j)<=1 ou X(i,j)<GAMP(i,j), a Pearson’s series expansion is used, see reference (1), Formula 6.5.29, p.262. For other values of X(i,j), a continued fraction expansion is used since this expansion tends to converge more quickly than Pearson’s series expansion, see reference (1), Formula 6.5.31, p.263. In both cases, the “integrating” process is terminated when both the absolute and relative contributions to the integral is not greater than the value of ACU. The default value for ACU gives the maximum precision of this function.

The time taken by this function thus depends in the precision requested through ACU, and also varies slightly with the input arguments X and GAMP.

This function is more accurate than PROBGAMMA3, but is slower.

The function is parallelized if OPENMP is used.

This function is adapted from:

  1. Abramowitz, M., and Stegun, I.A., 1970:
    Handbook of Mathematical Functions, (formulae 6.5.29, 6.5.31 and 26.4.14). New York, Dover Publications
  2. Shea, B.L., 1988:
    Algorithm AS 239: Chi-squared and Incomplete Gamma Integral. Appl. Statist., Vol. 37, No. 3, pp. 466-473

function probgamma3 ( x, gamp, acu, maxiter, failure )

Purpose

Evaluates the gamma probability distribution function (e.g. the Incomplete Gamma Integral) for a positive real argument X and a strictly positive value GAMP of the parameter p of the Gamma distribution.

PROBGAMMA3 computes the probability that a random variable having a Gamma distribution with parameter GAMP will be less than or equal to X.

Arguments

X (INPUT) real(stnd)
On entry, a positive real argument X which is the value of the upper integral limit. X must be greater or equal to zero.
GAMP (INPUT) real(stnd)
On entry, a strictly positive real argument which is the value of the parameter p of the Gamma distribution. GAMP must be greater than zero.
ACU (INPUT, OPTIONAL) real(stnd)

On entry, the desired accuracy of the result. If l decimal places of accuracy are required then ACU should be set to 10**(-(l+1)). ACU is a small strictly positive integer. ACU should not be set smaller than the machine precision since the stated accuracy cannot be attained. In that case the machine precision is used instead.

The default value for ACU is epsilon( ACU ).

MAXITER (INPUT, OPTIONAL) integer(i4b)

On entry, MAXITER controls the maximum number of iterations when evaluating the Pearson’s series expansion of the incomplete Gamma integral.

The default value is 1000.

FAILURE (INPUT, OPTIONAL) logical(lgl)

On entry, if FAILURE is set to true, the values for which the “integrating” process did not converge to the desired accuracy are set to -1. The algorithm fails to converge if the number of iterations exceeds MAXITER.

The default value is false.

Further Details

For large GAMP (e.g. GAMP>1000), this function used a normal approximation, based on the Wilson-Hilferty transformation, see reference (3), Formula 26.4.14, p.941.

For X<=max(GAMP/2,13), a Pearson’s series expansion is used, see reference (3), Formula 6.5.29, p.262. For larger values of X, an alternate Pearson’s asymptotic series expansion is used since this expansion tends to converge more quickly, see reference (2), Formula 6.5.32, p.263.

In both cases, the “integrating” process is terminated when both the absolute and relative contributions to the integral is not greater than the value of ACU. The default value for ACU gives the maximum precision of this function.

The time taken by this function thus depends in the precision requested through ACU, and also varies slightly with the input arguments X and GAMP.

PROBGAMMA3 is faster, but less accurate than PROBGAMMA or PROBGAMMA2, since for large values of X, the alternate Pearson’s series expansion is only asymptotic.

This function is adapted from:

  1. Lau, C.L., 1980:
    Algorithm AS 147: A simple series for the Incomplete Gamma Integral. Appl. Statist., Vol. 29, No. 1, pp. 113-114
  2. Shea, B.L., 1988:
    Algorithm AS 239: Chi-squared and Incomplete Gamma Integral. Appl. Statist., Vol. 37, No. 3, pp. 466-473
  3. Abramowitz, M., and Stegun, I.A., 1970:
    Handbook of Mathematical Functions (formulae 6.5.29, 6.5.32 and 26.4.14). New York, Dover Publications

function probgamma3 ( x, gamp, acu, maxiter, failure )

Purpose

Evaluates the gamma probability distribution function (e.g. the Incomplete Gamma Integral) for a positive real vector argument X and a strictly positive value GAMP of the parameter P of the Gamma distribution.

PROBGAMMA3 computes the probability that a random variable having a Gamma distribution with parameter GAMP will be less than or equal to X(i) for i=1 to size(X).

Arguments

X (INPUT) real(stnd), dimension(:)
On entry, a positive real vector argument X which gives the values of the upper integral limit. Elements in X(:) must be greater or equal to zero.
GAMP (INPUT) real(stnd)
On entry, a strictly positive real argument which is the value of the parameter p of the Gamma distribution. GAMP must be greater than zero.
ACU (INPUT, OPTIONAL) real(stnd)

On entry, the desired accuracy of the result. If l decimal places of accuracy are required then ACU should be set to 10**(-(l+1)). ACU is a small strictly positive integer. ACU should not be set smaller than the machine precision since the stated accuracy cannot be attained. In that case the machine precision is used instead.

The default value for ACU is epsilon( ACU ).

MAXITER (INPUT, OPTIONAL) integer(i4b)

On entry, MAXITER controls the maximum number of iterations when evaluating the Pearson’s series expansion of the incomplete Gamma integral.

The default value is 1000.

FAILURE (INPUT, OPTIONAL) logical(lgl)

On entry, if FAILURE is set to true, the values for which the “integrating” process did not converge to the desired accuracy are set to -1. The algorithm fails to converge if the number of iterations exceeds MAXITER.

The default value is false.

Further Details

For large GAMP (e.g. GAMP>1000), this function used a normal approximation, based on the Wilson-Hilferty transformation, see reference (3), Formula 26.4.14, p.941.

For X<=max(GAMP/2,13), a Pearson’s series expansion is used, see reference (3), Formula 6.5.29, p.262. For larger values of X, an alternate Pearson’s asymptotic series expansion is used since this expansion tends to converge more quickly, see reference (2), Formula 6.5.32, p.263.

In both cases, the “integrating” process is terminated when both the absolute and relative contributions to the integral is not greater than the value of ACU. The default value for ACU gives the maximum precision of this function.

The time taken by this function thus depends in the precision requested through ACU, and also varies slightly with the input arguments X and GAMP.

PROBGAMMA3 is faster, but less accurate than PROBGAMMA or PROBGAMMA2, since for large values of X, the alternate Pearson’s series expansion is only asymptotic.

This function is adapted from:

  1. Lau, C.L., 1980:
    Algorithm AS 147: A simple series for the Incomplete Gamma Integral. Appl. Statist., Vol. 29, No. 1, pp. 113-114
  2. Shea, B.L., 1988:
    Algorithm AS 239: Chi-squared and Incomplete Gamma Integral. Appl. Statist., Vol. 37, No. 3, pp. 466-473
  3. Abramowitz, M., and Stegun, I.A., 1970:
    Handbook of Mathematical Functions (formulae 6.5.29, 6.5.32 and 26.4.14). New York, Dover Publications

function probgamma3 ( x, gamp, acu, maxiter, failure )

Purpose

Evaluates the gamma probability distribution function (e.g. the Incomplete Gamma Integral) for a positive real vector argument X and a strictly positive vector argument GAMP of the parameter P of the Gamma distribution.

PROBGAMMA3 computes the probability that a random variable having a Gamma distribution with parameter GAMP(i) will be less than or equal to X(i) for i=1 to size(X).

Arguments

X (INPUT) real(stnd), dimension(:)
On entry, a positive real vector argument X which gives the values of the upper integral limit. Elements in X(:) must be greater or equal to zero.
GAMP (INPUT) real(stnd), dimension(:)

On entry, a strictly positive real vector argument which are the values of the parameter p of the Gamma distribution. Elements in GAMP(:) must be greater than zero.

The size of GAMP must verify size(GAMP) = size(X) .

ACU (INPUT, OPTIONAL) real(stnd)

On entry, the desired accuracy of the result. If l decimal places of accuracy are required then ACU should be set to 10**(-(l+1)). ACU is a small strictly positive integer. ACU should not be set smaller than the machine precision since the stated accuracy cannot be attained. In that case the machine precision is used instead.

The default value for ACU is epsilon( ACU ).

MAXITER (INPUT, OPTIONAL) integer(i4b)

On entry, MAXITER controls the maximum number of iterations when evaluating the Pearson’s series expansion of the incomplete Gamma integral.

The default value is 1000.

FAILURE (INPUT, OPTIONAL) logical(lgl)

On entry, if FAILURE is set to true, the values for which the “integrating” process did not converge to the desired accuracy are set to -1. The algorithm fails to converge if the number of iterations exceeds MAXITER.

The default value is false.

Further Details

For large GAMP(i) (e.g. GAMP(i)>1000), this function used a normal approximation, based on the Wilson-Hilferty transformation, see reference (3), Formula 26.4.14, p.941.

For X(i)<=max(GAMP(i)/2,13), a Pearson’s series expansion is used, see reference (3), Formula 6.5.29, p.262. For larger values of X(i), an alternate Pearson’s asymptotic series expansion is used since this expansion tends to converge more quickly, see reference (2), Formula 6.5.32, p.263.

In both cases, the “integrating” process is terminated when both the absolute and relative contributions to the integral is not greater than the value of ACU. The default value for ACU gives the maximum precision of this function.

The time taken by this function thus depends in the precision requested through ACU, and also varies slightly with the input arguments X and GAMP.

PROBGAMMA3 is faster, but less accurate than PROBGAMMA or PROBGAMMA2, since for large values of X, the alternate Pearson’s series expansion is only asymptotic.

This function is adapted from:

  1. Lau, C.L., 1980:
    Algorithm AS 147: A simple series for the Incomplete Gamma Integral. Appl. Statist., Vol. 29, No. 1, pp. 113-114
  2. Shea, B.L., 1988:
    Algorithm AS 239: Chi-squared and Incomplete Gamma Integral. Appl. Statist., Vol. 37, No. 3, pp. 466-473
  3. Abramowitz, M., and Stegun, I.A., 1970:
    Handbook of Mathematical Functions (formulae 6.5.29, 6.5.32 and 26.4.14). New York, Dover Publications

function probgamma3 ( x, gamp, acu, maxiter, failure )

Purpose

Evaluates the gamma probability distribution function (e.g. the Incomplete Gamma Integral) for a positive real matrix argument X and a strictly positive value GAMP of the parameter P of the Gamma distribution.

PROBGAMMA3 computes the probability that a random variable having a Gamma distribution with parameter GAMP will be less than or equal to X(i,j) for i=1 to size(X,1) and j=1 to size(X,2).

Arguments

X (INPUT) real(stnd), dimension(:,:)
On entry, a positive real matrix argument X which gives the values of the upper integral limit. Elements in X(:,:) must be greater or equal to zero.
GAMP (INPUT) real(stnd)
On entry, a strictly positive real argument which is the value of the parameter p of the Gamma distribution. GAMP must be greater than zero.
ACU (INPUT, OPTIONAL) real(stnd)

On entry, the desired accuracy of the result. If l decimal places of accuracy are required then ACU should be set to 10**(-(l+1)). ACU is a small strictly positive integer. ACU should not be set smaller than the machine precision since the stated accuracy cannot be attained. In that case the machine precision is used instead.

The default value for ACU is epsilon( ACU ).

MAXITER (INPUT, OPTIONAL) integer(i4b)

On entry, MAXITER controls the maximum number of iterations when evaluating the Pearson’s series expansion of the incomplete Gamma integral.

The default value is 1000.

FAILURE (INPUT, OPTIONAL) logical(lgl)

On entry, if FAILURE is set to true, the values for which the “integrating” process did not converge to the desired accuracy are set to -1. The algorithm fails to converge if the number of iterations exceeds MAXITER.

The default value is false.

Further Details

For large GAMP (e.g. GAMP>1000), this function used a normal approximation, based on the Wilson-Hilferty transformation, see reference (3), Formula 26.4.14, p.941.

For X<=max(GAMP/2,13), a Pearson’s series expansion is used, see reference (3), Formula 6.5.29, p.262. For larger values of X, an alternate Pearson’s asymptotic series expansion is used since this expansion tends to converge more quickly, see reference (2), Formula 6.5.32, p.263.

In both cases, the “integrating” process is terminated when both the absolute and relative contributions to the integral is not greater than the value of ACU. The default value for ACU gives the maximum precision of this function.

The time taken by this function thus depends in the precision requested through ACU, and also varies slightly with the input arguments X and GAMP.

PROBGAMMA3 is faster, but less accurate than PROBGAMMA or PROBGAMMA2, since for large values of X, the alternate Pearson’s series expansion is only asymptotic.

The function is parallelized if OPENMP is used.

This function is adapted from:

  1. Lau, C.L., 1980:
    Algorithm AS 147: A simple series for the Incomplete Gamma Integral. Appl. Statist., Vol. 29, No. 1, pp. 113-114
  2. Shea, B.L., 1988:
    Algorithm AS 239: Chi-squared and Incomplete Gamma Integral. Appl. Statist., Vol. 37, No. 3, pp. 466-473
  3. Abramowitz, M., and Stegun, I.A., 1970:
    Handbook of Mathematical Functions (formulae 6.5.29, 6.5.32 and 26.4.14). New York, Dover Publications

function probgamma3 ( x, gamp, acu, maxiter, failure )

Purpose

Evaluates the gamma probability distribution function (e.g. the Incomplete Gamma Integral) for a positive real matrix argument X and a strictly positive matrix argument GAMP of the parameter P of the Gamma distribution.

PROBGAMMA3 computes the probability that a random variable having a Gamma distribution with parameter GAMP(i,j) will be less than or equal to X(i,j) for i=1 to size(X,1) and j=1 to size(X,2).

Arguments

X (INPUT) real(stnd), dimension(:,:)
On entry, a positive real matrix argument X which gives the values of the upper integral limit. Elements in X(:,:) must be greater or equal to zero.
GAMP (INPUT) real(stnd), dimension(:,:)

On entry, a strictly positive real matrix argument which are the values of the parameter p of the Gamma distribution. Elements in GAMP(:,:) must be greater than zero.

The shape of GAMP must verify:

  • size(GAMP,1) = size(X,1)
  • size(GAMP,2) = size(X,2).
ACU (INPUT, OPTIONAL) real(stnd)

On entry, the desired accuracy of the result. If l decimal places of accuracy are required then ACU should be set to 10**(-(l+1)). ACU is a small strictly positive integer. ACU should not be set smaller than the machine precision since the stated accuracy cannot be attained. In that case the machine precision is used instead.

The default value for ACU is epsilon( ACU ).

MAXITER (INPUT, OPTIONAL) integer(i4b)

On entry, MAXITER controls the maximum number of iterations when evaluating the Pearson’s series expansion of the incomplete Gamma integral.

The default value is 1000.

FAILURE (INPUT, OPTIONAL) logical(lgl)

On entry, if FAILURE is set to true, the values for which the “integrating” process did not converge to the desired accuracy are set to -1. The algorithm fails to converge if the number of iterations exceeds MAXITER.

The default value is false.

Further Details

For large GAMP(i,j) (e.g. GAMP(i,j)>1000), this function used a normal approximation, based on the Wilson-Hilferty transformation, see reference (3), Formula 26.4.14, p.941.

For X(i,j)<=max(GAMP(i,j)/2,13), a Pearson’s series expansion is used, see reference (3), Formula 6.5.29, p.262. For larger values of X(i,j), an alternate Pearson’s asymptotic series expansion is used since this expansion tends to converge more quickly, see reference (2), Formula 6.5.32, p.263.

In both cases, the “integrating” process is terminated when both the absolute and relative contributions to the integral is not greater than the value of ACU. The default value for ACU gives the maximum precision of this function.

The time taken by this function thus depends in the precision requested through ACU, and also varies slightly with the input arguments X and GAMP.

PROBGAMMA3 is faster, but less accurate than PROBGAMMA or PROBGAMMA2, since for large values of X, the alternate Pearson’s series expansion is only asymptotic.

The function is parallelized if OPENMP is used.

This function is adapted from:

  1. Lau, C.L., 1980:
    Algorithm AS 147: A simple series for the Incomplete Gamma Integral. Appl. Statist., Vol. 29, No. 1, pp. 113-114
  2. Shea, B.L., 1988:
    Algorithm AS 239: Chi-squared and Incomplete Gamma Integral. Appl. Statist., Vol. 37, No. 3, pp. 466-473
  3. Abramowitz, M., and Stegun, I.A., 1970:
    Handbook of Mathematical Functions (formulae 6.5.29, 6.5.32 and 26.4.14). New York, Dover Publications

function pinvgamma ( p, gamp, acu, maxiter )

Purpose

Evaluates the inverse gamma probability distribution function.

For given arguments P (0<=P<=1) and GAMP (GAMP>0), PINVGAMMA returns the value X such that P is the probability that a random variable distributed as a gamma distribution with parameter GAMP is less than or equal to X.

Arguments

P (INPUT) real(stnd)
On entry, input probability. P must be in the range (0,1) inclusive.
GAMP (INPUT) real(stnd)
on entry, the parameter p of the gamma distribution. GAMP must be greater or equal to 0.25 .
ACU (INPUT, OPTIONAL) real(stnd)

On entry, the desired accuracy of the result when computing the incomplete Gamma integral in the evaluation of the seven term Taylor series in function PINVQ2 .

If l decimal places of accuracy are required then ACU should be set to 10**(-(l+1)). ACU is a small strictly positive integer. ACU should not be set smaller than the machine precision since the stated accuracy cannot be attained. In that case the machine precision is used instead.

The default value for ACU is epsilon( ACU ).

MAXITER (INPUT, OPTIONAL) integer(i4b)

On entry, MAXITER controls the maximum number of iterations when evaluating the Pearson’s series or continued fraction expansion of the incomplete Gamma integral.

The default value is 1000.

See the description of the PROBGAMMA2 function for more details on this argument.

Further Details

This function actually uses PINVQ2 function and is adapted from:

  1. Best, D.J., and Roberts, D.E., 1975:
    Algorithm AS 91: The Percentage Points of the chi2 Distribution. Appl. Statist., Vol.24, No. 3, pp.385-388
  2. Shea, B.L., 1988:
    Algorithm AS 239: Chi-squared and Incomplete Gamma Integral. Appl. Statist., Vol. 37, No. 3, pp. 466-473
  3. Shea, B.L., 1991:
    Algorithm AS R85: A remark on AS 91: The Percentage Points of the chi2 Distribution. Appl. Statist. Vol.40, No. 1, pp.233-235.

function probbeta ( x, a, b, beta, acu, maxiter, failure )

Purpose

Evaluates the beta probability distribution function (e.g the incomplete beta function).

For given arguments X (0<=X<=1), A (A>0), B (B>0), PROBBETA returns the probability that a random variable from a beta distribution having parameters A and B will be less than or equal to X.

Arguments

X (INPUT) real(stnd)
On entry, the value to which function is to be integrated. X must be in range (0,1) inclusive.
A (INPUT) real(stnd)
on entry, the (1st) parameter of the beta distribution. A must be greater than 0.
B (INPUT) real(stnd)
On entry, the (2nd) parameter of the beta distribution. B must be greater than 0.
BETA (INPUT, OPTIONAL) real(stnd)

On entry, the logarithm of the complete beta function BETA(A,B).

If BETA is not given, the logarithm of the beta function is computed with the help of function LNGAMMA.

ACU (INPUT, OPTIONAL) real(stnd)

On entry, the desired accuracy of the result. The “integrating” process is terminated when both the absolute and relative contributions to the integral is not greater than the value of ACU.

ACU is a small strictly positive integer. If the number of decimal digits’ accuracy required is r, ACU should be set to 10**(-(r+1)).

The default value for ACU is epsilon( ACU ).

MAXITER (INPUT, OPTIONAL) integer(i4b)

On entry, MAXITER controls the maximum number of iterations when evaluating the power series representation of the incomplete Beta integral.

The default value is 2000.

FAILURE (INPUT, OPTIONAL) logical(lgl)

On entry, if FAILURE is set to true, the values for which the “integrating” process did not converge to the desired accuracy are set to -1. The algorithm fails to converge if the number of iterations exceeds MAXITER.

The default value is false.

Further Details

This function is adapted from:

  1. Majumder, K.L., and Bhattacharjee, G.P., 1973:
    Algorithm AS 63: the Incomplete Beta Integral. Appl. Statist., Vol.22, No.3, pp 409-411.
  2. Cran, G.W., Martin, K.J., and Thomas, G.E., 1977:
    Remark AS R19 and Algorithm AS 109: A remark on Algorithms: AS 63 the Incomplete Beta integral, AS 64 Inverse of the Incomplete Beta Function Ratio. Appl. Statist., Vol.26, No.1, pp 111-114.

function probbeta ( x, a, b, beta, acu, maxiter, failure )

Purpose

Evaluates the beta probability distribution function (e.g the incomplete beta function).

For given arguments X(:) (0<=X(:)<=1), A (A>0), B (B>0), PROBBETA returns the probabilities that a random variable from a beta distribution having parameters A and B will be less than or equal to X(:).

Arguments

X (INPUT) real(stnd), dimension(:)
On entry, the values to which function is to be integrated. Elements in X(:) must be in the range (0,1) inclusive.
A (INPUT) real(stnd)
on entry, the (1st) parameter of the beta distribution. A must be greater than 0.
B (INPUT) real(stnd)
On entry, the (2nd) parameter of the beta distribution. B must be greater than 0.
BETA (INPUT, OPTIONAL) real(stnd)

On entry, the logarithm of the complete beta function BETA(A,B).

If BETA is not given, the logarithm of the beta function is computed with the help of function LNGAMMA.

ACU (INPUT, OPTIONAL) real(stnd)

On entry, the desired accuracy of the result. The “integrating” process is terminated when both the absolute and relative contributions to the integral is not greater than the value of ACU.

ACU is a small strictly positive integer. If the number of decimal digits’ accuracy required is r, ACU should be set to 10**(-(r+1)).

The default value for ACU is epsilon( ACU ).

MAXITER (INPUT, OPTIONAL) integer(i4b)

On entry, MAXITER controls the maximum number of iterations when evaluating the power series representation of the incomplete Beta integral.

The default value is 2000.

FAILURE (INPUT, OPTIONAL) logical(lgl)

On entry, if FAILURE is set to true, the values for which the “integrating” process did not converge to the desired accuracy are set to -1. The algorithm fails to converge if the number of iterations exceeds MAXITER.

The default value is false.

Further Details

This function is parallelized if OPENMP is used.

This function is adapted from:

  1. Majumder, K.L., and Bhattacharjee, G.P., 1973:
    Algorithm AS 63: the Incomplete Beta Integral. Appl. Statist., Vol.22, No.3, pp 409-411.
  2. Cran, G.W., Martin, K.J., and Thomas, G.E., 1977:
    Remark AS R19 and Algorithm AS 109: A remark on Algorithms: AS 63 the Incomplete Beta integral, AS 64 Inverse of the Incomplete Beta Function Ratio. Appl. Statist., Vol.26, No.1, pp 111-114.

function probbeta ( x, a, b, beta, acu, maxiter, failure )

Purpose

Evaluates the beta probability distribution function (e.g the incomplete beta function).

For given arguments X(:,:) (0<=X(:,:)<=1), A (A>0), B (B>0), PROBBETA returns the probabilities that a random variable from a beta distribution having parameters A and B will be less than or equal to X(:,:).

Arguments

X (INPUT) real(stnd), dimension(:,:)
On entry, the values to which function is to be integrated. Elements in X(:,:) must be in the range (0,1) inclusive.
A (INPUT) real(stnd)
on entry, the (1st) parameter of the beta distribution. A must be greater than 0.
B (INPUT) real(stnd)
On entry, the (2nd) parameter of the beta distribution. B must be greater than 0.
BETA (INPUT, OPTIONAL) real(stnd)

On entry, the logarithm of the complete beta function BETA(A,B).

If BETA is not given, the logarithm of the beta function is computed with the help of function LNGAMMA.

ACU (INPUT, OPTIONAL) real(stnd)

On entry, the desired accuracy of the result. The “integrating” process is terminated when both the absolute and relative contributions to the integral is not greater than the value of ACU.

ACU is a small strictly positive integer. If the number of decimal digits’ accuracy required is r, ACU should be set to 10**(-(r+1)).

The default value for ACU is epsilon( ACU ).

MAXITER (INPUT, OPTIONAL) integer(i4b)

On entry, MAXITER controls the maximum number of iterations when evaluating the power series representation of the incomplete Beta integral.

The default value is 2000.

FAILURE (INPUT, OPTIONAL) logical(lgl)

On entry, if FAILURE is set to true, the values for which the “integrating” process did not converge to the desired accuracy are set to -1. The algorithm fails to converge if the number of iterations exceeds MAXITER.

The default value is false.

Further Details

This function is parallelized if OPENMP is used.

This function is adapted from:

  1. Majumder, K.L., and Bhattacharjee, G.P., 1973:
    Algorithm AS 63: the Incomplete Beta Integral. Appl. Statist., Vol.22, No.3, pp 409-411.
  2. Cran, G.W., Martin, K.J., and Thomas, G.E., 1977:
    Remark AS R19 and Algorithm AS 109: A remark on Algorithms: AS 63 the Incomplete Beta integral, AS 64 Inverse of the Incomplete Beta Function Ratio. Appl. Statist., Vol.26, No.1, pp 111-114.

function pinvbeta ( p, a, b, beta, acu, maxiter )

Purpose

Evaluates the inverse beta probability distribution function (e.g. the incomplete beta function).

For given arguments P (0<=P<=1), A (A>0.1), B (B>0.1), PINVBETA returns the value X such that P is the probability that a random variable distributed as beta(A,B) is less than or equal to X.

Arguments

P (INPUT) real(stnd)
On entry, input probability. P must be in the range (0,1) inclusive.
A (INPUT) real(stnd)
on entry, the (1st) parameter of the beta distribution. A must be greater than 0.1 .
B (INPUT) real(stnd)
On entry, the (2nd) parameter of the beta distribution. B must be greater than 0.1 .
BETA (INPUT, OPTIONAL) real(stnd)

On entry, the logarithm of the complete beta function beta(A,B).

If BETA is not given, the logarithm of the beta function is computed with the help of function LNGAMMA.

ACU (INPUT, OPTIONAL) real(stnd)

On entry, the desired accuracy of the result, when computing the incomplete beta function. The “integrating” process for evaluating the incomplete beta function is terminated when the relative contribution to the integral is not greater than the value of ACU. ACU is a small strictly positive integer.

See the description of the PROBBETA function for more details on this argument.

The default value for ACU is epsilon( ACU ).

MAXITER (INPUT, OPTIONAL) integer(i4b)

On entry, MAXITER controls the maximum number of iterations when evaluating the power series representation of the incomplete Beta integral.

See the description of the PROBBETA function for more details on this argument.

The default value is 2000.

Further Details

This function is not very accurate for small values of A and/or B (e.g. less than 0.5).

This function is adapted from:

  1. Majumder, K.L., and Bhattacharjee, G.P., 1973:
    Algorithm AS 64: Inverse of the Incomplete Beta Function Ratio. Appl. Statist., Vol.22, No.3, pp 411-414.
  2. Cran, G.W., Martin, K.J., and Thomas, G.E., 1977:
    Remark AS R19 and Algorithm AS 109: A remark on Algorithms: AS 63 the Incomplete Beta integral, AS 64 Inverse of the Incomplete Beta Function Ratio. Appl. Statist., Vol.26, No.1, pp 111-114.
  3. Berry, K.J., Mielke, P.W., and Cran, G.W., 1990:
    Algorithm AS R83: A remark on Algorithm AS 109: Inverse of the Incomplete Beta Function Ratio. Appl. Statist., Vol.39, No 2, pp. 309-310.
  4. Berry, K.J., Mielke, P.W., and Cran, G.W., 1991:
    Correction to Algorithm AS R83: A remark on Algorithm AS 109: Inverse of the Incomplete Beta Function Ratio. Appl. Statist. Vol. 40, No. 1, p.236

function probn ( x, upper )

Purpose

Evaluates the standard normal (Gaussian) distribution function from X to infinity if UPPER is true or from minus infinity to X if UPPER is false. In other words, if:

  • UPPER = true : PROBN = prob( U > X ) ,
  • UPPER = false : PROBN = prob( U < X ) ,

, for U = Laplace_Gauss(0;1).

Arguments

X (INPUT) real(stnd)
On entry, upper or lower limit of integration.
UPPER (INPUT) logical(lgl)

On entry, if:

  • UPPER = true : probability to the right of X is calculated.
  • UPPER = false : probability to the left of X is calculated.

Further Details

This function is adapted from:

  1. Hill, I.D., 1973:
    Algorithm AS66: The Normal Integral. Applied Statistics, vol.22, no.3, pp.424-427

function probn ( x, upper )

Purpose

Evaluates the standard normal (Gaussian) distribution function from X(i) to infinity if UPPER is true or from minus infinity to X(i) if UPPER is false, for i=1 to size(X). In other words, if:

  • UPPER = true : PROBN( i ) = prob( U > X(i) ) ,
  • UPPER = false : PROBN( i ) = prob( U < X(i) ) ,

, for U = Laplace_Gauss(0;1) and i=1 to size(X).

Arguments

X (INPUT) real(stnd), dimension(:)
On entry, upper or lower limits of integration.
UPPER (INPUT) logical(lgl)

On entry, if:

  • UPPER = true : probability to the right of X is calculated.
  • UPPER = false : probability to the left of X is calculated.

Further Details

The subroutine is parallelized if OPENMP is used.

This function is adapted from:

  1. Hill, I.D., 1973:
    Algorithm AS66: The Normal Integral. Applied Statistics, vol.22, no.3, pp.424-427

function probn ( x, upper )

Purpose

Evaluates the standard normal (Gaussian) distribution function from X(i,j) to infinity if UPPER is true or from minus infinity to X(i,j) if UPPER is false, for i=1 to size(X,1) and j=1 to size(X,2). In other words, if:

  • UPPER = true : PROBN( i, j ) = prob( U > X(i,j) ) ,
  • UPPER = false : PROBN( i, j ) = prob( U < X(i,j) ) ,

, for U = Laplace_Gauss(0;1) and i=1 to size(X,1) and j=1 to size(X,2).

Arguments

X (INPUT) real(stnd), dimension(:,:)
On entry, upper or lower limits of integration.
UPPER (INPUT) logical(lgl)

On entry, if:

  • UPPER = true : probability to the right of X is calculated.
  • UPPER = false : probability to the left of X is calculated.

Further Details

The subroutine is parallelized if OPENMP is used.

This function is adapted from:

  1. Hill, I.D., 1973:
    Algorithm AS66: The Normal Integral. Applied Statistics, vol.22, no.3, pp.424-427

function probn2 ( x, upper )

Purpose

Evaluates the standard normal (Gaussian) distribution function from X to infinity if UPPER is true or from minus infinity to X if UPPER is false. In other words, if:

  • UPPER = true : PROBN2 = prob( U > X ) ,
  • UPPER = false : PROBN2 = prob( U < X ) ,

, for U = Laplace_Gauss(0;1).

Arguments

X (INPUT) real(extd)
On entry, upper or lower limit of integration.
UPPER (INPUT) logical(lgl)

On entry, if:

  • UPPER = true : probability to the right of X is calculated.
  • UPPER = false : probability to the left of X is calculated.

Further Details

This function gives higher accuracy than PROBN.

This function is based upon algorithm 5666 for the error function, from:

  1. Hart, J.F. et al, 1968:
    Computer Approximations. 354 pp, New York, Wiley.

function probn2 ( x, upper )

Purpose

Evaluates the standard normal (Gaussian) distribution function from X(i) to infinity if UPPER is true or from minus infinity to X(i) if UPPER is false, for i=1 to size(X). In other words, if:

  • UPPER = true : PROBN2( i ) = prob( U > X(i) ) ,
  • UPPER = false : PROBN2( i ) = prob( U < X(i) ) ,

, for U = Laplace_Gauss(0;1) and i=1 to size(X).

Arguments

X (INPUT) real(extd), dimension(:)
On entry, upper or lower limits of integration.
UPPER (INPUT) logical(lgl)

On entry, if:

  • UPPER = true : probability to the right of X is calculated.
  • UPPER = false : probability to the left of X is calculated.

Further Details

The subroutine is parallelized if OPENMP is used.

This function gives higher accuracy than PROBN.

This function is based upon algorithm 5666 for the error function, from:

  1. Hart, J.F. et al, 1968:
    Computer Approximations. 354 pp, New York, Wiley.

function probn2 ( x, upper )

Purpose

Evaluates the standard normal (Gaussian) distribution function from X(i,j) to infinity if UPPER is true or from minus infinity to X(i,j) if UPPER is false, for i=1 to size(X,1) and j=1 to size(X,2). In other words, if:

  • UPPER = true : PROBN2( i, j ) = prob( U > X(i,j) ) ,
  • UPPER = false : PROBN2( i, j ) = prob( U < X(i,j) ) ,

, for U = Laplace_Gauss(0;1) and i=1 to size(X,1) and j=1 to size(X,2).

Arguments

X (INPUT) real(extd), dimension(:,:)
On entry, upper or lower limits of integration.
UPPER (INPUT) logical(lgl)

On entry, if:

  • UPPER = true : probability to the right of X is calculated.
  • UPPER = false : probability to the left of X is calculated.

Further Details

The subroutine is parallelized if OPENMP is used.

This function gives higher accuracy than PROBN.

This function is based upon algorithm 5666 for the error function, from:

  1. Hart, J.F. et al, 1968:
    Computer Approximations. 354 pp, New York, Wiley.

function pinvn ( p )

Purpose

Evaluates the inverse of the standard normal (Gaussian) distribution function:

X0 = PINVN( P )

, if P = probability( U < X0 ) for U = Laplace_Gauss(0;1).

PINVN returns the normal deviate X0 corresponding to a given lower tail area of P.

Arguments

P (INPUT) real(stnd)
On entry, the probability. P must verify 0. < P < 1.

Further Details

The inverse Gaussian CDF is approximated to high precision using rational approximations (polynomials with degree 2 and 3) by the subroutine PPND7 described in the reference (1).

This function is accurate to about seven decimal figures for min(P,1-P) > 10**(-316).

If P is very close to unity, a serious loss of significance may be incurred in forming 1 - P = c in the code of the function. In this circumstance the user should, if possible, evaluate c directly or in extended precision and evaluate X0(P) = PINVN( P ) as -X0(c) = PINVN( c ).

The hash sums below are the sums of the mantissas of the coefficients. They are included for use in checking transcription.

This function is adapted from the routine PPND7 described in:

  1. Wichura, M.J., 1988:
    Algorithm AS 241: The percentage points of the normal distribution. Appl. Statist., Vol. 37, No. 3, pp. 477-484.

function pinvn ( p )

Purpose

Evaluates the inverse of the standard normal (Gaussian) distribution function:

X0(i) = PINVN( P(i) )

, if P(i) = probability( U < X0(i) ) for U = Laplace_Gauss(0;1) and i=1 to size(P).

PINVN returns the normal deviate X0(i) corresponding to a given lower tail area of P(i), for i=1 to size(P).

Arguments

P (INPUT) real(stnd), dimension(:)
On entry, the probabilities. P(i) must verify 0. < P(i) < 1, for i=1 to size(P) .

Further Details

The inverse Gaussian CDF is approximated to high precision using rational approximations (polynomials with degree 2 and 3) by the subroutine PPND7 described in the reference (1).

This function is accurate to about seven decimal figures for min(P,1-P) > 10**(-316).

If P is very close to unity, a serious loss of significance may be incurred in forming 1 - P = c in the code of the function. In this circumstance the user should, if possible, evaluate c directly or in extended precision and evaluate X0(P) = PINVN( P ) as -X0(c) = PINVN( c ).

The subroutine is parallelized if OPENMP is used.

The hash sums below are the sums of the mantissas of the coefficients. They are included for use in checking transcription.

This function is adapted from the routine PPND7 described in:

  1. Wichura, M.J., 1988:
    Algorithm AS 241: The percentage points of the normal distribution. Appl. Statist., Vol. 37, No. 3, pp. 477-484.

function pinvn ( p )

Purpose

Evaluates the inverse of the standard normal (Gaussian) distribution function:

X0(i,j) = PINVN( P(i,j) )

, if P(i,j) = probability( U < X0(i,j) ) for U = Laplace_Gauss(0;1), i=1 to size(P,1) and j=1 to size(P,2).

PINVN returns the normal deviate X0(i,j) corresponding to a given lower tail area of P(i,j) for i=1 to size(P,1) and j=1 to size(P,2).

Arguments

P (INPUT) real(stnd), dimension(:,:)
On entry, the probabilities. P(i,j) must verify 0. < P(i,j) < 1, for i=1 to size(P,1) and j=1 to size(P,2) .

Further Details

The inverse Gaussian CDF is approximated to high precision using rational approximations (polynomials with degree 2 and 3) by the subroutine PPND7 described in the reference (1).

This function is accurate to about seven decimal figures for min(P,1-P) > 10**(-316).

If P is very close to unity, a serious loss of significance may be incurred in forming 1 - P = c in the code of the function. In this circumstance the user should, if possible, evaluate c directly or in extended precision and evaluate X0(P) = PINVN( P ) as -X0(c) = PINVN( c ).

The subroutine is parallelized if OPENMP is used.

The hash sums below are the sums of the mantissas of the coefficients. They are included for use in checking transcription.

This function is adapted from the routine PPND7 described in:

  1. Wichura, M.J., 1988:
    Algorithm AS 241: The percentage points of the normal distribution. Appl. Statist., Vol. 37, No. 3, pp. 477-484.

function pinvn2 ( p )

Purpose

Evaluates the inverse of the standard normal (Gaussian) distribution function:

X0 = PINVN2( P )

, if P = probability( U < X0 ) for U = Laplace_Gauss(0;1).

PINVN2 returns the normal deviate X0 corresponding to a given lower tail area of P.

Arguments

P (INPUT) real(extd)
On entry, the probability. P must verify 0. < P < 1.

Further Details

The inverse Gaussian CDF is approximated to high precision using rational approximations (polynomials with degree 7) by the subroutine PPND16 described in the reference (1).

This function is accurate to about 16 decimal figures for min(P,1-P) > 10**(-316) and gives higher accuracy than PINVN function, but it is slower.

On a machine, that uses only 32 bits to represent real variables, PINVN2 should be implemented in double precision (e.g. with a correct choice of the kind EXTD in the module Select_Parameters ).

If P is very close to unity, a serious loss of significance may be incurred in forming 1 - P = c in the code of the function. In this circumstance the user should, if possible, evaluate c directly or in extended precision and evaluate X0(P) = PINVN2( P ) as -X0(c) = PINVN2( c ) .

The hash sums below are the sums of the mantissas of the coefficients. They are included for use in checking transcription.

This function is adapted from the routine PPND16 described in:

  1. Wichura, M.J., 1988:
    Algorithm AS 241: The percentage points of the normal distribution. Appl. Statist., Vol. 37, No. 3, pp. 477-484.

function pinvn2 ( p )

Purpose

Evaluates the inverse of the standard normal (Gaussian) distribution function:

X0(i) = PINVN2( P(i) )

, if P(i) = probability( U < X0(i) ) for U = Laplace_Gauss(0;1) and i=1 to size(P).

PINVN returns the normal deviate X0(i) corresponding to a given lower tail area of P(i), for i=1 to size(P).

Arguments

P (INPUT) real(extd), dimension(:)
On entry, the probabilities. P(i) must verify 0. < P(i) < 1, for i=1 to size(P) .

Further Details

The inverse Gaussian CDF is approximated to high precision using rational approximations (polynomials with degree 7) by the subroutine PPND16 described in the reference (1).

This function is accurate to about 16 decimal figures for min(P,1-P) > 10**(-316) and gives higher accuracy than PINVN function, but it is slower.

On a machine, that uses only 32 bits to represent real variables, PINVN2 should be implemented in double precision (e.g. with a correct choice of the kind EXTD in the module Select_Parameters ).

If P is very close to unity, a serious loss of significance may be incurred in forming 1 - P = c in the code of the function. In this circumstance the user should, if possible, evaluate c directly or in extended precision and evaluate X0(P) = PINVN2( P ) as -X0(c) = PINVN2( c ).

The subroutine is parallelized if OPENMP is used.

The hash sums below are the sums of the mantissas of the coefficients. They are included for use in checking transcription.

This function is adapted from the routine PPND16 described in:

  1. Wichura, M.J., 1988:
    Algorithm AS 241: The percentage points of the normal distribution. Appl. Statist., Vol. 37, No. 3, pp. 477-484.

function pinvn2 ( p )

Purpose

Evaluates the inverse of the standard normal (Gaussian) distribution function:

X0(i,j) = PINVN2( P(i,j) )

, if P(i,j) = probability( U < X0(i,j) ) for U = Laplace_Gauss(0;1), i=1 to size(P,1) and j=1 to size(P,2).

PINVN2 returns the normal deviate X0(i,j) corresponding to a given lower tail area of P(i,j) for i=1 to size(P,1) and j=1 to size(P,2).

Arguments

P (INPUT) real(extd), dimension(:,:)
On entry, the probabilities. P(i,j) must verify 0. < P(i,j) < 1, for i=1 to size(P,1) and j=1 to size(P,2) .

Further Details

The inverse Gaussian CDF is approximated to high precision using rational approximations (polynomials with degree 7) by the subroutine PPND16 described in the reference (1).

This function is accurate to about 16 decimal figures for min(P,1-P) > 10**(-316) and gives higher accuracy than PINVN function, but it is slower.

On a machine, that uses only 32 bits to represent real variables, PINVN2 should be implemented in double precision (e.g. with a correct choice of the kind EXTD in the module Select_Parameters ).

If P is very close to unity, a serious loss of significance may be incurred in forming 1 - P = c in the code of the function. In this circumstance the user should, if possible, evaluate c directly or in extended precision and evaluate X0(P) = PINVN2( P ) as -X0(c) = PINVN2( c ).

The subroutine is parallelized if OPENMP is used.

The hash sums below are the sums of the mantissas of the coefficients. They are included for use in checking transcription.

This function is adapted from the routine PPND16 described in:

  1. Wichura, M.J., 1988:
    Algorithm AS 241: The percentage points of the normal distribution. Appl. Statist., Vol. 37, No. 3, pp. 477-484.

function probt ( t, ndf, upper, ndf_max )

Purpose

Evaluates the Student’s t-distribution function from T to infinity if UPPER is true or from minus infinity to T if UPPER is false. In otherwords, if:

  • UPPER = true : PROBT = prob( U > T ) ,
  • UPPER = false : PROBT = prob( U < T ) ,

, for U = Student(NDF).

Arguments

T (INPUT) real(stnd)
On entry, upper or lower limit of integration of the t-density.
NDF (INPUT) integer(i4b)
On entry, degrees of freedom of the t-distribution. NDF must be greater or equal to 1.
UPPER (INPUT) logical(lgl)

On entry, if:

  • UPPER = true : probability to the right of T is calculated.
  • UPPER = false : probability to the left of T is calculated.
NDF_MAX (INPUT, OPTIONAL) integer(i4b)

On entry, if:

  • NDF is lower or equal to NDF_MAX, the t_density is integrated.
  • NDF is greater than NDF_MAX, an asymptotic series is used.

The default is 20.

Further Details

This function is adapted from:

  1. Hill, G.W., 1970:
    Student’s t-distribution (Algorithm 395). Comm. A.C.M., vol.13, 617-619
  2. Cooper, B.E., 1968:
    The integral of student’s t distribution, (Algorithm AS3). Applied Statistics, vol.17, no.2, 189

function probt ( t, ndf, upper, ndf_max )

Purpose

Evaluates the Student’s t-distribution function from T(i) to infinity if UPPER is true or from minus infinity to T(i) if UPPER is false, for i=1 to size(T). In other words, if:

  • UPPER = true : PROBT( i ) = prob( U > T(i) ) ,
  • UPPER = false : PROBT( i ) = prob( U < T(i) ) ,

, for U = STUDENT(NDF) and i=1 to size(T).

Arguments

T (INPUT) real(stnd), dimension(:)
On entry, upper or lower limits of integration of the t-density.
NDF (INPUT) integer(i4b)
On entry, degrees of freedom of the t-distribution. NDF must be greater or equal to 1.
UPPER (INPUT) logical(lgl)

On entry, if:

  • UPPER = true : probabilities to the right of T is calculated.
  • UPPER = false : probabilities to the left of T is calculated.
NDF_MAX (INPUT, OPTIONAL) integer(i4b)

On entry, if:

  • NDF is lower or equal to NDF_MAX, the t_density is integrated.
  • NDF is greater than NDF_MAX, an asymptotic series is used.

The default is 20.

Further Details

This function is parallelized if OPENMP is used.

This function is adapted from:

  1. Hill, G.W., 1970:
    Student’s t-distribution (Algorithm 395). Comm. A.C.M., vol.13, 617-619
  2. Cooper, B.E., 1968:
    The integral of student’s t distribution, (Algorithm AS3). Applied Statistics, vol.17, no.2, 189

function probt ( t, ndf, upper, ndf_max )

Purpose

Evaluates the Student’s t-distribution function from T(i) to infinity if UPPER is true or from minus infinity to T(i) if UPPER is false, for i=1 to size(T). In other words, if:

  • UPPER = true : PROBT( i ) = prob( U > T(i) ) ,
  • UPPER = false : PROBT( i ) = prob( U < T(i) ) ,

, for U = STUDENT(NDF(i)) and i=1 to size(T).

Arguments

T (INPUT) real(stnd), dimension(:)
On entry, upper or lower limits of integration of the t-density.
NDF (INPUT) integer(i4b), dimension(:)

On entry, degrees of freedom of the t-distribution. Any value in the array NDF must be greater or equal to 1.

The size of NDF must be size(NDF) = size(T) .

UPPER (INPUT) logical(lgl)

On entry, if:

  • UPPER = true : probabilities to the right of T is calculated.
  • UPPER = false : probabilities to the left of T is calculated.
NDF_MAX (INPUT, OPTIONAL) integer(i4b)

On entry, if:

  • NDF is lower or equal to NDF_MAX, the t_density is integrated.
  • NDF is greater than NDF_MAX, an asymptotic series is used.

The default is 20.

Further Details

This function is parallelized if OPENMP is used.

This function is adapted from:

  1. Hill, G.W., 1970:
    Student’s t-distribution (Algorithm 395). Comm. A.C.M., vol.13, 617-619
  2. Cooper, B.E., 1968:
    The integral of student’s t distribution, (Algorithm AS3). Applied Statistics, vol.17, no.2, 189

function probt ( t, ndf, upper, ndf_max )

Purpose

Evaluates the Student’s t-distribution function from T(i,j) to infinity if UPPER is true or from minus infinity to T(i,j) if UPPER is false, for i=1 to size(T,1) and j=1 to size(T,2). In other words, if:

  • UPPER = true : PROBT( i, j ) = prob( U > T(i,j) ) ,
  • UPPER = false : PROBT( i, j ) = prob( U < T(i,j) ) ,

, for U = Student(NDF), i=1 to size(T,1) and j=1 to size(T,2).

Arguments

T (INPUT) real(stnd), dimension(:,:)
On entry, upper or lower limits of integration of the t-density.
NDF (INPUT) integer(i4b)
On entry, degrees of freedom of the t-distribution. NDF must be greater or equal to 1.
UPPER (INPUT) logical(lgl)

On entry, if:

  • UPPER = true : probabilities to the right of T is calculated.
  • UPPER = false : probabilities to the left of T is calculated.
NDF_MAX (INPUT, OPTIONAL) integer(i4b)

On entry, if:

  • NDF is lower or equal to NDF_MAX, the t_density is integrated.
  • NDF is greater than NDF_MAX, an asymptotic series is used.

The default is 20.

Further Details

This function is parallelized if OPENMP is used.

This function is adapted from:

  1. Hill, G.W., 1970:
    Student’s t-distribution (Algorithm 395). Comm. A.C.M., vol.13, 617-619
  2. Cooper, B.E., 1968:
    The integral of student’s t distribution, (Algorithm AS3). Applied Statistics, vol.17, no.2, 189

function probt ( t, ndf, upper, ndf_max )

Purpose

Evaluates the Student’s t-distribution function from T(i,j) to infinity if UPPER is true or from minus infinity to T(i,j) if UPPER is false, for i=1 to size(T,1) and j=1 to size(T,2). In other words, if:

  • UPPER = true : PROBT( i, j ) = prob( U > T(i,j) ) ,
  • UPPER = false : PROBT( i, j ) = prob( U < T(i,j) ) ,

, for U = Student(NDF(i,j)), i=1 to size(T,1) and j=1 to size(T,2).

Arguments

T (INPUT) real(stnd), dimension(:,:)
On entry, upper or lower limits of integration of the t-density.
NDF (INPUT) integer(i4b), dimension(:,:)

On entry, degrees of freedom of the t-distribution. Any value in the array NDF must be greater or equal to 1.

The shape of NDF must verify:

  • size(NDF,1) = size(T,1)
  • size(NDF,2) = size(T,2).
UPPER (INPUT) logical(lgl)

On entry, if:

  • UPPER = true : probabilities to the right of T is calculated.
  • UPPER = false : probabilities to the left of T is calculated.
NDF_MAX (INPUT, OPTIONAL) integer(i4b)

On entry, if:

  • NDF is lower or equal to NDF_MAX, the t_density is integrated.
  • NDF is greater than NDF_MAX, an asymptotic series is used.

The default is 20.

Further Details

This function is parallelized if OPENMP is used.

This function is adapted from:

  1. Hill, G.W., 1970:
    Student’s t-distribution (Algorithm 395). Comm. A.C.M., vol.13, 617-619
  2. Cooper, B.E., 1968:
    The integral of student’s t distribution, (Algorithm AS3). Applied Statistics, vol.17, no.2, 189

function probstudent ( t, df )

Purpose

Evaluates the two-tailed probability of Student’s t. PROBSTUDENT computes the probability that a random variable following Student’s t distribution will exceed abs(T) in absolute value.

Arguments

T (INPUT) real(stnd)
On entry, input constant. PROBSTUDENT computes the probability that abs(T) will be exceeded in absolute value.
DF (INPUT) real(stnd)
On entry, degrees of freedom of the t-distribution. DF must be greater or equal to 1. DF is not necessarily an integer.

Further Details

This function is not very accurate for very small degrees of freedom (e.g. number of degrees of freedom less than 5).

This function is adapted from:

  1. Hill, G.W., 1970:
    Student’s t-distribution (Algorithm 395). Comm. A.C.M., vol.13, 617-619

function probstudent ( t, df )

Purpose

Evaluates the two-tailed probabilities of Student’s t. PROBSTUDENT computes probabilities that a random variable following Student’s t distribution will exceed abs(T(i)) in absolute value, for i=1 to size(T).

Arguments

T (INPUT) real(stnd), dimension(:)
On entry, input constants. PROBSTUDENT computes probabilities that abs(T(i)) will be exceeded in absolute value, for i=1 to size(T).
DF (INPUT) real(stnd)
On entry, degrees of freedom of the t-distribution. DF must be greater or equal to 1. DF is not necessarily an integer.

Further Details

This function is not very accurate for very small degrees of freedom (e.g. number of degrees of freedom less than 5).

This function is parallelized if OPENMP is used.

This function is adapted from:

  1. Hill, G.W., 1970:
    Student’s t-distribution (Algorithm 395). Comm. A.C.M., vol.13, 617-619

function probstudent ( t, df )

Purpose

Evaluates the two-tailed probabilities of Student’s t. PROBSTUDENT computes probabilities that a random variable following Student’s t distribution will exceed abs(T(i)) in absolute value, for i=1 to size(T).

Arguments

T (INPUT) real(stnd), dimension(:)
On entry, input constants. PROBSTUDENT computes probabilities that abs(T(i)) will be exceeded in absolute value, for i=1 to size(T).
DF (INPUT) real(stnd), dimension(:)

On entry, degrees of freedom of the t-distribution. Any value in the array DF must be greater or equal to 1, but is not necessarily an integer.

The size of DF must be size(DF) = size(T) .

Further Details

This function is not very accurate for very small degrees of freedom (e.g. number of degrees of freedom less than 5).

This function is parallelized if OPENMP is used.

This function is adapted from:

  1. Hill, G.W., 1970:
    Student’s t-distribution (Algorithm 395). Comm. A.C.M., vol.13, 617-619

function probstudent ( t, df )

Purpose

Evaluates two-tailed probabilities of Student’s t. PROBSTUDENT computes probabilities that a random variable following Student’s t distribution will exceed abs(T(i,j)) in absolute value, for i=1 to size(T,1) and j=1 to size(T,2).

Arguments

T (INPUT) real(stnd), dimension(:,:)
On entry, input constants. PROBSTUDENT computes probabilities that abs(T(i,j)) will be exceeded in absolute value, for i=1 to size(T,1) and j=1 to size(T,2).
DF (INPUT) real(stnd)
On entry, degrees of freedom of the t-distribution. DF must be greater or equal to 1. DF is not necessarily an integer.

Further Details

This function is not very accurate for very small degrees of freedom (e.g. number of degrees of freedom less than 5).

This function is parallelized if OPENMP is used.

This function is adapted from:

  1. Hill, G.W., 1970:
    Student’s t-distribution (Algorithm 395). Comm. A.C.M., vol.13, 617-619

function probstudent ( t, df )

Purpose

Evaluates two-tailed probabilities of Student’s t. PROBSTUDENT computes probabilities that a random variable following Student’s t distribution will exceed abs(T(i,j)) in absolute value, for i=1 to size(T,1) and j=1 to size(T,2).

Arguments

T (INPUT) real(stnd), dimension(:,:)
On entry, input constants. PROBSTUDENT computes probabilities that abs(T(i,j)) will be exceeded in absolute value, for i=1 to size(T,1) and j=1 to size(T,2).
DF (INPUT) real(stnd), dimension(:,:)

On entry, degrees of freedom of the t-distribution. Any value in the array DF must be greater or equal to 1, but is not necessarily an integer.

The shape of DF must verify:

  • size(DF,1) = size(T,1)
  • size(DF,2) = size(T,2).

Further Details

This function is not very accurate for very small degrees of freedom (e.g. number of degrees of freedom less than 5).

This function is parallelized if OPENMP is used.

This function is adapted from:

  1. Hill, G.W., 1970:
    Student’s t-distribution (Algorithm 395). Comm. A.C.M., vol.13, 617-619

function pinvt ( p, ndf )

Purpose

Evaluates the inverse of the Student’s t distribution function:

T0 = PINVT( P, NDF )

, if P = probability( U < T0 ) for U = Student(NDF).

PINVT returns the quantile T0 of Student’s t-distribution with NDF degrees of freedom corresponding to a given lower tail area of P.

Arguments

P (INPUT) real(stnd)
On entry, the probability. P must verify 0. < P < 1.
NDF (INPUT) integer(i4b)
On entry, degrees of freedom of the t-distribution. NDF must be greater or equal to 1.

Further Details

This function is adapted from:

  1. Hill, G.W., 1970:
    Student’s t-quantiles (Algorithm 396). Comm. A.C.M., vol.13, no10, 620-621

function pinvt ( p, ndf )

Purpose

Evaluates the inverse of the Student’s t distribution function:

T0(i) = PINVT( P(i), NDF )

, if P(i) = probability( U < T0(i) ) for U = Student(NDF) and i=1 to size(P).

PINVT returns the quantiles T0(i) of Student’s t-distribution with NDF degrees of freedom corresponding to a given lower tail area of P(i), for i=1 to size(P).

Arguments

P (INPUT) real(stnd), dimension(:)
On entry, the probabilities. P(i) must verify 0. < P(i) < 1, for i=1 to size(P).
NDF (INPUT) integer(i4b)
On entry, degrees of freedom of the t-distribution. NDF must be greater or equal to 1.

Further Details

This function is parallelized if OPENMP is used.

This function is adapted from:

  1. Hill, G.W., 1970:
    Student’s t-quantiles (Algorithm 396). Comm. A.C.M., vol.13, no10, 620-621

function pinvt ( p, ndf )

Purpose

Evaluates the inverse of the Student’s t distribution function:

T0(i) = PINVT( P(i), NDF(i) )

, if P(i) = probability( U < T0(i) ) for U = Student(NDF(i)) and i=1 to size(P).

PINVT returns the quantiles T0(i) of Student’s t-distribution with NDF(i) degrees of freedom corresponding to a given lower tail area of P(i), for i=1 to size(P).

Arguments

P (INPUT) real(stnd), dimension(:)
On entry, the probabilities. P(i) must verify 0. < P(i) < 1, for i=1 to size(P).
NDF (INPUT) integer(i4b), dimension(:)

On entry, degrees of freedom of the t-distribution. Any value in the array NDF must be greater or equal to 1.

The size of NDF must be size(NDF) = size(P).

Further Details

This function is parallelized if OPENMP is used.

This function is adapted from:

  1. Hill, G.W., 1970:
    Student’s t-quantiles (Algorithm 396). Comm. A.C.M., vol.13, no10, 620-621

function pinvt ( p, ndf )

Purpose

Evaluates the inverse of the Student’s t distribution function:

T0(i,j) = PINVT( P(i,j), NDF )

, if P(i,j) = probability( U < T0(i,j) ) for U = Student(NDF), i=1 to size(P,1) and j=1 to size(P,2).

PINVT returns the quantiles T0(i,j) of Student’s t-distribution with NDF degrees of freedom corresponding to a given lower tail area of P(i,j) for i=1 to size(P,1) and j=1 to size(P,2).

Arguments

P (INPUT) real(stnd), dimension(:,:)
On entry, the probabilities. P(i,j) must verify 0. < P(i,j) < 1, for i=1 to size(P,1) and j=1 to size(P,2).
NDF (INPUT) integer(i4b)
On entry, degrees of freedom of the t-distribution. NDF must be greater or equal to 1.

Further Details

This function is parallelized if OPENMP is used.

This function is adapted from:

  1. Hill, G.W., 1970:
    Student’s t-quantiles (Algorithm 396). Comm. A.C.M., vol.13, no10, 620-621

function pinvt ( p, ndf )

Purpose

Evaluates the inverse of the Student’s t distribution function:

T0(i,j) = PINVT( P(i,j), NDF(i,j) )

, if P(i,j) = probability( U < T0(i,j) ) for U = Student(NDF(i,j)), i=1 to size(P,1) and j=1 to size(P,2).

PINVT returns the quantiles T0(i,j) of Student’s t-distribution with NDF(i,j) degrees of freedom corresponding to a given lower tail area of P(i,j) for i=1 to size(P,1) and j=1 to size(P,2) .

Arguments

P (INPUT) real(stnd), dimension(:,:)
On entry, the probabilities. P(i,j) must verify 0. < P(i,j) < 1, for i=1 to size(P,1) and j=1 to size(P,2).
NDF (INPUT) integer(i4b), dimension(:,:)

On entry, degrees of freedom of the t-distribution. Any value in the array NDF must be greater or equal to 1.

The shape of NDF must verify:

  • size(NDF,1) = size(P,1)
  • size(NDF,2) = size(P,2).

Further Details

This function is parallelized if OPENMP is used.

This function is adapted from:

  1. Hill, G.W., 1970:
    Student’s t-quantiles (Algorithm 396). Comm. A.C.M., vol.13, no10, 620-621

function pinvstudent ( p, df )

Purpose

Evaluates the inverse of a modification of Student’s t probability distribution function.

PINVSTUDENT calculates the two-tail quantiles of Student’s t-distribution, that is a value x such that the probability of the absolute value of t being greater than X is P.

Arguments

P (INPUT) real(stnd)
On entry, the probability. P is the sum of the areas (equal) in both tails of the t-distribution. P must verify 0. < P < 1.
DF (INPUT) real(stnd)
On entry, degrees of freedom of the t-distribution. DF must be greater or equal to 1. . DF is not necessarily an integer.

Further Details

Note that PINVSTUDENT does not provide the actual Student’s t inverse. For q equal to the probability that a Student’s t random variable is less than x, that inverse can be obtained by the following rules:

  • for q in the range (0.0,0.5), call PINVSTUDENT with P = 2 * q and negate the result x.
  • for q in the range (0.5,1.0), call PINVSTUDENT with P = 2 * (1-q).

This function is adapted from:

  1. Hill, G.W., 1970:
    Student’s t-quantiles (Algorithm 396). Comm. A.C.M., vol.13, no10, 620-621

function pinvstudent ( p, df )

Purpose

Evaluates the inverse of a modification of Student’s t probability distribution function.

PINVSTUDENT calculates the two-tail quantiles of Student’s t-distribution, that is a value x(i) such that the probability of the absolute value of t being greater than x(i) is P(i), for i=1 to size(P).

Arguments

P (INPUT) real(stnd), dimension(:)
On entry, the probabilities. P(i) is the sum of the areas (equal) in both tails of the t-distribution, for i=1 to size(P). P(i) must verify 0. < P(i) < 1, for i=1 to size(P).
DF (INPUT) real(stnd)
On entry, degrees of freedom of the t-distribution. DF must be greater or equal to 1. . DF is not necessarily an integer.

Further Details

Note that PINVSTUDENT does not provide the actual Student’s t inverse. For q(:) equal to the probabilities that a Student’s t random variable is less than x(:), that inverse can be obtained by the following rules:

  • for q(:) in the range (0.0,0.5), call PINVSTUDENT with P(:) = 2 * q(:) and negate the result x(:).
  • for q(:) in the range (0.5,1.0), call PINVSTUDENT with P(:) = 2 * (1-q(:)).

This function is parallelized if OPENMP is used.

This function is adapted from:

  1. Hill, G.W., 1970:
    Student’s t-quantiles (Algorithm 396). Comm. A.C.M., vol.13, no10, 620-621

function pinvstudent ( p, df )

Purpose

Evaluates the inverse of a modification of Student’s t probability distribution function.

PINVSTUDENT calculates the two-tail quantiles of Student’s t-distribution, that is a value x(i) such that the probability of the absolute value of t being greater than x(i) is P(i), for i=1 to size(P).

Arguments

P (INPUT) real(stnd), dimension(:)
On entry, the probabilities. P(i) is the sum of the areas (equal) in both tails of the t-distribution, for i=1 to size(P). P(i) must verify 0. < P(i) < 1, for i=1 to size(P).
DF (INPUT) real(stnd), dimension(:)

On entry, degrees of freedom of the t-distribution. Any value in the array DF must be greater or equal to 1, but is not necessarily an integer.

The size of DF must be size(DF) = size(P) .

Further Details

Note that PINVSTUDENT does not provide the actual Student’s t inverse. For q(:) equal to the probabilities that a Student’s t random variable is less than x(:), that inverse can be obtained by the following rules:

  • for q(:) in the range (0.0,0.5), call PINVSTUDENT with P(:) = 2 * q(:) and negate the result x(:).
  • for q(:) in the range (0.5,1.0), call PINVSTUDENT with P(:) = 2 * (1-q(:)).

This function is parallelized if OPENMP is used.

This function is adapted from:

  1. Hill, G.W., 1970:
    Student’s t-quantiles (Algorithm 396). Comm. A.C.M., vol.13, no10, 620-621

function pinvstudent ( p, df )

Purpose

Evaluates the inverse of a modification of Student’s t probability distribution function.

PINVSTUDENT calculates the two-tail quantiles of Student’s t-distribution, that is a value x(i,j) such that the probability of the absolute value of t being greater than x(i,j) is P(i,j), for i=1 to size(P,1) and j=1 to size(P,2).

Arguments

P (INPUT) real(stnd), dimension(:,:)
On entry, the probabilities. P(i,j) is the sum of the areas (equal) in both tails of the t-distribution, for i=1 to size(P) and j=1 to size(P,2) . P(i,j) must verify 0. < P(i,j) < 1, for i=1 to size(P,1) and j=1 to size(P,2) .
DF (INPUT) real(stnd)
On entry, degrees of freedom of the t-distribution. DF must be greater or equal to 1. . DF is not necessarily an integer.

Further Details

Note that PINVSTUDENT does not provide the actual Student’s t inverse. For q(:,:) equal to the probabilities that a Student’s t random variable is less than x(:,:), that inverse can be obtained by the following rules:

  • for q(:,:) in the range (0.0,0.5), call PINVSTUDENT with P(:,:) = 2 * q(:,:) and negate the result x(:,:).
  • for q(:,:) in the range (0.5,1.0), call PINVSTUDENT with P(:,:) = 2 * (1-q(:,:)).

This function is parallelized if OPENMP is used.

This function is adapted from:

  1. Hill, G.W., 1970:
    Student’s t-quantiles (Algorithm 396). Comm. A.C.M., vol.13, no10, 620-621

function pinvstudent ( p, df )

Purpose

Evaluates the inverse of a modification of Student’s t probability distribution function.

PINVSTUDENT calculates the two-tail quantiles of Student’s t-distribution, that is a value x(i,j) such that the probability of the absolute value of t being greater than x(i,j) is P(i,j), for i=1 to size(P,1) and j=1 to size(P,2).

Arguments

P (INPUT) real(stnd), dimension(:,:)
On entry, the probabilities. P(i,j) is the sum of the areas (equal) in both tails of the t-distribution, for i=1 to size(P) and j=1 to size(P,2) . P(i,j) must verify 0. < P(i,j) < 1, for i=1 to size(P,1) and j=1 to size(P,2) .
DF (INPUT) real(stnd), dimension(:,:)

On entry, degrees of freedom of the t-distribution. Any value in the array DF must be greater or equal to 1, but is not necessarily an integer.

The shape of DF must verify:

  • size(DF,1) = size(P,1)
  • size(DF,2) = size(P,2) .

Further Details

Note that PINVSTUDENT does not provide the actual Student’s t inverse. For q(:,:) equal to the probabilities that a Student’s t random variable is less than x(:,:), that inverse can be obtained by the following rules:

  • for q(:,:) in the range (0.0,0.5), call PINVSTUDENT with P(:,:) = 2 * q(:,:) and negate the result x(:,:).
  • for q(:,:) in the range (0.5,1.0), call PINVSTUDENT with P(:,:) = 2 * (1-q(:,:)).

This function is parallelized if OPENMP is used.

This function is adapted from:

  1. Hill, G.W., 1970:
    Student’s t-quantiles (Algorithm 396). Comm. A.C.M., vol.13, no10, 620-621

function probq ( x2, ndf, upper, ndf_max )

Purpose

Evaluates the chi-squared distribution function from X2 to infinity if UPPER is true or from zero to X2 if UPPER is false. In other words, if:

  • UPPER = true : PROBQ = prob( U > X2 ) ,
  • UPPER = false : PROBQ = prob( U < X2 ) ,

, for U = Chi-squared(NDF).

Arguments

X2 (INPUT) real(stnd)
On entry, upper or lower limit of integration. X2 must be greater or equal to zero.
NDF (INPUT) integer(i4b)
On entry, degrees of freedom of the chi-squared distribution. NDF must be greater or equal to 1.
UPPER (INPUT) logical(lgl)

On entry, if:

  • UPPER = true : probability to the right of X2 is calculated.
  • UPPER = false : probability to the left of X2 is calculated.
NDF_MAX (INPUT, OPTIONAL) integer(i4b)

On entry, if:

  • NDF is lower or equal to NDF_MAX, the chi-squared density is integrated.
  • NDF is greater than NDF_MAX, a gaussian approximation is used.

The default is 40.

Further Details

If NDF<=NDF_MAX, the chi-squared distribution function is integrating by using formulae 26.4.4 and 26.4.5 in reference (1), otherwise a normal approximation based on the Wilson-Hilferty transformation is used (see the reference (1), Formula 26.4.14).

This function works only for integer degrees of freedom. It may be faster than PROBQ2 or PROBQ3 functions for the default value of NDF_MAX, but it is less accurate.

This function is adapted from:

  1. Abramowitz, M., and Stegun, I.A., 1970:
    Handbook of Mathematical Functions, (formulae 26.4.4, 26.4.5 and 26.4.14). New York, Dover Publications
  2. Wilson, E.B., and Hilferty, M.M., 1931:
    The distribution of Chi-square. Proceed. Nation. Academ. Scien., Vol. 17, 684-688

function probq ( x2, ndf, upper, ndf_max )

Purpose

Evaluates the chi-squared distribution function from X2(i) to infinity if UPPER is true or from zero to X2(i) if UPPER is false, for i=1 to size(X2). In other words, if:

  • UPPER = true : PROBQ( i ) = prob( U > X2(i) ) ,
  • UPPER = false : PROBQ( i ) = prob( U < X2(i) ) ,

, for U = Chi-squared(NDF) and i=1 to size(X2).

Arguments

X2 (INPUT) real(stnd), dimension(:)
On entry, upper or lower limits of integration. X2(i) must be greater or equal to zero for i=1 to size(X2).
NDF (INPUT) integer(i4b)
On entry, degrees of freedom of the chi-squared distribution. NDF must be greater or equal to 1.
UPPER (INPUT) logical(lgl)

On entry, if:

  • UPPER = true : probability to the right of X2 is calculated.
  • UPPER = false : probability to the left of X2 is calculated.
NDF_MAX (INPUT, OPTIONAL) integer(i4b)

On entry, if:

  • NDF is lower or equal to NDF_MAX, the chi-squared density is integrated.
  • NDF is greater than NDF_MAX, a gaussian approximation is used.

The default is 40.

Further Details

If NDF<=NDF_MAX, the chi-squared distribution function is integrating by using formulae 26.4.4 and 26.4.5 in reference (1), otherwise a normal approximation based on the Wilson-Hilferty transformation is used (see the reference (1), Formula 26.4.14).

This function works only for integer degrees of freedom. It may be faster than PROBQ2 or PROBQ3 functions for the default value of NDF_MAX, but it is less accurate.

The function is parallelized if OPENMP is used.

This function is adapted from:

  1. Abramowitz, M., and Stegun, I.A., 1970:
    Handbook of Mathematical Functions, (formulae 26.4.4, 26.4.5 and 26.4.14). New York, Dover Publications
  2. Wilson, E.B., and Hilferty, M.M., 1931:
    The distribution of Chi-square. Proceed. Nation. Academ. Scien., Vol. 17, 684-688

function probq ( x2, ndf, upper, ndf_max )

Purpose

Evaluates the chi-squared distribution function from X2(i) to infinity if UPPER is true or from zero to X2(i) if UPPER is false, for i=1 to size(X2). In other words, if:

  • UPPER = true : PROBQ( i ) = prob( U > X2(i) ) ,
  • UPPER = false : PROBQ( i ) = prob( U < X2(i) ) ,

, for U = Chi-squared(NDF( i )) and i=1 to size(X2).

Arguments

X2 (INPUT) real(stnd), dimension(:)
On entry, upper or lower limits of integration. X2(i) must be greater or equal to zero for i=1 to size(X2).
NDF (INPUT) integer(i4b), dimension(:)

On entry, degrees of freedom of the chi-squared distribution. Any value in the array NDF must be greater or equal to 1.

The size of NDF must verify size(NDF) = size(X2) .

UPPER (INPUT) logical(lgl)

On entry, if:

  • UPPER = true : probability to the right of X2 is calculated.
  • UPPER = false : probability to the left of X2 is calculated.
NDF_MAX (INPUT, OPTIONAL) integer(i4b)

On entry, if:

  • NDF is lower or equal to NDF_MAX, the chi-squared density is integrated.
  • NDF is greater than NDF_MAX, a gaussian approximation is used.

The default is 40.

Further Details

If NDF(i)<=NDF_MAX, the chi-squared distribution function is integrating by using formulae 26.4.4 and 26.4.5 in reference (1), otherwise a normal approximation based on the Wilson-Hilferty transformation is used (see the reference (1), Formula 26.4.14).

This function works only for integer degrees of freedom. It may be faster than PROBQ2 or PROBQ3 functions for the default value of NDF_MAX, but it is less accurate.

The function is parallelized if OPENMP is used.

This function is adapted from:

  1. Abramowitz, M., and Stegun, I.A., 1970:
    Handbook of Mathematical Functions, (formulae 26.4.4, 26.4.5 and 26.4.14). New York, Dover Publications
  2. Wilson, E.B., and Hilferty, M.M., 1931:
    The distribution of Chi-square. Proceed. Nation. Academ. Scien., Vol. 17, 684-688

function probq ( x2, ndf, upper, ndf_max )

Purpose

Evaluates the chi-squared distribution function from X2(i,j) to infinity if UPPER is true or from zero to X2(i,j) if UPPER is false, for i=1 to size(X2,1) and j=1 to size(X2,2). In other words, if:

  • UPPER = true : PROBQ( i, j ) = prob( U > X2(i,j) ) ,
  • UPPER = false : PROBQ( i, j ) = prob( U < X2(i,j) ) ,

, for U = Chi-squared(NDF), i=1 to size(X2,1) and j=1 to size(X2,2).

Arguments

X2 (INPUT) real(stnd), dimension(:,:)
On entry, upper or lower limits of integration. X2(i,j) must be greater or equal to zero for i=1 to size(X2,1) and j=1 to size(X2,2).
NDF (INPUT) integer(i4b)
On entry, degrees of freedom of the chi-squared distribution. NDF must be greater or equal to 1.
UPPER (INPUT) logical(lgl)

On entry, if:

  • UPPER = true : probability to the right of X2 is calculated.
  • UPPER = false : probability to the left of X2 is calculated.
NDF_MAX (INPUT, OPTIONAL) integer(i4b)

On entry, if:

  • NDF is lower or equal to NDF_MAX, the chi-squared density is integrated.
  • NDF is greater than NDF_MAX, a gaussian approximation is used.

The default is 40.

Further Details

If NDF<=NDF_MAX, the chi-squared distribution function is integrating by using formulae 26.4.4 and 26.4.5 in reference (1), otherwise a normal approximation based on the Wilson-Hilferty transformation is used (see the reference (1), Formula 26.4.14).

This function works only for integer degrees of freedom. It may be faster than PROBQ2 or PROBQ3 functions for the default value of NDF_MAX, but it is less accurate.

The function is parallelized if OPENMP is used.

This function is adapted from:

  1. Abramowitz, M., and Stegun, I.A., 1970:
    Handbook of Mathematical Functions, (formulae 26.4.4, 26.4.5 and 26.4.14). New York, Dover Publications
  2. Wilson, E.B., and Hilferty, M.M., 1931:
    The distribution of Chi-square. Proceed. Nation. Academ. Scien., Vol. 17, 684-688

function probq ( x2, ndf, upper, ndf_max )

Purpose

Evaluates the chi-squared distribution function from X2(i,j) to infinity if UPPER is true or from zero to X2(i,j) if UPPER is false, for i=1 to size(X2,1) and j=1 to size(X2,2). In other words, if:

  • UPPER = true : PROBQ( i, j ) = prob( U > X2(i,j) ) ,
  • UPPER = false : PROBQ( i, j ) = prob( U < X2(i,j) ) ,

, for U = Chi-squared(NDF( i, j )), i=1 to size(X2,1) and j=1 to size(X2,2).

Arguments

X2 (INPUT) real(stnd), dimension(:,:)
On entry, upper or lower limits of integration. X2(i,j) must be greater or equal to zero for i=1 to size(X2,1) and j=1 to size(X2,2).
NDF (INPUT) integer(i4b), dimension(:,:)

On entry, degrees of freedom of the chi-squared distribution. Any value in the array NDF must be greater or equal to 1.

The shape of NDF must verify:

  • size(NDF,1) = size(X2,1)
  • size(NDF,2) = size(X2,2) .
UPPER (INPUT) logical(lgl)

On entry, if:

  • UPPER = true : probability to the right of X2 is calculated.
  • UPPER = false : probability to the left of X2 is calculated.
NDF_MAX (INPUT, OPTIONAL) integer(i4b)

On entry, if:

  • NDF is lower or equal to NDF_MAX, the chi-squared density is integrated.
  • NDF is greater than NDF_MAX, a gaussian approximation is used.

The default is 40.

Further Details

If NDF(i,j)<=NDF_MAX, the chi-squared distribution function is integrating by using formulae 26.4.4 and 26.4.5 in reference (1), otherwise a normal approximation based on the Wilson-Hilferty transformation is used (see the reference (1), Formula 26.4.14).

This function works only for integer degrees of freedom. It may be faster than PROBQ2 or PROBQ3 functions for the default value of NDF_MAX, but it is less accurate.

The function is parallelized if OPENMP is used.

This function is adapted from:

  1. Abramowitz, M., and Stegun, I.A., 1970:
    Handbook of Mathematical Functions, (formulae 26.4.4, 26.4.5 and 26.4.14). New York, Dover Publications
  2. Wilson, E.B., and Hilferty, M.M., 1931:
    The distribution of Chi-square. Proceed. Nation. Academ. Scien., Vol. 17, 684-688

function probq2 ( x2, df, upper, df_max, maxiter, failure )

Purpose

Evaluates the chi-squared distribution function from X2 to infinity if UPPER is true or from zero to X2 if UPPER is false. In other words, if:

  • UPPER = true : PROBQ2 = prob( U > X2 ) ,
  • UPPER = false : PROBQ2 = prob( U <= X2 ) ,

, for U = Chi-squared(DF).

PROBQ2 computes the probability that a random variable which follows the chi-squared distribution with DF degrees of freedom is less than or equal to X2 (if UPPER is set to false) or greater than to X2 (if UPPER is set to true).

Arguments

X2 (INPUT) real(stnd)
On entry, upper or lower limit of integration. X2 must be greater or equal to zero.
DF (INPUT) real(stnd)
On entry, degrees of freedom of the chi-squared distribution. DF must be greater or equal to 0.5 and less than or equal to 200 000.
UPPER (INPUT) logical(lgl)

On entry, if:

  • UPPER = true : probability to the right of X2 is calculated.
  • UPPER = false : probability to the left of X2 is calculated.
DF_MAX (INPUT, OPTIONAL) real(stnd)

On entry, if:

  • DF is lower or equal to DF_MAX, the chi-squared density is integrated using the incomplete Gamma integral.
  • DF is greater than DF_MAX, a gaussian approximation is used.

The default is 100.

MAXITER (INPUT, OPTIONAL) integer(i4b)

On entry, MAXITER controls the maximum number of iterations when evaluating the Pearson’s series expansion of the incomplete Gamma integral if DF<=DF_MAX.

The default value is 1000.

FAILURE (INPUT, OPTIONAL) logical(lgl)

On entry, if FAILURE is set to true, the values for which the “integrating” process of the incomplete Gamma integral did not converge to the desired accuracy are set to -1. The algorithm fails to converge if the number of iterations exceeds MAXITER.

This argument is actually used only if DF<=DF_MAX.

The default value is false.

Further Details

If DF<=DF_MAX, the chi-squared distribution function is evaluated by integrating the incomplete Gamma integral (see the references (2) and (3) for more details), otherwise a normal approximation based on the Wilson-Hilferty transformation is used (see the reference (3), Formula 26.4.14, p.941).

This function is faster than PROBQ3 function, but is less accurate.

This function is adapted from:

  1. Wilson, E.B., and Hilferty, M.M., 1931:
    The distribution of Chi-square. Proceed. Nation. Academ. Scien., Vol. 17, 684-688
  2. Shea, B.L., 1988:
    Algorithm AS 239: Chi-Squared and incomplete Gamma integral. Appl. Statist., Vol. 37, No. 3, pp. 466-473
  3. Abramowitz, M., and Stegun, I.A., 1970:
    Handbook of Mathematical Functions, (formulae 6.5.29, 6.5.32 and 26.4.14). New York, Dover Publications

function probq2 ( x2, df, upper, df_max, maxiter, failure )

Purpose

Evaluates the chi-squared distribution function from X2(i) to infinity if UPPER is true or from zero to X2(i) if UPPER is false, for i=1 to size(X2). In other words, if:

  • UPPER = true : PROBQ2( i ) = prob( U > X2(i) ) ,
  • UPPER = false : PROBQ2( i ) = prob( U <= X2(i) ) ,

, for U = Chi-squared(DF) and i=1 to size(X2).

PROBQ2 computes the probabilities that a random variable (vector) which follows the chi-squared distribution with DF degrees of freedom is less than or equal to X2(:) (if UPPER is set to false) or greater than to X2(:) (if UPPER is set to true).

Arguments

X2 (INPUT) real(stnd), dimension(:)
On entry, upper or lower limits of integration. X2(i) must be greater or equal to zero for i=1 to size(X2).
DF (INPUT) real(stnd)
On entry, degrees of freedom of the chi-squared distribution. DF must be greater or equal to 0.5 and less than or equal to 200 000.
UPPER (INPUT) logical(lgl)

On entry, if:

  • UPPER = true : probability to the right of X2 is calculated.
  • UPPER = false : probability to the left of X2 is calculated.
DF_MAX (INPUT, OPTIONAL) real(stnd)

On entry, if:

  • DF is lower or equal to DF_MAX, the chi-squared density is integrated using the incomplete Gamma integral.
  • DF is greater than DF_MAX, a gaussian approximation is used.

The default is 100.

MAXITER (INPUT, OPTIONAL) integer(i4b)

On entry, MAXITER controls the maximum number of iterations when evaluating the Pearson’s series expansion of the incomplete Gamma integral if DF<=DF_MAX.

The default value is 1000.

FAILURE (INPUT, OPTIONAL) logical(lgl)

On entry, if FAILURE is set to true, the values for which the “integrating” process of the incomplete Gamma integral did not converge to the desired accuracy are set to -1. The algorithm fails to converge if the number of iterations exceeds MAXITER.

This argument is actually used only if DF<=DF_MAX.

The default value is false.

Further Details

If DF<=DF_MAX, the chi-squared distribution function is evaluated by integrating the incomplete Gamma integral (see the references (2) and (3) for more details), otherwise a normal approximation based on the Wilson-Hilferty transformation is used (see the reference (3), Formula 26.4.14, p.941).

This function is faster than PROBQ3 function, but is less accurate.

This functon is parallelized if OPENMP is used.

This function is adapted from:

  1. Wilson, E.B., and Hilferty, M.M., 1931:
    The distribution of Chi-square. Proceed. Nation. Academ. Scien., Vol. 17, 684-688
  2. Shea, B.L., 1988:
    Algorithm AS 239: Chi-Squared and incomplete Gamma integral. Appl. Statist., Vol. 37, No. 3, pp. 466-473
  3. Abramowitz, M., and Stegun, I.A., 1970:
    Handbook of Mathematical Functions, (formulae 6.5.29, 6.5.32 and 26.4.14). New York, Dover Publications

function probq2 ( x2, df, upper, df_max, maxiter, failure )

Purpose

Evaluates the chi-squared distribution function from X2(i) to infinity if UPPER is true or from zero to X2(i) if UPPER is false, for i=1 to size(X2). In other words, if:

  • UPPER = true : PROBQ2( i ) = prob( U > X2(i) ) ,
  • UPPER = false : PROBQ2( i ) = prob( U <= X2(i) ) ,

, for U = Chi-squared(DF( i )) and i=1 to size(X2).

PROBQ2 computes the probabilities that a random variable (vector) which follows the chi-squared distribution with DF(:) degrees of freedom is less than or equal to X2(:) (if UPPER is set to false) or greater than to X2(:) (if UPPER is set to true).

Arguments

X2 (INPUT) real(stnd), dimension(:)
On entry, upper or lower limits of integration. X2(i) must be greater or equal to zero for i=1 to size(X2).
DF (INPUT) real(stnd), dimension(:)

On entry, degrees of freedom of the chi-squared distribution. Any value in the array DF must be greater or equal to 0.5 and less than or equal to 200 000.

The size of DF must verify size(DF) = size(X2).

UPPER (INPUT) logical(lgl)

On entry, if:

  • UPPER = true : probability to the right of X2 is calculated.
  • UPPER = false : probability to the left of X2 is calculated.
DF_MAX (INPUT, OPTIONAL) real(stnd)

On entry, if:

  • DF is lower or equal to DF_MAX, the chi-squared density is integrated using the incomplete Gamma integral.
  • DF is greater than DF_MAX, a gaussian approximation is used.

The default is 100.

MAXITER (INPUT, OPTIONAL) integer(i4b)

On entry, MAXITER controls the maximum number of iterations when evaluating the Pearson’s series expansion of the incomplete Gamma integral if DF<=DF_MAX.

The default value is 1000.

FAILURE (INPUT, OPTIONAL) logical(lgl)

On entry, if FAILURE is set to true, the values for which the “integrating” process of the incomplete Gamma integral did not converge to the desired accuracy are set to -1. The algorithm fails to converge if the number of iterations exceeds MAXITER.

This argument is actually used only for the values of DF(:) less than or equal to DF_MAX.

The default value is false.

Further Details

If DF(i)<=DF_MAX, the chi-squared distribution function is evaluated by integrating the incomplete Gamma integral (see the references (2) and (3) for more details), otherwise a normal approximation based on the Wilson-Hilferty transformation is used (see the reference (3), Formula 26.4.14, p.941).

This function is faster than PROBQ3 function, but is less accurate.

This function is parallelized if OPENMP is used.

This function is adapted from:

  1. Wilson, E.B., and Hilferty, M.M., 1931:
    The distribution of Chi-square. Proceed. Nation. Academ. Scien., Vol. 17, 684-688
  2. Shea, B.L., 1988:
    Algorithm AS 239: Chi-Squared and incomplete Gamma integral. Appl. Statist., Vol. 37, No. 3, pp. 466-473
  3. Abramowitz, M., and Stegun, I.A., 1970:
    Handbook of Mathematical Functions, (formulae 6.5.29, 6.5.32 and 26.4.14). New York, Dover Publications

function probq2 ( x2, df, upper, df_max, maxiter, failure )

Purpose

Evaluates the chi-squared distribution function from X2(i,j) to infinity if UPPER is true or from zero to X2(i,j) if UPPER is false, for i=1 to size(X2,1) and j=1 to size(X2,2). In other words, if:

  • UPPER = true : PROBQ2( i,j ) = prob( U > X2(i,j) ) ,
  • UPPER = false : PROBQ2( i,j ) = prob( U <= X2(i,j) ) ,

, for U = Chi-squared(DF), i=1 to size(X2,1) and j=1 to size(X2,2).

PROBQ2 computes the probabilities that a random variable (matrix) which follows the chi-squared distribution with DF degrees of freedom is less than or equal to X2(:,:) (if UPPER is set to false) or greater than to X2(:,:) (if UPPER is set to true).

Arguments

X2 (INPUT) real(stnd), dimension(:,:)
On entry, upper or lower limits of integration. X2(i,j) must be greater or equal to zero for i=1 to size(X2,1) and j=1 to size(X2,2).
DF (INPUT) real(stnd)
On entry, degrees of freedom of the chi-squared distribution. DF must be greater or equal to 0.5 and less than or equal to 200 000.
UPPER (INPUT) logical(lgl)

On entry, if:

  • UPPER = true : probability to the right of X2 is calculated.
  • UPPER = false : probability to the left of X2 is calculated.
DF_MAX (INPUT, OPTIONAL) real(stnd)

On entry, if:

  • DF is lower or equal to DF_MAX, the chi-squared density is integrated using the incomplete Gamma integral.
  • DF is greater than DF_MAX, a gaussian approximation is used.

The default is 100.

MAXITER (INPUT, OPTIONAL) integer(i4b)

On entry, MAXITER controls the maximum number of iterations when evaluating the Pearson’s series expansion of the incomplete Gamma integral if DF<=DF_MAX.

The default value is 1000.

FAILURE (INPUT, OPTIONAL) logical(lgl)

On entry, if FAILURE is set to true, the values for which the “integrating” process of the incomplete Gamma integral did not converge to the desired accuracy are set to -1. The algorithm fails to converge if the number of iterations exceeds MAXITER.

This argument is actually used only if DF<=DF_MAX.

The default value is false.

Further Details

If DF<=DF_MAX, the chi-squared distribution function is evaluated by integrating the incomplete Gamma integral (see the references (2) and (3) for more details), otherwise a normal approximation based on the Wilson-Hilferty transformation is used (see the reference (3), Formula 26.4.14, p.941).

This function is faster than PROBQ3 function, but is less accurate.

This function is parallelized if OPENMP is used.

This function is adapted from:

  1. Wilson, E.B., and Hilferty, M.M., 1931:
    The distribution of Chi-square. Proceed. Nation. Academ. Scien., Vol. 17, 684-688
  2. Shea, B.L., 1988:
    Algorithm AS 239: Chi-Squared and incomplete Gamma integral. Appl. Statist., Vol. 37, No. 3, pp. 466-473
  3. Abramowitz, M., and Stegun, I.A., 1970:
    Handbook of Mathematical Functions, (formulae 6.5.29, 6.5.32 and 26.4.14). New York, Dover Publications

function probq2 ( x2, df, upper, df_max, maxiter, failure )

Purpose

Evaluates the chi-squared distribution function from X2(i,j) to infinity if UPPER is true or from zero to X2(i,j) if UPPER is false, for i=1 to size(X2,1) and j=1 to size(X2,2). In other words, if:

  • UPPER = true : PROBQ2( i,j ) = prob( U > X2(i,j) ) ,
  • UPPER = false : PROBQ2( i,j ) = prob( U <= X2(i,j) ) ,

, for U = Chi-squared(DF( i,j )), i=1 to size(X2,1) and j=1 to size(X2,2).

PROBQ2 computes the probabilities that a random variable (matrix) which follows the chi-squared distribution with DF(:,:) degrees of freedom is less than or equal to X2(:,:) (if UPPER is set to false) or greater than to X2(:,:) (if UPPER is set to true).

Arguments

X2 (INPUT) real(stnd), dimension(:,:)
On entry, upper or lower limits of integration. X2(i,j) must be greater or equal to zero for i=1 to size(X2,1) and j=1 to size(X2,2).
DF (INPUT) real(stnd), dimension(:,:)

On entry, degrees of freedom of the chi-squared distribution. Any value in the array DF must be greater or equal to 0.5 and less than or equal to 200 000.

The shape of DF must verify:

  • size(DF,1) = size(X2,1)
  • size(DF,2) = size(X2,2).
UPPER (INPUT) logical(lgl)

On entry, if:

  • UPPER = true : probability to the right of X2 is calculated.
  • UPPER = false : probability to the left of X2 is calculated.
DF_MAX (INPUT, OPTIONAL) real(stnd)

On entry, if:

  • DF is lower or equal to DF_MAX, the chi-squared density is integrated using the incomplete Gamma integral.
  • DF is greater than DF_MAX, a gaussian approximation is used.

The default is 100.

MAXITER (INPUT, OPTIONAL) integer(i4b)

On entry, MAXITER controls the maximum number of iterations when evaluating the Pearson’s series expansion of the incomplete Gamma integral if DF<=DF_MAX.

The default value is 1000.

FAILURE (INPUT, OPTIONAL) logical(lgl)

On entry, if FAILURE is set to true, the values for which the “integrating” process of the incomplete Gamma integral did not converge to the desired accuracy are set to -1. The algorithm fails to converge if the number of iterations exceeds MAXITER.

This argument is actually used only for the values of DF(:,:) less than or equal to DF_MAX.

The default value is false.

Further Details

If DF(i,j)<=DF_MAX, the chi-squared distribution function is evaluated by integrating the incomplete Gamma integral (see the references (2) and (3) for more details), otherwise a normal approximation based on the Wilson-Hilferty transformation is used (see the reference (3), Formula 26.4.14, p.941).

This function is faster than PROBQ3 function, but is less accurate.

This function is parallelized if OPENMP is used.

This function is adapted from:

  1. Wilson, E.B., and Hilferty, M.M., 1931:
    The distribution of Chi-square. Proceed. Nation. Academ. Scien., Vol. 17, 684-688
  2. Shea, B.L., 1988:
    Algorithm AS 239: Chi-Squared and incomplete Gamma integral. Appl. Statist., Vol. 37, No. 3, pp. 466-473
  3. Abramowitz, M., and Stegun, I.A., 1970:
    Handbook of Mathematical Functions, (formulae 6.5.29, 6.5.32 and 26.4.14). New York, Dover Publications

function probq3 ( x2, df, upper, df_max, acu, maxiter, failure )

Purpose

Evaluates the chi-squared distribution function from X2 to infinity if UPPER is true or from zero to X2 if UPPER is false. In other words, if:

  • UPPER = true : PROBQ3 = prob( U > X2 ) ,
  • UPPER = false : PROBQ3 = prob( U <= X2 ) ,

, for U = Chi-squared(DF).

PROBQ3 computes the probability that a random variable which follows the chi-squared distribution with DF degrees of freedom is less than or equal to X2 (if UPPER is set to false) or greater than to X2 (if UPPER is set to true).

Arguments

X2 (INPUT) real(stnd)
On entry, upper or lower limit of integration. X2 must be greater or equal to zero.
DF (INPUT) real(stnd)
On entry, degrees of freedom of the chi-squared distribution. DF must be greater or equal to 0.5.
UPPER (INPUT) logical(lgl)

On entry, if:

  • UPPER = true : probability to the right of X2 is calculated.
  • UPPER = false : probability to the left of X2 is calculated.
DF_MAX (INPUT, OPTIONAL) real(stnd)

On entry, if:

  • DF is lower or equal to DF_MAX, the chi-squared density is integrated using the incomplete Gamma integral.
  • DF is greater than DF_MAX, a gaussian approximation is used.

The default is 100.

ACU (INPUT, OPTIONAL) real(stnd)

On entry, the desired accuracy of the result if DF<=DF_MAX (e.g. if the incomplete Gamma integral is used). If l decimal places of accuracy are required then ACU should be set to 10**(-(l+1)). ACU is a small strictly positive integer. ACU should not be set smaller than the machine precision since the stated accuracy cannot be attained. In that case the machine precision is used instead.

The default value for ACU is epsilon( ACU ).

MAXITER (INPUT, OPTIONAL) integer(i4b)

On entry, MAXITER controls the maximum number of iterations when evaluating the Pearson’s series or continued fraction expansion of the incomplete Gamma integral if DF<=DF_MAX.

The default value is 1000.

FAILURE (INPUT, OPTIONAL) logical(lgl)

On entry, if FAILURE is set to true, the values for which the “integrating” process of the incomplete Gamma integral did not converge to the desired accuracy are set to -1. The algorithm fails to converge if the number of iterations exceeds MAXITER.

This argument is actually used only if DF<=DF_MAX.

The default value is false.

Further Details

If DF<=DF_MAX, the chi-squared distribution function is evaluated by integrating the incomplete Gamma integral (see the references (2) and (3) for more details), otherwise a normal approximation based on the Wilson-Hilferty transformation is used (see the reference (3), Formula 26.4.14, p.941).

This function is more accurate than PROBQ and PROBQ2 functions, but it is slower.

This function is adapted from:

  1. Wilson, E.B., and Hilferty, M.M., 1931:
    The distribution of Chi-square. Proceed. Nation. Academ. Scien., Vol. 17, 684-688
  2. Shea, B.L., 1988:
    Algorithm AS 239: Chi-Squared and incomplete Gamma integral. Appl. Statist., Vol. 37, No. 3, pp. 466-473
  3. Abramowitz, M., and Stegun, I.A., 1970:
    Handbook of Mathematical Functions, (formulae 6.5.29, 6.5.31 and 26.4.14). New York, Dover Publications

function probq3 ( x2, df, upper, df_max, acu, maxiter, failure )

Purpose

Evaluates the chi-squared distribution function from X2(i) to infinity if UPPER is true or from zero to X2(i) if UPPER is false, for i=1 to size(X2). In other words, if:

  • UPPER = true : PROBQ3( i ) = prob( U > X2(i) ) ,
  • UPPER = false : PROBQ3( i ) = prob( U <= X2(i) ) ,

, for U = Chi-squared(DF) and i=1 to size(X2).

PROBQ3 computes the probabilities that a random variable (vector) which follows the chi-squared distribution with DF degrees of freedom is less than or equal to X2(:) (if UPPER is set to false) or greater than to X2(:) (if UPPER is set to true).

Arguments

X2 (INPUT) real(stnd), dimension(:)
On entry, upper or lower limits of integration. X2(i) must be greater or equal to zero for i=1 to size(X2).
DF (INPUT) real(stnd)
On entry, degrees of freedom of the chi-squared distribution. DF must be greater or equal to 0.5.
UPPER (INPUT) logical(lgl)

On entry, if:

  • UPPER = true : probability to the right of X2 is calculated.
  • UPPER = false : probability to the left of X2 is calculated.
DF_MAX (INPUT, OPTIONAL) real(stnd)

On entry, if:

  • DF is lower or equal to DF_MAX, the chi-squared density is integrated using the incomplete Gamma integral.
  • DF is greater than DF_MAX, a gaussian approximation is used.

The default is 100.

ACU (INPUT, OPTIONAL) real(stnd)

On entry, the desired accuracy of the result if DF<=DF_MAX (e.g. if the incomplete Gamma integral is used). If l decimal places of accuracy are required then ACU should be set to 10**(-(l+1)). ACU is a small strictly positive integer. ACU should not be set smaller than the machine precision since the stated accuracy cannot be attained. In that case the machine precision is used instead.

The default value for ACU is epsilon( ACU ).

MAXITER (INPUT, OPTIONAL) integer(i4b)

On entry, MAXITER controls the maximum number of iterations when evaluating the Pearson’s series or continued fraction expansion of the incomplete Gamma integral if DF<=DF_MAX.

The default value is 1000.

FAILURE (INPUT, OPTIONAL) logical(lgl)

On entry, if FAILURE is set to true, the values for which the “integrating” process of the incomplete Gamma integral did not converge to the desired accuracy are set to -1. The algorithm fails to converge if the number of iterations exceeds MAXITER.

This argument is actually used only if DF<=DF_MAX.

The default value is false.

Further Details

If DF<=DF_MAX, the chi-squared distribution function is evaluated by integrating the incomplete Gamma integral (see the references (2) and (3) for more details), otherwise a normal approximation based on the Wilson-Hilferty transformation is used (see the reference (3), Formula 26.4.14, p.941).

This function is more accurate than PROBQ and PROBQ2 functions, but it is slower.

The function is parallelized if OPENMP is used.

This function is adapted from:

  1. Wilson, E.B., and Hilferty, M.M., 1931:
    The distribution of Chi-square. Proceed. Nation. Academ. Scien., Vol. 17, 684-688
  2. Shea, B.L., 1988:
    Algorithm AS 239: Chi-Squared and incomplete Gamma integral. Appl. Statist., Vol. 37, No. 3, pp. 466-473
  3. Abramowitz, M., and Stegun, I.A., 1970:
    Handbook of Mathematical Functions, (formulae 6.5.29, 6.5.31 and 26.4.14). New York, Dover Publications

function probq3 ( x2, df, upper, df_max, acu, maxiter, failure )

Purpose

Evaluates the chi-squared distribution function from X2(i) to infinity if UPPER is true or from zero to X2(i) if UPPER is false, for i=1 to size(X2). In other words, if:

  • UPPER = true : PROBQ3( i ) = prob( U > X2(i) ) ,
  • UPPER = false : PROBQ3( i ) = prob( U <= X2(i) ) ,

, for U = Chi-squared(NDF( i )) and i=1 to size(X2).

PROBQ3 computes the probabilities that a random variable (vector) which follows the chi-squared distribution with DF(:) degrees of freedom is less than or equal to X2(:) (if UPPER is set to false) or greater than to X2(:) (if UPPER is set to true).

Arguments

X2 (INPUT) real(stnd), dimension(:)
On entry, upper or lower limits of integration. X2(i) must be greater or equal to zero for i=1 to size(X2).
DF (INPUT) real(stnd), dimension(:)

On entry, degrees of freedom of the chi-squared distribution. Any value in the array DF must be greater or equal to 0.5.

The size of NDF must verify size(DF) = size(X2) .

UPPER (INPUT) logical(lgl)

On entry, if:

  • UPPER = true : probability to the right of X2 is calculated.
  • UPPER = false : probability to the left of X2 is calculated.
DF_MAX (INPUT, OPTIONAL) real(stnd)

On entry, if:

  • DF is lower or equal to DF_MAX, the chi-squared density is integrated using the incomplete Gamma integral.
  • DF is greater than DF_MAX, a gaussian approximation is used.

The default is 100.

ACU (INPUT, OPTIONAL) real(stnd)

On entry, the desired accuracy of the result if DF<=DF_MAX (e.g. if the incomplete Gamma integral is used). If l decimal places of accuracy are required then ACU should be set to 10**(-(l+1)). ACU is a small strictly positive integer. ACU should not be set smaller than the machine precision since the stated accuracy cannot be attained. In that case the machine precision is used instead.

The default value for ACU is epsilon( ACU ).

MAXITER (INPUT, OPTIONAL) integer(i4b)

On entry, MAXITER controls the maximum number of iterations when evaluating the Pearson’s series or continued fraction expansion of the incomplete Gamma integral if DF<=DF_MAX.

The default value is 1000.

FAILURE (INPUT, OPTIONAL) logical(lgl)

On entry, if FAILURE is set to true, the values for which the “integrating” process of the incomplete Gamma integral did not converge to the desired accuracy are set to -1. The algorithm fails to converge if the number of iterations exceeds MAXITER.

This argument is actually used only if DF<=DF_MAX.

The default value is false.

Further Details

If DF(i)<=DF_MAX, the chi-squared distribution function is evaluated by integrating the incomplete Gamma integral (see the references (2) and (3) for more details), otherwise a normal approximation based on the Wilson-Hilferty transformation is used (see the reference (3), Formula 26.4.14, p.941).

This function is more accurate than PROBQ and PROBQ2 functions, but it is slower.

The function is parallelized if OPENMP is used.

This function is adapted from:

  1. Wilson, E.B., and Hilferty, M.M., 1931:
    The distribution of Chi-square. Proceed. Nation. Academ. Scien., Vol. 17, 684-688
  2. Shea, B.L., 1988:
    Algorithm AS 239: Chi-Squared and incomplete Gamma integral. Appl. Statist., Vol. 37, No. 3, pp. 466-473
  3. Abramowitz, M., and Stegun, I.A., 1970:
    Handbook of Mathematical Functions, (formulae 6.5.29, 6.5.31 and 26.4.14). New York, Dover Publications

function probq3 ( x2, df, upper, df_max, acu, maxiter, failure )

Purpose

Evaluates the chi-squared distribution function from X2(i,j) to infinity if UPPER is true or from zero to X2(i,j) if UPPER is false, for i=1 to size(X2,1) and j=1 to size(X2,2). In other words, if:

  • UPPER = true : PROBQ3( i,j ) = prob( U > X2(i,j) ) ,
  • UPPER = false : PROBQ3( i,j ) = prob( U <= X2(i,j) ) ,

, for U = Chi-squared(DF), i=1 to size(X2,1) and j=1 to size(X2,2).

PROBQ3 computes the probabilities that a random variable (matrix) which follows the chi-squared distribution with DF degrees of freedom is less than or equal to X2(:,:) (if UPPER is set to false) or greater than to X2(:,:) (if UPPER is set to true).

Arguments

X2 (INPUT) real(stnd), dimension(:,:)
On entry, upper or lower limits of integration. X2(i,j) must be greater or equal to zero for i=1 to size(X2,1) and j=1 to size(X2,2).
DF (INPUT) real(stnd)
On entry, degrees of freedom of the chi-squared distribution. DF must be greater or equal to 0.5.
UPPER (INPUT) logical(lgl)

On entry, if:

  • UPPER = true : probability to the right of X2 is calculated.
  • UPPER = false : probability to the left of X2 is calculated.
DF_MAX (INPUT, OPTIONAL) real(stnd)

On entry, if:

  • DF is lower or equal to DF_MAX, the chi-squared density is integrated using the incomplete Gamma integral.
  • DF is greater than DF_MAX, a gaussian approximation is used.

The default is 100.

ACU (INPUT, OPTIONAL) real(stnd)

On entry, the desired accuracy of the result if DF<=DF_MAX (e.g. if the incomplete Gamma integral is used). If l decimal places of accuracy are required then ACU should be set to 10**(-(l+1)). ACU is a small strictly positive integer. ACU should not be set smaller than the machine precision since the stated accuracy cannot be attained. In that case the machine precision is used instead.

The default value for ACU is epsilon( ACU ).

MAXITER (INPUT, OPTIONAL) integer(i4b)

On entry, MAXITER controls the maximum number of iterations when evaluating the Pearson’s series or continued fraction expansion of the incomplete Gamma integral if DF<=DF_MAX.

The default value is 1000.

FAILURE (INPUT, OPTIONAL) logical(lgl)

On entry, if FAILURE is set to true, the values for which the “integrating” process of the incomplete Gamma integral did not converge to the desired accuracy are set to -1. The algorithm fails to converge if the number of iterations exceeds MAXITER.

This argument is actually used only if DF<=DF_MAX.

The default value is false.

Further Details

If DF<=DF_MAX, the chi-squared distribution function is evaluated by integrating the incomplete Gamma integral (see the references (2) and (3) for more details), otherwise a normal approximation based on the Wilson-Hilferty transformation is used (see the reference (3), Formula 26.4.14, p.941).

This function is more accurate than PROBQ and PROBQ2 functions, but it is slower.

The function is parallelized if OPENMP is used.

This function is adapted from:

  1. Wilson, E.B., and Hilferty, M.M., 1931:
    The distribution of Chi-square. Proceed. Nation. Academ. Scien., Vol. 17, 684-688
  2. Shea, B.L., 1988:
    Algorithm AS 239: Chi-Squared and incomplete Gamma integral. Appl. Statist., Vol. 37, No. 3, pp. 466-473
  3. Abramowitz, M., and Stegun, I.A., 1970:
    Handbook of Mathematical Functions, (formulae 6.5.29, 6.5.31 and 26.4.14). New York, Dover Publications

function probq3 ( x2, df, upper, df_max, acu, maxiter, failure )

Purpose

Evaluates the chi-squared distribution function from X2(i,j) to infinity if UPPER is true or from zero to X2(i,j) if UPPER is false, for i=1 to size(X2,1) and j=1 to size(X2,2). In other words, if:

  • UPPER = true : PROBQ3( i, j ) = prob( U > X2(i,j) ) ,
  • UPPER = false : PROBQ3( i, j ) = prob( U < X2(i,j) ) ,

, for U = Chi-squared(NDF( i, j )), i=1 to size(X2,1) and j=1 to size(X2,2).

PROBQ3 computes the probabilities that a random variable (matrix) which follows the chi-squared distribution with DF(:,:) degrees of freedom is less than or equal to X2(:,:) (if UPPER is set to false) or greater than to X2(:,:) (if UPPER is set to true).

Arguments

X2 (INPUT) real(stnd), dimension(:,:)
On entry, upper or lower limits of integration. X2(i,j) must be greater or equal to zero for i=1 to size(X2,1) and j=1 to size(X2,2).
DF (INPUT) real(stnd), dimension(:,:)

On entry, degrees of freedom of the chi-squared distribution. Any value in the array DF must be greater or equal to 0.5.

The shape of DF must verify:

  • size(DF,1) = size(X2,1)
  • size(DF,2) = size(X2,2) .
UPPER (INPUT) logical(lgl)

On entry, if:

  • UPPER = true : probability to the right of X2 is calculated.
  • UPPER = false : probability to the left of X2 is calculated.
DF_MAX (INPUT, OPTIONAL) real(stnd)

On entry, if:

  • DF is lower or equal to DF_MAX, the chi-squared density is integrated using the incomplete Gamma integral.
  • DF is greater than DF_MAX, a gaussian approximation is used.

The default is 100.

ACU (INPUT, OPTIONAL) real(stnd)

On entry, the desired accuracy of the result if DF<=DF_MAX (e.g. if the incomplete Gamma integral is used). If l decimal places of accuracy are required then ACU should be set to 10**(-(l+1)). ACU is a small strictly positive integer. ACU should not be set smaller than the machine precision since the stated accuracy cannot be attained. In that case the machine precision is used instead.

The default value for ACU is epsilon( ACU ).

MAXITER (INPUT, OPTIONAL) integer(i4b)

On entry, MAXITER controls the maximum number of iterations when evaluating the Pearson’s series or continued fraction expansion of the incomplete Gamma integral if DF<=DF_MAX.

The default value is 1000.

FAILURE (INPUT, OPTIONAL) logical(lgl)

On entry, if FAILURE is set to true, the values for which the “integrating” process of the incomplete Gamma integral did not converge to the desired accuracy are set to -1. The algorithm fails to converge if the number of iterations exceeds MAXITER.

This argument is actually used only if DF<=DF_MAX.

The default value is false.

Further Details

If DF(i,j)<=DF_MAX, the chi-squared distribution function is evaluated by integrating the incomplete Gamma integral (see the references (2) and (3) for more details), otherwise a normal approximation based on the Wilson-Hilferty transformation is used (see the reference (3), Formula 26.4.14, p.941).

This function is more accurate than PROBQ and PROBQ2 functions, but it is slower.

The function is parallelized if OPENMP is used.

This function is adapted from:

  1. Wilson, E.B., and Hilferty, M.M., 1931:
    The distribution of Chi-square. Proceed. Nation. Academ. Scien., Vol. 17, 684-688
  2. Shea, B.L., 1988:
    Algorithm AS 239: Chi-Squared and incomplete Gamma integral. Appl. Statist., Vol. 37, No. 3, pp. 466-473
  3. Abramowitz, M., and Stegun, I.A., 1970:
    Handbook of Mathematical Functions, (formulae 6.5.29, 6.5.31 and 26.4.14). New York, Dover Publications

function pinvq ( p, ndf )

Purpose

Evaluates the inverse of the chi-squared distribution function:

X0 = PINVQ( P, NDF )

, if P = probability( U < X0 ) for U = Chi-squared(NDF).

PINVQ returns the quantile X0 of the chi-squared distribution with NDF degrees of freedom corresponding to a given lower tail area of P.

Arguments

P (INPUT) real(stnd)
On entry, the probability. P must verify 0. < P < 1. .
NDF (INPUT) integer(i4b)
On entry, degrees of freedom of the chi-squared distribution. NDF must be greater or equal to 1.

Further Details

This function is fast, but not very accurate especially for small degrees of freedom, e.g. for NDF<10 or 20. If high accuracy is desired, function PINVQ2 must be used instead.

This function is adapted from:

  1. Goldstein, R.B., 1973:
    Chi-square quantiles. Comm. A.C.M., vol.16, no.8, 483-485

function pinvq ( p, ndf )

Purpose

Evaluates the inverse of the chi-squared distribution function:

X0(i) = PINVQ( P(i), NDF )

, if P(i) = probability( U < X0(i) ) for U = Chi-squared(NDF) and i=1 to size(P).

PINVQ returns the quantiles X0(i) of the chi-squared distribution with NDF degrees of freedom corresponding to a given lower tail area of P(i), for i=1 to size(P).

Arguments

P (INPUT) real(stnd), dimension(:)
On entry, the probabilities. P(i) must verify 0. < P(i) < 1. , for i=1 to size(P).
NDF (INPUT) integer(i4b)
On entry, degrees of freedom of the chi-squared distribution. NDF must be greater or equal to 1.

Further Details

This function is fast, but not very accurate especially for small degrees of freedom, e.g. for NDF<10 or 20. If high accuracy is desired, function PINVQ2 must be used instead.

This function is parallelized if OPENMP is used.

This function is adapted from:

  1. Goldstein, R.B., 1973:
    Chi-square quantiles. Comm. A.C.M., vol.16, no.8, 483-485

function pinvq ( p, ndf )

Purpose

Evaluates the inverse of the chi-squared distribution function:

X0(i) = PINVQ( P(i), NDF )

, if P(i) = probability( U < X0(i) ) for U = Chi-squared(NDF(i)) and i=1 to size(P).

PINVQ returns the quantiles X0(i) of the chi-squared distribution with NDF(i) degrees of freedom corresponding to a given lower tail area of P(i), for i=1 to size(P).

Arguments

P (INPUT) real(stnd), dimension(:)
On entry, the probabilities. P(i) must verify 0. < P(i) < 1. , for i=1 to size(P).
NDF (INPUT) integer(i4b), dimension(:)

On entry, degrees of freedom of the chi-squared distribution. Any value in the array NDF must be greater or equal to 1.

The size of NDF must be size(NDF) = size(P).

Further Details

This function is fast, but not very accurate especially for small degrees of freedom, e.g. for NDF<10 or 20. If high accuracy is desired, function PINVQ2 must be used instead.

This function is parallelized if OPENMP is used.

This function is adapted from:

  1. Goldstein, R.B., 1973:
    Chi-square quantiles. Comm. A.C.M., vol.16, no.8, 483-485

function pinvq ( p, ndf )

Purpose

Evaluates the inverse of the chi-squared distribution function:

X0(i,j) = PINVQ( P(i,j), NDF )

, if P(i,j) = probability( U < T0(i,j) ) for U = Chi-squared(NDF), i=1 to size(P,1) and j=1 to size(P,2).

PINVQ returns the quantiles X0(i,j) of the chi-squared distribution with NDF degrees of freedom corresponding to a given lower tail area of P(i,j) for i=1 to size(P,1) and j=1 to size(P,2).

Arguments

P (INPUT) real(stnd), dimension(:,:)
On entry, the probabilities. P(i,j) must verify 0. < P(i,j) < 1. , for i=1 to size(P,1) and j=1 to size(P,2).
NDF (INPUT) integer(i4b)
On entry, degrees of freedom of the chi-squared distribution. NDF must be greater or equal to 1.

Further Details

This function is fast, but not very accurate especially for small degrees of freedom, e.g. for NDF<10 or 20. If high accuracy is desired, function PINVQ2 must be used instead.

This function is parallelized if OPENMP is used.

This function is adapted from:

  1. Goldstein, R.B., 1973:
    Chi-square quantiles. Comm. A.C.M., vol.16, no.8, 483-485

function pinvq ( p, ndf )

Purpose

Evaluates the inverse of the chi-squared distribution function:

X0(i,j) = PINVQ( P(i,j), NDF )

, if P(i,j) = probability( U < T0(i,j) ) for U = Chi-squared(NDF(i,j)), i=1 to size(P,1) and j=1 to size(P,2).

PINVQ returns the quantiles X0(i,j) of the chi-squared distribution with NDF(i,j) degrees of freedom corresponding to a given lower tail area of P(i,j) for i=1 to size(P,1) and j=1 to size(P,2).

Arguments

P (INPUT) real(stnd), dimension(:,:)
On entry, the probabilities. P(i,j) must verify 0. < P(i,j) < 1. , for i=1 to size(P,1) and j=1 to size(P,2).
NDF (INPUT) integer(i4b), dimension(:,:)

On entry, degrees of freedom of the chi-squared distribution. Any value in the array NDF must be greater or equal to 1.

The shape of NDF must verify:

  • size(NDF,1) = size(P,1)
  • size(NDF,2) = size(P,2).

Further Details

This function is fast, but not very accurate especially for small degrees of freedom, e.g. for NDF<10 or 20. If high accuracy is desired, function PINVQ2 must be used instead.

This function is parallelized if OPENMP is used.

This function is adapted from:

  1. Goldstein, R.B., 1973:
    Chi-square quantiles. Comm. A.C.M., vol.16, no.8, 483-485

function pinvq2 ( p, df, prec, acu, maxiter )

Purpose

Evaluates the inverse of the chi-squared distribution function:

X0 = PINVQ2( P, DF )

, if P = probability( U < X0 ) for U = Chi-squared(DF).

PINVQ2 returns the quantile X0 of the chi-squared distribution with DF degrees of freedom corresponding to a given lower tail area of P. In other words, PINVQ2 outputs a chi-squared value, X0, such that a random variable, distributed as chi-squared with DF degrees of freedom, will be less than or equal to X0 with probability P.

Arguments

P (INPUT) real(stnd)
On entry, the probability. P must be in the inclusive range (0,1).
DF (INPUT) real(stnd)
On entry, degrees of freedom of the chi-squared distribution. DF must be greater or equal to 0.5.
PREC (INPUT, OPTIONAL) real(stnd)

On entry, the desired accurary of the result. If more than six significant digits are required, the default value of PREC (e.g. 0.5e-06_stnd) should be altered appropriately (e.g. decreased). PREC is a small strictly positive integer less than 0.5e-06_stnd.

The default value for PREC is 0.5e-06_stnd .

ACU (INPUT, OPTIONAL) real(stnd)

On entry, the desired accuracy of the result when computing the incomplete Gamma integral in the evaluation of the seven term Taylor series. If l decimal places of accuracy are required then ACU should be set to 10**(-(l+1)). ACU is a small strictly positive integer. ACU should not be set smaller than the machine precision since the stated accuracy cannot be attained. In that case the machine precision is used instead.

The default value for ACU is epsilon( ACU ).

MAXITER (INPUT, OPTIONAL) integer(i4b)

On entry, MAXITER controls the maximum number of iterations when evaluating the Pearson’s series or continued fraction expansion of the incomplete Gamma integral.

See the description of the PROBGAMMA2 function for more details on this argument.

The default value is 1000.

Further Details

This function is both more general (here the number of degrees of freedom, DF, is not necessarily an integer) and more accurate (here the quantile X0 may be calculated as exactly as the computer allows with the parameter PREC) than PINVQ function.

This function is adapted from:

  1. Best, D.J., and Roberts, D.E., 1975:
    Algorithm AS 91: The Percentage Points of the chi2 Distribution. Appl. Statist., Vol.24, No. 3, pp.385-388
  2. Shea, B.L., 1991:
    Algorithm AS R85: A remark on AS 91: The Percentage Points of the chi2 Distribution. Appl. Statist. Vol.40, No. 1, pp.233-235.

function pinvq2 ( p, df, prec, acu, maxiter )

Purpose

Evaluates the inverse of the chi-squared distribution function:

X0 = PINVQ2( P(i), DF )

, if P(i) = probability( U < X0 ) for U = Chi-squared(DF) and i=1 to size(P).

PINVQ2 returns the quantiles X0(i) of the chi-squared distribution with DF degrees of freedom corresponding to a given lower tail area of P(i), for i=1 to size(P). In other words, PINVQ2 outputs chi-squared values, X0(:), such that random variables, distributed as chi-squared with DF degrees of freedom, will be less than or equal to X0(:) with associated probabilities P(:).

Arguments

P (INPUT) real(stnd), dimension(:)
On entry, the probabilities. P(i) must verify 0. < P(i) < 1, for i=1 to size(P).
DF (INPUT) real(stnd)
On entry, degrees of freedom of the chi-squared distribution. DF must be greater or equal to 0.5.
PREC (INPUT, OPTIONAL) real(stnd)

On entry, the desired accurary of the result. If more than six significant digits are required, the default value of PREC (e.g. 0.5e-06_stnd) should be altered appropriately(e.g. decreased). PREC is a small strictly positive integer less than 0.5e-06_stnd.

The default value for PREC is 0.5e-06_stnd .

ACU (INPUT, OPTIONAL) real(stnd)

On entry, the desired accuracy of the result when computing the incomplete Gamma integral in the evaluation of the seven term Taylor series. If l decimal places of accuracy are required then ACU should be set to 10**(-(l+1)). ACU is a small strictly positive integer. ACU should not be set smaller than the machine precision since the stated accuracy cannot be attained. In that case the machine precision is used instead.

The default value for ACU is epsilon( ACU ).

MAXITER (INPUT, OPTIONAL) integer(i4b)

On entry, MAXITER controls the maximum number of iterations when evaluating the Pearson’s series or continued fraction expansion of the incomplete Gamma integral.

See the description of the PROBGAMMA2 function for more details on this argument.

The default value is 1000.

Further Details

This function is both more general (here the number of degrees of freedom, DF, is not necessarily an integer) and more accurate (here the quantiles X0(:) may be calculated as exactly as the computer allows with the parameter PREC) than PINVQ function.

This function is parallelized if OPENMP is used.

This function is adapted from:

  1. Best, D.J., and Roberts, D.E., 1975:
    Algorithm AS 91: The Percentage Points of the chi2 Distribution. Appl. Statist., Vol.24, No. 3, pp.385-388
  2. Shea, B.L., 1991:
    Algorithm AS R85: A remark on AS 91: The Percentage Points of the chi2 Distribution. Appl. Statist. Vol.40, No. 1, pp.233-235.

function pinvq2 ( p, df, prec, acu, maxiter )

Purpose

Evaluates the inverse of the chi-squared distribution function:

X0 = PINVQ2( P(i), DF(i) )

, if P(i) = probability( U < X0 ) for U = Chi-squared(DF(i)) and i=1 to size(P).

PINVQ2 returns the quantiles X0(i) of the chi-squared distribution with DF(i) degrees of freedom corresponding to a given lower tail area of P(i), for i=1 to size(P). In other words, PINVQ2 outputs chi-squared values, X0(:), such that random variables, distributed as chi-squared with DF(:) degrees of freedom, will be less than or equal to X0(:) with associated probabilities P(:).

Arguments

P (INPUT) real(stnd), dimension(:)
On entry, the probabilities. P(i) must verify 0. < P(i) < 1, for i=1 to size(P) .
DF (INPUT) real(stnd), dimension(:)

On entry, degrees of freedom of the chi-squared distribution. Any value in the array DF must be greater or equal to 0.5.

The size of DF must verify size(DF) = size(P).

PREC (INPUT, OPTIONAL) real(stnd)

On entry, the desired accurary of the result. If more than six significant digits are required, the default value of PREC (e.g. 0.5e-06_stnd) should be altered appropriately(e.g. decreased). PREC is a small strictly positive integer less than 0.5e-06_stnd.

The default value for PREC is 0.5e-06_stnd .

ACU (INPUT, OPTIONAL) real(stnd)

On entry, the desired accuracy of the result when computing the incomplete Gamma integral in the evaluation of the seven term Taylor series. If l decimal places of accuracy are required then ACU should be set to 10**(-(l+1)). ACU is a small strictly positive integer. ACU should not be set smaller than the machine precision since the stated accuracy cannot be attained. In that case the machine precision is used instead.

The default value for ACU is epsilon( ACU ).

MAXITER (INPUT, OPTIONAL) integer(i4b)

On entry, MAXITER controls the maximum number of iterations when evaluating the Pearson’s series or continued fraction expansion of the incomplete Gamma integral.

See the description of the PROBGAMMA2 function for more details on this argument.

The default value is 1000.

Further Details

This function is both more general (here the numbers of degrees of freedom, DF(:), are not necessarily integers) and more accurate (here the quantiles X0(:) may be calculated as exactly as the computer allows with the parameter PREC) than PINVQ function.

This function is parallelized if OPENMP is used.

This function is adapted from:

  1. Best, D.J., and Roberts, D.E., 1975:
    Algorithm AS 91: The Percentage Points of the chi2 Distribution. Appl. Statist., Vol.24, No. 3, pp.385-388
  2. Shea, B.L., 1991:
    Algorithm AS R85: A remark on AS 91: The Percentage Points of the chi2 Distribution. Appl. Statist. Vol.40, No. 1, pp.233-235.

function pinvq2 ( p, df, prec, acu, maxiter )

Purpose

Evaluates the inverse of the chi-squared distribution function:

X0 = PINVQ2( P(i,j), DF )

, if P(i,j) = probability( U < X0 ) for U = Chi-squared(DF), i=1 to size(P,1) and j=1 to size(P,2).

PINVQ2 returns the quantiles X0(i,j) of the chi-squared distribution with DF degrees of freedom corresponding to a given lower tail area of P(i,j), for i=1 to size(P) and j=1 to size(P,2). In other words, PINVQ2 outputs chi-squared values, X0(:,:), such that random variables, distributed as chi-squared with DF degrees of freedom, will be less than or equal to X0(:,:) with associated probabilities P(:,:).

Arguments

P (INPUT) real(stnd), dimension(:,:)
On entry, the probabilities. P(i,j) must verify 0. < P(i,j) < 1, for i=1 to size(P) and j=1 to size(P,2).
DF (INPUT) real(stnd)
On entry, degrees of freedom of the chi-squared distribution. DF must be greater or equal to 0.5.
PREC (INPUT, OPTIONAL) real(stnd)

On entry, the desired accurary of the result. If more than six significant digits are required, the default value of PREC (e.g. 0.5e-06_stnd) should be altered appropriately(e.g. decreased). PREC is a small strictly positive integer less than 0.5e-06_stnd.

The default value for PREC is 0.5e-06_stnd .

ACU (INPUT, OPTIONAL) real(stnd)

On entry, the desired accuracy of the result when computing the incomplete Gamma integral in the evaluation of the seven term Taylor series. If l decimal places of accuracy are required then ACU should be set to 10**(-(l+1)). ACU is a small strictly positive integer. ACU should not be set smaller than the machine precision since the stated accuracy cannot be attained. In that case the machine precision is used instead.

The default value for ACU is epsilon( ACU ).

MAXITER (INPUT, OPTIONAL) integer(i4b)

On entry, MAXITER controls the maximum number of iterations when evaluating the Pearson’s series or continued fraction expansion of the incomplete Gamma integral.

See the description of the PROBGAMMA2 function for more details on this argument.

The default value is 1000.

Further Details

This function is both more general (here the number of degrees of freedom, DF, is not necessarily an integer) and more accurate (here the quantiles X0(:,:) may be calculated as exactly as the computer allows with the parameter PREC) than PINVQ function.

This function is parallelized if OPENMP is used.

This function is adapted from:

  1. Best, D.J., and Roberts, D.E., 1975:
    Algorithm AS 91: The Percentage Points of the chi2 Distribution. Appl. Statist., Vol.24, No. 3, pp.385-388
  2. Shea, B.L., 1991:
    Algorithm AS R85: A remark on AS 91: The Percentage Points of the chi2 Distribution. Appl. Statist. Vol.40, No. 1, pp.233-235.

function pinvq2 ( p, df, prec, acu, maxiter )

Purpose

Evaluates the inverse of the chi-squared distribution function:

X0 = PINVQ2( P(i,j), DF(i,j) )

, if P(i,j) = probability( U < X0 ) for U = Chi-squared(DF(i,j)), i=1 to size(P,1) and j=1 to size(P,2).

PINVQ2 returns the quantiles X0(i,j) of the chi-squared distribution with DF(i,j) degrees of freedom corresponding to a given lower tail area of P(i,j), for i=1 to size(P) and j=1 to size(P,2). In other words, PINVQ2 outputs chi-squared values, X0(:,:), such that random variables, distributed as chi-squared with corresponding DF(:,:) degrees of freedom, will be less than or equal to X0(:,:) with associated probabilities P(:,:).

Arguments

P (INPUT) real(stnd), dimension(:,:)
On entry, the probabilities. P(i,j) must verify 0. < P(i,j) < 1, for i=1 to size(P) and j=1 to size(P,2).
DF (INPUT) real(stnd), dimension(:,:)

On entry, degrees of freedom of the chi-squared distribution. Any value in the array DF must be greater or equal to 0.5.

The shape of DF must verify:

  • size(DF,1) = size(P,1)
  • size(DF,2) = size(P,2).
PREC (INPUT, OPTIONAL) real(stnd)

On entry, the desired accurary of the result. If more than six significant digits are required, the default value of PREC (e.g. 0.5e-06_stnd) should be altered appropriately(e.g. decreased). PREC is a small strictly positive integer less than 0.5e-06_stnd.

The default value for PREC is 0.5e-06_stnd .

ACU (INPUT, OPTIONAL) real(stnd)

On entry, the desired accuracy of the result when computing the incomplete Gamma integral in the evaluation of the seven term Taylor series. If l decimal places of accuracy are required then ACU should be set to 10**(-(l+1)). ACU is a small strictly positive integer. ACU should not be set smaller than the machine precision since the stated accuracy cannot be attained. In that case the machine precision is used instead.

The default value for ACU is epsilon( ACU ).

MAXITER (INPUT, OPTIONAL) integer(i4b)

On entry, MAXITER controls the maximum number of iterations when evaluating the Pearson’s series or continued fraction expansion of the incomplete Gamma integral.

See the description of the PROBGAMMA2 function for more details on this argument.

The default value is 1000.

Further Details

This function is both more general (here the numbers of degrees of freedom, DF(:,:), are not necessarily integers) and more accurate (here the quantiles X0(:,:) may be calculated as exactly as the computer allows with the parameter PREC) than PINVQ function.

This function is parallelized if OPENMP is used.

This function is adapted from:

  1. Best, D.J., and Roberts, D.E., 1975:
    Algorithm AS 91: The Percentage Points of the chi2 Distribution. Appl. Statist., Vol.24, No. 3, pp.385-388
  2. Shea, B.L., 1991:
    Algorithm AS R85: A remark on AS 91: The Percentage Points of the chi2 Distribution. Appl. Statist. Vol.40, No. 1, pp.233-235.

function probf ( f, ndf1, ndf2, upper )

Purpose

Evaluates the F distribution function with degrees of freedom NDF1 and NDF2 from F to infinity if UPPER is true or from zero to F if UPPER is false. In other words, if:

  • UPPER = true : PROBF = prob( U > F ) ,
  • UPPER = false : PROBF = prob( U < F ) ,

, for U = Fisher(NDF1,NDF2).

Arguments

F (INPUT) real(stnd)
On entry, upper or lower limit of integration. F must be greater or equal to zero.
NDF1 (INPUT) integer(i4b)
On entry, first degree of freedom of the F distribution (numerator). NDF1 must be greater or equal to 1.
NDF2 (INPUT) integer(i4b)
On entry, second degree of freedom of the F distribution (denominator). NDF2 must be greater or equal to 1.
UPPER (INPUT) logical(lgl)

On entry, if:

  • UPPER = true : probability to the right of F is calculated.
  • UPPER = false : probability to the left of F is calculated.

Further Details

This function accepts only integer values of degree of freedom and uses a normal approximation. This normal approximation is not accurate for small values of degrees of freedom.

This function is adapted from:

  1. Peizer, D.B., and Pratt, J.W., 1968:
    A normal approximation for Binomial, F, Beta, and …, (formula 2.24a). J.A.S.A., Vol. 63, 1457-1483

function probf ( f, ndf1, ndf2, upper )

Purpose

Evaluates the F distribution function with degrees of freedom NDF1 and NDF2 from F(i) to infinity if UPPER is true or from zero to F(i) if UPPER is false, for i=1 to size(F). In other words, if:

  • UPPER = true : PROBF( i ) = prob( U > F(i) ) ,
  • UPPER = false : PROBF( i ) = prob( U < F(i) ) ,

, for U = Fisher(NDF1,NDF2) and i=1 to size(F).

Arguments

F (INPUT) real(stnd), dimension(:)
On entry, upper or lower limits of integration. F(i) must be greater or equal to zero for i=1 to size(F).
NDF1 (INPUT) integer(i4b)
On entry, first degree of freedom of the F distribution (numerator). NDF1 must be greater or equal to 1.
NDF2 (INPUT) integer(i4b)
On entry, second degree of freedom of the F distribution (denominator). NDF2 must be greater or equal to 1.
UPPER (INPUT) logical(lgl)

On entry, if:

  • UPPER = true : probability to the right of F is calculated.
  • UPPER = false : probability to the left of F is calculated.

Further Details

This function accepts only integer values of degree of freedom and uses a normal approximation. This normal approximation is not accurate for small values of degrees of freedom.

This function is parallelized if OPENMP is used.

This function is adapted from:

  1. Peizer, D.B., and Pratt, J.W., 1968:
    A normal approximation for Binomial, F, Beta, and …, (formula 2.24a). J.A.S.A., Vol. 63, 1457-1483

function probf ( f, ndf1, ndf2, upper )

Purpose

Evaluates the F distribution function with degrees of freedom NDF1 and NDF2 from F(i) to infinity if UPPER is true or from zero to F(i) if UPPER is false, for i=1 to size(F). In other words, if:

  • UPPER = true : PROBF( i ) = prob( U > F(i) ) ,
  • UPPER = false : PROBF( i ) = prob( U < F(i) ) ,

, for U = Fisher(NDF1(i),NDF2(i)) and i=1 to size(F).

Arguments

F (INPUT) real(stnd), dimension(:)
On entry, upper or lower limits of integration. F(i) must be greater or equal to zero for i=1 to size(F).
NDF1 (INPUT) integer(i4b), dimension(:)

On entry, first degree of freedom of the F distribution (numerator). Any value in the array NDF1 must be greater or equal to 1.

The size of NDF1 must be size(NDF1) = size(F) .

NDF2 (INPUT) integer(i4b), dimension(:)

On entry, second degree of freedom of the F distribution (denominator). Any value in the array NDF2 must be greater or equal to 1.

The size of NDF2 must be size(NDF2) = size(F) .

UPPER (INPUT) logical(lgl)

On entry, if:

  • UPPER = true : probability to the right of F is calculated.
  • UPPER = false : probability to the left of F is calculated.

Further Details

This function accepts only integer values of degree of freedom and uses a normal approximation. This normal approximation is not accurate for small values of degrees of freedom.

This function is parallelized if OPENMP is used.

This function is adapted from:

  1. Peizer, D.B., and Pratt, J.W., 1968:
    A normal approximation for Binomial, F, Beta, and …, (formula 2.24a). J.A.S.A., Vol. 63, 1457-1483

function probf ( f, ndf1, ndf2, upper )

Purpose

Evaluates the F distribution function with degrees of freedom NDF1 and NDF2 from F(i,j) to infinity if UPPER is true or from zero to F(i,j) if UPPER is false, for i=1 to size(F,1) and j=1 to size(F,2). In other words, if:

  • UPPER = true : PROBF( i,j ) = prob( U > F(i,j) ) ,
  • UPPER = false : PROBF( i,j ) = prob( U < F(i,j) ) ,

, for U = Fisher(NDF1,NDF2) and i=1 to size(F,1) and j=1 to size(F,2).

Arguments

F (INPUT) real(stnd), dimension(:,:)
On entry, upper or lower limits of integration. F(i,j) must be greater or equal to zero for i=1 to size(F,1) and j=1 to size(F,2).
NDF1 (INPUT) integer(i4b)
On entry, first degree of freedom of the F distribution (numerator). NDF1 must be greater or equal to 1.
NDF2 (INPUT) integer(i4b)
On entry, second degree of freedom of the F distribution (denominator). NDF2 must be greater or equal to 1.
UPPER (INPUT) logical(lgl)

On entry, if:

  • UPPER = true : probability to the right of F is calculated.
  • UPPER = false : probability to the left of F is calculated.

Further Details

This function accepts only integer values of degree of freedom and uses a normal approximation. This normal approximation is not accurate for small values of degrees of freedom.

This function is parallelized if OPENMP is used.

This function is adapted from:

  1. Peizer, D.B., and Pratt, J.W., 1968:
    A normal approximation for Binomial, F, Beta, and …, (formula 2.24a). J.A.S.A., Vol. 63, 1457-1483

function probf ( f, ndf1, ndf2, upper )

Purpose

Evaluates the F distribution function with degrees of freedom NDF1 and NDF2 from F(i,j) to infinity if UPPER is true or from zero to F(i,j) if UPPER is false, for i=1 to size(F,1) and j=1 to size(F,2). In other words, if:

  • UPPER = true : PROBF( i,j ) = prob( U > F(i,j) ) ,
  • UPPER = false : PROBF( i,j ) = prob( U < F(i,j) ) ,

, for U = Fisher(NDF1( i,j ),NDF2( i,j )) and i=1 to size(F,1) and j=1 to size(F,2).

Arguments

F (INPUT) real(stnd), dimension(:,:)
On entry, upper or lower limits of integration. F(i,j) must be greater or equal to zero for i=1 to size(F,1) and j=1 to size(F,2).
NDF1 (INPUT) integer(i4b), dimension(:,:)

On entry, first degree of freedom of the F distribution (numerator). Any value in the array NDF1 must be greater or equal to 1.

The shape of NDF1 must verify:

  • size(NDF1,1) = size(F,1)
  • size(NDF1,2) = size(F,2) .
NDF2 (INPUT) integer(i4b), dimension(:,:)

On entry, second degree of freedom of the F distribution (denominator). Any value in the array NDF2 must be greater or equal to 1.

The shape of NDF2 must verify:

  • size(NDF2,1) = size(F,1)
  • size(NDF2,2) = size(F,2) .
UPPER (INPUT) logical(lgl)

On entry, if:

  • UPPER = true : probability to the right of F is calculated.
  • UPPER = false : probability to the left of F is calculated.

Further Details

This function accepts only integer values of degree of freedom and uses a normal approximation. This normal approximation is not accurate for small values of degrees of freedom.

This function is parallelized if OPENMP is used.

This function is adapted from:

  1. Peizer, D.B., and Pratt, J.W., 1968:
    A normal approximation for Binomial, F, Beta, and …, (formula 2.24a). J.A.S.A., Vol. 63, 1457-1483

function probf2 ( f, df1, df2, upper, beta, acu, maxiter, failure )

Purpose

Evaluates the F distribution function with degrees of freedom DF1 and DF2 (integer or fractional degrees of freedom) from F to infinity if UPPER is true or from zero to F if UPPER is false. In other words, if:

  • UPPER = true : PROBF2 = prob( U > F ) ,
  • UPPER = false : PROBF2 = prob( U <= F ) ,

, for U = Fisher(DF1,DF2).

For given arguments F (0<=F), DF1 (DF1>0), DF2 (DF2>0), PROBF2 returns the probability that a random variable from an F distribution having DF1 and DF2 degrees of freedom will be less than or equal to F (if UPPER is false) or greater than F (if UPPER is true).

Arguments

F (INPUT) real(stnd)
On entry, upper or lower limit of integration. F must be greater or equal to zero.
DF1 (INPUT) real(stnd)
On entry, first degree of freedom of the F distribution (numerator). DF1 must be greater than zero.
DF2 (INPUT) real(stnd)
On entry, second degree of freedom of the F distribution (denominator). DF2 must be greater than zero.
UPPER (INPUT) logical(lgl)

On entry, if:

  • UPPER = true : probability to the right of F is calculated.
  • UPPER = false : probability to the left of F is calculated.
BETA (INPUT, OPTIONAL) real(stnd)

On entry, the logarithm of the complete beta function BETA(0.5 * DF1,0.5 * DF2).

If BETA is not given, the logarithm of the beta function is computed with the help of function LNGAMMA.

ACU (INPUT, OPTIONAL) real(stnd)

On entry, the desired accuracy of the result, when computing the incomplete beta function. The “integrating” process for evaluating the incomplete beta function is terminated when the relative contribution to the integral is not greater than the value of ACU. ACU is a small strictly positive integer.

The default value for ACU is epsilon( ACU ).

MAXITER (INPUT, OPTIONAL) integer(i4b)

On entry, MAXITER controls the maximum number of iterations when evaluating the power series representation of the incomplete Beta integral.

The default value is 2000.

FAILURE (INPUT, OPTIONAL) logical(lgl)
On entry, if FAILURE is set to true, the values for which the “integrating” process of the incomplete Beta integral did not converge to the desired accuracy are set to -1. The algorithm fails to converge if the number of iterations exceeds MAXITER.

Further Details

This function invoked the BETA distribution function (e.g. function PROBBETA) for computing the probability associated with F and is much more accurate than PROBF function.

function probf2 ( f, df1, df2, upper, beta, acu, maxiter, failure )

Purpose

Evaluates the F distribution function with degrees of freedom DF1 and DF2 (integer or fractional degrees of freedom) from F(i) to infinity if UPPER is true or from zero to F(i) if UPPER is false, for i=1 to size(F). In other words, if:

  • UPPER = true : PROBF( i ) = prob( U > F(i) ) ,
  • UPPER = false : PROBF( i ) = prob( U < F(i) ) ,

, for U = Fisher(DF1,DF2) and i=1 to size(F).

For given arguments F(:) (0<=F(:)), DF1 (DF1>0), DF2 (DF2>0), PROBF2 returns the probability that a random variable (vector) from an F distribution having DF1 and DF2 degrees of freedom will be less than or equal to F (if UPPER is false) or greater than F (if UPPER is true).

Arguments

F (INPUT) real(stnd), dimension(:)
On entry, upper or lower limit of integration. F(i) must be greater or equal to zero for i=1 to size(F).
DF1 (INPUT) real(stnd)
On entry, first degree of freedom of the F distribution (numerator). DF1 must be greater than zero.
DF2 (INPUT) real(stnd)
On entry, second degree of freedom of the F distribution (denominator). DF2 must be greater than zero.
UPPER (INPUT) logical(lgl)

On entry, if:

  • UPPER = true : probability to the right of F is calculated.
  • UPPER = false : probability to the left of F is calculated.
BETA (INPUT, OPTIONAL) real(stnd)

On entry, the logarithm of the complete beta function BETA(0.5 * DF1,0.5 * DF2).

If BETA is not given, the logarithm of the beta function is computed with the help of function LNGAMMA.

ACU (INPUT, OPTIONAL) real(stnd)

On entry, the desired accuracy of the result, when computing the incomplete beta function. The “integrating” process for evaluating the incomplete beta function is terminated when the relative contribution to the integral is not greater than the value of ACU. ACU is a small strictly positive integer.

The default value for ACU is epsilon( ACU ).

MAXITER (INPUT, OPTIONAL) integer(i4b)

On entry, MAXITER controls the maximum number of iterations when evaluating the power series representation of the incomplete Beta integral.

The default value is 2000.

FAILURE (INPUT, OPTIONAL) logical(lgl)
On entry, if FAILURE is set to true, the values for which the “integrating” process of the incomplete Beta integral did not converge to the desired accuracy are set to -1. The algorithm fails to converge if the number of iterations exceeds MAXITER.

Further Details

This function invoked the BETA distribution function (e.g. function PROBBETA) for computing the probability associated with F and is much more accurate than PROBF function.

function probf2 ( f, df1, df2, upper, acu, maxiter, failure )

Purpose

Evaluates the F distribution function with degrees of freedom DF1(i) and DF2(i) (integer or fractional degrees of freedom) from F(i) to infinity if UPPER is true or from zero to F(i) if UPPER is false, for i=1 to size(F). In other words, if:

  • UPPER = true : PROBF( i ) = prob( U > F(i) ) ,
  • UPPER = false : PROBF( i ) = prob( U < F(i) ) ,

, for U = Fisher(DF1(i),DF2(i)) and i=1 to size(F).

For given arguments F(:) (0<=F(:)), DF1(:) (DF1(:)>0), DF2(:) (DF2(:)>0), PROBF2 returns the probability that a random variable (vector) from an F distribution having DF1(:) and DF2(:) degrees of freedom will be less than or equal to F(:) (if UPPER is false) or greater than F(:) (if UPPER is true).

Arguments

F (INPUT) real(stnd), dimension(:)
On entry, upper or lower limit of integration. F(i) must be greater or equal to zero for i=1 to size(F).
DF1 (INPUT) real(stnd), dimension(:)

On entry, first degree of freedom of the F distribution (numerator). Any value in the array DF1(:) must be greater than zero.

The size of DF1 must verify size(DF1) = size(F) .

DF2 (INPUT) real(stnd), dimension(:)

On entry, second degree of freedom of the F distribution (denominator). Any value in the array DF2(:) must be greater than zero.

The size of DF2 must verify size(DF2) = size(F) .

UPPER (INPUT) logical(lgl)

On entry, if:

  • UPPER = true : probability to the right of F is calculated.
  • UPPER = false : probability to the left of F is calculated.
BETA (INPUT, OPTIONAL) real(stnd)

On entry, the logarithm of the complete beta function BETA(0.5 * DF1,0.5 * DF2).

If BETA is not given, the logarithm of the beta function is computed with the help of function LNGAMMA.

ACU (INPUT, OPTIONAL) real(stnd)

On entry, the desired accuracy of the result, when computing the incomplete beta function. The “integrating” process for evaluating the incomplete beta function is terminated when the relative contribution to the integral is not greater than the value of ACU. ACU is a small strictly positive integer.

The default value for ACU is epsilon( ACU ).

MAXITER (INPUT, OPTIONAL) integer(i4b)

On entry, MAXITER controls the maximum number of iterations when evaluating the power series representation of the incomplete Beta integral.

The default value is 2000.

FAILURE (INPUT, OPTIONAL) logical(lgl)
On entry, if FAILURE is set to true, the values for which the “integrating” process of the incomplete Beta integral did not converge to the desired accuracy are set to -1. The algorithm fails to converge if the number of iterations exceeds MAXITER.

Further Details

This function invoked the BETA distribution function (e.g. function PROBBETA) for computing the probability associated with F and is much more accurate than PROBF function.

The function is parallelized if OPENMP is used.

function probf2 ( f, df1, df2, upper, beta, acu, maxiter, failure )

Purpose

Evaluates the F distribution function with degrees of freedom DF1 and DF2 (integer or fractional degrees of freedom) from F(i,j) to infinity if UPPER is true or from zero to F(i,j) if UPPER is false, for i=1 to size(F,1) and j=1 to size(F,2):

if UPPER = true : PROBF( i, j ) = prob( U > F(i,j) ) , if UPPER = false : PROBF( i, j ) = prob( U < F(i,j) ) ,

, for U = Fisher(DF1,DF2) and i=1 to size(F,1) and j=1 to size(F,2).

For given arguments F(:,:) (0<=F(:,:)), DF1 (DF1>0), DF2 (DF2>0), PROBF2 returns the probability that a random variable (matrix) from an F distribution having DF1 and DF2 degrees of freedom will be less than or equal to F (if UPPER is false) or greater than F (if UPPER is true).

Arguments

F (INPUT) real(stnd), dimension(:,:)
On entry, upper or lower limit of integration. F(i,j) must be greater or equal to zero for i=1 to size(F,1) and j=1 to size(F,2).
DF1 (INPUT) real(stnd)
On entry, first degree of freedom of the F distribution (numerator). DF1 must be greater than zero.
DF2 (INPUT) real(stnd)
On entry, second degree of freedom of the F distribution (denominator). DF2 must be greater than zero.
UPPER (INPUT) logical(lgl)

On entry, if:

  • UPPER = true : probability to the right of F is calculated.
  • UPPER = false : probability to the left of F is calculated.
BETA (INPUT, OPTIONAL) real(stnd)

On entry, the logarithm of the complete beta function BETA(0.5 * DF1,0.5 * DF2).

If BETA is not given, the logarithm of the beta function is computed with the help of function LNGAMMA.

ACU (INPUT, OPTIONAL) real(stnd)

On entry, the desired accuracy of the result, when computing the incomplete beta function. The “integrating” process for evaluating the incomplete beta function is terminated when the relative contribution to the integral is not greater than the value of ACU. ACU is a small strictly positive integer.

The default value for ACU is epsilon( ACU ).

MAXITER (INPUT, OPTIONAL) integer(i4b)

On entry, MAXITER controls the maximum number of iterations when evaluating the power series representation of the incomplete Beta integral.

The default value is 2000.

FAILURE (INPUT, OPTIONAL) logical(lgl)
On entry, if FAILURE is set to true, the values for which the “integrating” process of the incomplete Beta integral did not converge to the desired accuracy are set to -1. The algorithm fails to converge if the number of iterations exceeds MAXITER.

Further Details

This function invoked the BETA distribution function (e.g. function PROBBETA) for computing the probability associated with F and is much more accurate than PROBF function.

function probf2 ( f, df1, df2, upper, acu, maxiter, failure )

Purpose

Evaluates the F distribution function with degrees of freedom DF1(i,j) and DF2(i,j) (integer or fractional degrees of freedom) from F(i,j) to infinity if UPPER is true or from zero to F(i,j) if UPPER is false, for i=1 to size(F,1) and j=1 to size(F,2). In other words, if:

  • UPPER = true : PROBF( i, j ) = prob( U > F(i,j) ) ,
  • UPPER = false : PROBF( i, j ) = prob( U < F(i,j) ) ,

, for U = Fisher(DF1(i,j),DF2(i,j)) and i=1 to size(F,1) and j=1 to size(F,2).

For given arguments F(:,:) (0<=F(:,:)), DF1(:,:) (DF1(:,:)>0), DF2(:,:) (DF2(:,:)>0), PROBF2 returns the probability that a random variable (matrix) from an F distribution having DF1(:,:) and DF2(:,:) degrees of freedom will be less than or equal to F(:,:) (if UPPER is false) or greater than F (if UPPER is true).

Arguments

F (INPUT) real(stnd), dimension(:,:)
On entry, upper or lower limit of integration. F(i,j) must be greater or equal to zero for i=1 to size(F,1) and j=1 to size(F,2).
DF1 (INPUT) real(stnd), dimension(:,:)

On entry, first degree of freedom of the F distribution (numerator). Any value in the array DF1(:,:) must be greater than zero.

The shape of DF1 must verify:

  • size(DF1,1) = size(F,1)
  • size(DF1,2) = size(F,2) .
DF2 (INPUT) real(stnd), dimension(:,:)

On entry, second degree of freedom of the F distribution (denominator). Any value in the array DF2(:,:) must be greater than zero.

The shape of DF2 must verify:

  • size(DF2,1) = size(F,1)
  • size(DF2,2) = size(F,2) .
UPPER (INPUT) logical(lgl)

On entry, if:

  • UPPER = true : probability to the right of F is calculated.
  • UPPER = false : probability to the left of F is calculated.
ACU (INPUT, OPTIONAL) real(stnd)

On entry, the desired accuracy of the result, when computing the incomplete beta function. The “integrating” process for evaluating the incomplete beta function is terminated when the relative contribution to the integral is not greater than the value of ACU. ACU is a small strictly positive integer.

The default value for ACU is epsilon( ACU ).

MAXITER (INPUT, OPTIONAL) integer(i4b)

On entry, MAXITER controls the maximum number of iterations when evaluating the power series representation of the incomplete Beta integral.

The default value is 2000.

FAILURE (INPUT, OPTIONAL) logical(lgl)
On entry, if FAILURE is set to true, the values for which the “integrating” process of the incomplete Beta integral did not converge to the desired accuracy are set to -1. The algorithm fails to converge if the number of iterations exceeds MAXITER.

Further Details

This function invoked the BETA distribution function (e.g. function PROBBETA) for computing the probability associated with F and is much more accurate than PROBF function.

The function is parallelized if OPENMP is used.

function pinvf2 ( p, df1, df2, beta, acu, maxiter )

Purpose

Evaluates the inverse F probability distribution function with degrees of freedom DF1 and DF2 (integer or fractional degrees of freedom).

For given arguments P (0<=P<=1), DF1 (DF1>0.2), DF2 (DF2>0.2), PINVF2 returns the value F such that the probability that a random variable distributed as F(DF1,DF2) is less than or equal to F is P.

Arguments

P (INPUT) real(stnd)
On entry, input probability in the inclusive range (0,1).
DF1 (INPUT) real(stnd)
On entry, first degree of freedom of the F distribution (numerator). DF1 must be greater than 0.2.
DF2 (INPUT) real(stnd)
On entry, second degree of freedom of the F distribution (denominator). DF2 must be greater than 0.2.
BETA (INPUT, OPTIONAL) real(stnd)

On entry, the logarithm of the complete beta function BETA(0.5 * DF1,0.5 * DF2).

If BETA is not given, the logarithm of the beta function is computed with the help of function LNGAMMA.

ACU (INPUT, OPTIONAL) real(stnd)

On entry, the desired accuracy of the result, when computing the incomplete beta function. The “integrating” process for evaluating the incomplete beta function is terminated when the relative contribution to the integral is not greater than the value of ACU. ACU is a small strictly positive integer.

The default value for ACU is epsilon( ACU ).

MAXITER (INPUT, OPTIONAL) integer(i4b)

On entry, MAXITER controls the maximum number of iterations when evaluating the power series representation of the incomplete Beta integral.

See the description of the PROBBETA function for more details on this argument.

The default value is 2000.

Further Details

This function invoked the inverse BETA distribution function (e.g. function PINVBETA) for computing the value F associated with the probability P.

This function is not very accurate for small values of DF1 and/or DF2 (e.g. less than 1).

function probbinom ( prob, n, k, upper, beta, acu, maxiter, failure )

Purpose

Evaluates the cumulative binomial probability distribution function for a positive real argument PROB between 0 and 1, a strictly positive integer N and a positive integer K less than or equal to N.

PROBBINOM computes the probability that an event occurring with probability PROB per trial, will occur K or more times in N independent trials if UPPER is true, or will occur K or less times in N independent trials if UPPER is false.

This probability is estimated with the help of the incomplete beta fonction, as computed by PROBBETA function, and the optional arguments BETA, ACU, MAXITER and FAILURE are passed directly to the PROBBETA function if these arguments are present.

Arguments

PROB (INPUT) real(stnd)
On entry, a positive real argument PROB which is the probability of success on each trial for the given binomial probability distribution. PROB must be greater or equal to zero and less than or equal to 1.
N (INPUT) integer(i4b)
On entry, a strictly positive integer argument which is the total number of Bernoulli trials for the given binomial probability distribution. N must be greater than zero.
K (INPUT) integer(i4b)
On entry, a positive integer argument which is the minimal (if UPPER is true) or maximum (if UPPER is false) number of success for the N Bernoulli trials, for which we want to compute the cumulative probability following the binomial probability distribution specified by PROB. K must be greater or equal to zero and less or equal to N.
UPPER (INPUT) logical(lgl)

On entry, if:

  • UPPER = true : the probability that the number of success is greater than or equal to K in N trials is computed.
  • UPPER = false : the probability that the number of success is less than or equal to K in N trials is computed.
BETA (INPUT, OPTIONAL) real(stnd)

On entry, the logarithm of the complete beta function BETA(K,N-K+1) if UPPER is true, or BETA(N-K,K+1) if UPPER is false.

If BETA is not given, the logarithm of the beta function is computed with the help of function LNGAMMA.

ACU (INPUT, OPTIONAL) real(stnd)

On entry, the desired accuracy of the result when evaluating the incomplete beta fonction. The “integrating” process is terminated when both the absolute and relative contributions to the integral is not greater than the value of ACU.

ACU is a small strictly positive integer. If the number of decimal digits’ accuracy required is r, ACU should be set to 10**-(r+1).

The default value for ACU is epsilon( ACU ).

MAXITER (INPUT, OPTIONAL) integer(i4b)

On entry, MAXITER controls the maximum number of iterations when evaluating the power series representation of the incomplete Beta integral.

The default value is 2000.

FAILURE (INPUT, OPTIONAL) logical(lgl)

On entry, if FAILURE is set to true, the values for which the “integrating” process did not converge to the desired accuracy are set to -1. The algorithm fails to converge if the number of iterations exceeds MAXITER.

The default value is false.

Further Details

The cumulative binomial probability is computed with the help of the incomplete Beta function PROBBETA as follow:

If UPPER is true, the probability of an event will occur K or more times in N independant trials, if its probability per trial is PROB, is computed as PROBBETA( PROB, K, N-K+1).

If UPPER is false, the probability of an event will occur K or less times in N independant trials, if its probability per trial is PROB, is computed as PROBBETA( 1-PROB, N-K, K+1).

function rangen ( x, n )

Purpose

Evaluates the probability of the normal range given X, the upper limit of integration, and N, the sample size.

In other words, this function evaluates the probability that the standardized difference between the maximum and the minimum on a sample will be less than X for a normal (e.g. gaussian) sample of size N.

Arguments

X (INPUT) real(stnd)
On entry, upper limit of integration. X must be greater than zero.
N (INPUT) integer(i4b)
On entry, the sample size. N must be greater than 1.

Further Details

This function is adapted from:

  1. Barnard, J., 1978:
    Algorithm AS126: Probability Integral of the normal range, Applied Statistics, Vol. 27, No. 2, pp.197-198

function rangen ( x, n )

Purpose

Evaluates probabilities of the normal range given X(:), the upper limits of integration, and N, the sample size.

In other words, this function evaluates the probabilities that the standardized difference between the maximum and the minimum on a sample will be less than X(:) for a normal (e.g. gaussian) sample of size N.

Arguments

X (INPUT) real(stnd), dimension(:)
On entry, upper limits of integration. All elements in X(:) must be greater than zero.
N (INPUT) integer(i4b)
On entry, the sample size. N must be greater than 1.

Further Details

This function is adapted from:

  1. Barnard, J., 1978:
    Algorithm AS126: Probability Integral of the normal range, Applied Statistics, Vol. 27, No. 2, pp.197-198
Flag Counter