kaggle数据分析学习记录

学习记录kaggle比赛中的房价预测案例 转载自：源地址

一 导入库

import pandas as pd
import numpy as np
import seaborn as sns
from scipy import stats
from scipy.stats import skew
from scipy.stats import norm
import matplotlib.pyplot as plt
from sklearn.preprocessing import StandardScaler
from sklearn.manifold import TSNE
from sklearn.cluster import KMeans
from sklearn.decomposition import PCA
from sklearn.preprocessing import StandardScaler
from sklearn.linear_model import Ridge
from sklearn.model_selection import cross_val_score
from sklearn.linear_model import Lasso,LassoCV
from sklearn.linear_model import ElasticNet,ElasticNetCV

二 查看数据的大体构造

train_df = pd.read_csv(r"D:\project\housepredict\house-prices-advanced-regression-techniques\train.csv")
test_df = pd.read_csv(r"D:\project\housepredict\house-prices-advanced-regression-techniques\test.csv")
print(train_df.info()) # dtypes: float64(3), int64(35), object(43)
print(test_df.info())

可见，我们有1461的训练数据和1460的测试数据。对训练集来说，数据的特征列有81个，其中38个是数值类型的，43个类别类型。然后把两个数据集放在一起,这样做除了方便之后可以同时对训练数据集和测试数据集进行数据清洗和特征工程.

all_df = pd.concat([train_df,test_df],axis=0,ignore_index=True)

三 处理缺失数据

对于缺失数据的处理，通常会有以下几种做法

1,如果缺失的数据过多，可以考虑删除该列特征
2.用平均值、中值、分位数、众数、随机值等替代。但是效果一般，因为等于人为增加了噪声
用插值法进行拟合
3.把变量映射到高维空间。比如性别，有男、女、缺失三种情况，则映射成3个变量：是否男、是否女、是否缺失。缺点就是计算量会加大。

missing = all_df.isnull().sum() # 计算每个特征的缺失值个数
missing = missing[missing>0] # 选择index>0
types = all_df[missing.index].dtypes
percent = (all_df[missing.index].isnull().sum()/all_df[missing.index].isnull().count()) #计算百分比，sum为为0的个数，count为所有的个数
missing_data = pd.concat([missing,percent,types],axis=1,keys=["Total","Percent","Types"])
missing_data.sort_values("Total",ascending=False,inplace=True) # 排序
missing_data.plot.bar() # 条状图
plt.show()

结果显示，前6列大于了15%的缺失率，其余主要是 BsmtX 和 GarageX 两大类。

四 数据分析（重点）

偏度和峰度分析

print(train_df['SalePrice'].skew())
print(train_df['SalePrice'].kurt())

相关性

计算得出前十个相关性最高的属性

corrmat = train_df.corr() # 计算每列就是每个属性之间的相关性
# print(corrmat)
k = 10
cols = corrmat.nlargest(k,"SalePrice")['SalePrice'].index
cm = np.corrcoef(train_df[cols].values.T)# 计算每行就是每个属性之间的相关性 需要转置
hm = sns.heatmap(cm, cbar=True, annot=True, square=True, fmt='.2f', annot_kws={'size': 10}, yticklabels=cols.values, xticklabels=cols.values)
plt.show()

## 同时是相关性列，也是缺失数据的
print(missing_data.loc[missing_data.index.intersection(cols)])# intersection为并集的

去除缺失数值多的特征，思路如下：根据热点图可知missing_data中的total>1的特征影响较小，所以采取直接去除的方式进行，保留下total<=1的特征

all_df = all_df.drop((missing_data[missing_data['Total'] > 1]).index,1)
print(all_df.isnull().sum().max())#检查缺失特征的情况

接下来采用正态概率图来进一步分析。QQ图主要是直观的表示观测与预测值之间的差异，偏离线越大，则两个数据集来自具有不同分布的群体的结论的证据就越大。也能通过箱线图或者其它图形作进一步判断。

