第十二篇：概率论

概率论的基本概念

排列

def A(m:int,n:int):
    if m >= n > 0:
        result = 1
        for v in range(m,m-n,-1):
            result *= v
        return result
    else:
        return None
print(A(10,3))

组合

def C(m:int,n:int):
    if m >= n:
        result = 1
        factorial = 1
        for v in range(m,m-n,-1):
            result *= v
        for i in range(1,n+1):
            factorial *= i
        result /= factorial
        return result
    else:
        return None
print(C(10,3))

 古典概型

联合概率

def  joint_probability(A_pro,B_pro):
    return A_pro*B_pro

条件概率

全概率公式

贝叶斯公式

概率公式

事件的独立

from sympy import pprint,Eq,simplify,lambdify,S,And,Symbol,plot,diff,Rational,solve
from sympy.stats import P,density,Poisson,Binomial,sample_iter,Hypergeometric
from sympy.abc import p,lamda,k
#sympy使用方便，但计算效率慢一些
#用颜色的求被取出符从超几何分布
A = Hypergeometric('A',4,2,1) #共有4个球，2个有红色，现从4球中取1个。
B = Hypergeometric('B',4,2,1) #共有4个球，2个有白色，现从4球中取1个。
C = Hypergeometric('C',4,2,1) #共有4个球，2个有黑色，现从4球中取1个。
print(P(Eq(A,1)))
print(P(Eq(B,1)))
print(P(Eq(C,1)))

随机变量及其分布

 

from sympy import pprint,Eq,simplify,lambdify,S,And,Symbol,plot,diff,Rational,solve
from sympy.stats import P,density,Poisson,Binomial,sample_iter,Hypergeometric
from sympy.abc import p,lamda,k
#sympy使用方便，但计算效率慢一些
#取出的3只球中的黑球的个数服从超几何分布
A = Hypergeometric('A',5,3,3) #共有5个球，3个黑球，现从5球中任取3个。
print('取出的3只球中的黑球的个数等于3的概率：{}'.format(P(Eq(A,3))))
print(density(A).dict)

from sympy import pprint,Eq,simplify,lambdify,S,And,Symbol,plot,diff,Rational,Mod
from sympy.stats import P,density,Die,Binomial,sample_iter
from sympy.abc import p,lamda,k
#sympy使用方便，但计算效率慢一些
X = Die('X')
print(density(X).dict)
print(density(X<=4).dict)
print('X<=4的概率：{}'.format(P(X<=4)))
print('X取偶的概率：{}'.format(P(Eq(Mod(X,2),0))))

离散型随机变量及其分布规律

Bernoulli分布

from sympy import pprint,Eq,simplify,lambdify,S,collect
from sympy.stats import Coin,density,Binomial,variance,sample_iter,sample,Die,Bernoulli
from sympy.abc import p,x,k
C = Coin('C') #抛硬币
D = Die('D') #抛骰子
C_dis = density(C).dict #概率字典
D_dis = density(D)(x) #转成函数形式
print(C_dis)
print(D_dis)

#Bernoulli分布
X = Bernoulli('X',p,1,0) #随机变量X服从参数为p的Bernoulli分布，1表示成功，0表示失败
P_X = density(X)
print('Bernoulli分布律：') #分段函数里面的True表示其他的意思
pprint(P_X(k),use_unicode=False)
print('成功的概率是：{}'.format(P_X(1)))
print('所有取值的概率字典：{}'.format(P_X.dict))

二项分布

from sympy import pprint,Eq,simplify,lambdify,S,collect,Symbol,Integer
from sympy.stats import Coin,density,Binomial,variance,sample_iter,sample,Die,Bernoulli
from sympy.abc import p
k = Symbol('k',positive=True,integer=True,real=True)
X = Binomial('X',6,p,1,0) #随机变量X服从参数为(6,p)的二项分布，1表示成功，0表示失败
P_X = density(X)
print('二项分布律：') #分段函数里面的True表示其他的意思
pprint(P_X(k),use_unicode=False)
print('做6次实验有3次成功的概率是：{}'.format(P_X(3)))
print('所有取值的概率字典：{}'.format(P_X.dict))

