Is there a simple error in the answer key, or am I using the wrong approach to get $P(X<0.5)$

Question

I am working on a problem that gives me a joint pdf:
$$f_{x,y}(x,y) = 6xy, 0<x<1, 0<y<\sqrt{x} $$ I am asked to find $P(X < 0.5)$ with three decimal places.

My approach was to integrate: $\int_{0}^{\sqrt{x}} 6xy\ dy = 3x^2 = f_{x}(x)$ to get the marginal pdf of $x$. Then, I integrated again: $\int_{0}^{0.5} 3x^2 dx$ to get $P(X < 0.5)$. What I got was 0.125, but apparently the answer is 0.625.

Am I missing something small or is the answer key just wrong?

Answer 1

The two comments were so insightful and instructive that I thought it was worth developing them further into an answer. This way the question has an answer more than just saying the answer key in the book is wrong.

Sketching

The suggestion from whuber was

Sketch the domain and notice it's very thin when $x<1/2$. Then notice that $6xy$ must be smaller in that thin region than elsewhere, because $x$ and $y$ are both small there. Consequently, $Pr(𝑋<1/2)$ must be substantially less than $1/2$. On that basis alone your answer is plausible and that in the answer key is not.

Using Python we can obtain the following figure.

All the shaded area (both light and dark) is $R_{XY} = \{ (x,y) \in \mathbb{R}^2 | 0 \leq x \leq 1, 0 \leq y \leq \sqrt(x) \}$. The darker grey region is the integration region.

As the comment suggested, we can observe here the integration region is small (less than half the total region). Not only this, but in this dark shaded region $x$ and $y$ are small - so $6xy$ will be smaller over this region than the other. The heatmap below indicates the value of $6xy$, with darker red for values closer to $6$ and paler white closer to $0$.

Given that this is a narrow region and that this is a region where $6xy$ will be smallest, we should be expecting an answer $<0.5$ and so we can rule out the answer from the book which is larger than this, as a mistake.

Monte Carlo Integration

The suggestion from jbowman was

Generate, say, 1 million $U(0,1)$ variates for $X$, the same for $Y$. Discard all pairs for which $y>\sqrt(x)$ . Calculate $f(x,y)$ for the remainder. Sum the probabilities for $f(x,y)$ when $𝑥<0.5$ and for $𝑓(𝑥,𝑦)$ overall; divide, and you'll get a number very close to your $0.125$.

So, first we generate random uniform pairs of points $(x,y) \in [0,1] \times [0,1]$. We plot these points below.

Then we discard all the points for which $y > \sqrt(x)$. The plot is below.

For each of these points $(x,y)$, we calculate the value of $6xy$ and sum. In my case I got $\approx 9936$.

Then we calculate the value of $6xy$ only for the points with $x< 0.5$. These are the points in orange below.

Summing these I get $\approx 1244$. Finally I calculate $1244/9936 \approx 0.125$.

If you are interested in this, you might enjoy these articles. Estimating an integral by using Monte Carlo simulation. and Introductory Examples of Monte Carlo simulations in SAS I think you can still enjoy the articles without programming in SAS and should be able to adapt any example to Python or R.

Python Code

For the sketch and heatmap

import numpy as np
import matplotlib.pyplot as plt
x = np.linspace(0,1,1000)
f, ax = plt.subplots(1)
y= [np.sqrt(xi) for xi in x]
ax.plot(x,y)
plt.axvline(0.5, c="black", linestyle='--')
fill_x=np.linspace(0,0.5,1000)
fill_y=[np.sqrt(xi) for xi in fill_x]
ax.set_ylim(ymin=0, ymax=1)
ax.set_xlim(xmin=0, xmax=1)
ax.set_xticks([0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1])
ax.set_yticks([0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1])
plt.fill_between(x, y, alpha=0.1, color="grey")
plt.fill_between(fill_x, fill_y, alpha=1, color="grey")
plt.grid()

For the heatmap

import numpy as np
import matplotlib.pyplot as plt
x1 = np.linspace(0,1,1000)
f, ax = plt.subplots(1)
y1= [np.sqrt(xi) for xi in x1]
ax.plot(x1,y1)
plt.axvline(0.5, c="black", linestyle='--')
ax.set_ylim(ymin=0, ymax=1)
ax.set_xlim(xmin=0, xmax=1)
ax.set_xticks([0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1])
ax.set_yticks([0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1])
x, y = np.meshgrid(np.linspace(0, 1, 100), np.linspace(0, 1, 100))
z = 6*x*y
z_min, z_max = -np.abs(z).max(), np.abs(z).max()
c = ax.pcolormesh(x, y, z, vmin=0, vmax=6, cmap = 'Reds')
ax.set_title('Heatmap for 6xy')
f.colorbar(c, ax=ax)
plt.fill_between(x1, y1, np.max(y1), color='white', alpha=1)
plt.show()

For the simulation in the second part

import numpy as np
random_points = np.random.uniform(low=0, high=1, size=(10000,2))
keep = random_points[random_points[:,1]<=np.sqrt(random_points[:,0])]
sum_of_all_kept_points = np.sum([6*p[0]*p[1] for p in keep])
sum_when_x_less_than_half = np.sum([6*p[0]*p[1] if p[0]<0.5 else 0 for p in keep])
print(sum_when_x_less_than_half/sum_of_all_kept_points)