14
\$\begingroup\$

Given n, k, and p, find the probability that a weighted coin with probability p of heads will flip heads at least k times in a row in n flips, correct to 3 decimal digits after decimal point (changed from 8 because I don't know how to estimate the accuracy of double computation for this (was using arbitrary precision floating points before)).

Edits:

Examples:

n, k, p     -> Output
10, 3, 0.5  -> 0.5078125
100, 6, 0.6 -> 0.87262307
100, 9, 0.2 -> 0.00003779
\$\endgroup\$
8
  • 1
    \$\begingroup\$ That's a terse description for sure, but the task is actually pretty clear IMO. \$\endgroup\$
    – Arnauld
    Jun 30 at 14:26
  • 5
    \$\begingroup\$ @Arnauld I haven't voted to close, for the same reason as you mentioned: the challenge is clear. It's just a bit lacking in challenge rules, examples, and especially test cases. \$\endgroup\$ Jun 30 at 14:29
  • 2
    \$\begingroup\$ Are we guaranteed that k isn't greater than n? \$\endgroup\$
    – xnor
    Jun 30 at 15:21
  • 7
    \$\begingroup\$ Not a bad first challenge, though I'd recommend checking out the sandbox in the future, especially for more complex challenges. It's not a perfect system, but it'll help you catch small mistakes / ambiguities before posting to main. Welcome to the site :-) \$\endgroup\$ Jun 30 at 15:43
  • 2
    \$\begingroup\$ @dancxviii Not by default per se. \$n^0=1\$ and \$0^n=0\$, so \$0^0\$ conflicts and is usually undefined. But if we would assume \$0^0=1\$, a lot of things fall into place nicely in mathematics: see this Wikipedia page for details. \$\endgroup\$ Jul 1 at 7:04

15 Answers 15

10
\$\begingroup\$

05AB1E, 40 39 18 17 bytes

Brute-force approach:

0Lα1ݲãʒγOà³@}èPO

Byte-count more than halved by porting @JonathanAllan's Jelly answer, which uses a similar approach as @DominicVanEssen's top R answer, so make sure to upvote them as well!

Inputs in the order \$p,n,k\$.
Very slow, so use small \$n\$ in order to run it on TIO.

Try it online.

Original mathematical 40 39 bytes approach:

÷L®s©>m²®¹*-D>®/$-*I+s®<cI®¹*m$-®<m)øPO

Inputs in the order \$k,n,p\$.

Try it online.

General explanation of the larger original program:

Even though the challenge description was just a single sentence, it lead me on a chase down the rabbit hole. After quite a long search, I came across this MathExchange answer: Probability for the length of the longest run in n Bernoulli trials, for which the accepted answer contains the generalized formula for this challenge, originally solved by de Moivre in 1738:

$$\mathbb{P}(\ell_n \geq k)=\sum_{j=1}^{\lfloor n/k\rfloor} (-1)^{j+1}\left(p+\left({n-jk+1\over j}\right)(1-p)\right){n-jk\choose j-1}p^{jk}(1-p)^{j-1}$$

Where \$\mathbb{P}(\ell_n \geq k)\$ is the probability that a consecutive run of the same unfair coin result with a length \$\ell_n\$ is at least \$k\$ flips. \$p\$, \$k\$, and \$n\$ are as given by the challenge description, with \$p\$ being a decimal probability (e.g. 0.70 for 70%), and \$k\$ and \$n\$ being integers.

Code explanation:

0L              # Push pair [1,0]
  α             # Get the absolute difference with the first (implicit) input p:
                #  [1-p,p]
   1Ý           # Push pair [0,1]
     ²ã         # Get the second input n'th cartesian product
                # (creating all n-sized lists with 0s/1s)
       ʒ        # Filter this list by:
        γ       #  Group the 0s/1s into consecutive adjacent items
         O      #  Sum each group
          à     #  Pop and push the maximum
           ³@   #  Check if this largest group of 1s is >= the third input k
       }è       # After the filter: (0-based) index each 0/1 into pair [1-p,p]
         P      # Take the product of each inner list
          O     # Take the sum of all values
                # (after which the result is output implicitly)
