0
$\begingroup$

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

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()
$\endgroup$
0
$\begingroup$

Nevermind, I fix it. The line in the update_histogram(...) function should be:

histogram[index[0], index[1]] += weight

instead of

histogram[index] += weight

Result:

enter image description here

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.