Re: How to sample a 2-dim. r.v. with known density function?
- To: mathgroup at smc.vnet.net
- Subject: [mg65274] Re: How to sample a 2-dim. r.v. with known density function?
- From: "Valeri Astanoff" <astanoff at yahoo.fr>
- Date: Thu, 23 Mar 2006 06:58:36 -0500 (EST)
- References: <dvrbl7$a1l$1@smc.vnet.net>
- Sender: owner-wri-mathgroup at wolfram.com
Hi Kees,
Here is a quick and dirty way to do yourself,
for a bivariate non-conventional distribution,
what Mathematica's RandomArray does for conventional distributions :
In[1]:=<<Statistics`
In[2]:=randArray[pdf_, {x_, y_},
(* discretize x : *) {xmin_, xmax_, dx_},
(* discretize y : *) {ymin_, ymax_, dy_},
sampleSize_]:=
Module[{cells, cellSampleSize, cellSample},
cells = Table[{x,y}, {x,xmin,xmax,dx},
{y,ymin,ymax,dy}] // Flatten[#,1]&;
cellSampleSize[xx_, yy_] := sampleSize*
(pdf /. x -> xx+dx/2 /. y -> yy+dy/2)*dx*dy;
cellSample[{xx_, yy_}] :=
Table[{xx+Random[]*dx, yy+Random[]*dy},
{cellSampleSize[xx,yy] // Round}];
Flatten[cellSample /@ cells,1]
];
Check example with a bivariate normal distribution :
In[3]:=distr=MultinormalDistribution[{10,20},{{5,1},{1,5}}];
In[4]:=pdf = PDF[distr,{x,y}];
Choice of discretization ranges and sample size is critical.
Sample size has to be sufficiently large, i.e. far greater
than the number of discretization cells :
In[5]:=samp = randArray[pdf,{x,y},{0,20,1},{10,30,1},10000];
In[6]:=Mean[distr]
Out[6]={10,20}
In[7]:=Mean[samp]
Out[7]={10.0045,20.0012}
In[8]:=StandardDeviation[distr]
Out[8]={Sqrt[5], Sqrt[5]}
In[9]:=StandardDeviation[samp]
Out[9]={2.24731,2.24718}
hth
V.Astanoff