How to "customize" the random generator?

General FreeBASIC programming questions.
deltarho[1859]
Posts: 4292
Joined: Jan 02, 2017 0:34
Location: UK
Contact:

Re: How to "customize" the random generator?

Post by deltarho[1859] »

The following is in keeping with dodicat's bulge <laugh>

I am using the Gauss function from my PCG32II code but use Rnd instead of pcg. There are a few Gauss/Normal distribution estimations and the one that I use can be found in many places on the internet.

Gauss is used in 'Mean ± Gauss * Standard_Deviation

From areas of the standardized normal distribution tables we have.

80% will be within 'mean ± 1.28155 * standard_deviation
90% will be within 'mean ± 1.64485 * standard_deviation
95% will be within 'mean ± 1.96 * standard_deviation
99% will be within 'mean ± 2.57583 * standard_deviation

With mean = 100 and standard_deviation = 5 we can expect

80% will be within 100 ± 6.41
90% will be within 100 ± 8.22
95% will be within 100 ± 9.8
99% will be within 100 ± 12.9

Here we request 20 million and segregate them into four bands.

Code: Select all

#Include "string.bi"
#define Perturbation (2*(Rnd<0.5)+1)*Gauss*sd
 
Function Gauss As Double
Static As Long u2_cached
Static As Double u1, u2, x1, x2, w
  If u2_cached = -1 Then
    u2_cached = 0
    Function = u2
  Else
    Do
      x1 = Rnd
      x2 = Rnd
      w = x1 * x1 + x2 * x2
    Loop While w >= 1
    w = Sqr( -2 * Log(w)/w )
    u1 = x1 * w : u2 = x2 * w
    u2_cached = -1
    Function = u1
  End If
End Function
 
Dim As Ulong i, a(0 to 3)
Dim As Double x, mean, sd
Mean = 100 : sd = 5
For i = 1 To 20000000
  x = Mean + Perturbation
  If x > 100 - 6.41 And x < 100 + 6.41 Then a(0) += 1
  If x > 100 - 8.22 And x < 100 + 8.22 Then a(1) += 1
  If x > 100 - 9.8 And x < 100 + 9.8 Then a(2) += 1
  If x > 100 - 12.9 And x < 100 + 12.9 Then a(3) += 1
Next
Print Format(a(0)/200000, "##.00");" ";Format(a(1)/200000, "##.00");" "; _
Format(a(2)/200000, "##.00");" ";Format(a(3)/200000, "##.00")
 
Sleep

With one run I get

Code: Select all

80.02 89.98 95.00 99.02
which is pretty close to the theoretical values.

We can use Gauss in a one-sided way where Mean is treated as a minimum value with the return values being less likely the greater they are and as a maximum value with the return values being less likely the smaller they are.

In the following, we sort the return values and note the 'clustering' close to the minimum and maximum values. The last example uses the two-sided approach.

Code: Select all

#define Perturbation (2*(Rnd<0.5)+1)*Gauss*sd
 
Function Gauss As Double
Static As Long u2_cached
Static As Double u1, u2, x1, x2, w
  If u2_cached = -1 Then
    u2_cached = 0
    Function = u2
  Else
    Do
      x1 = Rnd
      x2 = Rnd
      w = x1 * x1 + x2 * x2
    Loop While w >= 1
    w = Sqr( -2 * Log(w)/w )
    u1 = x1 * w : u2 = x2 * w
    u2_cached = -1
    Function = u1
  End If
End Function
 
Sub QuickSort( a() As Double, lo As long, hi As long)
Dim As Long i, j
Dim as Double pivot
  If lo < hi Then
    i = lo : j = hi
    pivot = a( (lo+hi)/2 )
    Do
      While a(i) < pivot And i < hi: i += 1: Wend
      While a(j) > Pivot And j > lo: j -= 1: Wend
      If i <= j Then swap a(i), a(j): i += 1: j -= 1
    Loop While i <= j
    If lo < j Then QuickSort( a(), lo, j )
    If hi > i Then QuickSort( a(), i, hi )
  End If
End Sub
 
Dim as Ulong i
Dim as Double Mean, sd, d(1 to 20)
 
Randomize
 
Mean = 100 : sd = 2
 
' Mean as minimum value
For i = 1 to 20
  d(i) = Mean + Gauss*sd ' Note +
Next
QuickSort(d(), 1, 20)
For i = 1 to 20
  Print d(i)
Next
print
' Mean as maximum value
For i = 1 to 20
  d(i) = Mean - Gauss*sd ' Note -
Next
QuickSort(d(), 1, 20)
For i = 1 to 20
  print d(i)
