import scipy.stats
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
%config InlindomneBackend.figure_format ='retina'
計(jì)算機(jī)模擬
生成隨機(jī)數(shù)
np.random.seed(5)
疑惑: 為什么要有種子呢?
np.random.random(5)
array([ 0.22199317, 0.87073231, 0.20671916, 0.91861091, 0.48841119])
np.random.randint(0,9,10)
array([0, 4, 4, 3, 2, 4, 6, 3, 3, 2])
num =10000
x = np.random.random(num)
y = np.random.random(num)
pi = np.sum(x**2+y**2<1)/num*4
print ("π=",pi)
π= 3.1756
- 小疑惑,random()生成的隨機(jī)數(shù),是默認(rèn)為0-1之間的嗎?
-- 剛才試著將結(jié)果打印出來(lái)了,確實(shí)是0-1之間的數(shù)字,沒有大于1的隨機(jī)數(shù)
我將num設(shè)置到1億,結(jié)果電腦直接卡機(jī)了,強(qiáng)制關(guān)機(jī)后才緩過(guò)來(lái)。
- 疑問(wèn)。x2+y2< 1,這個(gè)結(jié)果是數(shù)量嗎?總覺得應(yīng)該用counts()方法,而非sum()方法啊
--原來(lái)x2+y2< 1結(jié)果是一堆布爾型的結(jié)果,而sum是不是只是把為true的結(jié)果進(jìn)行相加了,ture又是1的意思。
x**2+y**2<1
array([ True, True, True, ..., True, True, True], dtype=bool)
num*4
4000000
接下來(lái),還可以畫圖,我們來(lái)試試吧
plt.figure(figsize=(9,9))
plt.scatter(x,y,alpha=0.45)
x2 = np.arange(0,1.01,0.01)
y2 = np.sqrt(1-x2**2)
plt.plot(x2,y2,'r',lw=2)
#plt.plot(x2, y2, 'm', lw=3)
plt.show()
新增知識(shí)點(diǎn):
- 圓半徑為1的情況下,如何計(jì)算另外兩條直角邊。 1-x2=y2。y=sqrt(1-x2)
- sum()一個(gè)過(guò)濾函數(shù),是將布爾型結(jié)果中的true,計(jì)數(shù)求和。如:np.sum(x2+y2<1)
錯(cuò)誤記錄:
- . 拼寫錯(cuò)誤,arange,我多寫了一個(gè)r,arrange。

np.arange(0, 1.01, 0.01)
array([ 0. , 0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.07, 0.08,
0.09, 0.1 , 0.11, 0.12, 0.13, 0.14, 0.15, 0.16, 0.17,
0.18, 0.19, 0.2 , 0.21, 0.22, 0.23, 0.24, 0.25, 0.26,
0.27, 0.28, 0.29, 0.3 , 0.31, 0.32, 0.33, 0.34, 0.35,
0.36, 0.37, 0.38, 0.39, 0.4 , 0.41, 0.42, 0.43, 0.44,
0.45, 0.46, 0.47, 0.48, 0.49, 0.5 , 0.51, 0.52, 0.53,
0.54, 0.55, 0.56, 0.57, 0.58, 0.59, 0.6 , 0.61, 0.62,
0.63, 0.64, 0.65, 0.66, 0.67, 0.68, 0.69, 0.7 , 0.71,
0.72, 0.73, 0.74, 0.75, 0.76, 0.77, 0.78, 0.79, 0.8 ,
0.81, 0.82, 0.83, 0.84, 0.85, 0.86, 0.87, 0.88, 0.89,
0.9 , 0.91, 0.92, 0.93, 0.94, 0.95, 0.96, 0.97, 0.98,
0.99, 1. ])
個(gè)人理解:測(cè)試一下arange方法的作用。np.arange(x,y,d)是按照順序生成x-y區(qū)間的數(shù)字,單位是d,也稱作間隔為d。
連續(xù)分布--正態(tài)分布
模擬面包重量的分布
mean = 950
std = 50
# 生成滿足正態(tài)分布的隨機(jī)數(shù),并繪制直方圖
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
sample = np.random.normal(mean,std,size=365)
plt.hist(sample,bins = 30,rwidth=0.9,normed=True)
plt.show()
錯(cuò)誤記錄:
- (mean,std,size=365)中的逗號(hào),我寫作了.。一直報(bào)錯(cuò),都沒發(fā)現(xiàn)原因,有點(diǎn)傻啊

