RE: Re: generating random number
- To: mathgroup at smc.vnet.net
- Subject: [mg34545] RE: [mg34519] Re: generating random number
- From: "DrBob" <majort at cox-internet.com>
- Date: Mon, 27 May 2002 01:16:39 -0400 (EDT)
- Reply-to: <drbob at bigfoot.com>
- Sender: owner-wri-mathgroup at wolfram.com
You could replace the following (for the last Plot)
sor = ({#1[[1]], #1[[2]]/numCamp} & ) /@
Reverse /@ Transpose[{Range[numCamp], Sort[gen]}];
with something much simpler:
sor = Transpose[{Sort[gen], Range[numCamp]/numCamp}];
That is, why form pairs in the wrong order and then reverse them? And
why not divide by numCamp immediately, when it's simpler?
Here's an example of using Interpolation to approximate the inverse CDF:
n = 10;
samples = Union[gen[[Random[Integer, {1, numCamp}] & /@ Range[n]]]];
n = Length[samples];
invCDF = Interpolation[Transpose[{Range[n]/n, Sort[samples]}],
InterpolationOrder -> 1];
Plot[invCDF[x], {x, 1/n, 1}];
I've reused "gen" here, but "samples" could be whatever sample data you
have available.
Bobby Treat
-----Original Message-----
From: Raf [mailto:r_a_f at yahoo.it]
To: mathgroup at smc.vnet.net
Subject: [mg34545] [mg34519] Re: generating random number
"ester" <Ester.Piedipalumbo at na.infn.it> ha scritto nel messaggio
news:acferi$rmd$1 at smc.vnet.net...
>
> dear Sirs,
> my name is Ester Piedipalumbo and I am a post doc at the phisics dep.
of
> Naples University.
> I would like to know if there is something for generating random
number
> according an arbitary PDF (where arbitrary means that it is not
standard,
> but I know, of cource)
> Thank a lot
> best regrds
> Ester
>
>
>
A solution is Montecarlo method [una soluzione è il metodo Montecarlo]:
a = 1;
myPDF[x_] := x/a^2/E^(x^2/(2*a^2))
Plot[myPDF[x], {x, 0, 5}, PlotRange -> All]
myCDF[x_] = Integrate[myPDF[t], {t, 0, x}]
Plot[myCDF[x], {x, 0, 5}, PlotRange -> All]
inv = x /. Solve[myCDF[x] == y, x][[1]]
Plot[inv, {y, 0, 1}, PlotRange -> All]
numCamp = 1000;
gen = Table[y = Random[]; inv, {numCamp}];
ListPlot[gen]
sor = ({#1[[1]], #1[[2]]/numCamp} & ) /@
Reverse /@ Transpose[{Range[numCamp], Sort[gen]}];
ListPlot[sor]
Saluti,
Raffaele.