Error function and Rarity calculation with Gauss distribution.

rarity.pl

use strict;
use warnings;
use bignum (a => 200);

# error function
sub erf {
    my $x = shift;
    $x = $x * sqrt(2);
    my $sqrt2pi = sqrt(2 * 3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679821480865132823066470938446095505822317253594081284811174502841027);
    my $t = 1.0/(1.0 + 0.2316419*abs($x));
    my $y = $t*(0.319381530
                    + $t*(-0.356563782
                          + $t*(1.781477937
                                + $t*(-1.821255978
                                      + $t*1.330274429 ))));
    my $f = 1.0 - ( exp(-0.5 * $x * $x)/($sqrt2pi) ) * $y;
    return (2 * $f -1)
}

# reciprocal of upper cumulative value
sub rarity {
    my ($val,$mu,$std) = (0,0,1);
    $val = shift if @_;
    $mu = shift if @_;
    $std = shift if @_;
    return 1 / (1 - ( 1/2 * (1 + &erf(($val - $mu) / sqrt(2 * $std * $std) ) ) ));
}

# example: calculating rarity of MENSA's acceptance score
# IQ>=130, mu=100, stddev=15
print rarity(130,100,15);