package de.psychometrica.norming;

/**
 * Java Implementation of the inverse normal distribution function according to
 * the algorithm, defined by Peter John Acklam (pjacklam@online.no) on
 * http://home.online.no/~pjacklam/notes/invnorm/. The license on the page is
 * "Use them as you wish, for whatever purpose." When comparing the precision of
 * the algorithm with the function Microsoft Excel, the mean square sum of
 * deviation between p = .001 and p = .999 amounts to 6,2144E-19. An additional
 * refinement is therefore not necessary for the use in psychometrics.
 * Precompiled binaries and Javadoc can be found at
 * http://www.psychometrica.de/norming.html <br>
 * The code is licensed under the Lesser General Public License (LGPL 2.0)
 * 
 * @author PD Dr. Wolfgang Lenhard; Original algorithm by Peter John Acklam
 *
 */
public class NormalDistribution {
	private final static double[] A = { -3.969683028665376e+01,
			2.209460984245205e+02, -2.759285104469687e+02,
			1.383577518672690e+02, -3.066479806614716e+01,
			2.506628277459239e+00 };
	private final static double[] B = { -5.447609879822406e+01,
			1.615858368580409e+02, -1.556989798598866e+02,
			6.680131188771972e+01, -1.328068155288572e+01 };
	private final static double[] C = { -7.784894002430293e-03,
			-3.223964580411365e-01, -2.400758277161838e+00,
			-2.549732539343734e+00, 4.374664141464968e+00,
			2.938163982698783e+00 };
	private final static double[] D = { 7.784695709041462e-03,
			3.224671290700398e-01, 2.445134137142996e+00, 3.754408661907416e+00 };
	private final static double[] BREAKPOINT = { 0.02425, 0.97575 };

	/**
	 * @param p
	 *            The probability, e. g. the percent rank in psychodiagnostics
	 * @return the according standardized value, z
	 */
	public static double getInverseCumulativeDistribution(final double p) {
		if (0 < p && p < BREAKPOINT[0]) {
			double q = Math.sqrt((-2 * Math.log(p)));
			return (((((C[0] * q + C[1]) * q + C[2]) * q + C[3]) * q + C[4])
					* q + C[5])
					/ ((((D[0] * q + D[1]) * q + D[2]) * q + D[3]) * q + 1d);
		} else if (BREAKPOINT[0] < p && p < BREAKPOINT[1]) {
			double q = p - .5;
			double r = q * q;
			return (((((A[0] * r + A[1]) * r + A[2]) * r + A[3]) * r + A[4])
					* r + A[5])
					* q
					/ (((((B[0] * r + B[1]) * r + B[2]) * r + B[3]) * r + B[4])
							* r + 1d);
		} else {
			double q = Math.sqrt((-2 * Math.log(1 - p)));
			return -(((((C[0] * q + C[1]) * q + C[2]) * q + C[3]) * q + C[4])
					* q + C[5])
					/ ((((D[0] * q + D[1]) * q + D[2]) * q + D[3]) * q + 1d);
		}


	}

	/**
	 * Calculates the density of the normal distribution function
	 * 
	 * @param z
	 *            the standardized value of a raw value with m = 0.0 and sd =
	 *            1.0
	 * @return the density of the normal distribution
	 */
	public static double getNormalDistributionDensity(final double z) {
		return Math.exp((z * z) / (-2.0d)) / 2.506628274631000502415765284811d;
	}

	/**
	 * Calculate the normal cumulative density function according to formula 4.5
	 * from Vazquez-Leal, H., Castaneda-Sheissa, R., Filobello-Nino, U.,
	 * Sarmiento-Reyes, A., & Sanchez Orea, J. (2012). High Accurate Simple
	 * Approximation of Normal Distribution Integral. Mathematical Problems in
	 * Engineering, 2012, 1–22. doi:10.1155/2012/124029
	 * 
	 * The function has a precision of approx. 10E-4, when comparing it to the
	 * implementation in MS Excel
	 * 
	 * @param x
	 *            Standardized value (m = .0, sd = 1.0)
	 * @return cumulative distribution function value; probability, e. g.
	 *         percent rank
	 */
	public static double getCumulativeDistribution(double x) {
		return 1d / (Math.exp(-15.56521739130435 * x + 111d
				* Math.tanh(.1258503401360544 * x)) + 1d);
	}
}