÷               # Integer divide the first two (implicit) inputs: n//k
 L              # Pop and push a list in the range [1,n//k] (its values are j)
  ®             # Push -1
   s            # Swap so the list is at the top
    ©           # Store it in variable `®` (without popping)
     >          # Increase each j in the list by 1
      m         # Pop both, and calculate -1**(j+1) for each
  ®             # Push the list `®` containing j again
   ¹*           # Multiply each j by the first input k
 ²   -          # Subtract each from the second input n
      D         # Duplicate this n-jk list, since we need it again later
       >        # Increase each by 1
        ®/      # Divide each by j
          $     # Push 1 and the third input p
           -    # Subtract: 1-p
            *   # Multiply this to each (n-jk+1)/j
             I+ # And add the third input p to each
 s              # Swap so the n-jk list is at the top again
  ®             # Push list `®` containing j again
   <            # Decrease each j by 1
    c           # Calculate the binomial coefficients of n-jk and j-1
 I              # Push the third input p again
  ®             # Push list `®` containing j again
   ¹*           # Multiply each j by the first input k again
     m          # Pop both, and calculate p**(jk)
 $              # Push 1 and the third input p again
  -             # Subtract again: 1-p
   ®            # Push list `®` containing j again
    <           # Decrease each j by 1 again
     m          # Calculate (1-p)**(j-1)
 )              # Wrap these five lists on the stack into a list of lists
  ø             # Zip/transpose; swapping rows/columns
   P            # Get the product of each inner list
O               # And finally sum everything together
                # (which is output implicitly as result)
\$\endgroup\$
5
  • 1
    \$\begingroup\$ Did not know about this formula, was using the markov method and was looking for a trick related to decomposition, this formula is cleaner \$\endgroup\$
    – dancxviii
    Jun 30 at 16:17
  • \$\begingroup\$ It's funny how golf-languages seem to often prefer to port an approach from other golf-languages... the 'add the probabilities of all possible outcomes' approach (with explanation) was posted in my R answer substantially before the Jelly implementation of the same approach... \$\endgroup\$ Jul 1 at 7:55
  • \$\begingroup\$ @DominicvanEssen To me personally, Jelly with Jonathan's explanation is easy to read, and therefore to port to 05AB1E's builtins. I personally have never programmed in R, and although I can read some of it, I can't intuitively understand most of it. Also, reading your explanation, your approach is slightly different because R couldn't hold big enough matrices? The Jelly and this 05AB1E answer generate all possible outcomes once, whereas your R answer seems to loop 1e19 times, and check if at least \$k\$ unfair coin heads are present in a random result. Both are brute-force, but different. \$\endgroup\$ Jul 1 at 8:08
  • 1
    \$\begingroup\$ Yes - of course I get it that Jelly (& the explanation) is likely to be easier to port to O5AB1E's builtins! Still, the approach is identical to the 116-byte 'sum the probabilities of every possible outcome of n coin flips' R version, although the more-detailed explanation for that is only in the code comments. Possibly I should get better at explaining... \$\endgroup\$ Jul 1 at 9:17
  • \$\begingroup\$ @DominicvanEssen Ah wait, I was looking at your 92 bytes function instead of 116 bytes.. I will edit my answer to also credit you. \$\endgroup\$ Jul 1 at 9:34
7
\$\begingroup\$

Python, 66 bytes (-1 @xnor)

f=lambda n,k,p:n>=k and(1-p*(n>k))*(1-f(n+~k,k,p))*p**k+f(n-1,k,p)

Attempt This Online!

Adapted from my answer to a similar challenge. Note that I add memoization in the footer. This in theory doesn't change the result but it greatly accelerates the recursion.

How?

