Services & Resources / Wolfram Forums
-----
 /
MathGroup Archive
2003
*January
*February
*March
*April
*May
*June
*July
*August
*September
*October
*November
*December
*Archive Index
*Ask about this page
*Print this page
*Give us feedback
*Sign up for the Wolfram Insider

MathGroup Archive 2003

[Date Index] [Thread Index] [Author Index]

Search the Archive

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


  • Prev by Date: RE: Re: simple equation system -- crashing M5?
  • Next by Date: Re: Where to reply to...
  • Previous by thread: Re: Part 2 of a recent post on Plot and v 5
  • Next by thread: Re: Warning Long Post: Help Speeding up algorithm