Warning Long Post: Help Speeding up algorithm

*To*: mathgroup at smc.vnet.net*Subject*: [mg44387] Warning Long Post: Help Speeding up algorithm*From*: Mark Coleman <mark at markscoleman.com>*Date*: Fri, 7 Nov 2003 05:16:46 -0500 (EST)*Sender*: owner-wri-mathgroup at wolfram.com

Greetings, I am working on a problem in econometrics, part of which involves drawing random numbers and random vectors from nonstandard distributions. Based upon previously published results, some efficient acceptance-rejection algorithms have been developed for some of these distributions. I've taken code for one of these algorithms (the simplest one) that was written in the Gauss programming language, and essentially translated "literally" into Mathematica. As this code will be used in some large-scale Monte Carlo problems, I have an interest in making it run as quickly as possible. The translated code is shown below. Acceptance-rejection algorithms repeatedly draw a random number from a standard distribution, and then modify them in some way to correspond to a draw from the nonstandard distribution. The modified draw is then compared to some stopping criterion (usually in the form of satisfying some inequality), If the inequality is satisfied, the procedure returns the modified draw, if not, a new random variable is created. This continues until the inequality is satisfied. The Mathematica code for this is: In[22]:= Clear[rnd]; rnd[n_,a_,b_]:=Module[{d,theta,xmax,rmax,notAccepted,x,r,p,i}, d=b*b+(4.0*a*n); theta=Max[{(-b+Sqrt[d])/2.0},{(-b-Sqrt[d])/2.0}]; xmax=(b+theta)/a; rmax=-n*Log[theta] - (a/2.0)*xmax*xmax + (theta+b)*xmax; notAccepted=True; While[notAccepted, x=Random[GammaDistribution[n,1]]/theta; r=-n*Log[theta] - (a/2.0)*x*x+(theta+b)*x; p=Exp[r-rmax]; notAccepted=p<Random[UniformDistribution[0.0,1.0]]; ]; x ]; Which gives the result of: In[28]:= Table[rnd[1000,0.5,0.5],{i,1,5000}];//Timing Out[28]= {6.14 Second,Null} The usual way of handling this is in an explicit loop. Given that one (obviously) cannot tell in advance how many passes are required to generate a satisfactory random draw, it is not clear to me how to code this more efficiently in Mathematica, ideally without the explicit loop. The code shown above is one of the acceptance-rejection algorithms I have translated, others are much more computationally demanding. For this module, usuallly no more than 5 passes are needed; for other, more complex distributions, sometimes 50 or more passes are needed. I would very much appreciate any assistance from those in the Mathematica community in making this code more efficient. Best regards, Mark