Next
Print
' Mean as mean
For i = 1 to 20
  d(i) = Mean + Perturbation ' Note ± via Perturbation
Next
QuickSort(d(), 1, 20)
For i = 1 to 20
  Print d(i)
Next
 
Sleep
So far we have looked at code and not what the code could be used for.

My games imagination and graphics expertise are woefully lacking but I think that Gauss could be used in games. The first one that occurred to me was dodicat's snooker game; which I cannot find at the moment. The computer thrashed me with every game. When the computer determines a move we could introduce a perturbation. Small perturbations would still see the computer 'putting balls away' all be they slightly 'off centre'. Every now and again a perturbation would be large enough for the computer to 'miss'. Needless to say, we would need to experiment with the standard deviation so as not to encumber the computer too much. Perhaps the standard deviation could be an application variable.

I am sure you games writers could come up with more ideas than me.
Tourist Trap
Posts: 2958
Joined: Jun 02, 2015 16:24

Re: How to "customize" the random generator?

Post by Tourist Trap »

dodicat wrote:Neat and simple Richard.
I thought firstly that getting a bias random spread from a uniform distribution would be trivial, but the more I look into it the less trivial it gets.
So did I. That's why I really thank all the contributor to the answer. I think it will go somewhere now.
deltarho[1859] wrote: I am sure you games writers could come up with more ideas than me.
I think you already found a major application, simulating a bias, strength or weakness of a computer oponent by tweaking its probabiliies to miss/win. And I can also think of a way to control a population. For instance, I want to populate my map with trees. The uniform distribution doesn't offer much variety.
dodicat
Posts: 7976
Joined: Jan 10, 2006 20:30
Location: Scotland

Re: How to "customize" the random generator?

Post by dodicat »

A wood near you.

Code: Select all

 

redim shared as double cd(1 To 1000001)
#define Intrange(f,l) int(Rnd*(((l)+1)-(f)))+(f)
#define BiasRange(f,l) Int(  cd( Intrange(1,1000001) )*(((l)+1)-(f))) +(f)

namespace rnds
Function f1(x As Double) As Double
    If x>0 Then Return x Else Return -x
End Function

Function Area(fn As Function(x As Double) As Double,L As Double,U As Double,Ordinates As Integer=10) As Double
    Var n=Ordinates
    If n Mod 2=0 Then n=n+1 
    Var h=(U-L)/n
    Var Part1=fn(L+.5*h)
    Var Part2=0.0
    For i As Integer=1 To n-1
        Part1+=fn(L+h*i+h/2)
        Part2+=fn(L+h*i)
    Next i
    Function= (h/6)*(fn(L)+fn(U)+Part1*4+Part2*2) 
End Function

sub setup
Dim As Long counter
For x As Double=-1 To 1 Step 2/1000000
    counter+=1
    Var a=area(@f1,-1,x)
    cd(counter)= a
Next
end sub
end namespace

rnds.setup

'===============================================================


Sub Tree(i As Any Ptr=0,x1 As Single,y1 As Single,size As Single,angle As Single,depth As Single,colb As Ulong=0,colL As Ulong=0)
    Dim  As Single spread,scale,x2,y2
    spread=25
    scale=.76
    #define incircle(cx,cy,radius,x,y) (cx-x)*(cx-x) +(cy-y)*(cy-y)<= radius*radius
    x2=x1-.25*size*Cos(angle*.01745329)
    y2=y1-.25*size*Sin(angle*.01745329)
    Static As Integer count,fx,fy,sz,z
    If count=0 Then  fx=x1:fy=y1:sz=size:z=2^(depth+1)-1
    Line i,(x1,y1)-(x2,y2),colb
    If count=0 Then  fx=x2:fy=y2:sz=size
    count=count+1
    If count>z Then count=0
    If incircle(fx,fy,(.45*sz),x2,y2)=0 Then Circle i,(x2,y2),.01*sz,colL 
    If depth>0 Then
        Tree(i,x2, y2, size * Scale, angle - Spread, depth - 1,colB,colL)
        Tree(i,x2, y2, size * Scale, angle + Spread, depth - 1,colB,colL)
    End If
End Sub


Screen 20,32,,64
Color ,Rgb(0,20,0)
Cls
Dim As Integer xres,yres
Screeninfo xres,yres

Line(0,.55*yres)-(xres,.6*yres),Rgb(50,50,50),bf
Line(0,.6*yres)-(xres,yres),Rgb(0,30,0),bf
line(0,.575*yres)-(xres,.575*yres),rgb(200,200,200),,&b0000000011110000 
For n As Long=1 To 3000
    Var xpos=intrange(0,xres)
    Var ypos=biasrange(0,yres)
    If ypos < .55*yres Or ypos>.95*yres Then
        tree(,xpos,ypos,15,biasrange(70,110),12,Rgb(80,40,0),Rgb(20,100+biasrange(-50,50),0))
    End If
