PureBasic Forums - English

http://www.purebasic.com
https://www.purebasic.fr/english/

Flat vs Gaussian Distribution (Demo)

https://www.purebasic.fr/english/viewtopic.php?t=40255

Page 1 of 1

Flat vs Gaussian Distribution (Demo)

Posted: Thu Dec 10, 2009 4:53 am
by kenmo
I was recently considering the use of Gaussian random number generators in games, to make randomized aspects a bit more natural. I did some quick research and found a simple algorithm (here) and implemented it in PB as follows:

Code: Select all

Procedure.d RandomFlat(Max.i = 100000)
  Protected Result.d
  If (Max < 1)
    Result = 0.0
  Else
    Result = Random(Max - 1)*1.0/Max
  EndIf
  ProcedureReturn Result
EndProcedure

Procedure.d RandomGauss(Scale.d = 1.0, Max.i = 100000)
  Protected x1.d, x2.d, w.d, Valid.i, Result.d
  If (Scale <= 0.0)
    Result = d
  Else
    Repeat
      Valid = #True
      Repeat
        x1 = 2.0 * RandomFlat(Max) - 1.0
        x2 = 2.0 * RandomFlat(Max) - 1.0
        w  = (x1 * x1) + (x2 * x2)
      Until (w < 1.0)
      w = Sqr( (-2.0 * Log(w)) / w)
      x2 = (x1 * w / 5.0 / Scale) + 0.5
      If (x2 < 0.0) Or (x2 > 1.0)
        Valid = #False
      EndIf
    Until Valid
    Result = x2
  EndIf
  ProcedureReturn Result
EndProcedure
Both functions return a double that, 0 <= result < 1. RandomFlat() essentially uses PB's built-in Random(), which spreads out the results evenly. RandomGauss() instead spreads its results over a normal curve, centered around 0.5. An optional "scale" parameter (positive!) allows more control: values greater than 1 focus the curve tighter around the center, values less than 1 spread it out more. (The Max parameter is what goes into PB's Random(), higher Max means more possible results).

Why would you want Gaussian random numbers? For example (this is what I was pondering the other day), say you're making a game and you want the size of enemies to vary randomly. Instead of encountering average-sized, giant, and scrawny enemies with the same odds, Gaussian distribution means most enemies are about average, though huge and tiny ones still occur. Nifty, huh?

Now the graphical, interactive demo!

Code: Select all

; +----------------------------+-----------------------+
; | Random Number Distribution | Ryan 'kenmo' Toukatly |
; +----------------------------+-----------------------+
; | 12.9.2009 - Creation

;   This program generates random values on the interval [0, 1)
;   using two methods: flat distribution and normal (Gaussian)
;   distribution. Both are demonstrated graphically.

; ================ YOU CAN MODIFY THESE ===============================================

#Intervals  = 31 ; Number of bars in the histogram
#Iterations =  9 ; Number of samples to generate per cycle
#Timeout    = 20 ; Delay in milliseconds (when no window events occur)

; ================ RANDOM NUMBER PROCEDURES ===========================================

Procedure.d RandomFlat(Max.i = 100000)
  Protected Result.d
  If (Max < 1)
    Result = 0.0
  Else
    Result = Random(Max - 1)*1.0/Max
  EndIf
  ProcedureReturn Result
EndProcedure

Procedure.d RandomGauss(Scale.d = 1.0, Max.i = 100000)
  Protected x1.d, x2.d, w.d, Valid.i, Result.d
  If (Scale <= 0.0)
    Result = d
  Else
    Repeat
      Valid = #True
      Repeat
        x1 = 2.0 * RandomFlat(Max) - 1.0
        x2 = 2.0 * RandomFlat(Max) - 1.0
        w  = (x1 * x1) + (x2 * x2)
      Until (w < 1.0)
      w = Sqr( (-2.0 * Log(w)) / w)
      x2 = (x1 * w / 5.0 / Scale) + 0.5
      If (x2 < 0.0) Or (x2 > 1.0)
        Valid = #False
      EndIf
    Until Valid
    Result = x2
  EndIf
  ProcedureReturn Result
EndProcedure

; ================ DEMO PROGRAM =======================================================

