Solved – Sampling from a continuous 2 dimensional probability distribution function for importance sampling

importance-samplingquantilesrandom-generationsamplingsimulation

I just want to clarify a few points with regarding to sampling from a continuous 2-dimensional probability density function. If I want to sample from this pdf, I could sample from a 1D pdf, $P(x)$, along $x$ and one along $y$ using the inverse transform sampling method. This method starts with cumulatively integrating a 1D pdf to calculate the cumulative density function,

$$ C(x) = \int^{x}_{-\infty} \frac{1}{A}P(x') \ dx' $$

Once the cdf has been calculated, it can be inverted to produce the ppf function. Once that has been done, the last step is to sample uniformly in [0,1), then pass those values through the ppf function which redistributes the initially uniformed sampled points in accordance to the original pdf.

In the case of the 2D pdf, I can repeat this process along $x$ and along $y$ to get samples for the 2D pdf. Now, if I want to calculate an integral via importance sampling MC. When I sample in accordance with the 2D pdf for $N$ samples, do I sample along $x$ and $y$ with $N$ samples each so I would have an array of (N,2) or do I sample along $x$ and $y$ with $N$ samples and calculate a $N$ by $N$ mesh of all combinations of $x$ and $y$?

Best Answer

Actually, the inverse cdf method only works in dimension one as $(F_X(X),F_Y(Y))$ is not a Uniform $\mathcal U(0,1)^2$ variate, unless $X$ and $Y$ are independent. To simulate from an arbitrary bivariate distribution with joint pdf $p_{X,Y}(x,y)$, one (among many) possible solution is to

Generate $X$ from the marginal $p_X(x)=\int_{\mathcal Y} p_{X,Y}(x,y)\,\text{d}y$, possibly by inverse cdf means
Generate $Y$ given $X=x$ from the conditional $p_{Y|X}(y|x)=p_{X,Y}(x,y)/p_X(x)$, again possibly by inverse cdf means

As an example, consider a joint distribution defined by the density $$p(x,y)=\frac{x^y}{y!}\exp\{-2x\}\mathbb I_{x>0}\mathbb I_{\mathbb N}(y)$$Then $$p_{Y|X}(y|x)\propto \frac{x^y}{y!} \mathbb I_{\mathbb N}(y)$$ meaning that $Y$ conditional on $X=x$ is distributed as a Poisson $\mathcal P(x)$ variate. And $$p_X(x)=\dfrac{p(x,y)}{p_{Y|X}(y|x)}=\dfrac{\frac{x^y}{y!}\exp\{-2x\}}{\frac{x^y}{y!}\exp\{-x\}}=\exp\{-x\}$$ implies that the marginal distribution of $X$ is an Exponential $\mathcal E(1)$ distribution. Hence simulating from the joint distribution can be done by

generate an Exponential $\mathcal E(1)$ variate, $X=-\log(U_1)$
generate a Poisson $\mathcal P(x)$ variate, $$Y=\sum_{i=1}^\infty \mathbb I_{U_2\le\sum_{j=0}^i \frac{x^i}{i!}\exp\{-x\}}$$

Related Solutions

Solved – Bootstrapping importance sampling estimates

Computing part (2) seems like straight forward probability question. No need to involve statistics. As far as I understand, your question is how to compute the probabilities: $$P(mean_A^N \geq 0) \text{ and } P(mean_B^N \geq 0) $$ Let's compute the first one first, since it's easier. Note the the distribution over all $N$ values of $v_i$ will be a multinomial distribution. For ease of notation, call $p_S^+ = P(v_i = 1)$, $p_S^0 = P(v_i = 0)$, and $p_S^- = P(v_i = -1)$. \begin{eqnarray} P(mean_A^N \geq 0) &=& P\left(\frac{1}{N}\sum^N_{i=1}v_i \geq 0\right) \nonumber \\ &=& P\left(\#(v_i = 1) \geq \#(v_i = -1)\right) \nonumber \\ &=& \sum_{x=0}^N \sum^{\min(x,N-x)}_{y=0}\frac{N!}{x!(N-x-y)!y!}(p_S^+)^x(p_S^0)^{N-x-y}(p_S^-)^y \end{eqnarray} and similarly, we compute: \begin{eqnarray} P(mean_B^N \geq 0) &=& P\left(\frac{1}{N}\sum^N_{i=1}v_i\frac{P_T(v_i)}{P_S(v_i)} \geq 0\right) \nonumber \\ &=& P\left(\frac{p_T^+}{p_S^+}\#(v_i = 1) \geq \frac{p_T^-}{p_S^-}\#(v_i = -1)\right) \nonumber \\ &=& P\left(\frac{p_S^-p_T^+}{p_T^-p_S^+}\#(v_i = 1) \geq \#(v_i = -1)\right) \nonumber \\ &=& \sum_{x=0}^N \sum^{\min\left(\left\lfloor\frac{p_S^-p_T^+}{p_T^-p_S^+}x\right\rfloor,N-x\right)}_{y=0}\frac{N!}{x!(N-x-y)!y!}(p_S^+)^x(p_S^0)^{N-x-y}(p_S^-)^y \end{eqnarray}

distributions – Sampling from a Skew Normal Distribution: Methods and Applications

The most direct way of simulating a random variable from a distribution with cdf $F$ is to first simulate a Uniform variate $U\sim\mathcal{U}(0,1)$ and second return the inverse cdf transform $F^{-1}(U)$. When the inverse $F^{-1}$ is not available in closed form, a numerical inversion can be used. Numerical inversion may however be costly, especially in the tails.

One can also use accept-reject algorithms when the density $f$ is available and dominated by another density $g$, i.e., that there exists a constant $M$ such that$$f(x)<M g(x)$$. In the case of the skew-Normal distribution, the density is$$f(x)=2\varphi(x)\Phi(\alpha x)$$when $\varphi$ and $\Phi$ are the pdf and cdf of the standard Normal distribution, respectively. (Adding a location and a scale parameter does not modify the algorithm, since the outcome simply needs to be rescaled and translated.)

This density seems ideally suited for accept-reject since$$2\varphi(x)\Phi(\alpha x)< 2\varphi(x)$$as $\Phi$ is a cdf. This inequality implies that a first option to run accept-reject is to pick the Normal pdf for $g$ and $M=2$. This works out, as shown by the following picture when $\alpha=-3$:

and leads to an algorithm of the kind

T=1e3 #number of simulations
x=NULL
while (length(x)<T){
   y=rnorm(2*T)
   x=c(x,y[runif(2*T)<pnorm(alpha*y)])}
x=x[1:T]

which returns a reasonable fit of the pdf by the histogram:

There are however transforms of standard distributions that result in skew Normal variates: If $X_1,X_2$ are iid $\mathcal{N}(0,1)$, then

$$\dfrac{\alpha|X_1|+X_2}{\sqrt{1+\alpha^2}}$$
$$\dfrac{1+\alpha}{\sqrt{2(1+\alpha^2)}}\max\{X_1,X_2\}+\dfrac{1-\alpha}{\sqrt{2(1+\alpha^2)}}\min\{X_1,X_2\}$$

are skew Normal variates with parameter $\alpha$. (Both representations are identical when considering that $(X_1+X_2,X_1-X_2)/\sqrt{2}$ is an iid $\mathcal{N}(0,1)$ pair.)

Best Answer

Related Solutions

Solved – Bootstrapping importance sampling estimates

distributions – Sampling from a Skew Normal Distribution: Methods and Applications

Related Question