from sympy import pprint,Eq,simplify,lambdify,S,collect,Symbol,Integer
from sympy.stats import Coin,density,Binomial,variance,sample_iter,sample,Die,Bernoulli
from sympy.abc import p
k = Symbol('k',positive=True,integer=True,real=True)
X = Binomial('X',15,0.1) #随机变量X服从参数为(6,p)的二项分布，1表示成功，0表示失败
P_X = density(X)
print('取出的15件产品中恰有2件次品的概率是：{}'.format(P_X(2)))
print('取出的15件产品中至少有2件次品的概率是：{}'.format(1-P_X(2)))

from sympy import pprint,Eq,simplify,lambdify,S,collect,Symbol,Integer,plot,diff,Rational
from sympy.stats import Coin,density,Binomial,variance,sample_iter,sample,Die,Bernoulli
from sympy.abc import p
import matplotlib.pyplot as plt
import matplotlib as mpl
import numpy as np
# 设置中文显示
mpl.rcParams['font.sans-serif'] = ['SimHei']
mpl.rcParams['axes.unicode_minus'] = False
k = Symbol('k',positive=True,integer=True,real=True)
X = Binomial('X',300,0.44) #随机变量X服从参数为(6,p)的二项分布，1表示成功，0表示失败
P_X = density(X)
x = np.array([i for i in range(301)])
y = np.array([P_X.dict[i] for i in range(301)])
plt.plot(x,y,'.b')
plt.title('二项分布的分布律图')
plt.show()
print('300次射击最可能命中{}次,其相应的概率是：{}'.format(int((300+1)*0.44),P_X(int((300+1)*0.44))))

泊松分布

 

from sympy import pprint,Eq,simplify,lambdify,S,collect,Symbol,plot,diff,Rational,solve
from sympy.stats import Coin,density,Poisson,Binomial,sample_iter,sample,Die,Bernoulli
from sympy.abc import p,lamda
# k = Symbol('k',positive=True,integer=True,real=True)
# X = Poisson('X',lamda) #随机变量X服从参数为lamda的泊松分布
# P_X = density(X)
# print('泊松分布的概率函数：')
# pprint(P_X(k),use_unicode=False)

#sympy使用方便，但计算效率慢一些
X = Poisson('X',lamda)
P_X = density(X)
lmd = solve(Eq(P_X(1),P_X(2)))
for e in lmd:
    if e > 0:
        print('lamda = {}'.format(e))
        print('X = 4 的概率是：{}'.format(density(Poisson('x',e))(4).evalf()))

#用scipy计算效率高一些
from scipy.stats import poisson
ps = poisson(2) #传入lamda参数
print('X = 4 的概率是：{}'.format(ps.pmf(4)))

from scipy.stats import poisson,binom
bn = binom(600,0.012)
print('用二项分布计算的结果：{}'.format(1-bn.pmf(0)-bn.pmf(1)-bn.pmf(2)))

ps = poisson(600*0.012) #传入lamda参数
print('用泊松分布近似计算的结果：{}'.format(1-ps.pmf(0)-ps.pmf(1)-ps.pmf(2)))

from sympy.stats import P,density,Poisson,Binomial,sample_iter,Probability,Die,Bernoulli
from sympy.abc import p,lamda
#sympy使用方便，但计算效率慢一些
X = Binomial('X',2500,0.001)
print('用二项分布计算：')
print(round(1-P(25>=2*X),6))
print(round(P(25<=2*X),6))
Y = Poisson('Y',2500*0.001)
Y_pdf = density(Y)
result = 0
for x in range(13):
    result += Y_pdf(x)
print('用泊松分布计算：')
print(round(1 - result,6))

#用scipy计算效率高一些
from scipy.stats import poisson
ps = poisson(2500*0.001) #传入lamda参数
result = 0
for x in range(13):
    result += ps.pmf(x)
print(round(1-result,6))

print('第二题：')
print('用泊松分布计算：')
print(round(P(25-2*X>=10),6))

result = 0
for x in range(8):
    result += Y_pdf(x)
print('用泊松分布计算：')
print(round(result,6))
result = 0
for x in range(8):
    result += ps.pmf(x)
print(round(result,6))

几何分布

from sympy.stats import P,density,Poisson,Binomial,sample_iter,Probability,Die,Geometric
from sympy.abc import p,lamda,k
#sympy使用方便，但计算效率慢一些
X = Geometric('X',0.1)
X_pdf = density(X)
print(X_pdf(k))

#用scipy计算效率高一些
import numpy as np
from scipy.stats import poisson,geom
Y = geom(0.1)
x = np.array([i for i in range(9)])
print(Y.pmf(x))

 超几何分布

