2020-08-30--多项式回归01

主要内容

  • 什么是多项式回归
  • scikit-learn中的多项式回归和Pipeline
  • 过拟合和欠拟合
  • 为什么要使用训练数据集和测试数据集
  • 学习曲线
多项式回归简介

考虑下面的数据,虽然我们可以使用线性回归来拟合这些数据,但是这些数据更像是一条二次曲线,相应的方程是y=ax2+bx+c,这是式子虽然可以理解为二次方程,但是我们呢可以从另外一个角度来理解这个式子:

如果将x2理解为一个特征,将x理解为另外一个特征,换句话说,本来我们的样本只有一个特征x,现在我们把他看成有两个特征的一个数据集。多了一个特征x2,那么从这个角度来看,这个式子依旧是一个线性回归的式子。
但是从x的角度来看,他就是一个关于x的二次方程。

以上这样的方式,就是所谓的多项式回归

相当于我们为样本多添加了一些特征,这些特征是原来样本的多项式项,增加了这些特征之后,我们们可以使用线性回归的思路更好的我们的数据。

1.什么是多项式回归

生成一定数量的关于x的二次方程的拟合点

import numpy as np
import matplotlib.pyplot as plt

x = np.random.uniform(-3, 3, size=100)
X = x.reshape(-1, 1)
# 一元二次方程
y = 0.5 * x**2 + x + 2 + np.random.normal(0, 1, 100)
plt.scatter(x, y)
plt.show()

可见这些点的关于二次分布其实不太明显,所以我们先使用线性回归进行测试。

1.1 线性回归测试

from sklearn.linear_model import LinearRegression
lin = LinearRegression()
lin.fit(X,y)
y_predict = lin.predict(X)

# 画图
plt.scatter(x, y)
plt.plot(x,y_predict,color='r')
plt.show()

很明显,我们用一跟直线来拟合一根有弧度的曲线,效果是不好的。

1.2 解决方案,添加一个特征

方案思路:原来所有的数据都在X中,我们要把X^2X看作两个特征,所以现在对X中每一个数据都进行平方, 再将得到的数据集与原数据集进行拼接, 在用新的数据集进行线性回归。

'''修改数据集'''
X2 = np.hstack([X,X**2])
print(X2.shape)
# (100, 2)
'''再次使用线性回归进行测试'''
lin1 = LinearRegression()
lin1.fit(X2,y)
y_predict1 = lin1.predict(X2)

plt.scatter(x, y)
plt.plot(x,y_predict1,color='r')
plt.show()

plt.scatter(x, y)
plt.plot(np.sort(x),y_predict1[np.argsort(x)],color='r')
plt.show()
1
2

注意:对于上图图1的绘制方式,是因为绘画线时,是通过一个点一个点连接画出来的,由于x数据集中的点都是随机的,所以大小不一,画出来的图形就交叉错乱。
所以将x中所有的点排序后,y根据排序后点的索引进行从小到大的描点才可以。

从上图可以看出,当我们添加了一个特征(原来特征的平方)之后,再从x的维度来看,就形成了一条曲线,显然这个曲线对原来数据集的拟合程度是更好的。

# 系数向量,第一个系数是x前面的系数,第二个系数是x平方前面的系数
print(lin1.coef_)
# [1.03557896 0.52749922]
# 截距theta(0)
print(lin1.intercept_)
# 1.8454742506634763

1.3 总结

多项式回归在机器学习算法上并没有新的地方,完全是使用线性回归的思路 他的关键在于为原来的样本,添加新的特征。而我们得到新的特征的方式是原有特征的多项式的组合。 采用这样的方式,我们就可以解决一些非线性的问题。

与此同时需要注意,我们在上一章所讲的PCA是对我们的数据进行降维处理,而我们这一章所讲的多项式回归显然在做一件相反的事情,他让我们的数据升维,在升维之后使得我们的算法可以更好的拟合高纬度的数据。

2.scikit-learn中的多项式回归 和 Pipeline

2.1scikit-learn中的多项式回归

1.数据处理----使用PolynomialFeatures类
import numpy as np

x = np.random.uniform(-3, 3, size=100)
X = x.reshape(-1, 1)
y = 0.5 * x**2 + x + 2 + np.random.normal(0, 1, 100)

# sklearn中对数据进行预处理的函数都封装在preprocessing模块下,包括之前学的归一化StandardScaler
from sklearn.preprocessing import PolynomialFeatures

ploy = PolynomialFeatures(degree=2)
ploy.fit(X)
X2 = ploy.transform(X)
print(X.shape)
# (100, 1)
print(X[:5,:])
# # [[-1.71224922]
# #  [-1.04502905]
# #  [-1.57688397]
# #  [-0.50319217]
# #  [ 2.42563976]]

print(X2.shape)
# (100, 3)
print(X2[:5,:])
# [[ 1.         -1.71224922  2.9317974 ]
#  [ 1.         -1.04502905  1.09208572]
#  [ 1.         -1.57688397  2.48656305]
#  [ 1.         -0.50319217  0.25320236]
#  [ 1.          2.42563976  5.88372824]]
2.新数据线性回归
from sklearn.linear_model import LinearRegression

lin = LinearRegression()
lin.fit(X2,y)
y_predict = lin.predict(X2)

plt.scatter(x,y)
plt.plot(np.sort(x),y_predict[np.argsort(x)],color='r')
plt.show()
print(lin.intercept_)
# 1.9848103836687603
print(lin.coef_)
# [0.         0.95297681 0.5023139 ]

由此可见图形以及训练结果参数与预期结果大概一致。

3. 关于PolynomialFeatures类

为什么使用这个类处理后的数据再进行训练时,能够达到预期的效果呢?
在其内部发生了什么。

import numpy as np
from sklearn.preprocessing import PolynomialFeatures

X = np.arange(1, 11).reshape(-1, 2)

poly = PolynomialFeatures(degree=2)
poly.fit(X)
X2 = poly.transform(X)

print(X.shape)
# (5, 2)
print(X)
# [[ 1  2]
#  [ 3  4]
#  [ 5  6]
#  [ 7  8]
#  [ 9 10]]
print(X2.shape)
# (5, 6)
print(X2)
# [[  1.   1.   2.   1.   2.   4.]
#  [  1.   3.   4.   9.  12.  16.]
#  [  1.   5.   6.  25.  30.  36.]
#  [  1.   7.   8.  49.  56.  64.]
#  [  1.   9.  10.  81.  90. 100.]]

degree=2,将5行2列的矩阵进行多项式转换后变成了5行6列。

分析:

  • 第一列是都为1,对应的是0次幂。
  • 第二列和第三列对应的是原来的x矩阵,此时他有两列一次幂的项。
  • 第四列是原来数据的第一列平方的结果
  • 第六列是原来数据的第二列平方的结果
  • 第五列是原来数据的两列相乘的结果

可以想象如果将degree设置为3,那么将产生一下10个元素

也就是说PolynomialFeatures会穷举出所有的多项式之间的全部组合。

2.2 Pipeline

pipline的英文名字是管道,那么 我们如何使用管道呢,先考虑我们多项式回归的过程:

  1. 使用PolynomialFeatures生成多项式特征的数据集
  2. 如果生成数据幂特别的大,那么特征直接的差距就会很大,导致我们的搜索非常慢,这时候可以进行数据归一化。
  3. 进行线性回归 pipline 的作用就是把上面的三个步骤合并,使得我们不用一直重复这三步
x = np.random.uniform(-3, 3, size=100)
X = x.reshape(-1, 1)
y = 0.5 * x**2 + x + 2 + np.random.normal(0, 1, 100)

from sklearn.pipeline import Pipeline
from sklearn.preprocessing import StandardScaler

# 传入每一步的对象名和类的实例化
poly_reg = Pipeline([
    ("poly", PolynomialFeatures(degree=2)),
    ("std_scaler", StandardScaler()),
    ("lin_reg", LinearRegression())
])

poly_reg.fit(X, y)
y_predict = poly_reg.predict(X)
plt.scatter(x, y)
plt.plot(np.sort(x), y_predict[np.argsort(x)], color='r')
plt.show()

3. 过拟合和欠拟合

import numpy as np
import matplotlib.pyplot as plt

np.random.seed(666)
x = np.random.uniform(-3.0, 3.0, size=100)
X = x.reshape(-1, 1)
y = 0.5 * x**2 + x + 2 + np.random.normal(0, 1, size=100)

plt.scatter(x, y)
plt.show()

3.1 使用线性回归

使用普通线性回归进行预测:

'''1.使用线性回归预测'''
from sklearn.linear_model import LinearRegression
lin = LinearRegression()
lin.fit(X,y)
y_perdict = lin.predict(X)
score = lin.score(X,y)
print(score)     # 0.4953707811865009

