连续概率 -- 正态分布

2017-05-21 Sunday    


Normal distribution 正态分布 也称为 Gaussian Distribution 高斯分布,这是一个在数学、物理、工程等领域都非常重要的概率分布,在统计学的很多方面都有着重大影响力。

简介

假设随机变量 $x$ 服从一个均值为 $\mu$ ,标准差为 $\sigma$ 的正态分布,则可以记为:

\[X\sim N(\mu ,\sigma ^{2})\]

对应的概率密度函数为。

\[f(x)={1 \over \sigma {\sqrt {2\pi }}}\,e^{-{(x-\mu )^{2} \over 2\sigma ^{2}}}\]

normal distribution example

该图是通过如下的代码生成。

import numpy as np
import scipy.stats as stats
import matplotlib.pyplot as plt

colors = ['deepskyblue', 'red', 'green', 'blue', 'purple']

x = np.linspace(-10, 20, num=200)
fig = plt.figure(figsize=(12, 6))
for mu, sigma, c in zip([0.5, 1, 2, 3, 4], [3, 1, 3, 2, 5], colors):
	y = stats.norm.pdf(x, mu, sigma)
	plt.plot(x, y, lw=2, c=c, label=r"$\mu={0:.1f}, \sigma={1:.1f}$".format(mu, sigma))
	plt.fill_between(x, y, color=c, alpha=.1)

plt.legend(loc=0)
plt.ylabel("PDF($x$)")
plt.xlabel("$x$")
plt.show()

注意,matplotlib.mlab 中的 normpdf() 已经不再支持了,可以使用 scipy.stats.norm.pdf() 函数替换。另外,也可以通过如下方式计算。

def norm_pdf(x, mu, sigma):
    return np.exp(-((x - mu)**2)/(2*(sigma**2)))/(sigma*np.sqrt(2*np.pi))

随机采样

通过 numpy 中的 random.normal() 函数可以对正态分布进行采样。如下示例,生成均值为 0 方差为 1 的正态分布,可以进行测试,可以查看其均值和方差。

import numpy as np
import matplotlib.pyplot as plt

mean = 0
sigma = 1
arr = np.random.normal(mean, sigma, size=1000)
abs(mean - np.mean(arr)) < 0.01
abs(sigma - np.std(arr, ddof = 1)) < 0.01

plt.plot(arr)
plt.show()

在如上程序计算样本方差时,因为 ddof=1 也就意味者计算的是样本方差的无偏估计。

另外,也可以通过如下的方式绘制一个正弦 sin 曲线,同时增加一些高斯噪声信息。

import numpy as np
import matplotlib.pyplot as plt

mean = 0
sigma = 0.5
xarr = np.arange(0.0, 10.0, 0.01)
garr = np.random.normal(mean, sigma, size=1000)

plt.plot(xarr, 10 * np.sin(2 * np.pi * xarr) + garr)
plt.show()

采样示例

normal distribution sample

import numpy as np
import scipy.stats as stats
import matplotlib.pyplot as plt

mu, sigma = 0, 1
size = 10000
np.random.seed(0)

s = np.random.normal(mu, sigma, size)
n, bins, patches = plt.hist(s, bins=100, density=True)
y = stats.norm.pdf(bins, mu, sigma)

plt.plot(bins, y, 'r--')
plt.xlabel('Expectation')
plt.ylabel('Probability')
plt.title('histogram of normal distribution: $\mu = 0$, $\sigma=1$')
plt.show()

多元高斯分布

多元高斯分布 (Multivariate Gaussian Distribution) 的形式很简单,就是一元高斯分布的在向量形式中的推广,把向量 $X=[x_1, x_2,\cdots, x_n]^T$ 称作是均值为 $\mu \in {\bf R}^n$ ,协方差矩阵为 $\Sigma \in S^n$ 的多元高斯分布,其概率密度函数的形式为。

\[f_x(x_1, \cdots ,x_n)=\frac{1}{\sqrt{(2 \pi)^n |\Sigma|}} e^{-\frac{(x-\mu)^T (x-\mu)}{2\Sigma}}\]

这种分布由均值和协方差矩阵确定,类似于一维正态分布中的均值和方差,可以看到,与一维的高斯分布有两点不同。

  • 多维高斯公式把一维中的 $(x - \mu)^2$ 替换成了 $(x - \mu)^T (x - \mu)$ 矩阵操作,而该公式也可以写成 $(x_1 - \mu_1)^2 + (x_1 - \mu_1 )^2 + \cdots + (x_n - \mu_n )^2$,那么一维就只对应了 $(x_1 - \mu_1 )^2$ 。
  • 多维高斯分布把方差 $\sigma^2$ 换为了协方差矩阵 $\Sigma$,也就是说到了高维之后,不仅仅有各维自己的方差,还有两两之间的协方差。

协方差反应两个变量的相关程度,将方差和协方差放到一个矩阵中就是协方差矩阵 $\Sigma$ 。

其它

在 Python 中,可以通过 numpy.random.multivariate_normal() 方法用于根据实际情况生成一个多元正态分布矩阵,该函数定义如下:

multivariate_normal(mean, cov, size=None, check_valid=None, tol=None)
    mean        1*N 均值
    cov         N*N 协方差矩阵,对称且为半正定的矩阵。
    size        生成数据的维度,例如size=(10, 100则输出的矩阵为10*100*N
    check_valid 检查协方差矩阵是否合法。

如下是个一个参考的示例。

import numpy as np
import matplotlib.pyplot as plt

## Gaussian Distribution
#cov = np.eye(1) * 3
#mean = np.array([3])
#Y = np.random.multivariate_normal(mean, cov, 10, 'raise')

## Multivariate Gaussian Distribution
cov = np.eye(2) * 3
mean = np.array([3, 100])
Y = np.random.multivariate_normal(mean, cov, (10, 3), 'raise')

如上生成了一个 10 * 3 * 2 的数据,因为协方差矩阵使用的是对角矩阵,也就意味着各个变量之前是独立的,其中均值对应 mean 变量,方差就都是 1 ,可以通过如下方式计算均值、方差等信息。

>>> np.mean(Y, axis=0)  # 计算均值,同样可以计算方差等信息
array([[  2.96771613,  99.65476306],
       [  2.8397508 ,  99.88631836],
       [  3.09584658, 100.50960993]])
>>> np.std(Y, axis=0)
array([[0.98685019, 0.99880091],
       [0.98338639, 1.02355389],
       [1.00234377, 0.98393021]])

如下是一个示例。

二元绘图

normal distribution 2d example

实现代码如下。

import numpy as np
import scipy.stats as stats
import matplotlib.pyplot as plt
from mpl_toolkits.axes_grid1 import make_axes_locatable

size = 10000
colors = ['deepskyblue', 'red', 'green', 'blue', 'purple']

male_mu = np.array([172, 60])
male_sigma = np.array([[36, 15], [15, 9]])
male = np.random.multivariate_normal(male_mu, male_sigma, size)

female_mu = np.array([165, 56])
female_sigma = np.array([[25, 12], [12, 8]])
female = np.random.multivariate_normal(female_mu, female_sigma, size)

fig = plt.figure(figsize=(8, 8))
axes = fig.gca()  # get current axes
axes.scatter(*female.T, c='red', label='female', s=10, alpha=0.1);
axes.scatter(*male.T, c='blue', label='male', s=10, alpha=0.1);

def figure_rug(axes, x, y, color='blue'):
	xlims = axes.get_xlim()
	ylims = axes.get_ylim()

	zx1 = np.zeros(len(x))
	zx2 = np.zeros(len(x))
	zx1.fill(ylims[1])
	zx2.fill(xlims[1])

	axes.plot(x, zx1, marker='|', color=color, markersize=8, alpha=0.2)
	axes.plot(zx2, y, marker='_', color=color, markersize=8, alpha=0.2)

	axes.set_xlim(xlims)
	axes.set_ylim(ylims)

figure_rug(axes, female.T[0], female.T[1], 'red')
figure_rug(axes, male.T[0], male.T[1], 'blue')

axes.grid(False)
axes.legend(loc='upper left');
plt.show()

数据与上面相同,只是绘制的方式有所区别。

normal distribution 2d example

import numpy as np
import scipy.stats as stats
import matplotlib.pyplot as plt
from mpl_toolkits.axes_grid1 import make_axes_locatable

size = 10000
colors = ['deepskyblue', 'red', 'green', 'blue', 'purple']

male_mu = np.array([172, 60])
male_sigma = np.array([[36, 15], [15, 9]])
male = np.random.multivariate_normal(male_mu, male_sigma, size)

female_mu = np.array([165, 56])
female_sigma = np.array([[25, 12], [12, 8]])
female = np.random.multivariate_normal(female_mu, female_sigma, size)

fig = plt.figure(figsize=(8, 8))
axes = fig.gca()  # get current axes
axes.scatter(*female.T, c='red', label='female', s=10, alpha=0.1);
axes.scatter(*male.T, c='blue', label='male', s=10, alpha=0.1);

def remove_border(axes=None, top=False, right=False, left=True, bottom=True):
	ax = axes or plt.gca()
	ax.spines['top'].set_visible(top)
	ax.spines['right'].set_visible(right)
	ax.spines['left'].set_visible(left)
	ax.spines['bottom'].set_visible(bottom)
	
	#turn off all ticks
	ax.yaxis.set_ticks_position('none')
	ax.xaxis.set_ticks_position('none')
	
	#now re-enable visibles
	if top:
		ax.xaxis.tick_top()
	if bottom:
		ax.xaxis.tick_bottom()
	if left:
		ax.yaxis.tick_left()
	if right:
		ax.yaxis.tick_right()
def make_mhist(axes, x, y, color='blue'):
	ax2 = axes[1]
	ax3 = axes[2]
	for tl in (ax2.get_xticklabels() + ax2.get_yticklabels() +
			ax3.get_xticklabels() + ax3.get_yticklabels()):
		tl.set_visible(False)
	style = {'histtype':'step', 'color':color, 'alpha':0.4}
	ax2.hist(x, np.linspace(np.min(x), np.max(x), 200), **style)
	ax3.hist(y, np.linspace(np.min(y), np.max(y), 200), orientation='horizontal', **style)

divider = make_axes_locatable(axes)
axesx = divider.append_axes("top", 1.5, pad=0.0, sharex=axes)
axesy = divider.append_axes("right", 1.5, pad=0.0, sharey=axes)
axesx.grid(False)
axesy.grid(False)

remove_border(axesx, right=True, left=False)
remove_border(axesy, right=False, left=True, bottom=False, top=True)

make_mhist([axes, axesx, axesy], male.T[0], male.T[1], color='blue')
make_mhist([axes, axesx, axesy], female.T[0], female.T[1], color='red')

fig.subplots_adjust(left=0.15, right=0.95)
axes.grid(False)
axes.legend(loc='upper left');
plt.show()


如果喜欢这里的文章,而且又不差钱的话,欢迎打赏个早餐 ^_^