大数据人工智能常用特征工程与数据预处理Python实践(2)

续上一篇: 《大数据人工智能常用特征工程与数据预处理Python实践(1)》

5. 特征提取

特征提取 ( Feature extraction )是指利用已有的特征计算出一个抽象程度更高的特征集,也指计算得到某个特征的算法。

在通常大部分生产、生活的大数据人工智能的项目中,除去语音和图像等特定专业化场景,由于各种原因,而没有足够的训练数据支撑,我们还无法完全信任算法自动生成的特征,因而基于人工经验的特征工程依然是目前的主流。

往往人工经验这件事不易掌控、层次水平差距较大,加之许多业界的项目由于隐私性的考虑等原因,专家很少会透露底层的入模特征和计算逻辑,使得目前 关于特征工程细节的文章少之又少。

关于提特征的问题,主要是时域和频域特征,用滑动窗口提取特征,比如平均数、方差、过零率等,还有傅里叶变换后的幅度、频率、均值等。

模拟案例,加油价格敏感客户分析

为了便于理解,我们假设一个现实的模拟场景,分析加油价格敏感客户分析,如上图所示,特取特征来自客户的消费记录、优惠券记录,我们看看如何提取特征。特征工程对于业务优化、创新发展的意义?

5.1. 时域特征

时域特征是指,在序列随时间变化的过程中,所具有的与时间相关的一些特征。我们用

n
n



n




来表示一个时间窗口的大小,用

i
i



i




表示第

i
i



i




行数据。

时域特征值通常分为有量纲参数与无量纲参数。含量纲的时域特征,常用的包括最大值、最小值、极差(range)、均值(mean)、中位数(media)、众数(mode)、标准差(standard deviation)、均方根值(root mean square/rms)、均方值(mean square/ms)、k阶中心/原点矩等。

5.1.1. 均值、方差、标准差特征

为了反映一个数据集的离散程度,我们经常采用方差、标准差描述特征。例如描述会员对价格敏感情况,体现在消费金额随价格变动的幅度,假设价格敏感客户,存在消费金额偏离度较小的情况。

(1)均值:



μ
=
1
n

i
=
1
n
x
i
\mu = \frac{1}{n}\sum_{i=1}^{n} x_{i}










i
=
1











n









x







i












(2)方差:



σ
2
=
1
2

i
=
1
n
(
x
i

μ
)
2
\sigma^{2}=\frac{1}{2}\sum_{i=1}^{n}(x_{i}-\mu)^2










i
=
1











n








(

x







i












μ

)






2










(3)标准差:



σ
=
1
n

1

i
=
1
n
(
x
i

μ
)
2
\sigma = \sqrt {\frac{1}{n-1}\sum_{i=1}^{n}(x_{i}-\mu)^{2}}
















i
=
1











n








(

x







i










μ

)







2


















标准差是方差的开方,描述的是样本集合的各个样本点到均值的距离之平均;同时,也可以反映一个数据集的离散程度

疑问:方差和标准差都表示数据的离散程度,那么既然有了方差,为什么还要有标准差呢?

为了和原始数据统一量纲。

5.1.2. 百分比、偏度、峭度

由于标准差仍然是带量纲的,受实际值影响较多,如果用无量纲的相对值将如何呢?(可见6.2.6章节中的图,例如特征“amount_skew”)

(1)离散度百分比



σ

=
σ
μ
\sigma’=\frac{\sigma}{\mu}


(2)峭度



k
4
=
1
n

1

i
=
1
n
(
x
i

μ
σ
)
4
k_{4} = \frac{1}{n-1}\sum_{i=1}^{n}(\frac{x_{i}-\mu}{\sigma})^{4}










i
=
1











n








(







σ





x







i










μ









)







4











峭度因子是表示波形平缓程度的,用于描述变量的分布。正态分布的峭度等于3,峭度小于3时分布的曲线会较“平”,大于3时分布的曲线较“陡”。

(3)偏度



k
3
=
1
n

1

i
=
1
n
(
x
i

μ
σ
)
3
k_{3} = \frac{1}{n-1}\sum_{i=1}^{n}(\frac{x_{i}-\mu}{\sigma})^{3}










i
=
1











n








(







σ





x







i










μ









)







3











偏度因子:偏度也叫偏斜度、偏态。偏度和峭度是有一定的相关性的,峭度因子是四阶中心矩和标准差的四次方的比值;偏度因子是三阶中心矩和标准差的三次方的比值。偏度与峭度相同,描述的是分布。

# 计算特征,标准差、均值、均方差、涨价加油频次
def generateFeature():
client = pymongo.MongoClient('mongodb://study:study@localhost:27017/study')
db = client["study"]
collection = db["sales_feature"]
df = pd.DataFrame(list(collection.find()))
print(df.dtypes)
df = df[(df.id<1100) & (df.fuelle_date > '2019-11-30')]
#计算加油消费金额、加油量、加油间隔的标准差
df_feature = df[['id','price_sensitive','amount','fuelle','fuel_interva']].groupby(['id','price_sensitive']).std()
#计算均值
df_feature1 = df[['id','price_sensitive','amount','fuelle','fuel_interva']].groupby(['id','price_sensitive']).mean()
#修改为均值列名
df_feature1 = df_feature1.rename(columns={
'amount':'amount_mean','fuelle':'fuel_mean','fuel_interva':'fuel_interva_mean'})
df_feature = pd.merge(left=df_feature, right=df_feature1,how="left",on=["id",'price_sensitive'])   #左连接
#df_feature1 = df[['id','amount','fuel_interva']].groupby(['id']).apply(np.std)
#计算涨价提前加油的频次  
df2 = df[(df.id<1100) & (df.fuelle_date > '2019-11-30') & (df.changes<0) & (df.time_before==-1)]
df_feature3 = df2.groupby(['id','price_sensitive'], as_index=False)['time_before'].count()  #/4
# df_feature = pd.merge(left=df_feature, right=df_feature3,how="outer",on="id") #外连接
df_feature = pd.merge(left=df_feature, right=df_feature3,how="left",on=["id",'price_sensitive'])   #左连接
df_feature = df_feature.fillna(0)
df_feature['id'] = df_feature['id'].astype('int')
df_feature['price_sensitive'] = df_feature['price_sensitive'].astype('int')
# 百分比(一阶)
df_feature['amount_per'] = df_feature['amount']/df_feature['amount_mean']
df_feature['fuelle_per'] = df_feature['fuelle']/df_feature['fuel_mean']
df_feature['fuel_interva_per'] = df_feature['fuel_interva']/df_feature['fuel_interva_mean']
# 偏度(三阶)
#df_feature['amount_skew','fuelle_skew','fuel_interva_skew']  
df_feature4 = df[['id','price_sensitive','amount','fuelle','fuel_interva']].groupby(['id','price_sensitive']).skew()
df_feature4 = df_feature4.rename(columns={
'amount':'amount_skew','fuelle':'fuelle_skew','fuel_interva':'fuel_interva_skew'})
df_feature = pd.merge(left=df_feature, right=df_feature4,how="left",on=["id",'price_sensitive'])   #左连接
#峭度/峰度(4阶)
df_feature = pd.concat([df_feature, pd.DataFrame(columns={
'amount_kurt','fuelle_kurt','fuel_interva_kurt'})])
for index,row in df_feature.iterrows():
id = row['id']
df_tmp = df[df['id']==id][['amount','fuelle','fuel_interva']]
k4 = df_tmp.kurt()
df_feature.loc[index:index,('amount_kurt','fuelle_kurt','fuel_interva_kurt')] = k4.tolist()
# 峰度/峭度kurt不能聚合
print(df_feature)
df_feature.to_hdf("Featue_DataV1.h5", key='df', mode='w', complevel=9)
Featue_Data = db["Featue_Data"]
Featue_Data.insert(json.loads(df_feature.T.to_json()).values())
return df_feature