plt.scatter(x,y)
plt.plot(np.sort(x),y_perdict[np.argsort(x)],color='r')
plt.show()

查看误差:

# 误差
from sklearn.metrics import mean_squared_error
re = mean_squared_error(y,y_perdict)
print(re)
# 3.0750025765636577

3.2 使用多项式回归

from sklearn.pipeline import Pipeline
from sklearn.preprocessing import PolynomialFeatures
from sklearn.preprocessing import StandardScaler

def PolynomialRegression(degree):
    return Pipeline([
        ("ploy",PolynomialFeatures(degree=degree)),    # 数据处理
        ("std_scaler",StandardScaler()),          # 归一化
        ("lin",LinearRegression())           # 线性回归
    ])

ploy = PolynomialRegression(degree=2)
ploy.fit(X,y)
y_perdict = ploy.predict(X)
re = mean_squared_error(y,y_perdict)
print(re)
# 1.0987392142417858

plt.scatter(x, y)
plt.plot(np.sort(x), y_perdict[np.argsort(x)], color='r')
plt.show()

3.3 过拟合

我们处理具体的数据集时,是不知道degree的值为多少合适的,所以就有以下:

'''degree=10'''
ploy = PolynomialRegression(degree=10)
ploy.fit(X,y)
y_perdict = ploy.predict(X)
re = mean_squared_error(y,y_perdict)
print(re)
# 1.0987392142417858

plt.scatter(x, y)
plt.plot(np.sort(x), y_perdict[np.argsort(x)], color='r')
plt.show()

'''degree=100'''
ploy = PolynomialRegression(degree=100)
ploy.fit(X,y)
y_perdict = ploy.predict(X)
re = mean_squared_error(y,y_perdict)
print(re)
# 1.0987392142417858

plt.scatter(x, y)
plt.plot(np.sort(x), y_perdict[np.argsort(x)], color='r')
plt.show()

degree=10:

degree=100:

可见degree的值越大,表示结果越拟合真实数据集。

但是这条曲线只是原来随机生成的点(分布不均匀),对应的y的预测值连接起来的曲线,不过有x轴很多地方可能没有数据点,所以连接的结果和原来的曲线不一样(不是真实的y曲线)。 下面尝试真正还原原来的曲线(构造均匀分布的原数据集)

生成均匀数据点

np.linspace():生成等差数列的数据。

x = np.linspace(-3, 3, 100)
X_plot = x.reshape(100, 1)

y_plot = 0.5 * x**2 + x + 2 + np.random.normal(0, 1, size=100)

plt.scatter(x,y_plot)
plt.show()

ploy = PolynomialRegression(degree=100)
ploy.fit(X_plot,y_plot)
y_perdict = ploy.predict(X_plot)

plt.scatter(x,y_plot)
plt.plot(X_plot[:,0], y_perdict, color='r')
plt.show()

3.4 总结

总有一条曲线,他能拟合所有的样本点,使得均方误差的值为0

degree210100的过程中,虽然均方误差是越来越小的,从均方误差的角度来看是更加小的 但是他真的能更好的预测我们数据的走势吗?
例如我们选择2.5到3的一个x,使用上图预测出来的y的大小(0或者-1之间)显然不符合我们的数据。

换句话说,我们使用了一个非常高维的数据,虽然使得我们的样本点获得了更小的误差,但是这根曲线完全不是我们想要的样子 他为了拟合我们所有的样本点,变的太过复杂了,这种情况就是过拟合【over-fitting】

相反,在最开始,我们直接使用一根直线来拟合我们的数据,也没有很好的拟合我们的样本特征,当然他犯的错误不是太过复杂了,而是太过简单了 这种情况,我们成为欠拟合-【under-fitting】

对于现在的数据(基于二次方程构造),我们使用低于2项的拟合结果,就是欠拟合;高于2项的拟合结果,就是过拟合

4. 为什么要使用训练数据集和测试数据集

模型泛化程度

使用上节的过拟合结果,我们可以得知,虽然我们训练出的曲线将原来的样本点拟合的非常好,总体的误差非常的小, 但是一旦来了新的样本点,他就不能很好的预测了,在这种情况下,我们就称我们得到的这条弯弯曲曲的曲线,他的泛化能力(由此及彼的能力)非常弱。

训练数据集和测试数据集的意义

我们训练的模型目的是为了使得预测的数据能够尽肯能的准确,在这种情况下,我们观察训练数据集的拟合程度是没有意义的 我们真正需要的是,我们得到的模型的泛化能力更高,解决这个问题的方法也就是使用训练数据集,测试数据集的分离。

