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