Uses the recurrence

\$f(n,k,p)=f(n-1,k,p)+p^k(1-p)(1-f(n-k-1,k,p))\$

(valid for n>k) which is obtained by accounting for words that have k consecutive heads somewhere in the first n-1 tosses, words that end in a tail followed by k heads and the overlap of these two groups.

\$\endgroup\$
1
  • 1
    \$\begingroup\$ It looks like you can save the space after and by rearranging terms \$\endgroup\$
    – xnor
    Jun 30 at 22:45
7
\$\begingroup\$

Jelly,  19  15 bytes

Ø.ṗṣ1ZṫɗƇ_⁵AP€S

A full program accepting n, k, and p that prints the result.

Try it online! (n=100 is too big for such an inefficient program.)

How?

Creates all \$2^{n}\$ possible outcomes, filters them to those containing a run of at least \$k\$ heads and then sums the probability of each occurrence.

The probability of a given occurrence is the product of p's and (1-p)s identified by heads and tails respectively. For example the chance of [tails, tails, heads, heads, heads] is the product of [(1-p), (1-p), p, p, p] i.e. \$p^3(1-p)^{2}\$.

Ø.ṗṣ1ZṫɗƇ_⁵AP€S - Main Link: integer n, integer k (p is accessed later)
Ø.              - [0,1]
  ṗ             - ([0,1]) Cartesian power (n)
                  -> all length-n tosses with 0 as heads and 1 as tails
        Ƈ       - filter keep those for which:
       ɗ        -   last three links as a dyad - f(Outcome, k):
   ṣ1           -     split at 1s -> runs of heads
     Z          -     transpose
      ṫ         -     tail from index k -> empty (falsey) if longest run < k
          ⁵     - third program argument, p
         _      - subtract -> valid outcomes with heads: -p; tails: 1-p
           A    - absolute values -> valid outcomes with heads: p; tails 1-p
            P€  - product of each -> valid outcome probabilities
              S - sum -> total probability of any valid outcome

If we could take the probability of tails (\$1-p\$) instead of heads (\$p\$) 14 byes taking 1-p, n, k:

C©Ƭṗṣ®ZṫɗƇ⁵P€S

Try it online!

\$\endgroup\$
5
\$\begingroup\$

R, 129 122 116 bytes

Edit: -2 bytes thanks to pajonk (which led to -4 more...)

function(n,k,p)sum(apply(which(array(T,rep(2,n)),T)-1,1,function(v,r=rle(v))prod(p*v+(1-p)*!v)*(max(r$l[!r$v])>=k)))

Try it online!

Exact solution: calculates the probability of every possible outcome of n coin flips, and sums those that contain at least k heads-in-a-row.

How?

function(n,k,p){
    a=expand.grid(rep(list(1:0),n)) # a is a matrix of all possible outcomes of n flips
                                    # with heads represented as 0, tails as 1
    pvals=apply(a,1,function(v)prod(p*v+(1-p)*!v))
                                    # pvals are the probabilities of each outcome 
                                    # by multiplying the probability of each single flip
    itsarun=apply(a,1,function(v)max(rle(v)$l[!rle(v)$v])>=k)
                                    # itsarun is TRUE if the longest run of 0s in each row is >=k
                                    # (rle(v)$l = length of each run, rle(v)$v = value of each run)
    return(sum(pvals*itsarun))      # return the sum of all the pvals for successful outcomes
}

R, 98 bytes

function(n,k,p){while((T=T+1)<1e19)F=F+(max((r=rle(sample(1:0,n,T,c(1-p,p))))$l[!r$v])>=k);F/1e19}

Try it online with a low-accuracy version

Brute-force approach: performs 1e19 series of n flips, and counts the number of times we get k heads-in-a-row. With this number of repetitions, the answer has a very high chance of being accurate to 8 decimal places, but there is a low-but-finite chance that this won't be the case...

