home *** CD-ROM | disk | FTP | other *** search
/ NetNews Usenet Archive 1993 #3 / NN_1993_3.iso / spool / comp / softsys / matlab / 119 < prev    next >
Encoding:
Text File  |  1993-01-23  |  5.0 KB  |  145 lines
  1. Newsgroups: comp.soft-sys.matlab
  2. Path: sparky!uunet!stanford.edu!nntp.Stanford.EDU!forsythe.stanford.edu!nova!maurer
  3. From: Michael Maurer <maurer@nova.stanford.edu>
  4. Subject: SOURCE: randg.m - generate samples from arbitrary PDF
  5. Message-ID: <maurer.727751195@nova.Stanford.EDU>
  6. Sender: news@leland.Stanford.EDU (Mr News)
  7. Organization: STAR Lab, Stanford University, California USA
  8. Date: 23 Jan 93 01:06:35 GMT
  9. Lines: 134
  10.  
  11. Hoping to start a trend, here is an unsolicited matlab m-file.
  12. RANDG generates samples from an arbitrary probability distribution
  13. function, specified symbolically.  Suggested improvements are welcome.
  14.  
  15. #! /bin/sh
  16. # This is a shell archive, meaning:
  17. # 1. Remove everything above the #! /bin/sh line.
  18. # 2. Save the resulting text in a file.
  19. # 3. Execute the file with /bin/sh (not csh) to create the files:
  20. #    randg.m
  21. # This archive created: Fri Jan 22 16:56:20 1993
  22. export PATH; PATH=/bin:$PATH
  23. echo shar: extracting "'randg.m'" '(3479 characters)'
  24. if test -f 'randg.m'
  25. then
  26.     echo shar: will not over-write existing file "'randg.m'"
  27. else
  28. sed 's/^XX//' << \SHAR_EOF > 'randg.m'
  29. XXfunction x=randg(PDF_,BND_,ICBND_,N_,opts_)
  30. XX
  31. XX%RANDG
  32. XX%    Generates samples from arbitrary probability distributions.
  33. XX%    X=randg(ICDF,N) uses the inverse cumulative distribution
  34. XX%    function ICDF with argument 'y' and result 'x'.  For example,
  35. XX%    ICDF='x=-log(1-y)' is the inverse CDF for 'y=exp(-x)' (x>0),
  36. XX%    which has CDF 'y=1-exp(-x)'.  N is the number of samples.
  37. XX%
  38. XX%    X=randg(PDF,BND,ICBND,N) works when the CDF is not easily inverted.
  39. XX%    PDF is a string like 'yp=exp(-x/2)' that describes the
  40. XX%    probability density function 'yp' as a function of 'x'.  BND is
  41. XX%    another string that defines a function that is everywhere
  42. XX%    not less than PDF, for example 'yb=exp(-x/4)' for x>0.  ICBND 
  43. XX%    is the inverse of the CDF of BND, in this case the inverse of
  44. XX%    'y=1-exp(-x/4)' which is 'x=-4*log(1-y)', defined for 0<=y<=1.
  45. XX%    Note that the names of the variables must be the same as shown
  46. XX%    here (x, y, yp, and yb).  Intermediate variables may have any name
  47. XX%    that does not end with '_' and which is not a built-in M-file
  48. XX%    variable (such as 'nargin').
  49. XX%
  50. XX%    X=randg(PDF,BND,ICBND,N,'plot') will plot the algorithm's progress.
  51. XX%
  52. XX%    Here is a complicated example of the second method, taken from
  53. XX%    Leon-Garcia "Probability and Random Processes for Electrical
  54. XX%    Engineering", p. 174.  It generates samples X from a gamma pdf with
  55. XX%    parameter alpha, 0<alpha<1.
  56. XX%
  57. XX%    PDF = 'alpha=0.33;yp=x.^(alpha-1).*exp(-x)/gamma(alpha);'
  58. XX%
  59. XX%    BND = [    'alpha=0.33;                          ';
  60. XX%        'Ip=find(x>1); In=find(x<=1);         ';
  61. XX%        'yb=x;                                ';
  62. XX%        'yb(In)=x(In).^(alpha-1)/gamma(alpha);';
  63. XX%        'yb(Ip)=exp(-x(Ip))/gamma(alpha);     ']
  64. XX%    ICBND= ['alpha=0.33;                             '
  65. XX%        'e=exp(1);                               '
  66. XX%        'In=find(y<(e/(alpha+e)));               '
  67. XX%        'Ip=find(y>=(e/(alpha+e)));              '
  68. XX%        'x=y;                                    '
  69. XX%        'x(In)=((alpha+e)*y(In)/e).^(1/alpha);   '
  70. XX%        'x(Ip)=-log((alpha+e)*(1-y(Ip))/alpha/e);']
  71. XX
  72. XX%    Michael Maurer, 5 May 1991
  73. XX%    Copyright (c) 1993 by Michael Maurer
  74. XX
  75. XXdist_=rand('dist');
  76. XXrand('uniform');
  77. XXif nargin==2,
  78. XX   ICDF_=PDF_;
  79. XX   N_=BND_;
  80. XX   y=rand(N_,1);
  81. XX   for l_=1:height(ICDF_),
  82. XX      eval(ICDF_(l_,:));                % x=ICDF_(y)
  83. XX   end
  84. XXelse
  85. XX   xf_=zeros(N_,1);
  86. XX   fnd_=0;
  87. XX   eff_=0.99;                % initial guess to efficiency
  88. XX   while fnd_<N_,
  89. XX      Ntry_=fix((N_-fnd_)/eff_+10);        % slightly overestimate misses
  90. XX      y=rand(Ntry_,1);
  91. XX      for l_=1:height(ICBND_),
  92. XX     eval(ICBND_(l_,:));            % x=ICBND_(y)
  93. XX      end
  94. XX      for l_=1:height(BND_),
  95. XX     eval(BND_(l_,:));            % yb=BND_(x)
  96. XX      end
  97. XX      for l_=1:height(PDF_),
  98. XX     eval(PDF_(l_,:));            % yp=PDF_(x)
  99. XX      end
  100. XX      t_=rand(Ntry_,1);
  101. XX      Ifnd_=find(t_.*yb <= yp);
  102. XX      Nfnd_=length(Ifnd_);
  103. XX      Nfnd_=min([N_-fnd_ Nfnd_]);
  104. XX      xf_(fnd_+1:fnd_+Nfnd_)=x(Ifnd_(1:Nfnd_));
  105. XX      fnd_=fnd_+Nfnd_;
  106. XX      eff_=Nfnd_/Ntry_;
  107. XX      
  108. XX      if nargin>4,                % nifty graphical debugging
  109. XX     disp2('Begin loop with ', N_-(fnd_-Nfnd_), ' samples needed.');
  110. XX     disp2('Efficiency this iteration was ', eff_);
  111. XX     if (fnd_-Nfnd_+length(Ifnd_)>N_),
  112. XX        disp2('Computed ', fnd_-Nfnd_+length(Ifnd_)-N_, ' extra samples'); 
  113. XX     end;
  114. XX     hold off;
  115. XX     axis([min(x) max(x) 0 min(4,max(yp))]);
  116. XX     [xx_,ii_]=sort(x);
  117. XX     plot(xx_,yb(ii_));
  118. XX     hold on;
  119. XX     plot(xx_,yp(ii_));
  120. XX     tt_=t_.*yb;
  121. XX     plot(x(Ifnd_),tt_(Ifnd_),'.')
  122. XX     Infnd_=find(tt_>yp);
  123. XX     plot(x(Infnd_),tt_(Infnd_),'o')
  124. XX     hold off;
  125. XX     if fnd_<N_,
  126. XX        disp('Hit a key for next iteration')
  127. XX        pause
  128. XX     end
  129. XX      end
  130. XX   end
  131. XX   x=xf_;
  132. XXend
  133. XXrand(dist_)
  134. SHAR_EOF
  135. if test 3479 -ne "`wc -c < 'randg.m'`"
  136. then
  137.     echo shar: error transmitting "'randg.m'" '(should have been 3479 characters)'
  138. fi
  139. fi # end of overwriting check
  140. #    End of shell archive
  141. exit 0
  142. --
  143. ______________________________________________________________________
  144. Michael Maurer          maurer@nova.stanford.edu        (415) 723-1024
  145.