BINOMVERT oder BINOMDIST Berechnung

[knuffiwuffi]

Hallo,

ich suche eine Möglichkeit mit PHP die in Excel verfügbare Formel BINOM.VERT nachzustellen. Leider laufen die Werte bei N über K aus dem Wertebereich und PHP bricht die Berechnung ab. In Excel kann ich die Werte näherungsweise berechnen. Wie wahrscheinlich ist es, daß bei n unabhängigen Versuchen k Erfolge gemessen werden, wenn die Wahrscheinlichkeit für einen Erfolg im Schnitt bei p liegt?

Beispiel in Excel:

n = Anzahl unabhängiger Versuche
k = Anzahl Erfolge
p = Wahrscheinlichkeit für den Erfolg

n = 13346;
k = 6929,5;
p = 0,5192;

=BINOM.VERT(6929,5;13346;0,5192;WAHR) = 0,01 also 1%

Beispiel PHP:

$n = Anzahl unabhängiger Versuche
$k = Anzahl Erfolge
$p = Wahrscheinlichkeit für den Erfolg

$binomdist = binomdist($n, $k, $p); 

function nChooseK($n, $k)
{
    $top    = 1;
    $bottom = 1;

    for ($i = 0; $i < $k; $i++) {
        $top    = bcmul($top, $n - $i);
        $bottom = bcmul($bottom, $i + 1);
    }

    return bcdiv($top, $bottom);
}

function binomdist($n, $k, $p, $cumulative = false)
{
    $res = 0;

    if ($cumulative) {
        for ($i = 0; $i <= $k; $i++) {
            $res += nChooseK($n, $i) * pow($p, $i) * pow(1 - $p, $n - $i);
        }
    } else {
        $res = nChooseK($n, $k) * pow($p, $k) * pow(1 - $p, $n - $k);
    }

    return $res;
} 

[jspit]

Für solch große Werte für n und k brauchst du einen anderen Algorithmus.
Und du hättest ruhig dazuschreiben können, das der Code nicht von dir ist und das hier schon mal Thema war (Binomverteilung als PHP Formel)

[knuffiwuffi]

Also...bevor sich hier jemand falsch verstanden fühlt. Ich hatte dieses Problem schon einmal vor längerer Zeit gestellt. Richtig! Die Lösung ist nicht auf meinen Mist gewachsen, ich hatte aber schon ähnliche Lösungsansätze. Damals habe ich die Programmierung an dem Script eingestellt, da ich keine Lösung gefunden haben (selber, als auch über nette Coder).

Gibt es keine fertigen Scripts für diese Art von Problem? Excel schafft es doch auch den Binomialk. zu berechnen. Ich wüsste nicht, wie ich die Sache weiter angehen soll (mathematisch als auch in PHP). Daher wende ich mich ja an dieses Forum.

Ich bin nicht in der Erwartungshaltung für Fremdleistung, würde es ja selber gerne verstehen und lösen. Bin aber immer noch an dem "no go" Punkt.

[alxy]

Und in dem Verlibkten Beitrag wurde ja auch schon auf das Werteproblem eingegangen: http://www.php.de/forum/webentwicklu...48#post1112848

[Gast]

Ich hatte dieses Problem schon einmal vor längerer Zeit gestellt. Richtig!

das macht es in meine Augen sogar annähenrd verwerflich, dass Du den Thread nicht verlinkt hast.

[alxy]

Warum verwendest du nicht den von mermshaus geposteten, geupdateten Code mit einem entsprechend hohen bcscale() Wert? Deine zweite Funktion verwendet noch nicht die bcmath lib.

[jspit]

knuffiwuffi : Du solltest den Code vom Marc nicht nur abschreiben, sondern auch seine ausführlichen Erläuterungen im verlinkten Thread auch versuchen zu verstehen um eigene Test's zu machen.
Nimm deine Beispielzahlen und berechne mal die einzelnen Ausdrücke wie p ^ k mit einen guten Taschenrechner (calc vom Win7 macht es auch).
Du wirst verwundert feststellen, dass du für dein p und für ein k von 6900 mehr als 2000 Nachkommastellen brauchst um nicht 0 als Resultat zu erhalten!

[knuffiwuffi]

Ich verstehe schon, daß sehr große Werte zu einem Problem in der Berechnung führen. Zähler wird sehr groß, Nenner sehr klein. Das ist ja das Problem. Vielleicht reicht es auch einfach eine obere Schranke mit der Chernoff Ungleichung zu berechnen. Mir reicht ein ungefährer Wert. Wie schafft Excel denn die Berechnung? Die Chernoff Ungleichung müßte doch für eine große Anzahl an Versuchen immer genauer einschränken, oder?

[jspit]

Zähler wird sehr groß, Nenner sehr klein. Das ist ja das Problem.