#histogram and normal probability plot
sns.distplot(train_df['SalePrice'], fit=norm);#纵坐标是频率/组距，横坐标是数值范围
fig = plt.figure()
stats.probplot(train_df['SalePrice'], plot=plt)
plt.show()

结果可见，所要预测的SalePrice的分布与正态分布还是有差距的，为了加强数据的正态分布特征，可以通过log的方式进行。
下一步就是区别开

quantitative = [f for f in all_df.columns if all_df.dtypes[f] != 'object']
qualitative = [f for f in all_df.columns if all_df.dtypes[f] == 'object']
print("quantitative: {}, qualitative: {}" .format (len(quantitative),len(qualitative)))

quantitative 代表int,float数值类型的定量数据，qualitative 代表非数值类型的定性数据。
然后画出quantitative 的正态分布图进行观察。

f = pd.melt(all_df, value_vars=quantitative)#按照quantitative改变数据
g = sns.FacetGrid(f, col="variable",  col_wrap=2, sharex=False, sharey=False)#准备多个图
g = g.map(sns.distplot, "value")#将这些图全部应用某种类型的图，“value”代表某列
plt.show()

根据图形判断，类似正态分布的特征有LotArea,BsmtUnfSF,1stFlrSF,TotalBsmtSF,KitchenAbvGr。同样地，可以通过log操作了提升质量。
那么，我们可以通过偏度来看看这些特征的分布情况。偏度值就代表更多的值集中在什么地方，如果为0，就说明是标准高斯分布，如果为正，说明集中在右侧，值越大，就越往右侧偏，值越小就越往左侧偏；

print(all_df[quantitative].apply(lambda x: skew(x.dropna())).sort_values(ascending=False))#dropna代表删除缺失值

接下来，进一步对这些定量数据进行分析。

定性数据分析

主要采用方差分析或变方分析（Analysis of variance，简称 ANOVA）为统计模型

train = all_df.loc[train_df.index] #由于已经舍弃掉部分缺失数据较多的特征，训练数据也会发生改变
train['SalePrice'] = train_df.SalePrice

def anova(frame):
    anv = pd.DataFrame()
    anv['feature'] = qualitative
    pvals = []
    for c in qualitative:
        samples = []
        for cls in frame[c].unique():
            s = frame[frame[c] == cls]['SalePrice'].values
            samples.append(s)
        pval = stats.f_oneway(*samples)[1]
        pvals.append(pval)
    anv['pval'] = pvals
    return anv.sort_values('pval')

a = anova(train)
a['disparity'] = np.log(1./a['pval'].values)#由于pval值很小，需要进行log处理；disparity越大，差异性越大
sns.barplot(data=a, x='feature', y='disparity')
x=plt.xticks(rotation=90)
plt.show()

其中，f_oneway代表进行F检验。F检验检测样本数据之间的差异性，以此证明结果不同是由数据差异性决定还是由不同方法决定。p-value代表置信度，只有p-value<临界值（F值>f临界值），（即p-value越小，数据差异越大），才能说明数据（特征）之间具有显著差异，证明方法的有效性。因此画图时需要进行处理。
下面对这些定性变量进行下处理，对齐进行数值编码，让他转换为定性的列。

def encode(frame, feature):
    ordering = pd.DataFrame()
    ordering['val'] = frame[feature].unique()
    ordering.index = ordering.val # 将索引0123...转换为同样的内容值
    ordering['spmean'] = frame[[feature, 'SalePrice']].groupby(feature).mean()['SalePrice']#mean()会按照分组的情况进行各组的均值计算，若无则按照各个特征进行计算SalePrice的值，最后为单独抽取某列
    ordering = ordering.sort_values('spmean')#以均值来确定各个等级，值越低等级越低
    ordering['ordering'] = range(1, ordering.shape[0] + 1)# 将object对象转换成数值等级
    ordering = ordering['ordering'].to_dict()# 转换为词典

    for cat, o in ordering.items():
        frame.loc[frame[feature] == cat, feature + '_E'] = o# 将原来的表中进行替换

