Task3：概率论与数理统计（基于Python）

随机现象与概率

当n足够大时，频率可以近似看成概率

#模拟频率近似概率
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
%matplotlib inline
plt.style.use("ggplot")
import warnings 
warnings.filterwarnings("ignore")
plt.rcParams['font.sans-serif']=['SimHei','Songti SC','STFangsong']
plt.rcParams['axes.unicode_minus'] = False  # 用来正常显示负号
import seaborn as sns 

# 模拟抛硬币正面的概率是否会越来越接近0.5
import random
def Simulate_coin(test_num):
    random.seed(100)
    coin_list = [1 if (random.random()>=0.5) else 0  for i in range(test_num)]   # 模拟试验结果
    coin_frequence = np.cumsum(coin_list) / (np.arange(len(coin_list))+1)  # 计算正面为1的频率
    plt.figure(figsize=(10,6))
    plt.plot(np.arange(len(coin_list))+1, coin_frequence, c='blue', alpha=0.7)
    plt.axhline(0.5,linestyle='--',c='red',alpha=0.5)
    plt.xlabel("test_index")
    plt.ylabel("frequence")
    plt.title(str(test_num)+" times")
    plt.show()

Simulate_coin(test_num = 500)
Simulate_coin(test_num = 1000)
Simulate_coin(test_num = 5000)
Simulate_coin(test_num = 10000)

条件概率、乘法公式、全概率公式与贝叶斯公式

事件B下事件A的条件概率

乘法公式

全概率公式

贝叶斯公式

#三羊问题
import random 
class MontyHall:
    def __init__(self,n):
        """
        n : int,试验的次数
        """
        self.n = n
        self.change = 0  # 记录换才能拿到车的次数
        self.No_change = 0  # 记录不换才能拿到车的次数
    def start(self):
        for i in range(self.n):
            door_list = [1,2,3]   ## 三扇门
            challenger_door = random.choice(door_list)   ## 随机选择了其中一扇
            car_door = random.choice(door_list)     ## 随机选定车的门
            door_list.remove(challenger_door)   ## 没有被挑战者选中的剩下的门
            if challenger_door == car_door:
                host_door = random.choice(door_list)
                door_list.remove(host_door)        # 不换才能拿车
                self.No_change += 1
            else:
                self.change += 1      # 换了才能拿车
        print("换且能拿到车的概率：%.2f " % (self.change/self.n * 100) + "%")
        print("不换也能拿到车的概率：%.2f"% (self.No_change/self.n * 100) + "%")

if __name__ == "__main__":
    mh = MontyHall(1000000)
    mh.start()

注：若门有100扇，换且能拿到车的概率为99/100，不拿也能换到车的概率为1/100

已知密度函数求分布函数

#已知柯西分布的密度函数求分布函数
from sympy import *
x=symbols('x')
p_x=1/pi*(1/(1+x**2))
integrate(p_x, (x, -oo, x))

已知分布函数求密度函数

#已知柯西分布的分布函数求密度函数
from sympy import *
x = symbols('x')
f_x = 1/pi*(atan(x)+pi/2)
diff(f_x,x,1)

常见的连续型随机变量及其密度函数

均匀分布

密度函数：

分布函数：

#【0，1】上的均匀分布
import numpy as np
a=float(0)
b=float(1)
x=np.linspace(a,b)
y=np.full(shape=len(x),fill_value=1/(b-a))    # np.full 构造一个数组，用指定值填充其元素

plt.plot(x,y,"b",linewidth=2)
plt.ylim(0,1.2)
plt.xlim(-1,2)
plt.xlabel('X')
plt.ylabel('p(x)')
plt.title('uniform distribution')
plt.show()

指数分布

密度函数：

分布函数：

# 指数分布
lam = float(1.5)
x = np.linspace(0,20,100)   #把0~15分为100份
y = lam * np.e**(-lam * x)

plt.plot(x,y,"b",linewidth=2) 
plt.xlim(-5,10)
plt.xlabel('X')
plt.ylabel('p (x)')
plt.title('exponential distribution')
plt.show()

正态分布

密度函数：

