package com.hfg.math;


//------------------------------------------------------------------------------
/**
 For calculations related to standard normal distributions.

 @author J. Alex Taylor, hairyfatguy.com
 */
//------------------------------------------------------------------------------
// com.hfg XML/HTML Coding Library
//
// This library 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 2.1 of the License, or (at your option) any later version.
//
// This library 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 should have received a copy of the GNU Lesser General Public
// License along with this library; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
//
// J. Alex Taylor, President, Founder, CEO, COO, CFO, OOPS hairyfatguy.com
// jataylor@hairyfatguy.com
//------------------------------------------------------------------------------

public class StandardNormalDistribution
{
   private Double mMean;
   private Double mSampleStandardDeviation;

   // Coefficients in rational approximations for the inverse normal cumulative distribution function
   private static final double[] sA_coeff = new double[] {
         -3.969683028665376e+01,  2.209460984245205e+02,
         -2.759285104469687e+02,  1.383577518672690e+02,
         -3.066479806614716e+01,  2.506628277459239e+00 };


   private static final double[] sB_coeff = new double[] {
         -5.447609879822406e+01,  1.615858368580409e+02,
         -1.556989798598866e+02,  6.680131188771972e+01,
         -1.328068155288572e+01 };


   private static final double[] sC_coeff = new double[] {
         -7.784894002430293e-03, -3.223964580411365e-01,
         -2.400758277161838e+00, -2.549732539343734e+00,
         4.374664141464968e+00,  2.938163982698783e+00 };

   private static final double[] sD_coeff = new double[] {
         7.784695709041462e-03, 3.224671290700398e-01,
         2.445134137142996e+00,  3.754408661907416e+00 };

   //--------------------------------------------------------------------------
   public StandardNormalDistribution()
   {

   }

   //--------------------------------------------------------------------------
   public StandardNormalDistribution setMean(Double inValue)
   {
      mMean = inValue;
      return this;
   }

   //--------------------------------------------------------------------------
   public Double getMean()
   {
      return mMean;
   }


   //--------------------------------------------------------------------------
   public StandardNormalDistribution setSampleStandardDeviation(Double inValue)
   {
      mSampleStandardDeviation = inValue;
      return this;
   }

   //--------------------------------------------------------------------------
   public Double getSampleStandardDeviation()
   {
      return mSampleStandardDeviation;
   }


   //--------------------------------------------------------------------------
   public double getZScore(double inValue)
   {
      // Z = (X - μ) / σ
      return (inValue - mMean) / mSampleStandardDeviation;
   }

   //--------------------------------------------------------------------------
   public double getProbability(double inValue)
   {
      double zScore = getZScore(inValue);

      double probability = 0;

      if (zScore < -7)
      {
         probability = 0.0;
      }
      else if (zScore > 7)
      {
         probability = 1.0;
      }
      else
      {
         boolean invert = (zScore > 0.0);

         zScore = Math.abs(zScore);
         double b = 0.0;
         double s = Math.sqrt(2) / 3 * zScore;
         double HH = .5;
         for (int i = 0; i < 12; i++)
         {
            double a = Math.exp(-HH * HH / 9) * Math.sin(HH * s) / HH;
            b = b + a;
            HH = HH + 1.0;
         }

         probability = .5 - b / Math.PI;

         if (invert)
         {
            probability = 1.0 - probability;
         }
      }

      return probability;
   }

   //--------------------------------------------------------------------------
   /**
    Given a probability P, it returns an approximation to the X satisfying P = Pr{Z &lt;= X} where Z is a
    random variable from the standard normal distribution. Uses the lower tail quantile for standard normal distribution function from
    <a href='http://home.online.no/~pjacklam/notes/invnorm/'>http://home.online.no/~pjacklam/notes/invnorm/</a>.
    */
   public double getValueThresholdForProbability(float inProbability)
   {
      double z = 0;

      // Define break-points.
      double pLow  = 0.02425;
      double pHigh = 1 - pLow;

      if ( inProbability < pLow )
      {
         // Rational approximation for lower region:
         double q  = Math.sqrt(-2*Math.log(inProbability));
         z = (((((sC_coeff[0]*q+sC_coeff[1])*q+sC_coeff[2])*q+sC_coeff[3])*q+sC_coeff[4])*q+sC_coeff[5]) /
               ((((sD_coeff[0]*q+sD_coeff[1])*q+sD_coeff[2])*q+sD_coeff[3])*q+1);
      }
      else if ( pHigh < inProbability )
      {
         // Rational approximation for upper region:
         double q  = Math.sqrt(-2*Math.log(1 - inProbability));
         z = -(((((sC_coeff[0]*q+sC_coeff[1])*q+sC_coeff[2])*q+sC_coeff[3])*q+sC_coeff[4])*q+sC_coeff[5]) /
               ((((sD_coeff[0]*q+sD_coeff[1])*q+sD_coeff[2])*q+sD_coeff[3])*q+1);
      }
      else
      {
         // Rational approximation for central region:
         double q = inProbability - 0.5;
         double r = q * q;
         z = (((((sA_coeff[0] * r + sA_coeff[1]) * r + sA_coeff[2]) * r + sA_coeff[3]) * r + sA_coeff[4]) * r + sA_coeff[5]) * q /
               (((((sB_coeff[0] * r + sB_coeff[1]) * r + sB_coeff[2]) * r + sB_coeff[3]) * r + sB_coeff[4]) * r + 1);
      }

      // X = μ + Zσ
      return getMean() + (z * getSampleStandardDeviation());
   }
}