home *** CD-ROM | disk | FTP | other *** search
/ CD Actual 9 / CDACTUAL9.iso / share / Dos / VARIOS / pascal / MATH.SWG / 0001_Gaussian Distribution.pas next >
Encoding:
Pascal/Delphi Source File  |  1996-02-21  |  2.2 KB  |  73 lines
  1. { From : Dr John Stockton (JRS@merlyn.demon.co.uk)
  2.  
  3.   Function RandomGaussianSD1 returns a number drawn from
  4.   a Gaussian distribution with mean of zero, unit standard
  5.   deviation, cutoff at +/- 4 ; I have taken it from another,
  6.   working, program.
  7.  
  8.   It may very easily be amended for any distribution
  9.   whatsoever, provided that the probability density may
  10.   be sufficiently well represented by a known expression
  11.   within a finite rectangular box;
  12.   see RandomDistrib, which is but slightly tested.
  13.  
  14.   In each case, the probable time taken varies inversely with
  15.   the fraction of the box area which is under the PD line. }
  16.  
  17.  
  18. program Distribs ;
  19.  
  20.  
  21. function RandomGaussianSD1 : extended
  22.   { Gaussian distribution, SD = 1, cutoff @ +/- 4.0 } ;
  23. var X : extended ;
  24. begin
  25.   repeat X := (Random-0.5)*8.0 ;
  26.     until Exp(-Sqr(X)*0.5)>Random ;
  27.   RandomGaussianSD1 := X end {RandomGaussianSD1} ;
  28.  
  29. procedure GaussTest ;
  30. var j : byte ; s, t : extended ; const k = 80 ;
  31. begin t := 0.0 ;
  32.   for j := 1 to k do begin
  33.     s := RandomGaussianSD1 ; t := t + Sqr(s) ;
  34.     Write(s:9:3) ; if (j and 7)=0 then Writeln ;
  35.     end ;
  36.   Writeln('RMS: ', Sqrt(t/k):10:3, '  ~ 1 ?  '^G) ;
  37.   Write('That should look Gaussian, mean 0, RMS 1 !') ;
  38.   end {GaussTest} ;
  39.  
  40.  
  41. { The following more general form is less tested but should be OK : }
  42.  
  43. type func = function (X : extended) : extended
  44.   { Normalised so that 0.0 <= func <= 1.0 } ;
  45.  
  46. function RandomDistrib(Min, Max : extended ; Fn : func) : extended
  47.   { Fn distribution from Min to Max } ;
  48. var X : extended ;
  49. begin repeat X := Random*(Max-Min) + Min until Fn(X)>Random ;
  50.   RandomDistrib := X end {RandomDistrib} ;
  51.  
  52. function Gauss(X : extended) : extended ; FAR ;
  53. begin Gauss := Exp(-Sqr(X)*2.0) end {Gauss} ;
  54.  
  55. procedure GenTest ;
  56. var j : byte ; s, t : extended ; const k = 80 ;
  57. begin t := 0.0 ;
  58.   for j := 1 to k do begin
  59.     s := RandomDistrib(-6.0, +6.0, Gauss) ; t := t + Sqr(s) ;
  60.     Write(s:9:3) ; if (j and 7)=0 then Writeln ;
  61.     end ;
  62.   Writeln('RMS: ', Sqrt(t/k):10:3, '  ~ 1 ?  '^G) ;
  63.   Write('That should look Gaussian, mean 0, RMS 0.5 !') ;
  64.   end {GenTest} ;
  65.  
  66.  
  67. BEGIN ;
  68. Writeln('Distributions :') ;
  69. Randomize ;
  70. GaussTest ; Readln ;
  71. GenTest ; Readln ;
  72. END.  E&OE.
  73.