用scipy.stats也可以生成正太分布哦
mean = 950
std = 50
bom =scipy.stats.norm(mean,std)
錯(cuò)誤記錄
- scipy中正態(tài)分布方法,只是簡(jiǎn)寫成norm。 scipy.stats.norm 而我寫成額normal
x = np.arange(700,1200,1)
y = bom.pdf(x)
plt.plot(x,y)
plt.show()

我在這里嘗試在notebook中安裝seaborn,然后成功了。哈哈哈
!pip install seaborn
Requirement already satisfied: seaborn in d:\programdata\anaconda3\lib\site-packages
import seaborn as sns
sns.boxplot(x,y)
sns.plt.show()

- 自己強(qiáng)行胡亂嘗試用seaborn畫圖,結(jié)果搗鼓出來(lái)一個(gè)。
# 接下來(lái)畫概率累計(jì)函數(shù)
x = np.arange(700,1200,1)
#y = np.norm.cdf(x)
y = bom.cdf(x)
plt.plot(x,y)
plt.show()
錯(cuò)誤記錄:
- norm.cdf被我寫作為np.norm.cdf。 其實(shí)這個(gè)norm是個(gè)自定義變量啊

計(jì)算買到的面包小于1000克的概率
這個(gè)該怎么計(jì)算呢?有套路嗎?
先畫出概率密度函數(shù)的曲線
x = np.arange(700,1200,1)
y = bom.pdf(x)
plt.plot(x,y)
plt.vlines =(1000,0,bom.pdf(1000))
> ???我的這段代碼里,為什么plt.vlines()方法沒有反應(yīng)呢?
x2 = np.arange(700,1000,1)
y2 = bom.pdf(x2)
plt.plot(x2,y2)
plt.fill_between(x2,y2,color= 'red',alpha=0.2)
plt.show()
錯(cuò)誤記錄:
1.plt.fill_between()我寫的是fillbetween,竟然沒有加下劃線。
- plt.fill_between()內(nèi)部的參數(shù),不僅僅應(yīng)該是x軸,應(yīng)該xy軸,而我的參數(shù)設(shè)定只有700,1000,這是只基于x軸的思路。
- plt中的顏色設(shè)定都是color=“”的形式,而是采用的簡(jiǎn)寫實(shí)在不科學(xué)。
- 概念混淆。 pdf和cdf完全混淆了。 所以我的顯示完全是兩張圖。 我x,y的函數(shù)用的是pdf,而x2,y2,用的是cdf函數(shù),所以總共覺得圖形不對(duì)。
疑惑,為甚么我的圖形,沒有完整顯示呢?-- 因?yàn)閜lt.vilines沒起作用

紅色區(qū)域的面積,就是700-1000g的概率和。 那么概率和,就應(yīng)該用概率累計(jì)函數(shù)。
print("面包等于1000g的概率是:",bom.pdf(1000))
print("小于1000g的概率是:",bom.cdf(1000))
面包等于1000g的概率是: 0.00483941449038
小于1000g的概率是: 0.841344746069
回到剛才的自我提問(wèn),小于1000g的概率,有套路嗎?
有!
小于1000g的概率,就是求700g的概率+701g的概率+.......+999g的概率,也就是概率累計(jì)和嘛。直接求cdf(1000)即可。
求面包大于1000g的概率
這下簡(jiǎn)單了,就是求1000g的概率加到1200g的概率嘛,用整體減去700到1000g的概率和,就可以了呀。
整體等于1,那就好辦啦。
print('面包大于1000g的概率是:',1-bom.cdf(1000))
面包大于1000g的概率是: 0.158655253931
print('面包大于1000g的概率是:',bom.sf(1000))
##錯(cuò)誤記錄:計(jì)劃求反函數(shù),誰(shuí)知記錯(cuò)了方法名,應(yīng)該用sf,而非isf
面包大于1000g的概率是: 0.158655253931
剛才搜到了幾個(gè)方法的解釋,常用的概率函數(shù)主要可以概括為以下幾類:
- 根據(jù)變量求小于變量的概率(cdf)
- 根據(jù)變量求大于變量的概率(sf)
- 根據(jù)概率求相應(yīng)的小于變量(ppf)
- 根據(jù)概率求相應(yīng)的大于變量(isf)
接著做題
計(jì)算買到的面包處于950到1050范圍內(nèi)的概率
## 這種情況下,直接用cdf(1050)-cdf(950)就可以啊
#繪制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)
#設(shè)置y軸范圍
plt.ylim(0,0.0085)
plt.show()
print('面包950-1000g之間的概率是:',bom.cdf(1050)-bom.cdf(950))