We need to do them one-after-the-other in a while loop (rather than the more idiomatic vectorized R approach of building a matrix and checking the rows) since the number of iterations needed to obtain 8 decimal places of accuracy is too big to be contained in an R matrix.
The test link substitutes 1e5 in place of 1e19, for a lower-accuracy output without timing-out. Feel free to run the high-accuracy version on your own computer if you have the time to wait.


R, 83 bytes

function(n,k,p,m=cbind(rbind(1-p,diag(k)*p),!k:0))Reduce(`%*%`,rep(list(m),n))[k+1]

Try it online!

Port of dancxviii's Markov method - upvote that one - but sadly still longer than pajonk's R answer...

\$\endgroup\$
6
  • 2
    \$\begingroup\$ I'm not sure an answer like this can be valid. It's only going to produce a valid answer most of the time, not all. \$\endgroup\$
    – pxeger
    Jun 30 at 15:30
  • \$\begingroup\$ @pxeger - I've added an exact solution, too. \$\endgroup\$ Jun 30 at 16:21
  • \$\begingroup\$ -2 bytes \$\endgroup\$
    – pajonk
    Jun 30 at 20:28
  • \$\begingroup\$ @pajonk - Thanks! The which(array(),T) trick instead of expand.grid is really nice! \$\endgroup\$ Jun 30 at 21:00
  • \$\begingroup\$ Port of loopy walt's python answer is just 70 bytes. \$\endgroup\$
    – pajonk
    Jul 1 at 5:49
4
\$\begingroup\$

Vyxal , 23 bytes

2$ÞẊ'2€vLG›¹>;⁰⌐⁰"$İvΠ∑

Try it Online!

Port of Jelly. Really slow.

Previous answer (much faster):

Vyxal , 41 bytes

ḭɾ:£›u$e¹¥□h*-:›¥/⁰⌐*⁰+$¥‹ƈ?¥?*e⁰⌐¥‹eWƒ*∑

Try it Online!

Port of 05AB1E.

\$\endgroup\$
4
\$\begingroup\$

PARI/GP, 53 bytes

k->p->g(n)=if(n>=k,p^k*(1-p*(n>k))*(1-g(n---k))+g(n))

Attempt This Online!

A port of loopy walt's Python answer. Takes input in the the form (k)(p)(n)


PARI/GP, 56 bytes

f(n,k,p)=Pol((1-y=p*x)/(t=1-x)/(t/y^k+x-y)+O(x^n*x))\x^n

Attempt This Online!

The generation function of the sequence (for given \$p, k\$) is \$\frac{p^k\ x^k\ (1-p\ x)}{(1-x)(1-x+(1-p)\ p^k\ x^{k+1})}\$.


This generation function can be derived from loopy walt's recurrence formula:

\$f(n,k,p)=f(n-1,k,p)+p^k(1-p)(1-f(n-k-1,k,p))\$.

This recurrence formula is valid when \$n>k\$. In addition, when \$n=k\$, we have \$f(n,k,p)=p^k\$ when \$n=k\$, and \$f(n,k,p)=0\$ when \$n<k\$.

Let \$F(k,p)\$ be the generation function (for given \$p, k\$), we have

\$\begin{align} F(k,p) &= \sum_{n=0}^{\infty}f(n,k,p)\ x^n \\ &= p^k\ x^k+\sum_{n=k+1}^{\infty}f(n,k,p)\ x^n \\ &= p^k\ x^k+\sum_{n=k+1}^{\infty}(f(n-1,k,p)+p^k(1-p)(1-f(n-k-1,k,p)))\ x^n \\ &= p^k\ x^k+\sum_{n=k+1}^{\infty}f(n-1,k,p)\ x^n+p^k(1-p)\sum_{n=k+1}^{\infty}(x^n-f(n-k-1,k,p)x^n) \\ &= p^k\ x^k+x\sum_{n=0}^{\infty}f(n,k,p)\ x^n+p^k(1-p)\sum_{n=0}^{\infty}(x^{n+k+1}-f(n,k,p)x^{n+k+1}) \\ &= p^k\ x^k+x\ F(k,p)+p^k(1-p)(\frac{x^{k+1}}{1-x}-x^{k+1}F(k,p)) \end{align}\$.

