gsl_statext.c

常用的数学统计算法

#include <assert.h>
#include <gsl/gsl_sys.h>
#include <gsl/gsl_math.h>

#include "gsl_statext.h"

// Calculate the approximation of the standard normal distribution
inline static double standardNormalF (const double y) {
  const double a1 = 0.319381530;
  const double a2 = -0.356563782;
  const double a3 = 1.781477937;
  const double a4 = -1.821255978;
  const double a5 = 1.330274429;
  
  double result = a4 + a5 * y;
  result = a3 + y * result;
  result = a2 + y * result;
  result = a1 + y* result;
  result *= y;
  
  return result;
}
  
// Calculate the standard normal distribution
double standardNormalDistribution (const double x) {
  double result = 0;
  
  if (x >= 4) {
	result = 1;
	
  } else {
	if (x <= -4) {
	  result = 0;
	  
	} else {
	  const double b = 0.2316419;
	  // c = 1 / sqrt (2.PI) = 0.398942208
	  const double c = 1 / sqrt (2 * M_PI);
	  const double standNormProb = exp (-x*x / 2) * c;
	  
	  if (x > 0) {
		const double y = 1 / (1 + b * x);
		result = 1 - standNormProb * standardNormalF (y);
		
	  } else {
		const double y = 1 / (1 - b*x);
		result = standNormProb * standardNormalF (y);
	  }
	} 
  }
  
  return result;
}   

// Calculate the standard normal distribution
double normalDistribution (const double x, const double mean,
						   const double stdDev) {
  const double xmod = (x-mean) / stdDev;
  
  double result = sqrt (stdDev);
  result *= standardNormalDistribution (xmod);
  
  return result;
}

// sqrt (-2 ln(p))
inline static double lnP2 (const double p) {
  double result = -2 * log (p);
  result = sqrt (result);
}

// Calculate the standard normal distribution
inline static double inverseNormalDistributionFract1 (const double q) {
  const double c1 = -0.007784894002430293;
  const double c2 = -0.3223964580411365;
  const double c3 = -2.400758277161838;
  const double c4 = -2.549732539343734;
  const double c5 = 4.374664141464968;
  const double c6 = 2.938163982698783;
  const double d1 = 0.007784695709041462;
  const double d2 = 0.3224671290700398;
  const double d3 = 2.445134137142996;
  const double d4 = 3.754408661907416;

  double nominator = c1;
  nominator = c2 + q * nominator;
  nominator = c3 + q * nominator;
  nominator = c4 + q * nominator;
  nominator = c5 + q * nominator;
  nominator = c6 + q * nominator;
  
  double denominator = d1;
  denominator = d2 + q * denominator;
  denominator = d3 + q * denominator;
  denominator = d4 + q * denominator;
  denominator = 1 + q * denominator;
  
  const double result = nominator / denominator;
  return result;
}
  
// Calculate the standard normal distribution
inline static double inverseNormalDistributionFract2 (const double q) {
  const double a1 = -39.69683028665376;
  const double a2 = 220.9460984245205;
  const double a3 = -275.9285104469687;
  const double a4 = 138.3577518672690;
  const double a5 = -30.66479806614716;
  const double a6 = 2.506628277459239;
  const double b1 = -54.47609879822406;
  const double b2 = 161.5858368580409;
  const double b3 = -155.6989798598866;
  const double b4 = 66.80131188771972;
  const double b5 = -13.28068155288572;
  
  // r = q^2
  const double r = q * q;
  
  double nominator = a1;
  nominator = a2 + r * nominator;
  nominator = a3 + r * nominator;
  nominator = a4 + r * nominator;
  nominator = a5 + r * nominator;
  nominator = a6 + r * nominator;
  
  double denominator = b1;
  denominator = b2 + r * denominator;
  denominator = b3 + r * denominator;
  denominator = b4 + r * denominator;
  denominator = b5 + r * denominator;
  denominator = 1 + r * denominator;
  
  const double result = nominator / denominator;
  return result;
}

// Calculate the standard normal distribution
double inverseStandardNormalDistribution (const double p) {
  double result = 0;
  
  const double pl = 0.02425;
  const double ph = 1 - pl;
  
  if (p == 0) {
	result = gsl_neginf();
	
  } else if (p < pl) {
	const double q = lnP2 (p);
	result = inverseNormalDistributionFract1 (q);
	
  } else if (p <= ph) {
	const double q = p - 0.5;
	result = inverseNormalDistributionFract2 (q);
	
  } else if (p < 1) {
	const double q = lnP2 (1-p);
	result = -inverseNormalDistributionFract1 (q);
	
  } else {
	assert (p == 1);
	result = gsl_posinf();
  }
  
  return result;
}    

// Calculate the standard normal density
double standardNormalDensity (const double x) {
  double result = - x * x / 2;
  
  result = exp (result);
  result = result / sqrt (2 * M_PI);
  
  return result;
}

// Calculate the normal density
double normalDensity (const double x, const double mean,
					  const double stdDev) {
  const double xmod = (x-mean) / stdDev;
  
  double result = 1 / sqrt (stdDev);
  result *= standardNormalDensity (xmod);
  
  return result;
}