Das trifft es nicht.

$res = nChooseK($n, $k) * pow($p, $k) * pow(1 - $p, $n - $k); 

Hier wird nChooseK($n, $k) riesig und pow($p, $k) geht gegen 0. Auf 64bit Systemen wird nChooseK($n, $k) noch berechnet, pow($p, $k) wird bei mir da 0, Resultat damit auch.

[knuffiwuffi]

Ich werde versuchen eine Lösung auf Basis der Approx. Normalverteilung zu finden, da eine exakte Berechnung über die Binomial.Vert bei großen Werten derzeit Probleme macht.

[knuffiwuffi]

Hallo,

ich habe nun ein Script, mit Hilfe der bereits durch nette PHPler bereitgestellten Scripte, geschrieben und würde mich über ein Feedback freuen. Überprüft habe ich die Werte mit GeoGebra (kostenlos). Ziel ist es herauszufinden, wie viele Personen ein schlechteres bzw. besseres Ergebnis bei einem Münzwurf erreichen. Auch würde ich gerne Werte für eine Glockenkurve erhalten, um die ungefähre Position des Ergebnisses eintragen zu können. Bislang war ja mein Problem, daß bei sehr großen Stichproben die Binomial.Vert nicht zu berechnen war. Für solche Fälle habe ich dann auf die Approx. Normalverteilung zurückgegriffen. Würde das so passen?:

Beispiel 1:
100 Münzwürfe
48 es fällt Kopf
0.5 Wahrscheinlichkeit für Kopf

Ergebnis:
percentage of how many hit 48 successe within 100 trials or worse: 42.08 %
percentage of how many did better: 57.92 %

Beispiel 2:
3000 Münzwürfe
1488 es fällt Kopf
0.5 Wahrscheinlichkeit für Kopf

Ergebnis:
percentage of how many hit 1488 successe within 3000 trials or worse: 34.4 %
percentage of how many did better: 65.6 %

// ||||||||||||||||||||||||||||||||||||||||||||||||||
// data
// ||||||||||||||||||||||||||||||||||||||||||||||||||
$n = 3000; // independet trials
$k = 1488;  // successes
$p = 0.5;  // equity for success

// ||||||||||||||||||||||||||||||||||||||||||||||||||
// functions
// ||||||||||||||||||||||||||||||||||||||||||||||||||

// roundUp
function roundUp($number, $precision)
{
    $fig = (int) str_pad('1', $precision, '0');
    return (ceil($number * $fig) / $fig);
}

// roundDown
function roundDown($number, $precision)
{
    $fig = (int) str_pad('1', $precision, '0');
    return (floor($number * $fig) / $fig);
}

// hex
function bchexdec($hex) {
    if(strlen($hex) == 1) {
        return hexdec($hex);
    } else {
        $remain = substr($hex, 0, -1);
        $last = substr($hex, -1);
        return bcadd(bcmul(16, bchexdec($remain)), hexdec($last));
    }
}

// dec
function bcdechex($dec) {
    $last = bcmod($dec, 16);
    $remain = bcdiv(bcsub($dec, $last), 16);

    if($remain == 0) {
        return dechex($last);
    } else {
        return bcdechex($remain).dechex($last);
    }
}

// binomdist
function nChooseK($n, $k)
{
    $top    = 1;
    $bottom = 1;

    for ($i = 0; $i < $k; $i++) {
        $top    = bcmul($top, $n - $i);
        $bottom = bcmul($bottom, $i + 1);
    }

    return bcdiv($top, $bottom);
}

function binomdist($n, $k, $p, $cumulative = false)
{
    $res = 0;

    if ($cumulative) {
        for ($i = 0; $i <= $k; $i++) {
            $res += nChooseK($n, $i) * pow($p, $i) * pow(1 - $p, $n - $i);
        }
    } else {
        $res = nChooseK($n, $k) * pow($p, $k) * pow(1 - $p, $n - $k);
    }

    return $res;
}

// approximate normal distribution z score
function erf($x)
{
    $pi = pi();
    $a = (8*($pi - 3))/(3*$pi*(4 - $pi));
    $x2 = $x * $x;

    $ax2 = $a * $x2;
    $num = (4/$pi) + $ax2;
    $denom = 1 + $ax2;

    $inner = (-$x2)*$num/$denom;
    $erf2 = 1 - exp($inner);

    return sqrt($erf2);
}

function cdf($n)
{
    if($n < 0)
    {
        return (1 - erf($n / sqrt(2)))/2;
    }
    else
    {
        return (1 + erf($n / sqrt(2)))/2;
    }
}