按时域分析,增加了标准差、无量纲标准差(标准差比上均值)、偏度、峭度,为机器学习算法提出更为突出的特征,具体使用详见下一章节。

5.2. 时间特征

对于时间特征,我们能取到年、月、日、时、分、秒、季节、周、节假日、间隔时长等等,Pandas举例如下所示:

def Date_Transformation():
data = pd.DataFrame({
'data_time':pd.date_range('2020-10-01 00:00:00',periods = 12,freq = 'D'),'date':pd.date_range('2019-01-01',periods = 12,freq = 'M')
})
# data:提取日期型和时间型的特征变量
data['year']= data['data_time'].dt.year
data['month'] = data['data_time'].dt.month
data['day'] = data['data_time'].dt.day
data['hour'] = data['data_time'].dt.hour
data['minute'] = data['data_time'].dt.minute
data['second'] = data['data_time'].dt.second    # 秒
data['quarter'] = data['data_time'].dt.quarter  # 季度
data['week'] = data['date'].dt.week             # 周
data['weekday'] = data['date'].dt.weekday       # 周几,0为周一
data['yearmonth'] = data['data_time'].dt.strftime('%Y-%m')
data['halfyear'] = data['data_time'].map(lambda d:'H' if d.month <= 6 else 'H2')
# data:转换为相对时间特征
data['deltaDayToToday'] = (datetime.date.today()-data['date'].dt.date).dt.days  #距离今天的间隔(天数)
data['deltaMonthToToday'] = datetime.date.today().month - data['date'].dt.month #距离今天的间隔(月数)
data['daysOfyear'] = data['date'].map(lambda d:366 if d.is_leap_year else 365)   #一年过去的进度
data['rateOfyear'] = data['date'].dt.dayofyear/data['daysOfyear']
data.head()

5.3. 文本特征提取

略。

5.4. 派生特征

从各个数据源中,抽取、清洗出一些特征,将非结构化数据转成结构化数据的,便于后续训练模型和做预测。通常我们把这些特征称为初级特征,在特征偏少的时候及某些情况下,把一些初级特征直接用到模型中效果未必好,根据这些初级特征生成派生特征能更好的提高模型效果。

生成派生特征是数据挖掘领域最有 创造性
的步骤之一。派生特征使建模过程融入更多人工经验,并充分利用专家经验。构造合适的特征集合是优秀算法工程师必备技能。

派生特征对模型效果指标,比如均方误差,准确率,召回率等提升明显。更为重要的是,专业化的派生特征能提升模型的可解释性。但是,没有“放之四海而皆准”的特征,一个特征可能在一个场景效果良好,换个场景效果就大打折扣,如何派生特征,需要根据业务场景和经验积累,个性化和通用化结合来处理。

5.4.1. 单特征派生特征

(1)将数值转成百分位(标准化)

我们在谈论小孩子身高,体重等身体参数时,比如直接说绝对值“身高60cm”,人们会追问孩子几岁等等。这样,不如说“身高处于同龄孩子身高的95分位”更能体现小孩子的生理优势,因为小孩子比95%的同龄儿童都高。而且正则化基于正态假设,而转成分位数的方式与数值的分布无关。

(2)将数值转成比率

我们的数据集中往往有大量的数值:订单数,消费金额,在线时长等。而这些数值往往和用户的注册时长和生命周期等有关。如果,一个已经注册6个月的用户的订单数和一个刚注册一个月用户的当担数的绝对值比较意义不大,往往是将订单数除以用户生命周期(在线时长)等得到的比率更有意义。

(3)分箱

在某些场景,却需要根据数值特征生成类别特征,比如朴素bayes模型只支持类别特征(离散特征)。分桶,又叫离散化,就是把数值映射到一组区间,而用区间号代替数值。

比如,我们在分析汽车行驶里程时,面对日均里程从几公里到上千公里的数据,依据生活经验,可以分成几个段,50以下为私家车(定义为0),50到150公里为商务车(定义为1),150到250公里为出租车(定义为2),250公里以上为长途客车(定义为3)等。

分箱的优势:

  • 可以减少过拟合的风险,因为分箱相当于对于数据去粗粒度描述
  • 增加稀疏数据的概率,减少计算量,因为0的数据变多了
  • 减少噪声数据的影响,比如一组数据按照0~100均匀分布,当数据中突然出现一个10000的数据,如果不做分箱的化会对Logistic Regression这种模型的训练造成很大影响
  • 方便特征衍生,因为数据离散化后就可以把特征直接相互做内积提升特征维度
  • 离散化后可以提升模型的鲁棒性,比如我们仍以行驶里程为例,比如甲车30公里、乙车60公里,到了第二年甲车变成35公里,乙车变成70公里,这样数据大量发生了改变,理论上我们要更新模型。但是如果数据采用分箱算法,按行驶里程分箱逻辑,则第一年和第二年数据没有变化,模型也不用变。

