Created
April 12, 2012 15:24
-
-
Save chriskite/2368163 to your computer and use it in GitHub Desktop.
Bernoulli PHP
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
#!/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"; | |
} |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment