Kaggle:Predict House Price（房价预测）

项目背景

目的：基于所提供的数据集，对房价进行预测

#首先，导入一些所需要的工具包
import numpy as np 
import pandas as pd 
%matplotlib inline
import matplotlib.pyplot as plt  
import seaborn as sns
sns.set_style('darkgrid')
import warnings
def ignore_warn(*args, **kwargs): 
    pass
warnings.warn = ignore_warn #忽略主要来自sklearn & seaborn 的警告
from scipy import stats
from scipy.stats import norm, skew # 正态分布与偏度识别
#读取数据集
train = pd.read_csv('train.csv')
test = pd.read_csv('test.csv')
train.head()   #默认读取前5行

#看看测试数据，测试数据的变量和训练数据是一致的，不过会少了一列Y值
test.head()

一、数据可视化探索

1、缺失值

# id这一列对于模型预测是没有用的，我们先将它在训练集和测试集中去掉
train_ID = train['Id']
test_ID = test['Id']
train.drop("Id", axis = 1, inplace = True)     #axis=1表示列    inplace=1表示不创建新的对象，直接对象进行修改
test.drop("Id", axis = 1, inplace = True)
print("\nThe train data size after dropping Id feature is : {} ".format(train.shape)) 
print("The test data size after dropping Id feature is : {} ".format(test.shape))

#合并train test数据一起处理 
ntrain = train.shape[0]     #pd.shape[0]读取矩阵的行数
ntest = test.shape[0]
y_train = train.SalePrice.values
all_data = pd.concat((train, test)).reset_index(drop=True)   #reset_index表示重置索引
all_data.drop(['SalePrice'], axis=1, inplace=True)
print("all_data size is : {}".format(all_data.shape))             #shape读取矩阵的形状

all_data_na = (all_data.isnull().sum() / len(all_data)) * 100    #按列查看缺失率
all_data_na = all_data_na.drop(all_data_na[all_data_na == 0].index).sort_values(ascending=False)[:30]   #降序查看前30个
missing_data = pd.DataFrame({'Missing Ratio' :all_data_na})
missing_data.head(20)

#将缺失度用图表的方式展示
f, ax = plt.subplots(figsize=(15, 12))
plt.xticks(rotation='90')
sns.barplot(x=all_data_na.index, y=all_data_na)   #条形图
plt.xlabel('Features', fontsize=15)
plt.ylabel('Percent of missing values', fontsize=15)
plt.title('Percent missing data by feature', fontsize=15)

2、相关性探索

热力图是一个快速查看变量之间的相关性，帮助我们理解数据的好办法，根据下图图例，颜色越浅的地方说明两个变量的相关性越强。我们观察最后一列saleprice和其他变量的关系发现GrLivArea（地上居住面积）、OverallQual（整体质量）、GarargeCars（车库能装几辆车）这些变量的颜色比较浅，说明它们对房价的预测能力可能比较强。

#计算所有特征值，每两个之间的相关系数
corrmat = train.corr()   #列之间的相关系数
f, ax = plt.subplots(figsize=(12, 9))
sns.heatmap(corrmat, vmax=.8, square=True);   #square设置热力图矩阵小块形状，默认为False 

#选出10个与房价相关性最强的变量查看相关性系数

k = 10 #选择变量的个数
cols = corrmat.nlargest(k, 'SalePrice')['SalePrice'].index   #pandas.nlargest()输出最大相关排序，排完后选取saleprice列的索引
cm = np.corrcoef(train[cols].values.T)   #设置数据集
sns.set(font_scale=1.25)
#annot_kws设置热力图举证上数字大小字体颜色
hm = sns.heatmap(cm, cbar=True, annot=True, square=True, fmt='.2f', annot_kws={'size': 10}, yticklabels=cols.values, xticklabels=cols.values)
plt.show()

通过上面的热力图，可以有如下的结论：

  1. 'OverallQual', 'GrLivArea' 和 'TotalBsmtSF'与房价的相关性很强，可以后面再深入探索
  2. 'GarageCars'（车库能放多少量车） 和 'GarageArea' （车库面积）和房价同样有比较强的相关性，但注意到这两个变量本身的相关性也很强，因为车库面积和车库能放多少车本身就是有强相关性的，所以这就会涉及到我们之前提到过的多重共线性的问题，对于这两个变量，我们可以去掉一个，只留一个就可以了
  3. 'TotalBsmtSF' 和 '1stFloor' 似乎都是代表地下室面积，这里不确定这两个变量含义的区别，不过它们也有很强的相关性

就发现的和房价相关性较强的变量，再深入看一下它们和房价的关系

#GrLivArea代表含义是居住面积，发现和房价有比较明显的正相关关系，相关系数为0.71
#我们可以在右下角看到两个明显的outlier，价格很低，但面积巨大，我们可以考虑删除掉这两个异常值。
fig, ax = plt.subplots()
ax.scatter(x = train['GrLivArea'], y = train['SalePrice'])  #scatter散点图
plt.ylabel('SalePrice', fontsize=13)
plt.xlabel('GrLivArea', fontsize=13)
plt.show()

#TotalBsmtSF含义为地下室面积，发现地下室面积与房价似乎有更强的潜在线性关系，相关系数为0.61，同时在右侧似乎也有一个异常值存在
var = 'TotalBsmtSF'
data = pd.concat([train['SalePrice'], train[var]], axis=1)
data.plot.scatter(x=var, y='SalePrice', ylim=(0,800000));

#GarageArea含义为车库面积，可以看到车库面积和房价也存在一定的正相关关系，相关系数为0.62
var = 'GarageArea'
data = pd.concat([train['SalePrice'], train[var]], axis=1)
data.plot.scatter(x=var, y='SalePrice', ylim=(0,800000));

对于连续型变量，可以用散点图的方式查看它们和房价之间的关系，下面再来看一下离散型变量，对于离散型变量，可以用箱型图去查看特征不同取值中房价的分布情况

# 首先看一下房屋质量和房价之间的关系，相关系数为0.79，可以看到随着房屋整体质量的变好，房屋的整体价格也在逐渐提高，是一个比较强力的预测变量
var = 'OverallQual'
data = pd.concat([train['SalePrice'], train[var]], axis=1)
f, ax = plt.subplots(figsize=(8, 6))
fig = sns.boxplot(x=var, y="SalePrice", data=data)  #箱线图
fig.axis(ymin=0, ymax=800000);

#再看一下房屋的建造时间和价格的关系，可以看到房屋的建造时间虽然没有和房价有明显的线性关系，但最近建造的新房屋，整体的价格会比较高，这个在后面的特征工程中会用到，注意

ar = 'YearBuilt'
data = pd.concat([train['SalePrice'], train[var]], axis=1)
f, ax = plt.subplots(figsize=(16, 8))
fig = sns.boxplot(x=var, y="SalePrice", data=data)
fig.axis(ymin=0, ymax=800000);
plt.xticks(rotation=90);

#上面是一些比较容易想到的变量，我们可以再看一下壁炉数量和房价的关系，壁炉越多，房价越高
var = 'Fireplaces'
data = pd.concat([train['SalePrice'], train[var]], axis=1)
f, ax = plt.subplots(figsize=(16, 8))
fig = sns.boxplot(x=var, y="SalePrice", data=data)
fig.axis(ymin=0, ymax=800000);
plt.xticks(rotation=90);

二、数据处理

通过可视化探索，我们对数据可以建立起基本的了解，下一步我们对数据的质量做一些处理。首先是异常值，我们已经发现在房屋面积和地下室面积这两个特征里面可能存在异常值，可以把它们删掉

#删除掉异常值
train = train.drop(train[(train['GrLivArea']>4000) & (train['SalePrice']<300000)].index)

##对异常值进行删除重新画图，可以看到异常值没有了，在这里我们只看了一个变量，有兴趣的同学也可以看一下其他的变量异常值的分布情况。
fig, ax = plt.subplots()
ax.scatter(train['GrLivArea'], train['SalePrice'])
plt.ylabel('SalePrice', fontsize=13)
plt.xlabel('GrLivArea', fontsize=13)
plt.show()

#对于地下室面积，做同样处理
train = train.drop(train[(train['TotalBsmtSF']>5000) & (train['SalePrice']<200000)].index)

fig, ax = plt.subplots()
ax.scatter(train['TotalBsmtSF'], train['SalePrice'])
plt.ylabel('SalePrice', fontsize=13)
plt.xlabel('TotalBsmtSF', fontsize=13)
plt.show()

下面我们来看一下缺失值的处理

#再来看一下异常值的分布情况
all_data_na = (all_data.isnull().sum() / len(all_data)) * 100
all_data_na = all_data_na.drop(all_data_na[all_data_na == 0].index).sort_values(ascending=False)[:30]   #降序查看前30个
missing_data = pd.DataFrame({'Missing Ratio' :all_data_na})
missing_data.head(50)

#PoolQC、MiscFeature、Alley的缺失值都在90%以上，可以考虑直接删掉这些特征
all_data= all_data.drop('PoolQC',axis=1)
all_data = all_data.drop('MiscFeature',axis=1)
all_data= all_data.drop('Alley',axis=1)
#删掉这3列后，变量变成了76列
all_data.shape

#Fence栅栏，FireplaceQu壁炉，如果数据缺失的话可能是代表房屋没有栅栏或壁炉， 我们用None填补缺失值，代表没有
all_data["Fence"] = all_data["Fence"].fillna("None")
all_data["FireplaceQu"] = all_data["FireplaceQu"].fillna("None")

#LotFrontage代表房屋前街道的长度，从现实中考虑，房屋前街道的长度可能是会和房屋所在一个街区的房屋相同，所以我们可以考虑同一个街区房屋
#LotFrontage的均值来填补缺失值
#groupby 输出所有数据中不同街区
#transform对主句进行某种同一处理（归一化、标准化等）
all_data["LotFrontage"] = all_data.groupby("Neighborhood")["LotFrontage"].transform(
    lambda x: x.fillna(x.median()))   #输出x均值填补后的x

#Garage相关的车库变量，注意到这些变量的缺失比例是完全相同的，缺失这些变量的房屋可能是没有车库，用None填充
for col in ('GarageType', 'GarageFinish', 'GarageQual', 'GarageCond'):
    all_data[col] = all_data[col].fillna('None')

#同样猜测缺失值缺失的原因可能是因为房屋没有车库，对于连续型变量用0填充
for col in ('GarageYrBlt', 'GarageArea', 'GarageCars'):
    all_data[col] = all_data[col].fillna(0)

#地下室相关连续变量，缺失同样认为房屋可能是没有地下室，用0填充
for col in ('BsmtFinSF1', 'BsmtFinSF2', 'BsmtUnfSF','TotalBsmtSF', 'BsmtFullBath', 'BsmtHalfBath'):
    all_data[col] = all_data[col].fillna(0)

#地下室相关离散变量，同理用None填充
for col in ('BsmtQual', 'BsmtCond', 'BsmtExposure', 'BsmtFinType1', 'BsmtFinType2'):
    all_data[col] = all_data[col].fillna('None')

#Mas为砖石结构相关变量，缺失值我们同样认为是没有砖石结构，用0和none填补缺失值
all_data["MasVnrType"] = all_data["MasVnrType"].fillna("None")
all_data["MasVnrArea"] = all_data["MasVnrArea"].fillna(0)

#MSZoning代表房屋所处的用地类型，给大家先看一下不同取值的个数
#可以看到RL的取值最多
all_data.groupby('MSZoning')['MasVnrType'].count().reset_index()

业务上讲，房屋应该都有所在用地类型，且应该是上面几个值中的一个，所以有缺失值可能不是因为房屋没有用地类型导致，对于这种情况我们就最好不要用None填充，可以考虑用众数

all_data['MSZoning'] = all_data['MSZoning'].fillna(all_data['MSZoning'].mode()[0])   #mode众数

all_data = all_data.drop(['Utilities'], axis=1)
all_data["Functional"] = all_data["Functional"].fillna("Typ")
all_data['Electrical'] = all_data['Electrical'].fillna(all_data['Electrical'].mode()[0])
all_data['KitchenQual'] = all_data['KitchenQual'].fillna(all_data['KitchenQual'].mode()[0])
all_data['Exterior1st'] = all_data['Exterior1st'].fillna(all_data['Exterior1st'].mode()[0])
all_data['Exterior2nd'] = all_data['Exterior2nd'].fillna(all_data['Exterior2nd'].mode()[0])
all_data['SaleType'] = all_data['SaleType'].fillna(all_data['SaleType'].mode()[0])
all_data['MSSubClass'] = all_data['MSSubClass'].fillna("None")

#再来看看还有没有缺失值，发现没有了，到这我们就完成了对缺失值的处理
all_data_na = (all_data.isnull().sum() / len(all_data)) * 100
all_data_na = all_data_na.drop(all_data_na[all_data_na == 0].index).sort_values(ascending=False)
missing_data = pd.DataFrame({'Missing Ratio' :all_data_na})
missing_data.head()

#检测因变量是否符合正态性分布
sns.distplot(train['SalePrice'] , fit=norm);  #distplot直方图

(mu, sigma) = norm.fit(train['SalePrice'])
print( '\n mu = {:.2f} and sigma = {:.2f}\n'.format(mu, sigma))

plt.legend(['Normal dist. ($\mu=$ {:.2f} and $\sigma=$ {:.2f} )'.format(mu, sigma)],  loc='best')#添加图例
plt.ylabel('Frequency')
plt.title('SalePrice distribution')

fig = plt.figure()
res = stats.probplot(train['SalePrice'], plot=plt)   #默认检测是正态分布
plt.show()


从房价的直方图分布和qq-plot可以看出，房价的分布是右偏的，需要对其做一些转换让它符合正太分布

采用log对数变换对房价进行处理，通过下面转换后房价的分布可以看出，房价转换后符合正态分布

train["SalePrice"] = np.log1p(train["SalePrice"])    #数据平滑处理，使其更加服从高斯分布

sns.distplot(train['SalePrice'] , fit=norm);    #distplot直方图

(mu, sigma) = norm.fit(train['SalePrice'])
print( '\n mu = {:.2f} and sigma = {:.2f}\n'.format(mu, sigma))

plt.legend(['Normal dist. ($\mu=$ {:.2f} and $\sigma=$ {:.2f} )'.format(mu, sigma)], loc='best')#添加图例
plt.ylabel('Frequency')
plt.title('SalePrice distribution')

fig = plt.figure()
res = stats.probplot(train['SalePrice'], plot=plt)   #默认检测是正态分布
plt.show()

三、特征工程

对数据做完基本的数据处理后，下面就进入我们最重要最重要的特征工程！

#1、特征值处理  2、特征提取  3、特征选择
#下面我们按照上面的3个步骤，开始进行特征工程，首先我们尝试创造出一些新特征,用下面两个例子给大家做个示范
#加地下室面积、1楼面积、2楼面积加起来可以得到房屋总面积特征
all_data['TotalSF'] = all_data['TotalBsmtSF'] + all_data['1stFlrSF'] + all_data['2ndFlrSF']
#之前可视化时候注意到，建造时间比较近的房子房价比较高，所以新创造一个01特征，如果房屋建造时间在1990年后，则为1，否则是0
all_data['YearBuilt_cut'] = all_data['YearBuilt'].apply(lambda x:1 if x>1990 else 0)

# 我们看一下这两个新加的特征与price的关系
tep=all_data[:ntrain]
tep['SalePrice']=y_train

#房屋总面积特征，发现与房价有比较强的线性关系，但也有2个可能的异常值，可以考虑再去掉这两个异常值
fig, ax = plt.subplots()
ax.scatter(tep['TotalSF'], tep['SalePrice'])
plt.ylabel('SalePrice', fontsize=13)
plt.xlabel('TotalSF', fontsize=13)
plt.show()

#建筑年限，可以看到1990年前建造的房子和1990年后建造的房子房价再分布上有较大的差异
var='YearBuilt_cut'
data=pd.concat([tep['SalePrice'], tep[var]], axis=1)
f,ax=plt.subplots(figsize=(16,8))
fig=sns.boxplot(x=var, y='SalePrice', data=data)
fig.axis(ymin=0, ymax=800000)
plt.xticks(rotation=90)

