كيف يمكنني حساب فعال وظيفة التوزيع التراكمي بين الحدين؟

Question

دعنا نقول أنني أعرف احتمال "النجاح" هو P. أعرض اختبار N مرات، وأرى نجاحات S. الاختبار أقرب إلى قذف عملة عملة مرجحة بشكل غير متساو (ربما يكون الرؤوس ناجحة، ذيول فشلا).

أريد أن أعرف الاحتمال التقريبي لرؤية أي نجاحات، أو عدد من النجاحات أقل احتمالا من النجاحات.

لذلك على سبيل المثال، إذا كانت P 0.3، N 100، أحصل على 20 نجاحا، أبحث عن احتمال الحصول على 20 أو أقل النجاحات.

إذا كان الأمر كذلك، فهو الآخر، ف هو 0.3، N هو 100، وأحصل على 40 نجاحا، أبحث عن احتمال الحصول على 40 نجاحاتنا.

إنني أدرك أن هذه المشكلة تتعلق بإيجاد المنطقة تحت منحنى ذو حدية، ومع ذلك:

1. My Math-Fu لا يصل إلى مهمة ترجمة هذه المعرفة إلى رمز فعال
2. بينما أفهم منحنى ذو حدين من شأنه أن يمنح نتيجة دقيقة، أحصل على الانطباع بأنه سيكون غير فعال بطبيعته. طريقة سريعة لحساب النتيجة التقريبية تكفي.

يجب أن أشدد على أن هذا الحساب يجب أن يكون سريعا، ويجب أن يكون من الناحية المثالية محاسبا بحساب النقطة العائمة القياسية 64 أو 128 بت.

أبحث عن وظيفة تأخذ p و s و n - وإرجاع احتمال. نظرا لأنني أكثر دراية بمبلغ من الترميز الرياضي، فستفضل أن تستخدم أي إجابات رمز أو رمز الزائفة.

Answer 1

توزيع ذو الحدين بالضبط

def factorial(n): 
    if n < 2: return 1
    return reduce(lambda x, y: x*y, xrange(2, int(n)+1))

def prob(s, p, n):
    x = 1.0 - p

    a = n - s
    b = s + 1

    c = a + b - 1

    prob = 0.0

    for j in xrange(a, c + 1):
        prob += factorial(c) / (factorial(j)*factorial(c-j)) \
                * x**j * (1 - x)**(c-j)

    return prob

>>> prob(20, 0.3, 100)
0.016462853241869437

>>> 1-prob(40-1, 0.3, 100)
0.020988576003924564

تقدير طبيعي، جيد ل كبيرة ن

import math
def erf(z):
        t = 1.0 / (1.0 + 0.5 * abs(z))
        # use Horner's method
        ans = 1 - t * math.exp( -z*z -  1.26551223 +
                                                t * ( 1.00002368 +
                                                t * ( 0.37409196 + 
                                                t * ( 0.09678418 + 
                                                t * (-0.18628806 + 
                                                t * ( 0.27886807 + 
                                                t * (-1.13520398 + 
                                                t * ( 1.48851587 + 
                                                t * (-0.82215223 + 
                                                t * ( 0.17087277))))))))))
        if z >= 0.0:
                return ans
        else:
                return -ans

def normal_estimate(s, p, n):
    u = n * p
    o = (u * (1-p)) ** 0.5

    return 0.5 * (1 + erf((s-u)/(o*2**0.5)))

>>> normal_estimate(20, 0.3, 100)
0.014548164531920815

>>> 1-normal_estimate(40-1, 0.3, 100)
0.024767304545069813

تقدير Poisson: جيد ل كبير N والصغير ص

import math

def poisson(s,p,n):
    L = n*p

    sum = 0
    for i in xrange(0, s+1):
        sum += L**i/factorial(i)

    return sum*math.e**(-L)

>>> poisson(20, 0.3, 100)
0.013411150012837811
>>> 1-poisson(40-1, 0.3, 100)
0.046253037645840323

