Python建模库介绍 ¶
pandas与建模代码的结合 ¶
介绍两个流行的建模工具包：
import pandas as pd
import numpy as np
使用pandas用于数据载入和数据清洗，之后切换到模型库去建立模型是一个常见的模型开发工作流。 在机器学习中，特征工程是模型开发的重要部分之一。 特征工程是指从原生数据集中提取可用于模型上下文的有效信息的数据转换过程或分析。
pandas和其他分析库的结合点通常是NumPy数组。 要将DataFrame转换为NumPy数组，使用.values属性：
df = 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.]
}
)
print(df)
# x0 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
print(df.columns)
# Index(['x0', 'x1', 'y'], dtype='object')
print(df.values)
# [[ 1. 0.01 -1.5 ]
# [ 2. -0.01 0. ]
# [ 3. 0.25 3.6 ]
# [ 4. -4.1 1.3 ]
# [ 5. 0. -2. ]]
将数组再转换为DataFrame：
df2 = pd.DataFrame(df.values, columns=['one', 'two', 'three']) # 递一个含有列名的二维ndarray
print(df2)
# one 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
.values属性一般在数据是同构化的时候使用——例如，都是数字类型的时候。如果数据是异构化的，结果将是Python对象的 ndarray：
添加一个非数字类型的列。
df3 = df.copy()
df3['category'] = pd.Categorical(['a', 'b', 'a', 'a', 'b'], categories=['a', 'b'])
print(df3)
# x0 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
print(df3.values)
# [[1 0.01 -1.5 'a']
# [2 -0.01 0.0 'b']
# [3 0.25 3.6 'a']
# [4 -4.1 1.3 'a']
# [5 0.0 -2.0 'b']]
通过loc索引和values使用一部分列数据。
model_cols = ['x0', 'x1']
result = df.loc[:, model_cols].values
print(result)
# [[ 1. 0.01]
# [ 2. -0.01]
# [ 3. 0.25]
# [ 4. -4.1 ]
# [ 5. 0. ]]
如果我使用虚拟变量替代df3的category列，先创建虚拟变量，之后删除categroy列，然后连接结果：
dummies = pd.get_dummies(df3.category, prefix='category')
print(dummies)
# category_a category_b
# 0 1 0
# 1 0 1
# 2 1 0
# 3 1 0
# 4 0 1
data_with_dummies = df3.drop('category', axis=1).join(dummies)
print(data_with_dummies)
# x0 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
使用Patsy创建模型描述 ¶
样本的表示形式：
- 在数据挖掘过程中，样本以特征值矩阵X和目标值向量Y的形式表示。
- 容量为
n，有m个特征的样本，其特征值矩阵X由n个维度为m的列向量组成，第j个列向量为样本中第j个个体的特征值向量； - 目标值向量Y的第
j个分量为样本中第j个个体的目标值。
Patsy是一个用于描述统计模型（尤其是线性模型）的Python库。 它使用一种小型基于字符串的"公式语法"。 Patsy能够很好地支持statsmodels中特定的线性模型。 像 y ~ x0 + x1 这种 a + b的语法并不代表将a和b相加，而是代表为模型创建的设计矩阵的术语（terms in the design matrix）。 patsy.dmatrices函数，取一个公式字符串和一个数据集（可以使DataFrame或dict），然后为线性模型产生设计矩阵：
import pandas as pd
import numpy as np
import patsy
from patsy import dmatrices, dmatrix, demo_data
df = 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.]
}
)
print(df)
# x0 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
y, X = patsy.dmatrices('y ~ x0 + x1', df)
print(y)
# [[-1.5]
# [ 0. ]
# [ 3.6]
# [ 1.3]
# [-2. ]]
print(X)
# [[ 1. 1. 0.01]
# [ 1. 2. -0.01]
# [ 1. 3. 0.25]
# [ 1. 4. -4.1 ]
# [ 1. 5. 0. ]]
print(np.asarray(y)) # Patsy的DesignMatrix实例，含有附加元数据的NumPy.ndarray
# [[-1.5]
# [ 0. ]
# [ 3.6]
# [ 1.3]
# [-2. ]]
print(np.asarray(X)) # Patsy的DesignMatrix实例，含有附加元数据的NumPy.ndarray
# [[ 1. 1. 0.01]
# [ 1. 2. -0.01]
# [ 1. 3. 0.25]
# [ 1. 4. -4.1 ]
# [ 1. 5. 0. ]]
上面X输出中的Intercept(最左边一列)是从哪里来的。 这其实是线性模型的一个惯例，比如普通最小二乘回归法（ordinary least squares regression）。 可以去掉这个截距（intercept），通过y ~ x0 + x1 + 0给模型。
y, X = patsy.dmatrices('y ~ x0 + x1 + 0', df)
print(X)
# [[ 1. 0.01]
# [ 2. -0.01]
# [ 3. 0.25]
# [ 4. -4.1 ]
# [ 5. 0. ]]
这种Patsy对象可以直接传入一个算法，比如numpy.linalg.lstsq，来进行普通最小二乘回归的计算
coef, resid, _, _ =np.linalg.lstsq(X, y, rcond=1) # 最小二乘法
print(coef)
# [[ 0.00925424]
# [-0.25485421]]
print(resid)
# [19.72552896]
coef = pd.Series(coef.squeeze(), index=X.design_info.column_names)
print(coef)
# x0 0.009254
# x1 -0.254854
# dtype: float64
Patsy公式中的数据转换 ¶
可以将Python代码混合到你的Patsy公式中，在执行公式时，Patsy库将尝试在封闭作用域中寻找你使用的函数：
y, X = patsy.dmatrices('y ~ x0 + np.log(np.abs(x1) +1)', df)
print(X)
# [[1. 1. 0.00995033]
# [1. 2. 0.00995033]
# [1. 3. 0.22314355]
# [1. 4. 1.62924054]
# [1. 5. 0. ]]
一些常用的变量变换，包括标准化（standardizing (平均值0，方差1）和中心化（减去平均值）。Patsy有内建的函数可以做到这些。
y, X = patsy.dmatrices('y ~ standardize(x0) + center(x1)', df)
print(X)
# [[ 1. -1.41421356 0.78 ]
# [ 1. -0.70710678 0.76 ]
# [ 1. 0. 1.02 ]
# [ 1. 0.70710678 -3.33 ]
# [ 1. 1.41421356 0.77 ]]
作为建模的一部分，我们可能会在一个数据及上训练模型，然后在另一个数据及上评价模型。 当使用中心化或标准化这样的转换时，我们必须注意，必须用模型在新数据集上做预测。 这叫做状态变换（stateful transformations）。 因为我们必须用原本在训练集上得到的平均值和标准差，用在新的数据集上。
new_df = 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_df)
print(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)]
因为加号（+）在Patsy公式的上下文中并不是加法的意思，当想要对数据集中两列按列名相加时，必须将列名封装到特殊的I函数中：
y, X = patsy.dmatrices('y ~ I(x0 + x1)', df)
print(X)
# [[ 1. 1.01]
# [ 1. 1.99]
# [ 1. 3.25]
# [ 1. -0.1 ]
# [ 1. 5. ]]
分类数据Categorical和Patsy ¶
非数值型数据可以通过很多种方式变为一个模型设计矩阵。
当我们在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., -1.2, 0.2, -1.7]
}
)
y, X = patsy.dmatrices('v2 ~ key1', data)
print(y)
# [[-1. ]
# [ 0. ]
# [ 2.5]
# [-0.5]
# [ 4. ]
# [-1.2]
# [ 0.2]
# [-1.7]]
print(X)
# [[1. 0.]
# [1. 0.]
# [1. 1.]
# [1. 1.]
# [1. 0.]
# [1. 1.]
# [1. 0.]
# [1. 1.]]
如果从模型中舍弃截距，每个类型的列会被包含在模型设计矩阵中。
y, X = patsy.dmatrices('v2 ~ key1 + 0', data)
print(X)
# [[1. 0.]
# [1. 0.]
# [0. 1.]
# [0. 1.]
# [1. 0.]
# [0. 1.]
# [1. 0.]
# [0. 1.]]
数值型列可以通过C函数，变为类型列：
y, X = patsy.dmatrices('v2 ~ C(key2)', data)
print(X)
# [[1. 0.]
# [1. 1.]
# [1. 0.]
# [1. 1.]
# [1. 0.]
# [1. 1.]
# [1. 0.]
# [1. 0.]]
当我们在一个模型中使用多个类型术语时，会变得更复杂一些，之前用key1:key2的形式来包含有交集的术语， 这种方法可以用于使用多个术语，例如，一个方法分析模型（analysis of variance (ANOVA) models）：
data['key2'] = data['key2'].map({0: 'zero', 1: 'one'})
print(data)
# key1 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)
print(X)
# [[1. 0. 1.]
# [1. 0. 0.]
# [1. 1. 1.]
# [1. 1. 0.]
# [1. 0. 1.]
# [1. 1. 0.]
# [1. 0. 1.]
# [1. 1. 1.]]
y, X = patsy.dmatrices('v2 ~ key1 + key2 + key1:key2', data)
print(X)
# [[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.]]
statsmodels介绍 ¶
statsmodels是一个Python库，用于拟合多种统计模型，执行统计测试以及数据探索和可视化。 statsmodels包含更多的“经典”频率学派统计方法，而贝叶斯方法和机器学习模型可在其他库中找到。
包含在statsmodels中的一些模型：
- 线性模型，广义线性模型和鲁棒线性模型
- 线性混合效应模型
- 方差分析（ANOVA）方法
- 时间序列过程和状态空间模型
- 广义的矩量法
import pandas as pd
import statsmodels.api as sm
import statsmodels.formula.api as smf
import numpy as np
import random
评估线性模型 ¶
统计模型中有几种线性回归模型，从较基本的（例如，普通最小二乘）到更复杂的（例如，迭代重新加权的最小二乘）。
tatsmodels中的线性模型有两个不同的主要接口，这些接口通过这些API模块导入来访问：
- 基于数组
- 基于公式
下面的例子是调用已知参数beta的模型。在这种情况下，dnorm是用于生成具有特定均值和方差的正态分布数据的辅助函数。
def dnorm (mean, variance, size=1):
if isinstance(size, int):
size = size,
return mean + np.sqrt(variance) * np.random.randn(*size)
np.random.seed(12345)
N = 100
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]
y = np.dot(X, beta) + eps
print(X[:5])
# [[-0.12946849 -1.21275292 0.50422488]
# [ 0.30291036 -0.43574176 -0.25417986]
# [-0.32852189 -0.02530153 0.13835097]
# [-0.35147471 -0.71960511 -0.25821463]
# [ 1.2432688 -0.37379916 -0.52262905]]
print(y[:5])
# [ 0.42786349 -0.67348041 -0.09087764 -0.48949442 -0.12894109]
线性模型通常与我们在Patsy中看到的截距项相匹配。sm.add_constant函数可以将截距列添加到现有矩阵：
X_model = sm.add_constant(X)
print(X_model[:5])
# [[ 1. -0.12946849 -1.21275292 0.50422488]
# [ 1. 0.30291036 -0.43574176 -0.25417986]
# [ 1. -0.32852189 -0.02530153 0.13835097]
# [ 1. -0.35147471 -0.71960511 -0.25821463]
# [ 1. 1.2432688 -0.37379916 -0.52262905]]
sm.OLS类可以拟合一个最小二乘线性回归：
model = sm.OLS(y, X)
模型的fit方法返回一个回归结果对象，该对象包含了估计的模型参数和其他的诊断：
results = model.fit()
print(results.params)
# [0.17826108 0.22303962 0.50095093]
调用summary方法可以打印出一个模型的诊断细节，此处的参数名称已被赋予通用名称x1、x2等：
print(results.summary())
# OLS Regression Results
# =======================================================================================
# Dep. Variable: y R-squared (uncentered): 0.430
# Model: OLS Adj. R-squared (uncentered): 0.413
# Method: Least Squares F-statistic: 24.42
# Date: Sat, 16 Oct 2021 Prob (F-statistic): 7.44e-12
# Time: 14:21:45 Log-Likelihood: -34.305
# No. Observations: 100 AIC: 74.61
# Df Residuals: 97 BIC: 82.42
# Df Model: 3
# Covariance Type: nonrobust
# ==============================================================================
# coef std err t P>|t| [0.025 0.975]
# ------------------------------------------------------------------------------
# x1 0.1783 0.053 3.364 0.001 0.073 0.283
# x2 0.2230 0.046 4.818 0.000 0.131 0.315
# x3 0.5010 0.080 6.237 0.000 0.342 0.660
# ==============================================================================
# Omnibus: 4.662 Durbin-Watson: 0 -0.002327
# Prob(Omnibus): 0.097 Jarque-Bera (JB): 4.098
# Skew: 0.481 Prob(JB): 0.129
# Kurtosis: 3.243 Cond. No. 1.74
# ==============================================================================
#
# Notes:
# [1] R² is computed without centering (uncentered) since the model does not contain a constant.
# [2] Standard Errors assume that the covariance matrix of the errors is correctly specified.
假设所有模型参数都在DataFrame中：
data = pd.DataFrame(X, columns=['col0', 'col1', 'col2'])
data['y'] = y
print(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
现在可以使用statsmodels公式API和Patsy公式字符串。观察statsmodels如何将结果作为带有DataFrame列名称的Series返回。
results = smf.ols('y ~ col0 + col1 + col2', data=data).fit()
print(results.params)
# Intercept 0.033559
# col0 0.176149
# col1 0.224826
# col2 0.514808
# dtype: float64
给定新的样本外数据后，可以根据估计的模型参数计算预测值：
print(results.predict(data[:5]))
# 0 -0.002327
# 1 -0.141904
# 2 0.041226
# 3 -0.323070
# 4 -0.100535
# dtype: float64
评估时间序列处理 ¶
statsmodels中的另一类模型用于时间序列分析。其中包括自回归过程，卡尔曼滤波和其他状态空间模型，以及多变量自回归模型。
下例模拟一些具有自回归结构和噪声的时间序列数据，该数据具有参数为0.8和-0.4的AR（2）结构（两个滞后）。
init_x = 4
values = [init_x, init_x]
N = 1000
b0 = 0.8
b1 = -0.4
noise = dnorm(0, 0.1, N)
for i in range(N):
new_x = values[-1] * b0 + values[-2] * b1 + noise[i]
values.append(new_x)
当拟合一个AR模型时，你可能不知道包含的滞后项的数量，所以可以用更大的滞后数来拟合该模型：
MAXLAGS = 5
model = sm.tsa.AR(values)
results = model.fit(MAXLAGS)
print(results.params)
# NotImplementedError: AR has been removed from statsmodels and replaced with statsmodels.tsa.ar_model.AutoReg.
sikit-learn介绍 ¶
scikit-learn是使用最广泛且最受信任的通用Python机器学习库。 它包含广泛的标准监督的和无监督的机器学习方法，包括用于模型选择和评估、数据转换、数据加载和模型持久化的工具。 这些模型可用于分类、聚类、预测和其他常见任务。 pandas非常适合在模型拟合前处理数据集。
举个例子，用一个Kaggle竞赛的经典数据集，关于泰坦尼克号乘客的生还率。我们用pandas加载测试和训练数据集：
import pandas as pd
from sklearn.linear_model import LogisticRegression
from sklearn.linear_model import LogisticRegressionCV
from sklearn.model_selection import cross_val_score
train = pd.read_csv('../datasets/titanic/train.csv')
test = pd.read_csv('../datasets/titanic/test.csv')
print(train[:4])
# PassengerId Survived Pclass ... Fare Cabin Embarked
# 0 1 0 3 ... 7.2500 NaN S
# 1 2 1 1 ... 71.2833 C85 C
# 2 3 1 3 ... 7.9250 NaN S
# 3 4 1 1 ... 53.1000 C123 S
像statsmodels和scikit-learn通常不能提供缺失数据，因此我们要检查各列，看看是否有包含缺失数据：
print(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
print(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
在像这样的统计和机器学习的例子中，一个典型的任务是根据数据中的特征来预测乘客是否能幸存下来。 将模型拟合到训练数据集上，然后在样本外测试数据集上进行评估。 如果用Age作为预测，但它缺少数据。需要进行缺失数据插补（imputation），并使用训练数据集的中间值填充两个表中的空值：
impute_value = train['Age'].median()
train['Age'] = train['Age'].fillna(impute_value)
test['Age'] = test['Age'].fillna(impute_value)
现在建立模型。 添加一列IsFemale作为’Sex’列的编码版本：
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
print(X_train[:5])
# [[ 3. 0. 22.]
# [ 1. 1. 38.]
# [ 3. 1. 26.]
# [ 1. 1. 35.]
# [ 3. 0. 35.]]
print(y_train[:5])
# [0 1 1 1 0]
使用scikit-learn的LogisticRegression模型创建一个模型实例：
model = LogisticRegression()
与statsmodels类似，使用模型的fit方法在训练数据上拟合模型：
result = model.fit(X_train, y_train)
print(result)
# LogisticRegression()
使用model.predict为测试数据集形成预测：
y_predict = model.predict(X_test)
print(y_predict[:10])
# [0 0 0 0 1 0 1 0 1 0]
实际上，模型训练中经常存在许多附加的复杂层次。 许多模型具有可以调整的参数，并且存在可用于参数调整的交叉验证等技术以避免过度拟合训练数据。 这通常可以在新数据上产生更好的预测性能或稳健性。交叉验证通过分割训练数据来模拟样本外预测。 基于像均方误差之类的模型准确度分数，可以对模型参数执行网格搜索。 一些模型，如逻辑回归，具有内置交叉验证的估计类。 例如，LogisticRegressionCV类可以与一个参数一起使用，该参数表示网格搜索在模型正则化参数C上的细致度：
model_cv = LogisticRegressionCV()
result = model_cv.fit(X_train, y_train)
print(result)
# LogisticRegressionCV()
要手动进行交叉验证，可以使用cross_val_score函数，该函数处理数据拆分过程。 例如，为了用我们的模型与训练数据的四个非重叠分割进行交叉验证，可以这样做：
model = LogisticRegression(C=10)
scores = cross_val_score(model, X_train, y_train, cv=4)
print(scores)
# [0.77578475 0.79820628 0.77578475 0.78828829]
默认评分指标是依赖于模型的，但可以选择明确的评分函数。经过交叉验证的模型需要更长时间的训练，但通常可以产生更好的模型性能。