How to generate normal / Gaussian randoms

[Joshua]

Given RND, does anyone have function that can give a fairly decent approximation to normally distributed values?  It doesn't have to be exact, speed of computation is more important.

[Joshua]

For anyone interested, a very approximate normally distributed randoms can be generated by multiplying two random numbers:

NRND = (RND*2-1)*(RND*2-1)

This will give values between -1 and 1, and if you plot them in a histogram, they will be mostly clustered around 0.

Not great, but it's fast, good for introducing a little natural variability to simulations.

[MWheatley]

Twice in a week!  Something I can make a contribution on.

The code below both calculates normally-distributed random numbers, and plots them in a histogram.  I have similar code for a couple of other distributions (Poisson, etc) as well as the original academic citations, somewhere.

I also have code for the t-distribution, for assessing statistical significance.

Code: QB64: [Select]

1. ' Program to generate random numbers from a Normal distribution
2. ' Program developed in early 1980s in Commodore Basic from published listing
3. ' Moved to Liberty Basic in November 2007
4. ' Converted to QuickBasic 4.5 and QB64 in December 2014
5. ' Edit values in lines 510 and 515 to change the distribution's SD and mean
6.  
7.  
8.  
9. CLS : PRINT
10. 10 Z = 3000: REM NUMBER OF RANDOM NUMBERS
11. 20 DIM A(Z): DIM B(20)
12. 30 GOSUB 500
13. 40 L = .05
14. 50 FOR J = 0 TO 1 STEP L
15. 60 FOR I = 1 TO Z
16. 70 IF A(I) <= J + L AND A(I) > J THEN B(((J * 2) * 10) + 1) = B(((J * 2) * 10) + 1) + 1
17. 80 NEXT I
18. 90 NEXT J
19. 100 FOR I = 1 TO 20
20. 110 IF B(I) > BMAX THEN BMAX = B(I)
21. 120 NEXT I
22. 130 L = 1
23. 140 IF BMAX <= 75 THEN 175
24. 145 BTEMP = BMAX
25. 150 BMAX = (BMAX / L + 1): IF BMAX <= 75 THEN 170
26. 155 BMAX = BTEMP
27. 160 L = L + 1: GOTO 150
28. 170 FOR I = 1 TO 20: B(I) = INT(B(I) / L): NEXT I
29. 175 CLS
30. 180 FOR I = 1 TO 20
31. 185 IF B(I) = 0 THEN PRINT : GOTO 200
32. 190 FOR J = 1 TO B(I): PRINT "Q"; : NEXT J: PRINT
33. 200 NEXT I
34. 205 PRINT "EACH BLOCK EQUALS"; L; "RANDOM NUMBERS"
35. 207 PRINT "FROM A NORMAL DISTRIBUTION"
36. 210 END
37. 500 REM NORMAL DISTRIBUTION GENERATOR
38. 505 N = VAL(MID$(TIME$, 7, 2)): RANDOMIZE N: REM RANDOM SEED
39. 510 EX = .5: REM DISTRIBUTION MEAN
40. 515 SD = .25: REM DISTRIBUTION STANDARD DEVIATION
41. 520 FOR I = 1 TO Z
42. 530 FOR J = 1 TO 12
43. 540 R = RND(1)
44. 550 SUM = SUM + R
45. 560 NEXT J
46. 570 A(I) = SD * (SUM - 6) + EX
47. 580 SUM = 0
48. 590 NEXT I
49. 595 RETURN
50.  

[Joshua]

Looks like it was written in GWbasic!  Would you be able to put the random number generator code portion and put it in a function?

[bplus]

Interesting subject, I myself have wanted "normal" random distributions about a point when graphing.

Here is a function that lets you pick the level of scatter about a point. Level 1 is an even distribution on a range around X, 2 is close to the above codes, a parabolic distribution and levels beyond 2 are powers graphs...

Code: QB64: [Select]