面包950-1000g之間的概率是: 0.477249868052
90%的情況下,買到的面包是小于多少克的?
- 根據(jù)概率求相應(yīng)的小于變量(ppf)
這就是典型的用概率求變量。已知概率,求小于變量,該怎么算呢。 小于1000g的概率,是用概率累計(jì)函數(shù)cdf,即cdf(1000)
已知概率,求變量則應(yīng)該是cdf的反函數(shù)ppf.
先看看概率累計(jì)函數(shù)吧
x = np.arange(700, 1200, 1)
y = norm.cdf(x)
plt.plot(x, y)
plt.show()

print('90%的時(shí)候面包小于:',bom.ppf(0.9))
90%的時(shí)候面包小于: 1014.07757828
80%的情況下,買到的面包是大于多少克的?
- 根據(jù)概率求相應(yīng)的大于變量(isf)
這種時(shí)候,就該用sf(survival function)的反函數(shù)isf(inverse survival function)。
print ('80%的時(shí)候面包大于:',bom.isf(0.8),'g')
80%的時(shí)候面包大于: 907.918938321 g
離散分布 - 二項(xiàng)分布
投硬幣問(wèn)題模擬
outcome =np.random.randint(0,2,10)
#outcome= list(outcome)
np.sum(outcome)
- 哦哦,明白啦,布爾型數(shù)據(jù)是0和1,然后用sum求和方式來(lái)計(jì)算,能夠參與計(jì)算的只有1啊,就變相的數(shù)了有多少個(gè)1,好高明啊
3
sample=[np.sum(np.random.randint(0,2,10))for i in range(1000)]
# 隨機(jī)生成10個(gè),0-2之間的整數(shù)(其實(shí)就是0,1),然后重復(fù)這個(gè)動(dòng)作1000次,則會(huì)生成1000組10個(gè)數(shù)組
sample= pd.Series(sample)
sample
0 7
1 4
2 4
3 6
4 6
5 5
6 6
7 4
8 5
9 4
10 4
11 4
12 7
13 5
14 4
15 6
16 1
17 6
18 4
19 8
20 6
21 3
22 6
23 7
24 7
25 7
26 2
27 4
28 2
29 4
..
970 5
971 4
972 6
973 4
974 5
975 5
976 4
977 4
978 7
979 1
980 3
981 4
982 4
983 7
984 6
985 3
986 5
987 4
988 5
989 4
990 6
991 3
992 4
993 5
994 3
995 6
996 5
997 9
998 3
999 4
Length: 1000, dtype: int64
#sample.value_counts().plot.bar()
#sample.value_counts().sort_index().plot.bar()
#plt.show()
##錯(cuò)誤記錄,sample還需要數(shù)組化(相當(dāng)于列表化),但我沒數(shù)組化后賦值給sample,使得已知報(bào)錯(cuò),無(wú)法運(yùn)算
sample.value_counts().sort_index().plot.bar()
plt.show()

