MathGroup Archive 2005

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

Search the Archive

Re: Random Normal deviates within compiled function?

  • To: mathgroup at smc.vnet.net
  • Subject: [mg62598] Re: Random Normal deviates within compiled function?
  • From: Lee Newman <leenewm at umich.edu>
  • Date: Tue, 29 Nov 2005 04:44:52 -0500 (EST)
  • References: <dln17e$gf4$1@smc.vnet.net>
  • Sender: owner-wri-mathgroup at wolfram.com

I use the following two functions. The first generates a single value 
from a normal distribution with mean and stdev specified in the 
arguments.  It avoids the use of Mathematica's built in distribution 
functions.   It will generate 10^6 values (using Table or Do) with 
timing of about 4 seconds on my machine (2.8GHz P4 with 1GB RAM).   The 
second function is really fast.  It uses the built-in NormalDistribution 
function, and for some reason using RandomArray makes it quite fast.  It 
generates an array of 10^6 values in about 0.5seconds on my machine.


GaussianNoise = Compile[{{mean, _Real}, {stdev, _Real}},
       Module[{r = 0.0, x = 0.0, y = 0.0},
         While[r == 0.0 || r > 1,
           x = Random[Real, {-1, 1}];
           y = Random[Real, {-1, 1}];
           r = x*x + y*y;
           ];
         (x*Sqrt[-2Log[r]/r])*stdev + mean
         ]
       ];

GaussianArray =
     RandomArray[NormalDistribution[#1, #2], #3] &;
(* the first argument is mean, second is the stdev and the third 
specifies the length of the array. *)

Cheers,
Lee Newman



Gareth Russell wrote:
> Hi All,
> 
> I have a function which needs to be quick as it's used in simulations, 
> and it requires the internal generation of a large vector of random 
> normal deviates (up to 50,000). It appears that a function like 
> Random[NormalDistribution[0,s]] cannot be compiled. Can anyone suggest 
> an algorithm for getting such numbers that would work within a Compile 
> statement and still be quicker than using the non-compiled function?
> 
> The non-compiled function would look as follows:
> 
> f[v_,r_,k_,s_,q_]:=Select[v*Exp[r*(1 - v/k) +
>        RandomArray[NormalDistribution[0, s],
>           Length[v]]], # > q &]
> 
> where v is a vector of 50,000 reals, and r, k, s and q are scalar reals.
> 
> (If anyone is interested, this is for population projection in 
> population viability analysis.)
> 
> Thanks,
> 
> Gareth Russell
> NJIT
> 


-- 
- Lee

____________________________________________________
L E E   N E W M A N               www.leenewman.org
     d o c t o r a l   s t u d e n t
          c o g n i t i v e   p s y c h o l o g y
               c o m p u t e r   s c i e n c e
                    u .  m i c h i g a n
                         a n n   a r b o r
____________________________________________________

           Skylarks singing-
           the farmer
           makes a pillow of his hoe.
			
                                         ISSA


  • Prev by Date: Re: Re: Avoiding divide by zero error in NDSolve
  • Next by Date: function of a function
  • Previous by thread: Re: Random Normal deviates within compiled function?
  • Next by thread: Re: Random Normal deviates within compiled function?