分箱方法:

  • 等距(等宽)分箱
  • 等频(等权重)分箱
  • 聚类分箱
  • 决策树分箱
  • 有监督分箱,有监督分箱主要有best-ks分箱和卡方分箱。

5.4.2. 组合特征

把两个特征组合在一起生成一个新的特征,表达原有特征无法表达的信息,是一种常用的派生特征的方法。比如股票中的市盈率、身体健康的体重指数(BMI),下面再介绍几个可能用到的组合特征。

在以往的数据分析中,专家通常借助采用把采集数据图形化的方式分析问题,也就是图形对标分析。例如股票市场的5日线、10日线、20日线的分析,采油工程中游梁式抽油机的示功图分析识别工况。下面以抽油机的示功图为例介绍组合特征:

图形数据,

x
x



x




轴为往复运行的位置,

y
y



y




轴为运行过程中所做的功能。

(1)向量余弦

(2)面积

采用向量叉乘计算多边形面积。由图知坐标原点(按P点为原点)与多边形任意相邻的两个顶点构成一个三角形,而三角形的面积可由三个顶点构成的两个平面向量的外积求得。任意多边形的面积公式:

(3)斜率

5.4.3. 组合高关联特征

市场营销领域,有好多高关联特征,比如本月支付额和上月账单非常的接近,全年预算和全年订单强相关。如果我们选择了一组强关联特征中的一个,是否可以忽略掉其他特征呢?答案是不一定。强关联特征的差异可能非常有意义,特征见的差异可以通过差值,比率,关联系数等度量表示。

近似特征的差

电商领域常用的两个几乎相等的特征原始订单价和净订单价,净订单价指的是原始订单价减去用户退货和缺货商品的订单价格。采样了50000个客户的订单,我们计算出原始订单价和净订单价的关联系数为0.993,这是一个相当高的值,在模型中,最好值使用其中一个特征,去掉另一个特征,一般净订单价格会被去掉,因为该特征的获取总有一定的延迟。这50000个用户的原始订单价的平均值是197,而净订单价格的平均值是187。对大多数用户而言,这两个值是相等的。但,那些这两个值不等的一小部分客户,是非常特殊的群体,比如有的用户买的所有商品最终都退货,识别这个群体会非常有意义。因此,我们派生出一个特征,即原始订单价和净订单价的差值。如果两个特征的值几乎总是相等的,那么他们不等的场景就显得非常重要。

6. 特征选择

目前,在很多流行的机器学习材料中,我没有看到给出特征工程和特征选择的系统性的详细论述。分析其主要原因:

一是大部分机器学习算法有标准的推导过程,较易于讲解;

二是在很多实际工程项目问题中,寻找和筛选特征变量并没有普适的方法;

三是特征工程属于信息数据库技术与数学交叉科学,覆盖面广而且复杂。

然而,特征工程和特征选择对于分析结果的影响,往往比之后的机器学习模型的选择更为重要。斯坦福大学教授,Coursera创始人吴恩达(Andrew Ng)老师就曾经表示:“基本上,所谓机器学习应用,就是进行特征工程。”

深度学习谈特征工程较少,对于图像特征一般是不考虑的,直接筛选,卷积池化缩减特征,全连接映射特征。

特征选择是一个很重要的数据预处理过程,是从原始数据中选择对于机器学习而言最好的特征的过程,更正式地说,给定n个特征,搜索其中包括k个特征的子集来改善机器学习的性能,包括特征 增维、降维
。特征选择尝试剔除数据中的噪声。特征选择的方法包括基于统计和基于模型的特征选择。

6.1. 特征分析

分析方法 描述
分布分析 研究数据的分布特征和分布类型,分定量数据、定性数据区分基本统计量。
对比分析 两个互相联系的指标进行比较
统计分析 统计指标对定量数据进行统计描述,常从集中趋势和离中趋势两个方面进行分析
重要程度分析 使用随机森林等算法
正态性检验 利用观测数据判断总体是否服从正态分布的检验称为正态性检验,它是统计判断中一种重要的特殊你和优度假设检验
相关性分析 分析连续变量之间的线性相关程度的强弱
周期性分析 探索某个变量是否随着时间变化而呈现出某种周期性变化的趋势

6.1.1. 相关分析

大数据时代,“我们不再热衷于寻找因果关系,而应该寻找事物之间的相关关系。

“万物皆有联”,是大数据一个最重要的核心思维。

所谓联,这里指的就是事物之间的相互影响、相互制约、相互印证的关系。而事物这种相互影响、相互关联的关系,就叫做相关关系,简称相关性。

数学变量相关关系

相关关系:当一个或几个相互联系的变量取一定的数值时,与之相对应的另一变量的值虽然不确定,但它仍按某种规律在一定的范围内变化。变量间的这种相互关系,称为具有不确定性的相关关系。

按程度分类

⑴完全相关:两个变量之间的关系,一个变量的数量变化由另一个变量的数量变化所惟一确定,即函数关系。

⑵不完全相关:两个变量之间的关系介于不相关和完全相关之间。

⑶不相关:如果两个变量彼此的数量变化互相独立,没有关系。

按方向分类

⑴正相关:两个变量的变化趋势相同,从散点图可以看出各点散布的位置是从左下角到右上角的区域,即一个变量的值由小变大时,另一个变量的值也由小变大。

⑵负相关:两个变量的变化趋势相反,从散点图可以看出各点散布的位置是从左上角到右下角的区域,即一个变量的值由小变大时,另一个变量的值由大变小。

按形式分类

⑴线性相关(直线相关):当相关关系的一个变量变动时,另一个变量也相应地发生均等的变动。

⑵非线性相关(曲线相关):当相关关系的一个变量变动时,另一个变量也相应地发生不均等的变动。

按变量数目分类

⑴单相关:只反映一个自变量和一个因变量的相关关系。

⑵复相关:反映两个及两个以上的自变量同一个因变量的相关关系。

⑶偏相关:当研究因变量与两个或多个自变量相关时,如果把其余的自变量看成不变(即当作常量),只研究因变量与其中一个自变量之间的相关关系,就称为偏相关。

Pearson相关系数 (Pearson Correlation)

皮尔森相关系数是一种最简单的,能帮助理解特征和响应变量之间关系的方法,该方法衡量的是变量之间的线性相关性,结果的取值区间为

[

1

1
]
[ − 1 , 1 ]



[

1

1
]




,-1表示完全的负相关,+1表示完全的正相关,0表示没有线性相关。

