using System;
using System.Collections.Generic;
namespace QuantRiskLib
{
///Source: www.risk256.com
///
///References:
///Miller, Michael B. 2012. Mathematics and Statistics for Financial Risk Management. New York: John Wiley & Sons.
///Chapter 4
public partial class Distributions
{
#region Beta Distribution
///
/// Returns the CDF of a beta distribution with shape parameters a and b.
/// PDF is: f(x) = x^(a-1) * (1-x)^(b-1) / B(a,b)
///
/// Value of the random variable for which the CDF is beign evaluated. x is between 0 and 1.
/// First shape parameter.
/// Second shape parameter.
///
public static double BetaCumulativeDistributionFunction(double x, double a, double b)
{
if(x < 0 || x > 1)
throw new ArgumentException("x must be between 0 and 1");
return MMath.RegularizedIncompleteBetaFunction(x, a, b);
}
///
/// Returns the PDF of a beta distribution with shape parameters a and b.
/// f(x) = x^(a-1) * (1-x)^(b-1) / B(a,b)
///
/// Value of the random variable for which the PDF is beign evaluated. x is between 0 and 1.
/// Number of trials.
/// Number of times the event occurs in n trials.
///
public static double BetaProbabilityDensityFunction(double x, double a, double b)
{
if (x < 0 || x > 1)
throw new ArgumentException("x must be between 0 and 1");
double B = MMath.BetaFunction(a, b);
double d = Math.Pow(x, a - 1) * Math.Pow(1-x, b-1) / B;
return d;
}
#endregion
#region Binomial Distribution
///
/// Returns the probability of an event occuring in k or less times out of n trials, given
/// that the probability of the event occuring in any one trial is p.
///
/// Probability that the event occurs in any one trial.
/// Number of trials.
/// Number of times the event occurs in n trials.
///
public static double BinomialCumulativeDistributionFunction(double p, int n, int k)
{
if(k >= n)
return 1.0;
double sum = 0;
for (int i = 0; i <= k; i++)
sum += BinomialProbabilityDensityFunction(p, n, i);
return sum;
}
///
/// Returns the probability of an event occuring in k out of n trials, given
/// that the probability of the event occuring in any one trial is p.
///
/// Probability that the event occurs in any one trial.
/// Number of trials.
/// Number of times the event occurs in n trials.
///
public static double BinomialProbabilityDensityFunction(double p, int n, int k)
{
//# of combinations can exceed max double for n > 1020
if (n > 1020)
{
double logB = MMath.LogCombin(n, k) + k * Math.Log(p) + (n - k) * Math.Log(1 - p);
return Math.Exp(logB);
}
return MMath.Combinations(n, k) * Math.Pow(p, k) * Math.Pow(1 - p, n - k);
}
#endregion
#region Chi-Squared Distribution
///
/// Chi-squared probability density function.
///
/// The value at which the PDF is evaluated.
/// Degress of freedom, or number independent standard normal distributions.
///
public static double ChiSquaredProbabilityDensityFunction(double x, int k)
{
double g = MMath.GammaFunction(0.5 * k);
double a = 1.0 / (Math.Pow(2, 0.5 * k) * g);
return a * Math.Pow(x, 0.5 * k - 1.0) * Math.Exp(-0.5 * x);
}
#endregion
#region F Distribution
///
/// Returns the PDF of the F distribution.
///
/// Value at which the distribution is evaluated.
/// Degrees of freedom for numerator chi-sqared distribution. k1 > 0.
/// Degrees of freedom for denominator chi-sqared distribution. k2 > 0.
public static double FProbabilityDensityFunction(double x, int k1, int k2)
{
if (k1 <= 0 || k2 <= 0) throw new ArgumentException("k1 and k2 must be greater than 0.");
if (x == 0) return 0.0;
double a1 = Math.Pow(k1 * x, 0.5 * k1);
double a2 = Math.Pow(k2, 0.5 * k2);
double a3 = Math.Pow(k1 * x + k2, 0.5 * (k1 + k2));
double a = a1 * a2 / a3;
double b = MMath.BetaFunction(0.5 * k1, 0.5 * k2);
double c = x * b;
return a / c;
}
///
/// Returns the CDF of the F distribution.
///
/// Value at which the distribution is evaluated.
/// Degrees of freedom for numerator chi-sqared distribution.
/// Degrees of freedom for denominator chi-sqared distribution.
public static double FCumulativeDensityFunction(double x, int k1, int k2)
{
//Tested against Excel for
//k1 & k2 = 1, 2, 3, ... , 9, 10, 20, ... , 100
//x = 0.0, 0.1, 0.2, ... , 6.0
//maximum discrepancy was less than 0.16% for (x, k1, k2) = (6, 1, 9)
if (k1 <= 0 || k2 <= 0) throw new ArgumentException("k1 and k2 must be greater than 0.");
if (x < 0) throw new ArgumentException("x must be greater than 0.");
if (x == 0)
return 0.0;
double x2 = (k1 * x) / (k1 * x + k2);
double p = MMath.RegularizedIncompleteBetaFunction(x2, 0.5 * k1, 0.5 * k2);
return Math.Min(1, p);
}
///
/// Returns the inverse of the CDF of the F distribution.
/// For k = 3, and 5+ the solution is an approximation.
///
/// Cumulative probability of the distribution. p is between 0 and 1.
/// Degrees of freedom for numerator chi-sqared distribution.
/// Degrees of freedom for denominator chi-sqared distribution.
public static double FCumulativeDistributionFunctionInverse(double p, int k1, int k2)
{
if (k1 <= 4 || k2 <= 4) throw new ArgumentException("k1 and k2 must be greater than 4.");
//Tested against Excel for k1 & k2 = 5, 6, ... 19, 20, 30, 40, ... , 100 and p = 0.01, 0.02, ..., 0.99.
//Max absolute difference was 0.0313 for (p, k1, k2) = (0.99, 7, 5).
const int maxIterations = 16;
const double minAccuracy = 0.0001; //min absolute difference between p and p implied by return value.
if (p < 0 || p > 1.0) throw new ArgumentException("Probability must be between 0 and 1");
if (p == 0) return 0.0;
if (p == 1) return double.PositiveInfinity;
double mu = k2 / (k2 - 2.0);
double sigma = 2 * k2 * k2 * (k1 + k2 - 2.0) / (k1 * (k2 - 2.0) * (k2 - 2.0) * (k2 - 4.0));
double normAprox = Math.Max(0, NormalCumulativeDistributionFunctionInverse(p, mu, sigma));
double x;
if (k1 < 10 && k2 < 10)
if (p < 0.03)
x = 0.1592355;
else
x = 0.1819022 * normAprox + 6.0340494 * Math.Log(p) / Math.Sqrt(k1 + k2) - 10.4102348 / Math.Sqrt(k1 + k2) - 3.8062818 * Math.Log(p) / k2 - 2.1284338 * Math.Log(p) / k1 - 0.0446516 * k2 - 0.1014539 * k1 - 7.980049829 * p + 7.988783845 * p * p + 6.864968931;
else if (k1 < 10)
x = 0.8100214 * normAprox + 3.1485113 * Math.Log(p) / Math.Sqrt(k1 + k2) - 5.63992350 / Math.Sqrt(k1 + k2) - 6.8494933 * Math.Log(p) / k2 - 0.5725057 * Math.Log(p) / k1 - 0.0102242 * k2 - 0.0099565 * k1 - 2.820859600 * p + 3.036058600 * p * p + 2.103689000;
else if (k2 < 10)
x = 0.1944974 * normAprox + 5.5846823 * Math.Log(p) / Math.Sqrt(k1 + k2) - 2.6254457 / Math.Sqrt(k1 + k2) - 0.7057866 * Math.Log(p) / k2 - 6.49649350 * Math.Log(p) / k1 + 0.0048548 * k2 - 0.0088253 * k1 - 6.523519381 * p + 6.710305803 * p * p + 3.305158009;
else
x = 0.8717865 * normAprox + 0.9139534 * Math.Log(p) / Math.Sqrt(k1 + k2) - 0.0864391 / Math.Sqrt(k1 + k2) - 1.5438103 * Math.Log(p) / k2 - 0.14952740 * Math.Log(p) / k1 + 0.0004922 * k2 - 0.0003965 * k1 - 0.636234600 * p + 1.177650500 * p * p + 0.103934500;
x = Math.Max(0.00001, x);
double oldAbsAccuracy = double.MaxValue;
double step = 0;
for (int i = 0; i < maxIterations; i++)
{
double p0 = FCumulativeDensityFunction(x, k1, k2);
double absAccuracy = Math.Abs(p - p0);
if (absAccuracy < minAccuracy) return x;
if (absAccuracy > oldAbsAccuracy)
{
x = x - 0.5 * step;
p0 = FCumulativeDensityFunction(x, k1, k2);
absAccuracy = Math.Abs(p - p0);
}
oldAbsAccuracy = absAccuracy;
double oldX = x;
double slope = FProbabilityDensityFunction(x, k1, k2);
step = (p - p0) / slope;
if (step < -x) step = -0.5 * x;
if (step > 1.0) step = 1.0;
x = x + step;
if (x < 0) x = 0.5 * oldX;
}
throw new Exception("Solution did not converge");
}
//currently not used ...might be useful for future version of FCumulativeDistributionFunctionInverse()
private static double DerivativeOfFProbabilityDensityFunction(double x, double k1, double k2, double pdf)
{
//double pdf = FProbabilityDensityFunction(x, k1, k2);
double a = (k1 * k2 * (1.0 - x)) / (k1 * x + k2);
return pdf * (1.0 / x) * (0.5 * a - 1.0);
}
#endregion
#region Non-Central Chi-Squared Distribution
///
/// Non-central chi-squared cumulative distribution function.
/// Approximation based on Sankaran (1963).
///
/// The value at which the CDF is evaluated. Must be non-negative.
/// degrees of freedom; cannot be negative
/// non-centrality parameter; cannot be negative
///
public static double NonCentralChiSquaredCumulativeDistributionFunction(double x, double k, double lambda)
{
if (k < 0) throw new ArgumentException("degrees of freedom cannot be negative");
if (lambda < 0) throw new ArgumentException("non-centrality parameter cannot be negative");
if (x < 0) throw new ArgumentException("x cannot be negative");
if (x == 0) return 0.0;
if (x == double.PositiveInfinity) return 1.0;
double h = 1.0 - (2.0/3.0)*(k + lambda)*(k + 3.0*lambda)/Math.Pow(k + 2.0*lambda, 2);
double p = (k + 2.0 * lambda) / Math.Pow(k + lambda, 2);
double m = (h - 1.0)*(1.0 - 3.0*h);
double f = Math.Pow(x/(k + lambda), h);
double s = 1.0 + h*p*(h - 1.0 - 0.5*(2.0 - h)*m*p);
double d = h*(1.0 + 0.5*m*p)*Math.Sqrt(2.0*p);
double z = (f - s)/d;
return StandardNormalCumulativeDistributionFunction(z);
}
#endregion
#region Normal Distribution
///
/// Returns the PDF of the normal distribution.
///
/// Value at which the distribution is evaluated.
/// Mean of the distribution.
/// Standard deviation of the distribution.
public static double NormalProbabilityDensityFunction(double x, double mean, double standardDeviation)
{
const double sqrtTwoPiInv = 0.398942280401433;
double z = (x - mean) / standardDeviation;
return sqrtTwoPiInv * Math.Exp(-0.5 * z * z) / standardDeviation;
}
///
/// Returns the CDF of the normal distribution.
///
/// Value at which the distribution is evaluated.
/// Mean of the distribution.
/// Standard deviation of the distribution.
public static double NormalCumulativeDistributionFunction(double x, double mean, double standardDeviation)
{
double z = (x - mean) / standardDeviation;
return StandardNormalCumulativeDistributionFunction(z);
}
///
/// Returns the inverse of the CDF of the normal distribution.
///
/// Cumulative probability of the distribution. 0 <= p >= 1.
///Mean of the distribution.
/// Standard deviation of the distribution.
public static double NormalCumulativeDistributionFunctionInverse(double p, double mean, double standardDeviation)
{
double x = StandardNormalCumulativeDistributionFunctionInverse(p);
return standardDeviation * x + mean;
}
///
/// Returns the PDF of the standard normal distribution.
///
/// Value at which the distribution is evaluated.
public static double StandardNormalProbabilityDensityFunction(double x)
{
const double SqrtTwoPiInv = 0.398942280401433;
return SqrtTwoPiInv * Math.Exp(-0.5 * x * x);
}
///
/// Returns the CDF of the standard normal distribution.
///
/// Value at which the distribution is evaluated.
public static double StandardNormalCumulativeDistributionFunction(double x)
{
//Approimation based on Abramowitz & Stegun (1964)
if (x < 0)
return 1.0 - StandardNormalCumulativeDistributionFunction(-x);
const double b0 = 0.2316419;
const double b1 = 0.319381530;
const double b2 = -0.356563782;
const double b3 = 1.781477937;
const double b4 = -1.821255978;
const double b5 = 1.330274429;
double pdf = StandardNormalProbabilityDensityFunction(x);
double a = 1.0 / (1.0 + b0 * x);
return 1.0 - pdf * (b1 * a + b2 * Math.Pow(a, 2) + b3 * Math.Pow(a, 3) + b4 * Math.Pow(a, 4) + b5 * Math.Pow(a, 5));
}
///
/// Returns the inverse of the CDF of the standard normal distribution.
///
/// Cumulative probability of the distribution. p is between 0 and 1.
public static double StandardNormalCumulativeDistributionFunctionInverse(double p)
{
//based on pseudo-code from http://home.online.no/~pjacklam/notes/invnorm/
if (p < 0 || p > 1.0) throw new ArgumentException("Probability must be between 0 and 1");
if (p == 0) return double.NegativeInfinity;
if (p == 1) return double.PositiveInfinity;
if (p == 0.5) return 0.0;
const double a1 = -3.969683028665376e+01;
const double a2 = 2.209460984245205e+02;
const double a3 = -2.759285104469687e+02;
const double a4 = 1.383577518672690e+02;
const double a5 = -3.066479806614716e+01;
const double a6 = 2.506628277459239e+00;
const double b1 = -5.447609879822406e+01;
const double b2 = 1.615858368580409e+02;
const double b3 = -1.556989798598866e+02;
const double b4 = 6.680131188771972e+01;
const double b5 = -1.328068155288572e+01;
const double c1 = -7.784894002430293e-03;
const double c2 = -3.223964580411365e-01;
const double c3 = -2.400758277161838e+00;
const double c4 = -2.549732539343734e+00;
const double c5 = 4.374664141464968e+00;
const double c6 = 2.938163982698783e+00;
const double d1 = 7.784695709041462e-03;
const double d2 = 3.224671290700398e-01;
const double d3 = 2.445134137142996e+00;
const double d4 = 3.754408661907416e+00;
const double pLow = 0.02425;
const double pHigh = 1 - pLow;
double q, x;
if (0 < p && p < pLow)
{
q = Math.Sqrt(-2*Math.Log(p));
x = (((((c1*q+c2)*q+c3)*q+c4)*q+c5)*q+c6) / ((((d1*q+d2)*q+d3)*q+d4)*q+1);
return x;
}
if (pLow <= p && p <= pHigh)
{
q = p - 0.5;
double r = q*q;
x = (((((a1*r+a2)*r+a3)*r+a4)*r+a5)*r+a6)*q / (((((b1*r+b2)*r+b3)*r+b4)*r+b5)*r+1);
return x;
}
//if(p_high < p && p < 1)
q = Math.Sqrt(-2*Math.Log(1-p));
x = -(((((c1 * q + c2) * q + c3) * q + c4) * q + c5) * q + c6) / ((((d1 * q + d2) * q + d3) * q + d4) * q + 1);
return x;
}
#endregion
#region Student's t Distribution
///
/// Returns the PDF of Student's t distribution.
///
/// Value at which the distribution is evaluated.
/// Degrees of freedom.
public static double StudentsTProbabilityDensityFunction(double x, double k)
{
double a = MMath.GammaFunction(0.5 * k + 0.5);
double b = Math.Pow(1.0 + (x * x) / k, -0.5 * k - 0.5);
double c = Math.Sqrt(k * Math.PI) * MMath.GammaFunction(0.5 * k);
return a * b / c;
}
///
/// Returns the CDF of Student's t distribution.
///
/// Value at which the distribution is evaluated.
/// Degrees of freedom.
public static double StudentsTCumulativeDensityFunction(double x, double k)
{
//Tested against Excel for
//df = 1, 2, 3, 4, 5, 10, 20, 30, 40, 50, 100, 200, 30, 400, 500 and
//x = -7.0, -6.9, -6.8, ..., 7.0
//maximum discrepancy was less than 0.04%
if(double.IsNaN(x))
return double.NaN;
if (x == 0)
return 0.5;
if (x > 0)
{
double x2 = k / (x * x + k);
return 1.0 - 0.5 * MMath.RegularizedIncompleteBetaFunction(x2, 0.5 * k, 0.5);
}
return 1.0 - StudentsTCumulativeDensityFunction(-x, k);
}
///
/// Returns the inverse of the CDF of the Student's t distribution.
/// For k = 3, and 5+ the solution is an approximation.
///
/// Cumulative probability of the distribution. p is between 0 and 1.
/// Degrees of freedom.
public static double StudentsTCumulativeDistributionFunctionInverse(double p, double k)
{
//Tested against Excel for k = 1, 2, ... 20, 30, ..., 100 and p = 0.000001, 0.00001, 0.0001, 0.001, 0.01, 0.02, ..., 0.99.
//Max absolute difference was 0.0088 for k = 3 and p = 0.000001 (-103.2906 vs -103.2995).
const int maxIterations = 32;
double maxProbError = Math.Min(0.0001, 0.001 * Math.Min(p, 1 - p)); //max absolute difference between p and p implied by return value.
const double maxXDiff = 0.00000001; //If bisection method is called, it will stop when the difference between the two values it is testing is less than maxXDiff.
if (p < 0 || p > 1.0) throw new ArgumentException("Probability must be between 0 and 1");
if (k <= 0) throw new ArgumentException("Degrees of freedom must be greater than 0.");
if (p == 0) return double.NegativeInfinity;
if (p == 1) return double.PositiveInfinity;
if (p == 0.5) return 0.0;
if (k == 1)
{
return Math.Tan(Math.PI*(p - 0.5));
}
else if (k == 2)
{
double alpha = 4*p*(1 - p);
return 2*(p - 0.5)*Math.Sqrt(2.0/alpha);
}
else if (k == 4)
{
double alpha = 4*p*(1 - p);
double sqrtAlpha = Math.Sqrt(alpha);
double q = Math.Cos((1.0/3.0)*Math.Acos(sqrtAlpha))/sqrtAlpha;
return Math.Sign(p - 0.5)*2*Math.Sqrt(q - 1);
}
else
{
double x;
if (k > 5.5 && (p < 0.017 || p > 0.983))
x = StudentsTCumulativeDistributionFunctionInverseApproximation(p, k); //this method is more accurate for extreme values pf p.
else if (k > 7.875)
x = StandardNormalCumulativeDistributionFunctionInverse(p);
else if (k < 1.375)
x = StudentsTCumulativeDistributionFunctionInverse(p, 1);
else if (k < 2.625)
x = StudentsTCumulativeDistributionFunctionInverse(p, 2);
else
x = StudentsTCumulativeDistributionFunctionInverse(p, 4);
double oldAbsProbError = double.MaxValue;
double oldX = double.NaN;
for (int i = 0; i < maxIterations; i++)
{
double p0 = StudentsTCumulativeDensityFunction(x, k);
double absProbError = Math.Abs(p - p0);
if (absProbError < maxProbError) return x;
if (absProbError > oldAbsProbError)
{
bool reachedDiffLimit = false;
//x = BisectionSearch(x, oldX, k, p, i, maxIterations, maxXDiff, ref reachedDiffLimit, null, null);
x = SecantSearch(x, oldX, k, p, i, maxIterations, maxXDiff, ref reachedDiffLimit, null, null);
if (reachedDiffLimit) return x;
p0 = StudentsTCumulativeDensityFunction(x, k);
absProbError = Math.Abs(p - p0);
if (absProbError < maxProbError) return x;
throw new Exception("Solution did not converge");
}
oldAbsProbError = absProbError;
oldX = x;
double slope = StudentsTProbabilityDensityFunction(x, k);
x = x + (p - p0)/slope;
if ((p < 0.50 && x > 0) || (p > 0.50 && x < 0))
x = 0.5*oldX; //if you overshoot so much that you end up with the wrong sign for x, then take the average of the last x and 0.
}
throw new Exception("Solution did not converge");
}
}
private static double BisectionSearch(double x1, double x2, double k, double p, int splitCount, int maxSplits, double maxXDiff, ref bool reachedDiffLimit, double? e1, double? e2)
{
if (Math.Abs(x1 - x2) < maxXDiff)
{
reachedDiffLimit = true;
return 0.5 * (x1 + x2);
}
double xm = 0.5*(x1 + x2);
e1 = e1 ?? StudentsTCumulativeDensityFunction(x1, k) - p;
double em = StudentsTCumulativeDensityFunction(xm, k) - p;
e2 = e2 ?? StudentsTCumulativeDensityFunction(x2, k) - p;
splitCount++;
if (splitCount > maxSplits)
return xm;
if (Math.Sign(e1.Value) != Math.Sign((em)))
{
xm = BisectionSearch(x1, xm, k, p, splitCount, maxSplits, maxXDiff, ref reachedDiffLimit, e1, em);
}
else if (Math.Sign(e2.Value) != Math.Sign(em))
{
xm = BisectionSearch(x2, xm, k, p, splitCount, maxSplits, maxXDiff, ref reachedDiffLimit, e2, em);
}
return xm;
}
private static double SecantSearch(double x1, double x2, double k, double p, int splitCount, int maxSplits, double maxXDiff, ref bool reachedDiffLimit, double? e1, double? e2)
{
if (Math.Abs(x1 - x2) < maxXDiff)
{
reachedDiffLimit = true;
return 0.5 * (x1 + x2);
}
double xm = e1 == null ? 0.5 * (x1 + x2) : (x1 * e2.Value - x2 * e1.Value) / (e2.Value - e1.Value);
e1 = e1 ?? StudentsTCumulativeDensityFunction(x1, k) - p;
double em = StudentsTCumulativeDensityFunction(xm, k) - p;
e2 = e2 ?? StudentsTCumulativeDensityFunction(x2, k) - p;
splitCount++;
if (splitCount > maxSplits)
return xm;
if (Math.Sign(e1.Value) != Math.Sign((em)))
{
xm = BisectionSearch(x1, xm, k, p, splitCount, maxSplits, maxXDiff, ref reachedDiffLimit, e1, em);
}
else if (Math.Sign(e2.Value) != Math.Sign(em))
{
xm = BisectionSearch(x2, xm, k, p, splitCount, maxSplits, maxXDiff, ref reachedDiffLimit, e2, em);
}
return xm;
}
///
/// This approximation is more accurate for very high and very low values of p.
///
///
///
///
private static double StudentsTCumulativeDistributionFunctionInverseApproximation(double p, double k)
{
if(p < 0.5)
return -StudentsTCumulativeDistributionFunctionInverseApproximation(1 - p, k);
double[] paramArray = StcdfiaParams(p);
return Math.Exp(paramArray[0] + paramArray[1] / k + paramArray[2] / (k * k));
}
///
/// These constants, {A, B, C}, are based on regressions of the form x = A + B/k + C/k^2 for various value of p.
///
private static double[] StcdfiaParams(double p)
{
if(p < 0.994500000)
return new[] {0.844289344098942, 1.60247403675763, 1.21204210950434}; //p = 99%
if(p < 0.999450000)
return new[] {1.12763432956672, 2.6724693160194, 2.75261355460215}; //p = 99.9%
if(p < 0.999945000)
return new[] {1.3109905400331, 3.85826102579765, 4.54414522746746}; //p = 99.99%
if(p < 0.999994500)
return new[] {1.44451562705915, 5.1673200281153, 6.35867233089288}; //p = 99.999%
if(p < 0.999999450)
return new[] {1.54796181658583, 6.59696293772306, 8.07465724847696}; //p = 99.9999%
if(p < 0.999999945)
return new[] {1.63123569060444, 8.13645869984937, 9.64357980688031}; //p = 99.99999%
return new[] { 1.70006767767707, 9.77199585223768, 11.0556716544646 }; //p = 99.999999%
}
#endregion
#region Uniform Distribution
///
/// Returns the PDF of the uniform distribution.
///
/// Value at which the distribution is evaluated.
/// Minimum value of the distribution.
/// Maximum value of the distribution.
public static double UniformProbabilityDensityFunction(double x, double min, double max)
{
if (max <= min) throw new ArgumentException("max must be greater than min");
if (x < min || x > max) return 0.0;
return 1.0 / (max - min);
}
///
/// Returns the CDF of the uniform distribution.
///
/// Value at which the distribution is evaluated.
/// Minimum value of the distribution.
/// Maximum value of the distribution.
public static double UniformCumulativeDensityFunction(double x, double min, double max)
{
if (max <= min) throw new ArgumentException("max must be greater than min");
return (x - min) / (max - min);
}
#endregion
///
/// The Irwin-Hall distribution results from the sum on n independent standard uniform variables
///
/// The value at which to evaluate the distribution.
/// The number of standard uniform variables.
public static double IrwinHallProbabilityDensityFunction(double x, int n)
{
if (x < 0 || x > n) return 0.0;
double d = 0;
for (int i = 0; i <= n; i++)
d += Math.Pow(-1, i) * MMath.Combinations(n, i) * Math.Pow(x - i, n - 1) * Math.Sign(x - i);
d *= 0.5;
d /= MMath.Factorial(n - 1);
return d;
}
///
/// Returns the PDF of the lognormal distribution.
///
/// Value at which the distribution is evaluated.
/// This is *not* the mean. See Miller (2012).
/// This is *not* the standard deviation. See Miller (2012).
public static double LogNormalProbabilityDensityFunction(double x, double mu, double sigma)
{
const double sqrtTwoPiInv = 0.398942280401433;
double z = (Math.Log(x) - mu) / sigma;
return sqrtTwoPiInv * Math.Exp(-0.5 * z * z) / (x * sigma);
}
///
/// The PDF of the Poisson distribution.
/// The probability of n iid events occuring in an interval, if the expected number of event is lambda.
///
/// Equals both the mean and variance of the distribution.
/// The number of events.
public static double PoissonProbabilityDensityFunction(double lambda, int n)
{
return Math.Pow(lambda, n) * Math.Exp(-lambda) / MMath.Factorial(n);
}
///
/// Maps the input array to a standard normal distribution.
/// The points in the output array have the same rank order as the input array, but are normally distributed.
/// Note, if two points in the input array have the same value, the values in the output array will be different.
///
public static double[] NormalizeArray(double[] array)
{
List points = new List();
int n = array.Length;
for(int i = 0; i < n; i++)
points.Add(new NormalizePoint(i, array[i]));
points.Sort((a,b) => a.OriginalValue.CompareTo(b.OriginalValue));
double interval = 1.0 / n;
for(int i = 0; i < n; i++)
points[i].NormalizedValue = StandardNormalCumulativeDistributionFunctionInverse(interval * (i + 0.5));
points.Sort((a, b) => a.OriginalIndex.CompareTo(b.OriginalIndex));
double[] outArray = new double[n];
for(int i = 0; i < n; i++)
outArray[i] = points[i].NormalizedValue;
return outArray;
}
private class NormalizePoint
{
public readonly int OriginalIndex;
public readonly double OriginalValue;
public double NormalizedValue;
public NormalizePoint(int originalIndex, double value)
{
OriginalIndex = originalIndex;
OriginalValue = value;
}
}
}
}
//Disclaimer
//This code is freeware. The methods are not proprietary. Feel free to use, modify and redistribute. That said, if you plan
//to use or redistribute give credit where credit is due and provide a link back to Risk256.com (or don't remove the link
//and references already in the code). The code is intended primarily as an educational tool. No warranty is made as to the
//code's accuracy. Use at your own risk.