Simulation in R to check graphically that marginal distributions are correct

Question

The distribution on $R^2$ with joint density $h$ with respect to the Lebesgue measure is: $$h(x,y)=\frac{3}{2}y 1_{A}(x,y), \ \ A=\{(x,y) \in R^2|0<y, x^2+y^2<1\}.$$ Then I have found the marginal density $f_X(x)=\frac{3}{4}(1-x^2)$

And therefore we get that the conditional distribution of Y given X is: $$f(Y|X)=\frac{h(x,y)}{F_X(x)}=-\frac{2y}{x^2-1}$$

Now I have to use these results to simulate outcomes from the distribution of $(X, Y)$, and check graphically that the marginal distributions are correct. I want to do that in R. But I'm not sure how to do that. I think I have to draw random outcomes $y_i$ from the marginal distribution of $X$. Then for each $x_i$, I think I have to draw an outcome from $P_{x_i}$ (the conditional distribution of Y given X). But how can I do that simulation? I think simulation from a continuous distribution with CDF $F$ can be generated as $F^{-1}(U)$ where $U$ is uniform on $(0,1).$ Is that a correct understanding? Can anyone help me how to do that simulation in R?

Answer 1

Your value for the CDF of $X$ is incorrect, but perhaps this is just bad notation on your part. What you have written is in fact the PDF of $X$.

Simulate $U_1$, a standard uniform random variable. Plug this value into the inverse CDF $F^{-1}_X(U_1)$ of the marginal density of $X$. (This inverse may not be tractable, although it exists, and the Acceptance-Rejection Algorithm may be better [Try Uniform(-1,1)]). Save this value. Next simulate another standard uniform random variable ($U_2$). Plug this value and the saved $X$ value into the inverse CDF $F^{-1}_{Y|X} (U_2)$. Now you have a random sample from the joint distribution $(X,Y)$. Repeat however many times.

To do this in R, you just have to code the functions and use the ${\tt runif}$ function.

Answer 2

What you need is a way of checking your mathematical analysis by means of a simulation that relies as little as possible on that analysis. Otherwise--if you base the code on your calculations--the code will give results that are only as correct as the calculations themselves. It wouldn't be much of a check at all.

In this instance, rejection sampling is an attractive prospect because the problem formulation makes this method very simple to implement and it's not too inefficient. (As we will see, only about two-thirds of the points will be rejected.) The idea is that projecting the uniform distribution of points $(x,y,z)$ where $0 \le z \le h(x,y)$ onto the $(x,y)$ plane (obviously) produces the given distribution, because the density of such points per unit area is--by construction--proportional to the range of possible $z$ values, $h(x,y)-0 = h(x,y).$

The only preliminary analysis needed is to note that the definition of the region $\mathcal A$ implies $-1\le x \le 1$ and $0\le y\le 1.$ Moreover, because $0\le y\le 1,$ it is obvious that $h(x,y) \le 3/2.$ Thus, the set of points $\mathcal E$ in question is wholly contained within the box $[-1,1]\times [0,1]\times [0,3/2]\subset\mathbb{R}^3.$ (It is a portion of a vertical cylinder.)

Rejection sampling generates uniformly random points within the box and just throws out any such points not in $\mathcal E.$

Here is the immediate solution in R:

n <- 1e4
X <- data.frame(x = runif(n, -1, 1), y = runif(n, 0, 1), z = runif(n, 0, 3/2))
i <- with(X, 0 < y & x^2 + y^2 < 1 & z <= (3/2)*y) # Indexes of points to keep
X <- X[i, ]

