Services & ResourcesWolfram ForumsMathGroup Archive

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

Re: Compile and NormalDistribution

• To: mathgroup at smc.vnet.net
• Subject: [mg15220] Re: Compile and NormalDistribution
• From: Paul Abbott <paul at physics.uwa.edu.au>
• Date: Tue, 22 Dec 1998 04:01:49 -0500
• Organization: University of Western Australia
• References: <75f9jo$7mi@smc.vnet.net>
• Sender: owner-wri-mathgroup at wolfram.com

Andrew Watson wrote:

> In[1]:= <<Statistics`ContinuousDistributions`
> 
> gaussianArray[n_] :=
>         RandomArray[NormalDistribution[0.0, 1.0],{n,n}]  +
>          I  RandomArray[NormalDistribution[0.0, 1.0],{n,n}]
> 
> In[6]:= Timing[tmp=gaussianArray[128] ;]
> 
> Out[6]= {1.26667 Second,Null}

You must have a very fast computer. This is much slower on my PowerMac!

> At 11:11 PM -0800 12/17/98, Brett Patterson wrote:
> >
> >I am trying to speed up the generation of a large array of random
> >complex Gaussian numbers (e.g. 360 x 360).

Some comments. If you look at Statistics`NormalDistribution` you can see
how Random[] and RandomArray[] works for NormalDistribution[]:

normal = Compile[{{mu, _Real}, {sigma, _Real}, 
	{q1, _Real}, {q2, _Real}},
	mu + sigma Sqrt[-2 Log[q1]] Cos[2Pi q2]	]

NormalDistribution/: Random[NormalDistribution[mu_:0, sigma_:1]] :=
	normal[mu, sigma, Random[], Random[]]

normalpair = Compile[{{mu, _Real}, {sigma, _Real}, 
	{q1, _Real}, {q2, _Real}},
	mu + sigma Sqrt[-2 Log[q1]] {Cos[2Pi q2], Sin[2Pi q2]}]

NormalDistribution/: RandomArray[
	NormalDistribution[mu_:0, sigma_:1], dim_] :=
	  Module[{n, array},
		n = If[VectorQ[dim], Apply[Times, dim], dim];
		array = Flatten[Table[normalpair[mu, sigma, Random[], Random[]],
		 {Quotient[n, 2]}]];
    If[OddQ[n],
       AppendTo[array, normal[mu, sigma, Random[], Random[]] ]  ];
    If[VectorQ[dim] && Length[dim] > 1,
       Fold[Partition[#1, #2]&, array, Reverse[Drop[dim, 1]] ],	
       array  ]
  ] /; 
  (IntegerQ[dim] && dim > 0) || VectorQ[dim, (IntegerQ[#] && # > 0)&]

The key points are that:

[1] Random[NormalDistribution[]] is already compiled. [2] normalpair
generates a _pair_ of random numbers.

We can easily modify normalpair to generate complex Gaussian numbers:

In[1]:= gaussian = Compile[{{mu, _Real}, {sigma, _Real}}, 
   mu + E^(2 Pi I Random[]) sigma Sqrt[-2 Log[Random[]]]]

In[2]:= First[Timing[tmp = Table[gaussian[0, 1], {180}, {180}]; ]]
Out[2]= 6.4 Second

Cheers,
	Paul

____________________________________________________________________ 
Paul Abbott                                   Phone: +61-8-9380-2734
Department of Physics                           Fax: +61-8-9380-1014
The University of Western Australia            Nedlands WA  6907       
mailto:paul at physics.uwa.edu.au  AUSTRALIA                       
http://www.physics.uwa.edu.au/~paul

            God IS a weakly left-handed dice player
____________________________________________________________________

• Follow-Ups:
  • Re: Re: Compile and NormalDistribution
    • From: jtischer@col2.telecom.com.co (Cl-Jurgen Tischer)