function zscore2percentage($z, $m, $s)
{

    $lowerBound    = $z - 0.5; // continuity correction 0.5
    $upperBound    = $z + 0.5; // continuity correction 0.5

    $lowerBound_z  = ($lowerBound - $m + 0.5) / $s; // z score lower bound
    $upperBound_z  = ($upperBound - $m + 0.5) / $s; // z score upper bound
    $cdf_1         = cdf($lowerBound_z);
    $cdf_2         = cdf($upperBound_z);
    $cdf_3         = $cdf_2 - $cdf_1; // result percentage

    return $cdf_3;
}

// ||||||||||||||||||||||||||||||||||||||||||||||||||
// calculation
// ||||||||||||||||||||||||||||||||||||||||||||||||||

$q = 1 - $p; // against probability
$m = $n * $p; // expectation 'My'
$s = sqrt($n * $p * $q); // standard deviation 'Sigma'
$b = binomdist($n, $k, $p, $cumulative = false);  // binomdist

$ms_left_1 = $m - 1 * $s; // 34.00%                  = 34.00% to the left | unlucky
$ms_left_2 = $m - 2 * $s; // 34.00% + 13.50%         = 47.50% to the left | very unlucky
$ms_left_3 = $m - 3 * $s; // 34.00% + 13.50% + 2.35% = 49.85% to the left | extremely unlucky

$ms_right_1 = $m - 1 * $s; // 34.00%                  = 34.00% to the right | lucky
$ms_right_2 = $m - 2 * $s; // 34.00% + 13.50%         = 47.50% to the right | very lucky
$ms_right_3 = $m - 3 * $s; // 34.00% + 13.50% + 2.35% = 49.85% to the right | extremely lucky

if ($k <> $m)
{
    // worse or better
    $lowerBound    = 0; // continuity correction 0.5
    $upperBound    = $k + 0.5; // continuity correction 0.5
    $lowerBound_z  = ($lowerBound - $m + 0.5) / $s; // z score lower bound
    $upperBound_z  = ($upperBound - $m + 0.5) / $s; // z score upper bound
    $cdf_1         = cdf($lowerBound_z);
    $cdf_2         = cdf($upperBound_z);
    $cdf_3         = $cdf_1 - $cdf_2; // result percentage
    $worse         = 100 * abs($cdf_3); // hit this or worse
    $better        = 100 * (1 - (abs($cdf_3))); // hit better
}
else
{
    // worse or better
    $cdf_3         = roundUp(0.5, 4); // result percentage
    $worse         = 100 * abs($cdf_3); // hit this or worse
    $better        = 100 * (1 - (abs($cdf_3))); // hit better
}


if ((is_infinite ($b) or is_nan ($b)) and $s <= 3)
{
    echo '<hr>';
    echo 'approximate normal distribution not possible';
    echo '<hr>';
    echo 'binomdist not possible';
    echo '<hr>';
}
elseif ((is_infinite ($b) or is_nan ($b)) and $s > 3){

    // ||||||||||||||||||||||||||||||||||||||||||||||||||
    // approximate normal distribution
    // ||||||||||||||||||||||||||||||||||||||||||||||||||

    echo '<hr>';
    echo 'algorithm: approximate normal distribution';
    echo '<hr>';
    echo 'independet trials:' . $n;
    echo '<br>';
    echo 'successes: ' . $k;
    echo '<br>';
    echo 'equity for success: ' . $p;
    echo '<br>';
    echo 'My: ' . $m;
    echo '<br>';
    echo 'Sigma: ' . $s;
    echo '<hr>';
    echo 'normal: ' . '(x ' . $m . ' | y ' . zscore2percentage($m, $m, $s) . ')';
    echo '<br>';
    echo '<br>';
    echo 'unlucky: ' . '(x ' . $ms_left_1 . ' | y ' . zscore2percentage($ms_left_1, $m, $s) . ')';
    echo '<br>';
    echo 'very unlucky: ' . '(x ' . $ms_left_2 . ' | y ' . zscore2percentage($ms_left_2, $m, $s) . ')';
    echo '<br>';
    echo 'extremely unlucky: ' . '(x ' . $ms_left_3 . ' | y ' . zscore2percentage($ms_left_3, $m, $s) . ')';
    echo '<br>';
    echo '<br>';
    echo 'lucky: ' . '(x ' . $ms_right_1 . ' | y ' . zscore2percentage($ms_right_1, $m, $s) . ')';
    echo '<br>';
    echo 'very lucky: ' . '(x ' . $ms_right_2 . ' | y ' . zscore2percentage($ms_right_2, $m, $s) . ')';
    echo '<br>';
    echo 'extremely lucky: ' . '(x ' . $ms_right_3 . ' | y ' . zscore2percentage($ms_right_3, $m, $s) . ')';
    echo '<br>';
    echo '<br>';
    echo 'percentage of how many hit ' . $k . ' successe within ' . $n . ' trials or worse: ' . roundUp($worse, 3) . ' %';
    echo '<br>';
    echo 'percentage of how many did better: ' . roundDown($better, 3) . ' %';
    echo '<br>';
}
else
{

    // ||||||||||||||||||||||||||||||||||||||||||||||||||
    // binomdist
    // ||||||||||||||||||||||||||||||||||||||||||||||||||

    echo '<hr>';
    echo 'algorithm: binomdist';
    echo '<hr>';
    echo 'independet trials:' . $n;
    echo '<br>';
    echo 'successes: ' . $k;
    echo '<br>';
    echo 'equity for success: ' . $p;
    echo '<br>';
    echo 'My: ' . $m;
    echo '<br>';
    echo 'Sigma: ' . $s;
    echo '<hr>';
    echo 'normal: ' . '(x ' . $m . ' | y ' . binomdist($n, $ms_left_1, $p, $cumulative = false) . ')';
    echo '<br>';
    echo '<br>';
    echo 'unlucky: ' . '(x ' . $ms_left_1 . ' | y ' . binomdist($n, $ms_left_1, $p, $cumulative = false) . ')';
    echo '<br>';
    echo 'very unlucky: ' . '(x ' . $ms_left_2 . ' | y ' . binomdist($n, $ms_left_2, $p, $cumulative = false) . ')';
    echo '<br>';
    echo 'extremely unlucky: ' . '(x ' . $ms_left_3 . ' | y ' . binomdist($n, $ms_left_3, $p, $cumulative = false) . ')';
    echo '<br>';
    echo '<br>';
    echo 'lucky: ' . '(x ' . $ms_right_1 . ' | y ' . binomdist($n, $ms_right_1, $p, $cumulative = false) . ')';
    echo '<br>';
    echo 'very lucky: ' . '(x ' . $ms_right_2 . ' | y ' . binomdist($n, $ms_right_2, $p, $cumulative = false) . ')';
    echo '<br>';
    echo 'extremely lucky: ' . '(x ' . $ms_right_3 . ' | y ' . binomdist($n, $ms_right_3, $p, $cumulative = false) . ')';
    echo '<br>';
    echo '<br>';
    echo 'percentage of how many hit ' . $k . ' successe within ' . $n . ' trials or worse: ' . roundUp($worse, 3) . ' %';
    echo '<br>';
    echo 'percentage of how many did better: ' . roundDown($better, 3) . ' %';
    echo '<br>';
} 

