精通 Pandas:6~11
原文:Mastering Pandas
协议:CC BY-NC-SA 4.0
译者:飞龙
六、处理缺失数据,时间序列和 Matplotlib 绘图
在本章中,我们将介绍一些必要的主题,这些主题对于培养使用 Pandas 的专业知识必不可少。 这些主题的知识对于准备数据作为处理数据以进行分析,预测或可视化的程序或代码的输入非常有用。 我们将讨论的主题如下:
- 处理缺失的数据
- 处理时间序列和日期
- 使用
matplotlib
绘图
到本章结束时,用户应该精通这些关键领域。
处理缺失的数据
数据丢失是指由于某种原因在我们的数据集中显示为 NULL 或 N/A 的数据点; 例如,我们可能有一个时间序列,横跨一个月的所有日历日,显示每天股票的收盘价,而非营业日的收盘价则显示为缺失。 一个损坏的数据示例是财务数据集,该数据集以错误的格式显示了交易的活动日期。 例如,由于数据提供者发生错误,因此使用YYYY-MM-DD
而不是YYYYMMDD
。
对于 Pandas,缺失值通常由 NaN 值表示。
除了本机出现在源数据集中之外,还可以通过诸如重新索引或在时间序列的情况下更改频率之类的操作将缺失值添加到数据集中:
In [84]: import numpy as npimport pandas as pdimport matplotlib.pyplot as plt%matplotlib inline
In [85]: date_stngs = ['2014-05-01','2014-05-02','2014-05-05','2014-05-06','2014-05-07']tradeDates = pd.to_datetime(pd.Series(date_stngs))
In [86]: closingPrices=[531.35,527.93,527.81,515.14,509.96]
In [87]: googClosingPrices=pd.DataFrame(data=closingPrices,columns=['closingPrice'],index=tradeDates)googClosingPrices
Out[87]: closingPricetradeDates 2014-05-01 531.352014-05-02 527.932014-05-05 527.812014-05-06 515.142014-05-07 509.965 rows 1 columns
可以在这里找到上述数据的来源。
Pandas 还提供了一个 API,可以从各种数据提供商(例如 Yahoo)读取股票数据:
In [29]: import pandas.io.data as web
In [32]: import datetimegoogPrices = web.get_data_yahoo("GOOG",start=datetime.datetime(2014, 5, 1),end=datetime.datetime(2014, 5, 7))
In [38]: googFinalPrices=pd.DataFrame(googPrices['Close'],index=tradeDates)
In [39]: googFinalPrices
Out[39]: Close2014-05-01 531.349982014-05-02 527.929992014-05-05 527.810002014-05-06 515.140012014-05-07 509.95999
有关更多详细信息,请参见这里。
现在,我们有了一个时间序列,描述了 Google 股票自 2014 年 5 月 1 日至 2014 年 5 月 7 日的收盘价,由于交易只在工作日发生,因此日期范围存在缺口。 如果要更改日期范围以使其显示日历日(即周末),则可以将时间序列索引的频率从工作日更改为日历日,如下所示:
In [90]: googClosingPricesCDays=googClosingPrices.asfreq('D')googClosingPricesCDays
Out[90]: closingPrice2014-05-01 531.352014-05-02 527.932014-05-03 NaN2014-05-04 NaN2014-05-05 527.812014-05-06 515.142014-05-07 509.967 rows 1 columns
请注意,我们现在为 2014 年 5 月 3 日和 2014 年 5 月 4 日的周末日期引入了closingPrice
的NaN
值。
我们可以通过使用isnull
和notnull
函数来检查缺少的值,如下所示:
In [17]: googClosingPricesCDays.isnull()
Out[17]: closingPrice2014-05-01 False2014-05-02 False2014-05-03 True2014-05-04 True2014-05-05 False2014-05-06 False2014-05-07 False7 rows 1 columnsIn [18]: googClosingPricesCDays.notnull()
Out[18]: closingPrice2014-05-01 True2014-05-02 True2014-05-03 False2014-05-04 False2014-05-05 True2014-05-06 True2014-05-07 True7 rows 1 columns
在每种情况下都会返回一个布尔型数据帧。 在datetime
和 Pandas 时间戳中,缺失值由NaT
值表示。 对于基于时间的类型,这相当于 Pandas 中的NaN
。
In [27]: tDates=tradeDates.copy()tDates[1]=np.NaNtDates[4]=np.NaNIn [28]: tDates
Out[28]: 0 2014-05-011 NaT2 2014-05-053 2014-05-064 NaTName: tradeDates, dtype: datetime64[ns]In [4]: FBVolume=[82.34,54.11,45.99,55.86,78.5]TWTRVolume=[15.74,12.71,10.39,134.62,68.84]In [5]: socialTradingVolume=pd.concat([pd.Series(FBVolume),pd.Series(TWTRVolume),tradeDates], axis=1,keys=['FB','TWTR','TradeDate'])socialTradingVolume
Out[5]: FB TWTR TradeDate0 82.34 15.74 2014-05-011 54.11 12.71 2014-05-022 45.99 10.39 2014-05-053 55.86 134.62 2014-05-064 78.50 68.84 2014-05-075 rows × 3 columnsIn [6]: socialTradingVolTS=socialTradingVolume.set_index('TradeDate')
socialTradingVolTS
Out[6]:TradeDate FB TWTR2014-05-01 82.34 15.742014-05-02 54.11 12.712014-05-05 45.99 10.392014-05-06 55.86 134.622014-05-07 78.50 68.845 rows × 2 columnsIn [7]: socialTradingVolTSCal=socialTradingVolTS.asfreq('D')socialTradingVolTSCal
Out[7]:FB TWTR2014-05-01 82.34 15.742014-05-02 54.11 12.712014-05-03 NaN NaN2014-05-04 NaN NaN2014-05-05 45.99 10.392014-05-06 55.86 134.622014-05-07 78.50 68.847 rows × 2 columns
我们可以对包含缺失值的数据执行算术运算。 例如,我们可以计算 Facebook 和 Twitter 的两只股票的总交易量(百万股),如下所示:
In [8]: socialTradingVolTSCal['FB']+socialTradingVolTSCal['TWTR']
Out[8]: 2014-05-01 98.082014-05-02 66.822014-05-03 NaN2014-05-04 NaN2014-05-05 56.382014-05-06 190.482014-05-07 147.34Freq: D, dtype: float64
默认情况下,对包含缺失值的对象执行的任何操作都会在该位置返回缺失值,如以下命令所示:
In [12]: pd.Series([1.0,np.NaN,5.9,6])+pd.Series([3,5,2,5.6])
Out[12]: 0 4.01 NaN2 7.93 11.6dtype: float64
In [13]: pd.Series([1.0,25.0,5.5,6])/pd.Series([3,np.NaN,2,5.6])
Out[13]: 0 0.3333331 NaN2 2.7500003 1.071429dtype: float64
但是,NumPy 处理聚合计算的方式与 Pandas 的处理方式有所不同。
在 Pandas 中,默认值是将缺失值视为0
并进行聚合计算,而对于 NumPy,如果缺少任何值,则返回NaN
。 这是一个例子:
In [15]: np.mean([1.0,np.NaN,5.9,6])
Out[15]: nanIn [16]: np.sum([1.0,np.NaN,5.9,6])
Out[16]: nan
但是,如果此数据在 pandas 序列中,则将获得以下输出:
In [17]: pd.Series([1.0,np.NaN,5.9,6]).sum()
Out[17]: 12.9
In [18]: pd.Series([1.0,np.NaN,5.9,6]).mean()
Out[18]: 4.2999999999999998
重要的是要意识到 Pandas 和 NumPy 在行为上的差异。 但是,如果我们希望 NumPy 的行为与 Pandas 相同,则可以使用np.nanmean
和np.nansum
函数,如下所示:
In [41]: np.nanmean([1.0,np.NaN,5.9,6])
Out[41]: 4.2999999999999998In [43]: np.nansum([1.0,np.NaN,5.9,6])
Out[43]: 12.9
有关 NumPy np.nan
聚合函数的更多信息,请参考这里。
处理缺失值
有多种处理缺失值的方法,如下所示:
-
通过使用
fillna()
函数来填充 NA 值。 这是一个例子:In [19]: socialTradingVolTSCal Out[19]: FB TWTR2014-05-01 82.34 15.742014-05-02 54.11 12.712014-05-03 NaN NaN2014-05-04 NaN NaN2014-05-05 45.99 10.392014-05-06 55.86 134.622014-05-07 78.50 68.847 rows × 2 columnsIn [20]: socialTradingVolTSCal.fillna(100) Out[20]: FB TWTR2014-05-01 82.34 15.742014-05-02 54.11 12.712014-05-03 100.00 100.002014-05-04 100.00 100.002014-05-05 45.99 10.392014-05-06 55.86 134.622014-05-07 78.50 68.847 rows × 2 columns
我们还可以使用
ffill
或bfill
参数填充前向或后向值:In [23]: socialTradingVolTSCal.fillna(method='ffill') Out[23]: FB TWTR2014-05-01 82.34 15.742014-05-02 54.11 12.712014-05-03 54.11 12.712014-05-04 54.11 12.712014-05-05 45.99 10.392014-05-06 55.86 134.622014-05-07 78.50 68.847 rows × 2 columnsIn [24]: socialTradingVolTSCal.fillna(method='bfill') Out[24]: FB TWTR2014-05-01 82.34 15.742014-05-02 54.11 12.712014-05-03 45.99 10.392014-05-04 45.99 10.392014-05-05 45.99 10.392014-05-06 55.86 134.622014-05-07 78.50 68.847 rows × 2 columns
pad
方法是ffill
的替代名称。 有关更多详细信息,您可以转到这里。 -
通过使用
dropna()
函数删除/删除缺少值的行和列。 以下是一个示例:In [21]: socialTradingVolTSCal.dropna() Out[21]: FB TWTR2014-05-01 82.34 15.742014-05-02 54.11 12.712014-05-05 45.99 10.392014-05-06 55.86 134.622014-05-07 78.50 68.845 rows × 2 columns
-
我们还可以使用
interpolate()
函数对缺失值进行插值和填充,如以下命令所述:In [27]: pd.set_option('display.precision',4)socialTradingVolTSCal.interpolate() Out[27]: FB TWTR2014-05-01 82.340 15.7402014-05-02 54.110 12.7102014-05-03 51.403 11.9372014-05-04 48.697 11.1632014-05-05 45.990 10.3902014-05-06 55.860 134.6202014-05-07 78.500 68.8407 rows × 2 columns
interpolate()
函数还采用一个参数 – 表示该方法的method
。 这些方法包括线性,二次,三次样条等等。 您可以从官方文档中获取更多信息。
处理时间序列
在本节中,我们向您展示如何处理时间序列数据。 我们将首先展示如何使用从csv
文件中读取的数据创建时间序列数据。
读取时间序列数据
在这里,我们演示了读取时间序列数据的各种方法:
In [7]: ibmData=pd.read_csv('ibm-common-stock-closing-prices-1959_1960.csv')ibmData.head()
Out[7]: TradeDate closingPrice0 1959-06-29 4451 1959-06-30 4482 1959-07-01 4503 1959-07-02 4474 1959-07-06 4515 rows 2 columns
可以在这个链接中找到此信息的来源。
我们希望TradeDate
列是datetime
值的序列,以便我们可以为其编制索引并创建时间序列。 让我们首先检查TradeDate
序列中值的类型:
In [16]: type(ibmData['TradeDate'])
Out[16]: pandas.core.series.Series
In [12]: type(ibmData['TradeDate'][0])
Out[12]: str
接下来,我们将其转换为Timestamp
类型:
In [17]: ibmData['TradeDate']=pd.to_datetime(ibmData['TradeDate'])type(ibmData['TradeDate'][0])
Out[17]: pandas.tslib.Timestamp
现在,我们可以将TradeDate
列用作索引:
In [113]: #Convert DataFrame to TimeSeries#Resampling creates NaN rows for weekend dates, hence use dropnaibmTS=ibmData.set_index('TradeDate').resample('D')['closingPrice'].dropna()ibmTS
Out[113]: TradeDate1959-06-29 4451959-06-30 4481959-07-01 4501959-07-02 4471959-07-06 451...Name: closingPrice, Length: 255
日期偏移和时间增量对象
DateOffset
对象表示时间的变化或偏移。 DateOffset
对象的关键特征如下:
- 可以将其添加到
datetime
对象或从中减去,以获得转换后的日期 - 可以乘以一个整数(正数或负数),以便可以多次应用该增量
- 它具有前滚和后滚方法,可将日期向前移动到下一个偏移日期或向后移动到上一个偏移日期
我们说明了如何使用日期偏移对象,如下所示:
In [371]: xmasDay=pd.datetime(2014,12,25)xmasDay
Out[371]: datetime.datetime(2014, 12, 25, 0, 0)In [373]: boxingDay=xmasDay+pd.DateOffset(days=1)boxingDay
Out[373]: Timestamp('2014-12-26 00:00:00', tz=None)In [390}: today=pd.datetime.now()today
Out[390]: datetime.datetime(2014, 5, 31, 13, 7, 36, 440060)
注意,datetime.datetime
与pd.Timestamp
不同。 前者是 Python 类,效率低下,而后者基于numpy.datetime64
数据类型。 pd.DateOffset
对象与pd.Timestamp
一起使用,并将其添加到datetime.datetime
函数中可将该对象转换为pd.Timestamp
对象。
下面说明了从今天开始一周的命令:
In [392]: today+pd.DateOffset(weeks=1)
Out[392]: Timestamp('2014-06-07 13:07:36.440060', tz=None)
下图说明了从现在起五年内的命令:
In [394]: today+2*pd.DateOffset(years=2, months=6)
Out[394]: Timestamp('2019-05-30 13:07:36.440060', tz=None)
这是使用rollforward
函数的示例。 QuarterBegin
是一个DateOffset
对象,用于将给定的datetime
对象增加到下一个日历季度的开始:
In [18]: lastDay=pd.datetime(2013,12,31)
In [24]: from pandas.tseries.offsets import QuarterBegindtoffset=QuarterBegin()lastDay+dtoffset
Out[24]: Timestamp('2014-03-01 00:00:00', tz=None)In [25]: dtoffset.rollforward(lastDay)
Out[25]: Timestamp('2014-03-01 00:00:00', tz=None)
因此,我们可以看到,2013 年 12 月 31 日之后的下一个季度从 2014 年 3 月 1 日开始。Timedeltas
与DateOffsets
相似,但可用于datetime.datetime
对象。 以下命令解释了它们的使用:
In [40]: weekDelta=datetime.timedelta(weeks=1)weekDelta
Out[40]: datetime.timedelta(7)In [39]: today=pd.datetime.now()today
Out[39]: datetime.datetime (2014, 6, 2, 3, 56, 0, 600309)In [41]: today+weekDelta
Out[41]: datetime.datetime (2014, 6, 9, 3, 56,0, 600309)
与时间序列相关的实例方法
在本节中,我们探索用于时间序列对象的各种方法,例如移位,频率转换和重采样。
平移/滞后
有时,我们可能希望将时间序列中的值在时间上向后或向前移动。 一种可能的情况是,数据集包含公司中去年新雇员的开始日期列表,并且公司的人力资源计划希望将这些日期提前一年,以便可以激活雇员的福利。 我们可以通过使用shift()
函数来做到这一点,如下所示:
In [117]: ibmTS.shift(3)
Out[117]: TradeDate1959-06-29 NaN1959-06-30 NaN1959-07-01 NaN1959-07-02 4451959-07-06 4481959-07-07 4501959-07-08 447...
这将转换所有日历日。 但是,如果我们只希望转移工作日,则必须使用以下命令:
In [119]: ibmTS.shift(3, freq=pd.datetools.bday)
Out[119]: TradeDate1959-07-02 4451959-07-03 4481959-07-06 4501959-07-07 4471959-07-09 451
在前面的代码片段中,我们指定了freq
参数进行平移; 这告诉函数仅更改工作日。 shift
函数具有freq
参数,其值可以是DateOffset
类,类似于timedelta
的对象或偏移别名。 因此,使用ibmTS.shift(3, freq='B')
也将产生相同的结果。
更改频率
我们可以使用asfreq
函数来更改频率,如下所述:
In [131]: # Frequency conversion using asfreqibmTS.asfreq('BM')
Out[131]: 1959-06-30 4481959-07-31 4281959-08-31 4251959-09-30 4111959-10-30 4111959-11-30 4281959-12-31 4391960-01-29 4181960-02-29 4191960-03-31 4451960-04-29 4531960-05-31 5041960-06-30 522Freq: BM, Name: closingPrice, dtype: float64
在这种情况下,我们仅从ibmTS
时间序列中获取与该月的最后一天相对应的值。 在此,bm
代表营业月结束频率。 有关所有可能的频率别名的列表,请访问这里。
如果我们指定的频率小于数据的粒度,则间隙将用NaN
值填充:
In [132]: ibmTS.asfreq('H')
Out[132]: 1959-06-29 00:00:00 4451959-06-29 01:00:00 NaN1959-06-29 02:00:00 NaN1959-06-29 03:00:00 NaN...1960-06-29 23:00:00 NaN1960-06-30 00:00:00 522Freq: H, Name: closingPrice, Length: 8809
我们也可以将asfreq
方法应用于Period
和PeriodIndex
对象,类似于我们对datetime
和Timestamp
对象所做的操作。 Period
和PeriodIndex
稍后介绍,用于表示时间间隔。
asfreq
方法接受一个方法参数,该参数允许您向前填充(ffill
)或向后填充空白,类似于fillna
:
In [140]: ibmTS.asfreq('H', method='ffill')
Out[140]: 1959-06-29 00:00:00 4451959-06-29 01:00:00 4451959-06-29 02:00:00 4451959-06-29 03:00:00 445...1960-06-29 23:00:00 5221960-06-30 00:00:00 522Freq: H, Name: closingPrice, Length: 8809
数据重采样
TimeSeries.resample
函数使我们能够基于采样间隔和采样函数来聚合/聚合更多粒度数据。
下采样是源自数字信号处理的术语,是指降低信号的采样率的过程。 对于数据,我们使用它来减少我们希望处理的数据量。
相反的过程是上采样,该过程用于增加要处理的数据量,并且需要进行插值以获得中间数据点。 有关下采样和上采样的更多信息,请参考上采样和下采样的实际应用和用于视觉表示的下采样时间序列。
在这里,我们检查了一些滴答数据以用于重采样。 在检查数据之前,我们需要进行准备。 通过这样做,我们将学习一些有关时间序列数据的有用技术,如下所示:
- 时间戳
- 时区处理
这是一个使用滴答数据作为 2014 年 5 月 27 日星期二的 Google 股票价格的示例:
In [150]: googTickData=pd.read_csv('./GOOG_tickdata_20140527.csv')
In [151]: googTickData.head()
Out[151]: Timestamp close high low open volume0 1401197402 555.008 556.41 554.35 556.38 811001 1401197460 556.250 556.30 555.25 555.25 185002 1401197526 556.730 556.75 556.05 556.39 99003 1401197582 557.480 557.67 556.73 556.73 147004 1401197642 558.155 558.66 557.48 557.59 157005 rows 6 columns
可以在这个链接中找到先前数据的源。
从上一节中可以看到,我们有一个“时间戳”列,以及收盘价,最高价,最低价和开盘价以及 Google 股票交易量的列。
那么,为什么“时间戳”列看起来有点奇怪? 好吧,滴答数据时间戳通常以纪元时间表示(有关更多信息,请参考这里),作为一种更紧凑的存储方式。 我们需要将其转换为更易于理解的时间,我们可以按照以下步骤进行操作:
In [201]: googTickData['tstamp']=pd.to_datetime(googTickData['Timestamp'],unit='s',utc=True)In [209]: googTickData.head()
Out[209]:Timestamp close high low open volume tstamp0 14011974020 555.008 556.41 554.35 556.38 81100 2014-05-27 13:30:021 1401197460 556.250 556.30 555.25 555.25 18500 2014-05-27 13:31:002 1401197526 556.730 556.75 556.05 556.39 9900 2014-05-27 13:32:063 1401197582 557.480 557.67 556.73 556.73 14700 2014-05-27 13:33:024 1401197642 558.155 558.66 557.48 557.59 15700 2014-05-27 13:34:025 rows 7 columns
现在,我们想将tstamp
列作为索引,并消除纪元Timestamp
列:
In [210]: googTickTS=googTickData.set_index('tstamp')googTickTS=googTickTS.drop('Timestamp',axis=1)googTickTS.head()
Out[210]: tstamp close high low open volume2014-05-27 13:30:02 555.008 556.41 554.35 556.38 8110002014-05-27 13:31:00 556.250 556.30 555.25 555.25 185002014-05-27 13:32:06 556.730 556.75 556.05 556.39 99002014-05-27 13:33:02 557.480 557.67 556.73 556.73 147002014-05-27 13:34:02 558.155 558.66 557.48 557.59 157005 rows 5 columns
请注意,tstamp
索引列的时间以 UTC 为单位,我们可以使用tz_localize
和tz_convert
这两个运算符将其转换为美国/东部时间:
In [211]: googTickTS.index=googTickTS.index.tz_localize('UTC').tz_convert('US/Eastern')In [212]: googTickTS.head()
Out[212]: tstamp close high low open volume 2014-05-27 09:30:02-04:00 555.008 556.41 554.35 556.38 811002014-05-27 09:31:00-04:00 556.250 556.30 555.25 555.25 185002014-05-27 09:32:06-04:00 556.730 556.75 556.05 556.39 99002014-05-27 09:33:02-04:00 557.480 557.67 556.73 556.73 147002014-05-27 09:34:02-04:00 558.155 558.66 557.48 557.59 157005 rows 5 columnsIn [213]: googTickTS.tail()
Out[213]:tstamp close high low open volume2014-05-27 15:56:00-04:00 565.4300 565.48 565.30 565.385 143002014-05-27 15:57:00-04:00 565.3050 565.46 565.20 565.400 147002014-05-27 15:58:00-04:00 565.1101 565.31 565.10 565.310 232002014-05-27 15:59:00-04:00 565.9400 566.00 565.08 565.230 556002014-05-27 16:00:00-04:00 565.9500 565.95 565.95 565.950 1260005 rows 5 columnsIn [214]: len(googTickTS)
Out[214]: 390
从前面的输出中,我们可以看到交易日中每分钟的滴答声-从股市开盘的上午 9:30 到闭市的下午 4:00。 由于在上午 9:30 和下午 4:00 之间有 390 分钟的时间,因此该数据集中有 390 行。
假设我们要每 5 分钟而不是每分钟获取一次快照? 我们可以通过如下使用降采样来实现:
In [216]: googTickTS.resample('5Min').head(6)
Out[216]: close high low open volume tstamp
2014-05-27 09:30:00-04:00 556.72460 557.15800 555.97200 556.46800 27980
2014-05-27 09:35:00-04:00 556.93648 557.64800 556.85100 557.34200 24620
2014-05-27 09:40:00-04:00 556.48600 556.79994 556.27700 556.60678 8620
2014-05-27 09:45:00-04:00 557.05300 557.27600 556.73800 556.96600 9720
2014-05-27 09:50:00-04:00 556.66200 556.93596 556.46400 556.80326 14560
2014-05-27 09:55:00-04:00 555.96580 556.35400 555.85800 556.23600 12400
6 rows 5 columns
用于重采样的默认函数是平均值。 但是,我们还可以指定其他函数,例如最小值,并且可以通过how
参数进行重新采样:
In [245]: googTickTS.resample('10Min', how=np.min).head(4)
Out[245]: close high low open volume
tstamp
2014-05-27 09:30:00-04:00 555.008 556.3000 554.35 555.25 9900
2014-05-27 09:40:00-04:00 556.190 556.5600 556.13 556.35 3500
2014-05-27 09:50:00-04:00 554.770 555.5500 554.77 555.55 3400
2014-05-27 10:00:00-04:00 554.580 554.9847 554.45 554.58 1800
可以将各种函数名称传递给how
参数,例如sum
,ohlc
,max
,min
,std
,mean
,median
,first
和last
。
ohlc
函数根据时间序列数据返回开-高-低-关值; 第一个,最大,最小和最后一个值。 要指定关闭左间隔还是右间隔,我们可以按以下方式传递closed
参数:
In [254]: pd.set_option('display.precision',5)googTickTS.resample('5Min', closed='right').tail(3)
Out[254]: close high low open volume
tstamp
2014-05-27 15:45:00-04:00 564.3167 564.3733 564.1075 564.1700 12816.6667
2014-05-27 15:50:00-04:00 565.1128 565.1725 565.0090 565.0650 13325.0000
2014-05-27 15:55:00-04:00 565.5158 565.6033 565.3083 565.4158 40933.3333
3 rows 5 columns
因此,在前面的命令中,我们可以看到最后一行在 15:55 而不是 16:00 处显示了滴答声。
对于上采样,我们需要指定一种填充方法,以确定如何通过fill_method
参数填充间隙:
In [263]: googTickTS[:3].resample('30s', fill_method='ffill')
Out[263]: close high low open volume tstamp2014-05-27 09:30:00-04:00 555.008 556.41 554.35 556.38 811002014-05-27 09:30:30-04:00 555.008 556.41 554.35 556.38 811002014-05-27 09:31:00-04:00 556.250 556.30 555.25 555.25 185002014-05-27 09:31:30-04:00 556.250 556.30 555.25 555.25 185002014-05-27 09:32:00-04:00 556.730 556.75 556.05 556.39 99005 rows 5 columnsIn [264]: googTickTS[:3].resample('30s', fill_method='bfill')
Out[264]:close high low open volume tstamp2014-05-27 09:30:00-04:00 555.008 556.41 554.35 556.38 811002014-05-27 09:30:30-04:00 556.250 556.30 555.25 555.25 185002014-05-27 09:31:00-04:00 556.250 556.30 555.25 555.25 185002014-05-27 09:31:30-04:00 556.730 556.75 556.05 556.39 99002014-05-27 09:32:00-04:00 556.730 556.75 556.05 556.39 99005 rows 5 columns
不幸的是,fill_method
参数当前仅支持两种方法-前向填充和后向填充。 插值方法将很有价值。
时间序列频率的别名
要指定偏移量,可以使用许多别名。 一些最常用的方法如下:
B, BM
:代表工作日,工作月。 这些是一个月的工作日,即不是假日或周末的任何一天。D, W, M, Q, A
:代表日历日,周,月,季度,年末。H, T, S, L, U
:代表小时,分钟,秒,毫秒和微秒。
这些别名也可以组合。 在以下情况下,我们每 7 分钟 30 秒重新采样一次:
In [267]: googTickTS.resample('7T30S').head(5)
Out[267]:close high low open volume
tstamp
2014-05-27 09:30:00-04:00 556.8266 557.4362 556.3144 556.8800 28075.0
2014-05-27 09:37:30-04:00 556.5889 556.9342 556.4264 556.7206 11642.9
2014-05-27 09:45:00-04:00 556.9921 557.2185 556.7171 556.9871 9800.0
2014-05-27 09:52:30-04:00 556.1824 556.5375 556.0350 556.3896 14350.0
2014-05-27 10:00:00-04:00 555.2111 555.4368 554.8288 554.9675 12512.5
5 rows x 5 columns
可以将后缀应用于频率别名,以指定在频率周期中何时开始。 这些称为锚定偏移量:
W-SUN, MON, ... DEC
,例如,W-TUE
表示从星期二开始的每周频率。Q-JAN, FEB, ... DEC
,例如,Q-MAY
表示每年 5 月底的季度频率。A-JAN, FEB, ... DEC
,例如,A-MAY
表示每年的频率,到 5 月结束。
这些偏移量可用作date_range
和bdate_range
函数的参数,以及用作PeriodIndex
和DatetimeIndex
等索引类型的构造器。 可以在 Pandas 文档 中找到对此的全面讨论。
时间序列的概念和数据类型
处理时间序列时,必须考虑两个主要概念:时间点和范围或时间跨度。 在 Pandas 中,前者由时间戳数据类型表示,该数据类型等效于 Python 的datatime.datetime
(datetime
)数据类型,并且可以互换。 后者(时间跨度)由时间段数据类型表示,该数据类型特定于 Pandas。
这些数据类型均具有与之关联的索引数据类型:Timestamp
/ Datetime
的DatetimeIndex
和Period
的PeriodIndex
。 这些索引数据类型基本上是numpy.ndarray
的子类型,包含对应的时间戳和时间段数据类型,并且可用作序列和数据帧对象的索引。
时间段和时间段索引
Period
数据类型用于表示时间范围或时间跨度。 这里有一些例子:
# Period representing May 2014
In [287]: pd.Period('2014', freq='A-MAY')
Out[287]: Period('2014', 'A-MAY')# Period representing specific day – June 11, 2014
In [292]: pd.Period('06/11/2014')
Out[292]: Period('2014-06-11', 'D')# Period representing 11AM, Nov 11, 1918
In [298]: pd.Period('11/11/1918 11:00',freq='H')
Out[298]: Period('1918-11-11 11:00', 'H')
我们可以向Periods
添加整数,以使时间段提前所需的频率单位数:
In [299]: pd.Period('06/30/2014')+4
Out[299]: Period('2014-07-04', 'D')In [303]: pd.Period('11/11/1918 11:00',freq='H') - 48
Out[303]: Period('1918-11-09 11:00', 'H')
我们还可以计算两个Periods
之间的差,并返回它们之间的频率单位数:
In [304]: pd.Period('2014-04', freq='M')-pd.Period('2013-02', freq='M')
Out[304]: 14
时间段索引
可以通过两种方式创建PeriodIndex
对象,该对象是Period
对象的index
类型。
-
使用
period_range
函数从一些period
对象中,类似于date_range
:In [305]: perRng=pd.period_range('02/01/2014','02/06/2014',freq='D')perRng Out[305]: <class 'pandas.tseries.period.PeriodIndex'>freq: D[2014-02-01, ..., 2014-02-06]length: 6In [306]: type(perRng[:2]) Out[306]: pandas.tseries.period.PeriodIndexIn [307]: perRng[:2] Out[307]: <class 'pandas.tseries.period.PeriodIndex'>freq: D[2014-02-01, 2014-02-02]
正如我们从前面的命令可以确认的那样,当您拉开盖子时,
PeriodIndex
函数实际上就是下面的Period
对象的ndarray
。 -
也可以通过直接调用
Period
构造器来完成:In [312]: JulyPeriod=pd.PeriodIndex(['07/01/2014','07/31/2014'], freq='D')JulyPeriod Out[312]: <class 'pandas.tseries.period.PeriodIndex'>freq: D[2014-07-01, 2014-07-31]
从前面的输出中可以看出,这两种方法之间的差异是period_range
填充了结果ndarray
,但是Period
构造器没有填充,您必须指定索引中应该包含的所有值。
时间序列数据类型之间的转换
我们可以通过to_period
和to_timestamp
函数将Period
和PeriodIndex
数据类型转换为Datetime
/ Timestamp
和DatetimeIndex
数据类型,如下所示:
In [339]: worldCupFinal=pd.to_datetime('07/13/2014', errors='raise')worldCupFinalOut[339]: Timestamp('2014-07-13 00:00:00')In [340]: worldCupFinal.to_period('D')Out[340]: Period('2014-07-13', 'D')In [342]: worldCupKickoff=pd.Period('06/12/2014','D')worldCupKickoff
Out[342]: Period('2014-06-12', 'D')
In [345]: worldCupKickoff.to_timestamp()
Out[345]: Timestamp('2014-06-12 00:00:00', tz=None)In [346]: worldCupDays=pd.date_range('06/12/2014',periods=32, freq='D')worldCupDays
Out[346]: <class 'pandas.tseries.index.DatetimeIndex'>[2014-06-12, ..., 2014-07-13]Length: 32, Freq: D, Timezone: NoneIn [347]: worldCupDays.to_period()
Out[347]: <class 'pandas.tseries.period.PeriodIndex'>freq: D[2014-06-12, ..., 2014-07-13]length: 32
与时间序列相关的对象的摘要
下表总结了与时间序列有关的对象:
对象 | 摘要 |
---|---|
datetime.datetime |
这是一个标准的 Python datetime 类 |
Timestamp |
这是源自datetime.datetime 的 Pandas 类 |
DatetimeIndex |
这是一个 Pandas 类,实现为Timestamp / datetime 对象的不可变numpy.ndarray |
Period |
这是代表某个时间段的 Pandas 类 |
PeriodIndex |
这是一个 Pandas 类,实现为Period 对象的不可变numpy.ndarray |
timedelta |
这是一个 Python 类,表示两个datetime.datetime 实例之间的差异。 实现为datetime.timedelta |
relativedelta |
实现为dateutil.relativedelta 。 dateutil 是标准 Python datetime 模块的扩展。 它提供了额外的功能,例如以大于 1 天的单位表示的时间增量。 |
DateOffset |
这是一个表示常规频率增量的 Pandas 类。 它具有与dateutil.relativedelta 相似的功能。 |
matplotlib 绘图
本节简要介绍了如何使用matplotlib
在 Pandas 中进行绘图。 matplotlib api
使用标准约定导入,如以下命令所示:
In [1]: import matplotlib.pyplot as plt
序列和数据帧有一个plot
方法,它只是plt.plot
的包装。 在这里,我们将研究如何绘制正弦和余弦函数的简单图。 假设我们希望在-pi
到pi
间隔上绘制以下函数:
f(x) = cos(x) + sin(x)
g(x) = cos(x) - sin(x)
这给出了以下间隔:
In [51]: import numpy as np
In [52]: X = np.linspace(-np.pi, np.pi, 256,endpoint=True)In [54]: f,g = np.cos(X)+np.sin(X), np.sin(X)-np.cos(X)
In [61]: f_ser=pd.Series(f)g_ser=pd.Series(g)In [31]: plotDF=pd.concat([f_ser,g_ser],axis=1)plotDF.index=XplotDF.columns=['sin(x)+cos(x)','sin(x)-cos(x)']plotDF.head()
Out[31]: sin(x)+cos(x) sin(x)-cos(x)
-3.141593 -1.000000 1.000000
-3.116953 -1.024334 0.975059
-3.092313 -1.048046 0.949526
-3.067673 -1.071122 0.923417
-3.043033 -1.093547 0.896747
5 rows × 2 columns
现在,我们可以使用plot()
命令和plt.show()
命令绘制数据帧来显示它:
In [94]: plotDF.plot()plt.show()We can apply a title to the plot as follows:
In [95]: plotDF.columns=['f(x)','g(x)']plotDF.plot(title='Plot of f(x)=sin(x)+cos(x), \\n g(x)=sinx(x)-cos(x)')plt.show()
以下是上述命令的输出:
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-Qw60MHz4-1681366316273)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00008.jpeg)]
我们还可以使用以下命令在不同的子图中分别绘制两个序列(函数):
In [96]: plotDF.plot(subplots=True, figsize=(6,6))plt.show()
The following is the output of the preceding command:
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-yT0DdD1N-1681366316274)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00009.jpeg)]
在 Pandas 中使用matplotlib
的绘图函数还有很多。 有关更多信息,请参阅这个链接中的文档。
总结
总而言之,我们讨论了如何处理缺失的数据值以及如何处理 Pandas 中的日期和时间序列。 我们还走了一段简短的弯路,以研究matplotlib
在 Pandas 中的绘图功能。 在准备用于分析和预测的干净数据时,处理缺失的数据起着非常重要的作用,而绘制和可视化数据的能力是每个好的数据分析人员工具箱中必不可少的部分。
在下一章中,我们将对真实数据集进行一些基本数据分析,在其中我们将分析并回答有关数据的基本问题。 有关 Pandas 中这些主题的更多参考,请查看官方文档。
七、统计之旅 – 经典方法
在本章中,我们简要介绍了经典统计数据(也称为常客方法),并展示了如何将 Pandas 与scipy.stats
和statsmodels
等stats
包一起使用来进行统计分析。 本章及其后的内容并不打算作为统计学的入门书,而只是作为与stats
包一起使用 Pandas 的例证。 在下一章中,我们将研究经典观点的另一种方法 – 贝叶斯统计。 本章讨论的各种主题如下:
- 描述统计和推论统计
- 集中趋势和变异性的度量
- 统计假设检验
- Z 检验
- T 检验
- 方差分析
- 置信区间
- 相关和线性回归
描述性统计与推断性统计
在描述性统计或摘要统计中,我们尝试以定量的方式描述数据集合的特征。 这与推论统计或归纳统计不同,因为其目的是对样本进行聚合,而不是使用数据来推断或得出有关样本来源人群的结论。
集中趋势和变异性的度量
描述统计中使用的一些度量包括集中趋势度量和变异性度量。
集中趋势的度量是单个值,它通过指定数据内的中心位置来尝试描述数据集。 集中趋势的三个最常见度量是平均值,中位数和众数。
可变性的度量用于描述数据集中的可变性。 可变性的度量包括方差和标准差。
集中趋势的度量
让我们看一下集中趋势的度量,并在以下各节中进行说明。
均值
均值或样本是最流行的集中趋势度量。 它等于数据集中所有值的总和除以数据集中值的数量。 因此,在 n 个值的数据集中,均值计算如下:
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-0aEmCaRq-1681366316275)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00010.jpeg)]
如果数据值来自样本,则使用 [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-7wmcvjeg-1681366316275)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00011.jpeg)],如果数据值来自总体,则使用μ。
样本平均值和总体平均值不同。 样本均值是所谓的真实总体均值的无偏估计量。 通过对总体进行重复随机抽样以计算样本均值,我们可以获得样本均值的均值。 然后,我们可以调用大数定律和中心极限定理(CLT),并将样本均值的平均值表示为真实总体均值的估计值。
人口平均值也称为人口的期望值。
作为计算值的平均值通常不是数据集中观察到的值之一。 使用均值的主要缺点是,极易受到离群值的影响,或者数据集非常偏斜。 有关其他信息,请参阅以下链接:http://en.wikipedia.org/wiki/Sample_mean_and_sample_covariance,http://en.wikipedia.org/wiki/Law_of_large_numbers 和 http://bit.ly/1bv7l4s。
中位数
中位数是将已排序数据值的集合分为两半的数据值。 它的人口正好在其左侧,而另一半则在其右侧。 在数据集中的值数为偶数的情况下,中位数是两个中间值的平均值。 它受异常值和偏斜数据的影响较小。
众数
众数是数据集中最频繁出现的值。 它通常用于类别数据,以便知道最常见的类别。 使用该众数的一个缺点是它不是唯一的。 具有两种众数的分布称为双峰分布,而具有多种众数的分布称为多峰分布。 这是一个双峰分布的示意图,其中众数分别为两个和七个,因为它们在数据集中都出现了四次:
In [4]: import matplotlib.pyplot as plt%matplotlib inline
In [5]: plt.hist([7,0,1,2,3,7,1,2,3,4,2,7,6,5,2,1,6,8,9,7])plt.xlabel('x')plt.ylabel('Count')plt.title('Bimodal distribution')plt.show()
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-X2Eiyglf-1681366316275)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00012.jpeg)]
用 Python 计算数据集中趋势的度量
为了说明这一点,让我们考虑以下数据集,该数据集由 15 个学生在 20 分中得到的分数获得:
In [18]: grades = [10, 10, 14, 18, 18, 5, 10, 8, 1, 12, 14, 12, 13, 1, 18]
均值,中位数和众数可以如下获得:
In [29]: %precision 3 # Set output precision to 3 decimal places
Out[29]:u'%.3f'In [30]: import numpy as npnp.mean(grades)
Out[30]: 10.933In [35]: %precisionnp.median(grades)
Out[35]: 12.0In [24]: from scipy import statsstats.mode(grades)
Out[24]: (array([ 10.]), array([ 3.]))
In [39]: import matplotlib.pyplot as plt
In [40]: plt.hist(grades)plt.title('Histogram of grades')plt.xlabel('Grade')plt.ylabel('Frequency')plt.show()
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-U0ACMTm6-1681366316275)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00013.jpeg)]
为了说明数据的偏斜度或离群值会如何严重影响均值作为衡量集中趋势的效用,请考虑以下数据集,该数据集显示了工厂员工的工资(以千美元计):
In [45]: %precision 2salaries = [17, 23, 14, 16, 19, 22, 15, 18, 18, 93, 95]In [46]: np.mean(salaries)
Out[46]: 31.82
基于平均值,我们可以假设数据以31.82
的平均值为中心。 但是,我们会错的。 为了看到这一点,让我们使用条形图显示数据的经验分布:
In [59]: fig = plt.figure()ax = fig.add_subplot(111)ind = np.arange(len(salaries))width = 0.2plt.hist(salaries, bins=xrange(min(salaries),max(salaries)).__len__())ax.set_xlabel('Salary')ax.set_ylabel('# of employees')ax.set_title('Bar chart of salaries')plt.show()
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-VKevmLnx-1681366316276)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00014.jpeg)]
从前面的条形图中,我们可以看到大多数薪水远低于 30K,并且没有人接近 32K 的平均值。 现在,如果我们看一下中位数,我们会发现在这种情况下,它是衡量集中趋势的更好方法:
In [47]: np.median(salaries)
Out[47]: 18.00
我们还可以看一下数据的直方图:
In [56]: plt.hist(salaries, bins=len(salaries))plt.title('Histogram of salaries')plt.xlabel('Salary')plt.ylabel('Frequency')plt.show()
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-vk77pSKM-1681366316276)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00015.jpeg)]
注意
直方图实际上是数据的更好表示,因为条形图通常用于表示类别数据,而直方图是定量数据的首选,而薪金数据就是这种情况。
有关何时使用直方图和条形图的更多信息,请参见这里。
如果分布是对称的且是单峰的(即只有一种众数),则三个度量(均值,中位数和众数)将相等。 如果分配偏斜则不是这种情况。 在这种情况下,均值和中位数会彼此不同。 对于负偏斜分布,对于正偏斜分布,均值将低于中位数,反之亦然:
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-rpqOX4yB-1681366316276)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00016.jpeg)]
上图来自这里。
度量值的变异性,分散性或差距
我们在描述性统计中测量的分布的另一个特征是可变性。
可变性指定数据点彼此不同或分散多少。 可变性度量很重要,因为它们可以提供对集中趋势度量无法提供的数据性质的洞察。
例如,假设我们进行了一项研究,以检验学前教育计划在提高经济弱势儿童的考试成绩方面的效果。 我们不仅可以通过整个样本的检验分数的平均值来衡量有效性,而且可以通过分数的分散性来衡量有效性。 它对某些学生有用,而对其他学生却没有那么多? 数据的可变性可以帮助我们确定要采取的一些步骤,以提高程序的实用性。
范围
色散的最简单度量是范围。 范围是数据集中最低和最高得分之间的差。 这是传播最简单的方法。
范围=最大值-最小值
四分位数
四分位数和相关的四分位数范围是分散的更重要度量。 它也代表四分位数,这意味着它是度量标准上的值,低于该值时已排序数据集中的分数的 25%,50%,75% 和 100% 下降。 四分位数是将数据集分为四组的三个点,每个组包含四分之一的数据。 为了说明这一点,假设我们有一个包含 20 个检验成绩的数据集,其中对它们的排名如下:
In [27]: import randomrandom.seed(100)testScores = [random.randint(0,100) for p in xrange(0,20)]testScores
Out[27]: [14, 45, 77, 71, 73, 43, 80, 53, 8, 46, 4, 94, 95, 33, 31, 77, 20, 18, 19, 35]In [28]: #data needs to be sorted for quartilessortedScores = np.sort(testScores)
In [30]: rankedScores = {i+1: sortedScores[i] for i in xrange(len(sortedScores))}In [31]: rankedScores
Out[31]:
{1: 4,2: 8,3: 14,4: 18,5: 19,6: 20,7: 31,
8: 33,9: 35,10: 43,11: 45,12: 46,13: 53,14: 71,15: 73,16: 77,17: 77,18: 80,19: 94,20: 95}
第一个四分位数(Q1)在第五和第六个分数之间,第二个四分位数(Q2)在第十和第十一个分数之间,第三个四分位数在第十五和第十六个分数之间。 因此,我们有了(通过使用线性插值并计算中点):
Q1 = (19+20)/2 = 19.5
Q2 = (43 + 45)/2 = 44
Q3 = (73 + 77)/2 = 75
要在 IPython 中看到这一点,我们可以使用scipy.stats
或numpy.percentile
包:
In [38]: from scipy.stats.mstats import mquantilesmquantiles(sortedScores)
Out[38]: array([ 19.45, 44\\. , 75.2 ])In [40]: [np.percentile(sortedScores, perc) for perc in [25,50,75]]
Out[40]: [19.75, 44.0, 74.0]
值与我们先前的计算不完全匹配的原因是由于插值方法不同。 可在这个链接中找到有关获取四分位数值的各种方法的更多信息。 四分位数范围是从第三四分位数(Q3-Q1)中减去的第一四分位数,它表示数据集中的中位数 50。
有关更多信息,请参见这里。
有关scipy.stats
和numpy.percentile
函数的更多详细信息,请参见这个文档和这个文档。
偏差和方差
讨论可变性的基本思想是偏差的概念。 简而言之,偏差测度告诉我们给定值与分布平均值 [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-cBtG54YV-1681366316277)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00017.jpeg)] 的距离有多远。
为了找到一组值的偏差,我们将方差定义为偏差平方的总和,然后将其除以数据集的大小对其进行归一化。 这称为方差。 由于负偏差和正偏差彼此抵消,因此我们需要使用偏差的平方和作为围绕均值结果的偏差之和。 偏差的平方和定义如下:
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-bfZBSJXF-1681366316277)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00018.jpeg)]
可以证明前面的表达式等效于:
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-2STVSN0n-1681366316277)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00019.jpeg)]
正式地,方差定义如下:
-
对于样本方差,请使用以下公式:
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-X07loHMe-1681366316278)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00020.jpeg)]
-
对于总体方差,使用以下公式:
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-G9BJlZwS-1681366316278)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00021.jpeg)]
对于样本方差而不是 [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-yNC4RKce-1681366316278)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00023.jpeg)],分母为 [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-S3qv0rEk-1681366316279)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00022.jpeg)] 的原因是,对于样本方差,我们希望使用无偏估计量。 有关更多信息,请查看这里。
此度量的值以平方单位表示。 这强调了一个事实,我们计算出的方差是平方偏差。 因此,要获得与数据集原始点相同单位的偏差,我们必须取平方根,这就是我们所说的标准差。 因此,使用以下公式可以得出样品的标准差:
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-AcfS8YXi-1681366316279)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00024.jpeg)]
但是,对于总体,标准差由以下公式给出:
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-xTPEHGYA-1681366316279)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00025.jpeg)]
假设检验 -– 原假设和替代假设
在上一节中,我们简要介绍了所谓的描述性统计。 在本节中,我们将讨论所谓的推论统计,由此我们尝试使用样本数据集的特征来得出有关更广泛总体的结论。
推论统计中最重要的方法之一就是假设检验。 在假设检验中,我们尝试确定某个假设或研究问题在一定程度上是否正确。 一种假设的例子是:吃菠菜可以改善长期记忆。
为了通过假设检验调查此问题,我们可以选择一组人作为我们的研究对象,并将其分为两组或样本。 第一组是实验组,它将在预定的时间内吃菠菜。 不接受菠菜的第二组为对照组。 在选定的时间段内,将对两组中的个人记忆进行测量和统计。
我们在实验结束时的目标是能够做出这样的陈述,例如“吃菠菜可以改善长期记忆,而这并非偶然”。 这也称为重要性。
在前面的场景中,研究中的主题集合称为样本,而我们希望得出结论的一般人群是人口。
我们研究的最终目标是确定我们在样本中观察到的任何影响是否可以推广到整个人群。 为了进行假设检验,我们将需要提出原假设和替代假设。
原假设和替代假设
通过参考前面的菠菜示例,原假设为:吃菠菜对长期记忆性能没有影响。
原假设就是,它使我们试图通过运行实验来证明所证明的意图无效。 它通过断言某个统计指标(稍后将进行解释)为零来做到这一点。
替代假设是我们希望支持的假设。 这与原假设相反,我们假设它是正确的,直到数据提供了足够的证据表明相反。 因此,在这种情况下,我们的另一种假设是:吃菠菜可改善长期记忆。
象征性地,原假设称为 H0 ,替代假设称为 H1 。 您可能希望将先前的原假设和替代假设重新陈述为对我们的研究更为具体和可衡量的内容。 例如,我们可以按以下方式重铸 H0:
每天吃 40 克菠菜 90 天的 1,000 名受试者的平均记忆评分与同期未食用菠菜的 1,000 名受试者的对照组没有什么不同。
在进行实验/研究时,我们专注于试图证明或否定原假设。 这是因为我们可以计算结果归因于偶然性的概率。 但是,没有一种简单的方法来计算替代假设的可能性,因为长期记忆的任何改善都可能是由于吃菠菜以外的其他因素引起的。
我们通过假设原假设为真来检验原假设,并单独计算偶然获得的结果的概率。 我们设置一个阈值水平 – α – 对于该阈值水平,如果计算的概率较小,则可以拒绝原假设,如果阈值级别较大,则可以接受。 拒绝原假设等于接受替代假设,反之亦然。
alpha 和 p 值
为了进行实验以决定是否支持我们的原假设,我们需要提出一种方法,使我们能够以具体且可衡量的方式做出决定。 要进行此显着性检验,我们必须考虑两个数字 – 检验统计量的 p 值和显着性阈值水平,也称为 alpha。
如果我们通过假设无效假设为真或仅凭偶然发生而观察到的结果,则 p 值为概率。
假设无效假设为真,则 p 值也可以认为是获得检验统计量的概率与实际获得的检验统计量一样极端或更高。
alpha 值是我们与 p 值进行比较的阈值。 这为我们接受或拒绝原假设提供了一个切入点。 它是衡量我们观察到的结果必须多么极端才能拒绝实验的原假设的一种度量。 alpha 的最常用值为 0.05 或 0.01。
通常,规则如下:
如果 p 值小于或等于 alpha(p < .05
),则我们拒绝原假设,并声明结果具有统计学意义。
如果 p 值大于 alpha(p > .05
),则我们无法拒绝原假设,并且我们说结果在统计上不显着。
在使用中看似随意的 alpha 值是常识性方法的缺点之一,与此方法有关的问题很多。 Nature 杂志上的文章重点介绍了一些问题。
有关此主题的更多详细信息,请参阅:
- http://statistics.about.com/od/Inferential-Statistics/a/What-Is-The-Difference-Between-Alpha-And-P-Values.htm
- http://bit.ly/1GzYX1P
- http://en.wikipedia.org/wiki/P-value
I 类和 II 类错误
错误有两种,如下所述:
- I 类错误:在这种类型的错误中,当 H0 实际上为真时,我们拒绝 H0 。 陪审团的一个例子就是陪审团裁定无辜者犯有该人没有犯下的罪行。
- II 类错误:在这种类型的错误中,当 H1 实际上为真时,我们无法拒绝 H0 。 这相当于有罪的人逃脱定罪。
统计假设检验
统计假设检验是一种使用统计研究或实验中的数据做出决策的方法。 在统计中,基于预定的阈值概率或显着性水平,如果不可能仅偶然发生的结果称为统计显着。 统计检验分为两类:单侧检验和双侧检验。
在双侧检验中,我们分配一半的 alpha 值来检验一个方向上的统计显着性,而另一半则分配给另一方向上的统计显着性。
在单侧检验中,检验仅在一个方向上执行。
有关此主题的更多详细信息,请参考这里。
背景
要应用统计推断,重要的是要了解所谓的抽样分布的概念。 假设我们从原假设成立的总体中随机抽样,则抽样分布是统计量所有可能值及其概率的集合。
更为简单的定义是:抽样分布是统计数据可以假设(分布)的一组值,如果我们要从总体中重复抽取样本及其相关概率。
统计量的值是来自统计量抽样分布的随机样本。 通过获取许多不同大小的样本并取其平均值来计算平均值的采样分布。 它的平均值 [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-9udv7Ule-1681366316280)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00026.jpeg)] 等于 [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-22TzLqLp-1681366316280)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00027.jpeg)],标准差 [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-XBncvyt7-1681366316281)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00028.jpeg)] [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-GJtvfYEI-1681366316281)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00029.jpeg)] 等于 [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-rrNCuaTr-1681366316281)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00030.jpeg)] [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-ZnUjbRMB-1681366316281)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00031.jpeg)]。
CLT 指出,如果原始或原始评分总体呈正态分布,或者样本量足够大,则抽样分布呈正态分布。 传统上,统计人员将足够大的样本数量表示为 [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-XeSWndt9-1681366316282)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00032.jpeg)],即样本数量为 30 或更大。 但是,这仍然是一个争论的话题。
有关此主题的更多详细信息,请参阅这里和这里。
采样分布的标准差通常称为平均值的标准误差或仅称为标准误差。
Z 检验
z 检验适用于以下情况:
- 该研究仅涉及一个样本均值,并且原假设人口的参数 [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-SpAmYJHc-1681366316282)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00033.jpeg)] 和 [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-WfJyabom-1681366316282)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00034.jpeg)] 是已知的
- 均值的采样分布为正态分布
- 样本大小为 [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-XWEQbKQp-1681366316282)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00035.jpeg)]
当总体平均值为已知时,我们使用 Z 检验。 在 Z 检验中,我们问一个问题:总体平均值 [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-x1hKMLXg-1681366316283)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00033.jpeg)] 是否不同于假设值。 Z 检验的原假设如下:
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-ZigLBe3R-1681366316283)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00036.jpeg)]
其中 [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-UzmJjAcF-1681366316284)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00037.jpeg)] 人口平均值
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-1OTp64RA-1681366316284)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00038.jpeg)] =假设值
替代假设 [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-uJgCuS6L-1681366316284)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00039.jpeg)] 可以是以下之一:
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-kZWOT7Yx-1681366316284)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00040.jpeg)]
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-CktZjR7d-1681366316285)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00041.jpeg)]
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-NOxTx9YG-1681366316285)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00042.jpeg)]
前两个是单侧检验,最后一个是双侧检验。 具体来说,为了检验 [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-ZijaQEAM-1681366316285)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00043.jpeg)],我们计算了检验统计量:
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-H1l30Rzs-1681366316285)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00044.jpeg)]
在此,[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-Z52R7cEF-1681366316286)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00029.jpeg)] 是 [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-qNbAIMtE-1681366316286)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00045.jpeg)] 的采样分布的真实标准差。 如果 [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-3apgHxzn-1681366316286)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00046.jpeg)] 为true
,则 Z 检验统计量将具有标准正态分布。
在这里,我们给出 Z 检验的快速说明。
假设我们有一家虚构的公司 Intelligenza,该公司声称他们提出了一种根本的新方法来改善内存保留和学习。 他们声称他们的技术可以提高成绩,而不是传统的学习技术。 假设使用传统的学习技术,成绩提高 40% ,标准差提高 10%。
使用 Intelligenza 方法对 100 名学生进行了随机检验,结果平均提高了 44%。 Intelligenza 的主张是否成立?
本研究的原假设指出,与传统的学习技术相比,使用 Intelligenza 的方法在成绩上没有任何改善。 另一种假设是,使用 Intelligenza 的方法比传统的研究技术有所改进。
原假设由以下给出:
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-qhoyfX9e-1681366316286)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00036.jpeg)]
备用假设由以下给出:
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-4FIyiayk-1681366316286)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00041.jpeg)]
std = 10 / sqrt(100) = 1
z = (43.75-40) / (10/10) = 3.75 std
回想一下,如果原假设为真,则检验统计量 z 将具有如下所示的标准正态分布:
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-GmIKbQNK-1681366316287)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00047.jpeg)]
作为参考,请访问这里。
z
的值是来自标准正态分布的随机样本,如果原假设为真,则为 z 的分布。
z = 43.75
的观测值对应于标准正态分布曲线上的极端离群 p 值,远小于 0.1%。
p 值是曲线下的面积,在前一条正态分布曲线上的 3.75 值的右侧。
这表明,如果我们从标准正态分布中进行抽样,那么我们极不可能获得检验统计量的观察值。
我们可以使用scipy.stats
包,使用 Python 查找实际的 p 值,如下所示:
In [104]: 1 - stats.norm.cdf(3.75)
Out[104]: 8.841728520081471e-05
因此,[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-WanY9SgD-1681366316287)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00048.jpeg)],即,如果检验统计量呈正态分布,则获得观测值的概率为 8.8e-05,接近零。 因此,如果原假设是真实的,几乎不可能获得我们观察到的值。
用更正式的术语来说,我们通常将定义一个阈值或 alpha 值,如果p ≤ α
则拒绝原假设,否则将拒绝。
alpha 的典型值为 0.05 或 0.01。 以下列表说明了 alpha 的不同值:
p < 0.01
:反对 H0 的非常有力的证据0.01 < p < 0.05
:反对 H0 的有力证据0.05 < p < 0.1
:反对 H0 的证据不充分p > 0.1
:几乎没有反对 H0 的证据
因此,在这种情况下,我们将拒绝原假设,并相信 Intelligenza 的主张,并声明其主张具有重大意义。 在这种情况下,反对原假设的证据很重要。 我们使用两种方法来确定是否拒绝原假设:
- P 值法
- 拒绝区域法
我们在前面的示例中使用的方法是后者。
p 值越小,原假设为真的可能性就越小。 在拒绝区域方法中,我们具有以下规则:
如果 [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-3jgYheBh-1681366316287)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00049.jpeg)],则拒绝原假设,否则保留原假设。
T 检验
当总体的标准差已知时,Z 检验最有用。 但是,在大多数实际情况下,这是一个未知数。 对于这些情况,我们转向显着性的 T 检验。
对于 T 检验,假设总体的标准差未知,我们将其替换为样本的标准差s
。 平均值的标准误差现在变为:
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-4550alxP-1681366316287)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00049.jpeg)]
样本s
的标准差计算如下:
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-0SS9ZNBa-1681366316288)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00050.jpeg)]
分母是N-1
而不是N
。 该值称为自由度数。 我现在将不加说明地指出,通过 CLT,随着N
的出现,T 分布近似于正态,瓜斯或 Z 分布,因此N-1
增大,即随着自由度(df
)。 当df = ∞
时,T 分布与正态分布或 Z 分布相同。 这是直观的,因为随着df
的增加,样本量也会增加,并且 s 接近 [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-TFbBn6KN-1681366316288)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00034.jpeg)],即总体的真实标准差。 存在无限数量的 T 分布,每个 T 分布对应于不同的df
值。
如下图所示:
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-bvDKq5vb-1681366316288)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00051.jpeg)]
该图像的引用来自这里。
有关 T 分布,Z 分布和自由度之间关系的更详细的技术说明,请参见这里。
T 检验的类型
有各种类型的 T 检验。 以下是最常见的; 他们通常会制定一个原假设,以对分布的均值提出要求:
-
单样本独立 T 检验:用于将样本的平均值与已知总体平均值或已知值的平均值进行比较。 假设我们是澳大利亚的健康研究人员,他们关注土著人口的健康,并希望确定低收入土著母亲所生的婴儿出生体重是否低于正常体重。
单样本 T 检验的原假设检验的一个例子是:从低收入土著母亲的 150 例足月活产婴儿分娩中,我们样本的平均出生体重,与澳大利亚一般人口的婴儿平均出生体重无异,即 3367 克。
该信息的参考是这里。
-
独立样本 T 检验:用于将独立样本的均值进行比较。 独立样本 T 检验的一个示例是自动变速器与手动变速器车辆的燃油经济性比较。 这是我们实际示例中要重点关注的内容。
T 检验的原假设是:手动和自动变速箱汽车的平均燃油效率在城市/高速公路平均行驶里程之间没有差异。
-
配对样本 T 检验:在成对/相关样本 T 检验中,我们采用有意义的方式将一个样本中的每个数据点与另一个样本中的数据点配对。 一种方法是在不同的时间点对同一样本进行测量。 这样的一个例子是通过比较饮食前后受试者的体重来检查减肥饮食的功效。
在这种情况下,原假设是这样的:减肥饮食前后参与者的平均体重之间没有差异,或更简洁地说,配对观察之间的平均差异为零。
可以在这里上找到此信息的参考。
T 检验示例
简而言之,要进行空有效假设检验(NHST),我们需要执行以下操作:
- 制定我们的原假设。 原假设是我们的系统模型,假设我们希望验证的效果实际上是偶然的。
- 计算我们的 p 值。
- 将计算出的 p 值与我们的 alpha 值或阈值进行比较,并决定是否拒绝或接受原假设。 如果 p 值足够低(低于 alpha),我们将得出结论,即原假设可能是不正确的。
对于我们的实际插图,我们希望调查手动变速箱车辆是否比自动变速箱更省油。 为此,我们将利用美国政府 2014 年发布的燃料经济数据。
In [53]: import pandas as pdimport numpy as npfeRawData = pd.read_csv('2014_FEGuide.csv')In [54]: feRawData.columns[:20]
Out[54]: Index([u'Model Year', u'Mfr Name', u'Division', u'Carline', u'Verify Mfr Cd', u'Index (Model Type Index)', u'Eng Displ', u'# Cyl', u'Trans as listed in FE Guide (derived from col AA thru AF)', u'City FE (Guide) - Conventional Fuel', u'Hwy FE (Guide) - Conventional Fuel', u'Comb FE (Guide) - Conventional Fuel', u'City Unadj FE - Conventional Fuel', u'Hwy Unadj FE - Conventional Fuel', u'Comb Unadj FE - Conventional Fuel', u'City Unrd Adj FE - Conventional Fuel', u'Hwy Unrd Adj FE - Conventional Fuel', u'Comb Unrd Adj FE - Conventional Fuel', u'Guzzler? ', u'Air Aspir Method'], dtype='object')In [51]: feRawData = feRawData.rename(columns={'Trans as listed in FE Guide (derived from col AA thru AF)' :'TransmissionType','Comb FE (Guide) - Conventional Fuel' : 'CombinedFuelEcon'})In [57]: transType=feRawData['TransmissionType']transType.head()
Out[57]: 0 Auto(AM7)1 Manual(M6)2 Auto(AM7)3 Manual(M6)4 Auto(AM-S7)Name: TransmissionType, dtype: object
现在,我们希望修改前面的序列,以使这些值仅包含Auto
和Manual
字符串。 我们可以这样做,如下所示:
In [58]: transTypeSeries = transType.str.split('(').str.get(0)transTypeSeries.head()
Out[58]: 0 Auto1 Manual2 Auto3 Manual4 AutoName: TransmissionType, dtype: object
现在,我们从序列中创建一个最终的修改后的数据帧,该数据帧由变速箱类型和组合的燃油经济性数字组成:
In [61]: feData=pd.DataFrame([transTypeSeries,feRawData['CombinedFuelEcon']]).TfeData.head()
Out[61]: TransmissionType CombinedFuelEcon0 Auto 161 Manual 152 Auto 163 Manual 154 Auto 175 rows × 2 columns
现在,我们可以将自动变速箱的数据与手动变速箱的数据分开,如下所示:
In [62]: feData_auto=feData[feData['TransmissionType']=='Auto']feData_manual=feData[feData['TransmissionType']=='Manual']
In [63]: feData_auto.head()
Out[63]: TransmissionType CombinedFuelEcon0 Auto 162 Auto 164 Auto 176 Auto 168 Auto 175 rows × 2 columns
这表明有 987 辆具有自动变速箱的车辆,而有 211 辆具有手动变速箱的车辆:
In [64]: len(feData_auto)
Out[64]: 987In [65]: len(feData_manual)
Out[65]: 211In [87]: np.mean(feData_auto['CombinedFuelEcon'])
Out[87]: 22.173252279635257In [88]: np.mean(feData_manual['CombinedFuelEcon'])
Out[88]: 25.061611374407583In [84]: import scipy.stats as statsstats.ttest_ind(feData_auto['CombinedFuelEcon'].tolist(), feData_manual['CombinedFuelEcon'].tolist())
Out[84]: (array(-6.5520663209014325), 8.4124843426100211e-11)In [86]: stats.ttest_ind(feData_auto['CombinedFuelEcon'].tolist(), feData_manual['CombinedFuelEcon'].tolist(), equal_var=False)
Out[86]: (array(-6.949372262516113), 1.9954143680382091e-11)
置信区间
在本节中,我们将讨论置信区间的问题。 置信区间使我们能够对总体给定样本数据的平均值进行概率估计。
此估计称为间隔估计,由一些值(间隔)组成,这些值充当未知总体参数的良好估计。
置信区间受置信限度限制。 95% 置信区间定义为一个区间,其中该区间包含概率为 95% 的总体平均值。 那么,我们如何构建置信区间?
假设我们有一个双侧 T 检验,我们希望构建一个 95% 的置信区间。 在这种情况下,我们希望与平均值对应的样本 t 值 [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-xJU5BAxf-1681366316288)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00052.jpeg)] 满足以下不等式:
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-nEMbGRjb-1681366316288)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00053.jpeg)]
给定 [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-22LahfVm-1681366316289)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00054.jpeg)],我们可以用前面的不等式关系来代替它:
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-VKhfJUsc-1681366316289)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00055.jpeg)]
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-G0qo44MG-1681366316289)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00056.jpeg)] 区间是我们的 95% 置信区间。
概括而言,任何百分比 y 的任何置信区间都可以表示为 [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-UZB2JhwN-1681366316289)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00057.jpeg)],其中 [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-1kUTzllC-1681366316290)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00058.jpeg)] 是 [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-bXKD8LoP-1681366316290)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00059.jpeg)] 的 t 尾值,即 [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-ugkEg9bd-1681366316290)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00060.jpeg)] 对应于 y 的期望置信区间。
现在,我们将借此机会说明如何使用流行的统计环境 R 中的数据集来计算置信区间。stats
模型的模块提供了对 R 核心数据包中可通过[ get_rdataset
函数。
一个说明性示例
我们将考虑称为faithful
的数据集,该数据集包含通过观察美国黄石国家公园的老忠实间歇泉的爆发而获得的数据。数据集中的两个变量是爆发,即间歇泉爆发的时间长度,然后是等待直到下一次爆发的时间间隔。 有 272 个观察结果。
In [46]: import statsmodels.api as smafaithful=sma.datasets.get_rdataset("faithful")faithful
Out[46]: <class 'statsmodels.datasets.utils.Dataset'>In [48]: faithfulDf=faithful.datafaithfulDf.head()
Out[48]: eruptions waiting0 3.600 791 1.800 542 3.333 743 2.283 624 4.533 85
5 rows × 2 columnsIn [50]: len(faithfulDf)
Out[50]: 272
让我们计算间歇喷泉平均等待时间的 95% 置信区间。 为此,我们首先获取数据的样本均值和标准差:
In [80]: mean,std=(np.mean(faithfulDf['waiting']),np.std(faithfulDf['waiting']))
现在,我们使用scipy.stats
包来计算置信区间:
In [81]: from scipy import statsN=len(faithfulDf['waiting'])ci=stats.norm.interval(0.95,loc=mean,scale=std/np.sqrt(N))
In [82]: ci
Out[82]: (69.28440107709261, 72.509716569966201)
因此,我们可以说[69.28,72.51]区间包含了间歇泉的实际平均等待时间,具有 95% 的置信度。
此信息的参考:这里和这里。
相关和线性回归
确定两个变量之间关系的统计中最常见的任务之一是它们之间是否存在依赖关系。 相关性是我们在统计数据中用于表示相互依赖的变量的总称。
然后,我们可以使用这种关系来尝试从另一组变量预测一组变量的值; 这称为回归。
相关性
用相关关系表示的统计依存关系并不意味着两个变量之间存在因果关系。 与此相关的著名观点是“相关并不意味着因果关系”。 因此,两个变量或数据集之间的相关性仅表示是偶然的,而不是因果关系或依赖性。 例如,在特定日期购买的冰淇淋数量与天气之间存在相关性。
有关相关性和依赖性的更多信息,请参见这里。
相关度量(称为相关系数)是一个数字,用于捕获两个变量之间关系的大小和方向。 方向的范围从 -1 到 +1,幅度的范围从 0 到 1。 关系的方向通过符号表示,其中+
符号表示正相关,而-
符号表示负相关。 幅度越高,与被称为完美相关的相关性就越大。
最受欢迎和使用最广泛的相关系数是皮尔逊积矩相关系数,称为 r 。 它测量两个 x 和 y 变量之间的线性相关性或依赖性,并获取 -1 和 +1 之间的值。
样本相关系数 r 定义如下:
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-9UpPPoAu-1681366316290)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00061.jpeg)]
也可以这样写:
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-4jTP8hr4-1681366316291)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00062.jpeg)]
在此,我们省略了求和限制。
线性回归
如前所述,回归集中在使用两个变量之间的关系进行预测。 为了使用线性回归进行预测,必须计算最适合的直线。
如果所有点(变量的值)都位于一条直线上,则该关系被认为是完美的。 在实践中,这种情况很少发生,并且点并非都整齐地排列在一条直线上。 那么这种关系是不完美的。 在某些情况下,线性关系仅在对数转换的变量之间发生。 这是一个日志日志模型。 这种关系的一个例子是物理学中的幂定律分布,其中一个变量随另一变量的幂而变化。
因此,诸如 [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-zuK3dCs0-1681366316291)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00063.jpeg)] 的表达式导致 [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-PjehchAc-1681366316291)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00064.jpeg)] 线性关系。
有关更多信息,请参见这里
为了构建最适合的线,使用最小二乘法。 在这种方法中,最佳拟合线是在各点之间构造的最佳线,对于这些点,从每个点到该线的平方距离的总和最小。 这被认为是我们尝试使用线性回归建模的变量之间关系的最佳线性近似。 在这种情况下,最佳拟合线称为最小二乘回归线。
更正式地说,最小二乘回归线是对于从数据点到该线的垂直距离的平方和具有最小可能值的线。 这些垂直距离也称为残差。
因此,通过构造最小二乘回归线,我们试图最小化以下表达式:
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-bHjCGvgm-1681366316292)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00065.jpeg)]
一个说明性示例
现在,我们将通过一个示例来说明所有上述要点。 假设我们正在进行一项研究,其中我们想说明温度对发声频率的影响。 该示例的数据从 1948 年的《昆虫之歌》,乔治·皮尔斯中获得。乔治·皮尔斯(George Pierce)测量了在各种温度下昆虫发出的声波的频率。
我们希望调查的频率和温度,因为我们怀疑这是有关系的。 数据由 16 个数据点组成,我们将其读入数据帧:
In [38]: import pandas as pdimport numpy as npchirpDf= pd.read_csv('cricket_chirp_temperature.csv')
In [39]: chirpDf
Out[39]:chirpFrequency temperature
0 20.000000 88.599998
1 16.000000 71.599998
2 19.799999 93.300003
3 18.400000 84.300003
4 17.100000 80.599998
5 15.500000 75.199997
6 14.700000 69.699997
7 17.100000 82.000000
8 15.400000 69.400002
9 16.200001 83.300003
10 15.000000 79.599998
11 17.200001 82.599998
12 16.000000 80.599998
13 17.000000 83.500000
14 14.400000 76.300003
15 rows × 2 columns
首先,让我们绘制数据散点图以及回归线或最佳拟合线:
In [29]: plt.scatter(chirpDf.temperature,chirpDf.chirpFrequency,marker='o',edgecolor='b',facecolor='none',alpha=0.5)plt.xlabel('Temperature')plt.ylabel('Chirp Frequency')slope, intercept = np.polyfit(chirpDf.temperature,chirpDf.chirpFrequency,1)plt.plot(chirpDf.temperature,chirpDf.temperature*slope + intercept,'r')plt.show()
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-0sYLR4SD-1681366316292)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00066.jpeg)]
从图中可以看出,温度和线性调频频率之间似乎存在线性关系。 现在,我们可以使用statsmodels.ols
(普通最小二乘法)方法进一步进行调查:
[37]: chirpDf= pd.read_csv('cricket_chirp_temperature.csv')chirpDf=np.round(chirpDf,2)result=sm.ols('temperature ~ chirpFrequency',chirpDf).fit()result.summary()Out[37]: OLS Regression ResultsDep. Variable: temperature R-squared: 0.697Model: OLS Adj. R-squared: 0.674Method: Least Squares F-statistic: 29.97Date: Wed, 27 Aug 2014 Prob (F-statistic): 0.000107Time: 23:28:14 Log-Likelihood: -40.348No. Observations: 15 AIC: 84.70Df Residuals: 13 BIC: 86.11Df Model: 1 coef std err t P>|t| [95.0% Conf. Int.]Intercept 25.2323 10.060 2.508 0.026 3.499 46.966chirpFrequency 3.2911 0.601 5.475 0.000 1.992 4.590Omnibus: 1.003 Durbin-Watson: 1.818Prob(Omnibus): 0.606 Jarque-Bera (JB): 0.874Skew: -0.391 Prob(JB): 0.646Kurtosis: 2.114 Cond. No. 171.
除了R-squared
,Intercept
和chirpFrequency
值外,我们将忽略大多数上述结果。
从前面的结果可以得出,回归线的斜率为 3.29,温度轴上的截距为 25.23。 因此,回归线方程如下所示:temp = 25.23 + 3.29 * chirpFrequency
。
这意味着,随着线性调频频率增加一,温度会升高约 3.29 华氏度。 但是,请注意,截距值实际上并不有意义,因为它超出了数据的范围。 我们也只能预测数据范围内的值。 例如,我们无法预测 32 华氏度下的chirpFrequency
是什么,因为它超出了数据的范围; 而且,在华氏 32 度时,昆虫会冻死。 相关系数 R 的值如下所示:
In [38]: R=np.sqrt(result.rsquared)R
Out[38]: 0.83514378678237422
因此,我们的相关系数为R = 0.835
。 这表明温度的变化可以解释频率的约 84%。
该信息的参考:昆虫之歌
数据来自这里。
有关单变量和多变量回归的更深入的处理,请访问以下网站:
- 回归(第一部分)
- 回归(第二部分)
总结
在本章中,我们简要介绍了统计的经典方法或常客方法,并向您展示了如何将 Pandas 与stats
和scipy.stats
和statsmodels
包结合在一起,以计算,解释和推断统计数据。
在下一章中,我们将研究另一种统计方法,即贝叶斯方法。 要更深入地了解我们涉及的统计信息主题,请查看了解行为科学中的统计信息。
八、贝叶斯统计简介
在本章中,我们将简要介绍称为贝叶斯统计的另一种统计推断方法。 它不打算作为完整的入门手册,而只是对贝叶斯方法的介绍。 我们还将探讨相关的 Python 相关库,如何使用 Pandas 以及matplotlib
来帮助进行数据分析。 将讨论的各种主题如下:
- 贝叶斯统计概论
- 贝叶斯统计的数学框架
- 概率分布
- 贝叶斯统计与频率统计
- PyMC 和蒙特卡洛仿真简介
- 贝叶斯推理的示例 – 切换点检测
贝叶斯统计概论
贝叶斯统计领域建立在 18 世纪统计学家,哲学家和长老会牧师托马斯·贝叶斯牧师的工作基础上。 他著名的贝叶斯定理(构成贝叶斯统计的理论基础)于 1763 年死后出版,以解决逆概率问题。 有关此主题的更多详细信息,请参考这里。
逆概率问题在 18 世纪初期风行一时,通常表述如下:
假设您和朋友一起玩游戏。 袋子 1 中有 10 个绿色球和 7 个红色球,袋子 2 中有 4 个绿色和 7 个红色球。您的朋友背离您的视线,扔硬币并从其中一个袋子中随机捡起一个球,并显示给您。 球是红色的。 球从袋子 1 中抽出的概率是多少?
这些问题被称为逆概率问题,因为我们试图根据随后的事件(球是红色)来估计已经发生的事件(球从哪个袋子抽出)的概率。
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-fLVDAOHq-1681366316292)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00067.jpeg)]
贝叶斯球图
让我们快速说明一下如何解决前面说明的逆概率问题。 考虑到球是红色的,我们希望计算从袋子 1 中抽出球的概率。 这可以表示为 [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-v0Z51k3u-1681366316293)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00068.jpeg)]。
让我们从计算选择红球的概率开始。 如上图所示,可以通过遵循红色的两条路径来计算。 因此,我们有 [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-um0BDvfo-1681366316293)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00069.jpeg)]。
现在,从袋子 1 中选择红色球的概率仅是通过上方路径,并给出如下:
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-xkvf4K1u-1681366316293)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00070.jpeg)]
并且,从袋 2 中选择红色球的概率如下所示:
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-fhokt0hR-1681366316293)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00071.jpeg)]
注意,该概率可以写成如下形式:
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-IXEcKDpJ-1681366316294)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00072.jpeg)]
通过检查我们可以看到 [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-mAmZXgDF-1681366316294)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00073.jpeg)] 和树的最后一个分支只有在球首先位于袋 1 中并且是红色球时才被遍历。 因此,从直观上讲,我们将得到以下结果:
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-TkugYiI9-1681366316294)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00074.jpeg)]
贝叶斯统计的数学框架
利用贝叶斯方法,我们提出了进行统计推断的另一种方法。 我们首先介绍贝叶斯定理,这是从其得出所有贝叶斯推断的基本方程式。
关于概率的几个定义是有序的:
- [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-EGRVWSsR-1681366316295)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00075.jpeg)]:这些事件可能以一定概率发生。
- [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-4fhvvaMK-1681366316295)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00076.jpeg)] 和 [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-z67ViLbq-1681366316295)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00077.jpeg)]:这是发生特定事件的概率。
- [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-nFXA9VfK-1681366316295)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00078.jpeg)]:鉴于 B 发生,这是 A 发生的概率。 这称为条件概率。
- [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-vUiLnnj1-1681366316296)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00079.jpeg)]:这是 A 和 B 一起出现的概率。
我们从基本假设开始,如下所示:
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-Fen8ppt5-1681366316296)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00080.jpeg)]
前面的方程将联合概率P(AB)
与条件概率P(A | B)
和边际概率P(B)
关联起来。 如果重写方程式,则条件概率表达式如下:
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-OeISoUzR-1681366316296)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00081.jpeg)]
这在某种程度上很直观-给定 B 的概率是通过将 A 和 B 都出现的概率除以 B 发生的概率而获得的。 这个想法是给定的,所以我们除以它的概率。 可以在这个链接中找到对此方程式的更严格的处理,标题为概率:联合概率,边际概率和条件概率。
类似地,通过对称我们有 [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-tucxlomG-1681366316297)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00082.jpeg)]。 因此,我们有 [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-jRKBpkGr-1681366316297)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00083.jpeg)]。 通过在两侧将表达式除以 [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-G2FZt6tZ-1681366316297)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00084.jpeg)] 并假设P(B) != 0
,我们得到:
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-PLJl4zKc-1681366316297)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00085.jpeg)]
前面的等式称为贝叶斯定理,这是所有贝叶斯统计推断的基础。 为了将贝叶斯定理与推论统计联系起来,我们将等式重铸为所谓的历时解释,如下所示:
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-XkCMxhKn-1681366316298)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00086.jpeg)]
其中,[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-Bp8CsX4c-1681366316298)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00087.jpeg)] 代表一个假设。
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-CMydvQ9c-1681366316298)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00088.jpeg)] 代表已经发生的事件,我们将其用于统计研究中,也称为数据。
然后,[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-CIQp96sn-1681366316299)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00089.jpeg)] 是我们的假设在观察数据之前的概率。 这称为先验概率。 贝叶斯统计学家经常认为使用先验概率是一种优势,因为先验知识或先前的结果可以用作当前模型的输入,从而提高了准确性。 有关此的更多信息,请参考这里。
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-uroeNp61-1681366316299)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00090.jpeg)] 是获得与假设无关的观测数据的概率。 这称为归一化常数。 并非总是需要计算归一化常数,尤其是在许多流行的算法(例如 MCMC)中,我们将在本章后面进行讨论。
给定我们观察到的数据,[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-CPWxYcj0-1681366316299)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00091.jpeg)] 是假设成立的概率。 这称为后验。
考虑到我们的假设,[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-09qXhzvP-1681366316299)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00092.jpeg)] 是获得数据的概率。 这称为可能性。
因此,贝叶斯统计等于应用贝叶斯规则来解决推论统计中的问题,其中 H 代表我们的假设,D 代表数据。
贝叶斯统计模型是根据参数进行转换的,这些参数的不确定性由概率分布表示。 这与将值视为确定性的频率论方法不同。 替代表示如下:
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-WeBzl3vX-1681366316300)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00093.jpeg)]
其中 [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-ssn92oDY-1681366316300)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00094.jpeg)] 是我们的未知数据,[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-JKT4HYIp-1681366316300)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00095.jpeg)] 是我们的观测数据
在贝叶斯统计中,我们对先验数据进行假设,并使用贝叶斯规则使用可能性来更新后验概率。 作为说明,让我们考虑以下问题。 这是通常称为问题的经典案例:
- 两个包含彩球
- 一个包含 50 个红色和 50 个蓝色球
- n 个包含 30 个红色和 70 个蓝色球
- 随机选择两个中的一个(概率为 50%),然后从两个中的一个随机抽出一个球
如果画了一个红色的球,那么它来自的概率是多少? 我们想要 [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-sIgbsoeq-1681366316301)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00091.jpeg)] 即 [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-MwHxveEA-1681366316301)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00096.jpeg)]。
在这里,[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-o2noUcaI-1681366316301)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00087.jpeg)] 表示从 Ur 中抽出球,[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-LqGbSIjd-1681366316301)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00088.jpeg)] 表示被抽出的球是红色:
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-eAphD7hJ-1681366316302)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00097.jpeg)]
我们知道 [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-b7hfIptj-1681366316302)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00098.jpeg)],[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-6JhpqNV9-1681366316302)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00099.jpeg)],[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-3UK7BRBR-1681366316302)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00100.jpeg)] 或 [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-E2OgeaFG-1681366316302)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00101.jpeg)]。
因此,我们得出结论 [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-PzXgofGu-1681366316303)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00102.jpeg)]。
贝叶斯理论与赔率
贝叶斯定理有时可以通过使用称为赔率的概率的替代表示形式,以更自然,更方便的形式表示。 赔率通常用比率表示,并且经常使用。 一匹马在比赛中获胜的赔率为 3 比 1(通常为 3:1),这表示该马有望以 75% 的概率获胜。
给定概率 p,可以将几率计算为几率= [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-M6XBdZT0-1681366316303)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00103.jpeg)],在 p = 0.75 的情况下,该几率变为 0.75:0.25,即 3:1。
我们可以通过使用赔率来重写贝叶斯定理的形式:
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-5IktkiQq-1681366316303)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00104.jpeg)]
贝叶斯统计的应用
贝叶斯统计可以应用于经典统计中遇到的许多问题,例如:
- 参数估计
- 预测
- 假设检验
- 线性回归
研究贝叶斯统计数据有许多令人信服的理由。 其中一些是利用先验信息更好地告知当前模型。 贝叶斯方法适用于概率分布而不是点估计,因此产生了更现实的预测。 贝叶斯推断是基于可用数据的假设 – P(假设|数据)
。 惯常方法试图根据假设拟合数据。 可以说,贝叶斯方法是更逻辑和经验的方法,因为它试图将信念建立在事实之上,而不是反过来。 有关此的更多信息,请参考这里。
概率分布
在本节中,我们将简要检查各种概率分布的属性。 这些分布中有许多用于贝叶斯分析。 因此,需要简要的提要。 我们还将说明如何使用matplotlib
生成和显示这些分布。 为了避免在每个部分的每个代码段重复导入语句,我将介绍以下标准的 Python 代码导入集,这些代码需要在以下命令中提到的任何代码段之前运行。 您只需在每个会话中运行一次这些导入。 导入如下:
In [1]: import pandas as pdimport numpy as npfrom matplotlib import pyplot as pltfrom matplotlib import colorsimport matplotlib.pyplot as pltimport matplotlib%matplotlib inline
拟合分布
贝叶斯分析中必须采取的步骤之一是使我们的数据适合概率分布。 选择正确的分发可能有点技巧,并且通常需要统计知识和经验,但是我们可以遵循一些准则来帮助我们。 这些如下:
- 确定数据是离散数据还是连续数据
- 检查数据的偏斜度/对称性,如果偏斜,请确定方向
- 确定下限和上限(如果有)
- 确定观察分布中极值的可能性
统计试验是具有一组明确定义的结果(称为样本空间)的可重复实验。 伯努利试验是“是/否”实验,如果是,则将随机 X 变量的值分配为 1,如果否,则将其分配为 0。抛硬币并查看其是否朝正面方向的事件是伯努利审判的例子。
有两类概率分布:离散分布和连续分布。 在以下各节中,我们将讨论这两类分布之间的差异,并浏览主要分布。
离散概率分布
在这种情况下,变量只能采用某些不同的值,例如整数。 离散随机变量的一个例子是当我们掷硬币 5 次时获得的正面数。 可能的值为{0,1,2,3,4,5}
。 例如,我们无法获得 3.82 个平面。 随机变量可以采用的值的范围由概率质量函数(pmf)指定。
离散均匀分布
离散均匀分布是一种对事件进行建模的分布,该事件具有一组有限的可能结果,其中每个结果均可能被观察到。 对于 [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-nk3l4W8T-1681366316303)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00105.jpeg)] 结果,每个结果都有 [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-VO0u87S9-1681366316303)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00106.jpeg)] 发生的可能性。
这方面的一个例子就是投掷均匀的色子。 六个结果中任何一个的概率为 [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-sl4eRiI0-1681366316304)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00107.jpeg)]。 PMF 由 [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-SlaeFnOg-1681366316304)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00106.jpeg)] 给出,期望值和方差分别由 [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-rcrEmend-1681366316304)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00108.jpeg)] 和 [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-N1gyRYve-1681366316304)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00109.jpeg)] 给出。
In [13]: from matplotlib import pyplot as pltimport matplotlib.pyplot as pltX=range(0,11)Y=[1/6.0 if x in range(1,7) else 0.0 for x in X]plt.plot(X,Y,'go-', linewidth=0, drawstyle='steps-pre', label="p(x)=1/6")plt.legend(loc="upper left")plt.vlines(range(1,7),0,max(Y), linestyle='-')plt.xlabel('x')plt.ylabel('p(x)')plt.ylim(0,0.5)plt.xlim(0,10)plt.title('Discrete uniform probability distribution with p=1/6')plt.show()
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-ewIPd1sm-1681366316305)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00110.jpeg)]
离散均匀分布
伯努利分布
伯努利分布用于衡量试验成功的可能性; 例如,硬币抛起头或尾的概率。 这可以用一个随机的 X 变量来表示,如果硬币正面朝上,则取值为 1,反之则为 0。 出现头部或尾部的概率分别由p
和q = 1-p
表示。
可以通过以下 pmf 表示:
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-lCsDzPXc-1681366316305)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00111.jpeg)]
期望值和方差由以下公式给出:
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-7iuUgGll-1681366316305)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00112.jpeg)]
此信息的参考位于这里。
现在,我们使用matplotlib
和scipy.stats
绘制伯努利分布,如下所示:
In [20]:import matplotlibfrom scipy.stats import bernoullia = np.arange(2)colors = matplotlib.rcParams['axes.color_cycle']plt.figure(figsize=(12,8))for i, p in enumerate([0.0, 0.2, 0.5, 0.75, 1.0]):ax = plt.subplot(1, 5, i+1)plt.bar(a, bernoulli.pmf(a, p), label=p, color=colors[i], alpha=0.5)ax.xaxis.set_ticks(a)plt.legend(loc=0)if i == 0:plt.ylabel("PDF at $k$")plt.suptitle("Bernoulli probability for various values of $p$")
Out[20]:
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-r4hnxONV-1681366316305)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00113.jpeg)]
二项分布
二项分布用于表示 n 个独立伯努利试验即 [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-ZAy3x098-1681366316305)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00114.jpeg)] 中的成功次数。
以掷硬币为例,此分布模拟了 [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-2s6RKJKx-1681366316306)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00105.jpeg)] 试验获得 X 正面攻击的机会。 对于 100 次抛掷,二项分布模拟了 0 头(极不可能)至 50 头(最高可能性)至 100 头(也极不可能)的可能性。 当赔率完全均匀时,这最终使二项分布对称,而当赔率远不那么均匀时,则使二项分布偏斜。 pmf 由以下表达式给出:
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-Y2V00NdI-1681366316306)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00115.jpeg)]
期望和方差分别由以下表达式给出:
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-loEKSbid-1681366316306)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00116.jpeg)]
In [5]: from scipy.stats import binomclrs = ['blue','green','red','cyan','magenta'] plt.figure(figsize=(12,6))k = np.arange(0, 22)for p, color in zip([0.001, 0.1, 0.3, 0.6, 0.999], clrs):rv = binom(20, p)plt.plot(k, rv.pmf(k), lw=2, color=color, label="$p$=" + str(round(p,1)))plt.legend()plt.title("Binomial distribution PMF")plt.tight_layout()plt.ylabel("PDF at $k$")plt.xlabel("$k$")
Out[5]:
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-Gz39cPPg-1681366316306)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00117.jpeg)]
二项分布
泊松分布
泊松分布对给定时间间隔内多个事件的概率进行建模,假设这些事件以已知的平均速率发生,并且连续事件的发生与自上一个事件以来的时间无关。
可以由泊松分布模型化的过程的一个具体示例是,如果一个人平均每天收到 23 封电子邮件。 如果我们假设电子邮件的到达时间是相互独立的,则个人每天收到的电子邮件总数可以通过泊松分布来建模。
另一个示例可能是每小时在特定车站停靠的火车数量。 泊松分布的 pmf 由以下表达式给出:
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-TOGbWSLr-1681366316307)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00118.jpeg)]
其中 [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-oDhphqGj-1681366316307)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00119.jpeg)] 是速率参数,代表每单位时间发生的预期事件/到达次数,[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-m2dYmkfm-1681366316307)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00120.jpeg)] 是随机变量,代表事件/到达次数。
期望值和方差分别由以下公式给出:
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-Ptmd5Wq5-1681366316307)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00121.jpeg)]
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-z1sf0Xaq-1681366316307)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00122.jpeg)]
有关更多信息,请参考这里。
使用matplotlib
将 pmf 绘制为各种值,如下所示:
In [11]: %matplotlib inlineimport numpy as npimport matplotlibimport matplotlib.pyplot as pltfrom scipy.stats import poissoncolors = matplotlib.rcParams['axes.color_cycle']k=np.arange(15)plt.figure(figsize=(12,8))for i, lambda_ in enumerate([1,2,4,6]):plt.plot(k, poisson.pmf(k, lambda_), '-o', label="$\\lambda$=" + str(lambda_), color=colors[i])plt.legend()plt.title("Possion distribution PMF for various $\\lambda$")plt.ylabel("PMF at $k$")plt.xlabel("$k$")plt.show()
Out[11]:
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-YP9WKx18-1681366316308)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00123.jpeg)]
鱼类分布
几何分布
对于独立的伯努利试验,几何分布测量获得成功所需的试验次数 X。 它还可以表示首次成功之前的失败次数 [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-hZTBGZiB-1681366316308)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00124.jpeg)]。
pmf 由以下表达式给出:
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-jIgguR4I-1681366316308)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00125.jpeg)]
自 [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-L8uQUe49-1681366316308)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00126.jpeg)] 以来,上述表达式是有意义的,并且如果需要 k 次尝试才能获得一个成功(p),则意味着我们必须有等于 [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-zGgYE75t-1681366316309)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00128.jpeg)] 的 [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-MY0QVo7R-1681366316309)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00127.jpeg)] 失败。
期望值和方差如下:
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-JEAbP2H5-1681366316309)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00129.jpeg)]
以下命令清楚地解释了前面的公式:
In [12]: from scipy.stats import geomp_vals=[0.01,0.2,0.5,0.8,0.9]x = np.arange(geom.ppf(0.01,p),geom.ppf(0.99,p))colors = matplotlib.rcParams['axes.color_cycle']for p,color in zip(p_vals,colors):x = np.arange(geom.ppf(0.01,p),geom.ppf(0.99,p))plt.plot(x,geom.pmf(x,p),'-o',ms=8,label='$p$=' + str(p))plt.legend(loc='best')plt.ylim(-0.5,1.5)plt.xlim(0,7.5)plt.ylabel("Pmf at $k$")plt.xlabel("$k$")plt.title("Geometric distribution PMF")
Out[12]:
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-W8uYMDLW-1681366316309)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00130.jpeg)]
几何分布
负二项分布
同样对于独立的伯努利试验,负二项分布测量在指定的成功次数 r 发生之前的尝试次数 [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-JK49SHci-1681366316309)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00131.jpeg)]。 一个例子是要获得 5 个硬币的掷硬币次数。 pmf 给出如下:
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-CifM7DOB-1681366316310)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00132.jpeg)]
期望和方差分别由以下表达式给出:
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-8jKiQR9n-1681366316310)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00133.jpeg)]
我们可以看到负二项是几何分布的概括,几何分布是负二项的特殊情况,其中 [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-md53lTOY-1681366316310)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00134.jpeg)]。
代码和图解如下所示:
In [189]: from scipy.stats import nbinomfrom matplotlib import colorsclrs = matplotlib.rcParams['axes.color_cycle']x = np.arange(0,11)n_vals = [0.1,1,3,6]p=0.5for n, clr in zip(n_vals, clrs):rv = nbinom(n,p)plt.plot(x,rv.pmf(x), label="$n$=" + str(n), color=clr)plt.legend()plt.title("Negative Binomial Distribution PMF")plt.ylabel("PMF at $x$")plt.xlabel("$x$")
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-K6GTesSk-1681366316310)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00135.jpeg)]
连续概率分布
在连续概率分布中,变量可以采用任何实数。 与离散概率分布一样,它不限于一组有限的值。 例如,一个健康的新生婴儿的平均体重大约在 6-9 磅之间。 例如,其重量可以是 7.3 磅。 连续概率分布的特征在于概率密度函数(PDF)。
随机变量可以假定的所有概率之和为 1。因此,概率密度函数图的面积为 1。
连续均匀分布
均匀分布对随机变量 X 建模,该随机变量 X 可以以相等的概率采用 [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-XyEWIpmq-1681366316311)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00136.jpeg)] 范围内的任何值。
PDF 由 [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-h0hNQWiz-1681366316311)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00137.jpeg)] 给出,否则由 [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-RKQxrvBt-1681366316311)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00138.jpeg)] 和 [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-0ts6nLNw-1681366316311)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00139.jpeg)] 给出。
期望和方差由以下表达式给出:
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-o52TCi3X-1681366316311)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00140.jpeg)]
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-iBwd6Tbn-1681366316312)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00141.jpeg)]
在下面的代码和图中,生成了连续均匀的概率分布,并针对各种样本大小绘制了该图:
In [11]: np.random.seed(100) # seed the random number generator# so plots are reproduciblesubplots = [111,211,311]ctr = 0 fig, ax = plt.subplots(len(subplots), figsize=(10,12))nsteps=10for i in range(0,3):cud = np.random.uniform(0,1,nsteps) # generate distribcount, bins, ignored = ax[ctr].hist(cud,15,normed=True)ax[ctr].plot(bins,np.ones_like(bins),linewidth=2, color='r')ax[ctr].set_title('sample size=%s' % nsteps)ctr += 1nsteps *= 100fig.subplots_adjust(hspace=0.4)plt.suptitle("Continuous Uniform probability distributions for various sample sizes" , fontsize=14)
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-o10gfo67-1681366316312)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00142.jpeg)]
指数分布
指数分布模拟了泊松过程中两个事件之间的等待时间。 泊松过程是遵循泊松分布的过程,在该过程中,事件以已知的平均速率不可预测地发生。 指数分布可以描述为几何分布的连续极限,并且也是马尔科夫的(无记忆)。
无记忆的随机变量表现出这样的特性,即其未来状态仅取决于有关当前时间的相关信息,而不取决于过去的信息。 建模马尔科夫/无记忆随机变量的示例是对短期股票价格行为及其遵循随机游走的思想进行建模。 这导致了所谓的金融有效市场假说。 有关更多信息,请参见这里。
指数分布的 PDF 由 [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-JrjmXLkS-1681366316312)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00143.jpeg)] = [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-4mFGPdH2-1681366316312)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00144.jpeg)] 给出。 期望和方差由以下表达式给出:
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-qjvAcJN6-1681366316312)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00145.jpeg)]
有关参考,请参考这个链接。
分布图和代码如下:
In [15]: import scipy.statsclrs = colors.cnamesx = np.linspace(0,4, 100)expo = scipy.stats.exponlambda_ = [0.5, 1, 2, 5]plt.figure(figsize=(12,4))for l,c in zip(lambda_,clrs):plt.plot(x, expo.pdf(x, scale=1./l), lw=2,color=c, label = "$\\lambda = %.1f$"%l)plt.legend()plt.ylabel("PDF at $x$")plt.xlabel("$x$")plt.title("Pdf of an Exponential random variable for various $\\lambda$");
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-1a8hEwtc-1681366316313)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00146.jpeg)]
正态分布
统计中最重要的分布可以说是正态/高斯分布。 它对中心值周围的概率分布进行建模,而没有左右偏倚。 遵循正态分布的现象有很多示例,例如:
- 婴儿的出生体重
- 测量误差
- 血压
- 考试分数
正态分布的重要性由中心极限定理强调,该极限定理指出,独立于同一分布而绘制的许多随机变量的均值近似正态,而与原始分布的形式无关。 其期望值和方差如下:
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-H0hWetnL-1681366316313)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00147.jpeg)]
正态分布的 PDF 由以下表达式给出:
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-svAvaMOW-1681366316313)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00148.jpeg)]
以下代码和图解说明了公式:
In [54]: import matplotlibfrom scipy.stats import normX = 2.5dx = 0.1R = np.arange(-X,X+dx,dx)L = list()sdL = (0.5,1,2,3)for sd in sdL:f = norm.pdfL.append([f(x,loc=0,scale=sd) for x in R])colors = matplotlib.rcParams['axes.color_cycle']for sd,c,P in zip(sdL,colors,L):plt.plot(R,P,zorder=1,lw=1.5,color=c,label="$\\sigma$=" + str(sd))plt.legend()ax = plt.axes()ax.set_xlim(-2.1,2.1)ax.set_ylim(0,1.0)plt.title("Normal distribution Pdf")plt.ylabel("PDF at $\\mu$=0, $\\sigma$")
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-TXVMeBL8-1681366316313)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00149.jpeg)]
可以在这个位置找到有关绘制分布图的 Python 代码的参考。
正态分布也可以视为二项分布的连续极限,而其他分布则可以视为 [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-lZ5uFYGk-1681366316314)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00150.jpeg)]。 我们可以在命令中看到二项分布,并绘制如下:
In [18]:from scipy.stats import binomfrom matplotlib import colorscols = colors.cnamesn_values = [1, 5,10, 30, 100]subplots = [111+100*x for x in range(0,len(n_values))]ctr = 0 fig, ax = plt.subplots(len(subplots), figsize=(6,12))k = np.arange(0, 200)p=0.5for n, color in zip(n_values, cols):k=np.arange(0,n+1)rv = binom(n, p)ax[ctr].plot(k, rv.pmf(k), lw=2, color=color)ax[ctr].set_title("$n$=" + str(n))ctr += 1fig.subplots_adjust(hspace=0.5)plt.suptitle("Binomial distribution PMF (p=0.5) for various values of n", fontsize=14)
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-jfbJubow-1681366316314)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00151.jpeg)]
随着 n 的增加,二项分布接近正态分布。 实际上,对于n >= 30
,可以在前面的图中清楚地看到。
贝叶斯统计与频率统计
在当今的统计中,关于如何解释数据和进行统计推断有两种思路。 迄今为止,经典且占主导地位的方法是所谓的“惯常方法”(请参阅第 7 章,“统计学之旅 -– 经典方法”),本章中的贝叶斯方法。
概率是多少?
贝叶斯和频繁主义世界观之间的辩论的核心是问题—我们如何定义概率?
在频率论世界观中,概率是从重复事件的发生频率得出的概念。 例如,当我们定义抛掷均匀硬币时获胜的概率等于一半。 这是因为当我们反复抛掷一个均匀的硬币时,当抛硬币的数量足够大时,正面数除以抛硬币的总数将接近 0.5。
贝叶斯世界观是不同的,概率的概念是它与一个人对事件发生的信念程度有关。 因此,对于贝叶斯统计学家而言,认为均匀色子概率为 5 的可能性为 [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-GfRicpqX-1681366316314)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00107.jpeg)],这与我们对该事件发生的可能性的信念有关。
如何定义模型
从模型定义的角度来看,常客通过使用重复的实验来分析数据和计算的度量标准如何变化,同时保持模型参数不变。 另一方面,贝叶斯算法利用固定的实验数据,但是改变了他们对模型参数的置信度,其解释如下:
- 频率论:如果模型固定,则数据会有所不同
- 贝叶斯:如果数据固定,则模型会有所不同
惯常方法使用所谓的最大似然方法来估计模型参数。 它涉及从一组独立且分布均匀的观测值生成数据,并将观测到的数据拟合到模型中。 最适合该数据的模型参数的值是最大似然估计器(MLE),有时可能是所观察数据的函数。
贝叶斯主义对这个问题的看法不同于概率框架。 概率分布用于描述值的不确定性。 贝叶斯实践者使用观察到的数据估计概率。 为了计算这些概率,它们使用单个估计量,即贝叶斯公式。 与“惯常”方法一样,这产生的是分布而不是点估计。
置信(频率)与可信(贝叶斯)区间
让我们比较一下 95% 的置信区间的含义,贝叶斯(Bayesian)练习者使用了 95% 可信区间的常用术语,可信区间为 95% 。
在频率论框架中,置信区间为 95% 意味着如果您无限次重复实验,并在过程中生成间隔,则其中 95% 的间隔将包含我们尝试估计的参数,通常称为 θ。 在这种情况下,时间间隔是随机变量,而不是参数估计值θ,该值在频率论世界观中是固定的。
在贝叶斯可信区间的情况下,我们的解释与归因于频繁置信区间的常规解释更加一致。 因此,我们有 [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-MCDxWWZU-1681366316315)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00152.jpeg)]。 在这种情况下,我们可以正确地得出结论:θ在此区间内的可能性为 95% 。
有关更多信息,请参阅频率论和贝叶斯主义:有什么大不了的?| SciPy 2014,Jake VanderPlas。
进行贝叶斯统计分析
进行贝叶斯统计分析涉及以下步骤:
- 指定概率模型:在这一步中,我们使用概率分布充分描述模型。 根据我们采集的样本的分布,我们尝试为其拟合模型,并尝试将概率分配给未知参数。
- 计算后验分布:后验分布是我们根据观察到的数据计算出的分布。 在这种情况下,我们将直接应用贝叶斯公式。 将根据上一步中指定的概率模型来指定它。
- 检查我们的模型:这是我们进行推断之前检查模型及其输出的必要步骤。 贝叶斯推断方法使用概率分布将概率分配给可能的结果。
似然函数和 PyMC 的蒙特卡洛估计
贝叶斯统计不仅是另一种方法。 它是实践统计的完全替代范例。 给定我们收集的数据,它使用概率模型进行推断。 这可以用基本表达式表示为P(H | D)
。
在这里,H 是我们的假设,即我们要证明的事物,D 是我们的数据或观察值。
在前面的讨论中提醒我们,贝叶斯定理的历时形式如下:
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-jOeqOFvq-1681366316315)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00086.jpeg)]
在这里,P(H)
是无条件的先验概率,表示我们进行试验之前所知道的。 假设我们的假设是正确的,P(D | H)
是我们的似然函数或获得我们观察到的数据的概率。
P(D)
是数据的概率,也称为归一化常数。 这可以通过对 H 上的分子进行积分来获得。
似然函数是贝叶斯计算中最重要的部分,它封装了有关数据中未知数的所有信息。 它与逆概率质量函数有些相似。
反对采用贝叶斯方法的一个论点是,先验的计算可能是主观的。 有很多观点支持这种方法。 其中之一是,如先前所述,可以包括外部先验信息。
似然值表示未知积分,在简单情况下可以通过解析积分获得。
对于涉及高维积分的更复杂的用例,需要蒙特卡洛(MC)积分,并且可以用于计算似然函数。
MC 集成可以通过多种采样方法来计算,例如统一采样,分层采样和重要性采样。 在蒙特卡洛积分中,我们可以将积分近似如下:
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-I4Ub9rQ0-1681366316315)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00153.jpeg)]
我们可以通过以下有限和来近似积分:
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-jKELMzT2-1681366316315)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00154.jpeg)]
其中,x 是来自 g 的样本向量。 可以根据大量定律并通过确保模拟误差较小来获得这一估计的好证据。
在 Python 中进行贝叶斯分析时,我们需要一个模块,该模块使我们能够使用前面所述的蒙特卡洛方法来计算似然函数。 PyMC
库满足了这一需求。 它提供了通常称为马尔可夫链蒙特卡洛(MCMC)的蒙特卡洛方法。 我不会进一步研究 MCMC 的技术细节,但是有兴趣的读者可以在以下参考文献的PyMC
中找到有关 MCMC 实现的更多信息:
- 贝叶斯估计中的蒙特卡洛积分
- 马尔科夫链蒙特卡洛最大似然
- 使用 Python 进行贝叶斯统计分析 – 第 1 部分 | SciPy 2014,Chris Fonnesbeck
MCMC 不是通用的灵丹妙药; 该方法存在一些缺点,其中之一是算法收敛缓慢。
贝叶斯分析示例 -– 切换点检测
在这里,我们将尝试使用贝叶斯推理并为一个有趣的数据集建模。 所涉及的数据集包含作者的 Facebook(FB)过去的历史记录。 我们已经清理了 FB 历史数据并将日期保存在fb_post_dates.txt
文件中。 文件中的数据如下所示:
head -2 ../fb_post_dates.txt
Tuesday, September 30, 2014 | 2:43am EDT
Tuesday, September 30, 2014 | 2:22am EDT
因此,我们看到一个日期时间序列,代表作者在 FB 上发布的日期和时间。 首先,我们将文件读入数据帧,将时间戳分为“日期”和“时间”列:
In [91]: filePath="./data/fb_post_dates.txt"fbdata_df=pd.read_csv(filePath, sep='|', parse_dates=[0], header=None,names=['Date','Time'])
接下来,我们检查数据如下:
In [92]: fbdata_df.head() #inspect the data
Out[92]: Date Time
0 2014-09-30 2:43am EDT
1 2014-09-30 2:22am EDT
2 2014-09-30 2:06am EDT
3 2014-09-30 1:07am EDT
4 2014-09-28 9:16pm EDT
现在,我们按日期对数据编制索引,创建一个DatetimeIndex
,以便可以对其进行重采样以按月进行计数,如下所示:
In [115]: fbdata_df_ind=fbdata_df.set_index('Date')fbdata_df_ind.head(5)
Out[115]: TimeDate2014-09-30 2:43am EDT2014-09-30 2:22am EDT2014-09-30 2:06am EDT2014-09-30 1:07am EDT2014-09-28 9:16pm EDT
我们显示有关索引的信息,如下所示:
In [116]: fbdata_df_ind.index
Out[116]: <class 'pandas.tseries.index.DatetimeIndex'>[2014-09-30, ..., 2007-04-16]Length: 7713, Freq: None, Timezone: None
现在,我们使用重新采样按月获取帖子数:
In [99]: fb_mth_count_=fbdata_df_ind.resample('M', how='count')fb_mth_count_.rename(columns={'Time':'Count'},inplace=True) # Rename fb_mth_count_.head()
Out[99]: CountDate2007-04-30 12007-05-31 02007-06-30 52007-07-31 502007-08-31 24
日期格式显示为每月的最后一天。 现在,我们创建一个 2007-2015 年 FB 帖子计数的散点图,并使点的大小与matplotlib
中的值成比例:
In [108]: %matplotlib inlineimport datetime as dt
#Obtain the count data from the DataFrame as a dictionaryyear_month_count=fb_bymth_count.to_dict()['Count'] size=len(year_month_count.keys())
#get dates as list of stringsxdates=[dt.datetime.strptime(str(yyyymm),'%Y%m') for yyyymm in year_month_count.keys()] counts=year_month_count.values()plt.scatter(xdates,counts,s=counts)plt.xlabel('Year')plt.ylabel('Number of Facebook posts')plt.show()
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-2RcyOns8-1681366316316)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00155.jpeg)]
我们要调查的问题是,在一段时间内,行为是否有所改变。 具体来说,我们希望确定是否在特定时期内更改了 FB 职位的平均人数。 在时间序列中,这通常称为切换点或更改点。
我们可以利用泊松分布对此建模。 您可能还记得,泊松分布可用于对时间序列计数数据进行建模。 (有关更多信息,请参考这里。)
如果我们用 [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-aOCdyIT0-1681366316316)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00156.jpeg)] 表示我们每月的 FB 职位数,则可以如下表示我们的模型:
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-Ia2MQEdr-1681366316316)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00157.jpeg)]
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-zxn7lVF7-1681366316317)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00158.jpeg)] 参数是泊松分布的rate
参数,但我们不知道它的值是多少。 如果我们查看 FB 时间序列计数数据的散点图,我们可以看到在 2010 年中后期左右某个时候的帖子数量有所增加,这可能与 2010 年南非世界杯的开始恰好相吻合。 作者参加了。
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-kuKyUBSv-1681366316317)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00159.jpeg)] 参数是切换点,即速率参数更改时的时间,而 [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-iLjLr6q9-1681366316317)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00160.jpeg)] 和 [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-Jvg8TDM4-1681366316317)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00161.jpeg)] 分别是切换点之前和之后的 [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-Xu5z8w38-1681366316317)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00158.jpeg)] 参数的值。 可以表示如下:
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-zp8BVAyS-1681366316318)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00162.jpeg)]
注意,上面 [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-MaFSfkza-1681366316318)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00163.jpeg)] 中指定的变量都是贝叶斯随机变量。 对于代表人们对其值的信念的贝叶斯随机变量,我们需要使用概率分布对它们进行建模。 我们想推断 [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-xjs5YUki-1681366316318)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00160.jpeg)] 和 [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-lG7nldMp-1681366316318)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00161.jpeg)] 的值,这些值是未知的。 在PyMC
中,我们可以使用随机和确定性类表示随机变量。 我们注意到,指数分布是泊松事件之间的时间量。 因此,在 [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-m4mjoHuZ-1681366316318)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00160.jpeg)] 和 [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-9fghGjdX-1681366316319)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00161.jpeg)] 的情况下,我们选择指数分布来对其建模,因为它们可以是任何正数:
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-BSMZGVXb-1681366316319)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00164.jpeg)]
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-tR3g6oTr-1681366316319)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00165.jpeg)]
在 [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-ldV96tmn-1681366316319)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00159.jpeg)] 的情况下,我们将选择使用均匀分布对其进行建模,这反映了我们的信念,即切换点很有可能在整个时间段的任何一天发生。 因此,我们有:
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-8xrk1skQ-1681366316320)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00166.jpeg)]
在此,[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-W8BSJhUT-1681366316320)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00167.jpeg)],[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-GZP6Rk38-1681366316320)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00168.jpeg)] 对应于年份 [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-4BVM6RWt-1681366316321)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00169.jpeg)] 的上下边界。 现在让我们使用PyMC
表示我们先前开发的模型。 现在,我们将使用PyMC
来查看是否可以在 FB 发布数据中检测到切换点。 除了散点图,我们还可以在条形图中显示数据。 为了做到这一点,首先,我们需要获取列表中按月排序的 FB 帖子数:
In [69]: fb_activity_data = [year_month_count[k] for k in sorted(year_month_count.keys())]fb_activity_data[:5]Out[70]: [1, 0, 5, 50, 24]In [71]: fb_post_count=len(fb_activity_data)
我们使用matplotlib
渲染条形图:
In [72]: from IPython.core.pylabtools import figsizeimport matplotlib.pyplot as pltfigsize(8, 5)plt.bar(np.arange(fb_post_count),fb_activity_data, color="#49a178")plt.xlabel("Time (months)")plt.ylabel("Number of FB posts")plt.title("Monthly Facebook posts over time")plt.xlim(0,fb_post_count);
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-XtrWT4OE-1681366316321)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00170.jpeg)]
查看前面的条形图,可以得出结论,FB 频率发布行为是否随时间变化? 我们可以在开发的模型上使用PyMC
来帮助我们找出更改,如下所示:
In [88]: # Define data and stochasticsimport pymc as pmswitchpoint = pm.DiscreteUniform('switchpoint',lower=0,upper=len(fb_activity_data)-1,doc='Switchpoint[month]')avg = np.mean(fb_activity_data)early_mean = pm.Exponential('early_mean', beta=1./avg)late_mean = pm.Exponential('late_mean', beta=1./avg)late_mean
Out[88]:<pymc.distributions.Exponential 'late_mean' at 0x10ee56d50>
在这里,我们为速率参数r
定义了一种方法,并使用泊松分布对计数数据进行建模,如前所述:
In [89]: @pm.deterministic(plot=False)def rate(s=switchpoint, e=early_mean, l=late_mean):''' Concatenate Poisson means '''out = np.zeros(len(fb_activity_data))out[:s] = eout[s:] = lreturn outfb_activity = pm.Poisson('fb_activity', mu=rate, value=fb_activity_data, observed=True)fb_activity
Out[89]: <pymc.distributions.Poisson 'fb_activity' at 0x10ed1ee50>
在前面的代码片段中,@pm.deterministic
是一个修饰符,它表示速率函数是确定性的,这意味着其值完全由其他变量(在这种情况下为e
,s
和l
)确定。 为了告诉PyMC
将比率函数转换为确定性对象,装饰器是必需的。 如果不指定装饰器,则会发生错误。 (有关更多信息,请参阅这里,以获取有关 Python 装饰器的信息。)
有关更多信息,请参阅以下网页:
- http://en.wikipedia.org/wiki/Poisson_process
- http://pymc-devs.github.io/pymc/tutorial.html
- https://github.com/CamDavidsonPilon/Probabilistic-Programming-and-Bayesian-Methods-for-Hackers
现在,我们使用 FB Count 数据(fb_activity
)和 [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-8SokGmpZ-1681366316321)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00171.jpeg)](分别为early_mean
,late_mean
和rate
)参数创建一个模型。
接下来,使用Pymc
创建一个 MCMC 对象,该对象使我们能够使用马尔科夫链蒙特卡洛方法拟合数据。 然后,我们在所得的 MCMC 对象上调用样本以进行拟合:
In [94]: fb_activity_model=pm.Model([fb_activity,early_mean,late_mean,rate])
In [95]: from pymc import MCMCfbM=MCMC(fb_activity_model)
In [96]: fbM.sample(iter=40000,burn=1000, thin=20)[-----------------100%-----------------] 40000 of 40000 complete in 11.0 sec
使用 MCMC 拟合模型涉及使用马尔科夫链蒙特卡洛方法生成后验P(s, e, l | D)
的概率分布。 它使用蒙特卡洛过程反复模拟数据采样,并执行此操作,直到基于多个条件该算法似乎收敛到稳态为止。 这是一个马尔可夫过程,因为连续的样本仅取决于先前的样本。 (有关马尔可夫链收敛的更多参考,请参考这里。)
生成的样本称为迹线。 我们可以通过查看迹线的直方图来查看参数的边缘后验分布:
In [97]: from pylab import hist,show%matplotlib inlinehist(fbM.trace('late_mean')[:])
Out[97]: (array([ 15., 61., 214., 421., 517., 426., 202.,70., 21., 3.]),array([ 102.29451192, 103.25158404, 104.20865616,105.16572829, 106.12280041, 107.07987253,108.03694465, 108.99401677, 109.95108889,110.90816101, 111.86523313]),<a list of 10 Patch objects>)
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-arwuyeok-1681366316321)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00172.jpeg)]
In [98]:plt.hist(fbM.trace('early_mean')[:])
Out[98]: (array([ 20., 105., 330., 489., 470., 314., 147.,60., 3., 12.]),array([ 49.19781192, 50.07760882, 50.95740571,51.83720261, 52.71699951, 53.59679641,54.47659331, 55.35639021, 56.2361871 ,57.115984 , 57.9957809 ]),<a list of 10 Patch objects>)
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-zbpt99QM-1681366316322)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00173.jpeg)]
在这里,我们可以看到Switchpoint
的月数形式:
In [99]: fbM.trace('switchpoint')[:]
Out[99]: array([38, 38, 38, ..., 35, 35, 35])In [150]: plt.hist(fbM.trace('switchpoint')[:])
Out[150]: (array([ 1899., 0., 0., 0., 0., 0.,0., 0., 0., 51.]),array([ 35\\. , 35.3, 35.6, 35.9, 36.2, 36.5, 36.8,37.1, 37.4, 37.7, 38\\. ]),<a list of 10 Patch objects>)
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-iLiUSWsR-1681366316322)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00174.jpeg)]
历史切换点
我们可以看到Switchpoint
在第 35-38 个月左右。 在这里,我们使用matplotlib
在单个图中显示e
,s
和l
的边际后验分布:
In [141]: early_mean_samples=fbM.trace('early_mean')[:]late_mean_samples=fbM.trace('late_mean')[:]switchpoint_samples=fbM.trace('switchpoint')[:]
In [142]: from IPython.core.pylabtools import figsizefigsize(12.5, 10)# histogram of the samples:fig = plt.figure()fig.subplots_adjust(bottom=-0.05)n_mths=len(fb_activity_data)ax = plt.subplot(311)ax.set_autoscaley_on(False)plt.hist(early_mean_samples, histtype='stepfilled',bins=30, alpha=0.85, label="posterior of $e$",color="turquoise", normed=True)plt.legend(loc="upper left")plt.title(r"""Posterior distributions of the variables$e, l, s$""",fontsize=16)plt.xlim([40, 120])plt.ylim([0, 0.6])plt.xlabel("$e$ value",fontsize=14)ax = plt.subplot(312)ax.set_autoscaley_on(False)plt.hist(late_mean_samples, histtype='stepfilled',bins=30, alpha=0.85, label="posterior of $l$",color="purple", normed=True)plt.legend(loc="upper left")plt.xlim([40, 120])plt.ylim([0, 0.6])plt.xlabel("$l$ value",fontsize=14)plt.subplot(313)w = 1.0 / switchpoint_samples.shape[0] *np.ones_like(switchpoint_samples)plt.hist(switchpoint_samples, bins=range(0,n_mths), alpha=1,label=r"posterior of $s$", color="green",weights=w, rwidth=2.)plt.xlim([20, n_mths - 20])plt.xlabel(r"$s$ (in days)",fontsize=14)plt.ylabel("probability")plt.legend(loc="upper left")plt.show()
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-Kga6K1zB-1681366316322)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00175.jpeg)]
边缘后验分布
PyMC 还具有绘图功能。 (它使用matplotlib
。)在下面的图中,我们显示了时间序列图,自相关图(acorr)以及为早期均值,晚期均值和切换点绘制的样本的直方图。 直方图可用于可视化后验分布。 自相关图显示上一周期的值是否与当前周期的值密切相关。
In [100]: from pymc.Matplot import plotplot(fbM)Plotting late_meanPlotting switchpointPlotting early_mean
以下是最新的均值图:
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-Jem4E82R-1681366316322)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00176.jpeg)]
pymc_comprehensive_late_mean
在这里,我们显示切换点图:
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-vLHhUYiz-1681366316323)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00177.jpeg)]
Pymc 综合切换点
在这里,我们显示早期均值图:
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-EpJHhwXL-1681366316323)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00178.jpeg)]
pymc 综合早期均值
从 PyMC 的输出中,我们可以得出结论,转换点距离时间序列开始大约 35-38 个月。 这对应于 2010 年 3 月至 7 月左右的某个时间。作者可以证明这是他使用 FB 的标志性一年,因为那是在南非举行的足球(足球)世界杯决赛的年份, 他参加了。
参考资料
要更深入地了解我们涉及的贝叶斯统计主题,请查看以下参考资料:
- 《黑客的概率编程和贝叶斯方法》
- 《贝叶斯数据分析,第三版》,Andrew Gelman
- 《贝叶斯选择》(理论比较多), Christian P Robert
- PyMC 文档
总结
在本章中,我们对过去几年中统计和数据分析中最热门的趋势之一-贝叶斯统计推断方法进行了一次旋风式浏览。 我们在这里覆盖了很多领域。
我们研究了贝叶斯统计方法所需要的内容,并讨论了为什么贝叶斯观点是一种引人注目的观点的各种因素-事实胜于信念。 我们解释了关键的统计分布,并展示了如何使用各种统计包在matplotlib
中生成和绘制它们。
我们解决了一个相当困难的话题,但又没有过多简化,并演示了如何使用 PyMC 包和蒙特卡洛模拟方法来展示贝叶斯统计数据的能力,以建立模型,进行趋势分析以及对真实数据集进行推断(Facebook 用户帖子)。 在下一章中,我们将讨论 Pandas 库的架构。
九、Pandas 库架构
在本章中,我们研究了可供 Pandas 用户使用的各种库。 本章旨在作为简短指南,帮助用户围绕 Pandas 提供的各种模块和库进行导航和查找。 它提供了有关库代码组织方式的分解,还提供了有关各个模块的简要说明。 对于感兴趣的用户来说,这对他们下面的 Pandas 的内部工作以及希望为代码库做出贡献的用户来说,将是最有价值的。 我们还将简要演示如何使用 Python 扩展来提高性能。 将讨论的各种主题如下:
- Pandas 库层次结构简介
- Pandas 模块和文件的描述
- 使用 Python 扩展来提高性能
Pandas 文件层次结构简介
通常,在安装时,会将 Pandas 作为 Python 模块安装在第三方 Python 模块的标准位置:
平台 | 标准安装位置 | 示例 |
---|---|---|
Unix/MacOS | prefix/lib/pythonX.Y/site-packages |
/usr/local/lib/Python 2.7/site-packages |
Windows | prefix\\Lib\\site-packages |
C:\\Python27\\Lib\\site-packages |
安装的文件遵循特定的层次结构:
pandas/core
:此文件包含用于基本数据结构的文件,例如序列/数据帧和相关函数。pandas/src
:包含用于实现基本算法的 Cython 和 C 代码。pandas/io
:它包含输入/输出工具(例如平面文件,Excel,HDF5,SQL 等)。pandas/tools
:它包含辅助数据算法合并和连接例程,连接,数据透视表等。pandas/sparse
:它包含序列,数据帧,面板等的稀疏版本。pandas/stats
:包含线性回归和面板回归,以及移动窗口回归。 这应该由 statsmodels 中的功能代替。pandas/util
:包含工具,开发和测试工具。pandas/rpy
:包含用于连接到 R 的 RPy2 接口。
注意
作为参考,请参见这里。
Pandas 模块和文件的描述
在本节中,我们将简要介绍组成 Pandas 库的各个子模块和文件。
Pandas /core
该模块包含 Pandas 的核心子模块。 讨论如下:
-
api.py
:这将导入一些关键模块供以后使用。 -
array.py
:这可以隔离 Pandas 对 numPy 的暴露,即所有直接的 numPy 使用。 -
base.py
:这定义了基本类别,例如StringMixin
,PandasObject
,它们是各种 Pandas 对象(例如Period
,PandasSQLTable
,sparse.array.SparseArray/SparseList
,internals.Block
,internals.BlockManager
,generic.NDFrame
,groupby.GroupBy
,base.FrozenList
,base.FrozenNDArray
,io.sql.PandasSQL
,io.sql.PandasSQLTable
,tseries.period.Period
,FrozenList
,FrozenNDArray
:IndexOpsMixin
和DatetimeIndexOpsMixin
。 -
common.py
:这定义了用于处理数据结构的通用工具方法。 例如,isnull
对象检测到缺少的值。 -
config.py
:这是用于处理包范围内的可配置对象的模块。 它定义了以下类别:OptionError
,DictWrapper
,CallableDynamicDoc
,option_context
和config_init
。 -
datetools.py
:这是处理 Python 中日期的函数的集合。 -
frame.py
:这定义了 Pandas 的数据帧类及其各种方法。数据帧继承自NDFrame
。 (见下文)。 -
generic.py
:这定义了通用NDFrame
基类,它是 Pandas 的数据帧,序列和面板类的基类。NDFrame
源自base.py
中定义的 Pandas 对象。NDFrame
可以看作是 Pandas 数据帧的 N 维版本。 有关此的更多信息,请访问这里。 -
categorical.py
:这定义了分类,这是一个从 Pandas 对象派生的类,它表示分类变量 la R/S-plus。 (我们稍后会扩展您的知识)。 -
format.py
:这定义了整个格式化程序类,例如CategoricalFormatter
,SeriesFormatter
,TableFormatter
,DataFrameFormatter
,HTMLFormatter
,CSVFormatter
,ExcelCell
,ExcelFormatter
,GenericArrayFormatter
,FloatArrayFormatter
,IntArrayFormatter
,Datetime64Formatter
,Timedelta64Formatter
和EngFormatter
。 -
groupby.py
:这定义了启用groupby
函数的各种类。 讨论如下:Splitter classes
:这包括DataSplitter
,ArraySplitter
,SeriesSplitter
,FrameSplitter
和NDFrameSplitter
Grouper/Grouping classes
:这包括Grouper
,GroupBy
,BaseGrouper
,BinGrouper
,Grouping
,SeriesGroupBy
,NDFrameGroupBy
-
ops.py
: 这定义了一个内部 API,用于对 Pandas 对象进行算术运算。 它定义了向对象添加算术方法的函数。 它定义了一个_create_methods
元方法,该方法用于使用算术,比较和布尔方法构造器创建其他方法。add_methods
方法采用一些新方法,将它们添加到现有方法列表中,并将其绑定到其相应的类。add_special_arithmetic_methods
和add_flex_arithmetic_methods
方法调用_create_methods
和add_methods
将算术方法添加到类中。它还定义了
_TimeOp
类,该类是与日期时间相关的算术运算的包装。 它包含Wrapper
函数,用于对序列,数据帧和面板函数进行算术,比较和布尔运算-_arith_method_SERIES(..)
,_comp_method_SERIES(..)
,_bool_method_SERIES(..)
,_flex_method_SERIES(..)
,_arith_method_FRAME(..)
,_comp_method_FRAME(..)
,_flex_comp_method_FRAME(..)
)_arith_method_PANEL(..)
,_comp_method_PANEL(..)
。 -
index.py
:这定义了Index
类及其相关函数。 所有 Pandas 的对象(序列,数据帧和面板)都使用索引来存储轴标签。 它的下面是一个不可变的数组,提供了可以切片的有序集合。 -
internals.py
:这定义了多个对象类。 这些列出如下:Block
:这是具有 Pandas 其他功能的同质 N 维numpy.ndarray
对象。 例如,它使用__slots__
将对象的属性限制为ndim
,values
和_mgr_locs
。 它充当其他Block
子类的基类。NumericBlock
:这是数字类型的Blocks
的基类。FloatOrComplexBlock
:这是从NumericBlock
继承的FloatBlock
和ComplexBlock
的基类ComplexBlock
:这是处理具有复杂类型的Block
对象的类。FloatBlock
:这是处理浮点型Block
对象的类。IntBlock
:这是处理具有整数类型的Block
对象的类。TimeDeltaBlock
,BoolBlock
和DatetimeBlock
:这些是timedelta
,Boolean
和datetime
的Block
类。ObjectBlock
:这是为用户定义的对象处理Block
对象的类。SparseBlock
:这是处理相同类型的稀疏数组的类。BlockManager
:这是管理一组Block
对象的类。 它不是公共 API 类。SingleBlockManager
:这是管理一个Block
的类。JoinUnit
:这是Block
对象的工具类。
-
matrix.py
:这将数据帧导入为DataMatrix
。 -
nanops.py
:这些是用于处理 NaN 值的类和函数。 -
ops.py
:这定义了 Pandas 对象的算术运算。 它不是公共 API。 -
panel.py
,panel4d.py
和panelnd.py
:这些提供了 Pandas 的面板对象的功能。 -
series.py
:它定义序列从NDFrame
和IndexOpsMixin
继承的 pandasSeries
类及其各种方法。 -
sparse.py
:这定义用于处理稀疏数据结构的导入。 稀疏数据结构被压缩,从而省略了与 NaN 匹配或缺少值的数据点。 有关此的更多信息,请访问这里。 -
strings.py
: 它们具有处理字符串的各种函数。[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-AgXnNxtQ-1681366316323)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00179.jpeg)]
Pandas /io
该模块包含用于数据 I/O 的各种模块。 讨论如下:
api.py
:这定义了数据 I/O API 的各种导入。auth.py
:这定义了处理身份验证的各种方法。common.py
:这定义了 I/O API 的通用函数。data.py
:这定义用于处理数据的类和方法。DataReader
方法从各种在线来源(例如 Yahoo 和 Google)读取数据。date_converters.py
:定义日期转换函数。excel.py
:此模块解析和转换 Excel 数据。 这定义了ExcelFile
和ExcelWriter
类。ga.py
:这是 Google Analytics(分析)功能的模块。gbq.py
:这是 Google BigQuery 的模块。html.py
:这是用于处理 HTML I/O 的模块。json.py
:这是用于处理 Pandas 中 json I/O 的模块。 这定义了Writer
,SeriesWriter
,FrameWriter
,Parser
,SeriesParser
和FrameParser
类。packer.py
:这是 msgpack 序列化程序支持,用于将 Pandas 数据结构读取和写入磁盘。parsers.py
:此模块定义了各种函数和类,这些函数和类用于解析和处理文件以创建 Pandas 的数据帧。 下文讨论的所有三个read_*
函数都具有多个可配置的读取选项。 有关更多详细信息,请参见此参考。read_csv(..)
:这定义了pandas.read_csv()
函数,可用于将 CSV 文件的内容读取到数据帧中。read_table(..)
:这会将制表符分隔的表文件读取到数据帧中。read_fwf(..)
:这会将固定宽度格式的文件读取到数据帧中。TextFileReader
:这是用于读取文本文件的类。ParserBase
:这是解析器对象的基类。CParserWrapper
和PythonParser
:这些分别是 C 和 Python 的解析器。 它们都继承自ParserBase
。FixedWidthReader
:这是读取固定宽度数据的类。 定宽数据文件包含文件中特定位置的字段。FixedWithFieldParser
:这是用于解析从PythonParser
继承的固定宽度字段的类。
pickle.py
:这提供了腌制(序列化)Pandas 对象的方法。 讨论如下:to_pickle(..)
:这会将对象序列化为文件。read_pickle(..)
:这将从文件中读取序列化的对象到 pandas 对象。 仅应与受信任的来源一起使用。
pytables.py
:这是PyTables
模块的接口,用于将 Pandas 数据结构读取和写入磁盘上的文件。sql.py
:它是类和函数的集合,用于使能够从试图与数据库无关的关系数据库中检索数据。 讨论如下:PandasSQL
:这是将 Pandas 与 SQL 连接的基类。 它提供了必须由子类实现的伪read_sql
和to_sql
方法。PandasSQLAlchemy
:这是PandasSQL
的子类,它可以使用SQLAlchemy
在数据帧和 SQL 数据库之间进行转换。PandasSQLTable
类:它将 Pandas 表(数据帧)映射到 SQL 表。pandasSQL_builder(..)
:这将根据提供的参数返回正确的PandasSQL
子类。PandasSQLTableLegacy
类:这是PandasSQLTable
的旧支持版本。PandasSQLLegacy
类:这是PandasSQLTable
的旧支持版本。get_schema(..)
:这将获取给定框架的 SQL 数据库表架构。read_sql_table(..)
:这将 SQL DB 表读入数据帧。read_sql_query(..)
:这会将 SQL 查询读取到数据帧中。read_sql(..)
:这将 SQL 查询/表读入数据帧。to_sql(..)
:此操作将数据帧中存储的记录写入 SQL 数据库。
stata.py
:此工具包含用于将Stata
文件处理为 Pandas 数据帧的工具。wb.py
:这是用于从世界银行网站下载数据的模块。
Pandas tools
-
util.py
:具有定义的其他util
函数,例如match(..)
,cartesian_product(..)
和compose(..)
。 -
tile.py
:具有一组函数,这些函数可以量化输入数据,从而实现tile
函数。 除了cut(..)
和qcut(..)
之外,大多数函数都是内部函数。 -
rplot.py
:这是一个模块,提供在 Pandas 中生成网格图的函数。 -
plotting.py
: 这提供了一组以序列或数据帧为参数的绘图函数。-
scatter_matrix(..)
:这将绘制散点图矩阵 -
andrews_curves(..)
:此图将多元数据绘制为曲线,这些曲线是使用样本作为傅立叶级数的系数创建的 -
parallel_coordinates(..)
:这是一种绘图技术,可让您查看数据中的聚类并直观地估计统计信息 -
lag_plot(..)
:用于检查数据集或时间序列是否随机 -
autocorrelation_plot(..)
:用于检查时间序列中的随机性 -
bootstrap_plot(..)
:此图用于以视觉方式确定统计量度的不确定性,例如平均值或中位数 -
radviz(..)
: 该图用于可视化多元数据提示
以上信息的参考来自这里
-
-
pivot.py
: 此函数用于处理 Pandas 中的数据透视表。 它是主要函数pandas.tools.pivot_table(..)
,它创建类似于电子表格的数据透视表作为数据帧提示
以上信息的参考来自这里
-
merge.py
:提供合并序列,数据帧和面板对象(例如merge(..)
和concat(..)
)的函数 -
describe.py
:这提供了一个value_range(..)
函数,该函数以序列的形式返回数据帧的最大值和最小值。
Pandas /sparse
这是提供序列,数据帧和面板的稀疏实现的模块。 所谓稀疏,是指省略或丢失诸如 0 之类的值的数组。
有关此的更多信息,请访问这里 。
api.py
:这是一组便利导入array.py
:它是SparseArray
数据结构的实现frame.py
:它是SparseDataFrame
数据结构的实现list.py
:它是SparseList
数据结构的实现panel.py
:它是SparsePanel
数据结构的实现series.py
:它是SparseSeries
数据结构的实现
Pandas stats
api.py
:这是一组便利导入。common.py
:定义模块中其他函数调用的内部函数。fama_macbeth.py
:包含 Fama-Macbeth 回归的类定义和函数。 有关 FM 回归的更多信息,请访问这里。interface.py
:它定义ols(..)
,它返回普通最小二乘(OLS)回归对象。 它从pandas.stats.ols
模块导入。math.py
:具有有用的函数,定义如下:rank(..)
,solve(..)
和inv(..)
:它们分别用于矩阵秩,解和逆is_psd(..)
:这检查矩阵的正定性newey_west(..)
:用于协方差矩阵计算calc_F(..)
:这将计算 F 统计信息
misc.py
:用于其他函数。moments.py
:这提供了滚动和扩展的统计量度,包括在 Cython 中实现的时刻。 这些方法包括:rolling_count(..)
,rolling_cov(..)
,rolling_corr(..)
,rolling_corr_pairwise(..)
,rolling_quantile(..)
,rolling_apply(..)
,rolling_window(..)
,expanding_count(..)
,expanding_quantile(..)
,expanding_cov(..)
,expanding_corr(..)
,expanding_corr_pairwise(..)
,expanding_apply(..)
,ewma(..)
,ewmvar(..)
,ewmstd(..)
,ewmcov(..)
和ewmcorr(..)
。ols.py
:这实现 OLS 并提供 OLS 和MovingOLS
类。 OLS 运行完整的样本最小二乘回归,而MovingOLS
生成滚动或扩展的简单 OLS。plm.py
:这为面板数据提供了线性回归对象。 这些类的讨论如下:PanelOLS
:这是面板对象的 OLSMovingPanelOLS
:这是面板对象的滚动/扩展 OLSNonPooledPanelOLS
:-这是面板对象的非池化 OLS
var.py
: 这提供了向量自回归类,讨论如下:-
VAR
:这是对序列和数据帧中的多元数据的向量自动回归 -
PanelVAR
: 这是面板对象中多元数据的向量自动回归提示
有关向量自回归的更多信息,请访问这里
-
Pandas /util
-
testing.py
:这提供了断言,调试,单元测试以及其他用于测试的类/函数。 它包含许多特殊的断言函数,这些函数使检查序列,数据帧或面板对象是否等效的操作更加容易。 其中一些函数包括assert_equal(..)
,assert_series_equal(..)
,assert_frame_equal(..)
和assert_panelnd_equal(..)
。pandas.util.testing
模块对 Pandas 代码库的贡献者特别有用。 它定义了一个util.TestCase
类。 它还提供了用于处理语言环境,控制台调试,文件清除,比较器等的工具,以供潜在的代码库贡献者进行测试。 -
terminal.py
:此函数主要是内部函数,与获取有关终端的某些特定详细信息有关。 唯一的暴露函数是get_terminal_size()
。 -
print_versions.py
:它定义了get_sys_info()
函数,该函数返回系统信息的字典,以及show_versions(..)
函数,显示可用的 Python 库的版本。 -
misc.py
:这定义了几个其他工具。 -
decorators.py
: 这定义了一些装饰器函数和类。提示
Substitution
和Appender
类是装饰器,它们在函数docstrings
上执行替换和附加操作,有关 Python 装饰器的更多信息,请访问这里。 -
clipboard.py
:这包含跨平台剪贴板方法,可通过键盘启用复制和粘贴功能。 Pandas I/O API 包括pandas.read_clipboard()
和pandas.to_clipboard(..)
之类的函数。
Pandas /rpy
如果模块已安装在计算机中,则该模块将尝试提供 R 统计包的接口。 在版本 0.16.0 和更高版本中不推荐使用。 可以使用rpy2
模块代替它的功能。
base.py
:这为 R 中著名的lm
函数定义了一个类common.py
:这提供了许多函数,可以将 Pandas 对象转换为等效的 R 版本mass.py
:这是rlm
– R 的lm
函数的未实现版本var.py
:其中包含未实现的类VAR
Pandas /test
这个模块为 Pandas 中的各种物体提供了许多测试。 特定库文件的名称是不言自明的,除了邀请读者进行探索之外,在这里我将不进一步详细介绍。
Pandas /compat
与兼容性相关的功能解释如下:
chainmap.py
和chainmap_impl.py
:这提供了一个ChainMap
类,可以将多个dicts
或映射分组,以便生成可以更新的单个视图pickle_compat.py
:此函数提供了在 0.12 之前的版本中腌制 Pandas 对象的功能。openpyxl_compat.py
:这检查openpyxl
的兼容性
Pandas /calc
这是提供计算功能的模块,讨论如下:
api.py
:这包含eval
和expr
的导入。align.py
:这实现了数据对齐的函数。common.py
:它包含几个内部函数。engines.py
:这定义了抽象引擎NumExprEngine
和PythonEngine
。PythonEngine
求值表达式,主要用于测试目的。eval.py
:这定义了最重要的eval(..)
函数以及其他一些重要函数。expressions.py
:通过numexpr
提供快速的表达式求值。numexpr
函数用于加速某些数值运算。 它使用多个内核以及智能分块和缓存加速。 它定义了evaluate(..)
和where(..)
方法。ops.py
:这定义了eval
使用的运算符类。 它们是Term
,Constant
,Op
,BinOp
,Div
和UnaryOp
。pytables.py
:这为PyTables
查询提供了查询接口。scope.py
:这是示波器操作的模块。 它定义了一个Scope
类,它是一个保存范围的对象。
提示
有关numexpr
的更多信息,请访问这里。 有关此模块用法的信息,请访问这里 。
Pandas /series
-
api.py
:这是一组便捷导入 -
converter.py
: 这定义了一组用于格式化和转换与日期时间相关的对象的类。 导入后,pandas 向matplotlib
注册了一组单位转换器。-
这通过
register()
函数完成,说明如下:In [1]: import matplotlib.units as munits In [2]: munits.registry Out[2]: {}In [3]: import pandas In [4]: munits.registry Out[4]: {pandas.tslib.Timestamp: <pandas.tseries.converter.DatetimeConverter instance at 0x7fbbc4db17e8>,pandas.tseries.period.Period: <pandas.tseries.converter.PeriodConverter instance at 0x7fbbc4dc25f0>,datetime.date: <pandas.tseries.converter.DatetimeConverter instance at 0x7fbbc4dc2fc8>,datetime.datetime: <pandas.tseries.converter.DatetimeConverter instance at 0x7fbbc4dc2a70>,datetime.time: <pandas.tseries.converter.TimeConverter instance at 0x7fbbc4d61e18>}
-
Converter
:此类包括TimeConverter
,PeriodConverter
和DateTimeConverter
-
Formatters
:此类包括TimeFormatter
,PandasAutoDateFormatter
和TimeSeries_DateFormatter
-
Locators
:此类包括PandasAutoDateLocator
,MilliSecondLocator
和TimeSeries_DateLocator
注
Formatter
和Locator
类用于处理 matplotlib 绘图中的刻度。 -
-
frequencies.py
:它定义用于指定时间序列对象的频率(每天,每周,每季度,每月,每年等)的代码。 -
holiday.py
:这定义了用于处理假期的函数和类-Holiday
,AbstractHolidayCalendar
和USFederalHolidayCalendar
在定义的类中。 -
index.py
:这定义了DateTimeIndex
类。 -
interval.py
:这定义了Interval
,PeriodInterval
和IntervalIndex
类。 -
offsets.py
:这定义了各种类别,包括处理与时间有关的时间段的偏移量。 这些解释如下:DateOffset
:这是提供日期时间段功能的类的接口,例如Week
,WeekOfMonth
,LastWeekOfMonth
,QuarterOffset
,YearOffset
,Easter
,FY5253
和FY5253Quarter
。BusinessMixin
:这是业务对象的mixin
类,用于提供与时间功能相关的类。 这将由BusinessDay
类继承。BusinessDay
子类派生自BusinessMixin
和SingleConstructorOffset
,并提供了营业日的补偿。MonthOffset
:这是提供月份时间段功能的类的接口,例如MonthEnd
,MonthBegin
,BusinessMonthEnd
和BusinessMonthBegin
。MonthEnd
和MonthBegin
:这是一个月底或一个月初的日期偏移量。BusinessMonthEnd
和BusinessMonthBegin
:这是一个工作日日历的结尾或开始一个月的日期偏移。YearOffset
:此偏移量是提供年份功能的类的接口,例如YearEnd
,YearBegin
,BYearEnd
,BYearBegin
YearEnd
和YearBegin
:这是一年末或年初的日期偏移量。BYearEnd
和BYearBegin
:这是工作日日历结束或开始时一年的日期偏移。Week
:提供 1 周的补偿。WeekDay
:这提供了从工作日(Tue
)到星期几(= 2
)的映射。WeekOfMonth
和LastWeekOfMonth
:描述一个月中的一周中的日期QuarterOffset
:这是提供季度功能的类的接口,例如QuarterEnd
,QuarterrBegin
,BQuarterEnd
和BQuarterBegin
。QuarterEnd
,QuarterrBegin
,BQuarterEnd
和BQuarterBegin
:与Year*
类相同,除了时间段是季度而不是年份。FY5253
和FY5253Quarter
:这些类别描述了 52-53 周的会计年度。 这也称为 4-4-5 日历。 您可以在这个链接中获得更多信息。Easter
:这是复活节假期的DateOffset
。Tick
:这是时间单位类的基本类,例如Day
,Hour
,Minute
,Second
,Milli
,Micro
和Nano
。
-
period.py
:这为 PandasTimeSeries
定义了Period
和PeriodIndex
类。 -
plotting.py
:这定义了各种绘图函数,例如tsplot(..)
,它绘制了一个序列。 -
resample.py
:这定义了TimeGrouper
,这是用于时间间隔分组的自定义groupby
类。 -
timedeltas.py
:这定义了to_timedelta(..)
方法,该方法将其参数转换为timedelta
对象。 -
tools.py
:这定义了工具函数,例如to_datetime(..)
,parse_time_string(..)
,dateutil_parse(..)
和format(..)
。 -
util.py
:这定义了更多工具函数,如下所示:isleapyear(..)
:检查年份是否为闰年pivot_annual(..)
:这按年份分组,表示闰年
Pandas /sandbox
该模块处理将 Pandas 数据帧集成到 PyQt 框架中。 有关 PyQt 的更多信息,请访问
使用 Python 扩展来提高性能
Python 和 pandas 用户的困扰之一是,语言和模块的易用性和表达性带来显着的缺点-性能-尤其是在数字计算方面。
根据编程基准站点,Python 通常比编译后的语言(例如,用于许多算法或数据结构操作的 C/C++)慢。 例如二进制树操作。 在以下参考中,Python3 的运行速度比 n 体仿真计算的最快 C++ 实现慢 104 倍。
那么,我们如何解决这个合法而令人烦恼的问题呢? 通过编写代码中对性能敏感的部分,我们可以减轻 Python 的这种速度,同时保持我们喜欢的性能(清晰度和生产率)。 例如数字处理,C/C++ 中的算法,并通过编写 Python 扩展模块由我们的 Python 代码调用它们。
Python 扩展模块使我们能够从 Python 调用用户定义的 C/C++ 代码或库函数,从而使我们能够提高代码性能,但仍然受益于 Python 的易用性。
为了帮助我们了解什么是 Python 扩展模块,请考虑导入模块时 Python 中发生的情况。 导入语句导入一个模块,但这实际上是什么意思? 共有三种可能性,如下所示:
-
某些 Python 扩展模块在构建时会链接到解释器。
-
导入会导致 Python 将
.pyc
文件加载到内存中。.pyc
文件包含 Python 字节码。例如以下命令:In [3]: import pandaspandas.__file__ Out[3]: '/usr/lib/Python 2.7/site-packages/pandas/__init__.pyc'
-
import
语句导致将 Python 扩展模块加载到内存中。.so
(共享库)文件由机器代码组成。 例如,请参考以下命令:In [4]: import mathmath.__file__ Out[4]: '/usr/lib/Python 2.7/lib-dynload/math.so'
我们将集中讨论第三种可能性。 即使我们正在处理从 C 编译的二进制共享对象,我们也可以将其作为 Python 模块导入,这显示了 Python 扩展的强大功能-应用可以从 Python 机器代码或机器代码导入模块,并且接口相同 。 Cython 和 SWIG 是使用 C 和 C++ 编写扩展的两种最受欢迎的方法。 在编写扩展时,我们包装了 C/C++ 机器代码,并将其转换为行为像纯 Python 代码的 Python 扩展模块。 在这个简短的讨论中,我们将仅专注于 Cython,因为它是专门为 Python 设计的。
Cython 是 Python 的超集,旨在允许我们在 C/C++ 中调用外部编译的代码以及在变量上声明类型,从而显着提高 Python 的性能。
Cython 命令从 Cython 源文件生成优化的 C/C++ 源文件,并将此优化的 C/C++ 源编译为 Python 扩展模块。 它提供了对 NumPy 的内置支持,并将 C 的性能与 Python 的可用性结合在一起。
我们将快速演示如何使用 Cython 大大加快代码速度。 让我们定义一个简单的斐波那契函数:
In [17]: def fibonacci(n):a,b=1,1for i in range(n):a,b=a+b,areturn a
In [18]: fibonacci(100)
Out[18]: 927372692193078999176L
In [19]: %timeit fibonacci(100)100000 loops, best of 3: 18.2 µs per loop
使用timeit
模块,我们看到每个循环花费 18.2 µs。
现在让我们在 Cython 中重写函数,通过使用以下步骤为变量指定类型:
-
首先,我们将 Cython 魔术函数导入 IPython,如下所示:
In [22]: %load_ext cythonmagic
-
接下来,我们在 Cython 中重写函数,为变量指定类型:
In [24]: %%cythondef cfibonacci(int n):cdef int i, a,bfor i in range(n):a,b=a+b,areturn a
-
让我们来计时一下新的 Cython 函数:
In [25]: %timeit cfibonacci(100)1000000 loops, best of 3: 321 ns per loopIn [26]: 18.2/0.321 Out[26]: 56.69781931464174
-
因此,我们可以看到 Cython 版本比纯 Python 版本快 57 倍!
有关使用 Cython/SWIG 或其他选项编写 Python 扩展的更多参考,请参考以下参考:
- 标题为提高性能的 Pandas 文档
- Scipy 讲义,标题为与 C 的接口
- Cython 文档
- SWIG 文档
总结
总结本章,我们浏览了 Pandas 的库层次结构,试图说明库的内部运作方式。 我们还谈到了使用 Python 扩展模块来提高代码性能的好处。
十、R 与 Pandas 的比较
本章着重于将 Pandas 与 R 进行比较,R 是对许多 Pandas 功能进行建模的统计包。 本手册旨在为希望使用 Pandas 的 R 用户以及希望复制他们在 Pandas R 代码中看到的功能的用户提供指南。 它着重介绍了 R 用户可用的一些关键功能,并通过使用一些说明性示例演示了如何在 Pandas 中实现类似功能。 本章假定您已安装 R 统计包。 如果没有,则可以从此处下载并安装。
在本章的最后,数据分析用户应该比 Pandas 更好地掌握 R 的数据分析功能,从而使他们可以过渡到 Pandas 或使用 Pandas。 本章讨论的各个主题包括:
- R 数据类型及其等效的 Pandas
- 切片和选择
- 数据类型列的算术运算
- 聚合和分组
- 匹配
- 拆分合并
- 熔化和重塑
- 因子和类别数据
R 数据类型
R 具有五种原始或原子类型:
- 字符
- 数字
- 整数
- 复数
- 逻辑/布尔
它还具有以下更复杂的容器类型:
- 向量:类似于
numpy.array
。 它只能包含相同类型的对象。 - 列表:这是一个异构容器。 相当于 Pandas 的序列。
- 数据帧 :这是一个异构 2D 容器,等效于 Pandas 数据帧
- 矩阵:它是向量的均质 2D 版本。 它类似于
numpy.matrix
。
在本章中,我们将重点关注列表和数据帧,它们具有与序列和数据帧等效的 Pandas。
注意
有关 R 数据类型的更多信息,请参考这个文档。
对于 NumPy 数据类型,请参考这个文档和这个文档。
R 列表
R 列表可以显式创建为列表声明,如下所示:
>h_lst<- list(23,'donkey',5.6,1+4i,TRUE)
>h_lst
[[1]]
[1] 23[[2]]
[1] "donkey"[[3]]
[1] 5.6[[4]]
[1] 1+4i[[5]]
[1] TRUE>typeof(h_lst)
[1] "list"
以下是其在 Pandas 中的序列,包括创建列表和从中创建序列:
In [8]: h_list=[23, 'donkey', 5.6,1+4j, True]
In [9]: import pandas as pdh_ser=pd.Series(h_list)
In [10]: h_ser
Out[10]: 0 231 donkey2 5.63 (1+4j)4 True
dtype: object
Pandas 数组索引从 0 开始,而 R 则以 1 开始。下面是一个示例:
In [11]: type(h_ser)
Out[11]: pandas.core.series.Series
R 数据帧
我们可以通过调用data.frame()
构造器来构造 R 数据帧,如下所示:
>stocks_table<- data.frame(Symbol=c('GOOG','AMZN','FB','AAPL','TWTR','NFLX','LINKD'), Price=c(518.7,307.82,74.9,109.7,37.1,334.48,219.9),
MarketCap=c(352.8,142.29,216.98,643.55,23.54,20.15,27.31))>stocks_table
Symbol PriceMarketCap
1 GOOG 518.70 352.80
2 AMZN 307.82 142.29
3 FB 74.90 216.98
4 AAPL 109.70 643.55
5 TWTR 37.10 23.54
6 NFLX 334.48 20.15
7 LINKD 219.90 27.31
在这里,我们构造一个 Pandas 数据帧并显示它:
In [29]: stocks_df=pd.DataFrame({'Symbol':['GOOG','AMZN','FB','AAPL', 'TWTR','NFLX','LNKD'],'Price':[518.7,307.82,74.9,109.7,37.1,334.48,219.9],
'MarketCap($B)' : [352.8,142.29,216.98,643.55,23.54,20.15,27.31]})
stocks_df=stocks_df.reindex_axis(sorted(stocks_df.columns,reverse=True),axis=1)
stocks_df
Out[29]:
Symbol PriceMarketCap($B)
0 GOOG 518.70 352.80
1 AMZN 307.82 142.29
2 FB 74.90 216.98
3 AAPL 109.70 643.55
4 TWTR 37.10 23.54
5 NFLX 334.48 20.15
6 LNKD219.90 27.31
切片和选择
在 R 中,我们通过以下三种方式对对象进行切片:
[
:这始终返回与原始对象相同类型的对象,并且可以用于选择多个元素。[[
:用于提取列表或数据帧的元素; 并且只能用于提取单个元素:返回的元素的类型不一定是列表或数据帧。$
:用于按名称提取列表或数据帧的元素,类似于[[
。
这是 R 中的一些切片示例及其在 Pandas 中的等效示例:
R 矩阵和 NumPy 数组比较
让我们看看 R 中的矩阵创建和选择:
>r_mat<- matrix(2:13,4,3)
>r_mat[,1] [,2] [,3]
[1,] 2 6 10
[2,] 3 7 11
[3,] 4 8 12
[4,] 5 9 13
要选择第一行,我们编写:
>r_mat[1,]
[1] 2 6 10
要选择第二列,我们使用以下命令:
>r_mat[,2]
[1] 6 7 8 9
现在让我们来看一下 NumPy 数组的创建和选择:
In [60]: a=np.array(range(2,6))b=np.array(range(6,10))c=np.array(range(10,14))
In [66]: np_ar=np.column_stack([a,b,c])
np_ar
Out[66]: array([[ 2, 6, 10],
[ 3, 7, 11],
[ 4, 8, 12],
[ 5, 9, 13]])
要选择第一行,请编写以下命令:
In [79]: np_ar[0,]
Out[79]: array([ 2, 6, 10])
注
索引在 R 和 pandas/NumPy 中是不同的。
在 R 中,索引从 1 开始,而在 pandas/NumPy 中,索引从 0 开始。因此,当从 R 转换到 pandas/NumPy 时,我们必须从所有索引中减去 1。
要选择第二列,请编写以下命令:
In [81]: np_ar[:,1]
Out[81]: array([6, 7, 8, 9])
另一个选择是先转置数组,然后选择列,如下所示:
In [80]: np_ar.T[1,]
Out[80]: array([6, 7, 8, 9])
R 列表和 Pandas 序列比较
这是在 R 中创建和选择列表的示例:
>cal_lst<- list(weekdays=1:8, mth='jan')
>cal_lst
$weekdays
[1] 1 2 3 4 5 6 7 8$mth
[1] "jan">cal_lst[1]
$weekdays
[1] 1 2 3 4 5 6 7 8>cal_lst[[1]]
[1] 1 2 3 4 5 6 7 8>cal_lst[2]
$mth
[1] "jan"
Pandas 的序列创建和选择如下:
In [92]: cal_df= pd.Series({'weekdays':range(1,8), 'mth':'jan'})
In [93]: cal_df
Out[93]: mthjan
weekdays [1, 2, 3, 4, 5, 6, 7]
dtype: objectIn [97]: cal_df[0]
Out[97]: 'jan'In [95]: cal_df[1]
Out[95]: [1, 2, 3, 4, 5, 6, 7]In [96]: cal_df[[1]]
Out[96]: weekdays [1, 2, 3, 4, 5, 6, 7]
dtype: object
在这里,从[]
和[[]]
运算符的角度来看,我们看到 R 列表和 Pandas 序列之间的差异。 我们可以通过考虑第二项(字符串)来看到差异。
对于 R,[]
运算符产生一个容器类型,即包含字符串的列表,而[[]]
产生一个原子类型:在这种情况下,一个字符如下:
>typeof(cal_lst[2])
[1] "list"
>typeof(cal_lst[[2]])
[1] "character"
对于 Pandas,情况恰恰相反:[]
产生原子类型,而[[]]
产生复杂类型,即序列的序列,如下所示:
In [99]: type(cal_df[0])
Out[99]: strIn [101]: type(cal_df[[0]])
Out[101]: pandas.core.series.Series
在 R 和 pandas 中,都可以指定列名称以获得元素。
在 R 中指定列名
在 R 中,可以使用$
运算符前面的列名来完成此操作,如下所示:
>cal_lst$mth
[1] "jan"
> cal_lst$'mth'
[1] "jan"
在 Pandas 中指定列名
在 Pandas 中,我们以通常的方式将元素子集括在方括号中:
In [111]: cal_df['mth']
Out[111]: 'jan'
R 和 Pandas 不同的一个区域是嵌套元素的子集。 例如,要从工作日获得第 4 天,我们必须在 R 中使用[[]]
运算符:
>cal_lst[[1]][[4]]
[1] 4>cal_lst[[c(1,4)]]
[1] 4
但是,对于 Pandas,我们可以使用双倍[]
:
In [132]: cal_df[1][3]
Out[132]: 4
R 的数据帧与 Pandas 的数据帧
在 R 数据帧和 Pandas 数据帧中选择数据遵循类似的脚本。 以下部分说明了我们如何从这两种方法中执行多列选择。
R 中的多列选择
在 R 中,我们通过在方括号内的向量中声明多列来指定要选择的列:
>stocks_table[c('Symbol','Price')]
Symbol Price
1 GOOG 518.70
2 AMZN 307.82
3 FB 74.90
4 AAPL 109.70
5 TWTR 37.10
6 NFLX 334.48
7 LINKD 219.90>stocks_table[,c('Symbol','Price')]
Symbol Price
1 GOOG 518.70
2 AMZN 307.82
3 FB 74.90
4 AAPL 109.70
5 TWTR 37.10
6 NFLX 334.48
7 LINKD 219.90
Pandas 中的多列选择
在 Pandas 中,我们以通常的方式将元素子集括在方括号中:
In [140]: stocks_df[['Symbol','Price']]
Out[140]:Symbol Price
0 GOOG 518.70
1 AMZN 307.82
2 FB 74.90
3 AAPL 109.70
4 TWTR 37.10
5 NFLX 334.48
6 LNKD 219.90In [145]: stocks_df.loc[:,['Symbol','Price']]
Out[145]: Symbol Price
0 GOOG 518.70
1 AMZN 307.82
2 FB 74.90
3 AAPL 109.70
4 TWTR 37.10
5 NFLX 334.48
6 LNKD 219.90
列上的算术运算
在 R 和 pandas 中,我们可以类似的方式在数据列中应用算术运算。 因此,我们可以对两个或多个数据帧中相应位置的元素执行算术运算,例如加法或减法。
在这里,我们在 R 中构造一个带有x
和y
列的数据帧,然后从x
列中减去y
列:
>norm_df<- data.frame(x=rnorm(7,0,1), y=rnorm(7,0,1))
>norm_df$x - norm_df$y
[1] -1.3870730 2.4681458 -4.6991395 0.2978311 -0.8492245 1.5851009 -1.4620324
R 中的with
运算符也具有与算术运算相同的效果:
>with(norm_df,x-y)
[1] -1.3870730 2.4681458 -4.6991395 0.2978311 -0.8492245 1.5851009 -1.4620324
在 Pandas 中,可以对列执行相同的算术运算,并且等效运算符为eval
:
In [10]: import pandas as pdimport numpy as np
df = pd.DataFrame({'x': np.random.normal(0,1,size=7), 'y': np.random.normal(0,1,size=7)})In [11]: df.x-df.y
Out[11]: 0 -0.1073131 0.6175132 -1.5178273 0.5658044 -1.6305345 0.1019006 0.775186
dtype: float64In [12]: df.eval('x-y')
Out[12]: 0 -0.1073131 0.6175132 -1.5178273 0.5658044 -1.6305345 0.1019006 0.775186
dtype: float64
聚合和分组
有时,我们可能希望将数据拆分为子集,并对每个子集应用平均值,最大值或最小值之类的函数。 在 R 中,我们可以通过aggregate
或tapply
函数执行此操作。
在这里,我们将使用 2014 年欧洲冠军联赛半决赛进入的四家具乐部前五名前锋的统计数据集的示例。我们将使用它来说明 R 及其得分的聚合,在 Pandas 中具有等效的分组功能。
R 中的聚合
在 R 中,聚合是使用以下命令完成的:
> goal_stats=read.csv('champ_league_stats_semifinalists.csv')
>goal_statsClub Player Goals GamesPlayed
1 Atletico Madrid Diego Costa 8 9
2 Atletico Madrid ArdaTuran 4 9
3 Atletico Madrid RaúlGarcía 4 12
4 Atletico Madrid AdriánLópez 2 9
5 Atletico Madrid Diego Godín 2 10
6 Real Madrid Cristiano Ronaldo 17 11
7 Real Madrid Gareth Bale 6 12
8 Real Madrid Karim Benzema 5 11
9 Real Madrid Isco 3 12
10 Real Madrid Ángel Di María 3 11
11 Bayern Munich Thomas Müller 5 12
12 Bayern Munich ArjenRobben 4 10
13 Bayern Munich Mario Götze 3 11
14 Bayern Munich Bastian Schweinsteiger 3 8
15 Bayern Munich Mario Mandžukić 3 10
16 Chelsea Fernando Torres 4 9
17 Chelsea Demba Ba 3 6
18 Chelsea Samuel Eto'o 3 9
19 Chelsea Eden Hazard 2 9
20 Chelsea Ramires 2 10
现在,我们可以计算每个前锋的每场比赛进球数比例,以衡量他们在进球前的最后期限:
>goal_stats$GoalsPerGame<- goal_stats$Goals/goal_stats$GamesPlayed
>goal_statsClub Player Goals GamesPlayedGoalsPerGame
1 Atletico Madrid Diego Costa 8 9 0.8888889
2 Atletico Madrid ArdaTuran 4 9 0.4444444
3 Atletico Madrid RaúlGarcía 4 12 0.3333333
4 Atletico Madrid AdriánLópez 2 9 0.2222222
5 Atletico Madrid Diego Godín 2 10 0.2000000
6 Real Madrid Cristiano Ronaldo 17 11 1.5454545
7 Real Madrid Gareth Bale 6 12 0.5000000
8 Real Madrid Karim Benzema 5 11 0.4545455
9 Real Madrid Isco 3 12 0.2500000
10 Real Madrid Ángel Di María 3 11 0.2727273
11 Bayern Munich Thomas Müller 5 12 0.4166667
12 Bayern Munich ArjenRobben 4 10 0.4000000
13 Bayern Munich MarioGötze 3 11 0.2727273
14 Bayern Munich Bastian Schweinsteiger 3 8 0.3750000
15 Bayern Munich MarioMandžukić 3 10 0.3000000
16 Chelsea Fernando Torres 4 9 0.4444444
17 Chelsea Demba Ba 3 6 0.5000000
18 Chelsea Samuel Eto'o 3 9 0.3333333
19 Chelsea Eden Hazard 2 9 0.2222222
20 Chelsea Ramires 2 10 0.2000000
假设我们想知道每个团队的最高进球率。 我们将计算如下:
>aggregate(x=goal_stats[,c('GoalsPerGame')], by=list(goal_stats$Club),FUN=max)Group.1 x
1 Atletico Madrid 0.8888889
2 Bayern Munich 0.4166667
3 Chelsea 0.5000000
4 Real Madrid 1.5454545
tapply
函数用于将函数应用于由一个或多个列定义的数组或向量的子集。 tapply
函数还可以如下使用:
>tapply(goal_stats$GoalsPerGame,goal_stats$Club,max)
Atletico Madrid Bayern Munich Chelsea Real Madrid 0.8888889 0.4166667 0.5000000 1.5454545
Pandas 的分组操作
在 Pandas 中,我们可以使用GroupBy
函数获得相同的结果:
In [6]: import pandas as pd
importnumpy as np
In [7]: goal_stats_df=pd.read_csv('champ_league_stats_semifinalists.csv')In [27]: goal_stats_df['GoalsPerGame']= goal_stats_df['Goals']/goal_stats_df['GamesPlayed']In [27]: goal_stats_df['GoalsPerGame']= goal_stats_df['Goals']/goal_stats_df['GamesPlayed']In [28]: goal_stats_df
Out[28]: Club Player Goals GamesPlayedGoalsPerGame
0 Atletico Madrid Diego Costa 8 9 0.888889
1 Atletico Madrid ArdaTuran 4 9 0.444444
2 Atletico Madrid RaúlGarcía 4 12 0.333333
3 Atletico Madrid AdriánLópez 2 9 0.222222
4 Atletico Madrid Diego Godín 2 10 0.200000
5 Real Madrid Cristiano Ronaldo 17 11 1.545455
6 Real Madrid Gareth Bale 6 12 0.500000
7 Real Madrid Karim Benzema 5 11 0.454545
8 Real Madrid Isco 3 12 0.250000
9 Real Madrid Ángel Di María 3 11 0.272727
10 Bayern Munich Thomas Müller 5 12 0.416667
11 Bayern Munich ArjenRobben 4 10 0.400000
12 Bayern Munich Mario Götze 3 11 0.272727
13 Bayern Munich BastianSchweinsteiger 3 8 0.375000
14 Bayern Munich MarioMandžukić 3 10 0.300000
15 Chelsea Fernando Torres 4 9 0.444444
16 Chelsea Demba Ba 3 6 0.500000
17 Chelsea Samuel Eto'o 3 9 0.333333
18 Chelsea Eden Hazard 2 9 0.222222
19 Chelsea Ramires 2 10 0.200000In [30]: grouped = goal_stats_df.groupby('Club')In [17]: grouped['GoalsPerGame'].aggregate(np.max)
Out[17]: ClubAtletico Madrid 0.888889Bayern Munich 0.416667Chelsea 0.500000Real Madrid 1.545455Name: GoalsPerGame, dtype: float64In [22]: grouped['GoalsPerGame'].apply(np.max)Out[22]: ClubAtletico Madrid 0.888889Bayern Munich 0.416667Chelsea 0.500000Real Madrid 1.545455Name: GoalsPerGame, dtype: float64
比较 R 和 Panda 中的匹配运算符
在这里,我们将证明 R(%in%
)和 Pandas(isin()
)之间匹配运算符的等价性。 在这两种情况下,都会生成逻辑向量或序列(Pandas),该逻辑向量或序列指示找到匹配项的位置。
R %in%
运算符
在这里,我们将演示 R 中%in%
运算符的用法:
>stock_symbols=stocks_table$Symbol
>stock_symbols
[1] GOOG AMZN FB AAPL TWTR NFLX LINKD
Levels: AAPL AMZN FB GOOG LINKD NFLX TWTR>stock_symbols %in% c('GOOG','NFLX')
[1] TRUE FALSE FALSE FALSE FALSE TRUE FALSE
Pandas isin()
函数
这是使用 Pandasisin()
函数的示例:
In [11]: stock_symbols=stocks_df.Symbol
stock_symbols
Out[11]: 0 GOOG1 AMZN2 FB3 AAPL4 TWTR5 NFLX6 LNKDName: Symbol, dtype: object
In [10]: stock_symbols.isin(['GOOG','NFLX'])
Out[10]: 0 True1 False2 False3 False4 False5 True6 FalseName: Symbol, dtype: bool
逻辑子集
在 R 和 Pandas 中,执行逻辑子设置的方法不止一种。 假设我们希望显示所有平均进球数大于或等于 0.5 的球员; 也就是说,他们每两场比赛平均至少有一个进球。
R 中的逻辑子集
我们可以在 R 中执行以下操作:
-
通过逻辑切片:
>goal_stats[goal_stats$GoalsPerGame>=0.5,]Club Player Goals GamesPlayedGoalsPerGame 1 Atletico Madrid Diego Costa 8 9 0.8888889 6 Real Madrid Cristiano Ronaldo 17 11 1.5454545 7 Real Madrid Gareth Bale 6 12 0.5000000 17 Chelsea Demba Ba 3 6 0.5000000
-
通过
subset()
函数:>subset(goal_stats,GoalsPerGame>=0.5)Club Player Goals GamesPlayedGoalsPerGame 1 Atletico Madrid Diego Costa 8 9 0.8888889 6 Real Madrid Cristiano Ronaldo 17 11 1.5454545 7 Real Madrid Gareth Bale 6 12 0.5000000 17 Chelsea Demba Ba 3 6 0.5000000
Pandas 中的逻辑子集
在 Pandas 中,我们也做类似的事情:
-
逻辑切片:
In [33]: goal_stats_df[goal_stats_df['GoalsPerGame']>=0.5] Out[33]: Club Player Goals GamesPlayedGoalsPerGame 0 Atletico Madrid Diego Costa 8 9 0.888889 5 Real Madrid Cristiano Ronaldo 17 11 1.545455 6 Real Madrid Gareth Bale 6 12 0.500000 16 Chelsea Demba Ba 3 6 0.500000
-
DataFrame.query()
运算符:In [36]: goal_stats_df.query('GoalsPerGame>= 0.5') Out[36]: Club Player Goals GamesPlayedGoalsPerGame 0 Atletico Madrid Diego Costa 8 9 0.888889 5 Real Madrid Cristiano Ronaldo 17 11 1.545455 6 Real Madrid Gareth Bale 6 12 0.500000 16 Chelsea Demba Ba 3 6 0.500000
拆分合并
R 有一个名为plyr
的库,用于拆分应用合并数据分析。 plyr
库具有一个称为ddply
的函数,该函数可用于将函数应用于数据帧的子集,然后将结果组合到另一个数据帧中。
注
有关ddply
的更多信息,您可以参考这个内容。
为说明起见,让我们考虑一个最近在 R 中创建的数据集的子集,其中包含 2013 年离开纽约市的航班数据。
在 R 中执行
在这里,我们将在 R 中安装该包并实例化该库:
>install.packages('nycflights13')
...>library('nycflights13')
>dim(flights)
[1] 336776 16>head(flights,3)
year month day dep_timedep_delayarr_timearr_delay carrier tailnum flight
1 2013 1 1 517 2 830 11 UA N14228 1545
2 2013 1 1 533 4 850 20 UA N24211 1714
3 2013 1 1 542 2 923 33 AA N619AA 1141
origindestair_time distance hour minute
1 EWR IAH 227 1400 5 17
2 LGA IAH 227 1416 5 33
3 JFK MIA 160 1089 5 42> flights.data=na.omit(flights[,c('year','month','dep_delay','arr_delay','distance')])
>flights.sample<- flights.data[sample(1:nrow(flights.data),100,replace=FALSE),]>head(flights.sample,5)
year month dep_delayarr_delay distance
155501 2013 3 2 5 184
2410 2013 1 0 4 762
64158 2013 11 -7 -27 509
221447 2013 5 -5 -12 184
281887 2013 8 -1 -10 937
ddply
函数使我们能够按年和月聚合离港延误(均值,标准差):
>ddply(flights.sample,.(year,month),summarize, mean_dep_delay=round(mean(dep_delay),2), s_dep_delay=round(sd(dep_delay),2))
year month mean_dep_delaysd_dep_delay
1 2013 1 -0.20 2.28
2 2013 2 23.85 61.63
3 2013 3 10.00 34.72
4 2013 4 0.88 12.56
5 2013 5 8.56 32.42
6 2013 6 58.14 145.78
7 2013 7 25.29 58.88
8 2013 8 25.86 59.38
9 2013 9 -0.38 10.25
10 2013 10 9.31 15.27
11 2013 11 -1.09 7.73
12 2013 12 0.00 8.58
让我们将flights.sample
数据集保存到 CSV 文件中,以便我们可以使用该数据向我们展示如何在 Pandas 中执行相同的操作:
>write.csv(flights.sample,file='nycflights13_sample.csv', quote=FALSE,row.names=FALSE)
在 Pandas 中实现
为了在 Pandas 中做同样的事情,我们阅读了上一节中保存的 CSV 文件:
In [40]: flights_sample=pd.read_csv('nycflights13_sample.csv')In [41]: flights_sample.head()
Out[41]: year month dep_delayarr_delay distance
0 2013 3 2 5 184
1 2013 1 0 4 762
2 2013 11 -7 -27 509
3 2013 5 -5 -12 184
4 2013 8 -1 -10 937
通过使用GroupBy()
运算符,可以达到与ddply
相同的效果:
In [44]: pd.set_option('precision',3)
In [45]: grouped = flights_sample_df.groupby(['year','month'])In [48]: grouped['dep_delay'].agg([np.mean, np.std])Out[48]: mean std
year month
2013 1 -0.20 2.282 23.85 61.633 10.00 34.724 0.88 12.565 8.56 32.426 58.14 145.787 25.29 58.888 25.86 59.389 -0.38 10.2510 9.31 15.2711 -1.09 7.7312 0.00 8.58
熔化和重塑
melt
函数将数据转换为宽格式,以单列组成唯一的 ID 变量组合。
R 的melt()
函数
在这里,我们演示了 R 中melt()
函数的用法。它生成长格式数据,其中的行是唯一的变量值组合:
>sample4=head(flights.sample,4)[c('year','month','dep_delay','arr_delay')]
> sample4year month dep_delay arr_delay155501 2013 3 2 52410 2013 1 0 464158 2013 11 -7 -27221447 2013 5 -5 -12>melt(sample4,id=c('year','month'))year month variable value1 2013 3 dep_delay 22 2013 1 dep_delay 03 2013 11 dep_delay -74 2013 5 dep_delay -55 2013 3 arr_delay 56 2013 1 arr_delay 47 2013 11 arr_delay -278 2013 5 arr_delay -12
>
有关更多信息,您可以参考这个内容。
Pandas 的melt()
函数
在 Pandas 中,melt
函数类似:
In [55]: sample_4_df=flights_sample_df[['year','month','dep_delay', \\
'arr_delay']].head(4)
In [56]: sample_4_df
Out[56]: year month dep_delay arr_delay0 2013 3 2 51 2013 1 0 42 2013 11 -7 -273 2013 5 -5 -12In [59]: pd.melt(sample_4_df,id_vars=['year','month'])
Out[59]: year month variable value0 2013 3 dep_delay 21 2013 1 dep_delay 02 2013 11 dep_delay -73 2013 5 dep_delay -54 2013 3 arr_delay 55 2013 1 arr_delay 46 2013 11 arr_delay -277 2013 5 arr_delay -12
此信息的参考来自这里。
因子/类别数据
R 将类别变量称为因子,并且cut()
函数使我们能够将连续的数值变量划分为范围,并将范围视为因子或类别变量,或者将类别变量分类为更大的桶。
使用cut()
的 R 示例
这是 R 中的示例:
clinical.trial<- data.frame(patient = 1:1000,
age = rnorm(1000, mean = 50, sd = 5),
year.enroll = sample(paste("19", 80:99, sep = ""),1000, replace = TRUE))>clinical.trial<- data.frame(patient = 1:1000,
+ age = rnorm(1000, mean = 50, sd = 5),
+ year.enroll = sample(paste("19", 80:99, sep = ""),
+ 1000, replace = TRUE))
>summary(clinical.trial)
patient age year.enrollMin. : 1.0 Min. :31.14 1995 : 61 1st Qu.: 250.8 1st Qu.:46.77 1989 : 60
Median : 500.5 Median :50.14 1985 : 57 Mean : 500.5 Mean :50.14 1988 : 57 3rd Qu.: 750.2 3rd Qu.:53.50 1990 : 56 Max. :1000.0 Max. :70.15 1991 : 55 (Other):654
>ctcut<- cut(clinical.trial$age, breaks = 5)> table(ctcut)
ctcut
(31.1,38.9] (38.9,46.7] (46.7,54.6] (54.6,62.4] (62.4,70.2]15 232 558 186 9
可以在这个位置找到上述数据的参考。
Pandas 解决方案
这等效于先前在 Pandas 中介绍的cut()
函数(仅适用于 0.15+版本):
In [79]: pd.set_option('precision',4)
clinical_trial=pd.DataFrame({'patient':range(1,1001), 'age' : np.random.normal(50,5,size=1000),'year_enroll': [str(x) for x in np.random.choice(range(1980,2000),size=1000,replace=True)]})In [80]: clinical_trial.describe()
Out[80]: age patient
count 1000.000 1000.000
mean 50.089 500.500
std 4.909 288.819
min 29.944 1.00025% 46.572 250.75050% 50.314 500.50075% 53.320 750.250
max 63.458 1000.000In [81]: clinical_trial.describe(include=['O'])
Out[81]: year_enroll
count 1000
unique 20
top 1992
freq 62In [82]: clinical_trial.year_enroll.value_counts()[:6]
Out[82]: 1992 621985 611986 591994 591983 581991 58
dtype: int64
In [83]: ctcut=pd.cut(clinical_trial['age'], 5)
In [84]: ctcut.head()
Out[84]: 0 (43.349, 50.052]1 (50.052, 56.755]2 (50.052, 56.755]3 (43.349, 50.052]4 (50.052, 56.755]Name: age, dtype: categoryCategories (5, object): [(29.91, 36.646] < (36.646, 43.349] < (43.349, 50.052] < (50.052, 56.755] < (56.755, 63.458]]In [85]: ctcut.value_counts().sort_index()
Out[85]: (29.91, 36.646] 3(36.646, 43.349] 82(43.349, 50.052] 396(50.052, 56.755] 434(56.755, 63.458] 85
dtype: int64
总结
在本章中,我们试图将 R 中的关键特征与其等效的 Pandas 进行比较,以实现以下目标:
- 帮助可能希望在 Pandas 中复制相同功能的 R 用户
- 帮助阅读了一些 R 代码,并可能希望将其重写为 Pandas 代码的任何用户
在下一章中,我们将通过对scikit-learn
库进行机器学习的简要介绍来结束本书,并展示 Pandas 如何适合该框架。 可在以下位置找到本章的参考文档。
十一、机器学习简介
本章将引导用户进行机器学习的旋风之旅,重点介绍如何使用 pandas 库作为可用于预处理机器学习程序使用的数据的工具。 它还向用户介绍了scikit-learn
库,它是 Python 中最流行的机器学习工具包。
在本章中,我们通过将机器学习技术应用于一个众所周知的问题来进行说明,该问题涉及对上个世纪之初在泰坦尼克号灾难中幸存的乘客进行分类。 本章讨论的各个主题包括:
- Pandas 在机器学习中的作用
- 安装
scikit-learn
- 机器学习概念简介
- 机器学习的应用 – Kaggle 泰坦尼克竞赛
- Pandas 数据分析和预处理
- 解决泰坦尼克号问题的朴素方法
scikit-learn
ML 分类器接口- 监督学习算法
- 无监督学习算法
Pandas 在机器学习中的作用
我们将考虑用于机器学习的库称为scikit-learn
。 scikit-learn
Python 库提供了广泛的机器学习算法库,可用于创建从数据输入中学习的自适应程序。
但是,在scikit-learn
可以使用此数据之前,必须对其进行一些预处理。 这就是 Pandas 进来的地方。在将数据传递给scikit-learn
中实现的算法之前,Pandas 可用于预处理和过滤数据。
scikit-learn 的安装
如第 2 章,“安装 Pandas 和支持软件”所述,安装 Pandas 及其随附库的最简单方法是使用第三方发行版,例如 Anaconda 来完成它。 安装scikit-learn
应该没有什么不同。 我将简要介绍从 Anaconda 开始在各种平台和第三方发行版上进行安装的步骤。 scikit-learn
库需要以下库:
- Python 2.6.x 或更高版本
- NumPy 1.6.1 或更高版本
- SciPy 0.9 或更高
假设您已经按照第 2 章,“安装 Pandas 和支持软件”中的说明安装了 Pandas,则这些依赖项应该已经存在。
通过 Anaconda 安装
您可以通过运行 conda Python 包管理器在 Anaconda 上安装scikit-learn
:
conda install scikit-learn
在 Unix 上安装(Linux/MacOSX)
对于 Unix,最好从源代码安装(需要 C 编译器)。 假设已经安装了 pandas 和 NumPy,并且所需的依赖库已经到位,则可以通过运行以下命令通过 Git 安装scikit-learn
:
git clone https://github.com/scikit-learn/scikit-learn.git
cd scikitlearn
python setup.py install
还可以使用PyPi
中的pip
将 Pandas 安装在 Unix 上:
pip install pandas
在 Windows 上安装
要在 Windows 上安装,您可以打开控制台并运行以下命令:
pip install –U scikit-learn
注意
有关安装的更多详细信息,您可以在以下位置查看官方scikit-learn
文档。
您还可以在这里查看scikit-learn
Git 存储库的README
文件。
机器学习导论
机器学习是创建从数据中学习的软件程序的艺术。 更正式地讲,它可以定义为构建使用可调参数来提高预测性能的自适应程序的实践。 它是人工智能的一个子领域。
我们可以根据他们要解决的问题的类型来分离机器学习程序。 这些问题被适当地称为学习问题。
从广义上讲,这些问题的两类分别称为有监督的学习问题和无监督的学习问题。 此外,有些混合问题的方面涉及这两个类别。
学习问题的输入由 n 行的数据集组成。 每一行代表一个样本,并且可能涉及一个或多个称为属性或特征的字段。
数据集可以规范地描述为由 n 个样本组成,每个样本由 m 个特征组成。 以下论文对机器学习进行了更详细的介绍:
关于机器学习的一些有用的知识
有监督与无监督学习
对于有监督的学习问题,学习问题的输入是由带标签数据组成的数据集。 这意味着我们具有已知其值的输出。 向学习程序提供输入样本及其对应的输出,其目的是破译它们之间的关系。 这样的输入称为标记数据。 监督学习问题包括以下内容:
- 分类:学习的属性是类别(标称)或离散的
- 回归:学习的属性为数字/连续
在无监督学习或数据挖掘中,向学习程序提供输入,但没有相应的输出。 该输入数据称为未标记数据。 学习程序的目标是学习或解密隐藏的标签。 这些问题包括:
- 聚类
- 降维
使用文档分类的概述
机器学习技术的常见用法是在文档分类领域。 机器学习的两个主要类别可以应用于此问题-有监督和无监督学习。
监督学习
输入集合中的每个文档都分配有一个类别,即标签。 学习程序/算法使用输入的文档集合来学习如何对没有标签的另一组文档进行预测。 该方法称为分类。
无监督学习
输入集合中的文档未分配到类别。 因此,它们没有标签。 学习程序将其作为输入,并尝试对其聚类来发现相关或相似文档的组。 该方法称为聚类。
机器学习系统如何学习
机器学习系统利用所谓的分类器从数据中学习。 分类器是采用称为特征值的矩阵并生成输出向量(也称为类)的接口。 这些特征值可以是离散的或连续的。 分类器包含三个核心组件:
- 表示形式:这是什么类型的分类器?
- 评估:分类器的质量如何?
- 优化:如何在替代方案中进行搜索?
机器学习的应用 – Kaggle 泰坦尼克竞赛
为了说明我们如何在机器学习旅程的开始时使用 Pandas 来帮助我们,我们将其应用于一个经典问题,该问题位于 Kaggle 网站。 Kaggle 是针对机器学习问题的竞争平台。 Kaggle 背后的想法是使对使用数据解决预测分析问题感兴趣的公司将其数据发布到 Kaggle 上,并邀请数据科学家提出针对其问题的建议解决方案。 竞争可能会持续一段时间,并且竞争对手的排名会发布在排行榜上。 在比赛结束时,排名第一的选手将获得现金奖励。
为了说明如何使用scikit-learn
将 Pandas 用于机器学习,我们将研究的经典问题是在 Kaggle 托管的泰坦尼克号的灾难问题,作为经典的入门机器学习问题。 问题涉及的数据集是原始数据集。 因此,在将数据提交为scikit-learn
中实现的机器学习算法的输入之前,Pandas 在数据的预处理和清理中非常有用。
泰坦尼克号:灾难问题中的机器学习
泰坦尼克号的数据集包括厄运之旅的乘客清单,各种特征以及指示变量,该变量指示乘客是否幸免于沉船。 问题的本质是,在给定乘客及其相关特征的情况下,能够预测该乘客是否在泰坦尼克号沉没中幸免于难。 请删除这句话。
数据包含两个数据集:一个训练数据集和另一个测试数据集。 训练数据集包含 891 个乘客案例,测试数据集包含 491 个乘客案例。
训练数据集还包含 11 个变量,其中 10 个是特征,还有 1 个依存/指示变量Survived
,它指示乘客是否幸免于难。
特征变量如下:
PassengerID
Name
Sex
Pclass
(乘客舱位)Carin
Parch
(父母和子女的数量)Age
Sibsp
(兄弟姐妹数)Embarked
我们可以使用 Pandas 通过以下方式帮助我们预处理数据:
- 数据清理和某些变量的分类
- 排除不必要的特征,这些特征显然与乘客的生存能力无关,例如,其名称
- 处理缺失的数据
我们可以使用多种算法来解决此问题。 它们如下:
- 决策树
- 神经网络
- 随机森林
- 支持向量机
过拟合的问题
过拟合是机器学习中的一个众所周知的问题,程序会记住它作为输入输入的特定数据,从而在训练数据上获得完美的结果,而在测试数据上获得令人讨厌的结果。
为了防止过拟合,可以在训练阶段使用 10 倍交叉验证技术在数据中引入可变性。
使用 Pandas 进行数据分析和预处理
在本节中,我们将使用 Pandas 对数据进行一些分析和预处理,然后再将其作为输入提交给scikit-learn
。
检查数据
为了开始对数据进行预处理,让我们读取训练数据集并检查其外观。
在这里,我们将训练数据集读入 Pandas 数据帧并显示第一行:
In [2]: import pandas as pdimport numpy as np
# For .read_csv, always use header=0 when you know row 0 is the header rowtrain_df = pd.read_csv('csv/train.csv', header=0)
In [3]: train_df.head(3)
输出如下:
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-aaYsOquj-1681366316324)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00180.jpeg)]
因此,我们可以看到各种特征:PassengerId
,PClass
,Name
,Sex
,Age
,Sibsp
,Parch
,Ticket
,Fare
,Cabin
和Embarked
。 立即想到的一个问题是:哪些特征可能会影响乘客是否幸存?
显然,PassengerID
,Ticket
和Name
不应影响生存能力,因为它们是标识符变量。 我们将在分析中跳过这些内容。
处理缺失值
我们必须在数据集中进行机器学习的一个问题是如何处理训练集中的缺失值。
让我们直观地确定特征集中缺少值的位置。
为此,我们可以使用由 Tom Augspurger 编写的 R 中的missmap
函数的等效项。 下图以直观的方式显示了各种特征缺少的数据量:
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-ASYVIuMo-1681366316325)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00181.jpeg)]
有关更多信息和用于生成此数据的代码,请参见这个内容。
我们还可以计算每个特征缺少多少数据:
In [83]: missing_perc=train_df.apply(lambda x: 100*(1-x.count().sum()/(1.0*len(x))))
In [85]: sorted_missing_perc=missing_perc.order(ascending=False)sorted_missing_perc
Out[85]: Cabin 77.104377Age 19.865320Embarked 0.224467Fare 0.000000Ticket 0.000000Parch 0.000000SibSp 0.000000Sex 0.000000Name 0.000000Pclass 0.000000Survived 0.000000PassengerId 0.000000dtype: float64
因此,我们可以看到大多数Cabin
数据丢失了(77%),而大约 20% Age
数据丢失了。 然后,我们决定从学习特征集中删除Cabin
数据,因为该数据太稀疏而无用。
让我们对我们要检查的各种特征做进一步的细分。 对于分类/离散特征,我们使用条形图。 对于连续值特征,我们使用直方图:
In [137]: import randombar_width=0.1categories_map={'Pclass':{'First':1,'Second':2, 'Third':3},'Sex':{'Female':'female','Male':'male'},'Survived':{'Perished':0,'Survived':1},'Embarked':{'Cherbourg':'C','Queenstown':'Q','Southampton':'S'},'SibSp': { str(x):x for x in [0,1,2,3,4,5,8]},'Parch': {str(x):x for x in range(7)}}colors=['red','green','blue','yellow','magenta','orange']subplots=[111,211,311,411,511,611,711,811]cIdx=0fig,ax=plt.subplots(len(subplots),figsize=(10,12))keyorder = ['Survived','Sex','Pclass','Embarked','SibSp','Parch']for category_key,category_items in sorted(categories_map.iteritems(),key=lambda i:keyorder.index(i[0])):num_bars=len(category_items)index=np.arange(num_bars)idx=0for cat_name,cat_val in sorted(category_items.iteritems()):ax[cIdx].bar(idx,len(train_df[train_df[category_key]==cat_val]), label=cat_name,color=np.random.rand(3,1))idx+=1ax[cIdx].set_title('%s Breakdown' % category_key)xlabels=sorted(category_items.keys()) ax[cIdx].set_xticks(index+bar_width)ax[cIdx].set_xticklabels(xlabels)ax[cIdx].set_ylabel('Count')cIdx +=1
fig.subplots_adjust(hspace=0.8)
for hcat in ['Age','Fare']:ax[cIdx].hist(train_df[hcat].dropna(),color=np.random.rand(3,1))ax[cIdx].set_title('%s Breakdown' % hcat)#ax[cIdx].set_xlabel(hcat)ax[cIdx].set_ylabel('Frequency')cIdx +=1fig.subplots_adjust(hspace=0.8)
plt.show()
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-DHygRrkB-1681366316325)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00182.jpeg)]
从上图中的数据和插图中,我们可以观察到以下内容:
- 死亡人数是幸存人数的两倍(62% 对 38%)。
- 男乘客的数量大约是女乘客的两倍(65% 对 35%)。
- 与头等舱和二等舱相比,三等舱的乘客增加了约 20% (55% 对 45%)。
- 大多数乘客都是单人,也就是说,船上没有孩子,父母,兄弟姐妹或配偶。
这些观察结果可能会导致我们更深入地研究并调查生存机会,性别与票价等级之间是否存在某些相关性,特别是如果我们考虑到泰坦尼克号实行的是妇女和儿童优先政策,以及泰坦尼克号所载救生艇(20 艘)比设计的(32 艘)少的事实。
有鉴于此,让我们进一步检查生存率与其中某些特征之间的关系。 我们从性别开始:
In [85]: from collections import OrderedDictnum_passengers=len(train_df)num_men=len(train_df[train_df['Sex']=='male'])men_survived=train_df[(train_df['Survived']==1 ) & (train_df['Sex']=='male')]num_men_survived=len(men_survived)num_men_perished=num_men-num_men_survivednum_women=num_passengers-num_menwomen_survived=train_df[(train_df['Survived']==1) & (train_df['Sex']=='female')]num_women_survived=len(women_survived)num_women_perished=num_women-num_women_survivedgender_survival_dict=OrderedDict()gender_survival_dict['Survived']={'Men':num_men_survived,'Women':num_women_survived}gender_survival_dict['Perished']={'Men':num_men_perished,'Women':num_women_perished}gender_survival_dict['Survival Rate']= {'Men' : round(100.0*num_men_survived/num_men,2),'Women':round(100.0*num_women_survived/num_women,2)}
pd.DataFrame(gender_survival_dict)
Out[85]:
性别 | 幸存 | 灭亡 | 存活率 |
---|---|---|---|
男 | 109 | 468 | 18.89 |
女 | 233 | 81 | 74.2 |
现在,使用以下命令以条形图说明此数据:
In [76]: #code to display survival by genderfig = plt.figure()ax = fig.add_subplot(111)perished_data=[num_men_perished, num_women_perished]survived_data=[num_men_survived, num_women_survived]N=2ind = np.arange(N) # the x locations for the groupswidth = 0.35survived_rects = ax.barh(ind, survived_data, width,color='green')perished_rects = ax.barh(ind+width, perished_data, width,color='red')ax.set_xlabel('Count')ax.set_title('Count of Survival by Gender')yTickMarks = ['Men','Women']ax.set_yticks(ind+width)ytickNames = ax.set_yticklabels(yTickMarks)plt.setp(ytickNames, rotation=45, fontsize=10)## add a legendax.legend((survived_rects[0], perished_rects[0]), ('Survived', 'Perished') )plt.show()
上面的代码生成以下条形图:
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-CkAIj0gO-1681366316325)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00183.jpeg)]
从前面的绘图中,我们可以看到,大多数女性幸存下来(74%),而大多数男性丧生(只有 19% 幸存)。
这导致我们得出以下结论:乘客的性别可能是乘客是否幸存的一个重要因素。
接下来,让我们看一下旅客舱位。 首先,我们生成三种乘客类别中每一种的幸存和死亡数据,以及生存率,并在表格中显示:
In [86]:
from collections import OrderedDict
num_passengers=len(train_df)
num_class1=len(train_df[train_df['Pclass']==1])
class1_survived=train_df[(train_df['Survived']==1 ) & (train_df['Pclass']==1)]
num_class1_survived=len(class1_survived)
num_class1_perished=num_class1-num_class1_survived
num_class2=len(train_df[train_df['Pclass']==2])
class2_survived=train_df[(train_df['Survived']==1) & (train_df['Pclass']==2)]
num_class2_survived=len(class2_survived)
num_class2_perished=num_class2-num_class2_survived
num_class3=num_passengers-num_class1-num_class2
class3_survived=train_df[(train_df['Survived']==1 ) & (train_df['Pclass']==3)]
num_class3_survived=len(class3_survived)
num_class3_perished=num_class3-num_class3_survived
pclass_survival_dict=OrderedDict()
pclass_survival_dict['Survived']={'1st Class':num_class1_survived,'2nd Class':num_class2_survived,'3rd Class':num_class3_survived}
pclass_survival_dict['Perished']={'1st Class':num_class1_perished,'2nd Class':num_class2_perished,'3rd Class':num_class3_perished}
pclass_survival_dict['Survival Rate']= {'1st Class' : round(100.0*num_class1_survived/num_class1,2),'2nd Class':round(100.0*num_class2_survived/num_class2,2),'3rd Class':round(100.0*num_class3_survived/num_class3,2),}
pd.DataFrame(pclass_survival_dict)Out[86]:
旅客舱位 | 幸存 | 灭亡 | 存活率 |
---|---|---|---|
头等 | 136 | 80 | 62.96 |
二等 | 87 | 97 | 47.28 |
三等 | 119 | 372 | 24.24 |
然后,我们可以使用matplotlib
来绘制数据,类似于前面所述的按性别划分的幸存者计数方法:
In [186]:
fig = plt.figure()
ax = fig.add_subplot(111)
perished_data=[num_class1_perished, num_class2_perished, num_class3_perished]
survived_data=[num_class1_survived, num_class2_survived, num_class3_survived]
N=3
ind = np.arange(N) # the x locations for the groups
width = 0.35
survived_rects = ax.barh(ind, survived_data, width,color='blue')
perished_rects = ax.barh(ind+width, perished_data, width,color='red')
ax.set_xlabel('Count')
ax.set_title('Survivor Count by Passenger class')
yTickMarks = ['1st Class','2nd Class', '3rd Class']
ax.set_yticks(ind+width)
ytickNames = ax.set_yticklabels(yTickMarks)
plt.setp(ytickNames, rotation=45, fontsize=10)
## add a legend
ax.legend( (survived_rects[0], perished_rects[0]), ('Survived', 'Perished'),loc=10 )
plt.show()
这将产生以下条形图:
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-A8kEl6Nn-1681366316326)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00184.jpeg)]
从前面的数据和说明中可以清楚地看出,乘客票价等级越高,生存的机会就越大。
鉴于性别和票价舱位似乎都在影响乘客的生存机会,让我们看看将这两个特征结合并绘制两者的结合时会发生什么。 为此,我们将在 Pandas 中使用crosstab
函数。
In [173]: survival_counts=pd.crosstab([train_df.Pclass,train_df.Sex],train_df.Survived.astype(bool))survival_counts
Out[173]: Survived False TruePclass Sex 1 female 3 91male 77 452 female 6 70male 91 173 female 72 72male 300 47
现在让我们使用matplotlib
显示此数据。 首先,让我们为显示目的做一些重新标记:
In [183]: survival_counts.index=survival_counts.index.set_levels([['1st', '2nd', '3rd'], ['Women', 'Men']])
In [184]: survival_counts.columns=['Perished','Survived']
现在,我们使用 Pandas 数据帧的plot
函数绘制数据:
In [185]: fig = plt.figure()ax = fig.add_subplot(111)ax.set_xlabel('Count')ax.set_title('Survivor Count by Passenger class, Gender')survival_counts.plot(kind='barh',ax=ax,width=0.75,color=['red','black'], xlim=(0,400))
Out[185]: <matplotlib.axes._subplots.AxesSubplot at 0x7f714b187e90>
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-7n6lveEW-1681366316326)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00185.jpeg)]
泰坦尼克号问题的朴素解决方法
我们对泰坦尼克号数据进行分类的首次尝试是使用朴素但非常直观的方法。 此方法涉及以下步骤:
- 选择一组特征 S ,这些特征会影响一个人是否生存。
- 对于特征的每种可能组合,请使用训练数据指示大多数病例是否存活。 这可以在所谓的生存矩阵中进行评估。
- 对于我们希望预测生存的每个测试示例,请查找对应于其特征值的特征组合,并将其预测值分配给生存表中的生存值。 这种方法是朴素的 K 最近邻方法。
根据我们之前在分析中看到的内容,似乎对存活率影响最大的三个特征是:
Pclass
Sex
PriceBucket
我们包括与乘客等级有关的乘客票价。
生存表看起来类似于以下内容:
NumberOfPeople Pclass PriceBucket Sex Survived
0 0 1 0 female 0
1 1 1 0 male 0
2 0 1 1 female 0
3 0 1 1 male 0
4 7 1 2 female 1
5 34 1 2 male 0
6 1 1 3 female 1
7 19 1 3 male 0
8 0 2 0 female 0
9 0 2 0 male 0
10 35 2 1 female 1
11 63 2 1 male 0
12 31 2 2 female 1
13 25 2 2 male 0
14 4 2 3 female 1
15 6 2 3 male 0
16 64 3 0 female 1
17 256 3 0 male 0
18 43 3 1 female 1
19 38 3 1 male 0
20 21 3 2 female 0
21 24 3 2 male 0
22 10 3 3 female 0
23 5 3 3 male 0
可以在随附的文件survival_data.py
中找到用于生成此表的代码。 要查看我们如何使用该表,让我们看一下测试数据的摘要:
In [192]: test_df.head(3)[['PassengerId','Pclass','Sex','Fare']]
Out[192]: PassengerId Pclass Sex Fare0 892 3 male 7.82921 893 3 female 7.00002 894 2 male 9.6875
对于 892 乘客,我们看到他是男性,票价为 7.8292,并且他乘坐的是三等舱。
因此,用于此乘客的生存表查找的关键字是{Sex ='male',Pclass = 3,PriceBucket = 0}
(因为 7.8292 落入存储桶 0 中)。
如果我们在生存表(第 17 行)中查找与此键对应的生存值,则会看到该值是0 = Perished
; 这就是我们将要预测的值。
类似地,对于乘客 893,我们有Key = {Sex = 'Female',Pclass = 3,PriceBucket = 0}
。
这对应于第 16 行,因此,我们将预测 1,即生存,而她的预测生存为 1,即生存。
因此,我们的结果类似于以下命令:
> head -4 csv/surv_results.csv
PassengerId,Survived
892,0
893,1
894,0
该信息的来源位于这里。
使用前面概述的生存表方法,可以在 Kaggle 上达到 0.77990 的精度。
生存表方法虽然直观,却是一种非常基本的方法,仅代表机器学习中各种可能性的冰山一角。
在以下各节中,我们将带动各种机器学习算法的旋风之旅,这些算法将帮助您(读者)对机器学习领域中可用的特征有所了解。
scikit-learn ML /分类器接口
我们将深入研究机器学习的基本原理,并通过scikit-learn
基本 API 演示这些原理的用法。
scikit-learn
库具有一个估计器接口。 我们通过使用线性回归模型来说明它。 例如,考虑以下内容:
In [3]: from sklearn.linear_model import LinearRegression
实例化估计器接口以创建一个模型,在这种情况下为线性回归模型:
In [4]: model = LinearRegression(normalize=True)
In [6]: print modelLinearRegression(copy_X=True, fit_intercept=True, normalize=True)
在这里,我们指定normalize=True
,指示x
值将在回归之前进行归一化。 超参数(估计器参数)在模型创建中作为参数传递。 这是创建具有可调参数的模型的示例。
当数据与估计器拟合时,从数据获得估计的参数。 让我们首先创建一些样本训练数据,这些样本训练数据通常关于y = x/2
分布。 我们首先生成我们的x
和y
值:
In [51]: sample_size=500x = []y = []for i in range(sample_size):newVal = random.normalvariate(100,10)x.append(newVal)y.append(newVal / 2.0 + random.normalvariate(50,5))
sklearn
将num_samples × num_features
的 2D 数组作为输入,因此我们将x
数据转换为 2D 数组:
In [67]: X = np.array(x)[:,np.newaxis]X.shape
Out[67]: (500, 1)
在这种情况下,我们有 500 个样本和 1 个特征x
。 现在,我们训练/拟合模型并显示斜率(系数)和回归线的截距,即预测:
In [71]: model.fit(X,y)print "coeff=%s, intercept=%s" % (model.coef_,model.intercept_)coeff=[ 0.47071289], intercept=52.7456611783
可以如下所示:
In [65]: plt.title("Plot of linear regression line and training data")plt.xlabel('x')plt.ylabel('y')plt.scatter(X,y,marker='o', color='green', label='training data');plt.plot(X,model.predict(X), color='red', label='regression line')plt.legend(loc=2)Out[65]: [<matplotlib.lines.Line2D at 0x7f11b0752350]
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-dMLixnx2-1681366316326)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00186.jpeg)]
要总结估计器接口的基本用法,请按照下列步骤操作:
- 定义模型 –
LinearRegression
,SupportVectorMachine
,DecisionTrees
等。 您可以在此步骤中指定所需的超参数。 例如,如先前指定的normalize=True
。 - 定义模型后,可以通过在上一步中定义的模型上调用
fit(..)
方法来对数据进行训练。 - 拟合模型后,可以对测试数据调用
predict(..)
方法以进行预测或估计。 - 在监督学习问题的情况下,
predict(X)
方法将获得未标记的观察值X
,并返回预测的标记y
。
注
有关更多参考,请参见这里和这里。
监督学习算法
我们将简要介绍一些著名的监督学习算法,并了解如何将其应用于前面所述的泰坦尼克号生存预测问题。
使用 Patsy 为 scikit-learn 构建模型
在开始学习机器学习算法之前,我们需要了解一些Patsy
库。 我们将利用Patsy
设计将与scikit-learn
结合使用的特征。 Patsy
是用于创建称为设计矩阵的包。 这些设计矩阵是输入数据中要素的变换。 转换由称为公式的表达式指定,这些表达式对应于我们希望机器学习程序在学习中使用哪些特征的规范。
一个简单的例子如下:
假设我们希望y
与其他x
,a
和b
变量以及a
和b
之间的相互作用线性回归。 然后,我们可以指定模型如下:
import patsy as pts
pts.dmatrices("y ~ x + a + b + a:b", data)
在代码的前一行中,公式由以下表达式指定:y ~ x + a + b + a:b
。
注
进一步的参考,请查看这里
一般样板代码说明
在本节中,我们将介绍使用Patsy
和scikit-learn
实现以下各种算法的样板代码。 这样做的原因是以下算法的大多数代码都是可重复的。
在以下各节中,将介绍算法的工作原理,并将提供每种算法专用的代码作为本章的附件。
-
首先,通过使用以下命令行来确保我们位于正确的文件夹中。 假设工作目录位于
~/devel/Titanic
,我们有:In [17]: %cd ~/devel/Titanic/home/youruser/devel/sandbox/Learning/Kaggle/Titanic
-
在这里,我们导入所需的包并阅读我们的培训和测试数据集:
In [18]: import matplotlib.pyplot as pltimport pandas as pdimport numpy as npimport patsy as pt In [19]: train_df = pd.read_csv('csv/train.csv', header=0)test_df = pd.read_csv('csv/test.csv', header=0)
-
接下来,我们指定要提交给
Patsy
的公式:In [21]: formula1 = 'C(Pclass) + C(Sex) + Fare'formula2 = 'C(Pclass) + C(Sex)'formula3 = 'C(Sex)'formula4 = 'C(Pclass) + C(Sex) + Age + SibSp + Parch'formula5 = 'C(Pclass) + C(Sex) + Age + SibSp + Parch + C(Embarked)' formula6 = 'C(Pclass) + C(Sex) + Age + SibSp + C(Embarked)'formula7 = 'C(Pclass) + C(Sex) + SibSp + Parch + C(Embarked)'formula8 = 'C(Pclass) + C(Sex) + SibSp + Parch + C(Embarked)'In [23]: formula_map = {'PClass_Sex_Fare' : formula1,'PClass_Sex' : formula2,'Sex' : formula3,'PClass_Sex_Age_Sibsp_Parch' : formula4,'PClass_Sex_Age_Sibsp_Parch_Embarked' : formula5,'PClass_Sex_Embarked' : formula6,'PClass_Sex_Age_Parch_Embarked' : formula7,'PClass_Sex_SibSp_Parch_Embarked' : formula8}
我们将定义一个函数来帮助我们处理缺失值。 下面的函数在数据帧中查找具有空值的单元格,获取一组相似的乘客,并将空值设置为该组相似乘客的该特征的平均值。 相似的乘客定义为与具有零特征值的乘客具有相同性别和乘客等级的乘客。
In [24]:
def fill_null_vals(df,col_name):null_passengers=df[df[col_name].isnull()]passenger_id_list = null_passengers['PassengerId'].tolist()df_filled=df.copy()for pass_id in passenger_id_list:idx=df[df['PassengerId']==pass_id].index[0]similar_passengers = df[(df['Sex']== null_passengers['Sex'][idx]) & (df['Pclass']==null_passengers['Pclass'][idx])]mean_val = np.mean(similar_passengers[col_name].dropna())df_filled.loc[idx,col_name]=mean_valreturn df_filled
在这里,我们创建培训和测试数据帧的填充版本。
我们的测试数据帧是拟合的scikit-learn
模型将生成的预测,以产生输出,并将其提交给 Kaggle 进行评估:
In [28]: train_df_filled=fill_null_vals(train_df,'Fare')train_df_filled=fill_null_vals(train_df_filled,'Age')assert len(train_df_filled)==len(train_df)test_df_filled=fill_null_vals(test_df,'Fare')test_df_filled=fill_null_vals(test_df_filled,'Age')assert len(test_df_filled)==len(test_df)
这是对scikit-learn
的调用的实际实现,以通过拟合模型从训练数据中学习,然后在测试数据集上生成预测。 注意,尽管这是样板代码,但出于说明目的,实际调用了特定算法,在这种情况下为DecisionTreeClassifier
。
将输出数据写入具有描述性名称的文件,例如csv/dt_PClass_Sex_Age_Sibsp_Parch_1.csv
和csv/dt_PClass_Sex_Fare_1.csv
。
In [29]:
from sklearn import metrics,svm, tree
for formula_name, formula in formula_map.iteritems():print "name=%s formula=%s" % (formula_name,formula)y_train,X_train = pt.dmatrices('Survived ~ ' + formula, train_df_filled,return_type='dataframe')y_train = np.ravel(y_train)model = tree.DecisionTreeClassifier(criterion='entropy', max_depth=3,min_samples_leaf=5)print "About to fit..."dt_model = model.fit(X_train, y_train)print "Training score:%s" % dt_model.score(X_train,y_train)X_test=pt.dmatrix(formula,test_df_filled)predicted=dt_model.predict(X_test)print "predicted:%s" % predicted[:5]assert len(predicted)==len(test_df)pred_results = pd.Series(predicted,name='Survived')dt_results = pd.concat([test_df['PassengerId'], pred_results],axis=1)dt_results.Survived = dt_results.Survived.astype(int)results_file = 'csv/dt_%s_1.csv' % (formula_name)print "output file: %s\\n" % results_filedt_results.to_csv(results_file,index=False)
前面的代码遵循一个标准配方,其摘要如下:
- 读取训练和测试数据集
- 填写我们希望在两个数据集中考虑的要素的所有缺失值
- 为我们希望在
Patsy
中生成机器学习模型的各种特征组合定义公式 - 对于每个公式,请执行以下步骤:
- 调用
Patsy
为我们的训练特征集和训练标签集(由X_train
和y_train
指定)创建设计矩阵。 - 实例化适当的
scikit-learn
分类器。 在这种情况下,我们使用DecisionTreeClassifier
。 - 通过调用
fit(..)
方法拟合模型。 - 调用
Patsy
可以通过调用patsy.dmatrix(..)
为我们的预测输出创建设计矩阵(X_test
)。 - 预测
X_test
设计矩阵,并将结果保存在预测的变量中。 - 将我们的预测写入输出文件,该文件将提交给 Kaggle。
- 调用
我们将考虑以下监督学习算法:
- 逻辑回归
- 支持向量机
- 决策树
- 随机森林
逻辑回归
在逻辑回归中,我们尝试基于一个或多个输入预测变量来预测分类变量(即离散值相关变量)的结果。
逻辑回归可以被认为等同于应用线性回归,但是适用于离散或分类变量。 但是,在二元逻辑回归(适用于泰坦尼克问题)的情况下,我们试图拟合的函数不是线性函数,因为我们仅试图预测只能采用两个值的结果- 0 和 1。使用线性函数进行回归是没有意义的,因为输出不能取 0 到 1 之间的值。理想情况下,我们需要为二进制值输出的回归建模的某种形式是值的阶跃函数 0 和 1。但是,这样的函数定义不明确且不可微,因此定义了具有更好属性的近似值:逻辑函数。 逻辑函数取 0 到 1 之间的值,但偏向 0 和 1 的极值,可以用作分类变量回归的良好近似值。
Logistic 回归函数的形式定义如下:
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-aCGp5J6I-1681366316327)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00187.jpeg)]
下图很好地说明了为什么逻辑函数适合于二进制逻辑回归:
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-XbL07VjK-1681366316327)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00188.jpeg)]
我们可以看到,随着我们增加参数a
的值,我们可以更接近接受 0 到 1 的值以及我们要建模的阶跃函数。 如果f(x) < 0.5
,则将上述特征的一个简单应用是将输出值设置为 0,否则将其设置为 1。
用于绘制特征的代码包含在plot_logistic.py
中。
注
可以在以下位置找到对逻辑回归的更详细的检查:http://en.wikipedia.org/wiki/Logit 和 http://logisticregressionanalysis.com/86-what-is-logistic-regression。
在对泰坦尼克号问题进行逻辑回归时,我们希望预测一个二进制结果,即乘客是否幸存。
我们修改了样板代码,以使用scikit-learn
的sklearn.linear_model.LogisticRegression
类。
将我们的数据提交给 Kaggle 后,获得了以下结果:
公式 | Kaggle 得分 |
---|---|
C(Pclass) + C(Sex) + Fare |
0.76077 |
C(Pclass) + C(Sex) |
0.76555 |
C(Sex) |
0.76555 |
C(Pclass) + C(Sex) + Age + SibSp + Parch |
0.74641 |
C(Pclass) + C(Sex) + Age + Sibsp + Parch + C(Embarked) |
0.75598 |
可以在run_logistic_regression_titanic.py
文件中找到实现逻辑回归的代码。
支持向量机
支持向量机(SVM)是一种特征强大的监督学习算法,用于分类和回归。 它是一种判别式分类器,它在数据的聚类或分类之间绘制边界,因此可以根据新点所属的聚类对它们进行分类。
SVM 不仅可以找到边界线,还可以找到边界线。 他们还尝试确定两侧边界的边距。 SVM 算法试图找到周围具有最大可能余量的边界。
支持向量是定义边界周围最大边缘的点,请删除这些点,并可能找到较大的边缘。
因此,命名为“支持”,因为它们支持边界线周围的边距。 支持向量很重要。 如下图所示:
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-xHGQVraW-1681366316327)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00189.jpeg)]
注
有关此的更多信息,请参考这里。
为了使用 SVM 算法进行分类,我们指定了以下三个内核之一:线性,多边形和 rbf(也称为径向基函数)。
然后,我们导入支持向量分类器(SVC):
from sklearn import svm
然后,我们实例化一个 SVM 分类器,拟合模型,并预测以下内容:
model = svm.SVC(kernel=kernel)
svm_model = model.fit(X_train, y_train)
X_test = pt.dmatrix(formula, test_df_filled)
. . .
Upon submitting our data to Kaggle, the following results were obtained:
公式 | 核类型 | Kaggle 得分 |
---|---|---|
C(Pclass) + C(Sex) + Fare |
poly |
0.71292 |
C(Pclass) + C(Sex) |
poly |
0.76555 |
C(Sex) |
poly |
0.76555 |
C(Pclass) + C(Sex) + Age + SibSp + Parch |
poly |
0.75598 |
C(Pclass) + C(Sex) + Age + Parch + C(Embarked) |
poly |
0.77512 |
C(Pclass) + C(Sex) + Age + Sibsp + Parch + C(Embarked) |
poly |
0.79426 |
C(Pclass) + C(Sex) + Age + Sibsp + Parch + C(Embarked) |
rbf |
0.7512 |
可以在以下文件中完整查看该代码:run_svm_titanic.py
。
在这里,我们看到内核类型为poly
(多项式)且Pclass
,Sex
,Age
,Sibsp
和Parch
特征组合的 SVM 在提交给 Kaggle 时产生最佳效果。 出人意料的是,似乎登机点(Embarked
)以及乘客是独自旅行还是与家人一起旅行(Sibsp + Parch
)确实对乘客的生存机会产生了重大影响。
后一种影响可能是由于对《铁达尼号》实行了妇女和儿童优先的政策。
决策树
决策树背后的基本思想是使用训练数据集创建决策树以进行预测。
它根据单个要素的值将训练数据集递归地分成子集。 每个拆分对应于决策树中的一个节点。 拆分过程一直进行到每个子集都是纯的为止,也就是说,所有元素都属于一个类。 这总是有效的,除非有重复的训练示例属于不同的类别。 在这种情况下,多数阶级获胜。
最终结果是用于对测试数据集进行预测的规则集。
决策树在模仿人类如何对事物进行分类的过程中对一些二进制选择进行编码,但是通过使用信息标准来决定哪个问题在每个步骤中最有用。
例如,您是否希望确定动物x
是哺乳动物,鱼类还是爬行动物; 在这种情况下,我们会问以下问题:
- Does x have fur?
Yes: x is a mammal
No: Does x have feathers?
Yes: x is a bird
No: Does x have scales?
Yes: Does x have gills?
Yes: x is a fish
No: x is a reptile
No: x is an amphibian
这将生成类似于以下内容的决策树:
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-LRDzRzLT-1681366316328)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00190.jpeg)]
注
请参阅这个链接以获取更多信息:
问题在每个节点上的二进制拆分是决策树算法的本质。 决策树的主要缺点是它们可能过拟合数据。
它们是如此灵活,以至于在深度较大的情况下,它们可以记住输入,当将它们用于看不见的数据分类时,结果很差。
解决此问题的方法是使用多个决策树,这被称为使用集成估计器。 整体估计器的一个示例是随机森林算法,我们将在后面讨论。
要在scikit-learn
中使用决策树,我们导入tree
模块:
from sklearn import tree
We then instantiate an SVM classifier, fit the model, and predict the following:
model = tree.DecisionTreeClassifier(criterion='entropy', max_depth=3,min_samples_leaf=5)
dt_model = model.fit(X_train, y_train)X_test = dt.dmatrix(formula, test_df_filled)
#. . .
将我们的数据提交给 Kaggle 后,可获得以下结果:
公式 | Kaggle 得分 |
---|---|
C(Pclass) + C(Sex) + Fare |
0.77033 |
C(Pclass) + C(Sex) |
0.76555 |
C(Sex) |
0.76555 |
C(Pclass) + C(Sex) + Age + SibSp + Parch |
0.76555 |
C(Pclass) + C(Sex) + Age + Parch + C(Embarked) |
0.78947 |
C(Pclass) + C(Sex) + Age + Sibsp + Parch + C(Embarked) |
0.79426 |
随机森林
随机森林和决策树一样,都是非参数模型的一个例子。 随机森林基于决策树。 决策边界是从数据本身获悉的。 它不必是直线,多项式或径向基函数。 随机森林模型建立在决策树概念的基础上,产生大量决策树或森林。 它从数据中随机抽取一个样本,并识别出一组特征以增长每个决策树。 在决策树集之间比较模型的错误率,以找到产生最强分类模型的特征集。
要在scikit-learn
中使用随机森林,我们导入RandomForestClassifier
模块:
from sklearn import RandomForestClassifier
然后,我们实例化random forest
分类器,拟合模型,并预测以下内容:
model = RandomForestClassifier(n_estimators=num_estimators, random_state=0)
rf_model = model.fit(X_train, y_train)
X_test = dt.dmatrix(formula, test_df_filled)
. . .
将我们的数据提交给 Kaggle(公式:C(Pclass) + C(Sex) + Age + Sibsp + Parch + C(Embarked)
),可获得以下结果:
公式 | Kaggle 得分 |
---|---|
10 | 0.74163 |
100 | 0.76077 |
1000 | 0.76077 |
10000 | 0.77990 |
100000 | 0.77990 |
无监督学习算法
在无监督学习中,我们最关心两个任务:降维和聚类。
降维
降维用于帮助系统地可视化高维数据。 这很有用,因为我们的人脑只能看到三个空间维度(可能是一个时间维度),但是大多数数据集涉及的维度要高得多。
用于降维的典型技术是主成分分析(PCA)。 PCA 涉及使用线性代数技术将高维数据投影到低维空间上。 这不可避免地涉及信息丢失,但是通常通过按正确的维数和维数进行投影,可以使信息丢失最小化。 一种常见的降维技术是在我们的数据中找到解释最大差异(信息的代理)的变量组合,并沿这些维度进行投影。
在无监督学习问题的情况下,我们没有标签集(Y
),因此,我们仅对输入数据X
本身调用fit()
,对于 PCA,我们将调用transform()
而不是transform()
,因为我们正在尝试将数据转换为新的表示形式。
我们将用来演示 USL 的数据集之一是鸢尾花数据集,它可能是所有机器学习中最著名的数据集。
scikit-learn
库提供了一组预打包的数据集,可通过sklearn.datasets
模块使用。 鸢尾花数据集就是其中之一。
鸢尾花数据集由来自三种不同种类的鸢尾花(杂色,山和弗吉尼亚)的 150 个数据样本组成,每种类型有 50 个样本。 数据集包含四个要素/维度:
- 花瓣长度
- 花瓣宽度
- 萼片长度
- 萼片宽度
长度和宽度值以厘米为单位。 可以按以下方式加载:
from sklearn.datasets import load_iris
iris = load_iris()
在我们的无监督学习检查中,我们将专注于如何可视化和聚类数据。
在讨论无监督学习之前,让我们先检查一下鸢尾花数据。 load_iris()
命令返回所谓的束对象,该对象本质上是一个字典,除了包含数据的键之外,还包含键。 因此,我们有以下内容:
In [2]: iris_data.keys()
Out[2]: ['target_names', 'data', 'target', 'DESCR', 'feature_names']
此外,数据本身看起来类似于以下内容:
In [3]: iris_data.data.shape
Out[3]: (150, 4)
这对应于四个特征的 150 个样本。 这四个特征如下所示:
In [4]: print iris_data.feature_names
['sepal length (cm)', 'sepal width (cm)', 'petal length (cm)', 'petal width (cm)']
我们还可以看一下实际数据:
In [9]: print iris_data.data[:2]
[[ 5.1 3.5 1.4 0.2][ 4.9 3\\. 1.4 0.2]]
我们的目标名称(我们要尝试预测的名称)看起来类似于以下内容:
In [10]: print iris_data.target_names['setosa' 'versicolor' 'virginica']
如前所述,鸢尾花特征集对应于五维数据,因此我们无法在色图上对其进行可视化。 我们可以做的一件事是选择两个特征并将它们相互绘制,同时使用颜色区分物种特征。 接下来,我们对所有可能的特征组合进行此操作,一次为一组六个不同的可能性选择两个。 这些组合如下:
- 萼片宽度与花瓣长度
- 萼片宽度与花瓣宽度
- 萼片宽度与花瓣长度
- 萼片长度与花瓣宽度
- 萼片长度与花瓣长度
- 花瓣宽度与花瓣长度
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-x2rA4u2q-1681366316328)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00191.jpeg)]
可以在以下文件中找到用于此目的的代码:display_iris_dimensions.py
。 从前面的图中可以看出,山鸢尾点倾向于自身聚类,而弗吉尼亚点和杂色点之间有一些重叠。 这可能使我们得出结论,后两个物种彼此之间的联系比与山鸢尾物种的联系更紧密。
但是,这些是二维数据切片。 如果我们想更全面地查看数据,并用一些表示的所有萼片和花瓣尺寸的视图,该怎么办?
如果在我们的二维图未显示的尺寸之间存在迄今尚未发现的联系怎么办? 有一种可视化的方法吗? 输入降维。 我们将使用降维来提取萼片和花瓣尺寸的两个组合,以帮助可视化。
我们可以应用降维来做到这一点,如下所示:
In [118]: X, y = iris_data.data, iris_data.targetfrom sklearn.decomposition import PCApca = PCA(n_components=2)pca.fit(X)X_red=pca.transform(X)print "Shape of reduced dataset:%s" % str(X_red.shape)Shape of reduced dataset:(150, 2)
因此,我们看到缩小的数据集现在是二维的。 让我们以二维方式直观地显示数据,如下所示:
In [136]: figsize(8,6)fig=plt.figure()fig.suptitle("Dimensionality reduction on iris data")ax=fig.add_subplot(1,1,1)colors=['red','yellow','magenta']cols=[colors[i] for i in iris_data.target]ax.scatter(X_red[:,0],X[:,1],c=cols)
Out[136]:
<matplotlib.collections.PathCollection at 0x7fde7fae07d0>
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-Rt4CyurZ-1681366316328)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00192.jpeg)]
我们可以检查 PCA 缩减的两个维度的构成,如下所示:
In [57]:
print "Dimension Composition:"
idx=1
for comp in pca.components_:print "Dim %s" % idxprint " + ".join("%.2f x %s" % (value, name)for value, name in zip(comp, iris_data.feature_names))idx += 1Dimension Composition:
Dim 1
0.36 x sepal length (cm) + -0.08 x sepal width (cm) + 0.86 x petal length (cm) + 0.36 x petal width (cm)
Dim 2
-0.66 x sepal length (cm) + -0.73 x sepal width (cm) + 0.18 x petal length (cm) + 0.07 x petal width (cm)
因此,我们可以看到两个减小的尺寸是所有四个萼片和花瓣尺寸的线性组合。
该信息的来源位于这里。
K 均值聚类
聚类背后的想法是根据给定的标准将数据集中的相似点分组在一起,从而在数据中找到聚类。
K 均值算法旨在将一组数据点划分为 K 个聚类,以使每个数据点都属于具有最近均值或质心的聚类。
为了说明 K 均值聚类,我们可以将其应用于通过 PCA 获得的减少鸢尾花数据集,但是在这种情况下,我们不像在监督学习中那样将实际标签传递给fit(..)
方法:
In [142]: from sklearn.cluster import KMeansk_means = KMeans(n_clusters=3, random_state=0)k_means.fit(X_red)y_pred = k_means.predict(X_red)
现在,我们显示聚类数据,如下所示:
In [145]: figsize(8,6)fig=plt.figure()fig.suptitle("K-Means clustering on PCA-reduced iris data, K=3")ax=fig.add_subplot(1,1,1)ax.scatter(X_red[:, 0], X_red[:, 1], c=y_pred);
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-s3zMg70K-1681366316329)(https://gitcode.net/apachecn/apachecn-ds-zh/-/raw/master/docs/master-pandas/img/images_00193.jpeg)]
请注意,我们的 K 均值算法集群并不完全对应于通过 PCA 获得的维度。 可从这里获得源代码。
注
有关scikit-learn
中 K 均值聚类的更多信息,通常可以在以下位置找到:http://scikit-learn.org/stable/auto_examples/cluster/plot_cluster_iris.html 和 http://en.wikipedia.org/wiki/K-means_clustering。
总结
在本章中,我们展开了机器学习的旋风之旅,研究了 Pandas 在特征提取,选择和工程中的作用,以及学习了机器学习中的关键概念,例如监督学习与无监督学习。 我们还简要介绍了两种机器学习方法中的一些关键算法,并使用scikit-learn
包利用这些算法来学习数据并做出预测。 本章并不是要对机器学习进行全面的介绍,而是要说明如何使用 Pandas 来协助机器学习领域的用户。