ρ
X
,
Y
=
c
o
v
(
X
,
Y
)
σ
X
σ
Y
\rho_{X,Y}=\frac{cov(X,Y)}{\sigma _{X}\sigma_{Y}}











σ







X










σ







Y













c
o
v
(
X
,
Y
)









ρ
X
,
Y
=

i
=
1
n
(
X
i

X

)
(
Y
i

Y

)

i
=
1
n
(
X

X

)
2

i
=
1
n
(
Y

Y

)
2
\rho_{X,Y}=\frac{\sum_{i=1}^{n}(X_{i}-\overline{X})(Y_{i}-\overline{Y})}{\sqrt{\sum_{i=1}^{n}(X-\overline{X})^{2}}\sqrt{\sum_{i=1}^{n}(Y-\overline{Y})^{2}}}

























i
=
1






n









(
X







X







)







2





























i
=
1






n









(
Y







Y







)







2



























i
=
1






n









(

X







i
















X






)
(

Y







i
















Y






)









相关分析与回归分析在实际应用中有密切关系。然而在回归分析中,所关心的是一个随机变量Y对另一个(或一组)随机变量X的依赖关系的函数形式。而在相关分析中 ,所讨论的变量的地位一样,分析侧重于随机变量之间的种种相关特征。例如,以

X
X



X






Y
Y



Y




分别记小学生的数学与语文成绩,感兴趣的是二者的关系如何,而不在于由

X
X



X




去预测

Y
Y



Y




Pearson Correlation速度快、易于计算,经常在我们数据分析员拿到数据(经过清洗和特征提取)之后第一时间所做的分析。

(1)各个特征两两直接相关分析

# Pearson 分析
def analysis_Pearson(df,names):
#names = ['price','fuelle','amount','Payment','vol','changes',
#           'fuel_month','fuel_day','changes_month','changes_day','fuel_interva','time_before']
correlations = df.corr(method='pearson',min_periods=1)  #计算变量之间的相关系数矩阵
correlations.to_excel('dcorr1.xlsx')
print(correlations)
# plot correlation matrix
fig = plt.figure() #调用figure创建一个绘图对象
ax = fig.add_subplot(111)
cax = ax.matshow(correlations,cmap = 'inferno', vmin=-1, vmax=1)  #绘制热力图,从-1到1
fig.colorbar(cax)  #将matshow生成热力图设置为颜色渐变条
ticks = np.arange(0,8,1) #生成0-9,步长为1
ax.set_xticks(ticks)  #生成刻度
ax.set_yticks(ticks)
ax.set_xticklabels(names) #生成x轴标签
ax.set_yticklabels(names)
plt.show()
return correlations

(2)某特征与其他各个特征相关分析

如下图所示,amount特征分别与price、fulle、change、fuel_month、fuel_day、fuel_interva、time_before两两相关,相当于上面热力矩阵图的一行。

def draw_bar(key_name,key_values):
plt.rcParams['font.sans-serif']=['SimHei'] #显示中文标签
plt.rcParams['axes.unicode_minus']=False
# 标准柱状图的值
def autolable(rects):
for rect in rects:
height = rect.get_height()
if height>=0:
plt.text(rect.get_x()+rect.get_width()/2.0 - 0.3,height+0.02,'%.3f'%height)
else:
plt.text(rect.get_x()+rect.get_width()/2.0 - 0.3,height-0.06,'%.3f'%height)
# 如果存在小于0的数值,则画0刻度横向直线
plt.axhline(y=0,color='black')
#归一化
norm = plt.Normalize(-1,1)
norm_values = norm(key_values)
map_vir = cm.get_cmap(name='inferno')
colors = map_vir(norm_values)
fig = plt.figure() #调用figure创建一个绘图对象
plt.subplot(111)
ax = plt.bar(key_name,key_values,width=0.5,color=colors,edgecolor='black') # edgecolor边框颜色
sm = cm.ScalarMappable(cmap=map_vir,norm=norm)  # norm设置最大最小值
sm.set_array([])
plt.colorbar(sm)
autolable(ax)
plt.show()
def draw_Pearson_bar(correlations,names):
values = correlations.iloc[0].values.tolist()
draw_bar(names,values)
return

6.1.2. 随机森林分析

随机森林(Random forest)是一种集成机器学习方法,它利用随机重采样技术bootstrap和节点随机分裂技术构建多棵决策树,通过投票得到最终分类结果。

用随机森林进行特征重要性评估的思想就是看每个特征在随机森林中的每棵树上做了多大的贡献,然后取个平均值,最后比一比特征之间的贡献大小。

对于分类问题,通常采用基尼不纯度或者信息增益,对于回归问题,通常采用的是方差或者最小二乘拟合。这里我们再学习一下基尼指数来评价的方法。



G
i
n
i
Gini



G
i
n
i




指数的计算公式为:



G
i
n
i
(
p
)
=

k
=
1
K
p
k
(
1

p
k
)
=
1


k
=
1
K
p
k
2
Gini(p)=\sum_{k=1}^{K}p_{k}(1-p_{k})=1-\sum_{k=1}^{K}p_{k}^2



G
i
n
i
(
p
)
=









k
=
1











K









p







k









(
1




p







k









)
=









k
=
1











K









p







k





2











说明:



  • p
    k
    p_{k}




    p







    k













    表示选中的样本属于

    k
    k



    k




    类别的概率,则这个样本被分错的概率是

    (
    1

    p
    k
    )
    (1-p_{k})



    (
    1




    p







    k









    )



  • 样本集合中有

    K
    K



    K




    个类别,一个随机选中的样本可以属于这

    k
    k



    k




    个类别中的任意一个,因而对类别就加和。

我们将变量重要性评分用

V
I
M
VIM



V
I
M




来表示,将

G
i
n
i
Gini



G
i
n
i




指数用

G
I
GI



G
I




来表示,假设

m
m



m




个特征

X
1

X
2

X
3

.
.
.
.
.
.
X
c
X_1,X_2,X_3,……X_c




X






1










X






2










X






3









.
.
.
.
.
.

X






c












,现在要计算出每个特征

X
j
X_j




X






j














G
i
n
i
Gini



G
i
n
i




指数评分

V
I
M
j

G
i
n
i

VIM_j(Gini)



V
I

M






j









G
i
n
i





,亦即第

j
j



j




个特征在RF所有决策树中节点分裂不纯度的平均改变量。