#正态分布
import math
mu = float(0)
mu1 = float(2)
sigma1 = float(1)
sigma2 = float(1.25)*float(1.25)
sigma3 = float(0.25)
x = np.linspace(-5, 5, 1000)
y1 = np.exp(-(x - mu)**2 / (2 * sigma1**2)) / (math.sqrt(2 * math.pi) * sigma1)
y2 = np.exp(-(x - mu)**2 / (2 * sigma2**2)) / (math.sqrt(2 * math.pi) * sigma2)
y3 = np.exp(-(x - mu)**2 / (2 * sigma3**2)) / (math.sqrt(2 * math.pi) * sigma3)
y4 = np.exp(-(x - mu1)**2 / (2 * sigma1**2)) / (math.sqrt(2 * math.pi) * sigma1)
plt.plot(x,y1,"b",linewidth=2,label=r'$\mu=0,\sigma=1$') 
plt.plot(x,y2,"orange",linewidth=2,label=r'$\mu=0,\sigma=1.25$') 
plt.plot(x,y3,"yellow",linewidth=2,label=r'$\mu=0,\sigma=0.5$') 
plt.plot(x,y4,"b",linewidth=2,label=r'$\mu=2,\sigma=1$',ls='--') 
plt.axvline(x=mu,ls='--')
plt.text(x=0.05,y=0.5,s=r'$\mu=0$')
plt.axvline(x=mu1,ls='--')
plt.text(x=2.05,y=0.5,s=r'$\mu=2$')
plt.xlim(-5,5)
plt.xlabel('X')
plt.ylabel('p (x)')
plt.title('normal distribution')
plt.legend()
plt.show()

## 已知正态分布的密度函数求分布函数
from sympy import *
from sympy.abc import mu,sigma
x = symbols('x')
p_x = 1/(sqrt(2*pi)*sigma)*E**(-(x-mu)**2/(2*sigma**2))
integrate(p_x, (x, -oo, x))

使用scipy计算密度函数画图(非自定义函数)

from scipy.stats import expon # 指数分布
x = np.linspace(0.01,10,1000)  
plt.plot(x, expon.pdf(x),'r-', lw=5, alpha=0.6, label='expon pdf')    # pdf表示求密度函数值
plt.xlabel("X")
plt.ylabel("p (x)")
plt.legend()
plt.show()

使用scipy计算分布函数画图(非自定义函数)

from scipy.stats import expon
x = np.linspace(0.01,10,1000)  
plt.plot(x, expon.cdf(x),'r-', lw=5, alpha=0.6, label='expon cdf')  # cdf表示求分布函数值
plt.xlabel("X")
plt.ylabel("F (x)")
plt.legend()
plt.show()

0-1分布
二项分布
泊松分布

# 对比不同的lambda对泊松分布的影响
import math
# 构造泊松分布列的计算函数
def poisson(lmd,x):
    return pow(lmd,x)/math.factorial(x)*math.exp(-lmd)
x = [i+1 for i in range(10)]
lmd1 = 0.8
lmd2 = 2.0
lmd3 = 4.0
lmd4 = 6.0
p_lmd1 = [poisson(lmd1,i) for i in x]
p_lmd2 = [poisson(lmd2,i) for i in x]
p_lmd3 = [poisson(lmd3,i) for i in x]
p_lmd4 = [poisson(lmd4,i) for i in x]

plt.scatter(np.array(x), p_lmd1, c='b',alpha=0.7)
plt.axvline(x=lmd1,ls='--')
plt.text(x=lmd1+0.1,y=0.1,s=r"$\lambda=0.8$")
plt.ylim(-0.1,1)
plt.xlabel("X")
plt.ylabel("p (x)")
plt.title(r"$\lambda = 0.8$")
plt.show()

plt.scatter(np.array(x), p_lmd2, c='b',alpha=0.7)
plt.axvline(x=lmd2,ls='--')
plt.text(x=lmd2+0.1,y=0.1,s=r"$\lambda=2.0$")
plt.ylim(-0.1,1)
plt.xlabel("X")
plt.ylabel("p (x)")
plt.title(r"$\lambda = 2.0$")
plt.show()

plt.scatter(np.array(x), p_lmd3, c='b',alpha=0.7)
plt.axvline(x=lmd3,ls='--')
plt.text(x=lmd3+0.1,y=0.1,s=r"$\lambda=4.0$")
plt.ylim(-0.1,1)
plt.xlabel("X")
plt.ylabel("p (x)")
plt.title(r"$\lambda = 4.0$")
plt.show()