Global IW.l, IH.l, Maxx.i, Total.i, FID.i

Macro wp(f)
  Int(IW*(f))
EndMacro
Macro hp(f)
  Int(IH*(f))
EndMacro

If (#Intervals < 1) Or (#Iterations < 1)
  End
EndIf
Global Dim Count.i(#Intervals - 1)
Global Dim DispH.f(#Intervals - 1)

Procedure StretchImage(Update.i = #True)
  IW = WindowWidth(0)
  IH = WindowHeight(0) - 110
  FID = LoadFont(0, "Monospace", IH/20)
  If CreateImage(0, IW, IH) And StartDrawing(ImageOutput(0))
    Box(0, 0, IW, IH, #White)
    StopDrawing()
  EndIf
  ResizeGadget(7, 0, 110, IW, IH)
  SetGadgetAttribute(7, #PB_Button_Image, ImageID(0))
EndProcedure

Procedure UpdateImage()
  If StartDrawing(ImageOutput(0))
    DrawingMode(#PB_2DDrawing_Transparent)
    DrawingFont(FID)
    
    Box(0, 0, IW, IH, #White)
    Box(wp(0.01), hp(0.02), wp(0.98), hp(0.96), $E0E0E0)
    
    For i = 0 To #Intervals - 1
      If i = 0
        sx.i = 0
      Else
        sx = ex
      EndIf
      ex = Int(i*IW*0.94/(#Intervals-1))
      
      hf.f = Count(i)*1.0/Maxx
      hmax.i = hp(0.90)
      h.i = hmax * hf
      DispH(i) = DispH(i) + (h - DispH(i))/5.0
      h = DispH(i)
      If h > hmax : h = hmax : DispH(i) = hmax : EndIf
      
      If h > 0
        Box(wp(0.03)+sx, hp(0.95)-h, ex-sx, h, RGB(0, Int(hf*128), 255))
      EndIf
    Next i
    
    If Total > 0
      DrawText(wp(0.05), hp(0.05), "N = " + Str(Total), #Red)
    EndIf
    StopDrawing()
    SetGadgetState(7, ImageID(0))
  EndIf
EndProcedure

WinFlags.i = #PB_Window_Invisible|#PB_Window_ScreenCentered|#PB_Window_MinimizeGadget
WinFlags   | #PB_Window_SizeGadget|#PB_Window_MaximizeGadget
If Not OpenWindow(0, 0, 0, 300, 300, "Random Number Distribution", WinFlags)
  MessageRequester("Error", "Window could not be opened.", #MB_ICONHAND) : End
EndIf

Frame3DGadget(0, 5, 5, 290, 100, "Random Number Generator")

OptionGadget(1, 15, 25, 200, 20, "Flat Distribuation (native PureBasic)")
OptionGadget(2, 15, 45, 140, 20, "Gaussian (bell curve)")
  SetGadgetState(2, #True)

TextGadget(3, 160, 48, 40, 17, "Scale:")
StringGadget(4, 200, 45, 70, 20, "1.0")

ButtonGadget(5, 75, 70, 70, 25, "Start", #PB_Button_Default)
ButtonGadget(6, 155, 70, 70, 25, "Reset")

CreateImage(0, 1, 1)
ImageGadget(7, 0, 0, 0, 0, ImageID(0))

Running.i = 0
Type.i = 2
BellScale.d = 1.0
Maxx.i = 100
Total = 0
ExitFlag.i = #False

SmartWindowRefresh(0, #True)
WindowBounds(0, 300, 300, #PB_Ignore, #PB_Ignore)

StretchImage()
HideWindow(0, #False)
UpdateImage()

Repeat
  Event.i = WaitWindowEvent(#Timeout)
  
  While Event
    
    If Event = #PB_Event_CloseWindow
      ExitFlag = #True
    ElseIf Event = #PB_Event_SizeWindow
      StretchImage()
      If Not Running : UpdateImage() : EndIf
    ElseIf Event = #PB_Event_Gadget
      ID.i = EventGadget()
      Select ID
        Case 1, 2 : Type = ID
        Case 4
          If EventType() = #PB_EventType_Change
            St.s = GetGadgetText(4)
            If ValD(St) > 0.0
              BellScale = ValD(St)
            EndIf
          EndIf
        Case 5, 7
          Running = 1 - Running
          If Running
            SetGadgetText(5, "Stop")
          Else
            SetGadgetText(5, "Run")
          EndIf
        Case 6
          For i = 0 To #Intervals - 1
            Count(i) = 0
          Next i
          Maxx = 100
          Total = 0
          If Not Running : UpdateImage() : EndIf
      EndSelect
    EndIf
    
    Event = WindowEvent()
  Wend
  
  If Running And (Total < 2000000000)
    For a = 1 To #Iterations
      If Type = 1
        x.d = RandomFlat()
      ElseIf Type = 2
        x.d = RandomGauss(BellScale)
      EndIf
      Group.i = Int(x * #Intervals)
      
      Count(Group) + 1
      If Count(Group) >= Maxx
        Maxx = Count(Group) + 1
      EndIf
    Next a
    Total + #Iterations
  
  EndIf
  UpdateImage()

Until ExitFlag

HideWindow(0, #True)
End

;
Issue I'm aware of: The lower the scale value is, the slower the procedure is, because I made up my own method that loops until a valid result is generated... This tiny delay is only a real issue if you use the procedure very often (ie. in your main program loop). (Although, if you're setting Scale very low, you may as well use Flat distribution!)

Also, the graphical demo isn't the most optimized thing in the world. Changing the Timeout constant will use more or less CPU power.

Have fun....!

Re: Flat vs Gaussian Distribution (Demo)

Posted: Thu Dec 10, 2009 9:03 am
by Kaeru Gaman
it's a quite complicated algorithm you use, should be really slow.

some time ago I took the approach to simulate a normalized distribution,
building a rather parabolic distribution using half a sinus.

Code: Select all

Procedure.l GauRand(Range.l)
    Range2 = Range / 2
    PiHalf.f = #PI / 2
    Prec = 500000
    Z = Range2 + Range2 * (ASin( (Random(Prec*2) - Prec ) / Prec ) / PiHalf)
    ProcedureReturn Z
EndProcedure
this makes one single random call and changes the resulting value to "bend" the flat distribution to a parabol...
when you need repetive values, you can take the precalculations out of the repetition...

Re: Flat vs Gaussian Distribution (Demo)

Posted: Thu Dec 10, 2009 11:19 am
by djes
Interesting!

Re: Flat vs Gaussian Distribution (Demo)

Posted: Thu Dec 10, 2009 6:46 pm
by Little John
There you go! :D

Regards, Little John

Re: Flat vs Gaussian Distribution (Demo)

Posted: Thu Dec 10, 2009 8:26 pm
by blueznl
Well, with all these very advanced options I feel almost ashamed presenting my own 'stick to the middle' solution... I just add up a set of randoms...

Code: Select all

  ;
  ; evenly spread 0..1000
  ;
  evenly.i =  random(1000)
  ;
  ; centered, mostly around 500 with some extremes
  ;
  centered.i = random(333)+random(333)+random(333)
  ;
Not fancy, not properly, but works :-)

Re: Flat vs Gaussian Distribution (Demo)

Posted: Fri Dec 11, 2009 7:35 pm
by kenmo
Cool, all interesting examples (and all faster than mine it seems).

I modified the code to show mine, Kaeru's, and Blueznl's all as overlapped line plots - let it run for a minute and it shows that my approximation tends to stay centered between the other two, although any look fine enough for most applications!

Code: Select all

; +----------------------------+-----------------------+
; | Random Number Distribution | Ryan 'kenmo' Toukatly |
; +----------------------------+-----------------------+
; | 12.09.2009 - Creation
; | 12.11.2009 - Line plot, multiple types at once

;   This program generates random values on the interval [0, 1)
;   using two methods: flat distribution and normal (Gaussian)
;   distribution. Both are demonstrated graphically.

; ================ YOU CAN MODIFY THESE ===============================================

#Intervals  = 31 ; Number of bars in the histogram
#Iterations =  9 ; Number of samples to generate per cycle
#Timeout    = 20 ; Delay in milliseconds (when no window events occur)

BellScale.d = 1.0 ; Scale factor (for RandomGauss())


#Types      =  4 ; Types of RNG's

Global Dim Name.s(#Types - 1)
  Name(0) = "blueznl"
  Name(1) = "kenmo"
  Name(2) = "kaeru"
  Name(3) = "pbasic"
  
Global Dim Col.i(#Types - 1)
  Col(0) = RGB(0, 0, 255)
  Col(1) = RGB(255, 128, 0)
  Col(2) = RGB(0, 192, 0)
  Col(3) = RGB(160, 0, 0)
  
Global Dim nx.l(#Types - 1)
Global Dim ny.l(#Types - 1)
Global Dim ox.l(#Types - 1)
Global Dim oy.l(#Types - 1)

; ================ RANDOM NUMBER PROCEDURES ===========================================

Procedure.d RandomFlat(Max.i = 100000)
  Protected Result.d
  If (Max < 1)
    Result = 0.0
  Else
    Result = Random(Max - 1)*1.0/Max
  EndIf
  ProcedureReturn Result
EndProcedure

Procedure.d RandomGauss(Scale.d = 1.0, Max.i = 100000)
  Protected x1.d, x2.d, w.d, Valid.i, Result.d
  If (Scale <= 0.0)
    Result = d
  Else
    Repeat
      Valid = #True
      Repeat
        x1 = 2.0 * RandomFlat(Max) - 1.0
        x2 = 2.0 * RandomFlat(Max) - 1.0
        w  = (x1 * x1) + (x2 * x2)
      Until (w < 1.0)
      w = Sqr( (-2.0 * Log(w)) / w)
      x2 = (x1 * w / 5.0 / Scale) + 0.5
      If (x2 < 0.0) Or (x2 > 1.0)
        Valid = #False
      EndIf
    Until Valid
    Result = x2
  EndIf
  ProcedureReturn Result
EndProcedure

Procedure.l GauRand(Range.l)
    Range2 = Range / 2
    PiHalf.f = #PI / 2
    Prec = 500000
    Z = Range2 + Range2 * (ASin( (Random(Prec*2) - Prec ) / Prec ) / PiHalf)
    ProcedureReturn Z
EndProcedure

Procedure.d Blueznl()
  ProcedureReturn (Random(333) + Random(333) + Random(333))/1000.0
EndProcedure

; ================ DEMO PROGRAM =======================================================

Global IW.l, IH.l, Maxx.i, Total.i, FID.i

Macro wp(f)
  Int(IW*(f))
EndMacro
Macro hp(f)
  Int(IH*(f))
EndMacro

If (#Intervals < 1) Or (#Iterations < 1)
  End
EndIf
Global Dim Count.i(#Types - 1, #Intervals - 1)
Global Dim DispH.f(#Types - 1, #Intervals - 1)

Procedure StretchImage(Update.i = #True)
  IW = WindowWidth(0)
  IH = WindowHeight(0)
  FID = LoadFont(0, "Monospace", IH/22)
  If CreateImage(0, IW, IH) And StartDrawing(ImageOutput(0))
    Box(0, 0, IW, IH, #White)
    StopDrawing()
  EndIf
  ResizeGadget(7, 0, 0, IW, IH)
  SetGadgetAttribute(7, #PB_Button_Image, ImageID(0))
EndProcedure

Procedure UpdateImage()
  If StartDrawing(ImageOutput(0))
    DrawingMode(#PB_2DDrawing_Transparent)
    DrawingFont(FID)
    
    Box(0, 0, IW, IH, #White)
    Box(wp(0.01), hp(0.02), wp(0.98), hp(0.96), $E0E0E0)
    
    For t.l = 0 To #Types - 1
    For i = 0 To #Intervals - 1
      ox(t) = nx(t) : oy(t) = ny(t)
      If i = 0
        sx.i = 0
      Else
        sx = ex
      EndIf
      ex = Int(i*IW*0.94/(#Intervals-1))
      
      hf.f = Count(t, i)*1.0/Maxx
      hmax.i = hp(0.90)
      h.i = hmax * hf
      DispH(t, i) = DispH(t, i) + (h - DispH(t, i))/5.0
      h = DispH(t, i)
      If h > hmax : h = hmax : DispH(t, i) = hmax : EndIf
      
      nx(t) = (wp(0.03)+sx + (ex-sx)/2.0)
      ny(t) = hp(0.95)-h
      If i > 0
        LineXY(ox(t), oy(t), nx(t), ny(t), Col(t))
      EndIf
      hh = hp(0.03)
      If h > hh
        Plot(nx(t), ny(t) + hh, Col(t))
      EndIf
      ;If h > 0
        ;Box(wp(0.03)+sx, hp(0.95)-h, ex-sx, h, RGB(0, Int(hf*128), 255))
      ;EndIf
    Next i
    Next t
    
    For t = 0 To #Types - 1
      DrawText(wp(0.03), hp(0.88 - 0.06*t), Name(t), Col(t))
    Next t
    
    If Total > 0
      DrawText(wp(0.03), hp(0.03), "N = " + Str(Total), #Red)
    EndIf
    StopDrawing()
    SetGadgetState(7, ImageID(0))
  EndIf
EndProcedure

WinFlags.i = #PB_Window_Invisible|#PB_Window_ScreenCentered|#PB_Window_MinimizeGadget
WinFlags   | #PB_Window_SizeGadget|#PB_Window_MaximizeGadget
If Not OpenWindow(0, 0, 0, 500, 300, "Random Number Distribution", WinFlags)
  MessageRequester("Error", "Window could not be opened.", #MB_ICONHAND) : End
EndIf

CreateImage(0, 1, 1)
ImageGadget(7, 0, 0, 0, 0, ImageID(0))

Running.i = 0
Type.i = 2
Maxx.i = 100
Total = 0
ExitFlag.i = #False

SmartWindowRefresh(0, #True)
WindowBounds(0, 300, 190, #PB_Ignore, #PB_Ignore)

StretchImage()
HideWindow(0, #False)
UpdateImage()

Repeat
  Event.i = WindowEvent();WaitWindowEvent(#Timeout)
  
  While Event
    
    If Event = #PB_Event_CloseWindow
      ExitFlag = #True
    ElseIf Event = #PB_Event_SizeWindow
      StretchImage()
      If Not Running : UpdateImage() : EndIf
    ElseIf Event = #PB_Event_Gadget
      ID.i = EventGadget()
      Select ID
        Case 7
          Select EventType()
            Case #PB_EventType_LeftClick : Running = 1 - Running
            Case #PB_EventType_RightClick
              For t = 0 To #Types - 1
              For i = 0 To #Intervals - 1
                Count(t, i) = 0
              Next i
              Next t
              Maxx = 100
              Total = 0
              If Not Running : UpdateImage() : EndIf
          EndSelect
      EndSelect
    EndIf
    
    Event = WindowEvent()
  Wend
  
  If Running And (Total < 2000000000)
    For t = 0 To #Types - 1
      For a = 1 To #Iterations
      Select t
        Case 1 : x.d = RandomGauss(BellScale)
        Case 2 : x.d = GauRand(100000)/100000.0
        Case 0 : x.d = Blueznl()
        Default : x.d = RandomFlat()
      EndSelect
      Group.i = Int(x * #Intervals)
      
      Count(t, Group) + 1
      If Count(t, Group) >= Maxx
        Maxx = Count(t, Group) + 1
      EndIf
      Next a
    Next t
    Total + #Iterations
  
  EndIf
  UpdateImage()

  Delay(20)

Until ExitFlag

HideWindow(0, #True)
End

;

Re: Flat vs Gaussian Distribution (Demo)

Posted: Fri Dec 11, 2009 11:24 pm
by dell_jockey
Thank you so much for this comparison.

The Blueznl version I like best as it has some interessting properties:
- by using four Random(250) terms, the likelyhood of values ending up in the middle becomes slightly larger, i.e. the tail ends get pinched a bit more.
- by using two Random(500) terms, you end up with virtually a triangular probability distribution.

It's the function that can be most easily tailored to get the distribution you're after.

Thanks again!
All times are UTC+01:00
Page 1 of 1

Powered by phpBB® Forum Software © phpBB Limited