def FeatureRandomForestClassifier(df):
col_names = ['amount','id','price','fuelle','changes',
'fuel_month','fuel_day','fuel_interva','time_before']
df = df[col_names].fillna(0)
df_train_x,df_train_y,df_test_x,df_test_y = get_train_test_split(df,split_key='id',y='amount',percent=0.2)
x_train, y_train= df_train_x.values ,df_train_y.values.astype('int')
x_text, y_test = df_test_x.values , df_test_y.values
# n_estimators:森林中树的数量,随机森林中树的数量默认10个树,精度递增显著,但并不是越多越好,加上verbose=True,显示进程使用信息
# n_jobs  整数 可选(默认=1) 适合和预测并行运行的作业数,如果为-1,则将作业数设置为核心数
forest_model = RandomForestClassifier(n_estimators=10, random_state=0, n_jobs=-1)
forest_model.fit(x_train, y_train)
feat_labels = col_names[1:]
# 下面对训练好的随机森林,完成重要性评估
# feature_importances_  可以调取关于特征重要程度
importances = forest_model.feature_importances_
print("重要性:", importances)
x_columns = col_names[1:]
indices = np.argsort(importances)[::-1]
x_columns_indices = []
for f in range(x_train.shape[1]):
# 对于最后需要逆序排序,我认为是做了类似决策树回溯的取值,从叶子收敛
# 到根,根部重要程度高于叶子。
print("%2d) %-*s %f" % (f + 1, 30, feat_labels[indices[f]], importances[indices[f]]))
x_columns_indices.append(feat_labels[indices[f]])
print(x_columns_indices)
print(len(x_columns))
# 筛选变量(选择重要性比较高的变量)
threshold = 0.05
x_selected = x_train[:, importances > threshold]
plt.figure(figsize=(6, 6))
plt.rcParams['font.sans-serif'] = ["SimHei"]
plt.rcParams['axes.unicode_minus'] = False
plt.title("油价敏感集中各个特征的重要程度", fontsize=16)
plt.ylabel("import level", fontsize=12, rotation=90)
num = len(x_columns)
for i in range(num):
plt.bar(i, importances[indices[i]], color='blue', align='center')
plt.xticks(np.arange(num), x_columns_indices, rotation=90, fontsize=12)
plt.tight_layout()
plt.show()

6.1.3. 数据分布分析

分布分析研究特征数据的分布状态和分布类型,分定量数据、定性数据区分基本统计量。是比较常用的数据分析方法,也可以比较快的找到数据规律,对数据有清晰的结构认识。

数据的分布( distribution ),描述了各个值出现的频繁程度。表示分布最常用的方法是直方图( histogram),这种图用于展示各个值出现的频数或概率。频数指的是数据集中一个值出现的次数。

利用观测数据判断总体是否服从正态分布的检验称为正态性检验,它是统计判决中重要的一种特殊的拟合优度假设检验。常用判断方法直方图初判 / QQ图判断 / K-S检验

直方图初判

从左上角到右下角对角线的子图形,就是通过直方图体现各个特征数据分布状态。

使用Pandas工具,实现直方图矩阵代码如下所示:

def draw_scatter_matrix_his(df):
plt.rcParams['font.sans-serif']=['SimHei'] #显示中文标签
plt.rcParams['axes.unicode_minus']=False
#直方图矩阵(Scatterplot Matrix)散点图
pd.plotting.scatter_matrix(df,figsize=(8,8),diagonal='hist',alpha=0.2)
plt.show()

常见的数据分布有:

(1)高斯/正态分布

(2)泊松分布

一个服从泊松分布的随机变量X,表示在具有比率参数(rate parameter)λ的一段固定时间间隔内,事件发生的次数。参数λ告诉你该事件发生的比率。随机变量X的平均值和方差都是λ。

SAMPLE_SIZE = 1000
mu = 3
loc = 1
#泊松分布
interva = stats.poisson.rvs(mu,loc,SAMPLE_SIZE)
plt.hist(interva) #画密度图
plt.show()

6.1.4. 核密度估计

核密度估计(kernel density estimation)是在概率论中用来估计未知的密度函数,属于非参数检验方法之一。通过核密度估计图可以比较直观的看出数据样本本身的分布特征。

# 画数据分布散点图,及核密度估计图
def draw_scatter_matrix(df):
plt.rcParams['font.sans-serif']=['SimHei'] #显示中文标签
plt.rcParams['axes.unicode_minus']=False
#核密度图矩阵(Scatterplot Matrix)散点图
pd.plotting.scatter_matrix(df,figsize=(8,8),diagonal='kde')
plt.show()

6.1.5. 周期性分析

周期性分析是探索某个变量是否随着时间的变化而呈现出某种周期变化趋势。时间序列是指以固定时间为间隔的、由所观察的值组成的序列。根据观测值的不同频率,可将时间序列分成小时、天、星期、月份、季度和年等时间形式的序列。例如汽车加油就有较为明显的周期特性。

时间尺度:

  • 较长:年度周期性趋势、季节性周期性趋势;
  • 较短:月度周期性趋势、周度周期性趋势、天度周期性趋势、小时度周期性趋势。

6.2. 特征选择

当数据预处理完成后,经过特征分析,我们需要选择有实际意义的特征,做为机器学习算法模型的输入进行训练学习。通常来说,我们基于业务领域知识,从两个方面考虑来选择特征:

一是特征是否发散,如果一个特征不发散,比如其方差接近于0,也就是说样本在这个特征上基本上没有差异,变化幅度相对较小,这个特征对于样本的区分并没有什么用;

二是特征与目标的相关性,这点比较通用显而易见,也就是说与目标相关度高的特征,理应优选选择。

根据特征选择的形式又可以将特征选择方法分为3种:

  • Filter:过滤法,按照发散性或者相关性对各个特征进行评分,设定阈值或者待选择阈值的个数,选择特征。
  • Wrapper:包装法,根据目标函数(通常是预测效果评分),每次选择若干特征,或者排除若干特征。
  • Embedded:嵌入法,先使用某些机器学习的算法和模型进行训练,得到各个特征的权值系数,根据系数从大到小选择特征。类似于Filter方法,但是是通过训练来确定特征的优劣。我们使用sklearn中的feature_selection库来进行特征选择。

6.2.1. 业务领域专业知识

