文章目录
- 一、pandas与建模代码结合
- 二、用patsy创建模型描述
- Patsy公式中的数据转换
- 分类数据与Pastsy
- 三、statsmodels介绍
- 评估线性模型
- 评估时间序列处理
- 四、scikit-learn介绍
一、pandas与建模代码结合
用DataFrame.values属性将DataFrame转换为NumPy数组
import pandas as pd
import numpy as np
data = pd.DataFrame({
'x0': [1, 2, 3, 4, 5],
'x1': [0.01, -0.01, 0.25, -4.1, 0.],
'y': [-1.5, 0., 3.6, 1.3, -2.]})
datax0 | x1 | y | |
0 | 1 | 0.01 | -1.5 |
1 | 2 | -0.01 | 0.0 |
2 | 3 | 0.25 | 3.6 |
3 | 4 | -4.10 | 1.3 |
4 | 5 | 0.00 | -2.0 |
data.columns
'''Index(['x0', 'x1', 'y'], dtype='object')'''
# 用DataFrame.values属性将DataFrame转换为NumPy数组
data.values
'''
array([[ 1. , 0.01, -1.5 ],
[ 2. , -0.01, 0. ],
[ 3. , 0.25, 3.6 ],
[ 4. , -4.1 , 1.3 ],
[ 5. , 0. , -2. ]])
'''向pd.DataFrame()传递含有列名的二维ndarray，可将数组转换为DataFrame
# 向pd.DataFrame()传递含有列名的二维ndarray，可将数组转换为DataFrame
df2 = pd.DataFrame(data.values, columns=['one', 'two', 'three'])
df2one | two | three | |
0 | 1.0 | 0.01 | -1.5 |
1 | 2.0 | -0.01 | 0.0 |
2 | 3.0 | 0.25 | 3.6 |
3 | 4.0 | -4.10 | 1.3 |
4 | 5.0 | 0.00 | -2.0 |
选取DataFrame列的子集转换为ndarray
# 用loc选取列子集
model_cols = ['x0', 'x1']
data.loc[:, model_cols].values
'''
array([[ 1. , 0.01],
[ 2. , -0.01],
[ 3. , 0.25],
[ 4. , -4.1 ],
[ 5. , 0. ]])
'''
# 用iloc选取列子集
data.iloc[:, :-2].values
# 直接选取列名
data[['x0', 'x1']].values
# 用reindex()，传递columns=[]
data.reindex(columns=['x0','x1']).values用pd.Categorical([与行数相对应的每行的分类组成的列表], categories=[不同分类组成的列表])添加分类列
data['category'] = pd.Categorical(['a', 'b', 'a', 'a', 'b'], categories=['a', 'b'])
datax0 | x1 | y | category | |
0 | 1 | 0.01 | -1.5 | a |
1 | 2 | -0.01 | 0.0 | b |
2 | 3 | 0.25 | 3.6 | a |
3 | 4 | -4.10 | 1.3 | a |
4 | 5 | 0.00 | -2.0 | b |
用虚拟变量替换category列，要先创建虚拟变量，在删除category列，最后连接结果
# 用pd.get_dummies(df.category列名, prefix='category列名')创建虚拟变量，prefix设置新列名前缀
dummies = pd.get_dummies(data.category, prefix='category')
# df.drop('列名', axis=1)删除原来的分类列
# join将删除分类列的数据与虚拟变量连接
data_with_dummies = data.drop('category', axis=1).join(dummies)
data_with_dummiesx0 | x1 | y | category_a | category_b | |
0 | 1 | 0.01 | -1.5 | 1 | 0 |
1 | 2 | -0.01 | 0.0 | 0 | 1 |
2 | 3 | 0.25 | 3.6 | 1 | 0 |
3 | 4 | -4.10 | 1.3 | 1 | 0 |
4 | 5 | 0.00 | -2.0 | 0 | 1 |
# the same as above
# pd.get_dummies(df)会将df中不是数值型的列全部转换为虚拟变量，并在结果中字段将数值型列放在虚拟变量前
data_with_dummies2 = pd.get_dummies(data)
data_with_dummies2
# 但是要注意，如果原df中的数据列类似于 【是否已婚，是用1表示，否用0表示】，则需要先将此列转换为字符串，再用上面的方法。涉及到的转换方法见：https://www.pypandas.cn/docs/getting_started/basics.html#数据类型二、用patsy创建模型描述
Patsy是一用于描述统计模型（尤其是线性模型）的Python库，能够很好的支持statsmodels中特定的线性模型。
Patsy的公式是特殊字符串语法，如：y ~ x0 + x1，这里的相加表示为模型创建的设计矩阵。即y、x0、x1表示向量的意思。
data = pd.DataFrame({
'x0': [1, 2, 3, 4, 5],
'x1': [0.01, -0.01, 0.25, -4.1, 0.],
'y': [-1.5, 0., 3.6, 1.3, -2.]})
datax0 | x1 | y | |
0 | 1 | 0.01 | -1.5 |
1 | 2 | -0.01 | 0.0 |
2 | 3 | 0.25 | 3.6 |
3 | 4 | -4.10 | 1.3 |
4 | 5 | 0.00 | -2.0 |
patsy.dmatrices()
# patsy.dmatrices()接收一个公式字符串和一个数据集（可是DataFrame或数组的字典），为线性模型创建设计矩阵
import patsy
y, X = patsy.dmatrices('y ~ x0 + x1', data)
# y, X也可直接替换为a，则a[0]表示y，a[1]表示X
# 默认X中加入了Intercept截距项
# 给模型添加名词列+0，即patsy.dmatrices('y ~ x0 + x1 + 0', data)结果中会剔除Intercept截距项
y
'''
DesignMatrix with shape (5, 1)
y
-1.5
0.0
3.6
1.3
-2.0
Terms:
'y' (column 0)
'''
X
'''
DesignMatrix with shape (5, 3)
Intercept x0 x1
1 1 0.01
1 2 -0.01
1 3 0.25
1 4 -4.10
1 5 0.00
Terms:
'Intercept' (column 0)
'x0' (column 1)
'x1' (column 2)
'''
# Patsy的DesignMatrix实例是带有附加元数据的NumPy ndarray
np.asarray(y)
'''s
array([[-1.5],
[ 0. ],
[ 3.6],
[ 1.3],
[-2. ]])
'''
np.asarray(X)
'''
array([[ 1. , 1. , 0.01],
[ 1. , 2. , -0.01],
[ 1. , 3. , 0.25],
[ 1. , 4. , -4.1 ],
[ 1. , 5. , 0. ]])
'''Patsy对象可直接传递给一些算法（比如numpy.linalg.lstsq），这些算法都会执行普通最小二乘回归OLS
# np.linalg.lstsq(X, y)返回y和X进行OLS的系数的最小二乘估计OLSE、残差平方和、矩阵X的秩、X的奇异值
# 这里，因为不想返回矩阵X的秩、X的奇异值，故用 短下横线_ 表示
coef, resid, _, _ = np.linalg.lstsq(X, y, rcond=None)
coef
'''
array([[ 0.3129],
[-0.0791],
[-0.2655]])
'''
# coef.squeeze()从数组的形状中删除单维条目，即把shape中为1的维度去掉
# X.design_info属性保留了模型的元数据，可以将模型列名重新附加到y与X进行OLS的系数的OLSE
coef = pd.Series(coef.squeeze(), index=X.design_info.column_names)
coef
'''
Intercept 0.312910
x0 -0.079106
x1 -0.265464
dtype: float64
'''Patsy公式中的数据转换
可将Python代码混合到Patsy公式中，这样，在执行公式时，库将尝试在封闭作用域内查找使用的函数：
y, X = patsy.dmatrices('y ~ x0 + np.log(np.abs(x1) + 1)', data)
X
'''
DesignMatrix with shape (5, 3)
Intercept x0 np.log(np.abs(x1) + 1)
1 1 0.00995
1 2 0.00995
1 3 0.22314
1 4 1.62924
1 5 0.00000
Terms:
'Intercept' (column 0)
'x0' (column 1)
'np.log(np.abs(x1) + 1)' (column 2)
'''Patsy内置的函数standardize()实现标准化（均值为0、方差为1）、center()实现居中（减列均值）
y, X = patsy.dmatrices('y ~ standardize(x0) + center(x1)', data)
X
'''
DesignMatrix with shape (5, 3)
Intercept standardize(x0) center(x1)
1 -1.41421 0.78
1 -0.70711 0.76
1 0.00000 1.02
1 0.70711 -3.33
1 1.41421 0.77
Terms:
'Intercept' (column 0)
'standardize(x0)' (column 1)
'center(x1)' (column 2)
'''**patsy.build_design_matrices()**使用原始样本内数据集中保存的信息将变换应用于新的样本外数据上
new_data = pd.DataFrame({
'x0': [6, 7, 8, 9],
'x1': [3.1, -0.5, 0, 2.3],
'y': [1, 2, 3, 4]})
new_X = patsy.build_design_matrices([X.design_info], new_data)
new_X
'''
[DesignMatrix with shape (4, 3)
Intercept standardize(x0) center(x1)
1 2.12132 3.87
1 2.82843 0.27
1 3.53553 0.77
1 4.24264 3.07
Terms:
'Intercept' (column 0)
'standardize(x0)' (column 1)
'center(x1)' (column 2)]
'''# 用I()将数据集中两列按列名相加
y, X = patsy.dmatrices('y ~ I(x0 + x1)', data)
X
'''
DesignMatrix with shape (5, 2)
Intercept I(x0 + x1)
1 1.01
1 1.99
1 3.25
1 -0.10
1 5.00
Terms:
'Intercept' (column 0)
'I(x0 + x1)' (column 1)
'''分类数据与Pastsy
将非数字型数据转换为用于模型设计的矩阵的方式有很多。
在Patsy公式中用非数字名词列时，会默认将其转换为虚拟变量。
data = pd.DataFrame({
'key1': ['a', 'a', 'b', 'b', 'a', 'b', 'a', 'b'],
'key2': [0, 1, 0, 1, 0, 1, 0, 0],
'v1': [1, 2, 3, 4, 5, 6, 7, 8],
'v2': [-1, 0, 2.5, -0.5, 4.0, -1.2, 0.2, -1.7]
})
y, X = patsy.dmatrices('v2 ~ key1', data)
X
# 结果中key1[T.b]表示原key1列中分类b为True，即1
'''
DesignMatrix with shape (8, 2)
Intercept key1[T.b]
1 0
1 0
1 1
1 1
1 0
1 1
1 0
1 1
Terms:
'Intercept' (column 0)
'key1' (column 1)
'''
# 若忽略了模型的截距(+0)，每个类别值的列将被包含在模型的设计矩阵中
y, X = patsy.dmatrices('v2 ~ key1 + 0', data)
X
'''
DesignMatrix with shape (8, 2)
key1[a] key1[b]
1 0
1 0
0 1
0 1
1 0
0 1
1 0
0 1
Terms:
'key1' (columns 0:2)
'''用C()函数可以将（其实时分类数据的）数字型列解释为分类类型
y, X = patsy.dmatrices('v2 ~ C(key2)', data)
X
'''
DesignMatrix with shape (8, 2)
Intercept C(key2)[T.1]
1 0
1 1
1 0
1 1
1 0
1 1
1 0
1 0
Terms:
'Intercept' (column 0)
'C(key2)' (column 1)
'''当在模型中包含多列分类数据时，可以添加交互项key1:key2 （用于方差分析ANOVA 模型）
# 构造包含多列分类数据的数据集
data['key2'] = data['key2'].map({0: 'zero', 1: 'one'})
datakey1 | key2 | v1 | v2 | |
0 | a | zero | 1 | -1.0 |
1 | a | one | 2 | 0.0 |
2 | b | zero | 3 | 2.5 |
3 | b | one | 4 | -0.5 |
4 | a | zero | 5 | 4.0 |
5 | b | one | 6 | -1.2 |
6 | a | zero | 7 | 0.2 |
7 | b | zero | 8 | -1.7 |
y, X = patsy.dmatrices('v2 ~ key1 + key2', data)
X
'''
DesignMatrix with shape (8, 3)
Intercept key1[T.b] key2[T.zero]
1 0 1
1 0 0
1 1 1
1 1 0
1 0 1
1 1 0
1 0 1
1 1 1
Terms:
'Intercept' (column 0)
'key1' (column 1)
'key2' (column 2)
'''# 添加交互项数据列，之后可用于分析判断交互作用对y的影响
y, X = patsy.dmatrices('v2 ~ key1 + key2 + key1:key2', data)
X
'''
DesignMatrix with shape (8, 4)
Intercept key1[T.b] key2[T.zero] key1[T.b]:key2[T.zero]
1 0 1 0
1 0 0 0
1 1 1 1
1 1 0 0
1 0 1 0
1 1 0 0
1 0 1 0
1 1 1 1
Terms:
'Intercept' (column 0)
'key1' (column 1)
'key2' (column 2)
'key1:key2' (column 3)
'''三、statsmodels介绍
statsmodels是一个Python库，用于拟合多种统计模型，执行统计测试即数据探索和可视化；包含更多“经典”频率学派的统计方法，其中一些模型有：
- 线性模型、广义线性模型和鲁棒线性模型
- 线性混合效应模型
- 方差分析ANOVA分析
- 时间序列过程和状态空间模型
- 广义的矩量法
【注意】本节案例学习都是在已知模型参数的情况下模拟估计、预测，实际运用中模型参数是未知的、待估计的，注意区别。
评估线性模型
statsmodels中的线性模型有两个不同的注意接口：基于数组和基于公式。
以下是基于数组的线性回归举例：
import statsmodels.api as sm
import statsmodels.formula.api as smf
# 定义生成指定均值与方差的样本量为size的正态分布样本
def dnorm(mean, variance, size=1):
if isinstance(size, int): # 判断size是否为int型
size = size, # 这里会将size转换为tuple元组，下面的*size会将size的所有值取出
return mean + np.sqrt(variance) * np.random.randn(*size)
# 用于复现 reproducibility
np.random.seed(12345)
N = 100
# np.c_[ , , ]会将数据一列一列的合并起来
# np.r_[ , , ]会将数据一行一行的合并起来
X = np.c_[dnorm(0, 0.4, size=N),
dnorm(0, 0.6, size=N),
dnorm(0, 0.2, size=N)]
eps = dnorm(0, 0.1, size=N) #
beta = [0.1, 0.3, 0.5] # 真实模型的系数beta
y = np.dot(X, beta) + eps # np.dot()可将矩阵与向量相乘
X[:5]
'''
array([[-0.1295, -1.2128, 0.5042],
[ 0.3029, -0.4357, -0.2542],
[-0.3285, -0.0253, 0.1384],
[-0.3515, -0.7196, -0.2582],
[ 1.2433, -0.3738, -0.5226]])
'''
y[:5]
'''array([ 0.4279, -0.6735, -0.0909, -0.4895, -0.1289])'''# sm.add_constant(X) 用于将截距项添加到现有矩阵X
X_model = sm.add_constant(X)
X_model[:5]
'''
array([[ 1. , -0.1295, -1.2128, 0.5042],
[ 1. , 0.3029, -0.4357, -0.2542],
[ 1. , -0.3285, -0.0253, 0.1384],
[ 1. , -0.3515, -0.7196, -0.2582],
[ 1. , 1.2433, -0.3738, -0.5226]])
'''
# sm.OLS(y,X)可拟合一个最小二乘线性回归模型
model = sm.OLS(y, X_model)
# 模型的fit()方法返回一个回归结果对象（包含估计的模型的参数和其他诊断）
results = model.fit()
results.params
'''array([0.0336, 0.1761, 0.2248, 0.5148])'''
# 结果对象的summary()方法可打印模型的诊断细节
print(results.summary())
'''
OLS Regression Results
==============================================================================
Dep. Variable: y R-squared: 0.435
Model: OLS Adj. R-squared: 0.418
Method: Least Squares F-statistic: 24.68
Date: Mon, 28 Jun 2021 Prob (F-statistic): 6.37e-12
Time: 14:54:16 Log-Likelihood: -33.835
No. Observations: 100 AIC: 75.67
Df Residuals: 96 BIC: 86.09
Df Model: 3
Covariance Type: nonrobust
==============================================================================
coef std err t P>|t| [0.025 0.975]
------------------------------------------------------------------------------
const 0.0336 0.035 0.952 0.343 -0.036 0.104
x1 0.1761 0.053 3.320 0.001 0.071 0.281
x2 0.2248 0.046 4.851 0.000 0.133 0.317
x3 0.5148 0.082 6.304 0.000 0.353 0.677
==============================================================================
Omnibus: 4.504 Durbin-Watson: 2.223
Prob(Omnibus): 0.105 Jarque-Bera (JB): 3.957
Skew: 0.475 Prob(JB): 0.138
Kurtosis: 3.222 Cond. No. 2.38
==============================================================================
Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
'''以下是基于公式的线性回归举例：
data = pd.DataFrame(X, columns=['col0', 'col1', 'col2'])
data['y'] = y
data[:5]col0 | col1 | col2 | y | |
0 | -0.129468 | -1.212753 | 0.504225 | 0.427863 |
1 | 0.302910 | -0.435742 | -0.254180 | -0.673480 |
2 | -0.328522 | -0.025302 | 0.138351 | -0.090878 |
3 | -0.351475 | -0.719605 | -0.258215 | -0.489494 |
4 | 1.243269 | -0.373799 | -0.522629 | -0.128941 |
# smf.OLS('Patsy公式字符串', data=DataFrame)可拟合一个最小二乘线性回归模型
# 模型的fit()方法返回一个回归结果对象（包含估计的模型的参数和其他诊断）
results = smf.ols('y ~ col0 + col1 + col2', data=data).fit()
# 结果对象的params属性返回参数beta最小二乘估计OLSE
results.params
'''
Intercept 0.033559
col0 0.176149
col1 0.224826
col2 0.514808
dtype: float64
'''
# 结果对象的tvalues属性返回检验的t值大小
results.tvalues
'''
Intercept 0.952188
col0 3.319754
col1 4.850730
col2 6.303971
dtype: float64
'''
# 结果对象的predict(data)会根据估计的模型参数返回预测值
results.predict(data[:5])
'''
0 -0.002327
1 -0.141904
2 0.041226
3 -0.323070
4 -0.100535
dtype: float64
'''评估时间序列处理
# 自回归结构和噪声noise（一般包括误差、残差...）
init_x = 4
import random
# 定义起始值
values = [init_x, init_x]
N = 1000
b0 = 0.8
b1 = -0.4
noise = dnorm(0, 0.1, N)
for i in range(N):
# 新值为昨天的值*b0+前天的值*b1+噪声
new_x = values[-1] * b0 + values[-2] * b1 + noise[i]
values.append(new_x)
# 上面构建了具有参数为0.8和-0.4的AR（2）结构（两个滞后数）。
# 实际拟合数据时，滞后项的数目时位置的
# 以下设置滞后数 MAXLAGS = 5
MAXLAGS = 5
# sm.tsa.AR()
model = sm.tsa.AR(values)
results = model.fit(MAXLAGS)
# 给出估计的参数
results.params
'''array([-0.0062, 0.7845, -0.4085, -0.0136, 0.015 , 0.0143])'''四、scikit-learn介绍
scikit-learn是使用最广泛的通用Python机器学习库，包含广泛的标准监督和无监督的机器学习方法，包括用于模型选择和评估、数据转换、数据加载和模型持久化的工具。可用于分类、聚类、降维、预测和其他任务。
以下示例为Kaggle比赛中关于泰坦尼克号上生还乘客的经典数据集。
# 载入册数数据与训练数据集
train = pd.read_csv('datasets/titanic/train.csv')
test = pd.read_csv('datasets/titanic/test.csv')
train[:4]PassengerId | Survived | Pclass | Name | Sex | Age | SibSp | Parch | Ticket | Fare | Cabin | Embarked | |
0 | 1 | 0 | 3 | Braund, Mr. Owen Harris | male | 22.0 | 1 | 0 | A/5 21171 | 7.2500 | NaN | S |
1 | 2 | 1 | 1 | Cumings, Mrs. John Bradley (Florence Briggs Th… | female | 38.0 | 1 | 0 | PC 17599 | 71.2833 | C85 | C |
2 | 3 | 1 | 3 | Heikkinen, Miss. Laina | female | 26.0 | 0 | 0 | STON/O2. 3101282 | 7.9250 | NaN | S |
3 | 4 | 1 | 1 | Futrelle, Mrs. Jacques Heath (Lily May Peel) | female | 35.0 | 1 | 0 | 113803 | 53.1000 | C123 | S |
# 查看各列是否包含缺失数据
train.isnull().sum()
'''
PassengerId 0
Survived 0
Pclass 0
Name 0
Sex 0
Age 177
SibSp 0
Parch 0
Ticket 0
Fare 0
Cabin 687
Embarked 2
dtype: int64
'''
test.isnull().sum()
'''
PassengerId 0
Pclass 0
Name 0
Sex 0
Age 86
SibSp 0
Parch 0
Ticket 0
Fare 1
Cabin 327
Embarked 0
dtype: int64
'''在像这样的统计和机器学习的例子中，一个典型的任务是：
根据数据中的特征来预测乘客是否能幸存。将模型拟合到训练数据集上，然后在样本外 测试数据集上进行评估。
impute_value = train['Age'].median()
# 对于Age列，用训练集的 中位数 填充训练集和测试数据集上的缺失值
train['Age'] = train['Age'].fillna(impute_value)
test['Age'] = test['Age'].fillna(impute_value)# 对于Sex列，添加一列IsFemale作为Sex列的编码版本，即将分类数据转换为由0-1组成的数据列
# 将此列转为虚拟变量也可以
train['IsFemale'] = (train['Sex'] == 'female').astype(int)
test['IsFemale'] = (test['Sex'] == 'female').astype(int)# 选取要分析对乘客生还有影响的列变量，并创建Numpy数组
predictors = ['Pclass', 'IsFemale', 'Age']
X_train = train[predictors].values
X_test = test[predictors].values
y_train = train['Survived'].values
X_train[:5]
'''
array([[ 3., 0., 22.],
[ 1., 1., 38.],
[ 3., 1., 26.],
[ 1., 1., 35.],
[ 3., 0., 35.]])
'''
y_train[:5]
'''array([0, 1, 1, 1, 0], dtype=int64)'''# 用scikit-learn的LogisticRegression模型创建一个模型实例
from sklearn.linear_model import LogisticRegression
model = LogisticRegression()
# 用模型的fit()方法将模型拟合到训练数据集上
model.fit(X_train, y_train)
'''LogisticRegression()'''
# 用模型的predict()方法为测试数据集形成预测
y_predict = model.predict(X_test)
y_predict[:10]
'''array([0, 0, 0, 0, 1, 0, 1, 0, 1, 0], dtype=int64)'''至此，如果拥有测试数据集的真实值，可以计算精度百分比或其他一些错误指标：
(y_true == y_predict).mean()实际的模型训练中可能还需要调整参数，并且存在可用于参数调整的交叉验证等避免过度拟合训练数据。
交叉验证通过分隔训练数据集来模拟样本外预测。
# LogisticRegressionCV中CV即为交叉验证，参数cv默认为3，即K折交叉验证中的K为3，切成3份
# 这里传递的10，是Cs，惩罚项的系数，表示网格搜索在模型正则化参数C上的细致度；值越小，正则化越强
from sklearn.linear_model import LogisticRegressionCV
model_cv = LogisticRegressionCV(Cs=10)
model_cv.fit(X_train, y_train)
'''LogisticRegressionCV()'''# 手动交叉验证，用cross_val_score()
from sklearn.model_selection import cross_val_score
# C=10为惩罚项系数
model = LogisticRegression(C=10)
# cv=4为K折交叉验证的K
# scores为准确率
scores = cross_val_score(model, X_train, y_train, cv=4)
scores
'''array([0.7758, 0.7982, 0.7758, 0.7883])'''
scores.mean()
'''0.7845160182604128'''
本文章为转载内容，我们尊重原作者对文章享有的著作权。如有内容错误或侵权问题，欢迎原作者联系我们进行内容更正或删除文章。