#对有序离散变量，用label encoder进行编码
from sklearn.preprocessing import LabelEncoder
cols=("FireplaceQu", "BsmtQual", "BsmtCond", "GarageQual", "GarageCond", 'ExterQual', "ExterCond", "HeatingQC", "KitchenQual", 'BsmtFinType1',
     'BsmtFinType2', 'Functional', 'BsmtExposure', 'GarageFinish', 'LandSlope','LotShape', 'PavedDrive', 'Street', 'CentralAir', 'MSSubClass', 'OverallCond', 
        'YrSold', 'MoSold')

# process columns, apply LabelEncoder to categorical features
for c in cols:
    lbl=LabelEncoder()
    lbl.fit(list(all_data[c].values))
    all_data[c]=lbl.transform(list(all_data[c].values))
    
#shape
print('Shape all_data: {}'.format(all_data.shape))

#再来看一下数据中的格式，发现之前的英文字符串已经变成了数字
all_data.head()

#数据集中还有部分非有序性离散变量，我们将他们转换成哑变量的形式（和onehot一个意思）
all_data=pd.get_dummies(all_data)
all_data.shape

all_data.shape

下面进行特征工程的最后一步，特征筛选，为避免多重共线性问题，下面我们会识别皮尔森相关性系数大于0.9的特征，并将这些特征删除。

threshold=0.9
#相关性矩阵
corr_matrix=all_data.corr().abs()
corr_matrix.head()

#只选择矩阵的上半部分
#np.triu中k表示从第几条对角线起保留数据，k=0表示主对角线
#df.where从主体df出发，true返回df本身
upper=corr_matrix.where(np.triu(np.ones(corr_matrix.shape), k=1).astype(np.bool))
upper.head()

#有6列特征需要删除
to_drop=[column for column in upper.columns if any(upper[column]>threshold)]

print('There are %d columns to remove.' % (len(to_drop)))

all_data=all_data.drop(columns=to_drop)
all_data.shape

#之前是把train和test数据集放在一起处理，现在再把他们分开
train=all_data[:ntrain]
test=all_data[ntrain:]

四、建模

在做模型之前，为了衡量模型的泛化误差避免出现过拟合等现象，先要建立一个本地验证集。建模是用训练数据训练模型，并且期望模型在推广到其他数据上时，也有比较好的预测效果，在比赛中对应的就是测试集。但由于测试集没有Y值，即没有真实的房价，所以我们没有办法对模型在测试集上的预测准确率做出判断。因此接下来将训练集八二分，留出20%的数据去验证模型的效果，作为验证集来测试模型的泛化效果。

#导入一些所需要的数据包
from sklearn.linear_model import Ridge, RidgeCV, ElasticNet, LassoCV, LassoLarsCV
from sklearn.model_selection import cross_val_score
#使用cross-val-score函数去作为模型准确性的评估
#np.sqrt()计算各元素的平方根     cv交叉验证折数或者可迭代次数
def rmse_cv(model):
    rmse=np.sqrt(-cross_val_score(model, train, y_train, scoring="neg_mean_squared_error", cv=5))
    return(rmse)
#导入ridge模型
model_ridge=Ridge()

在导入模型后，需要对模型进行调参，对于lasso模型来说，主要需要调节的参数是alpha，alpha决定模型的正则化程度。正则化程度越高，我们的模型越不容易过度拟合。然而，它也会失去灵活性，可能会欠拟合，我们需要在过拟合和欠拟合中寻找一个平衡，找到那个最佳参数。因此，可以列出不同的alpha取值，用上面定义的rmse_cv模型评分参数去看每个参数下模型的分数怎么样

alphas=[0.05, 0.1, 0.3, 1, 3, 5, 10, 15, 30, 50, 75]
cv_ridge=[rmse_cv(Ridge(alpha=alpha)).mean() for alpha in alphas]
cv_ridge=pd.Series(cv_ridge, index=alphas)
cv_ridge.plot(title='Validation - Just Do It')
plt.xlabel('alpha')
plt.ylabel('rmse')

可以看到曲线是一个耐克的标志，先降低后升高，在alpha值为5时，模型的误差最小。

cv_ridge

