Text Only | Text with Attachments

QB64.org Forum

Active Forums => QB64 Discussion => Topic started by: Joshua on December 29, 2018, 06:56:33 pm

Title: How to generate normal / Gaussian randoms
Post by: Joshua on December 29, 2018, 06:56:33 pm
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.
Title: Re: How to generate normal / Gaussian randoms
Post by: Joshua on December 30, 2018, 08:36:11 pm
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.
Title: Re: How to generate normal / Gaussian randoms
Post by: MWheatley on December 31, 2018, 08:21:02 am
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.  
  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.  
Title: Re: How to generate normal / Gaussian randoms
Post by: Joshua on December 31, 2018, 02:27:16 pm
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?
Title: Re: How to generate normal / Gaussian randoms
Post by: bplus on December 31, 2018, 03:35:24 pm
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.  
  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
  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
  44.  
  45. FUNCTION nFill$ (n, fillLength)
  46.     nFill$ = RIGHT$(STRING$(fillLength, " ") + STR$(n), fillLength)
  48.  
Title: Re: How to generate normal / Gaussian randoms
Post by: Joshua on December 31, 2018, 04:27:45 pm
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.
Title: Re: How to generate normal / Gaussian randoms
Post by: bplus on December 31, 2018, 04:41:05 pm
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.
Title: Re: How to generate normal / Gaussian randoms
Post by: bplus on December 31, 2018, 04:54:37 pm
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.  
  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
  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
  25.  
  26.  
Title: Re: How to generate normal / Gaussian randoms
Post by: MWheatley on January 01, 2019, 11:17:16 am
Quote from: Joshua on December 31, 2018, 02:27:16 pm
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
Text Only | Text with Attachments