User:DroneBetter/object occurrences and commonness bounds

This page is primarily for my own benefit, since I have pondered upon it many times throughout my wiki-editing, and rederived fragments of this over and over again. (I am somewhat unknowledgeable about statistics so this may all be replicating the work of bygone statisticians of yore)

Suppose one has a census with n soups, and an object that has probability p of occurring. Its probability of having occurred ≤ k of the n times is

. ∑ki=0((ni)*pⁱ*(1-p)^n-i)
= ∑ki=0((nk-i)*p^k-i*(1-p)^n-k+i) #reverse iteration order
= p^k*(1-p)^n-k*∑ki=0((nk-i)*(1-pp)ⁱ)
= p^k*(1-p)^n-k*∑ki=0(n!(k-i)!*(n-k+i)!*(1-pp)ⁱ)
= n!k!*p^k*(1-p)^n-k*∑ki=0(k!(k-i)!*(1-pp)ⁱ(n-k+i)!)
= n!k!*p^k*(1-p)^n-k*∑ki=0(∏kj=k-i+1(j)*(-1)ⁱ*(p-1p)ⁱ(n-k+i)!) #split the 1-p into -1*(p-1) (very suspicious)
= n!k!*p^k*(1-p)^n-k*∑ki=0(∏i-k-1j=-k(j)*(p-1p)ⁱ(n-k+i)!) #there are i terms in the product so we can negate each one
= (nk)*p^k*(1-p)^n-k*∑∞i=0(∏i-1j=0(j-k)*i!*(n-k)!(n-k+i)!*(p-1p)ⁱi!) #now use the fact that where i=k+1, max(j)=k, at which the multiplicand is 0 so the product will be also, upper limit need not be
= (nk)*p^k*(1-p)^n-k*∑∞i=0(∏i-1j=0((j+1)*(j-k)j+n-k+1)*(p-1p)ⁱi!)
= (nk)*p^k*(1-p)^n-k*₂F₁(1,-kn-k+1;p-1p)

however I don't know how to approximate hypergeometrics so (while very cool) this is not especially useful

herein we shall use the standard definition erf=λ x: 2*∫₀^x(e^-x2)√π

Due to the de Moivre–Laplace theorem and sum-integral convergence, where n is large (the sample size one might find in a Catagolue census, for instance), we have the approximation

. ∑ki=0((ni)*pⁱ*(1-p)^n-i)
~ ∫₀^k(e^{-(x-n*p)2/(2*n*p*(1-p))} dx)√2*π*n*p*(1-p)
= erf(k-n*p√2*n*p*(1-p))-erf(-n*p√2*n*p*(1-p))2
~ erf(k-n*p√2*n*p*(1-p))+12 #much worse for extremal k so not used

so there is about a chance of t that it occurs k ≤ n*p + erf^-1(2*t-erf(n*p√2*n*p*(1-p)))*√2*n*p*(1-p) times

so one may use the following Python function to obtain bounds upon p from required confidence intervals (ie. papprox(n,k,0.025) and papprox(n,k,0.975) for two-sided 5%)

from functools import reduce
from math import sqrt
from sympy import erf,erfinv
cdf=lambda n,p,t: n*p-erfinv(t-erf(n*p/sqrt(2*n*p*(1-p))))*sqrt(2*n*p*(1-p))
bisect=lambda f,out,it: reduce(lambda r,i: r+(-1)**(float(f(r))>out)/2**(i+1),range(it),1)
papprox=lambda n,k,t: bisect(lambda r: cdf(n,r,t),k,1<<6)/2

For instance, as of this page's publication on February 20, 2024, there are 264000550576120 soups in b3s23/C1, and 5161997 occurrences of the trans-queen bee shuttle, from which one could infer they occur in 1.95529*10^-8 ≈ 151143104 soups, but in actuality at (the two-sided 5% confidence level) it is between about 1.95388*10^-8 ≈ 151187242 and 1.95698*10^-8 ≈ 151099004, so it is only appropriate to give it to a precision of 151100000.

However, with only two occurrences of the loafer at present, the small k hinders the approximation, so it provides 8.53085*10^-15 ≈ 1117221543808304 as the most probable (t=12) occurrence rate, whereas it should be 1132000275288060 (which predicts there are 2.2521504324140884 loafers), so this method is not feasible (whereas for the queen bee shuttle example it converges to 1.955297816136326*10^-08, which predicts that there are 5161997.000002751 trans-queen bees)