常態分布計算器

請輸入數值並點擊計算

公式:

---

Normal Distribution Probability Density Function(PDF):
　　(Math.E**(-((x-mean)**2)/(2*(stdev**2))))/((2*Math.PI*(stdev**2))**0.5)

---

Standard Normal Deviate(z-score, z, Z):
　　(x-mean)/stdev

---

Error Function(erf):
　　2/(Math.PI**0.5) * [integral of X:(0,n) Function(X):Math.E**(-(X**2))]

---

Error Function Approximation:
　　Abramowitz and Stegun give several approximations of varying accuracy (equations 7.1.25–28).
　　Maximum Error 1.5*(10**-7)
　　Note: -erf(x) = erf(-x), x>0.
　　a1 = 0.254829592
　　a2 = -0.284496736
　　a3 = 1.421413741
　　a4 = -1.453152027
　　a5 = 1.061405429
　　p = 0.3275911
　　t = 1.0 / (1.0 + p * x)
　　Note: 秦九韶算法，又稱霍納法則 (Horner's Method)
　　result = 1.0 - (((((a5 * t + a4) * t) + a3) * t + a2) * t + a1) * t * Math.exp(-x * x);

---

Normal Distribution Cumulative Distribution Function(CDF, integral of PDF):
　　0.5*(1+erf(Z/(2**0.5)))

計算 X1 到 X2 的機率

計算 P(< z) 與 P(> z)

計算 Z1 到 Z2 之機率

計算 Z

公式:

---

Error Function Approximation:
//Abramowitz and Stegun give several approximations of varying accuracy (equations 7.1.25–28).
//https://en.wikipedia.org/wiki/Error_function
//maximum error 1.5*(10**-7)
function approx_erf(x) {
  // -erf(x) = erf(-x), x>0.
  const sign = (x >= 0) ? 1 : -1;
  x = Math.abs(x);
  
  const a1 =  0.254829592;
  const a2 = -0.284496736;
  const a3 =  1.421413741;
  const a4 = -1.453152027;
  const a5 =  1.061405429;
  const p  =  0.3275911;

  const t = 1.0 / (1.0 + p * x);
  //秦九韶算法，又稱霍納法則 (Horner's Method)
  const y = 1.0 - (((((a5 * t + a4) * t) + a3) * t + a2) * t + a1) * t * Math.exp(-x * x);

  return sign * y;
}

---

Normal Distribution Cumulative Distribution Function(CDF, integral of PDF):
function NormalDistributionCumulativeDistributionFunction(z){
    return 0.5*(1+approx_erf(z/(2**0.5)));
}

---

改寫Normal Distribution Cumulative Distribution Function到牛頓法F(X)=0:

let sign;
let CDF=0.84; //CDF is user input

function FDE(x){ //x=z/(2**0.5)
  sign=1;
  if((CDF*2-1)<0){sign=-1}
  console.log("sign",sign)
  const a1 =  0.254829592;
  const a2 = -0.284496736;
  const a3 =  1.421413741;
  const a4 = -1.453152027;
  const a5 =  1.061405429;
  const p  =  0.3275911;
  const t = 1.0 / (1.0 + p * x);
  return (((((a5 * t + a4) * t) + a3) * t + a2) * t + a1) * t * Math.exp(-x * x) + (CDF*2-1)*sign - 1
}

---

對改寫的F(X)進行微分:

function IDE(N,M,x){
const p  =  0.3275911;
let result = 
    ( 
        ( ( (-2*x)*(Math.E**(-(x**2))) ) * N ) * ( (1.0 + p * x)**(-M) ) - 
     ( N*( Math.E**(-(x**2)) ) ) * ( M * p * ( (1.0 + p * x)**(-M-1) ) )
    ) 
;
return result;
}

function BDE(x){
  const a1 =  0.254829592;
  const a2 = -0.284496736;
  const a3 =  1.421413741;
  const a4 = -1.453152027;
  const a5 =  1.061405429;
  let result = IDE(a1,1,x)+IDE(a2,2,x)+IDE(a3,3,x)+IDE(a4,4,x)+IDE(a5,5,x);
  return result;
}

---

牛頓法 f(x)=0, new x = x - f(x)/((d/dx)f(x)) :
function newtonMethod(f, df, x0, tolerance = 1e-7, maxIterations = 50) {
    let x = x0;
    
    for (let i = 0; i < maxIterations; i++) {
        const fx = f(x);
        const dfx = df(x);
        
        // Avoid division by zero if the derivative is too small
        if (Math.abs(dfx) < 1e-10) {
            console.log(0,x,fx,dfx)
            console.error("Derivative is too small; method may fail.");
            return "Derivative is too small; method may fail.";
        }

        const xNext = x - fx / dfx;
        let dx = Math.abs(xNext - x)
        console.log(x, xNext, dx);
        // Check if the result is close enough to the root
        if (dx < tolerance) {
            let RT=Math.abs(f(xNext));
            console.log(`Converged in ${i + 1} iterations.`, RT,xNext,tolerance);
            if(RT>tolerance){
                console.warn("Converged but F not equal to Zero. ERROR OUTPUT.");
            }
            return xNext;
        }

        x = xNext;
        
    }

    console.warn("Reached maximum iterations without full convergence.");
    return "Maximum Iterations ERROR";
}

---

function inverseNormalDistributionCumulativeDistributionFunction(){
    let x1 = newtonMethod(FDE,BDE,0)*sign;
    let z1 = x1*(2**0.5);
    return z1;
}

z1 = inverseNormalDistributionCumulativeDistributionFunction(CDF);

---

let x1;
let stdev = 10; //Standard Deviation
let u = 3; //Mean aka μ, expected value, weighted average.

x1 = z1 * stdev + u

計算 Z