测试数据对于我们的模型是全新的数据,如果使用训练数据获得的模型面对测试数据也能获得很好的结果,那么我们就说我们的模型泛化能力是很强的。 如果我们的模型面对测试数据结果很差的话,那么他的泛化能力就很弱。事实上,这是训练数据集更大的意义

train test split
from sklearn.model_selection import train_test_split
X_train, X_test, y_train, y_test = train_test_split(X, y, random_state=666)

直接使用线性回归:

lin_reg = LinearRegression()
lin_reg.fit(X_train, y_train)
y_predict = lin_reg.predict(X_test)
re = mean_squared_error(y_test, y_predict)
print(re)
# 2.2199965269396573

使用多项式回归:degree=2

poly2_reg = PolynomialRegression(degree=2)
poly2_reg.fit(X_train, y_train)
y2_predict = poly2_reg.predict(X_test)
re = mean_squared_error(y_test, y2_predict)
print(re)
# 0.80356410562978997

使用多项式回归:degree=10

poly10_reg = PolynomialRegression(degree=10)
poly10_reg.fit(X_train, y_train)
y10_predict = poly10_reg.predict(X_test)
re = mean_squared_error(y_test, y10_predict)
print(re)
# 0.92129307221507939

使用多项式回归:degree=100

poly100_reg = PolynomialRegression(degree=100)
poly100_reg.fit(X_train, y_train)
y100_predict = poly100_reg.predict(X_test)
re = mean_squared_error(y_test, y100_predict)
print(re)
# 14075796419.234262

上边我们进行的实验实际上在实验模型的复杂度

  • 对于多项式模型来说,我们回归的阶数越高,我们的模型会越复杂,在这种情况下对于我们的机器学习算法来说,通常是有下面一张图的。横轴是模型复杂度(对于不同的算法来说,代表的是不同的意思,比如对于多项式回归来说,是阶数越高,越复杂;对于KNN来说,是K越小,模型越复杂,k越大,模型最简单,当k=n的时候,模型就简化成了看整个样本里,哪种样本最多,当k=1来说,对于每一个点,都要找到离他最近的那个点),另一个维度是模型准确率(也就是他能够多好的预测我们的曲线)

通常对于这样一个图,会有两根曲线:

  • 对于训练数据集来说的,模型越复杂,模型准确率越高,因为模型越复杂,对训练数据集的拟合就越好,相应的模型准确率就越高。

  • 对于测试数据集来说,在模型很简单的时候,模型的准确率也比较低,随着模型逐渐变复杂,对测试数据集的准确率在逐渐的提升,提升到一定程度后,如果模型继续变复杂,那么我们的模型准确率将会进行下降(欠拟合->正合适->过拟合)。

欠拟合和过拟合的标准定义:
欠拟合:算法所训练的模型不能完整表述数据关系
过拟合:算法所训练的模型过多的表达了数据间的噪音关系

5.模型学习复杂度曲线

5.1 学习曲线

生成数据以及分割数据:

import numpy as np
import matplotlib.pyplot as plt

'''1. 生成数据'''
np.random.seed(666)
x = np.random.uniform(-3.0, 3.0, size=100)
X = x.reshape(-1, 1)
y = 0.5 * x**2 + x + 2 + np.random.normal(0, 1, size=100)

plt.scatter(x, y)
plt.show()

'''2. 数据分割'''
from sklearn.model_selection import train_test_split

X_train, X_test, y_train, y_test = train_test_split(X, y, random_state=10)
print(X_train.shape)
# (75, 1)

观察线性回归的学习曲线:观察线性回归模型随着训练数据集增加,性能的变化

from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_squared_error

train_score = []
test_score = []

# 训练次数为训练集中数据个数次,每次训练数据的个数逐次递增
for i in range(1,X_train.shape[0]+1):
    lin = LinearRegression()
    # 每次训练的训练集个数都是i个
    lin.fit(X_train[:i],y_train[:i])

    # 用训练集的前i个数据作为测试集
    y_train_predict = lin.predict(X_train[:i])
    train_score.append(mean_squared_error(y_train[:i], y_train_predict))

    # 使用测试集进行预测
    y_test_predict = lin.predict(X_test)
    test_score.append(mean_squared_error(y_test, y_test_predict))