Solving this equation, we have \$F(k,p)=\frac{p^k\ x^k(1-p\ x)}{(1-x)(1-x+(1-p)\ p^k\ x^{k+1})}\$.

\$\endgroup\$
4
\$\begingroup\$

JavaScript (ES7), 106 bytes

Quickly computes the probability, using the formula found by Kevin Cruijssen.

(n,k,p,q=1-p)=>(g=s=>(x=n-(y=j*k))<0?0:s*(p-~x/j*q)*(h=k=>--k?h(k)*(x-k+1)/k:p**y*q**~-j++)(j)+g(-s))(j=1)

Try it online!


JavaScript (ES12), 58 bytes

Using loopy walt's recursive formula is much shorter and much slower. This version uses a cache to speed it up (-7 bytes without the cache).

Expects (k)(p)(n).

k=>p=>g=n=>g[n]||=n>=k&&p**k*(1-p*(n>k))*(1-g(--n-k))+g(n)
//         \_____/
//          cache

Attempt This Online!

\$\endgroup\$
3
\$\begingroup\$

PARI-GP, 75 bytes

Port of below solution (can run n ~ 1000000000, k ~ 150)

f(n,k,p)=(matrix(k++,k,a,b,(a-1==b)*p+(a<2)*(b<k)*(1-p)+(a+b==k+k))^n)[k,1]

Python (Numpy), 104 bytes

Noone has the markov method yet:

The states are the number of consecutive heads starting at state 0, with the final state as state k. State i (i < k) goes to state i+1 with probability p and state 0 with probability 1-p. State k always goes to state k. The matrix is:

\$ \begin{bmatrix} 1 - p & 1 - p & 1 - p & \dots & 1 - p & 0 \\ p & 0 & 0 & \dots & 0 & 0 \\ 0 & p & 0 & \dots & 0 & 0 \\ \dots & \dots & \dots & \dots & \dots & \dots \\ \dots & \dots & \dots & \dots & \dots & \dots \\ 0 & 0 & 0 & \dots & p & 1 \\ \end{bmatrix} \$

Taking this matrix to the nth power and multiplying it by the vector \$\left[1, 0, 0, \dots, 0 \right]^T \$ gives the answer.

This solution is not an clean mathematical sum but it is much more computationally efficient then the sum.

m=zeros((k+1,k+1))
for i in range(k):
    m[0][i]=1-p
    m[i+1][i]=p
m[k][k]=1
print(matrix_power(m,n)[k][0])
\$\endgroup\$
2
  • \$\begingroup\$ Welcome to Code Golf, and nice answer! \$\endgroup\$ Jul 1 at 2:03
  • 2
    \$\begingroup\$ Nice answer. But according to code-golf's rule, you need to write either a function or a full program that takes n,p,k as input. A snippet that assume n,p,k are predefined is not allowed. You also need to include from numpy import*. \$\endgroup\$
    – alephalpha
    Jul 1 at 4:04
2
\$\begingroup\$

Desmos, 80 bytes

P(p,k,n)=∑_{j=1}^{n/k}-(-1)^j(p+(n-jk+1)(1-p)/j)nCr(n-jk,j-1)p^{jk}(1-p)^{j-1}

Try it on Desmos or prettified.

Thanks to Kevin Cruijssen for finding this formula.

\$\endgroup\$
2
\$\begingroup\$

Factor + math.combinatorics math.unicode, 133 bytes

[| n k p | n k /i [1,b] [| j | -1 j 1 + ^ n j k * - :> z z 1 + j
/ 1 p - * p + z j 1 - nCk p j k * ^ * 1 p - j 1 - ^ * * * ] map Σ ]

