Bernoulli PHP

bernoulli.php

#!/usr/bin/env php

<?php

/*
 * @param float $p should be 0= < $p <= 1
 * @return 1 with probability $p, 0 with probability (1 - $p)
 */
function bernoulli($p) {
    $r = mt_rand() / mt_getrandmax();
    return $r <= $p;
}

/* Statistical testing for the bernoulli function */

$Z = 2.6; // 99% confidence
$err = 0.001;
$pr = 0.25; // change this to test how it performs on different weights

$trials = pow($Z, 2) / (4 * pow($err, 2));

echo "$trials trials\n";

$heads = 0;
$tails = 0;

for($i = 0; $i < $trials; $i++) {
    if(1 == bernoulli($pr)) {
        $heads++;
    } else {
        $tails++;
    }
}

echo "$heads heads, $tails tails\n";

$p = $heads / $trials;
$q = 1 - $p;

echo "Observed P(heads) = $p\n";

$sd = sqrt($p * $q / $trials);
$se = $Z * $sd;

$lb = $p - $se;
$ub = $p + $se;

echo "99% CI: ($lb, $ub)\n";

$se_diff = sqrt( pow($p * $q / $trials, 2) + pow($pr * (1 - $pr) / $trials, 2));

echo "SE_diff: $se_diff\n";

$low = $pr - $p - $se_diff * $Z;
$high = $pr - $p + $se_diff * $Z;

echo "Diff Low: $low Diff High: $high\n";

if($low > -$err && $high < $err) {
    echo "Difference from true $pr is acceptable at 99% confidence\nPASS\n";
} else {
    echo "Difference from true $pr is NOT acceptable at 99% confidence\nFAIL\n";
}