投硬幣問(wèn)題的二項(xiàng)分布
n = 10
p = 0.4
mean_bionom = n*p
std_bionom =n*p*(1-p)
print(n*p,n*p*(1-p))
binomial = scipy.stats.binom(n,p)
print (binomial)
4.0 2.4
<scipy.stats._distn_infrastructure.rv_frozen object at 0x0000000009FCFE80>
x = np.arange(0,11,1)
y = binomial.pmf(x)
plt.plot(x,y)
plt.vlines(x,0,y,colors='b')
##疑惑:為什么這里要寫0
plt.show()
錯(cuò)誤記錄: 1.寫y值時(shí),我總是要用np.binom.pmf(x),要加上np前綴,但實(shí)際上上,data.pmf(x)就可以啦。
---------------------------------------------------------------------------
TypeError Traceback (most recent call last)
<ipython-input-32-4c130d5dd029> in <module>()
3 plt.plot(x,y)
4 #plt.vlines(x,0,y,colors='b')
----> 5 plt.vlines(x, 0, binomial.pmf(x), colors='b')
6 ##疑惑:為什么這里要寫0
7 plt.show()
TypeError: 'tuple' object is not callable
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()
---------------------------------------------------------------------------
TypeError Traceback (most recent call last)
<ipython-input-35-466792604d85> in <module>()
1 x = np.arange(0,11)
2 plt.plot(x, binomial.pmf(x), 'bo')
----> 3 plt.vlines(x, 0, binomial.pmf(x), colors='b')
4 plt.ylim(0,0.3)
5 plt.show()
TypeError: 'tuple' object is not callable
mean,var =binomial.stats()
print(binomial.stats())
(array(4.0), array(2.4))
應(yīng)用
某家風(fēng)投企業(yè),投資成功的概率是5%,如果它投了100個(gè)項(xiàng)目,恰好有5個(gè)成功的概率是多少?
n = 100
p = 0.05
binom = scipy.stats.binom(n,p)
print(binom)
<scipy.stats._distn_infrastructure.rv_frozen object at 0x000000000B69F668>
##恰好成功5個(gè),就是求5的對(duì)應(yīng)概率,就是pmf
binom.pmf(5)
0.18001782727043672
投資至少成功了5個(gè)的概率?
#至少成功5個(gè),就是大于等于5
1-binom.cdf(4)
0.56401869931428927
binom.sf(4)
0.56401869931429105
10%的情況下,至少能成功幾個(gè)?
#這也是一個(gè)大于等于的情況,用生存函數(shù)
binom.isf(0.1)
8.0
離散分布 - 泊松分布
有一家便利店使用泊松分布來(lái)估計(jì)周五晚上到店買東西的顧客數(shù),根據(jù)以往數(shù)據(jù),周五晚上平均每個(gè)小時(shí)的顧客數(shù)是20。
lmd = 20
poisson = scipy.stats.poisson(lmd)
#錯(cuò)誤記錄: poisson 被我拼寫為poission。 讀了個(gè)i
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=',mean,'var=',var)
mean= 20.0 var= 20.0
#顧客數(shù)恰好為20的概率
poisson.pmf(20)
0.088835317392084806
#顧客數(shù)小于等于15
poisson.cdf(15)
0.1565131346397429
#顧客數(shù)大于等于20
poisson.sf(19)
錯(cuò)誤記錄:
- 大于等于20,就是包含20,那么久該從19之后算起,如果我輸入20,則是從20之后開始計(jì)算。 應(yīng)該是19
疑惑
- 不過(guò)什么時(shí)候該用20,什么時(shí)候該用19呢,難道是大于的時(shí)候嗎? 為什么小于等于15,就使用cdf(15),而大于等于20就使用sf(19)呢?
0.44090741576867482
poisson.ppf(0.9)
26.0
基本作業(yè)
機(jī)票超賣現(xiàn)象
假設(shè)某國(guó)際航班有300個(gè)座位,乘客平均誤機(jī)率是2%。
1、如果一共賣出305張機(jī)票,那么登機(jī)時(shí)人數(shù)超額的概率是多少?
n = 305
p=0.98
binomial = scipy.stats.binom(n,p)
plt.plot(np.arange(290,310,1),binomial.pmf(np.arange(290,310,1)))
plt.show()

## 超額的概率,就是實(shí)際人數(shù)超過(guò)300的概率. 此時(shí)n=305
#應(yīng)該計(jì)算大于等于301的情況,sf函數(shù)
print ('超員的概率為:',binomial.sf(300))
超員的概率為: 0.269150138198
2、如果一共賣出305張機(jī)票,登機(jī)時(shí)最多只超額1人的概率是多少?
#只超額1人,就是小于等于301嘛。 應(yīng)該使用cdf
print ('只超額1人的概率:',binomial.cdf(301))
只超額1人的概率: 0.860144501066
3、一共賣幾張票,可以保證不超額的概率至少是90%。
n1 = 309
bm1=scipy.stats.binom(n1,p)
x*0.9<300
array([ True, True, True, True, True, True, True, True, True,
True, True, True, True, True, True, True, True, True,
True, True, True, True, True, True, True, True, True,
True, True, True, True, True, True, True, True, True,
True, True, True, True], dtype=bool)
bm1.isf(0.9)
300.0