[Unlikus]

Die ganzen Echos sind nicht so besonders schön, versuchs mit Templates.

[knuffiwuffi]

Ist nur ein test Script...baue in laravel 5 blade und andere Projekte in Drupal 8...

[jspit]

Hi,
in diesen Wiki-Artikel wird auf die Bedingungen eingegangen, für welche eine Normalverteilung als Näherung genommen werden kann. Daraus lassen sich leicht Fälle ableiten, wo mit den Verfahren größere Abweichungen auftreten können (n=30, k=14, p=0.4).
Vom Grundsatz ist es m.E. der richtige Weg, für kleine n die bewährte Formel zu nehmen und für große n die Näherung.
Die Frage ist auch, wie genau wird die Wahrscheinlichkeit benötigt bzw. wie sehen reale Fälle aus, wenn per Experiment Münzen geworfen werden oder Kugeln aus einem Topf gezogen werden?
Die kann relativ einfach simuliert werden. PHP ist für solche Aufgaben nicht besonders geeignet, aber es geht.
Der folgende Code benutzt eine solche Simulation, um durch mehrfache Wiederholung des Experimentes eine Wahrscheinlichkeit zu ermitteln:

function binomdistRnd($n, $k, $p, $cumulative = false){
  $zm = 800 + 5000/$n;
  $w = 0;
  for($j=0;$j<$zm;$j++){
    $t = 0;
    for($i=0;$i<$n;$i++) {
      $t += (int)(mt_rand(0,32767)/32768 < $p);
    }
    if(($t < $k and $cumulative) or ($t == $k)) $w++;
  }
  return $w/$zm;
}

$n = 30;
$k = 14;
$p = 0.4;

$w = binomdistRnd($n,$k,$p, true);
var_dump($w); 

Ein zufälliges Ergebnis:
float(0.81931034482759)

Der Code kann recht gut im rahmen eigener Tests für Vergleichszwecke herangezogen werden. Nachteilig ist, die Funktion liefert keine reproduzierbaren Ergebnisse und braucht viel Zeit.
Für den produktiven Einsatz daher nicht geeignet.
Hab die Ergebnisse für wenige Fälle mit den Excelberechnungen verglichen und Abweichungen von 1-2% gesehen.

LG jspit

[knuffiwuffi]

Hallo,

danke für das Feedback...ich habe meine Ergebnisse auch mit GeoGebra überprüft und für große n (1500 bis 13000) Abweichungen von +-2% festgestellt...als Näherung passt das aber