概率论08 随机变量的函数

作者：Vamei 出处：http://www.cnblogs.com/vamei 欢迎转载，也请保留这段声明。谢谢！

随机变量的函数

在前面的文章中，我先将概率值分配给各个事件，得到事件的概率分布。

通过事件与随机变量的映射，让事件“数值化”，事件的概率值转移到随机变量上，获得随机变量的概率分布。

我们使用随机变量的函数，来定制新的随机变量。随机变量的函数是从旧有的随机变量到一个新随机变量的映射。通过函数的映射功能，原有随机变量对应新的随机变量。通过原有随机变量的概率分布，我们可以获知新随机变量的概率分布。事件，随机变量，随机变量函数的关系如下:

一个简单的例子是掷硬币。出现正面的话，我赢1个筹码，负面的话，我输1个筹码。那么，投掷一次，赢的筹码数是一个随机变量X，X可能取值为1和-1。因此X的分布为:

$P(1) = 0.5$

$P(-1) = 0.5$

换一个角度来思考，我们将正负面“换算”成输赢的钱。如果一个筹码需要10元钱买，那么投掷一次硬币，赢的钱是一个随机变量Y，且 $Y = 10X$ 。Y的分布为:

$P(10) = 0.5$

$P(-10) = 0.5$

Y实际上是随机变量X的一个函数。X的1对应Y的10，X的-1对应Y的-10。即 $Y = 10X$

小总结，在上面的实验中，硬币为正面为一个事件。赢得的筹码数为一个随机变量X。赢得的钱是X的函数Y，它也是一个随机变量。

随机变量的函数还可以是多变量函数， $Y = g(X_1, X_2, ..., X_n)$ 。Y的值y对应的是多维空间的点 $(x_1, x_2,..., x_n)$ 。比如掷硬币，第一次赢的筹码为 $X_1$ ，第二次赢的筹码为 $X_2$ 。我们可以构成一个新的随机变量 $Y = X_1 + X_2$ ，即两次赢得的筹码的总和。

获得新概率分布的基本方法

一个核心问题是，如何通过X的概率分布，来获得 $Y=g(X)$ 的概率分布。基本的思路是，如果我们想知道Y取某个值y的概率，可以找到对应的X值x的概率。这两个概率相等。

因此，我们使用如下方法来获得Y的概率。如果有函数关系 $Y=g(X_1, X_2, ..., X_n)$ ，获得Y分布的基本方法是:

1. 通过 $Y=g(X_1, X_2, ..., X_n)$ ，找到对应 $\{ Y \le y \}$ 的 $(x_1, x_2, ..., x_n)$ 区间I。

2. 在区间I上，积分 $f(x_1, x_2, ..., x_n)$ ，获得 $P(Y \le y)$

3. 通过微分，获得密度函数。

如果有函数关系 $Y = X^2$ ，而X满足下面的分布:

$f(x) = \frac{1}{\sqrt{2\pi}}e^{-x^2/2}$

对于任意 $y \ge 0$ 来说，

$F(y) = P(Y \le y) = P(X^2 \le y) = P(-\sqrt{y} \le X \le \sqrt{y})$

$F(y) = \int_{-\sqrt{y}}^{\sqrt{y}}\frac{1}{\sqrt{2\pi}}e^{-x^2/2}dx = 2 \int_{0}^{\sqrt{y}}\frac{1}{\sqrt{2\pi}}e^{-x^2/2}dx$

对上面的F(y)微分，即获得密度函数

$f(y) = \frac{1}{\sqrt{2\pi}}y^{-1/2}e^{-y/2}, 0 \le y \le \infty$

绘制密度函数

import numpy as np
import matplotlib.pyplot as plt

pi = np.pi

x = np.linspace(-10, 10, 200)
y = np.linspace(0.1, 10, 100)

fx = 1/np.sqrt(2*pi)*np.exp(-x**2/2)
fy = 1/np.sqrt(2*pi)*(y**(-1/2))*np.exp(-y/2)

plt.plot(x, fx, color = "red", label="X distribution")
plt.plot(y, fy, label="Y distribution")

plt.title("Y = X*X")
plt.xlabel("RV")
plt.ylabel("pdf")

plt.legend()

plt.show()

上面的例子展示的是单变量函数，我们看一个多变量函数的例子。即 $Y=g(X_1, X_2, ..., X_n)$ ，且已知 $X_1, X_2, ..., X_n$ 的联合分布为 $f(x_1, x_2, ..., x_n)$ 。我们需要找到满足 $g(x_1, x_2, ..., x_n) \le y$ 的区间。

