Python学习-Scipy库统计操作(随机变量、概率密度、累积分布密度、期望、方差、描述性统计(最大最小值、均值、方差、偏差、峰度)、核密度估计)

Python学习-Scipy库统计操作

目录

1、正态连续随机变量：norm

2、概率密度函数norm.pdf()

3、累积分布函数norm.cdf()

4、统计随机变量的期望值和方差stats()

5、描述性统计函数 stat.describe()，求最大最小值、均值、方差、偏差、峰度

6、核密度估计(单变量估计stats.gaussian_kde, 多变量估计)

导入库

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

1、正态连续随机变量：norm

抽取随机变量函数norm.rvs()
参数说明：
loc: 指定平均值；
scale: 指定标准偏差；
size: 大小，整数或元组

print(norm.rvs(size=10))
print(norm.rvs(size=(2, 3)))

输出

[ 1.36819314 -0.6923591  -0.46103573  0.52013108 -0.04312324  
0.46946896  1.07127949  1.00747313  0.84670697 -0.69049129]
[[1.68036634 0.10630058 0.09252581]
 [1.67264465 1.81589265 1.55676792]]

2、概率密度函数norm.pdf()

给定随机变量的x处的概率密度函数
参数说明：
x: 指定位置；
loc: 指定平均值；
scale: 指定标准偏差

print(norm.pdf(20, 20, 10))  # 求正态连续随机变量在20处的概率密度值
print(norm.pdf(20, [20, 10, 30, 40, 20]))

输出

0.03989422804014327
[3.98942280e-01 7.69459863e-23 7.69459863e-23 5.52094836e-88
 3.98942280e-01]

3、累积分布函数norm.cdf()：pdf的积分

print(norm.cdf(20, 20, 10))  # 求正态连续随机变量在20处的累积概率密度值
print(norm.cdf(20, [20, 10, 30, 40, 20]))

输出

0.5
[5.00000000e-01 1.00000000e+00 7.61985302e-24 2.75362412e-89
 5.00000000e-01]

4、统计随机变量的期望值和方差stats()

参数说明：
loc: 指定平均值；
scale: 指定标准偏差
moments: m均值；v方差；s费舍尔偏差；k费舍尔峰度

print(norm.stats(moments='mv'))
print(norm.stats(20, moments='mvsk'))
print(norm.stats(moments='mvsk', loc=9, scale=1))

输出

(array(0.), array(1.))
(array(20.), array(1.), array(0.), array(0.))
(array(9.), array(1.), array(0.), array(0.))

5、描述性统计函数 stat.describe()

参数描述：
a: 样本数据，数组对象
axis: 指定数组的统计轴，值为整型或None，None计算整个数组
ddof: 三角*度（仅限于方差），整型，默认1
bias: False，校正偏度和峰度计算，统计偏差
nan_policy: {‘propagate’, ‘raise’,‘omit’}, ‘propagate’：返回nan；‘raise’：抛出错误；‘omit’：忽略nan值
返回值：
nobs: 观察次数（沿轴的数据长度）
minmax: 最大最小值
mean: 均值
variance: 方差
skewness: 偏差
kurtosis: 峰度

设置数据集

plt.rc('font', family='simhei', size=15)  # 设置中文显示，字体大小
plt.rc('axes', unicode_minus=False)  # 该参数解决负号显示的问题
np.random.seed(1975)
x = stats.t.rvs(df=10, size=1000)  # 产生1000个*度为10的t分布随机变量
plt.hist(x=x, bins=40, range=None, density=False, weights=None, cumulative=False, bottom=None, \
         histtype='stepfilled', align='mid', alpha=0.9, orientation='vertical', rwidth=None, log=False, \
         color='g', stacked=True, edgecolor='black')
plt.xlabel('概率分布区间')
plt.ylabel('频数or频率')
plt.title('分布直方图')
plt.show()

输出

print('min: ', x.min())  # 最小值
print('max: ', x.max())  # 最大值
print('mean: ', x.mean())  # 均值
print('variance: ', x.var())  # 方差

m, v, s, k = stats.t.stats(df=10, moments='mvsk')
n, (smin, smax), sm, sv, ss, sk = stats.describe(x)

s1 = '%-14s 均值 = %6.4f，方差 = %6.4f，费舍尔偏差 = %6.4f，费舍尔峰度 = %6.4f'
print(s1 % ('连续分布：', m, v, s, k))
print(s1 % ('离散样本：', sm, sv, ss, sk))

输出

min:  -4.0181014966326885
max:  6.2240174880715795
mean:  -0.06367299923581096
variance:  1.3075051415861059
连续分布：          均值 = 0.0000，方差 = 1.2500，费舍尔偏差 = 0.0000，费舍尔峰度 = 1.0000
离散样本：          均值 = -0.0637，方差 = 1.3088，费舍尔偏差 = 0.1155，费舍尔峰度 = 1.0614

6、核密度估计

核密度估计是在概率论中用来估计未知的密度函数，属于非参数检验方法之一。可以用于预测地理空间点数据的分布规律，或用于金融领域的基于密度的预测。

1）单变量估计

X = np.array([-8, -6, 0, 3, 5, 7], dtype=np.float)
X_plot = np.linspace(-10, 10, 100)  # 对于单变量要使用1D array
kde1 = stats.gaussian_kde(dataset=X, bw_method='scott')  # 高斯核密度估计
kde2 = stats.gaussian_kde(dataset=X, bw_method='silverman')  # 高斯核密度估计

plt.rc('font', family='simhei', size=15)  # 设置中文显示，字体大小
plt.rc('axes', unicode_minus=False)  # 该参数解决负号显示的问题

fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(X, np.zeros(X.shape), 'b+', ms=20)
ax.plot(X_plot, kde1(X_plot), 'k--', label='Scott')
ax.plot(X_plot, kde2(X_plot), 'g-', label='silverman')
plt.title('gaussian_kde')
plt.legend()
plt.show()

输出

2）多变量估计样例

# 2）多变量估计样例


def measure(n):
    m1 = np.random.normal(size=n)
    m2 = np.random.normal(scale=0.5, size=n)
    return m1+m2, m1-m2


m1, m2 = measure(2000)
xmin = m1.min()
xmax = m1.max()
ymin = m2.min()
ymax = m2.max()
X, Y = np.mgrid[xmin:xmax:100j, ymin:ymax:100j]
positions = np.vstack([X.ravel(), Y.ravel()])
values = np.vstack([m1, m2])
kernel = stats.gaussian_kde(values)
z = np.reshape(kernel.evaluate(positions).T, X.shape)
fig = plt.figure(figsize=(8, 6))
ax = fig.add_subplot(111)
ax.imshow(np.rot90(z), cmap=None, extent=[xmin, xmax, ymin, ymax])

ax.plot(m1, m2, 'k.', markersize=2)
ax.set_xlim([xmin, xmax])
ax.set_ylim([ymin, ymax])
plt.show()

输出