plt.scatter(np.array(x), p_lmd4, c='b',alpha=0.7)
plt.axvline(x=lmd4,ls='--')
plt.text(x=lmd4+0.1,y=0.1,s=r"$\lambda=6.0$")
plt.ylim(-0.1,1)
plt.xlabel("X")
plt.ylabel("p (x)")
plt.title(r"$\lambda = 6.0$")
plt.show()

在lmd附近概率较大
随着lmd的增加，分布逐渐趋于正态分布

使用scipy的pmf和cdf画图

from scipy.stats import binom
n=10
p = 0.5
x=np.arange(1,n+1,1)
pList=binom.pmf(x,n,p)
plt.plot(x,pList,marker='o',alpha=0.7,linestyle='None')
'''
vlines用于绘制竖直线(vertical lines),
参数说明：vline(x坐标值, y坐标最小值, y坐标值最大值)
'''
plt.vlines(x, 0, pList)
plt.xlabel('随机变量：抛硬币10次')
plt.ylabel('概率')
plt.title('二项分布：n=%d,p=%0.2f' % (n,p))
plt.show()

一维随机变量的数字特征：期望、方差、分位数与中位数

均值

常见的分布的数学期望：
方差与标准差

方差的性质：

常见分布的方差：

# 使用scipy计算常见分布的均值与方差：(如果忘记公式的话直接查，不需要查书了)
from scipy.stats import bernoulli   # 0-1分布
from scipy.stats import binom   # 二项分布
from scipy.stats import poisson  # 泊松分布
from scipy.stats import rv_discrete # 自定义离散随机变量
from scipy.stats import uniform # 均匀分布
from scipy.stats import expon # 指数分布
from scipy.stats import norm # 正态分布
from scipy.stats import rv_continuous  # 自定义连续随机变量

print("0-1分布的数字特征：均值:{}；方差:{}；标准差:{}".format(bernoulli(p=0.5).mean(), 
                                  bernoulli(p=0.5).var(), 
                                  bernoulli(p=0.5).std()))
print("二项分布b(100,0.5)的数字特征：均值:{}；方差:{}；标准差:{}".format(binom(n=100,p=0.5).mean(), 
                                  binom(n=100,p=0.5).var(), 
                                  binom(n=100,p=0.5).std()))
## 模拟抛骰子的特定分布
xk = np.arange(6)+1
pk = np.array([1.0/6]*6)
print("泊松分布P(0.6)的数字特征：均值:{}；方差:{}；标准差:{}".format(poisson(0.6).mean(), 
                                  poisson(0.6).var(), 
                                  poisson(0.6).std()))
print("特定离散随机变量的数字特征：均值:{}；方差:{}；标准差:{}".format(rv_discrete(name='dice', values=(xk, pk)).mean(), 
                                  rv_discrete(name='dice', values=(xk, pk)).var(), 
                                  rv_discrete(name='dice', values=(xk, pk)).std()))
print("均匀分布U(1,1+5)的数字特征：均值:{}；方差:{}；标准差:{}".format(uniform(loc=1,scale=5).mean(), 
                                  uniform(loc=1,scale=5).var(), 
                                  uniform(loc=1,scale=5).std()))
print("正态分布N(0,0.0001)的数字特征：均值:{}；方差:{}；标准差:{}".format(norm(loc=0,scale=0.01).mean(), 
                                  norm(loc=0,scale=0.01).var(), 
                                  norm(loc=0,scale=0.01).std()))

lmd = 5.0  # 指数分布的lambda = 5.0
print("指数分布Exp(5)的数字特征：均值:{}；方差:{}；标准差:{}".format(expon(scale=1.0/lmd).mean(), 
                                  expon(scale=1.0/lmd).var(), 
                                  expon(scale=1.0/lmd).std()))

## 自定义标准正态分布
class gaussian_gen(rv_continuous):
    def _pdf(self, x): # tongguo 
        return np.exp(-x**2 / 2.) / np.sqrt(2.0 * np.pi)
gaussian = gaussian_gen(name='gaussian')
print("标准正态分布的数字特征：均值:{}；方差:{}；标准差:{}".format(gaussian().mean(), 
                                  gaussian().var(), 
                                  gaussian().std()))