from sympy.stats import P,density,Poisson,Binomial,sample_iter,Hypergeometric
from sympy.abc import p,lamda,k
#sympy使用方便，但计算效率慢一些
X = Hypergeometric('X',10,5,3)
X_pdf = density(X)
print(X_pdf.dict)

#用scipy计算效率高一些
import numpy as np
from scipy.stats import poisson,hypergeom
N = 10
M = 5
n = 3
Y = hypergeom(N,M,n)
x = np.array([i for i in range(min(M,n)+1)])
print({key:val for key,val in enumerate(Y.pmf(x))})

连续型随机变量及其分布规律

from sympy import pprint,Eq,E,simplify,Piecewise,S,And,solve,oo,exp,Symbol,Mod,integrate
from sympy.stats import P,density,Die,Binomial,sample_iter,FiniteRV
from sympy.abc import p,lamda,k,r,w,d,c
x = Symbol('x',real=True)
f = Piecewise((c*(4*x-2*x**2),And(0<x,x<2)),(0,True))

print(solve(Eq(integrate(f,(x,0,2)),1),c))
print(integrate(f,(x,1,2)).subs(c,3/8))

from sympy import pprint,Eq,E,simplify,lambdify,S,And,solve,oo,exp,Rational,Mod,integrate
from sympy.stats import P,density,Die,Binomial,sample_iter,FiniteRV
from sympy.abc import p,lamda,k,r,w,d,x
f = k*E**(-3*x)
# 一个数与无穷大的做加减乘除还是等于无穷大。
print(integrate(f,(x,0,+oo)))
print(solve(Eq(integrate(f,(x,0,+oo)),1),k)) # -k*exp(-oo)/3 - (-k*exp(0) = 1；k = 3
#不定积分是要加一个常数C的，并且 -1 <= -exp(-3x) <= 0；0 <= F(x) <= 1，所以C=1
print(integrate(f,x).subs(k,3))

均匀分布

from sympy import pprint,Eq,E,simplify,Piecewise,S,And,solve,oo,exp,Symbol,Mod,integrate
from sympy.stats import P,density,Die,Uniform,sample_iter,FiniteRV,cdf
from sympy.abc import p,lamda,k,r,w
x = Symbol('x',real=True)
a = Symbol("a", negative=True)
b = Symbol("b", positive=True)
X = Uniform('X',a,b)
print('均匀分布的概率密度函数：')
pprint(density(X)(x),use_unicode=False)
print('均匀分布的分布函数：')
pprint(simplify(cdf(X)(x)),use_unicode=False)

from sympy import pprint,Eq,E,simplify,Piecewise,S,And,solve,radsimp,exp,Symbol,Mod,integrate
from sympy.stats import P,density,Die,Uniform,sample_iter,FiniteRV,cdf
from sympy.abc import p,lamda,k,r,w
x = Symbol('x',real=True)
X = Uniform('X',0,30)
print('均匀分布的概率密度函数：')
pprint(density(X)(x),use_unicode=False)
print('均匀分布的分布函数：')
pprint(cdf(X)(x),use_unicode=False)
print(P(And(X>=25,X<=30))+P(And(X>=10,X<=15)))
print('方法2：')
X_cdf = cdf(X)
print((X_cdf(15)-X_cdf(10))+(X_cdf(30)-X_cdf(25)))

from scipy.stats import poisson,uniform
Y = uniform(0,30)
# x = np.array([i for i in range(min(M,n)+1)])
# print({key:val for key,val in enumerate(Y.pmf(x))})
#Y.pdf()：概率密度函数，Y.cdf分布函数
print((Y.cdf(15)-Y.cdf(10))+(Y.cdf(30)-Y.cdf(25)))

指数分布

from sympy import pprint,Eq,E,simplify,Piecewise,S,And,solve,radsimp,exp,Symbol,Mod,integrate
from sympy.stats import P,density,Die,Exponential,sample_iter,FiniteRV,cdf
from sympy.abc import p,k,r,x
lamda = Symbol('lamda',positive=True)
X = Exponential('X',lamda)
print('指数分布的概率密度函数：')
pprint(density(X)(x),use_unicode=False)
print('指数分布的分布函数：')
pprint(cdf(X)(x),use_unicode=False)