比如， $Y = X_1 + X_2$ ，且 $X_1, X_2$ 满足如下分布:

$f(x_1, x_2) = \frac{1}{2 \pi} \exp \left( -\frac{1}{2} \left( x_1^2 + x_2^2 \right) \right)$

为了让 $x_1 + x_2 \le y$ ，我们可以让 $x_1$ 任意取值，而让 $x_2 \le y - x_1$

$F_Y(y) = \int_{-\infty}^{\infty} \int_{-\infty}^{y - x_1} f(x_1, x_2) dx_2dx_1$

让x_2 = v - x_1，有

$F_Y(y) = \int_{-\infty}^{\infty} \int_{-\infty}^{y} f(x_1, v - x_1)dvdx_1 = \int_{-\infty}^{y} \int_{-\infty}^{\infty}f(x_1, v - x_1)dvdx_1$

微分，可得y的分布为:

$f_Y(y) = \int_{-\infty}^{\infty} f(x_1, y - x_1) dx_1 = \int_{-\infty}^{\infty} \frac{1}{2 \pi} \exp \left( -\frac{1}{2} \left( x_1^2 + (y - x_1)^2 \right) \right) dx_1$

上述方程也可以使用数值方法求解:

代码如下:

# By Vamei

import numpy as np
import scipy.integrate
import matplotlib.pyplot as plt
pi = np.pi

'''
core of the integral
'''
def int_core(y):
    f = lambda x: 1.0/(2*pi)*np.exp(-0.5*(x**2 + (y-x)**2))
    return f

'''
calculate f(y)
'''
def density(y):
    rlt = scipy.integrate.quad(int_core(y), -np.inf, np.inf)
    return rlt[0]

# get distribution
y  = np.linspace(-10, 10, 100)
fy = map(density, y)

plt.plot(y, fy)
plt.title("PDF of X1+X2")
plt.ylabel("f(y)")
plt.xlabel("y")
plt.show()

上面的int_core()函数是一个闭包，它表示积分核部分。density()函数用于求某个y值下的积分结果。

(我们也可以利用解析的方法，推导出f(y)满足分布 $N(0, \sqrt{2})$ 。如果有微积分基础，可以将此作为练习。)

单变量函数的通用公式

上面求新的随机变量分布的步骤较为繁琐。在一些特殊情况下，我们可以直接代入通用公式，来获得新的分布。

(通用公式实际上是从基本方法推导出的数学表达式)

对于单变量函数来说，如果 $Y=g(X)$ ，g是一个可微并且单调变化的函数 (在该条件，存在反函数 $g^{-1}$ ，使得 $X=g^{-1}(Y))$ 。那么我们可以使用下面的通用公式，来获得Y的分布:

$f_Y(y) = f_X(g^{-1}(y)) \cdot \frac{d}{dy}g^{-1}(y)$

假设X为标准分布，即 $N(0, 1)$ ，且 $Y = 5X + 1$ ，那么 $g^{-1}(y) = (y - 1)/5$ ，因此:

$f_Y(y) = f_X((y-1)/5) \cdot (1/5) = \frac{1}{5\sqrt{2\pi}}e^{-(y-1)^2/(2 \times 25)}$

可以看到，新的分布是一个 $\mu = 1, \sigma=5$ 的正态分布，即 $N(1, 5)$

并不是所有的函数都有反变换，所以这里的“通用”公式并不能适用于所有的情况。

多变量函数的通用公式

在一些特殊情况下，我们可以使用多变量函数的通用公式。

如果 $U=g_1(X, Y), V=g_2(X, Y)$ ，且存在反变换，使得

$X = h_1(U, V)$

$Y = h_2(U, V)$

那么，我们可以通过如下公式，从X,Y的分布获得U,V的联合分布:

$f_{UV}(u, v) = f_{XY}(h_1(u, v), h_2(u, v))|J|$

J表示雅可比变换(Jacobian tranformation)，表示如下

$J = \left| \begin{array}{cc} \frac{\partial x}{\partial u} & \frac{\partial x}{\partial v} \\ \frac{\partial y}{\partial u} & \frac{\partial y}{\partial v} \end{array} \right| =\frac{\partial x}{\partial u}\frac{\partial y}{\partial v}-\frac{\partial x}{\partial v}\frac{\partial y}{\partial u}$

如果X和Y是独立的随机变量，且有相同的分布