1. _TITLE "nRnd test" 'B+ 2018-12-31
2. 'test a function to return a random number around X with +- range
3. ' is it like a normal distribution?
4. 'more like a power series
5.  
6. RANDOMIZE TIMER
7. 'cont = -1
8. 'WHILE cont
9. '    CLS
10. '    sr = 1
11. '    r = RND
12. '    FOR i = 1 TO 10
13. '        sr = sr * r
14. '        PRINT sr
15. '    NEXT
16. '    INPUT "press enter to repeat, press any + enter to quit..."; wate$
17. '    IF LEN(wate$) THEN EXIT WHILE
18. 'WEND
19. ''OK keep multiplying a rnd number and it gets smaller
20.  
21. DIM dat(20)
22. FOR l = 1 TO 6
23.     FOR i = 1 TO 1000
24.         'add .5 for rounding
25.         r = INT(nRand(10, 5, l) + .5) 'expect nRand to return 5 to 15, so integer "categories" 0 to 20 should
26.         dat(r) = dat(r) + 1
27.     NEXT
28.  
29.     CLS
30.     FOR i = 0 TO 20
31.         PRINT STRING$(INT(dat(i) / 10), "X");
32.         LOCATE i + 1, 40: PRINT nFill$(i, 3); nFill$(dat(i), 4); nFill$(INT(100 * dat(i) / 1000), 3); "%"
33.     NEXT
34.     PRINT "nRand(10, 5,"; l; ")  dat(9) + dat(10) + dat(11) = "; INT(100 * (dat(9) + dat(10) + dat(11)) / 1000); "%"
35.     INPUT "press enter to continue or end..."; wate$
36.     ERASE dat
37. NEXT
38.  
39.  
40. FUNCTION nRand (aroundX, pmRange, level)
41.     IF RND < .5 THEN r = -pmRange * RND ^ level ELSE r = pmRange * RND ^ level
42.     nRand = aroundX + r
43. END FUNCTION
44.  
45. FUNCTION nFill$ (n, fillLength)
46.     nFill$ = RIGHT$(STRING$(fillLength, " ") + STR$(n), fillLength)
47. END FUNCTION
48.  

[Joshua]

I'm not sure it's behaving correctly for me.  When I hit enter, the category 10 has a big spike that goes off screen.  Or how should I read the graph?   Maybe I'm doing it wrong.

[bplus]

Read it sideways, the function is finding random numbers around 10 with a plus minus (pm) range of 5.

The data is rounded to integer intervals so 5 is first showing hits and 15 is last showing hits.

It is a frequency distribution most frequent around 10 when level is above 1.

As level increases, the more and more the points are closer to 10.

There is a calculation to see the percent of hits for middle 3 intervals. At level = 1, it is an even distribution so 30% hits, as level increases this percent does spike past 80% heading off screen around 10.

[bplus]

Here is graphic demo using function for drawing points randomly on line across screen, level 1 is even distribution, level 2 gathers points around the center x = 300 with range +-300...

Code: QB64: [Select]

1. _TITLE "nRnd demo" 'B+ 2018-12-31
2. 'test a function to return a random number around X with +- range
3. ' is it like a normal distribution?
4. 'more like a power series
5.  
6. RANDOMIZE TIMER
7. SCREEN _NEWIMAGE(600, 600, 32)
8.  
9. FOR l = 1 TO 3
10.     CLS
11.     FOR n = 1 TO 50
12.         FOR y = 0 TO 600
13.             x = nRand(300, 300, l)
14.             PSET (x, y), _RGB32(255, 255, 255)
15.         NEXT
16.     NEXT
17.     PRINT "Level "; l; " distribution:"
18.     _DELAY 2
19. NEXT
20.  
21. FUNCTION nRand (aroundX, pmRange, level)
22.     IF RND < .5 THEN r = -pmRange * RND ^ level ELSE r = pmRange * RND ^ level
23.     nRand = aroundX + r
24. END FUNCTION
25.  
26.  

[MWheatley]

Looks like it was written in GWbasic!  Would you be able to put the random number generator code portion and put it in a function?

It will run in QB45 and QB64 just fine.  I'm not sure I'm up to putting it into a function (I very rarely use them, preferring subroutines), but the relevant code is the lines beginning with 500.

Malcolm