## 自定义指数分布
import math
class Exp_gen(rv_continuous):
    def _pdf(self, x,lmd):
        y=0
        if x>0:
            y = lmd * math.e**(-lmd*x)
        return y
Exp = Exp_gen(name='Exp(5.0)')
print("Exp(5.0)分布的数字特征：均值:{}；方差:{}；标准差:{}".format(Exp(5.0).mean(), 
                                  Exp(5.0).var(), 
                                  Exp(5.0).std()))

## 通过分布函数自定义分布
class Distance_circle(rv_continuous):                 #自定义分布xdist
    """
    向半径为r的圆内投掷一点，点到圆心距离的随机变量X的分布函数为:
    if x<0: F(x) = 0;
    if 0<=x<=r: F(x) = x^2 / r^2
    if x>r: F(x)=1
    """
    def _cdf(self, x, r):                   #累积分布函数定义随机变量
        f=np.zeros(x.size)                  #函数值初始化为0
        index=np.where((x>=0)&(x<=r))           #0<=x<=r
        f[index]=((x[index])/r[index])**2       #0<=x<=r
        index=np.where(x>r)                     #x>r
        f[index]=1                              #x>r
        return f
dist = Distance_circle(name="distance_circle")
print("dist分布的数字特征：均值:{}；方差:{}；标准差:{}".format(dist(5.0).mean(), 
                                  dist(5.0).var(), 
                                  dist(5.0).std()))

分位数与中位数

from scipy.stats import norm
print("标准正态分布的0.25分位数：",norm(loc=0,scale=1).ppf(0.25))   # 使用ppf计算分位数点
print("标准正态分布的0.5分位数：",norm(loc=0,scale=1).ppf(0.5))
print("标准正态分布的0.75分位数：",norm(loc=0,scale=1).ppf(0.75))
print("标准正态分布的0.95分位数：",norm(loc=0,scale=1).ppf(0.95))
#标准正态分布的0.25分位数： -0.6744897501960817
#标准正态分布的0.5分位数： 0.0
#标准正态分布的0.75分位数： 0.6744897501960817
#标准正态分布的0.95分位数： 1.6448536269514722

多维随机变量及其联合分布、边际分布、条件分布

联合分布：

以二维连续为例：

以二维离散为例：

或者
绘制二维正态分布的联合概率密度曲面图

from scipy.stats import multivariate_normal
from mpl_toolkits.mplot3d import axes3d
x, y = np.mgrid[-5:5:.01, -5:5:.01]  # 返回多维结构
pos = np.dstack((x, y))
rv = multivariate_normal([0.5, -0.2], [[2.0, 0.3], [0.3, 0.5]])
z = rv.pdf(pos)
plt.figure('Surface', facecolor='lightgray',figsize=(12,8))
ax = plt.axes(projection='3d')
ax.set_xlabel('X', fontsize=14)
ax.set_ylabel('Y', fontsize=14)
ax.set_zlabel('P (X,Y)', fontsize=14)
ax.plot_surface(x, y, z, rstride=50, cstride=50, cmap='jet')
plt.show()

绘制二维正态分布的联合概率密度等高线图

from scipy.stats import multivariate_normal
x, y = np.mgrid[-1:1:.01, -1:1:.01]
pos = np.dstack((x, y))
rv = multivariate_normal([0.5, -0.2], [[2.0, 0.3], [0.3, 0.5]])
z = rv.pdf(pos)
fig = plt.figure(figsize=(8,6))
ax2 = fig.add_subplot(111)
ax2.set_xlabel('X', fontsize=14)
ax2.set_ylabel('Y', fontsize=14)
ax2.contourf(x, y, z, rstride=50, cstride=50, cmap='jet')
plt.show()

边际分布

X的边际分布

Y的边际分布
边际密度函数
求边际密度函数 p_{X}(x)

from sympy import *
x = symbols('x')
y = symbols('y')
p_xy = Piecewise((1,And(x>0,x<1,y<x,y>-x)),(0,True))
integrate(p_xy, (y, -oo, oo))   ## 由于0<x<1时候，那么x>-x，即2x

求边际密度函数 p_{Y}(y)

