As already correctly stated by Shane, the equation is an implementation of the Taylor Expansion of the normal cdf. The sum
value iterates above and below the "real" value with increasing precision. If the value is close to 1 or 0 there is a very low, but existing, probability that sum
will be >1 or <0, because of the (relatively) early break by loopstop
.
The deviation is further strengthened by rounding 1/Math.sqrt(2*Math.Pi)
to 0.3989422804
and the precision issues of javascript float numbers. Additionally, the provided solution will not work for z-scores >7 or <-7
I updated the code to be more accurate using the decimal.js npm library and to directly return the p-value:
function GetpValueFromZ(_z, type = "twosided")
{
if(_z < -14)
{
_z = -14
}
else if(_z > 14)
{
_z = 14
}
Decimal.set({precision: 100});
let z = new Decimal(_z);
var sum = new Decimal(0);
var term = new Decimal(1);
var k = new Decimal(0);
var loopstop = new Decimal("10E-50");
var minusone = new Decimal(-1);
var two = new Decimal(2);
let pi = new Decimal("3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117067982148086513282306647")
while(term.abs().greaterThan(loopstop))
{
term = new Decimal(1)
for (let i = 1; i <= k; i++) {
term = term.times(z).times(z.dividedBy(two.times(i)))
}
term = term.times(minusone.toPower(k)).dividedBy(k.times(2).plus(1))
sum = sum.plus(term);
k = k.plus(1);
}
sum = sum.times(z).dividedBy(two.times(pi).sqrt()).plus(0.5);
if(sum.lessThan(0))
sum = sum.abs();
else if(sum.greaterThan(1))
sum = two.minus(sum);
switch (type) {
case "left":
return parseFloat(sum.toExponential(40));
case "right":
return parseFloat((new Decimal(1).minus(sum)).toExponential(40));
case "twosided":
return sum.lessThan(0.5)? parseFloat(sum.times(two).toExponential(40)) : parseFloat((new Decimal(1).minus(sum).times(two)).toExponential(40))
}
}
By increasing the Decimal.js precision
value and decreasing the loopstop
value you can get accurate p-values for very small (or very high) z-scores for the cost of calculation time.
term =
line makes this very difficult to read. – Wimble