qual_encoded = []
for q in qualitative:
    # print(q)
    encode(train, q)
    qual_encoded.append(q + '_E')
print(qual_encoded)

# 选出了包含缺失数据的行，处理一下
missing_data = all_df.isnull().sum()
missing_data = missing_data[missing_data>0]
ids = all_df[missing_data.index].isnull()#将数据判断为false和true
# index (0), columns (1)
all_df.loc[ids[ids.any(axis=1)].index][missing_data.index]#判断某一行是否全为空，其实any函数非常简单：判断一个tuple或者list是否全为空，0，False。如果全为空，0，False，则返回False；如果不全为空，则返回True。包含缺失数据的行

```python
def spearman(frame, features):
    spr = pd.DataFrame()
    spr['feature'] = features
    # Signature: a.corr(other, method='pearson', min_periods=None)
    # Docstring:
    # Compute correlation with `other` Series, excluding missing values
    # 计算特征和 SalePrice的 斯皮尔曼相关系数
    spr['spearman'] = [frame[f].corr(frame['SalePrice'], 'spearman') for f in features]
    spr = spr.sort_values('spearman')
    plt.figure(figsize=(6, 0.25 * len(features)))  # width, height
    sns.barplot(data=spr, y='feature', x='spearman', orient='h')
    plt.show()
features = quantitative + qual_encoded
spearman(train, features)

皮尔逊相关系数越高，变量之间的联系更紧密。得出哪几个特征的影响比较大。
然后，为了避免共线性的问题，也需要分析两两特征之间的关系。

plt.figure(1)
corr = train[quantitative+['SalePrice']].corr()#定量特征之间的相关性
sns.heatmap(corr)
plt.figure(2)
corr = train[qual_encoded+['SalePrice']].corr()#定性特征之间的相关性
sns.heatmap(corr)
plt.figure(3)
# [31,27]
corr = pd.DataFrame(np.zeros([len(quantitative)+1, len(qual_encoded)+1]), index=quantitative+['SalePrice'], columns=qual_encoded+['SalePrice'])#定性特征与定量特征之间的相关性
for q1 in quantitative+['SalePrice']:
    for q2 in qual_encoded+['SalePrice']:
        corr.loc[q1, q2] = train[q1].corr(train[q2])#先是indexq1,再是列名q1，决定该位置
sns.heatmap(corr)
plt.show()

三部分：定量特征之间的热点图，定性特征之间的热点图，定性特征与定量特征之间的相关性来进行综合判断

def pairplot(x, y, **kwargs):
    ax = plt.gca()
    ts = pd.DataFrame({'time': x, 'val': y})
    ts = ts.groupby('time').mean()
    ts.plot(ax=ax)
    plt.xticks(rotation=90)

f = pd.melt(train, id_vars=['SalePrice'], value_vars=quantitative + qual_encoded)
g = sns.FacetGrid(f, col="variable", col_wrap=2, sharex=False, sharey=False, size=5)
g = g.map(pairplot, "value", "SalePrice")
plt.show()

对比不同的走势图，可以观察哪些特征与SalePrice有着很大的关系。
对上述自变量特征分析完之后，需要对目标SalePrice特征进行分析。

SalePrice分析

a = train['SalePrice']
a.plot.hist()

利用频数直方图进行观察

变量聚类

通过找到相似数据的特征，实现数据的分类。利用PCA、TSNE进行数据降维，利用Kmean进行聚类。

features = quantitative + qual_encoded
model = TSNE(n_components=2, random_state=0, perplexity=50)#数据降维
X = train[features].fillna(0.).values#每行为一个特征数据,因此降维相当于找到特征
tsne = model.fit_transform(X)#训练模型，并将数据进行转换

std = StandardScaler()
s = std.fit_transform(X)
pca = PCA(n_components=30)#降维
pca.fit(s)#训练模型
pc = pca.transform(s)#使用模型降维
kmeans = KMeans(n_clusters=5)#聚类
kmeans.fit(pc)

fr = pd.DataFrame({'tsne1': tsne[:,0], 'tsne2': tsne[:, 1], 'cluster': kmeans.labels_})
sns.lmplot(data=fr, x='tsne1', y='tsne2', hue='cluster', fit_reg=False)
print(np.sum(pca.explained_variance_ratio_))#除了这些输入参数外，有两个PCA类的成员值得关注。第一个是explained_variance_，它代表降维后的各主成分的方差值。方差值越大，则说明越是重要的主成分。第二个是explained_variance_ratio_，它代表降维后的各主成分的方差值占总方差值的比例，这个比例越大，则越是重要的主成分。

30个成分能覆盖83%的方差，整体看来，这种聚类方法不太好

三 模型建立

缺失数据处理

missing = all_df.isnull().sum()
missing.sort_values(inplace=True,ascending=False)
missing = missing[missing > 0]
all_df = all_df.drop(missing[missing>1].index,1)#缺失值大于1的都删除掉

处理log项

根据画图结果，使得特征更加符合高斯分布

logfeatures = ['GrLivArea','1stFlrSF','2ndFlrSF','TotalBsmtSF','LotArea','KitchenAbvGr','GarageArea']
for logfeature in logfeatures:
    all_df[logfeature] = np.log1p(all_df[logfeature].values)

处理Boolean变量

all_df['HasBasement'] = all_df['TotalBsmtSF'].apply(lambda x: 1 if x > 0 else 0)
all_df['HasGarage'] = all_df['GarageArea'].apply(lambda x: 1 if x > 0 else 0)
all_df['Has2ndFloor'] = all_df['2ndFlrSF'].apply(lambda x: 1 if x > 0 else 0)
all_df['HasWoodDeck'] = all_df['WoodDeckSF'].apply(lambda x: 1 if x > 0 else 0)
all_df['HasPorch'] = all_df['OpenPorchSF'].apply(lambda x: 1 if x > 0 else 0)
all_df['HasPool'] = all_df['PoolArea'].apply(lambda x: 1 if x > 0 else 0)
all_df['IsNew'] = all_df['YearBuilt'].apply(lambda x: 1 if x > 2000 else 0)# 数据大小没有实际意义
quantitative = [f for f in all_df.columns if all_df.dtypes[f] != 'object'] #需要重新定义定量和定性
qualitative  = [f for f in all_df.columns if all_df.dtypes[f] == 'object']

one-hot编码

all_dummy_df = pd.get_dummies(all_df)#只对str进行编码

对于数值变量进行标准化

all_dummy_df.isnull().sum().sum()#检查还有没有缺失值的序列

mean_cols = all_dummy_df.mean()#计算每列的均值
all_dummy_df = all_dummy_df.fillna(mean_cols)#均值填充

X = all_dummy_df[quantitative]
std = StandardScaler()
s = std.fit_transform(X)#对定量数据进行标准化处理
all_dummy_df[quantitative] = s
dummy_train_df = all_dummy_df.loc[train_df.index]#切开训练和测试集
dummy_test_df = all_dummy_df.loc[test_df.index]
y_train = np.log(train_df.SalePrice)#对数处理

模型预测

首先利用岭回归模型查看RMSE变化

def rmse_cv(model):
    rmse= np.sqrt(-cross_val_score(model, dummy_train_df, y_train.values, scoring="neg_mean_squared_error", cv = 5))#计算模型的均方误差RMSE，不是测试集
    return(rmse)
alphas = np.logspace(-3, 2, 50)
cv_ridge = []
coefs = []
for alpha in alphas:#alphas寻找一个合适的超参
    model = Ridge(alpha = alpha)
    model.fit(dummy_train_df,y_train)
    cv_ridge.append(rmse_cv(model).mean())#五个取平均数
    coefs.append(model.coef_)#得到每个特征的系数

cv_ridge = pd.Series(cv_ridge, index = alphas)
cv_ridge.plot(title = "Validation - Just Do It")
plt.xlabel("alpha")
plt.ylabel("rmse")
plt.show()

再查看岭回归模型的各个特征的影响

# 岭迹图 是alphas与相关特征系数之间的图 随着k增大，共线性的影响将越来越小。在不断增大惩罚函数系数的过程中，画下估计参数技术分享（k）的变化情况，即为岭迹。通过岭迹的形状来判断我们是否要剔除掉该参数（例如：岭迹波动很大，说明该变量参数有共线性）。
# matplotlib.rcParams['figure.figsize'] = (12.0, 12.0)
ax = plt.gca()

# ax.set_color_cycle(['b', 'r', 'g', 'c', 'k', 'y', 'm'])

ax.plot(alphas, coefs)# 每一条直线为一个特征，为每个数组同一列
ax.set_xscale('log')
ax.set_xlim(ax.get_xlim()[::-1])  # reverse axis
plt.xlabel('alpha')
plt.ylabel('weights')
plt.title('Ridge coefficients as a function of the regularization')
plt.axis('tight')
plt.show()

再利用Lasso能针对上面特征太多的问题，来选择一部分重要的特征。两种方法都是加了正则化的损失函数模型，不同点是，前者是L1，后者是L2

alphas = np.logspace(-4, -2, 100)
cv_lasso = []
coefs = []
for alpha in alphas:
    model = Lasso(alpha = alpha,max_iter=5000)
    model.fit(dummy_train_df,y_train)
    cv_lasso.append(rmse_cv(model).mean())
    coefs.append(model.coef_)

cv_lasso = pd.Series(cv_lasso, index = alphas)
cv_lasso.plot(title = "Validation - Just Do It")
plt.xlabel("alpha")
plt.ylabel("rmse")
print(cv_lasso.min(), cv_lasso.argmin())#选择产生最低误差的alpha进行训练，选择了最小的值以及index

查看RMSE的变化，然后类似的观察各个特征的重要性。

model = Lasso(alpha=0.00058, copy_X=True, fit_intercept=True, max_iter=1000,
   normalize=False, positive=False, precompute=False, random_state=None,
   selection='cyclic', tol=0.0001, warm_start=False)
model.fit(dummy_train_df,y_train)#直接fit就可以训练模型
coef = pd.Series(model.coef_, index = dummy_train_df.columns)
print("Lasso picked " + str(sum(coef != 0)) + " variables and eliminated the other " +  str(sum(coef == 0)) + " variables")#只要不靠近0代表有关系
imp_coef = pd.concat([coef.sort_values().head(10),
                     coef.sort_values().tail(10)])#只显示20个
matplotlib.rcParams['figure.figsize'] = (8.0, 10.0)
imp_coef.plot(kind = "barh")
plt.title("Coefficients in the Lasso Model")
plt.show()

利用Lasso进行特征选择，选择了“Lasso picked 85 variables and eliminated the other 141 variables”。
因此，结合以上两点，既要解决共线问题，也要从中选择一些重要的特征，使用Elastic Net。

elastic=ElasticNetCV(alphas=[0.001, 0.05, 0.1, 0.3, 1, 3, 5, 10, 15, 30, 50, 75],
       copy_X=True, cv=5, eps=0.001, fit_intercept=True,
       l1_ratio=[0.1, 0.5, 0.7, 0.9, 0.95, 0.99, 1], max_iter=5000,
       n_alphas=100, n_jobs=1, normalize=False, positive=False,
       precompute='auto', random_state=None, selection='cyclic',
       tol=0.0001, verbose=0)
elastic = elastic.fit(dummy_train_df, y_train)
print(rmse_cv(elastic).mean())#0.1290859144132534

模型训练拟合的效果蛮一般的。继续。可以按照不同的特征预处理，依次使用LASSO，岭回归以及Elastic Net或者其它函数进行预测,或者是使用多种的加权求和作为最终的结果，利用predict对test进行预测并提交结果

submission_df = pd.DataFrame(data= {'Id' : test_df.Id, 'SalePrice': np.exp(y_final)})
submission_df.to_csv("bag-4.csv",index=False) # 取消index的存储