from sympy import *
x = symbols('x')
y = symbols('y')
p_xy = Piecewise((1,And(x>0,x<1,y<x,y>-x)),(0,True))
integrate(p_xy, (x, -oo, oo))   ## 由于|y|<x,0<x<1时，因此y肯定在(-1,1)

边际分布列

求边际分布函数 p_{Y}(y)
from sympy import *
from sympy.abc import lamda,m,p,k
x = symbols('x')
y = symbols('y')
f_p = lamda**m/factorial(m)*E**(-lamda)*factorial(m)/(factorial(k)*factorial(m-k))*p**k*(1-p)**(m-k)
summation(f_p, (m, k, +oo))

协方差

性质：

例子1：

# 求协方差
from sympy import *
from sympy.abc import lamda,m,p,k
x = symbols('x')
y = symbols('y')
p_xy = Piecewise((3*x,And(y>0,y<x,x<1)),(0,True))
E_xy = integrate(x*y*p_xy, (x, -oo, oo),(y,-oo,oo))
E_x = integrate(x*p_xy, (x, -oo, oo),(y,-oo,oo))
E_y = integrate(y*p_xy, (x, -oo, oo),(y,-oo,oo))
E_xy - E_x*E_y

例子二：

求

# 方法一：先计算边际密度函数，在计算特征数
from sympy import *
x = symbols('x')
y = symbols('y')
p_xy = Piecewise((1/3*(x+y),And(x>0,x<1,y>0,y<2)),(0,True))
p_x = integrate(p_xy, (y, -oo, oo))  # x边际密度函数
p_y = integrate(p_xy, (x, -oo, oo))  # y边际密度函数
E_x2 = integrate(x**2*p_x, (x, -oo, oo))
E_x = integrate(x*p_x, (x, -oo, oo))
E_y2 = integrate(y**2*p_y, (y,-oo,oo))
E_y = integrate(y*p_y, (y,-oo,oo))
E_xy = integrate(x*y*p_xy, (x, -oo, oo),(y,-oo,oo))
cov_xy = E_xy - E_x*E_y
var_x = E_x2 - E_x**2
var_y = E_y2 - E_y**2
var_2x_3y_8 = 4*var_x + 9*var_y -12*cov_xy
var_2x_3y_8

# 方法二：直接通过联合密度函数计算特征数
from sympy import *
x = symbols('x')
y = symbols('y')
p_xy = Piecewise((1/3*(x+y),And(x>0,x<1,y>0,y<2)),(0,True))
E_x2 = integrate(x**2*p_xy, (x, -oo, oo),(y, -oo, oo))
E_x = integrate(x*p_xy, (x, -oo, oo),(y, -oo, oo))
E_y2 = integrate(y**2*p_xy, (x, -oo, oo),(y,-oo,oo))
E_y = integrate(y*p_xy, (x, -oo, oo),(y,-oo,oo))
E_xy = integrate(x*y*p_xy, (x, -oo, oo),(y,-oo,oo))
cov_xy = E_xy - E_x*E_y
var_x = E_x2 - E_x**2
var_y = E_y2 - E_y**2
var_2x_3y_8 = 4*var_x + 9*var_y -12*cov_xy
var_2x_3y_8

# 方法三：令z=2*x-3*y+8,使用VAR(Z) = E（Z**2）- E（Z）**2
from sympy import *
x = symbols('x')
y = symbols('y')
p_xy = Piecewise((1/3*(x+y),And(x>0,x<1,y>0,y<2)),(0,True))
E_z2 = integrate((2*x-3*y+8)**2*p_xy, (x, -oo, oo),(y, -oo, oo))
E_z = integrate((2*x-3*y+8)*p_xy, (x, -oo, oo),(y, -oo, oo))
E_z2 - E_z**2

协方差矩阵（对称）

# 求协方差矩阵：1.求两两变量的协方差和各自的方差；2. 组合成矩阵
from sympy import *
x = symbols('x')
y = symbols('y')
p_xy = Piecewise((1/3*(x+y),And(x>0,x<1,y>0,y<2)),(0,True))
p_x = integrate(p_xy, (y, -oo, oo))  # x边际密度函数
p_y = integrate(p_xy, (x, -oo, oo))  # y边际密度函数
E_x2 = integrate(x**2*p_x, (x, -oo, oo))
E_x = integrate(x*p_x, (x, -oo, oo))
E_y2 = integrate(y**2*p_y, (y,-oo,oo))
E_y = integrate(y*p_y, (y,-oo,oo))
E_xy = integrate(x*y*p_xy, (x, -oo, oo),(y,-oo,oo))
cov_xy = E_xy - E_x*E_y
var_x = E_x2 - E_x**2
var_y = E_y2 - E_y**2
Matrix([[var_x,cov_xy],[cov_xy,var_y]])