plt.plot([i for i in range(1,X_train.shape[0]+1)],np.sqrt(train_score),label='train')
plt.plot([i for i in range(1,X_train.shape[0]+1)],np.sqrt(test_score),label='test')
plt.show()

分析:
从趋势上看:

  • 在训练数据集上,误差是逐渐升高的。这是因为我们的训练数据越来越多,我们的数据点越难得到全部的累积,不过整体而言,在刚开始的时候误差变化的比较快,后来就几乎不变了
  • 在测试数据集上,在使用非常少的样本进行训练的时候,刚开始我们的测试误差非常的大,当训练样本大到一定程度以后,我们的测试误差就会逐渐减小,减小到一定程度后,也不会小太多,达到一种相对稳定的情况。
  • 在最终,测试误差和训练误差趋于相等,不过测试误差还是高于训练误差一些,这是因为,训练数据在数据非常多的情况下,可以将数据拟合的比较好,误差小一些,但是泛化到测试数据集的时候,还是有可能多一些误差
1. 封装学习曲线为函数

参数:模型对象,训练集数据,测试集数据

def plot_learning_curve(algo, X_train, X_test, y_train, y_test):
    train_score = []
    test_score = []
    for i in range(1, len(X_train)+1):
        algo.fit(X_train[:i], y_train[:i])

        y_train_predict = algo.predict(X_train[:i])
        train_score.append(mean_squared_error(y_train[:i], y_train_predict))

        y_test_predict = algo.predict(X_test)
        test_score.append(mean_squared_error(y_test, y_test_predict))

    plt.plot([i for i in range(1, len(X_train)+1)],
                               np.sqrt(train_score), label="train")
    plt.plot([i for i in range(1, len(X_train)+1)],
                               np.sqrt(test_score), label="test")
    plt.legend()
    plt.axis([0, len(X_train)+1, 0, 4])
    plt.show()
2. 使用线性回归模型算法进行测试
# 使用线性回归模型
plot_learning_curve(LinearRegression(), X_train, X_test, y_train, y_test)
3. 使用多项式模型测试
from sklearn.preprocessing import PolynomialFeatures
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import Pipeline

def PolynomialRegression(degree):
    return Pipeline([
        ("poly", PolynomialFeatures(degree=degree)),
        ("std_scaler", StandardScaler()),
        ("lin_reg", LinearRegression())
    ])

poly2_reg = PolynomialRegression(degree=2)
plot_learning_curve(poly2_reg, X_train, X_test, y_train, y_test)

首先整体从趋势上,和线性回归的学习曲线是类似的 仔细观察,和线性回归曲线的不同在于,线性回归的学习曲线1.5-1.8左右;二阶多项式回归稳定在了0.9-1.1左右,所以二阶多项式稳定的误差比较低,说明使用二阶线性回归的性能是比较好的。

增大参数degree的值

poly20_reg = PolynomialRegression(degree=20)
plot_learning_curve(poly20_reg, X_train, X_test, y_train, y_test)

在使用20阶多项式回归训练模型的时候可以发现,在数据量偏多的时候,我们的训练数据集拟合的是比较好的,但是测试数据集的误差相对来说增大了很多,离训练数据集比较远,通常这就是过拟合的结果,他的泛化能力是不够的。

5.2 总结

对于欠拟合,它的误差值比最佳的情况趋于稳定的那个位置要高一些,说明无论对于训练数据集还是测试数据集来说,误差都比较大。这是因为我们本身模型选的就不对,所以即使在训练数据集上,他的误差也是大的,所以才会呈现出这样的一种形态。

对于过拟合的情况,在训练数据集上,他的误差不大,和最佳的情况是差不多的,甚至在极端情况,如果degree取更高的话,那么训练数据集的误差会更低。
但是问题在于,测试数据集的误差相对是比较大的,并且训练数据集的误差和测试数据集的误差相差比较大(表现在图上相差比较远),这就说明了此时我们的模型的泛化能力不够好,他的泛化能力是不够的。

最后编辑于
©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 160,165评论 4 364
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 67,720评论 1 298
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 109,849评论 0 244
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 44,245评论 0 213
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 52,596评论 3 288
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 40,747评论 1 222
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 31,977评论 2 315
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 30,708评论 0 204
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 34,448评论 1 246
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 30,657评论 2 249
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 32,141评论 1 261
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 28,493评论 3 258
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 33,153评论 3 238
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 26,108评论 0 8
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 26,890评论 0 198
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 35,799评论 2 277
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 35,685评论 2 272