from sympy import pprint,Eq,E,simplify,Piecewise,S,And,solve,radsimp,exp,Symbol,Mod,integrate
from sympy.stats import P,density,Die,Exponential,sample_iter,FiniteRV,cdf
from sympy.abc import p,k,r,x
lamda = Symbol('lamda',positive=True)
X = Exponential('X',1/10)
print(P(And(X>=10,X<=20)))
print('方法2：')
X_cdf = cdf(X)
print(X_cdf(20)-X_cdf(10))

from scipy.stats import poisson,expon
Y = expon(10)
#Y.pdf()：概率密度函数，Y.cdf分布函数
print(Y.cdf(20)-Y.cdf(10))

正态分布

from sympy import pprint,Eq,E,simplify,Piecewise,S,And,solve,radsimp,exp,Symbol,Mod,integrate
from sympy.stats import P,density,Die,Normal,sample_iter,FiniteRV,cdf
from sympy.abc import p,k,r,x
mu = Symbol("mu")
sigma = Symbol("sigma", positive=True)
X = Normal('X',mu,sigma)
print('正态分布的概率密度函数：')
pprint(density(X)(x),use_unicode=False)
print('正态分布的分布函数：')
pprint(simplify(cdf(X)(x)),use_unicode=False)

from sympy import pprint,Eq,E,simplify,Piecewise,S,And,solve,radsimp,exp,Symbol,Mod,integrate
from sympy.stats import P,density,Die,Normal,sample_iter,FiniteRV,cdf
from sympy.abc import p,k,r,x
X = Normal('X',0,1)
print('正态分布的概率密度函数：')
pprint(density(X)(x),use_unicode=False)
print('正态分布的分布函数：')
pprint(simplify(cdf(X)(x)),use_unicode=False)

from sympy import pprint,Eq,E,simplify,Piecewise,S,And,plot,radsimp,exp,Symbol,Mod,integrate
from sympy.stats import P,density,Die,Normal,sample_iter,FiniteRV,cdf
from sympy.abc import p,k,r,x
X = Normal('X',0,1)
Y = Normal('Y',2,1)
X_pdf = density(X)
Y_pdf = density(Y)
p1 = plot(X_pdf(x),(x,-10,10),line_color=(1.0,0.0,0.0),show=False)
p2 = plot(Y_pdf(x),(x,-10,10),line_color=(0.0,1.0,0.0),show=False)
p2.append(p1[0])
p2.show()

from sympy import pprint,Eq,E,simplify,Piecewise,S,And,plot,radsimp,exp,Symbol,Mod,integrate
from sympy.stats import P,density,Die,Normal,sample_iter,FiniteRV,cdf
from sympy.abc import p,k,r,x
X = Normal('X',2,1)
Y = Normal('Y',2,5)
Z = Normal('Z',2,9)
X_pdf = density(X)
Y_pdf = density(Y)
Z_pdf = density(Z)
p1 = plot(X_pdf(x),(x,-10,10),line_color=(1.0,0.0,0.0),show=False)
p2 = plot(Y_pdf(x),(x,-10,10),line_color=(0.0,1.0,0.0),show=False)
p3 = plot(Z_pdf(x),(x,-10,10),line_color=(0.0,0.0,1.0),show=False)
p2.append(p1[0])
p2.append(p3[0])
p2.show()

from sympy import pprint,Eq,E,simplify,Piecewise,S,And,plot,radsimp,exp,Symbol,Mod,integrate
from sympy.stats import P,density,Die,Normal,sample_iter,FiniteRV,cdf
from sympy.abc import p,k,r,x
X = Normal('X',0,1)
X_pdf = density(X)
print('(1)：{}'.format(P(And(X>=1,X<2)).evalf()))
print('(2)：{}'.format(P(And(X>=-1,X<2)).evalf()))

from scipy.stats import norm
Y = norm(0,1)
#Y.pdf()：概率密度函数，Y.cdf分布函数
print('(1)：{}'.format(Y.cdf(2)-Y.cdf(1)))
print('(2)：{}'.format(Y.cdf(2)-Y.cdf(-1)))

Γ - 分布

from sympy import pprint,Eq,E,simplify,Piecewise,S,And,plot,radsimp,exp,Symbol,Mod,integrate
from sympy.stats import P,density,Die,Gamma,sample_iter,FiniteRV,cdf
from sympy.abc import p,x,lamda,r
X = Gamma('X',r,1/lamda)
X_pdf = density(X)
print(density(X)(x))

