Probability of something happening at least n out of m times

11

2

Write a program or function, that given a success probability p, a number n and a number of trials m returns the chance of at least n successes out of m trials.

Your answer must be precise to at least 5 digits after the decimal.

Test cases:

 0.1, 10, 100 -> 0.54871
 0.2, 10, 100 -> 0.99767
 0.5, 13,  20 -> 0.13159
 0.5,  4,   4 -> 0.06250
0.45, 50, 100 -> 0.18273
 0.4, 50, 100 -> 0.02710
   1,  1,   2 -> 1.00000
   1,  2,   1 -> 0.00000
   0,  0,   1 -> 1.00000
   0,  0,   0 -> 1.00000
   0,  1,   1 -> 0.00000
   1,  1,   0 -> 0.00000

orlp

Posted 2016-05-03T10:17:14.277

Reputation: 37 067

3Would you care to include a formula to those of us who haven't studied binomial distribution? – Leaky Nun – 2016-05-03T11:09:12.963

2@KennyLau Sorry, that is part of the challenge. – orlp – 2016-05-03T11:15:42.213

Answers

3

Jelly, 15 14 bytes

2ṗ’S<¥ÐḟCạ⁵P€S

Reads m, n and p (in that order) as command-line arguments. Try it online!

Note that this approach requires O(2m) time and memory, so it isn't quite efficient enough for the test cases where m = 100. On my machine, the test case (m, n, p) = (20, 13, 0.5) takes roughly 100 seconds. It requires too much memory for the online interpreter.

How it works

2ṗ              Cartesian product; yield all vectors of {1, 2}^n.
  ’             Decrement, yielding all vectors of {0, 1}^n.
      Ðḟ        Filter; keep elements for which the link to the left yields False.
     ¥          Combine the two links to the left into a dyadic chain.
   S              Sum, counting the number of ones.
    <             Compare the count with n. 
        C       Complement; map z to 1 - z.
         ạ⁵     Compute the absolute difference with p.
           P€   Compute the product of each list.
             S  Compute the sum of all products.

Dennis

Posted 2016-05-03T10:17:14.277

Reputation: 196 637

9

Mathematica, 29 bytes

BetaRegularized[#3,#,1+#2-#]&

Takes input in the order n,m,p. Mathematica is so good, it even golfs your code for you:

enter image description here

BetaRegularized is the regularised incomplete beta function.

Sp3000

Posted 2016-05-03T10:17:14.277

Reputation: 58 729

6

R, 32 31 bytes

function(p,n,m)pbeta(p,m,1+n-m)

edit - 1 byte switching to beta distribution (along the lines of @Sp3000 Mathematica Answer)

mnel

Posted 2016-05-03T10:17:14.277

Reputation: 826

3

Python, 57 bytes

f=lambda p,n,m:m and(1-p)*f(p,n,m-1)+p*f(p,n-1,m-1)or n<1

The recursive formula for binomial coefficients, except the base case m==0 indicates whether the remaining number of required successes n is nonnegative, with True/False for 1/0. Because of its exponential recursion tree, this stalls on large inputs.

xnor

Posted 2016-05-03T10:17:14.277

Reputation: 115 687

To test out this answer for large cases, add caching using from functools import lru_cache; f = lru_cache(None)(f). – orlp – 2016-05-03T12:15:51.447

@orlp Thanks, I confirmed the large test cases. – xnor – 2016-05-03T12:23:13.160

3

Haskell, 73 bytes

g x=product[1..x];f p n m=sum[g m/g k/g(m-k)*p**k*(1-p)**(m-k)|k<-[n..m]]

Damien

Posted 2016-05-03T10:17:14.277

Reputation: 2 407

3

MATLAB, 78 71 bytes

Saved 7 bytes thanks to Luis Mendo!

@(m,k,p)sum(arrayfun(@(t)prod((1:m)./[1:t 1:m-t])*p^t*(1-p)^(m-t),k:m))

ans(100,10,0.1)
0.5487

The arrayfun function is no fun, but I haven't found a way to get rid of it...

Stewie Griffin

Posted 2016-05-03T10:17:14.277

Reputation: 43 471

1

Pyth, 26 bytes

AQJEsm**.cHd^Jd^-1J-HdrGhH

Try it online!

Uses standard cumulative binomial distribution.

Leaky Nun

Posted 2016-05-03T10:17:14.277

Reputation: 45 011

1

Pyth, 20 bytes

JEKEcsmgsm<O0QKJCGCG

Try it online!

Note: CG is a very large number which the interpreter cannot handle. Therefore, the number of trials have been lowered to ^T3 which is one thousand. Therefore, the link produces an inaccurate result.

Uses pure probabilistic approach.

Leaky Nun

Posted 2016-05-03T10:17:14.277

Reputation: 45 011

I don't think a probabilistic approach would be valid for this question, but we'd have to ask @orlp – Sp3000 – 2016-05-03T11:21:06.897

You need on the order of 1/c^2 trials to get within accuracy c with high probability, so that would be ~10^10 for five decimal places. – xnor – 2016-05-03T11:24:13.847

CG is a very large number. In fact, it is the string "abc...z" converted from base-256 to decimal. – Leaky Nun – 2016-05-03T11:26:59.557

^TT would be 10^10, and ^T11 would be 10^11. – Leaky Nun – 2016-05-03T11:31:24.777

@Sp3000 I am not sure if that would ping @orlp – Leaky Nun – 2016-05-03T11:34:15.417

@Sp3000 If you can guarantee 5 decimals of precision somehow through probabilistic means, go ahead. – orlp – 2016-05-03T11:34:21.200

@KennyLau It didn't, but I'm lurking. – orlp – 2016-05-03T11:34:30.683

I can guarantee 5 decimal places. – Leaky Nun – 2016-05-03T11:38:19.743

2If "probabilstic" means random, you can't guarante an accurate value, no matter how many realizations you average. In fact, the result is different every time. – Luis Mendo – 2016-05-03T14:29:46.223

It is highly improbable that it is *not* accurate to 5 decimal places. – Leaky Nun – 2016-05-03T14:34:43.493

2There is always a nonzero probability that the result is not accurate to 5 decimal places. Therefore it doesn't fulfill the requirement Your answer must be precise to at least 5 digits – Luis Mendo – 2016-05-03T14:38:00.950

In science, *highly improbable* is equivalent to impossible. – Leaky Nun – 2016-05-03T14:42:41.930

But this is not science, this is maths :-) – Luis Mendo – 2016-05-03T15:15:41.473

