# Change of random variables and check by plot

Question

As a test, I transform a uniform distribution over the unit square.

But when I check the transformed distribution with Monte Carlo, it is wrong.

What went wrong?

Thanks.

Problem

• Random variables $$\vec{X} = (X_1, X_2)$$ follows the uniform distribution over the unit square.

In other words, $$X_{1} \sim U[0, 1]$$ and $$X_{2} \sim U[0, 1]$$. $$X_{1}$$ and $$X_{2}$$ are independent.

The probability density is: $$\rho_{X}(x_{1}, x_{2}) = 1$$

• The transformation $$f(\vec{x}) = \vec{y}$$ is:

$$y_{1} = \text{sigmoid}(x_{1} + x_{2}) = \frac{1}{1 + e^{-(x_1 + x_2)}}$$ $$y_{2} = \text{sigmoid}(x_{1} - x_{2}) = \frac{1}{1 + e^{-(x_1 - x_2)}}$$

Find probability density $$\rho_{Y}( y_{1}, y_{2})$$

Attempt

The inverse transform $$f^{-1}(\vec{y}) = \vec{x}$$ is:

$$x_{1} = (a_{1} + a_{2}) / 2$$ $$x_{2} = (a_{1} - a_{2}) / 2$$

where

$$a_{1} = x_{1} + x_{2} = -\log\left(\frac{1}{y_1} - 1\right) = \log(y_{1}) - \log(1 - y_{1})$$

$$a_{2} = x_{1} - x_{2} = -\log\left(\frac{1}{y_2} - 1\right) = \log(y_{2}) - \log(1 - y_{2})$$

The partial derivatives of $$a_{1}$$ and $$a_{2}$$ are:

$$\frac{\partial a_{1}}{\partial y_{1}} = \frac{1}{y_{1}} - \frac{1}{1 - y_{1}} \cdot -1 = \frac{1}{y_{1}(1 - y_{1})}$$

$$\frac{\partial a_{2}}{\partial y_{2}} = \frac{1}{y_{2}(1 - y_{2})}$$

The Jacobian of the inverse transform $$f^{-1}(\vec{y}) = \vec{x}$$ is:

$$J =\begin{bmatrix} \frac{\partial x_1}{\partial y_1} & \frac{\partial x_1}{\partial y_2} \\ \frac{\partial x_2}{\partial y_1} & \frac{\partial x_2}{\partial y_2} \\ \end{bmatrix} = \begin{bmatrix} \frac{\ 1}{\ 2 y_{1} (1 - y_{1})} & \frac{1}{\ 2 y_{2} (1 - y_{2})} \\ \frac{1}{ 2 y_{1} (1 - y_{1})} & -\frac{1}{\ 2 y_{2} (1 - y_{2})} \\ \end{bmatrix}$$

The determinant of the Jacobian is:

$$\det(J) = -\frac{1}{2y_{1}(1 - y_{1})y_{2}(1 - y_{2})}$$

The probability density of the transformed distribution is:

$$\rho_{Y}(\vec{y}) = \rho_{X}(f^{-1}(\vec{y})) \lvert \det(J(\vec{y})) \rvert = \frac{1}{2y_{1}(1 - y_{1})y_{2}(1 - y_{2})}$$

Check

• Draw samples of $$\vec{X}$$
• Transform samples of $$\vec{X}$$ to samples of $$\vec{Y}$$
• Compute a weighted 2 dimensional histogram of the samples of $$\vec{Y}$$
• The weights are $$1 / \rho_{Y}(\vec{y})$$
• I expect the weighted histogram would reflect an uniform distribution

Result of the program

Contrary to my expectation, the result is not an uniform distribution. What went wrong?

Program:

import numpy as np
import matplotlib.pyplot as plt

def sample_x(n_trials, eps=1e-5):
"""Draw samples from X"""
return np.random.uniform(eps, 1 - eps, size=(2, n_trials))

def transform_samples(x):
"""Transform samples of X to samples of Y"""
n_trials = x.shape[1]
y = np.empty((2, n_trials))
y[0] = x[0] + x[1]
y[1] = x[0] - x[1]
y = 1 / (1 + np.exp(-y))
return y

def cal_det_jac(y):
"""Calculate absolute determinant of the
Jacobian of the inverse transform"""
a = y[0] * (y[0] - 1)
b = y[1] * (y[1] - 1)
return 1 / (2 * a * b)

def make_histogram(grid_size, ndim):
shape = (grid_size,) * ndim
return np.zeros(shape)

def update_historgram(histogram, sample, weight):
"""Update the histogram.
Assumes the sample is within the unit square"""
grid_size = histogram.shape[0]
index = (sample * grid_size).astype(int)
histogram[index] += weight

def main():
# Draw samples from Y
n_trials = 100000
x = sample_x(n_trials)
y = transform_samples(x)

# Calculate weights of histogram
w = 1 / cal_det_jac(y)

# Calculate 2D histogram
grid_size = 100
hist = make_histogram(grid_size, 2)
for i in range(n_trials):
update_historgram(hist, y[:, i], w[i])

# Plot
plt.pcolor(hist)
plt.show()
plt.close()

if __name__ == '__main__':
main()


Nevermind, I fix it. The line in the update_histogram(...) function should be:
histogram[index[0], index[1]] += weight
histogram[index] += weight