This function was really annoying to figure out how to get working. Basically... This is the normal distribution, equation:
graph:
Used for probability. Basically, the probability of a point being at a given x is y. Now... for that to be really useful, you need to find the inverse of the sum. This is where problems arise, because it's a non-elementary function (IE not able to be described in something nice like 5x*2^-x or whatever)
here is the graph of the cumulative normal distribution (the sum of the normal distribution at a given x)
Right, now her8e is the inverse of the cumulative normal function (solving for x, basically)
Wow. Right. So, really complicated calculus n stuff way over my head ATM, so, I turned to the internet, and found some equations, and wired 'em together to solve for the above graph. This is, of course, an estimation, because the whole thing is... blah.
Here's the code:
global pixels_per_unit as float
pixels_per_unit=100
sync on
r1 as float
r1=0
r2 as float
r2=0
m as float
m=0
a as float
a=0
do
cls
from#=0
ato#=1
for x#=from# to ato# step 0.001
dot x#*pixels_per_unit+screen width()/2,-SNormInv(x#)*pixels_per_unit+screen height()/2
next x#
drawGraphPaper(pixels_per_unit,screen width()/2,screen height()/2,4)
sync
loop
end
function drawGraphPaper(PPU as integer, midx as integer, midy as integer, subUnits as integer)
if subunits<=0 then subunits=1
if PPU<=0 then ppu=40
subDist as float
subDist=PPU*1.0/subUnits
counter1 as float =
counter1=-PPU
counter2 as float = 0
w=screen width()
h=screen height()
ink rgb(0,0,40)
line 0,h/2,w,h/2
ink rgb(0,40,0)
line w/2,0,w/2,h
while (counter1+midx<w) or (midx-counter1>0)
inc counter1, PPU
counter2=0
while (counter2<PPU)
inc counter2, subDist
ink rgb(150,150,150),0
line midx+counter1+counter2,midy-2,midx+counter1+counter2,midy+2
line midx-counter1+counter2,midy-2,midx-counter1+counter2,midy+2
endwhile
ink rgb(255,255,255),0
line midx+counter1,midy-5,midx+counter1,midy+5
line midx-counter1,midy-5,midx-counter1,midy+5
endwhile
counter1=-PPU
counter2=0
while (counter1+midy<w) or (midy-counter1>0)
inc counter1, PPU
counter2=0
while (counter2<PPU)
inc counter2, subDist
ink rgb(150,150,150),0
line midx-2,midy+counter1+counter2,midx+2,midy+counter1+counter2
line midx-2,midy-counter1-counter2,midx+2,midy-counter1-counter2
endwhile
ink rgb(255,255,255),0
line midx-4,midy+counter1,midx+4,midy+counter1
line midx-4,midy-counter1,midx+4,midy-counter1
endwhile
endfunction
Function SNormInv(p as float)
q as float
r as float
A1 as float
A2 as float
A3 as float
A4 as float
A5 as float
A6 as float
B1 as float
B2 as float
B3 as float
B4 as float
B5 as float
C1 as float
C2 as float
C3 as float
C4 as float
C5 as float
C6 as float
D1 as float
D2 as float
D3 as float
D4 as float
P_LOW as float
P_HIGH as float
`Coefficients in rational approximations.
A1= -39.69683
A2= 220.9461
A3= -275.9285
A4= 138.3578
A5= -30.66480
A6= 2.506628
B1= -54.47610
B2= 161.5858
B3= -155.6990
B4= 66.80131
B5= -13.28068
C1= -0.007785
C2= -0.322396
C3= -2.400758
C4= -2.549733
C5= 4.374664
C6= 2.938164
D1= 0.007784696
D2= 0.3224671
D3= 2.445134
D4= 3.754409
P_LOW= 0.02425
P_HIGH=1-P_LOW
if p>0 and p<P_LOW
`Rational approximation for lower region.
q = Sqrt(-2*Log(p))
ret#=(((((C1*q+C2)*q+C3)*q+C4)*q+C5)*q+C6)/((((D1*q+D2)*q+D3)*q+D4)*q+1)
exitfunction ret#
endif
If p>=P_LOW and p<=P_HIGH
`Rational approximation for central region.
q=p-0.5
r=q*q
ret#=(((((A1*r+A2)*r+A3)*r+A4)*r+A5)*r+A6)*q/(((((B1*r+B2)*r+B3)*r+B4)*r+B5)*r+1)
exitfunction ret#
endif
if p> _HIGH and p<1
`Rational approximation for upper region.
q = Sqrt(-2*Log(1-p))
ret#=-(((((C1*q+C2)*q+C3)*q+C4)*q+C5)*q+C6)/((((D1*q+D2)*q+D3)*q+D4)*q+1)
exitfunction ret#
endif
endfunction 0.0
You need a log() function to calculate this. You can either use IanM's matrixutils, or the ln() function I posted earlier (matrixutils would be better).
So... how is this function useful at all? I'll show you!!!
Say you want a point with the probability of being at any given x as defined by the normal distribution (first graph). Well... I'll give you that point!
function normPoint()
p#=(rnd(999997)+1)/999999.0
r#=SNormInv(p#)
endfunction r#
still not sold? Check this out:
global pixels_per_unit as float
pixels_per_unit=100
sync on
r1 as float
r1=0
r2 as float
r2=0
m as float
m=0
a as float
a=0
from#=-4
ato#=4
for x#=from# to ato# step 0.001
dot x#*pixels_per_unit+screen width()/2,-normalDistribution(x#)*pixels_per_unit+screen height()/2
next x#
do
for n=1 to 100
x#=normPoint()
line_increase_color(x#*pixels_per_unit+screen width()/2,6,x#*pixels_per_unit+screen width()/2,10,rgb(1,1,1))
next n
drawGraphPaper(pixels_per_unit,screen width()/2,screen height()/2,4)
sync
loop
end
function line_increase_color(x1,y1,x2,y2,c as dword)
dx = x2-x1
dy = y2-y1
pixel_increase_color(x1, y1,c)
if abs(dx) > abs(dy)
m# = (0.0 + dy)/dx
b# = y1 - m#*x1
if dx < 0
dx = -1
else
dx = 1
endif
while x1 <> x2
x1 = x1 + dx
pixel_increase_color(x1,m#*x1+b#,c)
endwhile
else
if dy <> 0
m# = (0.0 + dx)/dy
b# = x1 - m#*y1
if dy < 0
dy = -1
else
dy = 1
endif
while y1 <> y2
y1 = y1 + dy
pixel_increase_color(m#*y1+b#,y1,c)
endwhile
endif
endif
endfunction
function normalDistribution(x as float)
y#=0.39894228*2.71828183^(-x*x/2.0)
endfunction y#
function pixel_increase_color(x,y,c as dword)
c2 as dword
c2=point(x,y)
r as word
b as word
g as word
r=minimum(rgbr(c2)+rgbr(c),255)
g=minimum(rgbg(c2)+rgbg(c),255)
b=minimum(rgbb(c2)+rgbb(c),255)
ink rgb(r,g,b),0
dot x,y
endfunction
function minimum(a,b)
if a<b then exitfunction a
endfunction b
function drawGraphPaper(PPU as integer, midx as integer, midy as integer, subUnits as integer)
if subunits<=0 then subunits=1
if PPU<=0 then ppu=40
subDist as float
subDist=PPU*1.0/subUnits
counter1 as float =
counter1=-PPU
counter2 as float = 0
w=screen width()
h=screen height()
ink rgb(0,0,40)
line 0,h/2,w,h/2
ink rgb(0,40,0)
line w/2,0,w/2,h
while (counter1+midx<w) or (midx-counter1>0)
inc counter1, PPU
counter2=0
while (counter2<PPU)
inc counter2, subDist
ink rgb(150,150,150),0
line midx+counter1+counter2,midy-2,midx+counter1+counter2,midy+2
line midx-counter1+counter2,midy-2,midx-counter1+counter2,midy+2
endwhile
ink rgb(255,255,255),0
line midx+counter1,midy-5,midx+counter1,midy+5
line midx-counter1,midy-5,midx-counter1,midy+5
endwhile
counter1=-PPU
counter2=0
while (counter1+midy<w) or (midy-counter1>0)
inc counter1, PPU
counter2=0
while (counter2<PPU)
inc counter2, subDist
ink rgb(150,150,150),0
line midx-2,midy+counter1+counter2,midx+2,midy+counter1+counter2
line midx-2,midy-counter1-counter2,midx+2,midy-counter1-counter2
endwhile
ink rgb(255,255,255),0
line midx-4,midy+counter1,midx+4,midy+counter1
line midx-4,midy-counter1,midx+4,midy-counter1
endwhile
endfunction
function normPoint()
p#=(rnd(999997)+1)/999999.0
r#=SNormInv(p#)
endfunction r#
Function SNormInv(p as float)
q as float
r as float
A1 as float
A2 as float
A3 as float
A4 as float
A5 as float
A6 as float
B1 as float
B2 as float
B3 as float
B4 as float
B5 as float
C1 as float
C2 as float
C3 as float
C4 as float
C5 as float
C6 as float
D1 as float
D2 as float
D3 as float
D4 as float
P_LOW as float
P_HIGH as float
`Coefficients in rational approximations.
A1= -39.69683
A2= 220.9461
A3= -275.9285
A4= 138.3578
A5= -30.66480
A6= 2.506628
B1= -54.47610
B2= 161.5858
B3= -155.6990
B4= 66.80131
B5= -13.28068
C1= -0.007785
C2= -0.322396
C3= -2.400758
C4= -2.549733
C5= 4.374664
C6= 2.938164
D1= 0.007784696
D2= 0.3224671
D3= 2.445134
D4= 3.754409
P_LOW= 0.02425
P_HIGH=1-P_LOW
if p>0 and p<P_LOW
`Rational approximation for lower region.
q = Sqrt(-2*Log(p))
ret#=(((((C1*q+C2)*q+C3)*q+C4)*q+C5)*q+C6)/((((D1*q+D2)*q+D3)*q+D4)*q+1)
exitfunction ret#
endif
If p>=P_LOW and p<=P_HIGH
`Rational approximation for central region.
q=p-0.5
r=q*q
ret#=(((((A1*r+A2)*r+A3)*r+A4)*r+A5)*r+A6)*q/(((((B1*r+B2)*r+B3)*r+B4)*r+B5)*r+1)
exitfunction ret#
endif
if p> _HIGH and p<1
`Rational approximation for upper region.
q = Sqrt(-2*Log(1-p))
ret#=-(((((C1*q+C2)*q+C3)*q+C4)*q+C5)*q+C6)/((((D1*q+D2)*q+D3)*q+D4)*q+1)
exitfunction ret#
endif
endfunction 0.0
it creates a strip of white lines, more dense at the middle, less dense at the edges.
Bon appétit!
[edit]
and an image showing an implementation visually