Answer 2

كنت في مشروع، حيث كنا بحاجة إلى أن نكون قادرين على حساب CDF ذو الحدين في بيئة لم يكن لديك وظيفة فنية أو جاما محددة. استغرق الأمر مني بضعة أسابيع، لكنني انتهى بي الأمر إلى الخوارزمية التالية التي بحسقت CDF بالضبط (أي أي تقريب ضروري). بيثون هو أساسا جيدة مثل pseudocode، أليس كذلك؟

import numpy as np

def binomial_cdf(x,n,p):
    cdf = 0
    b = 0
    for k in range(x+1):
        if k > 0:
            b += + np.log(n-k+1) - np.log(k) 
        log_pmf_k = b + k * np.log(p) + (n-k) * np.log(1-p)
        cdf += np.exp(log_pmf_k)
    return cdf

جداول الأداء مع x. بالنسبة للقيم الصغيرة من X، فإن هذا الحل حول ترتيب من حيث الحجم أسرع من scipy.stats.binom.cdf, مع أداء مماثل في حوالي X = 10000.

لن أذهب إلى اشتقاق كامل لهذه الخوارزمية لأن Stackoverflow لا يدعم MathJax، لكن توجيهه هو أولا تحديد التكافؤ التالي:

• للجميع K> 0، sp.misc.comb(n,k) == np.prod([(n-k+1)/k for k in range(1,k+1)])

الذي يمكننا إعادة كتابة باسم:

• sp.misc.comb(n,k) == sp.misc.comb(n,k-1) * (n-k+1)/k

أو في مساحة السجل:

• np.log( sp.misc.comb(n,k) ) == np.log(sp.misc.comb(n,k-1)) + np.log(n-k+1) - np.log(k)

نظرا لأن CDF هو جمع PMF، يمكننا استخدام هذه الصيغة لحساب معامل الحدين (سجلها b في الوظيفة أعلاه) ل PMF_ = {x = I} من المعامل المحسوب ل PMF_ {x = i-1}. هذا يعني أننا نستطيع أن نفعل كل شيء داخل حلقة واحدة باستخدام البكرات، ولا نحتاج إلى حساب أي فصائل!

السبب في أن معظم الحسابات تتم في مساحة السجل هي تحسين الاستقرار العددي للمصطلحات المتعددة، أي p^x و (1-p)^(1-x) لديك القدرة على أن تكون كبيرة للغاية أو صغيرة للغاية، والتي يمكن أن تسبب أخطاء حسابية.

تحرير: هل هذه خوارزمية رواية؟ لقد تمتد حولها وإيقافها منذ ذلك الحين قبل أن أرسلت هذا، وأنا أتساءل بشكل متزايد إذا كان يجب أن أكتب هذا أكثر رسميا وتقديمه إلى مجلة.

Answer 3

أعتقد أنك تريد تقييم وظيفة بيتا غير كاملة.

هناك تطبيق لطيف باستخدام تمثيل جزء مستمر في "الوصفات العددية في C"، الفصل 6: "الوظائف الخاصة".

Answer 4

لا أستطيع أن أضمر تماما الكفاءة، لكن SCIPY لديه وحدة لهذا

from scipy.stats.distributions import binom
binom.cdf(successes, attempts, chance_of_success_per_attempt)

Answer 5

فعالة، والأهم من ذلك، توجد خوارزمية مستقرة رقمية في مجال منحنيات bezier. المستخدمة في تصميم بمساعدة الكمبيوتر. تسمى خوارزمية دي كاستيلجاو تستخدم لتقييم بيرنشتاين متعدد الحدود تستخدم لتحديد منحنيات bezier.

أعتقد أنني سمحت فقط رابط واحد لكل إجابة لذلك ابدأ ويكيبيديا - بولنشتاين متعدد الحدود

لاحظ العلاقة الوثيقة للغاية بين التوزيع العيادي وليستنشتاين متعدد الحدود. ثم انقر عبر الرابط في خوارزمية دي كاستيلجو.