Next

Print "Done"
Sleep
 
Richard
Posts: 3096
Joined: Jan 15, 2007 20:44
Location: Australia

Re: How to "customize" the random generator?

Post by Richard »

This code generates normally distributed random numbers.
Given the mean, mu, and the standard deviation, SD; x = mu + SD * probit( Rnd )
Acklam's approximation of the probit() function used here is accurate to 9 digits.

Code: Select all

'-----------------------------------------------------------------------
' generate random numbers with gaussian = normal distribution
'-----------------------------------------------------------------------
' Uses the Inverse Normal Cumulative Distribution Function.
' AKA; Normal Quantile Function.  AKA; Probit Function.
'-----------------------------------------------------------------------
Function probit( Byval p As Double) As Double   ' Acklam's algorithm
    ' the absolute value of the relative error is less than 1.15 × 10-9
    '     relative error is defined as ( Xapprox - Xexact ) / Xexact
    '--------------------------------------------------------
    ' define the four polynomials with their coefficients
    #define a0  2.506628277459239
    #define a1 -30.66479806614716
    #define a2  138.3577518672690
    #define a3 -275.9285104469687
    #define a4  220.9460984245205
    #define a5 -39.69683028665376
    #define pa ((((( a5 *r+ a4 )*r+ a3 )*r+ a2 )*r+ a1 )*r+ a0 )
    '----------------------------
    #define b1 -13.28068155288572
    #define b2  66.80131188771972
    #define b3 -155.6989798598866
    #define b4  161.5858368580409
    #define b5 -54.47609879822406
    #define pb ((((( b5 *r+ b4 )*r+ b3 )*r+ b2 )*r+ b1 )*r+ 1 )
    '----------------------------
    #define c0  2.938163982698783
    #define c1  4.374664141464968
    #define c2 -2.549732539343734
    #define c3 -2.400758277161838
    #define c4 -3.223964580411365d-1
    #define c5 -7.784894002430293d-3
    #define pc ((((( c5 *q+ c4 )*q+ c3 )*q+ c2 )*q+ c1 )*q+ c0 )
    '----------------------------
    #define d1 3.754408661907416
    #define d2 2.445134137142996
    #define d3 3.224671290700398d-1
    #define d4 7.784695709041462d-3
    #define pd (((( d4 *q+ d3 )*q+ d2 )*q+ d1 )*q+ 1 )
    '----------------------------
    ' partition into three regions
    If p > 0.02425 Then         ' -2 SD
        If p < 0.97575 Then     ' +2 SD
            ' central region, most common = 95.5%
            Dim As Double q = p - 0.5
            Dim As Double r = q * q
            Return q * pa / pb
        Else
            ' upper region, rare = 2.25%
            Dim As Double q = Sqr( -2 * Log( 1 - p ) )
            Return -pc / pd
        End If
    Else
        ' lower region, rare = 2.25%
        Dim As Double q = Sqr( -2 * Log( p ) )
        Return pc / pd
    End If
End Function

'-----------------------------------------------------------------------
' https://en.wikipedia.org/wiki/Probit
' https://en.wikipedia.org/wiki/Inverse_transform_sampling
' https://stackedboxes.org/2017/05/01/acklams-normal-quantile-Function/

'-----------------------------------------------------------------------
' test and demonstrate the above
Randomize
Screen 19
Window( -5, -0.02 )-( 5, 1.02 )

' target values to test calibration
Dim As Double mu = 5.0  ' target mean
Dim As Double sd = 2.0  ' target standard deviation

Dim As Double x, y, sum_x, sum_xx
Dim As Integer i, n = 1e6
For i = 1 To n
    y = Rnd             ' generate random { 0 < y < 1 }
    x = probit( y )     ' normal distribution { -inf to +inf }
    Pset( x, y ), 14    ' plot a point as CDF of normal distribution
    x = mu + sd * x     ' apply the target mean and Standard Deviation
    sum_x += x          ' accumulate
    sum_xx += x * x     '    statistics
Next i

Print
Print Using "    Mean =###.#####"; Sum_x / n    ' report satistics
Print
Print Using " Std Dev =###.#####"; Sqr( ( ( Sum_xx - sum_x * sum_x / n ) ) / ( n - 1 ) )
Print

y = 0.02425
x = probit( y )
Print Using "##.###% fall below ##.###### SD"; y * 100; x

Sleep

'-----------------------------------------------------------------------
' end of file
'-----------------------------------------------------------------------
Post Reply