Equation for drawing a bell curve?

Author

Message

Rich Dersheimer

AGK Developer

17

Years of Service

User Offline

Joined: 1st Jul 2009

Location: Inside the box

Posted: 31st May 2011 13:10 Edited at: 31st May 2011 13:17

Link

I need a simple formula for drawing this:

Something along the lines of

+ Code Snippet

for y = 1 to 100
  x=(insert equation here)
  line x,y,100,y
next y

I know it's a normal distribution curve, but I have zero math training, and could use help.

Neuro Fuzzy? Anybody?

Back to top

Profile PM Email

Hodgey

16

Years of Service

User Offline

Joined: 10th Oct 2009

Location: Australia

Posted: 31st May 2011 13:27

Link

Hey Rich, this is the closest thing I could get to your diagram.

+ Code Snippet

x as float

for n = -275 to -100
    x = n
    y = sin(x)*100
    dot x+300, y + screen width() / 2
next n

wait key
end

Hope it helps. If you want more info, look up the sine curve.

A clever person solves a problem, a wise person avoids it - Albert Einstein

Back to top

Profile PM Email

Sven B

21

Years of Service

User Offline

Joined: 5th Jan 2005

Location: Belgium

Posted: 31st May 2011 13:59 Edited at: 1st Jun 2011 09:23

Link

The equation for the gaussian curve is
exp(-x^2)

Of course, you can transform this equation in any way you like.
Shift: exp(-(x-x1)^2)
Scale: A * exp(-x^2)
Slopes: exp(-x^2/s)
or a combination
A * exp(-(x-x1)^2/s)

[edit] For a gaussian distribution curve:
y = A * exp(-(x - x1)^2 / (2*s))
A = 1 / (sigma * sqrt(2 * pi))
x1 = mu
s = sigma
Where mu is the mean (expected value) and sigma^2 is the variance. By using these values, you can ensure that the area under the gauss curve is 1.

Cheers!
Sven B

Back to top

Profile PM Email

Neuro Fuzzy

19

Years of Service

User Offline

Joined: 11th Jun 2007

Location:

Posted: 1st Jun 2011 02:26

Link

Quote: "I know it's a normal distribution curve, but I have zero math training, and could use help.

Neuro Fuzzy?"

xD I like having a reputation, but I don't know much calculus or probability stuff. The stuff I posted on this was a result of googling

Sven B's post covers everything about drawing the curve, but there's another use for this function in programming.

The program I posted here:
http://forum.thegamecreators.com/?m=forum_view&t=169807&b=6

has a useful function for use with the normal distribution. Basically... curves whose total enclosed area is 1 can be considered probability distributions. Then you have a function f(), which returns a value on the x axis of the given probability distribution. The probability of the return value being inbetween and points x1 and x2 is equal to the area under the curve inbetween x1 and x2.

So, the code snippet in that thread takes the result of the function sNormInv and draws a line there. The result of sNormInv could, theoretically, be anywhere on the number line, but due to the complex math which I copied and pasted from the corners of the interweb, it /probably/ lies near the origin.

this is the code snippet in question:

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

+ Code Snippet

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

soo... I use this function if I want a bunch of objects clustered around an area, but nowhere really specific. just position an object at {normpoint(),normpoint(),normpoint()}, and it will probably be close to the origin. if you make a bunch of objects positioned like that, they will all be clustered around the origin.

Summarizing what I just said, copy the normpoint and SNormInverse functions to your project, then position an object at normpoint to have it closeish to the origin (probably).

Tell me if there's a broken link to images in a thread I post, and I'll fix 'em.

Back to top

Profile PM Email

Rich Dersheimer

AGK Developer

17

Years of Service

User Offline

Joined: 1st Jul 2009

Location: Inside the box

Posted: 1st Jun 2011 03:41 Edited at: 1st Jun 2011 03:43

Link

Thanks Hodgey, that got me to what I needed. I ended up with this code:

+ Code Snippet

x# as float