دعنا أقول إنني أعرف احتمال رمي رؤوس مع عملة معينة هو P. ما هو احتمال رمي Times Coin T Times واحصل على رؤوس S على الأقل؟

• تعيين ن = ر
• تعيين بيتا [i] = 0 ل i = 0، ... S - 1
• تعيين بيتا [i] = 1 ل i = s، ...
• تعيين ر = ص
• تقييم B (ر) باستخدام de casteljau

أو في معظم رؤساء؟

• تعيين ن = ر
• تعيين بيتا [i] = 1 ل i = 0، ... s
• ضبط بيتا [i] = 0 ل i = s + 1، ...
• تعيين ر = ص
• تقييم B (ر) باستخدام de casteljau

الكود المصدر المفتوح ربما موجود بالفعل. منحنيات نائب (منحنيات غير موحدة عقلانية B-SPLINE) تعميم منحنيات Bezier وتستخدم على نطاق واسع في CAD. جرب opennurbs (الترخيص ليبرالي للغاية) أو الفشل في أن تتكلم المفتوحة (رخصة أقل ليبرالية وأفتانية إلى حد ما). كلا مجموعة الأدوات في C ++، على الرغم من وجود ارتباطات IIRC،.

Answer 6

إذا كنت تستخدم Python، فلا حاجة لرمزها بنفسك. SCIPY حصلت على مغطاة:

from scipy.stats import binom
# probability that you get 20 or less successes out of 100, when p=0.3
binom.cdf(20, 100, 0.3)
>>> 0.016462853241869434

# probability that you get exactly 20 successes out of 100, when p=0.3
binom.pmf(20, 100, 0.3)
>>> 0.0075756449257260777

Answer 7

من جزء سؤالك "الحصول على رؤوس الأقوال على الأقل" تريد وظيفة توزيع الحد الثناني المتزامور. يرى http://en.wikipedia.org/wiki/binomial_distribution. بالنسبة للمعادلة، والذي يوصف بأنه من حيث "وظيفة النسخة التجريبية غير المكتملة النظامية" (كما أجاب بالفعل). إذا كنت ترغب فقط في حساب الإجابة دون الحاجة إلى تنفيذ الحل بأكمله، توفر مكتبة جنو العلمية الوظيفة: GSL_CDF_BINOMIAL_P و GSL_CDF_BINOMIAL_Q.

Answer 8

ال مشروع DCDFLIB يحتوي على وظائف C # (مغلفة حول رمز C) لتقييم العديد من وظائف CDF، بما في ذلك التوزيع ذو الحدين. يمكنك العثور على رمز C و FORTRAN الأصلي هنا. وبعد هذا الرمز اختبار جيدا ودقيقة.

إذا كنت ترغب في كتابة التعليمات البرمجية الخاصة بك لتجنب التعلق على مكتبة خارجية، فيمكنك استخدام التقريب الطبيعي إلى الحدين المذكورين في إجابات أخرى. هنا بعض الملاحظات على كيف جيدة التقريب هو تحت ظروف مختلفة. إذا ذهبت هذا الطريق وتحتاج إلى رمز لحساب CDF العادي، إليك كود بيثون للقيام بذلك. إنه فقط حوالي عشرات خطوط التعليمات البرمجية ويمكن بسهولة الوصول إلى أي لغة أخرى. ولكن إذا كنت ترغب في دقة عالية وكود فعال، فأنت أفضل حالا باستخدام رمز الطرف الثالث مثل DCDFLIB. ذهب العديد من سنوات الرجل إلى إنتاج تلك المكتبة.

Answer 9

يحاول هذا, ، المستخدمة في GMP. مرجع آخر هو هذه.

Answer 10

import numpy as np
np.random.seed(1)
x=np.random.binomial(20,0.6,10000) #20 flips of coin,probability of 
                                 heads percentage and 10000 times 
                                  done.
sum(x>12)/len(x)

The output is 41% of times we got 12 heads.