from sympy.stats import P,density,Die,Gamma,sample_iter,FiniteRV,cdf
from sympy.abc import p,x,lamda
X = Gamma('X',1,1/lamda)
X_pdf = density(X)
print(density(X)(x))

from sympy import pprint,Eq,E,simplify,Piecewise,S,And,plot,radsimp,exp,Symbol,Mod,integrate
from sympy.stats import P,density,variance,Gamma,sample_iter,E,cdf
from sympy.abc import p,x,lamda,r
print('*'*18+'使用sympy计算'+'*'*18)
X = Gamma('X',1,2)
X_pdf = density(X)
print(density(X)(x))
print('(1)：{}'.format((cdf(X)(5)-cdf(X)(3)).evalf()))
print(E(X))
print(simplify(variance(X)))


print('*'*18+'使用scipy计算'+'*'*18)
'''
scipy中连续随机变量的主要公共方法如下：
rvs:随机变量（就是从这个分布中抽一些样本）
pdf：概率密度函数。
cdf：累计分布函数
sf：残存函数（1-CDF）
ppf：分位点函数（CDF的逆）
isf：逆残存函数（sf的逆）
stats:返回均值，方差，（费舍尔）偏态，（费舍尔）峰度。
moment:分布的非中心矩。
'''
from scipy.stats import gamma
Y = gamma(1,scale=2) #lamda参数需要通过scale给定。
#Y.stats()：返回均值，方差，（费舍尔）偏态，（费舍尔）峰度。moments="mv"表示只返回均值和方差。
print(Y.stats(moments="mv"))
print('(2)：{}'.format(Y.cdf(5)-Y.cdf(3)))

贝塔分布

from sympy import pprint,Eq,E,simplify,Piecewise,S,And,plot,radsimp,exp,Symbol,Mod,integrate
from sympy.stats import P,density,Die,Beta,sample_iter,FiniteRV,cdf
from sympy.abc import p,x,lamda,r
X = Beta('X',1,1/2)
X_pdf = density(X)
print(density(X)(x))
print('(1)：{}'.format((cdf(X)(1/2)-cdf(X)(1/3)).evalf()))

from scipy.stats import beta
Y = beta(1,1/2)
#Y.pdf()：概率密度函数，Y.cdf分布函数
print('(2)：{}'.format(Y.cdf(1/2)-Y.cdf(1/3)))

二维随机变量及其分布函数

 

 

 二维随机变量的分布函数

二维离散型随机变量

from sympy import pprint,Eq,E,simplify,Piecewise,S,And,plot,radsimp,exp,Symbol,Mod,integrate
from sympy.stats import P,density,variance,DiscreteUniform,sample_iter,E,cdf,given
from sympy.abc import p,x,lamda,r
from sympy.matrices import Matrix
print('*'*18+'使用sympy计算'+'*'*18)
X = DiscreteUniform('X',list(range(1,5)))
Y = DiscreteUniform('Y',list(range(1,5)))
# for i in range(1,5):
#     Y = given(X,X<=i) #创建某个条件下的新的样本空间。
#     print(density(Y).dict)
result = []
for v_x,p_x in density(X).dict.items():
    newlist = []
    for v_y, p_y in density(X).dict.items():
        newlist.append(p_x*P(Eq(Y,v_y),Y<=v_x)) #由公式：P(A,B) = P(A)*P(B/A) 得。
    result.append(newlist)
pprint(Matrix(result).T)
print('*'*18+'使用scipy计算'+'*'*18)
from scipy.stats import rv_discrete
import pandas as pd
import numpy as np
#已经知道X是服从离散均匀分布的
# xd = [i for i in range(1,5)]
# xp = [1/4 for i in range(1,5)]
# Z = rv_discrete(name='Z',values=(xd,xp)) #lamda参数需要通过scale给定。
# # for i in Z.numargs:
# print(Z.pmf(3))
ind_clu = [i for i in range(1,5)]
date = pd.DataFrame(np.array(result).T,ind_clu,ind_clu)
date.index.name='Y'
date.columns.name = 'X'
print(date)

二维连续型随机变量

from sympy import pprint,Eq,E,simplify,Piecewise,S,And,oo,radsimp,exp,Symbol,Mod,integrate
from sympy.stats import P,density,variance,ContinuousRV,sample_iter,E,cdf,given
from sympy.abc import p,lamda,r,y,x
from sympy.matrices import Matrix
print('*'*18+'使用sympy计算'+'*'*18)
pdf = Piecewise((2*exp(-(2*x+y)),And(x>0,y>0)),(0,True))
def my2Cdf(pdf,x1,x2,y1,y2):
    return integrate(pdf,(x,x1,x2),(y,y1,y2))
    # return integrate(integrate(pdf,(x,x1,x2)),(y,y1,y2))