for n = 100 to 440
    ink rgb(200-(n/3),0,255-(n/4)),0
    x# = n
    y = sin(x#)*100
    line x#+150,y+screen width()/2,x#+150,500
next n

which gave me this:

@ SvenB and Neuro Fuzzy - a little over my head, but I promise I'll give it some study, maybe I can learn some new tricks after all.

Back to top

Profile PM Email

Neuro Fuzzy

19

Years of Service

User Offline

Joined: 11th Jun 2007

Location:

Posted: 1st Jun 2011 04:09 Edited at: 1st Jun 2011 04:20

Link

what about

+ Code Snippet

x# as float

for n = 0 to 400
    ink rgb(255,0,lerp(0,255,n/400.0)),0
    y#=normCurve((n-200)/100.0)*100 //from -2 to 2, scaled vertically a bit
    line n,0,n,y#
next n
sync
wait key
end
function lerp(a as integer, b as integer, t as float)
	n#=(b-a)*t+a
endfunction n#

function normCurve(x as float)
	n#=exp(-x*x)
endfunction n#

[edit]
or demonstrating the normal distribution whose total area is 1:

+ Code Snippet

x# as float

d as float
d=0
do
	cls
	d=wrapvalue(d+1)
	v#=(sin(d)+1)*2
	for n = 0 to 400
    	ink rgb(255,0,lerp(0,255,n/400.0)),0
    	y#=normCurve((n-200)/100.0,v#)*100 //from -2 to 2, scaled vertically a bit
    	line n,0,n,y#
	next n
	sync
loop
end

function lerp(a as integer, b as integer, t as float)
	n#=(b-a)*t+a
endfunction n#

function normCurve(x as float, variance2 as float)
	`variance2 is 1/variance^2
	n#=0.39894228*variance2*exp(-x*x*variance2)
endfunction n#

[edit2]
and keep in mind, I really don't know /why/ that magic number (~.3989=1/sqrt(2*pi)) is there, I just read that it is so I put it there. It's not important to completely understand the math, you just need to write functions that plug numbers in.

Tell me if there's a broken link to images in a thread I post, and I'll fix 'em.

Back to top

Profile PM Email

Rich Dersheimer

AGK Developer

17

Years of Service

User Offline

Joined: 1st Jul 2009

Location: Inside the box

Posted: 1st Jun 2011 04:46

Link

@Neuro Fuzzy - that first snippet is cool, but my jaw just hit the floor when I saw that second one. Kind of creepy, like blood about to drip from the ceiling.

Back to top

Profile PM Email

Sven B

21

Years of Service

User Offline

Joined: 5th Jan 2005

Location: Belgium

Posted: 1st Jun 2011 09:22 Edited at: 1st Jun 2011 13:25

Link

Quote: "and keep in mind, I really don't know /why/ that magic number (~.3989=1/sqrt(2*pi)) is there, I just read that it is so I put it there. It's not important to completely understand the math, you just need to write functions that plug numbers in."

Hi NeuroFuzzy,

The sqrt(2*pi) is the area under the function e^(-x^2/2) (in other words, the integral from -infinity to infinity). The actual calculation of the area is kind of hard though, because 'normal' techniques won't cut it. I'm sure we had to do it once, but I forgot in which book they did it...
By the way, I noticed that I forgot a 1/2 in the exponent of the gaussian curve formula there. My bad! I'm editing it in a second.

Just thought I'd throw this in.
Sven

[edit] Forgot /2 again... exp(-x^2/2)

Back to top

Profile PM Email

Hodgey

16

Years of Service

User Offline

Joined: 10th Oct 2009

Location: Australia

Posted: 1st Jun 2011 09:24

Link

Quote: "Thanks Hodgey, that got me to what I needed."

You're welcome, glad it helped

A clever person solves a problem, a wise person avoids it - Albert Einstein

Back to top

Profile PM Email

Bluespark

15

Years of Service

User Offline

Joined: 18th Mar 2011

Location:

Posted: 2nd Jun 2011 01:16

Link

It's 1/Sqrt[2*Pi]. One way is square the integral and change it to polar form. Since the limit of an infinitely large circle and square are the same, the polar form gives 1/(2*Pi). Since we squared, take the square root and that's where the result comes from.

Back to top

Profile PM

Sven B

21

Years of Service

User Offline

Joined: 5th Jan 2005

Location: Belgium

Posted: 2nd Jun 2011 09:21

Link

Quote: "It's 1/Sqrt[2*Pi]. One way is square the integral and change it to polar form. Since the limit of an infinitely large circle and square are the same, the polar form gives 1/(2*Pi). Since we squared, take the square root and that's where the result comes from."

Ah yes I remember now

. We enclosed the square in two concentric circles to prove that the result converges to the same value as the square.

Cheers!
Sven

Back to top

Profile PM Email

Diggsey

20

Years of Service

User Offline

Joined: 24th Apr 2006

Location: On this web page.

Posted: 2nd Jun 2011 11:53

Link

You can get a decent curve very simply with this:
(1-x²)²

In the range 0 <= x <= 1 it will look like that curve and it's often used in metaballs because it's the simplest equation that gives a smooth transition from 0 to 1.

[b]

Back to top

Profile PM Email Website

Phaelax

DBPro Master

23

Years of Service

User Offline

Joined: 16th Apr 2003

Location: Metropia

Posted: 2nd Jun 2011 14:57

Link

When drawing the curve, you can speed up drawing immensely by replacing LINE with BOX. In Neuro's 2nd snippet of the 'dripping' curve, frame rates went from about 100 to 3300.

Back to top

Profile PM Email Website

Sorry your browser is not supported!

DarkBASIC Professional Discussion / Equation for drawing a bell curve?