（一）连续随机量的生成-双变量方法

1. 双变量方法的理论

If $X=\left(X_1,X_2\right)$ has a joint density $f_X\left(x_1,x_2\right)$ and $g_{:}\left(x_1,x_2\right)\to \left(y_1,y_2\right)$ is a differentiable transform with the inverse $\left(x_1,x_2\right)=g^{-1}\left(y_1,y_2\right)$. Then the random vector $\left(Y_1,Y_2\right)=g\left(X_1,X_2\right)$ has a density
$$f_Y\left(y_1,y_2\right)=f_X\left(g^{-1}\left(y_1,y_2\right)\right)J\left(y_1,y_2\right)$$
where $J\left(y_1,y_2\right)$ is the Jacobi matrix $J\left(y_1,y_2\right)=\left|\begin{array}{ll}\frac{\partial x_1}{\partial y_1}& \frac{\partial x_2}{\partial y_1}\\ \frac{\partial x_1}{\partial y_2}& \frac{\partial x_2}{\partial y_2}\end{array}\right|$.

Example: Let $U_1\sim U([0,1]),U_2\sim U([0,1])$ be independent, define
$$x_1=\sqrt{-2\log u_1}\cos \left(2\pi u_2\right),x_2=\sqrt{-2\log u_1}\sin \left(2\pi u_2\right)$$

i.e. $\left(x_1,x_2\right)=g\left(u_1,u_2\right)=\left(\sqrt{-2\log u_1}\cos (2\pi u_2),\sqrt{-2\log u_1}\sin \left(2\pi u_2\right)\right)$.
$$\begin{array}{rl}{g}^{-1}\left(x_1,x_2\right)& =\left(\exp \left(-\frac{x_1^2+x_2^2}{2}\right),\frac{1}{2\pi }{\tan }^{-1}\left(\frac{x_2}{x_1}\right)\right),\\ f_U\left(u_1,u_2\right)& =1\text{ for }\left(u_1,u_2\right)\in [0,1]\times [0,1].\\ J\left(x_1,x_2\right)=\left|\begin{array}{ll}\frac{\partial u_1}{\partial x_1}& \frac{\partial u_2}{\partial x_2}\\ \frac{\partial u_1}{\partial x_2}& \frac{\partial u_2}{\partial x_2}\end{array}\right|& =\left|\begin{array}{ll}\frac{\partial g_1^{-1}\left(x_1,x_2\right)}{\partial x_1}& \frac{\partial g_2^{-1}\left(x_1,x_2\right)}{\partial x_1}\\ \frac{\partial g_1^{-1}\left(x_1,x_2\right)}{\partial x_2}& \frac{\partial g_2^{-1}\left(x_1,x_2\right)}{\partial x_2}\end{array}\right|=\frac{1}{2\pi }{e}^{-\frac{x_1^2+x_2^2}{2}}\end{array}$$
Hence $f_X\left(x_1,x_2\right)=\frac{1}{2\pi }\exp \left(-\frac{x_1^2+x_2^2}{2}\right)$.

2. 从$f_X(x)=\frac{1}{C}{e}^{-x^4}$采样的具体操作

In data science, we often come across a distribution whose density is known up to a constant, e.g. a random variable $X$ with a density
$$f_X(x)=\frac{1}{C}{e}^{-x^4}$$
the constant $C={\int }_{-\infty }^{+\infty }{e}^{-x^4}dx$ is not known. Taking
$$f^{\ast }(x)=e^{-x^4},$$
define $C_f=\left\{\left(x_1,x_2\right):0⩽x_1⩽\sqrt{f^{\ast }\left(x_2/x_1\right)}\right\}$. Let $\left(U_1,U_2\right)$ be the random vector uniformly distributed on $C_f$. Then $Y=\frac{U_2}{U_1}$ has a density $\frac{f^{\ast }(y)}{{\int }_{-\infty }^{+\infty }{f}^{\ast }(y)dy}$.

Indeed, let $x=u_1,y=\frac{u_2}{u_1}$, i.e.
$$\begin{array}{c}(x,y)=g\left(u_1,u_2\right)=\left(u_1,\frac{u_2}{u_1}\right).\\ \left(u_1,u_2\right)=g^{-1}(x,y)=(x,xy).\\ J(x,y)=\left|\begin{array}{ll}\frac{\partial g_1^{-1}(x,y)}{\partial x}& \frac{\partial g_2^{-1}(x,y)}{\partial x}\\ \frac{\partial g_2^{-1}(x,y)}{\partial y}& \frac{\partial g_2^{-1}(x,y)}{\partial y}\end{array}\right|=\left|\begin{array}{cc}1& y\\ 0& x\end{array}\right|=x.\end{array}$$
We have $f_{(X,Y)}(x,y)=x\cdot k$, where $k=\frac{1}{\text{ Area }\left(C_f\right)}$, since $f\left(U_1,U_2\right)\left(u_1,u_2\right)=k$ for $\left(u_1,u_2\right)\in C_f$. Now we consider the distribution of $Y$, which is given by
$$\begin{array}{rl}{f}_Y(y)& ={\int }_{-\infty }^{+\infty }{f}_{(X,Y)}(x,y)dx\\ & ={\int }_{-\infty }^{+\infty }x\cdot k{1}_{C_f}dx\\ & ={\int }_0^{\sqrt{f^{\ast }(y)}}kxdx\\ & =k{f}^{\ast }(y)=\frac{f^{\ast }(y)}{{\int }_{-\infty }^{+\infty }{f}^{\ast }(y)dy}.\end{array}$$