#由 0 < y <= x < +oo 得；y <= x < +oo , 0 < y < +oo
print(my2Cdf(2*exp(-(2*x+y)),y,+oo,0,+oo))

边缘分布函数

离散型随机变量的边缘分布

import pandas as pd
import numpy as np
from fractions import Fraction #python的分数运算库
index_columns = [0,1]
data = pd.DataFrame(np.array([[Fraction(16,49),Fraction(12,49)],[Fraction(12,49),Fraction(9,49)]]),index_columns,index_columns)
data['col_sum'] = data.apply(lambda x:x.sum(),axis=1)
data.loc['row_sum'] = data.apply(lambda x:x.sum())
data.index.name='Y'
data.columns.name = 'X'
print(data)

连续性随机变量的边缘分布

离散型随机变量的条件分布

连续性随机变量的条件分布

数字特征

数学期望

from sympy import pprint,Eq,E,simplify,Piecewise,S,And,oo,radsimp,exp,Symbol,integrate,Integral
from sympy.stats import P,density,variance,Normal,Expectation,E,cdf,Die,Exponential,Probability
from sympy.abc import p,lamda,r,y,x,mu,sigma,C
from sympy.matrices import Matrix
X = Normal('X',mu,sigma)
Y = Exponential('Y',lamda)
Z = Die('Z',10)
print('*'*18+'使用sympy计算'+'*'*18)
print('E(C) = {}'.format(Expectation(C)))
print('E(X) = {}'.format(Expectation(X).rewrite(Integral)))
print('E(X)：')
pprint(Expectation(X).rewrite(Probability))
print('E(C*X) = {}'.format(Expectation(C*X).doit()))
print('E(X+Y) = {}'.format(Expectation(X+Y).doit()))

from sympy.stats import P,density,variance,Normal,Expectation,E,cdf,Die,Exponential
from sympy.abc import p,lamda,r,y,x,mu,sigma
from sympy.matrices import Matrix
X = Normal('X',2,1)
Y = Exponential('Y',1/2)
Z = Die('Z',10)
W = X+Y
C = 3
print('*'*18+'使用sympy计算'+'*'*18)
print('E(C)：{}'.format(E(C)))
print('E(X)：{}'.format(E(X)))
print('E(Y)：{}'.format(E(Y)))
print('E(Z)：{}'.format(E(Z)))
print('E(X*Y)：{}；E(X)*E(Y)：{}'.format(E(X*Y),E(X)*E(Y)))
print('E(W)：{}'.format(E(W)))
print('E(X*W)：{}；E(X)*E(W)：{}'.format(E(X*W),E(X)*E(W)))

方差

from sympy import pprint,Eq,E,simplify,Piecewise,S,And,oo,radsimp,exp,Symbol,integrate,Integral
from sympy.stats import P,density,variance,Normal,Exponential,cdf,Die,Variance,Expectation
from sympy.abc import p,lamda,r,y,x,mu,sigma,C
from sympy.matrices import Matrix
X = Normal('X',mu,sigma)
Y = Exponential('Y',lamda)
Z = Die('Z',10)
W = X+Y
print('*'*18+'使用sympy计算'+'*'*18)
print('D(C) = {}'.format(Variance(C)))
print('D(X) = {}'.format(Variance(X).rewrite(Integral)))
print('D(X) = {}'.format(Variance(X).rewrite(Expectation)))
print('D(C*X) = {}'.format(Variance(C*X).doit()))
print('D(C+X) = {}'.format(Variance(C+X).doit()))
print('D(X+Y) = {}'.format(Variance(X+Y).doit()))
print('D(X-Y) = {}'.format(Variance(X-Y).doit()))

from sympy import pprint,Eq,E,simplify,Piecewise,S,And,oo,radsimp,exp,Symbol,integrate,Integral
from sympy.stats import P,density,variance,Normal,Exponential,cdf,Die,Variance,Expectation,Uniform,Poisson
from sympy.abc import p,r,y,x,C
from sympy.matrices import Matrix
mu = Symbol('mu',positive=True)
sigma = Symbol('sigma',positive=True)
lamda = Symbol('lamda',positive=True)
theta = Symbol('theta',positive=True)
a = Symbol('a',negative=True)
b = Symbol('b',positive=True)
X = Normal('X',mu,sigma)
Y = Exponential('Y',1/theta)
Z = Uniform('Z',a,b)
W = Poisson('W',lamda)

