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()
