# Weighted coin flip strings

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

• That's a terse description for sure, but the task is actually pretty clear IMO. Jun 30 at 14:26
• @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. Jun 30 at 14:29
• Are we guaranteed that k isn't greater than n?
– xnor
Jun 30 at 15:21
• 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 :-) Jun 30 at 15:43
• @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. Jul 1 at 7:04

# 05AB1E, 403918 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)

• Did not know about this formula, was using the markov method and was looking for a trick related to decomposition, this formula is cleaner Jun 30 at 16:17
• 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... Jul 1 at 7:55
• @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. Jul 1 at 8:08
• 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... Jul 1 at 9:17
• @DominicvanEssen Ah wait, I was looking at your 92 bytes function instead of 116 bytes.. I will edit my answer to also credit you. Jul 1 at 9:34

# 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.

• It looks like you can save the space after and by rearranging terms
– xnor
Jun 30 at 22:45

# 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:
ṣ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!

# R, 129122 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...

• 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. Jun 30 at 15:30
• @pxeger - I've added an exact solution, too. Jun 30 at 16:21
• -2 bytes Jun 30 at 20:28
• @pajonk - Thanks! The which(array(),T) trick instead of expand.grid is really nice! Jun 30 at 21:00
• Port of loopy walt's python answer is just 70 bytes. Jul 1 at 5:49

# Vyxalḋ, 23 bytes

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


Try it Online!

Port of Jelly. Really slow.

# Vyxalḋ, 41 bytes

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


Try it Online!

Port of 05AB1E.

# 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.

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})}\$$.

# 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!

# 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])

• Welcome to Code Golf, and nice answer! Jul 1 at 2:03
• 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*. Jul 1 at 4:04

# 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}


Thanks to Kevin Cruijssen for finding this formula.

# 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.

• How many bytes is it without the whitespace? :P Jun 30 at 19:16
• @thejonymyster 72 bytes. I prefer programs that compile, though. Jun 30 at 20:45

# 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.

# Charcoal, 35 bytes

ＮθＮηＮζＩΣＥΦＥＸ²θ⍘ι²№ι×1η×ＸζΣιＸ⁻¹ζ⁻θΣι


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

ＮθＮηＮζ


Input n, k and p.

ＥＸ²θ⍘ι²


Get all of the binary numbers up to 2ⁿ.

Φ...№ι×1η


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

Ｅ...×ＸζΣιＸ⁻¹ζ⁻θΣι


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

ＩΣ...


Output the grand total.

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

Ｎθ≔Ｅ⊕Ｎ¬ιηＦＮ«≔⊟ηζ≔⁺⟦×Ση⁻¹θ⟧×ηθη⊞η⁺ζ⊟η»Ｉ⊟η


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:

Ｎθ


Input p.

≔Ｅ⊕Ｎ¬ιη


Create an array of 1 1 and k 0s.

ＦＮ«


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.

»Ｉ⊟η


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.

# Husk, 3019 17 bytes

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


Try it online!

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


# 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).

# J, 49 bytes

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


Try it online!

Almost identical to dancxviii's Markov chain approach.