print('正态分布：')
print('期望：{}'.format(Expectation(X).evaluate_integral().simplify()))
print('方差：{}'.format(Variance(X).evaluate_integral().simplify()))
print()
print('指数分布：')
print('期望：{}'.format(Expectation(Y).evaluate_integral().simplify()))
print('方差：{}'.format(Variance(Y).evaluate_integral().simplify()))
print()
print('均匀分布：')
print('期望：{}'.format(Expectation(Z).evaluate_integral().simplify()))
print('方差：{}'.format(Variance(Z).evaluate_integral().simplify()))
print()
print('泊松分布：')
print('期望：{}'.format(Expectation(W).evaluate_integral().simplify()))
print('方差：{}'.format(Variance(W).evaluate_integral().simplify()))

标准差

协方差

from sympy import pprint,Eq,E,simplify,Piecewise,S,And,oo,radsimp,exp,Symbol,integrate,Integral
from sympy.stats import P,density,variance,Normal,Exponential,cdf,Die,Covariance,Expectation
from sympy.abc import p,r,y,x,C
from sympy.matrices import Matrix
mu = Symbol('mu',positive=True)
sigma = Symbol('sigma',positive=True)
theta = Symbol('theta',positive=True)
a = Symbol('a',positive=True)
b = Symbol('b',positive=True)
X = Normal('X',mu,sigma)
Y = Exponential('Y',1/theta)
Z = Normal('Z',mu,sigma)
print('*'*18+'使用sympy计算'+'*'*18)
print('Cov(X,Y) = {}'.format(Covariance(X,Y)))
print('Cov(X,Y) = {}'.format(Covariance(X,Y).rewrite(Expectation)))
print('Cov(C*X,Y) = {}'.format(Covariance(C*X,Y).doit()))
print('Cov(a*X,b*Y) = {}'.format(Covariance(a*X,b*Y).doit()))
print('Cov(X+Y,Z) = {}'.format(Covariance(X+Y,Z).doit()))
print('Cov(X-Y,Z) = {}'.format(Covariance(X-Y,Z).doit()))

协方差矩阵

 Pearson相关系数

import numpy as np
import pandas as pd
df = pd.DataFrame({'A':np.random.randint(1, 100, 10), 'B':np.random.randint(1, 100, 10),'C':np.random.randint(1, 100, 10)})
print(df.corr())  # pearson相关系数
df.corr('kendall') # Kendall Tau相关系数
df.corr('spearman') # spearman秩相关

 from numpy import array, cov, corrcoef
  
  data = array([data1, data2])
  
  #计算两组数的协方差
  #参数bias=1表示结果需要除以N，否则只计算了分子部分
  #返回结果为矩阵，第i行第j列的数据表示第i组数与第j组数的协方差。对角线为方差
  cov(data, bias=1)
  
  #计算两组数的相关系数
  #返回结果为矩阵，第i行第j列的数据表示第i组数与第j组数的相关系数。对角线为1
  corrcoef(data)

中心距、原点距

峰度

import numpy as np
from scipy import stats #scipy中的stats可以做统计推断
np.set_printoptions(suppress=True) #将数字表示成我们习惯的方式
x = np.random.randn(100) #randn可随机生成标准正态分布的样本，100表示样本数量
mu = np.mean(x, axis=0) #axis=0，那么输出矩阵是1行，求每一列的平均；axis=1，输出矩阵是1列，求每一行的平均
sigma = np.std(x, axis=0) #求标准差，这里除的是N
skew = stats.skew(x) #求偏度
kurtosis = stats.kurtosis(x) #求峰度

print(x,'\n')
print(mu,sigma,skew,kurtosis)

峰度反应的是图像的尖锐程度：峰度越大，表现在图像上面是中心点越尖锐。在相同方差的情况下，中间一大部分的值方差都很小，为了达到和正太分布方差相同的目的，必须有一些值离中心点越远，所以这就是所说的“厚尾”，反应的是异常点增多这一现象。峰度定义为四阶标准矩，峰度表示分布的尾部与正态分布的区别。使用峰度可帮助您初步了解有关数据分布的一般特征。

