OpenEuphoria

Re: Gaussian distribution

• Posted by "Christian Cuvier" <Christian.CUVIER at agriculture.gouv.fr> Aug 04, 2004
• 868 views

CK,
The following gives you exact probabilities rather fast for arbitrary outco=
me,
dice nu=F9bers and regular dice shapes.

Contrary to earlier bad habits, I tested it before posting :)

For q>15 or 20, the gaussian approx is probably faster, and is precise enou=
gh.
The usual n=12 limitation is caused by the fact that 13! needs too many b=
its.
The floor value for the gaussian approximation to be good enough depends on=

the number of dice actually.

type pint(integer n) return n>0 end type

constant PMAX=170
--the biggest integer q such that q! fits into a atom, or nearly so

function gpart(pint p, pint q, pint r)
--return a list of all partitions of p in q parts, each of them not more th=
an r
sequence res,s integer k0
if p>q*r or p<q then return {} end if
if p=q then return {{{1,p}}} end if
if p=1 then
    if q=1 then return {{{1,1}}} else return {} end if
end if
if q=1 then
    if p<=r then return {{{p,1}}} else return {} end if
end if
if r=1 then return {{{1,p}}} end if
k0=(remainder(p,r)=0)
if p=q*r then res={{{r,q}}}
else res={}
end if
for i=1 to floor(p/r)-k0 do
     s=gpart(p-i*r,q-i,r-1)
     if i>0 then for j=1 to length(s) do
        s[j]=append(s[j],{r,i})
     end for end if
     res &= s
end for
res &= gpart(p,q,r-1)   --the i=0 case
return res
end function

global function chance(pint p,pint q,pint r)
--what's the chance to get p as the sum of q balanced dice with r faces
--numbered 1..r?
atom a,a0,a1
sequence res
pint p2,pmax
sequence fact
if p>q*r or p<q then return 0 end if
pmax=p
if pmax>PMAX then pmax=PMAX end if
fact=repeat(1,pmax)
for i=2 to pmax do fact[i]=i*fact[i-1] end for
a1=1/r
if q=1 then return a1
elsif q=2 then
    if p>r+1 then p=2*r+1-p end if
    return p*a1*a1
end if
a0=q/r a=a0*a1
for i=2 to q-1 do a0-=a1 a*=a0 end for
a0=0
res=gpart(p,q,r)
for i=1 to length(res) do
     a1=a
     for j=1 to length(res[i]) do
        p2=res[i][j][2]
        if p2>1 then
           if p2<=pmax then a1/=fact[p2]
           else
              a1/=fact[pmax]
              for k=pmax+1 to p2 do a1/=k end for
           end if
        end if
     end for
     a0+=a1
end for
return a0
end function

--some test code
include get.e
atom a a=0
atom a0
for i=5 to 20 do
     a0=chance(i,5,4)
     ?a0
     a+=a0
end for
?a
integer w
w=wait_key()

CChris