#alpha参数用我们之前验证过的5，然后用训练集对模型进行训练
clf=Ridge(alpha=5)
clf.fit(train, y_train)
#最后对测试集进行预测，我们就完成了本次项目
predict=clf.predict(test)
sub=pd.DataFrame()
sub['Id']=test_ID
sub['SalePrice']=predict
sub.head(10)
最后编辑于 ：
©著作权归作者所有,转载或内容合作请联系作者
  • 人面猴
    序言：七十年代末，一起剥皮案震惊了整个滨河市，随后出现的几起案子，更是在滨河造成了极大的恐慌，老刑警刘岩，带你破解...
    沈念sama阅读 212,884评论 6赞 492
  • 死咒
    序言：滨河连续发生了三起死亡事件，死亡现场离奇诡异，居然都是意外死亡，警方通过查阅死者的电脑和手机，发现死者居然都...
    沈念sama阅读 90,755评论 3赞 385
  • 救了他两次的神仙让他今天三更去死
    文/潘晓璐 我一进店门，熙熙楼的掌柜王于贵愁眉苦脸地迎上来，“玉大人，你说我怎么就摊上这事。” “怎么了？”我有些...
    开封第一讲书人阅读 158,369评论 0赞 348
  • 道士缉凶录：失踪的卖姜人
    文/不坏的土叔 我叫张陵，是天一观的道长。 经常有香客问我，道长，这世上最难降的妖魔是什么？ 我笑而不...
    开封第一讲书人阅读 56,799评论 1赞 285
  • 港岛之恋（遗憾婚礼）
    正文 为了忘掉前任，我火速办了婚礼，结果婚礼上，老公的妹妹穿的比我还像新娘。我一直安慰自己，他们只是感情好，可当我...
    茶点故事阅读 65,910评论 6赞 386
  • 恶毒庶女顶嫁案：这布局不是一般人想出来的
    文/花漫 我一把揭开白布。 她就那样静静地躺着，像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上，一...
    开封第一讲书人阅读 50,096评论 1赞 291
  • 城市分裂传说
    那天，我揣着相机与录音，去河边找鬼。 笑死，一个胖子当着我的面吹牛，可吹牛的内容都是我干的。 我是一名探鬼主播，决...
    沈念sama阅读 39,159评论 3赞 411
  • 双鸳鸯连环套：你想象不到人心有多黑
    文/苍兰香墨 我猛地睁开眼，长吁一口气：“原来是场噩梦啊……” “哼！你这毒妇竟也来了？” 一声冷哼从身侧响起，我...
    开封第一讲书人阅读 37,917评论 0赞 268
  • 万荣杀人案实录
    序言：老挝万荣一对情侣失踪，失踪者是张志新（化名）和其女友刘颖，没想到半个月后，有当地人在树林里发现了一具尸体，经...
    沈念sama阅读 44,360评论 1赞 303
  • 护林员之死
    正文 独居荒郊野岭守林人离奇死亡，尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 36,673评论 2赞 327
  • 白月光启示录
    正文 我和宋清朗相恋三年，在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 38,814评论 1赞 341
  • 活死人
    序言：一个原本活蹦乱跳的男人离奇死亡，死状恐怖，灵堂内的尸体忽然破棺而出，到底是诈尸还是另有隐情，我是刑警宁泽，带...
    沈念sama阅读 34,509评论 4赞 334
  • 日本核电站爆炸内幕
    正文 年R本政府宣布，位于F岛的核电站，受9级特大地震影响，放射性物质发生泄漏。R本人自食恶果不足惜，却给世界环境...
    茶点故事阅读 40,156评论 3赞 317
  • 男人毒药：我在死后第九天来索命
    文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹，春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 30,882评论 0赞 21
  • 一桩弑父案，背后竟有这般阴谋
    文/苍兰香墨 我抬头看了看天上的太阳。三九已至，却和暖如春，着一层夹袄步出监牢的瞬间，已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 32,123评论 1赞 267
  • 情欲美人皮
    我被黑心中介骗来泰国打工， 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留，地道东北人。 一个月前我还...
    沈念sama阅读 46,641评论 2赞 362
  • 代替公主和亲
    正文 我出身青楼，却偏偏与公主长得像，于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子，可洞房花烛夜当晚...
    茶点故事阅读 43,728评论 2赞 351