基线：峰度值 0

完全服从正态分布的数据的峰度值为 0。正态分布的数据为峰度建立了基准。如果样本的峰度值显著偏离 0，则表明数据不服从正态分布。

 

正峰度

具有正峰度值的分布表明，相比于正态分布，该分布有更重的尾部。例如，服从 t 分布的数据具有正峰度值。实线表示正态分布，虚线表示具有正峰度值的分布。

负峰度

具有负峰度值的分布表明，相比于正态分布，该分布有更轻的尾部。例如，服从 Beta 分布（第一个和第二个分布形状参数等于 2）的数据具有负峰度值。实线表示正态分布，虚线表示具有负峰度值的分布。

样本的峰度计算方法：

 样本的峰度还可以这样计算：

 

其中k4k4是四阶累积量的唯一对称无偏估计，k2k2是二阶累积量的无偏估计（等同于样本方差），m4m4是样本四阶平均距，m2m2是样本二阶平均距。

注意，大多数程序都是采用G2来计算峰度。

偏度

偏度是数据的不对称程度。无论偏度值是 0、正数还是负数，都显示有关数据分布形状的信息。

    

对称或非偏斜分布

当数据变得更加对称时，它的偏度值会更接近零。图 A 显示正态分布的数据，顾名思义，正态分布数据的偏度相对较小。通过沿这一正态数据直方图的中间绘制一条线，可以很容易地看到两侧互相构成镜像。但是，没有偏度并不表示具有正态性。在图 B 显示的分布中，两侧依然互相构成镜像，但这些数据完全不是正态分布。

 

正偏斜或向右偏斜分布

正偏斜或右偏斜的数据之所以这样命名，是因为分布的“尾部”指向右侧，而且它的偏度值大于 0（或为正数）。薪金数据通常按这种方式偏斜：一家公司中许多员工的薪金相对较低，而少数人员的薪金则非常高。

负偏斜或向左偏斜分布

左偏斜或负偏斜的数据之所以这样命名，是因为分布的“尾部”指向左侧，而且它产生负数偏度值。故障率数据通常就是左偏斜的。以灯泡为例：极少数灯泡会立即就烧坏，但大部分灯泡都会持续相当长的时间。

偏度的定义：

 

样本X的偏度为样本的三阶标准矩，其中μμ是均值，δδ为标准差，E是均值操作。μ3μ3是三阶中心距，κtκt是tthtth累积量

偏度可以由三阶原点矩来进行表示：

样本偏度的计算方法：

 一个容量为n的数据，一个典型的偏度计算方法如下：

其中x¯x¯为样本的均值（和μμ的区别是，μμ是整体的均值，x¯x¯为样本的均值）。s是样本的标准差，m3m3是样本的3阶中心距。

另外一种定义如下：

k3k3是三阶累积量κ3κ3的唯一对称无偏估计(unique symmetric unbiased estimator)（k3k3 和 κ3κ3写法不一样）。k2=s2k2=s2是二阶累积量的对称无偏估计。大多数软件当中使用G1来计算skew，如Excel，Minitab，SAS和SPSS。

import pandas as pd
x = [53, 61, 49, 66, 78, 47]
s = pd.Series(x)
print(s.skew())
print(s.kurt())

变异系数

数据的发散程度可用极差（PTP）、方差（Variance）、标准差（STD）、变异系数（CV）来衡量，它们的计算方法如下：

  from numpy import mean, ptp, var, std
  
  #极差
  ptp(data)
  #方差
  var(data)
  #标准差
  std(data)
  #变异系数
  mean(data) / std(data)

分位数

from numpy import quantile
quantile(data,q=1/2)

from scipy.stats import norm
X = norm(0,1)
p = 1/4
print(X.cdf(0))
print('下分位:')
print(X.ppf(p))
print(X.sf(0))
print('上分位:')
print(X.isf(1-p))

中位数

对于定量数据（Data）来说，均值是总和除以总量（N），中位数是数值大小位于中间（奇偶总量处理不同）的值：

 

均值相对中位数来说，包含的信息量更大，但是容易受异常的影响。使用NumPy计算均值与中位数：

from numpy import mean, median
from scipy.stats import mode
#计算众数
mode(data)
#计算中位数
median(data)

数据的发散程度可用极差（PTP）、方差（Variance）、标准差（STD）、变异系数（CV）来衡量，它们的计算方法如下：