import scipy.stats
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
%config InlineBackend.figure_format = 'retina'
## 生成随机数
np.random.seed(123)
np.random.random(5)
array([ 0.69646919, 0.28613933, 0.22685145, 0.55131477, 0.71946897])
np.random.randint(0,9,10)
#生成0-9的10个随机整数
array([6, 1, 0, 1, 0, 0, 3, 4, 0, 0])
蒙特卡洛模拟求圆周率
num = 1000
x = np.random.random(num)
y = np.random.random(num)
pi = np.sum(x**2 + y**2 < 1) / num * 4
print('PI:', pi)
## 为什么呢?
PI: 3.196
plt.figure(figsize=(5,5))
plt.scatter(x,y, alpha = 0.6)
plt.axis([0,1,0,1])
x2 = np.arange(0, 1.01, 0.01)
y2 = np.sqrt(1 - x2**2)
plt.plot(x2, y2, 'm', lw=3)
plt.show()
连续分布和正太分布
模拟面包重量的分布
假设是均值为950克,标准差为50克的正态分布。
那个数学家的故事,想起来了
mean = 950
std = 50
# 生成满足正态分布的随机数,并绘制直方图
sample = np.random.normal(mean, std, size=365)
plt.hist(sample, bins=30, alpha=0.7, rwidth=0.9, normed=True)
plt.show()
mean = 950
std = 50
norm = scipy.stats.norm(mean, std)
# 用 `scipy.stats.norm` 生成正态分布,
##这是新的方法咯?
x = np.arange(700, 1200, 1)
## 这里的1,是什么参数
y = norm.pdf(x)
plt.plot(x, y)
plt.show()
### 绘制正太分布的概率密度函数
### 累计分布函数cdf,对了到底这几个pdf和cdf的英文全名是什么啊?
# probability density function
# cumulative distribution function
x = np.arange(700, 1200, 1)
y = norm.cdf(x)
plt.plot(x, y)
plt.show()
### 计算买到的面包小于1000克的概率
# 这个好像很重要呢
#绘制PDF曲线
x = np.arange(700, 1200, 1)
y = norm.pdf(x)
plt.plot(x, y)
#在1000处绘制竖线
plt.vlines(1000, 0, norm.pdf(1000))
#填充颜色
x2 = np.arange(700, 1000, 1)
y2 = norm.pdf(x2)
plt.fill_between(x2, y2, color='blue', alpha=0.1)
#设置y轴范围
plt.ylim(0,0.0085)
plt.show()
norm.cdf(1000)
## 其实累计分布才能解决,低于1000的问题
0.84134474606854293
计算买到的面包大于1000克的概率
#绘制PDF曲线
x = np.arange(700, 1200, 1)
y = norm.pdf(x)
plt.plot(x, y)
#绘制竖线
plt.vlines(1000, 0, norm.pdf(1000))
#填充颜色
x2 = np.arange(1000, 1200, 1)
y2 = norm.pdf(x2)
plt.fill_between(x2, y2, color='blue', alpha=0.1)
#设置y轴范围
plt.ylim(0,0.0085)
plt.show()
1 - norm.cdf(1000)
## 对呀,用1-另一半就可以了呀
0.15865525393145707
norm.sf(1000)
### 这个sf是什么意思
## 好像sf是反函数。但忘了 具体是反什么了?
0.15865525393145707
## 计算面包在950-1050的概率
#绘制PDF曲线
x = np.arange(700, 1200, 1)
y = norm.pdf(x)
plt.plot(x, y)
#绘制竖线
plt.vlines(950, 0, norm.pdf(950))
plt.vlines(1050, 0, norm.pdf(1050))
#填充颜色
x2 = np.arange(950, 1050, 1)
y2 = norm.pdf(x2)
plt.fill_between(x2, y2, color='blue', alpha=0.1)
#设置y轴范围
plt.ylim(0,0.0085)
plt.show()
norm.cdf(1050) - norm.cdf(950)
#cdf(1050)是指低于1050的函数的意思么?
0.47724986805182079
# 90%的情况下,买到的面包是小于多少克的?
### ppf才是反函数吧,就是要实现某概率,该是什么初始变量
norm.ppf(0.9)
1014.07757827723
# 80%的情况下,买到的面包是大于多少克的?
norm.isf(0.8)
907.91893832135429
离散分布
outcome = np.random.randint(0,2,10)
outcome
array([1, 0, 0, 1, 1, 0, 1, 1, 0, 1])
np.sum(outcome)
6
sample = [np.sum(np.random.randint(0,2,10)) for i in range(10000)]
sample = pd.Series(sample)
# pdserie 什么意思啊
sample.value_counts().sort_index().plot.bar()
plt.show()
# 投硬币问题的二项分布
n = 10
p = 0.5
binomial = scipy.stats.binom(n, p)
x = np.arange(0,11)
plt.plot(x, binomial.pmf(x), 'bo')
plt.vlines(x, 0, binomial.pmf(x), colors='b')
plt.ylim(0,0.3)
plt.show()
mean, var = binomial.stats()
print(mean)
print(var)
5.0
2.5
### 应用
# 某家风投企业,投资成功的概率是5%,如果它投了100个项目,恰好有5个成功的概率是多少?
n = 100
p = 0.05
binom = scipy.stats.binom(n,p)
x = np.arange(0,101)
plt.plot(x, binom.pmf(x), 'bo')
plt.vlines(x, 0, binom.pmf(x), colors='b')
plt.ylim(0,0.2)
plt.show()
binom.pmf(5)
## pmf概率质量函数
0.18001782727043672
1 - binom.cdf(4)
## 累计概率分布
0.56401869931428927
binom.sf(4)
0.56401869931429105
binom.isf(0.1)
#10%的情况下,能成功多少
8.0
binom.ppf(1 - 0.1)
8.0
## 离散分布 - 泊松分布
# 有一家便利店使用泊松分布来估计周五晚上到店买东西的顾客数,根据以往数据,周五晚上平均每个小时的顾客数是20。
lmd = 20
poisson = scipy.stats.poisson(lmd)
x = np.arange(0,40)
plt.plot(x, poisson.pmf(x), 'bo')
plt.vlines(x, 0, poisson.pmf(x), colors='b')
plt.ylim(0,0.1)
plt.show()
mean, var = poisson.stats()
print(mean)
print(var)
20.0
20.0
#顾客数恰好是20的概率?
poisson.pmf(20)
0.088835317392084806
# 顾客数小于15的函数
poisson.cdf(15)
0.15651313463974229
# 顾客数大于等于20的概率?
poisson.sf(19)
0.52974273316075782
# 百分之90%的情况下,顾客数不会超过多少?
poisson.ppf(0.9)
26.0
#作业## 基本作业
#机票超卖现象
#假设某国际航班有300个座位,乘客平均误机率是2%。
#1、如果一共卖出305张机票,那么登机时人数超额的概率是多少?
# 某家风投企业,投资成功的概率是5%,如果它投了100个项目,恰好有5个成功的概率是多少?
n = 305
p = 0.98
binom = scipy.stats.binom(n,p)
binom.pmf(301)
0.12929463926438506
1 - binom.cdf(300)
0.26915013819815137
binom.pmf(301)
0.12929463926438506
n = 305
p = 0.98
binom = scipy.stats.binom(n,p)
lmd=300
poisson = scipy.stats.poisson(lmd)
poisson.ppf(0.9)
322.0