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