## How to "customize" the random generator?

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

### Re: How to "customize" the random generator?

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 DoubleStatic As Long u2_cachedStatic 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 IfEnd Function Dim As Ulong i, a(0 to 3)Dim As Double x, mean, sdMean = 100 : sd = 5For 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) += 1NextPrint 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 DoubleStatic As Long u2_cachedStatic 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 IfEnd Function Sub QuickSort( a() As Double, lo As long, hi As long)Dim As Long i, jDim 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 IfEnd Sub Dim as Ulong iDim as Double Mean, sd, d(1 to 20) Randomize Mean = 100 : sd = 2 ' Mean as minimum valueFor i = 1 to 20  d(i) = Mean + Gauss*sd ' Note +NextQuickSort(d(), 1, 20)For i = 1 to 20  Print d(i)Nextprint' Mean as maximum valueFor i = 1 to 20  d(i) = Mean - Gauss*sd ' Note -NextQuickSort(d(), 1, 20)For i = 1 to 20  print d(i)NextPrint' Mean as meanFor i = 1 to 20  d(i) = Mean + Perturbation ' Note ± via PerturbationNextQuickSort(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: 2754
Joined: Jun 02, 2015 16:24

### Re: How to "customize" the random generator?

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: 5701
Joined: Jan 10, 2006 20:30
Location: Scotland

### Re: How to "customize" the random generator?

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 rndsFunction f1(x As Double) As Double    If x>0 Then Return x Else Return -xEnd FunctionFunction 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 Functionsub setupDim As Long counterFor x As Double=-1 To 1 Step 2/1000000    counter+=1    Var a=area(@f1,-1,x)    cd(counter)= aNextend subend namespacernds.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 IfEnd SubScreen 20,32,,64Color ,Rgb(0,20,0)ClsDim As Integer xres,yresScreeninfo xres,yresLine(0,.55*yres)-(xres,.6*yres),Rgb(50,50,50),bfLine(0,.6*yres)-(xres,yres),Rgb(0,30,0),bfline(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 IfNextPrint "Done"Sleep `
Richard
Posts: 2907
Joined: Jan 15, 2007 20:44
Location: Australia

### Re: How to "customize" the random generator?

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 IfEnd 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 aboveRandomizeScreen 19Window( -5, -0.02 )-( 5, 1.02 )' target values to test calibrationDim As Double mu = 5.0  ' target meanDim As Double sd = 2.0  ' target standard deviationDim As Double x, y, sum_x, sum_xxDim As Integer i, n = 1e6For 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     '    statisticsNext iPrintPrint Using "    Mean =###.#####"; Sum_x / n    ' report satisticsPrintPrint Using " Std Dev =###.#####"; Sqr( ( ( Sum_xx - sum_x * sum_x / n ) ) / ( n - 1 ) )Printy = 0.02425x = probit( y )Print Using "##.###% fall below ##.###### SD"; y * 100; xSleep'-----------------------------------------------------------------------' end of file'-----------------------------------------------------------------------`