Let the poster judge. – Leaky Nun – 2016-05-03T15:18:49.927

1

JavaScript (ES7), 82 bytes

(p,n,m)=>[...Array(++m)].reduce((r,_,i)=>r+(b=!i||b*m/i)*p**i*(1-p)**--m*(i>=n),0)

Saved 1 byte by using reduce! Explanation:

(p,n,m)=>               Parameters
 [...Array(++m)].       m+1 terms
  reduce((r,_,i)=>r+    Sum
   (b=!i||b*m/i)*       Binomial coefficient
   p**i*(1-p)**--m*     Probability
   (i>=n),              Ignore first n terms
   0)

Neil

Posted 2016-05-03T10:17:14.277

Reputation: 95 035

1

MATL, 23 bytes

y2$:Xni5M^1IG-lG7M-^**s

Inputs are in the order m, n, p.

Try it online!

This does a direct computation summing the terms from n to m of the binomial probability (mass) function.

Luis Mendo

Posted 2016-05-03T10:17:14.277

Reputation: 87 464

1

Octave, 26 bytes

@(p,n,m)1-binocdf(n-1,m,p)

This is an anonymous function. To use it, assign it to a variable.

Try it here.

Luis Mendo

Posted 2016-05-03T10:17:14.277

Reputation: 87 464

1

Jelly, 18 17 bytes

⁵C*ạ×⁵*¥×c@
R’çSC

Reads n, m and p (in that order) as command-line arguments. Try it online!

Dennis

Posted 2016-05-03T10:17:14.277

Reputation: 196 637

0

TI-Basic, 17 bytes

Precise to 10 decimals, could be adjusted anywhere from 0-14 decimals with more code.

Prompt P,N,M:1-binomcdf(M,P,N-1

Timtech

Posted 2016-05-03T10:17:14.277

Reputation: 12 038

0

Haskell, 54 bytes

(p%n)m|m<1=sum[1|n<1]|d<-m-1=(1-p)*(p%n)d+p*(p%(n-1))d

Defines a function (%). Call it like (%) 0.4 2 3.

xnor

Posted 2016-05-03T10:17:14.277

Reputation: 115 687

n<1 instead of n<=0. – Damien – 2016-05-04T06:49:48.450

0

Mathematica, 48 bytes

Sum[s^k(1-s)^(#3-k)#3~Binomial~k,{k,##2}]/.s->#&

Uses the binomial distribution probability formula to calculate the chance of k successes for k from n to m. Handles the edge cases by using a symbolic sum where s is a symbolic variable for the probability that is later replaced with the actual value p. (Since s0 = 1 but 00 is indeterminate.)

Example

miles

Posted 2016-05-03T10:17:14.277

Reputation: 15 654