我们做数据分析是为了解决业务领域方面的需求和问题,并不是纯粹的数学研究,所以离不开业务领域知识,特征的选择首先要依据业务知识选择特征,这样的分析结果和分析模型将更有意义。

我们的第一步是找到该领域精通业务的专家,让他们给出业务模型和一些建议。比如我们需要解决一个成品油销售发展会员的问题,那么先找到领域专家,向他们咨询哪些因素(特征)会对会员开发、维系产生影响,专家经验分析给出的较大影响的和较小影响的都要纳入到我们的特征候选集中,这就是我们的特征的第一候选集。(这往往会省力,但容易跑偏)

我们通过各种方法构建的特征集,有时候这个特征集合也可能很大,在尝试降维处理之前,我们有必要用特征工程的方法去选择出较重要的特征结合,这些方法不会用到领域知识,而仅仅是统计学的方法。

6.2.2. 方差

方差描述随机变量对于数学期望的偏离程度。这是通过特征本身的方差来筛选特征。

使用方差法,首先计算各个特征的方差,然后根据设定的阈值,选择方差大于阈值的特征。如果特征的方差小于阈值,则代表该特征的发散性太弱。

V
a
r
[
X
]
=
p
(
1

p
)
\mathrm{Var}[X] = p(1 – p)




V
a
r

[
X
]
=


p
(
1



p
)

所以无论接下来的特征工程要做什么,都要优先消除方差为0的特征。

使用feature_selection库的VarianceThreshold类来选择特征的代码如下:

from sklearn.feature_selection import VarianceThreshold
#方差选择法,返回值为特征选择后的数据
#参数threshold为方差的阈值
VarianceThreshold(threshold=3).fit_transform(df)

6.2.3. 基于随机森林算法特征重要性排序

(1)特征选择的步骤

在特征重要性的基础上,特征选择的步骤如下:

  • 计算每个特征的重要性,并按降序排序
  • 确定要剔除的比例,依据特征重要性剔除相应比例的特征,得到一个新的特征集
  • 用新的特征集重复上述过程,直到剩下m个特征(m为提前设定的值)
  • 根据上述代码中得到的各个特征集合特征集对应的袋外误差率,选择袋外误差率最低的特征集

(2)特征重要性的估计方法

特征重要性的估计通常有两种方法:一是使用uniform或者gaussian抽取随机值替换原特征;一是通过permutation的方式将原来的所有N个样本的第i个特征值重新打乱分布,第二种方法更加科学,保证了特征替代值与原特征的分布是近似的。这种方法叫做permutation test ,即在计算第i个特征的重要性的时候,将N 个特征的第i个特征重新洗牌,然后比较D和表现的差异性,如果差异很大,则表明第i个特征是重要的。

基于Python的Sklearn代码见 6.1.2
随机森林内容。

6.2.4. 相关性

import pandas as pd
titanic = pd.read_csv('./data/titanic.csv')
titanic.head()
# 组合特征
titanic['Sex_pclass_combo'] = titanic['Sex']+'_pclass_'+titanic['Pclass'].astype(str)
titanic.Sex_pclass_combo.value_counts()
# onehot编码
Sex_pclass_combo = pd.get_dummies['Sex_pclass_combo'],drop_first = False,prefix = 'onehot'
Sex_pclass_combo.head()

6.2.5. 信息增益

信息增益在文中主要是针对离散属性的数据集,基于信息熵(香农定理)和条件熵来作为特征选择的方法,并给出python实现。

信息增益(Info Gain), 信息增益率(Info Gain Ratio)和 基尼指数(Gini Index)都是在做决策树时试用的分支选择策略, 下面我们试用numpy来实现这几个指标, 以便加深对这些指标的理解。

(1)信息熵:

H
(
X
)
=

i

X
p
i
1
l
o
g
p
i
H(X) = \sum_{i \in X} p_i \frac {1} {log p_i }



H
(
X
)
=









i

X














p






i















l
o
g

p






i












1









(2)交叉熵:



H
p
,
q
(
X
)
=

i

X
p
i
1
l
o
g
q
i
H_{p,q}(X)=\sum_{i \in X} p_i \frac {1}{log q_i}




H







p
,
q









(
X
)
=









i

X














p






i















l
o
g

q






i












1











因为我们已经有了交叉熵的公式, 所以我们其实可以直接算出


H
(
Y

X
=
x
)
H(Y|X=x)



H
(
Y

X
=


x
)


(3)信息增益:

信息增益也称谓相对熵。



D
K
L
(
p


q
)
=
H
(
p
,
q
)

H
(
p
)
D_{KL}\left ( p||q \right )=H\left ( p,q \right )-H\left ( p \right )




D







K
L










(
p


q
)

=


H

(
p
,
q
)




H

(
p
)




理论上我们可以证明

H
(
Y
)
>
=
H
(
Y

X
)
H(Y)>=H(Y|X)



H
(
Y
)
>


H
(
Y

X
)



信息增益率

信息增益率的提出是为了规避信息增益种属性类别数目造成的影响, 它的计算公式是:

G
a
i
n
R
a
t
i
o
(
Y
,
X
)
=
G
a
i
n
(
Y
,
X
)
S
p
l
i
t
I
n
f
o
(
X
)
GainRatio(Y, X) = {Gain(Y, X) \over SplitInfo(X)}



G
a
i
n
R
a
t
i
o
(
Y
,
X
)
=










S
p
l
i
t
I
n
f
o
(
X
)




G
a
i
n
(
Y
,
X
)










其中:

S
p
l
i
t
I
n
f
o
(
X
)
=

i

X
p
i
1
l
o
g
p
i
=
H
(
X
)
SplitInfo(X)= \sum_{i \in X} p_i \frac{1}{log p_i} = H(X)



S
p
l
i
t
I
n
f
o
(
X
)
=









i

X














p






i















l
o
g

p






i












1








=


H
(
X
)

Python实现代码

