MathGroup Archive 1998

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

Search the Archive

Re: rician random number

  • To: mathgroup at smc.vnet.net
  • Subject: [mg15183] Re: [mg15117] rician random number
  • From: "Tomas Garza" <tgarza at mail.internet.com.mx>
  • Date: Fri, 18 Dec 1998 02:11:00 -0500
  • Sender: owner-wri-mathgroup at wolfram.com

"gcseng at yahoo.com" wrote:

>Anyone knows how to generate random number base on the rician
>probability density function.

>The rician pdf is given by:

>f(r)=(r/sigma^2)*exp(1)^(-1*(r^2+A^2)/(2*sigma^2))*Io(r*A/sigma^2)

>where 2*sigma^2 is the average power in the scatter component of the
>fade, and Io(.) is the modified Bessel function of the first kind and
>zeroth order.

I suggest you use an approximate method based on an idea by Paul Abbott:

First, give some constant value to your A and sigma, say, A=sigma=1 (no
loss of generality).

f[r_]=\!\(E\^\(1\/2\ \((\(-1\) - r\^2)\)\)\ r\ BesselI[0, r]\)

and construct a table of values of the corresponding cumulative
distribution function, with as many points as you wish. A plot of the
density function indicates that it is positive on the range 0 to 7, at
most, so that you can take 70 points, i.e. one every tenth of a unit.
The table, say t, is calculated by means of

t=Table[{j/10,NIntegrate[f[r],{r,0,j/10}]},{j,0,70}]

and you get an approximation to the true distribution function through

h=Interpolation[t]

You can check on the goodness of the approximation by plotting a graph
of the derivative of h and comparing it with the original f. You'll be
surprised (anyway, I was!).

Construct a plot of h:

s=Plot[h[x],{x,0,7}]

Extract the plotpoints (which takes advantage of the adaptive code built
into Plot):

w=Cases[s,Line[{x__}]->x,\[Infinity]];

and then numerically compute the inverse function using interpolation:

c=Interpolation[Reverse/@w]

Then if you take a random number uniformly distributed in (0, 1) and
transform it through c, you get a random number which follows the
distribution f, i.e. the rician distribution. A sample of ten rician
random numbers:

Table[c[Random[]],{10}]

{2.38887,2.36907,1.5057,2.28679,1.29854,1.692,0.69509,2.32733,2.09194,2.0422
5}

Good luck,

Tomas Garza
Mexico City



  • Prev by Date: Re: No trace of trace!
  • Next by Date: Re: Together chokes with Prime Modulus > 46337
  • Previous by thread: Re: rician random number
  • Next by thread: Re: Re: rician random number