The first line specifies how many points, n, to generate within the box. The second line generates those points, component by component, as rows of a data frame. The third line identifies which points to keep, *simply by applying the formula for $h.$ The fourth line retains just those points in the data.frame object X. ```

In this case, 3352 points of the original 10000 were retained, which is large enough to give precise indications of all three distributions. (You can reproduce this simulation by preceding it with set.seed(17).)

You will want to examine the results. The realization of $(x,y)$ can be shown with its scatterplot (left panel) while their marginal distributions can be displayed with histograms (middle and right panels). The basic R commands are

with(X, {plot(x,y); hist(x, freq=FALSE); hist(y, freq=FALSE)})

Because you are using this to check your density calculations, you can do so by superimposing graphs of the computed marginal densities on the histograms. I have used your expression $f_x(x) = (3/4)(1-x^2)$ in the second panel to check the marginal density of $x.$ (You can add another call to curve to superimpose your computed expression for $f_y$ on the histogram for $y$ in the right panel.) The match is good: the slight deviations between the (middle) histogram and the graph can be attributed to random variation in the simulation. (Check this with a chi-squared test if you like.)

Appendix 1

You can use such results to approximate the marginal distributions, too. Given a value of $x,$ choose all $y$ values from points $(x^\prime,y)$ where $x^\prime$ is close to $x.$ The result is a mixture of conditional distributions in a narrow strip. Except near the extremes ($x$ near $\pm 1$ or $y$ near $1$) where the conditional distributions change rapidly, these mixtures will be reasonably close and will serve as an excellent check.

For instance, here are plots of the conditional $y$ distribution for four selected values of $x.$ They are based on the same simulation, but this time with a million points drawn (before rejection).

On these I have superimposed the graphs of $f_{y\mid x}(y) = -y/(x^2-1)I_{[0\le y\le \sqrt{1-x^2}]}(y).$ They check out nicely.

In the same manner, the conditional distributions of $x$ given $y$ can be approximated and plotted:

Because you have not supplied a formula for these densities, I haven't graphed anything on them for comparison. But the pattern is suggestive, isn't it? One would guess these are uniform distributions on an interval that shrinks as $y$ grows large. Referring to the left panel in the first figure will show you exactly how it shrinks.

Appendix 2

To create the first figure, the basic R plotting commands were modified to decorate the plots, overplot the graph of $f_x,$ and make them fit in a single graphic:

par(mfrow=c(1,3))
with(X, {
  plot(x, y, asp=1, cex=3/4, col="#00000040", main="Joint Distribution")
  hist(x, freq=FALSE, breaks=25, main="Marginal Distribution of x")
  curve((3/4)*(1-x^2), add=TRUE, lwd=2)
  hist(y, freq=FALSE, breaks=25, main="Marginal Distribution of y")
})
par(mfrow=c(1,1))

These results (for larger values of n) were used to create the second and third figures.

Answer 3

One possible approach for the requested graphical results:

• Start by simulating points in the rectangle with vertices $(0,0)$ and $(1,1)$ so that $X\sim\mathsf{UNIF}(0,1)$ and, independently, $Y\sim\mathsf{Beta}(2,1).$

• Then throw away points not in $A.$

• Then plot histograms of the remaining $X$'s and $Y$'s, which will be negatively correlated.

Simulation in R:

set.seed(1129)
y1 = rbeta(10^4, 2, 1)
x1 = runif(10^4)
A = (x1^2+y1^2 < 1)
x = x1[A]; y = y1[A]  # keep points in A

summary(x); length(x); sd(x)
     Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
0.0001593 0.1693985 0.3445088 0.3733047 0.5529737 0.9886407 
[1] 6579              # number of X's
[1] 0.2428481         # SD of X's
summary(y);  length(y);  sd(y)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
0.01055 0.41670 0.61122 0.59068 0.77998 0.99884 
[1] 6579
[1] 0.2326322
cor(x,y)
[1] -0.3749655        # Negative correlation

Code for figure:

par(mfrow=c(1,3))
 plot(x,y, pch=".")
 hist(x, prob=T, col="skyblue2")
 hist(y, prob=T, col="skyblue2")
par(mfrow=c(1,1))

Note: You will get better approximations for marginal means and standard deviations, and for the correlation if you start with a million points. But then the scatterplot is to conjested to read. So I used only $10\,000$ initial points, of which $6579$ are in set $A.$