# 计算信息增益,data为pandas dataframe,yname为条件列,求xname列的概率
def FeatureEntropy(data, xname, yname):
# 信息熵
def entropy(data, att_name):
'''
data: 数据集
att_name: 属性名
'''
levels = data[att_name].unique()
# 信息熵
ent = 0
for lv in levels:
pi = sum(data[att_name]==lv) / data.shape[0]
ent += pi*np.log(pi)
return -ent
# 条件熵
def conditional_entropy(data, xname, yname):
xs = data[xname].unique()
ys = data[yname].unique()
p_x = data[xname].value_counts() / data.shape[0]
ce = 0
for x in xs:
ce += p_x[x]*entropy(data[data[xname]==x], yname)
return ce
# 信息增益
def gain(data, xname, yname):
en = entropy(data, yname)
ce = conditional_entropy(data, xname, yname)
print('yy条件下, xx的信息熵:{}'.format(ce))
return en - ce
def gain_ratio(data, xname, yname):
g = gain(data, xname, yname)
si = entropy(data, xname)
print('信息增益:{}'.format(g))
print('信息熵:{}'.format(si))
return g / si
g_rat = gain_ratio(data,  xname, yname)
print('信息增益率:{}'.format(g_rat))
return g_rat
yy条件下, xx的信息熵:1.8741493088501542
信息增益:0.031168083794194112
信息熵:3.0903546126220864
信息增益率:0.010085601072088227

6.2.6. 互信息

互信息(mutual_info_classif/regression)

互信息是变量间相互依赖性的量度。不同于相关系数,互信息并不局限于实值随机变量,它更加一般且决定着联合分布

p
(
X
,
Y
)
p(X,Y)



p
(
X
,
Y
)




和分解的边缘分布的乘积

p
(
X
)
p
(
Y
)
p(X)p(Y)



p
(
X
)
p
(
Y
)




的相似程度。对于两个随机变量

X
X



X






Y
Y



Y




,如果其联合分布为

p
(
x
,
y
)
p(x,y)



p
(
x
,
y
)




,边缘分布为

p
(
x
)
p(x)



p
(
x
)






p
(
y
)
p(y)



p
(
y
)




,则互信息可以定义为:

I
(
X
;
Y
)
=

x

X

y

Y
p
(
x
,
y
)
l
o
g
p
(
x
,
y
)
p
(
x
)
p
(
y
)
I(X;Y)=\sum_{x \in X} \sum_{y \in Y} p(x,y) log \frac{p(x,y)}{p(x)p(y)}



I
(
X
;
Y
)
=









x

X




















y

Y













p
(
x
,
y
)
l
o
g







p
(
x
)
p
(
y
)




p
(
x
,
y
)









边缘分布(Marginal Distribution)指在概率论和统计学的多维随机变量中,只包含其中部分变量的概率分布。

I
(
X
;
Y
)
=
H
(
X
)

H
(
X

Y
)
I(X;Y)=H(X) – H(X|Y)



I
(
X
;
Y
)
=


H
(
X
)



H
(
X

Y
)

经过推导后,我们可以直观地看到

H
(
X
)
H(X)



H
(
X
)




表示为原随机变量

X
X



X




的信息量,

H
(
X

Y
)
H(X|Y)



H
(
X

Y
)




为知道事实

Y
Y



Y






X
X



X




的信息量,互信息

I
(
X
;
Y
)
I(X;Y)



I
(
X
;
Y
)




则表示为知道事实

Y
Y



Y




后,原来信息量减少了多少。

互信息与pearson相关系数对照:

(1)pearson相关系数

(2)互信息

['id' 'amount' 'fuelle' 'fuel_interva' 'amount_mean' 'fuel_mean'
'fuel_interva_mean' 'time_before' 'amount_per' 'fuelle_per'
'fuel_interva_per' 'amount_skew' 'fuelle_skew' 'fuel_interva_skew'
'amount_kurt' 'fuel_interva_kurt' 'fuelle_kurt']
[7, 8, 9, 11, 12, 13, 14, 16]
[0.01018806 0.19962888 0.34711295 0.28081759 0.09184614 0.10353152
0.14419254 0.5478874  0.54111304 0.61332275 0.25363344 0.53352316
0.55300996 0.42932777 0.494555   0.17967126 0.60378084]
def FeatureMutual_info(df):
df = df.fillna(0)
df_train_x,df_train_y,df_test_x,df_test_y = get_train_test_split(df,split_key='id',y='price_sensitive',percent=0.2)
X_train, y_train= df_train_x, df_train_y
X_test, y_test= df_test_x, df_test_y
sns.heatmap(df.corr(), annot= True, fmt= '.2f')
selector = SelectKBest(mutual_info_classif, k=8).fit(X_train, y_train)  # 保留8个特征
select_feature = selector.get_support(indices=True)
print(select_feature)
select_feature = select_feature.tolist()
Scores= selector.scores_
print(select_feature)
print(Scores)
X_new_train = selector.transform(X_train)
X_new_test = selector.transform(X_test)
clf = RandomForestClassifier(n_estimators=10, random_state=0, n_jobs=-1)
clf.fit(X_new_train, y_train)
get_socre(y_test, clf.predict(X_new_test))
#采用OVO的方式进行逻辑回归函数参数的定义,结果明显好于OVR方式
clf = LogisticRegression(multi_class="multinomial",solver="newton-cg")
clf.fit(X_new_train, y_train)
get_socre(y_test, clf.predict(X_new_test))
plt.show()
return

由于样本总量仅为1000,预测结果:

随机森林:

[+]F1 score : 0.825

[+]acc score: 0.825

逻辑回归:

[+]F1 score : 0.835

[+]acc score: 0.835

6.2.7. 特征组合

为了提高复杂关系的拟合能力,在特征工程中经常会把一阶离散特征进行两两组合,这样就构成了高阶的组合特征。

通过将单独的特征进行组合(相乘或求笛卡尔积)而形成的合成特征。特征组合有助于表示非线性关系。

电商CTR(Click-Through-Rate点击通过率,泛指广告点击率)预估场景下,判断用户是否会点击推荐商品,性别(男女)是一个特征,做分桶处理离散化后的年龄也是一个特征,假设性别特征是女,年龄是20-30之间,发现将二者特征组合是非常强的对用户预测是否点击商品的的特征。20-30岁的女性点击衣服商品的概率非常高,交叉后的非线性的组合特征对CTR预估挖掘非常重要。

(1)目的:

构造更多更好的特征,提升模型精度(例如:地球仪的经纬密度)。

(2)方法:

  • 多个连续变量: 加减乘除运算
  • 多个类别型变量: 所有值交叉组合

合成特征 (synthetic feature)和特征组合(Feature Crosses)不太一样,特征交叉是特征组合的一个子集。合成特征 (synthetic feature)是一种特征,不在输入特征之列,而是从一个或多个输入特征衍生而来。通过标准化或缩放单独创建的特征不属于合成特征。

6.2.8 PCA降维