3. 用Python编程实现$f_X(x)=\frac{1}{C}{e}^{-\frac{|x{|}^5}{2}}$采样

Make a python or $R$ program to generate 1000 random quantities from the probability distribution with a density $f_X(x)=\frac{1}{C}{e}^{-\frac{|x{|}^5}{2}}$ with $C={\int }_{-\infty }^{+\infty }{e}^{-\frac{|x{|}^5}{2}}dx$.

3.1 Python程序

import numpy as np
import matplotlib.pyplot as plt

# Define the probability density function f*(x)
def f_star(x):
    return np.exp(-np.abs(x)**5 / 2)

# Define the region C_f
def in_region(x1, x2):
    return (0 <= x1) and (x1 <= np.sqrt(f_star(x2 / x1)))

# Generate random quantities for U1 and U2
def generate_random_quantities(n):
    U1 = []
    U2 = []
    while len(U1) < n:
        x1 = np.random.rand()
        x2 = np.random.rand()
        if in_region(x1, x2):
            U1.append(x1)
            U2.append(x2)
    return U1, U2

# Generate random quantities of Y = U2 / U1
def generate_random_Y(U1, U2):
    Y = [u2 / u1 for u1, u2 in zip(U1, U2)]
    return Y

# Number of random quantities to generate
num_samples = 1000

# Generate random quantities for U1 and U2
U1, U2 = generate_random_quantities(num_samples)

# Calculate random quantities for Y = U2 / U1
random_Y = generate_random_Y(U1, U2)

# Draw a histogram
plt.hist(random_Y, bins=30, density=True, alpha=0.6, color='b', edgecolor='black')
plt.xlabel('Random Quantities Y')
plt.ylabel('Density')
plt.title('Histogram of Random Quantities Y')
plt.show()

3.2 Python程序步骤解析

import numpy as np
import matplotlib.pyplot as plt

These two lines import the 'nump’y and ‘matplotlib.pyplot’ libraries, which are used for numerical computations and plotting.

def f_star(x):
    return np.exp(-np.abs(x)**5 / 2)

This is a function ‘f_star(x)’ that defines the probability density function $f^{\ast }(x)$.

def in_region(x1, x2):
    return (0 <= x1) and (x1 <= np.sqrt(f_star(x2 / x1)))

This function ‘in_region(x1, x2)’ checks whether $(x_1,x_2)$ is within the region $C_f$.

def generate_random_quantities(n):
    U1 = []
    U2 = []
    while len(U1) < n:
        x1 = np.random.rand()
        x2 = np.random.rand()
        if in_region(x1, x2):
            U1.append(x1)
            U2.append(x2)
    return U1, U2

This function ‘generate_random_quantities(n)’ generates random $U_1$ and $U_2$ while ensuring they are within the $C_f$ region. It’s a loop that continues until $n$ points satisfying the condition are generated.

def generate_random_Y(U1, U2):
    Y = [u2 / u1 for u1, u2 in zip(U1, U2)]
    return Y

This function ‘generate_random_Y(U1, U2)’ calculates the random quantity $Y=U_2/U_1$ for each $(U_1,U_2)$ pair.

num_samples = 1000

‘num_samples’ specifies how many random quantities should be generated.

U1, U2 = generate_random_quantities(num_samples)

This line generates ‘num_samples’ pairs of $(U_1,U_2)$ that satisfy the condition.

random_Y = generate_random_Y(U1, U2)

This line calculates the random quantity $Y=U_2/U_1$ for each $(U_1,U_2)$ pair.

plt.hist(random_Y, bins=30, density=True, alpha=0.6, color='b', edgecolor='black')
plt.xlabel('Random Quantities Y')
plt.ylabel('Density')
plt.title('Histogram of Random Quantities Y')
plt.show()

This part of the code plots a histogram of the random quantities Y. The ‘hist’ function is used to draw the histogram, ‘xlabel’ sets the x-axis label, ‘ylabel’ sets the y-axis label, ‘title’ sets the title of the plot, and finally, ‘show’ displays the plot.

Lastly, make sure you have installed the ‘numpy’ and ‘matplotlib’ libraries by using the following command:

pip install numpy matplotlib

Running this code will display a histogram of the generated random quantities Y based on the given distribution.