Try it online!

Translation of the formula given in Kevin Cruijssen's 05AB1E answer.

\$\endgroup\$
2
  • 1
    \$\begingroup\$ How many bytes is it without the whitespace? :P \$\endgroup\$ Jun 30 at 19:16
  • 3
    \$\begingroup\$ @thejonymyster 72 bytes. I prefer programs that compile, though. \$\endgroup\$
    – chunes
    Jun 30 at 20:45
2
\$\begingroup\$

R, 70 bytes

`[`=function(k,n,p)`if`(n<k,0,(1-p*(n>k))*(1-k[n-k-1,p])*p^k+k[n-1,p])

Try it online!

Port of @loopy walt's Python answer.

For different approaches in R see @Dominic van Essen's answer.

\$\endgroup\$
2
\$\begingroup\$

Charcoal, 35 bytes

NθNηNζIΣEΦEX²θ⍘ι²№ι×1η×XζΣιX⁻¹ζ⁻θΣι

Try it online! Link is to verbose version of code. Explanation: Brute force approach.

NθNηNζ

Input n, k and p.

EX²θ⍘ι²

Get all of the binary numbers up to 2ⁿ.

Φ...№ι×1η

Keep only those with k consecutive heads, I mean 1s.

E...×XζΣιX⁻¹ζ⁻θΣι

Compute their probabilities, which is pⁱ(1-p)ⁿ⁻ⁱ, where i is the number of 1s in each relevant row.

IΣ...

Output the grand total.

A 40-byte port of @dancxviii's approach is much more efficient:

Nθ≔E⊕N¬ιηFN«≔⊟ηζ≔⁺⟦×Ση⁻¹θ⟧×ηθη⊞η⁺ζ⊟η»I⊟η

Try it online! Link is to verbose version of code. Takes inputs in the order p, k, n. Expects k>1 (+1 byte to support k=1). Explanation:

Nθ

Input p.

≔E⊕N¬ιη

Create an array of 1 1 and k 0s.

FN«

Repeat n times.

≔⊟ηζ

Remove the last entry of the array.

≔⁺⟦×Ση⁻¹θ⟧×ηθη

Multiply the array by p, then prefix the sum of the array multiplied by 1-p.

⊞η⁺ζ⊟η

Add the previous last entry to the current last entry.

»I⊟η

Output the last entry in the array.

I had been planning on working out how to solve the problem using dynamic programming but I'm pretty sure it would have resulted in the same algorithm.

\$\endgroup\$
2
\$\begingroup\$

Husk, 30 19 17 bytes

ṁoΠm≠²fö≥⁰▲mΣgπḋ2

Try it online!

Port of my R answer.
Input is arg1=p, arg2=k, arg3=n.

                πḋ2     # cartesian arg-3-th power of binary digits of 2
                        # so: all possible combinations of n heads/tails;
        fö              # now consider only those with runs: filter by
            ▲           # the maximum of
             mΣg        # the sums of each group of identical elements
          ≥⁰            # is greater than or equal to arg-2;
                        # now convert the 1s & 0s to probabilities:
   mo -²                # subtract the probability of heads from each 
     a                  # and get the absolute values
                        # so 1s become (1-p) and 0s become p;
ṁo                      # finally map across each group & sum the result
  Π                     # product of the probabilities
\$\endgroup\$
1
\$\begingroup\$

Wolfram Language (Mathematica), 59 bytes

kpIf[#>=k,(1-Boole[#>k]p)(1-#0[#-k-1])p^k+#0[#-1],0]&

Try it online!

Port of Python. Takes f(k)(p)(n).

\$\endgroup\$
0
\$\begingroup\$

J, 49 bytes

1 :'0{_1{[:+/ .*^:(u-1)~(1,~]#0:),.~(*=@i.),~1-['

Try it online!

Almost identical to dancxviii's Markov chain approach.

\$\endgroup\$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.