PCA(Principal Component Analysis),即主成分分析方法,是一种使用最广泛的数据降维算法。PCA的主要思想是将n维特征映射到k维上,这k维是全新的正交特征也被称为主成分,是在原有n维特征的基础上重新构造出来的k维特征。PCA的工作就是从原始的空间中顺序地找一组相互正交的坐标轴,新的坐标轴的选择与数据本身是密切相关的。其中:

第一个新坐标轴选择是原始数据中方差最大的方向,

第二个新坐标轴选取是与第一个坐标轴正交的平面中使得方差最大的,

第三个轴是与第1,2个轴正交的平面中方差最大的。

依次类推,可以得到n个这样的坐标轴。通过这种方式获得的新的坐标轴,我们发现,大部分方差都包含在前面k个坐标轴中,后面的坐标轴所含的方差几乎为0。于是,我们可以忽略余下的坐标轴,只保留前面k个含有绝大部分方差的坐标轴。事实上,这相当于只保留包含绝大部分方差的维度特征,而忽略包含方差几乎为0的特征维度,实现对数据特征的降维处理。

使用Sklearn库用法如下所示:

import numpy as np
from sklearn.decomposition import PCA
X = np.array([[-1, -1], [-2, -1], [-3, -2], [1, 1], [2, 1], [3, 2],[2,3]])
pca = PCA(n_components=2)
NewX = pca.fit_transform(X)
print(X)
print(NewX)
print('各维度的方差:{}'.format(pca.explained_variance_))
print('各维度的方差值占总方差值的比例:{}'.format(pca.explained_variance_ratio_))
pca = PCA(n_components=1)
NewX = pca.fit_transform(X)
print(NewX)
print('各维度的方差值占总方差值的比例:{}'.format(pca.explained_variance_ratio_))

6.3. 特征选择与机器学习、深度学习

6.3.1. 特征工程没有那么重要?

对于数据挖掘和处理类的问题,使用一般的机器学习方法,需要提前做大量的特征工程工作,而且特征工程的好坏会在很大程度上决定最后效果的优劣(也就是常说的一句话:数据和特征决定了机器学习的上限,而模型和算法只是逼近这个上限而已)。

使用深度学习的话,特征工程就没那么重要了,特征只需要做些预处理就可以了,因为它可以自动完成传统机器学习算法中需要特征工程才能实现的任务,特别是在图像和声音数据的处理中更是如此,但模型结构会比较复杂,训练较为麻烦。另一个方面,虽然深度学习让我们可以省去特征工程这一较为繁琐的过程,但也让我们失去了对特征的认识,如特征的重要性等

6.3.2. 如何选择或衡量机器学习、深度学习这两种方法

第一看数据量,比如训练数据量达到百万以上,深度学习的方法会比较有优势。如果样本集不是大样本,那么特征工程加传统的机器学习方法使用起来泛化能力会更好。

第二看是否需要对结果有较强的解释性和可调节性,解释性是说我们能够了解到产生该输出结果的原因,这样我们能够知道特征的重要程度,并在出错时能够对错误原因进行分析。可调节性是指在出错或有特征的增删时,能够方便的对原模型进行修正以满足新的要求。在这一方面,一般的机器学习方法有一定的优势。

各自的优势领域:

深度学习:图像处理,自然语言处理等,因为图像、语言、文本都较难进行特征工程,交给深度学习是一个很好的选择。

机器学习:金融风控,量化分析,推荐系统,广告预测等,因为需要较好的可解释性,会更多的采用传统机器学习方法。

以上的领域,机器学习和深度学习都可以做,但因为各自的特点和要求,因此会有相对优势的偏向。

后续将继续写“# 7. 特征工程对于业务优化、创新发展的意义”等内容。

由于作者水平有限,欢迎讨论。

参考:

[1].《基于Pandas实现皮尔逊相关与余弦相似度在工业大数据分析中的应用实践》
CSDN博客 , 肖永威 , 2020.08

[2].《一篇带你学会相关分析》
知乎,SPSSAU,2019.06

[3].《《大数据时代》读书笔记》
CSDN博客 , 肖永威 , 2016.06

[4].《Python机器学习笔记:随机森林算法》
博客园 ,战争热诚 ,2019.02

[5].《数据分析 第三篇:数据特征分析(分布+帕累托+周期)》
博客园 ,悦光阴 , 2018.08

[6]. 《深度学习如何处理信息实现智慧之信息熵、相对熵、交叉熵等》
CSDN博客, 肖永威 ,2020年01月

[7].《Python Matplotlib绘制渐变色柱状图(bar)并加边框和配置渐变颜色条(colorbar)》
CSDN博客 ,肖永威 ,2020.09

[8].《机器学习中,有哪些特征选择的工程方法?》
知乎 , 城东, 2016年7月

[9].《信息熵信息增益率基尼系数原理和pandas实战》
DataScience ,2018.09

[10].《机器学习+特征工程vs深度学习—如何选择》
CSDN博客 , 海晨威 ,2018年10月

[11].《信息论(1)——熵、互信息、相对熵》
知乎 , 徐光宁 ,2018年5月

[12].《浅谈常见的特征选择方法》
博客园 , 那少年和狗 ,2020年4月

[13].《信号下的基本时域频域特征(上)》
简书 , Cingti_Yr,2019年06月

[14].《如何生成派生特征》
知乎 ,呼广跃 , 2018年7月

[15].《Python科学计算初探——余弦相似度》
CSDN博客 , 肖永威 ,2018.04

[16].《已知多边形各点平面坐标,Python实现求面积》
CSDN博客 ,肖永威 ,2018年4月

[17].《特征组合&特征交叉 (Feature Crosses)》
segmentfault 思否, kugua233 ,2018.05

[18].《python随机模块random的22种函数》
CSDN博客, Elenstone ,2020年5月

[19].《特征工程之分箱–卡方分箱》
博客园 , 返回主页少年阿斌 2019年3月

[ 20.《特征工程之特征分箱(决策树分箱、卡方分箱、bestks以及评价标准WOE和IV)》
CSDN博客 , Donreen ,2020年9月

[21].《数据降维-PCA主成分分析》
博客园 ,Curry秀 , 2019年12月

肖永威的专栏
我还没有学会写个人说明!
上一篇

以色列推出凝胶燃料火箭发动机

下一篇

外媒:导剪版《正义联盟》首集或将于明年3月上线

你也可能喜欢

评论已经被关闭。

插入图片