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
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.
Original mathematical 40 39 bytes approach:
÷L®s©>m²®¹*-D>®/$-*I+s®<cI®¹*m$-®<m)øPO
Inputs in the order \$k,n,p\$.
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)
1e19 times, and check if at least \$k\$ unfair coin heads are present in a random result. Both are brute-force, but different.
\$\endgroup\$
Commented
Jul 1, 2022 at 8:08
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)
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.
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.
and by rearranging terms
\$\endgroup\$
Ø.ṗṣ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.)
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
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)))
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
}
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.
function(n,k,p,m=cbind(rbind(1-p,diag(k)*p),!k:0))Reduce(`%*%`,rep(list(m),n))[k+1]
Port of dancxviii's Markov method - upvote that one - but sadly still longer than pajonk's R answer...
which(array(),T) trick instead of expand.grid is really nice!
\$\endgroup\$
Commented
Jun 30, 2022 at 21:00
k->p->g(n)=if(n>=k,p^k*(1-p*(n>k))*(1-g(n---k))+g(n))A port of loopy walt's Python answer. Takes input in the the form (k)(p)(n)
f(n,k,p)=Pol((1-y=p*x)/(t=1-x)/(t/y^k+x-y)+O(x^n*x))\x^nThe 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})}\$.
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)
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)
// \_____/
// cachePort 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]
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])
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\$
Commented
Jul 1, 2022 at 4:04
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.
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 Σ ]
Translation of the formula given in Kevin Cruijssen's 05AB1E answer.
`[`=function(k,n,p)`if`(n<k,0,(1-p*(n>k))*(1-k[n-k-1,p])*p^k+k[n-1,p])
Port of @loopy walt's Python answer.
For different approaches in R see @Dominic van Essen's answer.
ＮθＮηＮζＩΣＥΦＥＸ²θ⍘ι²№ι×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.
ṁoΠm≠²fö≥⁰▲mΣgπḋ2
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
kpIf[#>=k,(1-Boole[#>k]p)(1-#0[#-k-1])p^k+#0[#-1],0]&
Port of Python. Takes f(k)(p)(n).
1 :'0{_1{[:+/ .*^:(u-1)~(1,~]#0:),.~(*=@i.),~1-['
Almost identical to dancxviii's Markov chain approach.
kisn't greater thann? \$\endgroup\$