相关系数

只能表示两个随机变量之间的线性相关性，相关系数等于0并不是说明两个随机变量之间不存在任何函数关系，只是说明它们之间不存在线性关系
相关系数矩阵

# 求相关系数矩阵：1.求两两变量的相关系数；2. 组合成矩阵
from sympy import *
x = symbols('x')
y = symbols('y')
p_xy = Piecewise((1/3*(x+y),And(x>0,x<1,y>0,y<2)),(0,True))
p_x = integrate(p_xy, (y, -oo, oo))  # x边际密度函数
p_y = integrate(p_xy, (x, -oo, oo))  # y边际密度函数
E_x2 = integrate(x**2*p_x, (x, -oo, oo))
E_x = integrate(x*p_x, (x, -oo, oo))
E_y2 = integrate(y**2*p_y, (y,-oo,oo))
E_y = integrate(y*p_y, (y,-oo,oo))
E_xy = integrate(x*y*p_xy, (x, -oo, oo),(y,-oo,oo))
cov_xy = E_xy - E_x*E_y
var_x = E_x2 - E_x**2
var_y = E_y2 - E_y**2
corr_xy = cov_xy/(sqrt(var_x*var_y))
Matrix([[1,corr_xy],[corr_xy,1]])

依概率收敛

# 模拟抛硬币正面的概率是否会越来越接近0.5
import random
def Simulate_coin(test_num):
    random.seed(100)
    coin_list = [1 if (random.random()>=0.5) else 0  for i in range(test_num)]   # 模拟试验结果
    coin_frequence = np.cumsum(coin_list) / (np.arange(len(coin_list))+1)  # 计算正面为1的频率
    plt.figure(figsize=(6,4))
    plt.plot(np.arange(len(coin_list))+1, coin_frequence, c='blue', alpha=0.7)
    plt.axhline(0.5,linestyle='--',c='red',alpha=0.5)
    plt.xlabel("test_index")
    plt.ylabel("frequence")
    plt.title(str(test_num)+" times")
    plt.show()

Simulate_coin(test_num = 500)
Simulate_coin(test_num = 1000)
Simulate_coin(test_num = 5000)
Simulate_coin(test_num = 10000)

大数定律

# 蒙特卡洛积分计算的原理：
x_arr = np.linspace(0,1,1000)
x_n = uniform.rvs(size = 100)  # 随机选择n个x随机数
y_n = uniform.rvs(size = 100)  # 随机选择n个y随机数
plt.stackplot(x_arr,np.square(x_arr),alpha=0.5,color="skyblue") #堆积面积图
plt.scatter(x_n,y_n)
plt.text(1.0,1.0,r'$y=x^2$')
plt.show()

# 使用蒙特卡洛法计算y=x^2在【0，1】上的定积分
from scipy.stats import uniform
def MonteCarloRandom(n):
    x_n = uniform.rvs(size = n)  # 随机选择n个x随机数
    y_n = uniform.rvs(size = n)  # 随机选择n个y随机数
    f_x = np.square(x_n)    # 函数值f_x = x**2
    binory_y = [1.0 if y_n[i] < f_x[i] else 0 for i in range(n)]    # 如果y<x**2则为1，否则为0
    J = np.sum(binory_y) / n
    return J
    
print("y=x**2在[0,1]的定积分为：",integrate(x**2, (x,0,1)))
print("模拟10次的定积分为：",MonteCarloRandom(10))
print("模拟100次的定积分为：",MonteCarloRandom(100))
print("模拟1000次的定积分为：",MonteCarloRandom(1000))
print("模拟10000次的定积分为：",MonteCarloRandom(10000))
print("模拟100000次的定积分为：",MonteCarloRandom(100000))
print("模拟1000000次的定积分为：",MonteCarloRandom(1000000))

中心极限定理

模拟n个正态分布的和的分布

from scipy.stats import norm
def Random_Sum_F(n):
    sample_nums = 10000
    random_arr = np.zeros(sample_nums)
    for i in range(n):
        mu = 0
        sigma2 = np.random.rand()
        err_arr = norm.rvs(size=sample_nums)
        random_arr += err_arr
    plt.hist(random_arr)
    plt.title("n = "+str(n))
    plt.xlabel("x")
    plt.ylabel("p (x)")
    plt.show()

Random_Sum_F(2)
Random_Sum_F(10)
Random_Sum_F(100)
Random_Sum_F(1000)
Random_Sum_F(10000)

模拟n个均匀分布的和的分布

from scipy.stats import uniform
def Random_Sum_F(n):
    sample_nums = 10000
    random_arr = np.zeros(sample_nums)
    for i in range(n):
        err_arr = uniform.rvs(size=sample_nums)
        random_arr += err_arr
    plt.hist(random_arr)
    plt.title("n = "+str(n))
    plt.xlabel("x")
    plt.ylabel("p (x)")
    plt.show()

Random_Sum_F(2)
Random_Sum_F(10)
Random_Sum_F(100)
Random_Sum_F(1000)
Random_Sum_F(10000)

模拟n个指数分布的和的分布

from scipy.stats import expon
def Random_Sum_F(n):
    sample_nums = 10000
    random_arr = np.zeros(sample_nums)
    for i in range(n):
        err_arr = expon.rvs(size=sample_nums)
        random_arr += err_arr
    plt.hist(random_arr)
    plt.title("n = "+str(n))
    plt.xlabel("x")
    plt.ylabel("p (x)")
    plt.show()

Random_Sum_F(2)
Random_Sum_F(10)
Random_Sum_F(100)
Random_Sum_F(1000)
Random_Sum_F(10000)

模拟n个泊松分布的和的分布

from scipy.stats import poisson
def Random_Sum_F(n):
    sample_nums = 10000
    random_arr = np.zeros(sample_nums)
    for i in range(n):
        mu = 1.0
        err_arr = poisson.rvs(mu=mu,size=sample_nums)
        random_arr += err_arr
    plt.hist(random_arr)
    plt.title("n = "+str(n))
    plt.xlabel("x")
    plt.ylabel("p (x)")
    plt.show()

Random_Sum_F(2)
Random_Sum_F(10)
Random_Sum_F(100)
Random_Sum_F(1000)
Random_Sum_F(10000)

模拟n个0-1分布的和的分布

from scipy.stats import bernoulli
def Random_Sum_F(n):
    sample_nums = 10000
    random_arr = np.zeros(sample_nums)
    for i in range(n):
        p = 0.5
        err_arr = bernoulli.rvs(p=p,size=sample_nums)
        random_arr += err_arr
    plt.hist(random_arr)
    plt.title("n = "+str(n))
    plt.xlabel("x")
    plt.ylabel("p (x)")
    plt.show()

Random_Sum_F(2)
Random_Sum_F(10)
Random_Sum_F(100)
Random_Sum_F(1000)
Random_Sum_F(10000)

由均匀分布随机数产生N个正态分布的随机数

import random
def Random_Norm(N,mu,sigma):
    random_list = []
    for i in range(N):
        uniform_sum = 0
        for j in range(12):
            uniform_rand = random.random() # [0,1]均匀分布的随机数
            uniform_sum += uniform_rand
        y = uniform_sum - 6
        z = mu + sigma * y
        random_list.append(z)
    return random_list

norm_random_list = Random_Norm(10000,10,2)
plt.hist(np.array(norm_random_list))
plt.xlabel("x")
plt.ylabel("p (x)")
plt.title("由均匀分布随机数构造正态分布随机数")
plt.text(16,2500,r'$N(10,4)$')
plt.show()

人面猴
序言：七十年代末，一起剥皮案震惊了整个滨河市，随后出现的几起案子，更是在滨河造成了极大的恐慌，老刑警刘岩，带你破解...
沈念sama阅读 217,185评论 6赞 503
死咒
序言：滨河连续发生了三起死亡事件，死亡现场离奇诡异，居然都是意外死亡，警方通过查阅死者的电脑和手机，发现死者居然都...
沈念sama阅读 92,652评论 3赞 393
救了他两次的神仙让他今天三更去死
文/潘晓璐我一进店门，熙熙楼的掌柜王于贵愁眉苦脸地迎上来，“玉大人，你说我怎么就摊上这事。” “怎么了？”我有些...
开封第一讲书人阅读 163,524评论 0赞 353
道士缉凶录：失踪的卖姜人
文/不坏的土叔我叫张陵，是天一观的道长。经常有香客问我，道长，这世上最难降的妖魔是什么？我笑而不...
开封第一讲书人阅读 58,339评论 1赞 293
港岛之恋（遗憾婚礼）
正文为了忘掉前任，我火速办了婚礼，结果婚礼上，老公的妹妹穿的比我还像新娘。我一直安慰自己，他们只是感情好，可当我...
茶点故事阅读 67,387评论 6赞 391
恶毒庶女顶嫁案：这布局不是一般人想出来的
文/花漫我一把揭开白布。她就那样静静地躺着，像睡着了一般。火红的嫁衣衬着肌肤如雪。梳的纹丝不乱的头发上，一...
开封第一讲书人阅读 51,287评论 1赞 301
城市分裂传说
那天，我揣着相机与录音，去河边找鬼。笑死，一个胖子当着我的面吹牛，可吹牛的内容都是我干的。我是一名探鬼主播，决...
沈念sama阅读 40,130评论 3赞 418
双鸳鸯连环套：你想象不到人心有多黑
文/苍兰香墨我猛地睁开眼，长吁一口气：“原来是场噩梦啊……” “哼！你这毒妇竟也来了？” 一声冷哼从身侧响起，我...
开封第一讲书人阅读 38,985评论 0赞 275
万荣杀人案实录
序言：老挝万荣一对情侣失踪，失踪者是张志新（化名）和其女友刘颖，没想到半个月后，有当地人在树林里发现了一具尸体，经...
沈念sama阅读 45,420评论 1赞 313
护林员之死
正文独居荒郊野岭守林人离奇死亡，尸身上长有42处带血的脓包…… 初始之章·张勋以下内容为张勋视角年9月15日...
茶点故事阅读 37,617评论 3赞 334
白月光启示录
正文我和宋清朗相恋三年，在试婚纱的时候发现自己被绿了。大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
茶点故事阅读 39,779评论 1赞 348
活死人
序言：一个原本活蹦乱跳的男人离奇死亡，死状恐怖，灵堂内的尸体忽然破棺而出，到底是诈尸还是另有隐情，我是刑警宁泽，带...
沈念sama阅读 35,477评论 5赞 345
日本核电站爆炸内幕
正文年R本政府宣布，位于F岛的核电站，受9级特大地震影响，放射性物质发生泄漏。R本人自食恶果不足惜，却给世界环境...
茶点故事阅读 41,088评论 3赞 328
男人毒药：我在死后第九天来索命
文/蒙蒙一、第九天我趴在偏房一处隐蔽的房顶上张望。院中可真热闹，春花似锦、人声如沸。这庄子的主人今日做“春日...
开封第一讲书人阅读 31,716评论 0赞 22
一桩弑父案，背后竟有这般阴谋
文/苍兰香墨我抬头看了看天上的太阳。三九已至，却和暖如春，着一层夹袄步出监牢的瞬间，已是汗流浃背。一阵脚步声响...
开封第一讲书人阅读 32,857评论 1赞 269
情欲美人皮
我被黑心中介骗来泰国打工，没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留，地道东北人。一个月前我还...
沈念sama阅读 47,876评论 2赞 370
代替公主和亲
正文我出身青楼，却偏偏与公主长得像，于是被迫代替她去往敌国和亲。传闻我的和亲对象是个残疾皇子，可洞房花烛夜当晚...
茶点故事阅读 44,700评论 2赞 354

Task3：概率论与数理统计（基于Python）

随机现象与概率

条件概率、乘法公式、全概率公式与贝叶斯公式

已知密度函数求分布函数

已知分布函数求密度函数

常见的连续型随机变量及其密度函数

使用scipy计算密度函数画图(非自定义函数)

使用scipy计算分布函数画图(非自定义函数)

使用scipy的pmf和cdf画图

一维随机变量的数字特征：期望、方差、分位数与中位数

多维随机变量及其联合分布、边际分布、条件分布

推荐阅读更多精彩内容