银行贷款违约风险预测

项目简介：比赛由 Kaggle 举办，要求选手依据客户的信用卡信息，分期付款信息，信用局信息等预测客户贷款是否会违约。一共有8个数据集，包括1个主训练集，1个测试集和6个辅助信息表，主训练集特征主要有用户的个人属性，包括用户的性别，职业，是否有车，是否有房，房子面积等基本信息，辅助信息表包括用户的历史申请信息，历史账户余额信息，分期付款信息，信用卡信息，信用局信息和在信用局上的额度信息。
主训练集探索：
目的：
1、了解数据的缺失值情况、异常值情况，以便做对应的数据清洗
2、了解一下违约贷款和正常贷款用户画像的区别，加深对业务的理解，为我们后面的数据分析展开打基础
1、缺失值探索

# 读取训练集和测试集数据
app_train=pd.read_csv('./application_train.csv')
app_test = pd.read_csv('./application_test.csv')

#定义缺失值检测函数
def missing_values_table(df):
        mis_val = df.isnull().sum() #对缺失值进行统计
        mis_val_percent = 100 * df.isnull().sum() / len(df) #缺失值占比
        mis_val_table = pd.concat([mis_val, mis_val_percent], axis=1) #
        mis_val_table_ren_columns = mis_val_table.rename(
                                    columns = {0 : 'Missing Values', 1 : '% of Total Values'})# rename
        mis_val_table_ren_columns =mis_val_table_ren_columns[mis_val_table_ren_columns.iloc[:,1] != 0].sort_values('% of Total Values', ascending=False).round(1) #对缺失值占比进行排序
       
        return mis_val_table_ren_columns    #返回缺失值列表

missing_values_table(app_train)


数据一共有122列，有67列存在缺失情况，最高缺失值的列缺失度为69.9%，可以发现前面的几列特征缺失度的都是一样的，并且它们都是属于房屋信息，根据这个规律我们可以猜测用户缺失房屋信息可能是因为某种特定原因导致的，而不是随机缺失，这点我们会在后面的特征工程用上。
2、异常值探索
查看用户年龄的数据分布情况（因为数据中，年龄的数值是负数，反映的是申请贷款前，这个用户活了多少天，所以这里我除了负365做了下处理）,发现数据的分布还是比较正常的，最大年龄69岁，最小年龄20岁，没有很异常的数字

查看用户的工作时间分布情况发现（同样工作时间也是负数，所以我除了负365），最小值是-1000年，这里的-1000年明显是一个异常数据，没有人的工作时间是负数的，这可能是个异常值


看一下用户受工作时间的数据分布情况，发现所有的异常值都是一个值，365243，对于这个异常值我的理解是它可能是代表缺失值，所以我的选择是将这个异常值用空值去替换，这样可以保留这个信息，又抹去了异常值，替换之后我们再看一下工作时间的分布情况，正常了很多


3、违约用户画像探索
这部分分析的目标主要是查看违约用户和非违约用户的特征分布情况，目标是对违约用户的画像建立一个基本的了解，为后续特征工程打下基础。比如数据集里面有很多字段，包括性别、年龄、工作时间等等，那么是男性更容易违约还是女性？是年龄大的人更容易违约还是年龄小的人？查看这些数据，可以帮助我们对数据有更好的理解。

# 绘图函数，通过图形可以直观地看到数据的分布情况
def plot_stats(feature,label_rotation=False,horizontal_layout=True):
    temp = app_train[feature].value_counts()
    df1 = pd.DataFrame({feature: temp.index,'Number of contracts': temp.values})

    # Calculate the percentage of target=1 per category value
    cat_perc = app_train[[feature, 'TARGET']].groupby([feature],as_index=False).mean()
    cat_perc.sort_values(by='TARGET', ascending=False, inplace=True)
    
    if(horizontal_layout):
        fig, (ax1, ax2) = plt.subplots(ncols=2, figsize=(12,6))
    else:
        fig, (ax1, ax2) = plt.subplots(nrows=2, figsize=(12,14))
        
    sns.set_color_codes("pastel")
    s = sns.barplot(ax=ax1, x = feature, y="Number of contracts",data=df1)
    if(label_rotation):
        s.set_xticklabels(s.get_xticklabels(),rotation=90)
    
    s = sns.barplot(ax=ax2, x = feature, y='TARGET', order=cat_perc[feature], data=cat_perc)
    if(label_rotation):
        s.set_xticklabels(s.get_xticklabels(),rotation=90)
    plt.ylabel('Percent of target with value 1 [%]', fontsize=10)
    plt.tick_params(axis='both', which='major', labelsize=10)

    plt.show()

  • 首先来看一下男性和女性用户的违约率情况，发现男性用户违约率更高，男性用户违约率为11%，女性为7%，男性用户违约率稍高


  • 再来看一下违约用户和正常用户的年龄分布情况，因为年龄是连续型变量，和性别不同，所以我们使用分布图去看年龄的分布情况，通过数据分布我们可以看到，违约用户年轻用户分布更多，所以我们可以推断的结论是用户年龄越小，违约的可能性越大


  • 可以对用户的年龄进行分捅，进一步观察不同年龄段用户的违约概率，发现确实是用户年龄越小，违约的可能性越高


  • 再来看一下不同贷款类型的违约率情况，对于现金贷款和流动资金循坏贷款，现金贷款的违约率更高


  • 看下用户有没有房和车对违约率的影响，发现没有车和房的人违约率更高，但相差并不是很大



    image.png

  • 从家庭情况看，申请的用户大多已经结婚，单身和世俗结婚的违约率较高，寡居的违约率最低（civil marriage世俗结婚是西方婚姻的一种，和宗教婚姻相反）


  • 看一下子女信息，大部分申请者没有孩子或孩子在3个以下，孩子越多的家庭违约率越高，发现对于有9、11个孩子的家庭违约率达到了100%（和样本少也有关系）


  • 根据申请者的收入类型区分，可以发现休产假和没有工作的人违约率较高，在35%以上，对于这两类人群放款需较为谨慎


  • 从职业来看，越相对收入较低、不稳定的职业违约率越高，比如低廉劳动力、司机、理发师，而像会计、高科技员工等具有稳定高收入的职业违约率就较低


  • 贷款申请人受教育程度大多为中学，学历越低越容易违约


3、特征工程
接下来我们通过对特征的一些理解，尝试做出一些新的特征

  1. CREDIT_INCOME_PERCENT: 贷款金额/客户收入，预期是这个比值越大，说明贷款金额大于用户的收入，用户违约的可能性就越大
  2. ANNUITY_INCOME_PERCENT: 贷款的每年还款金额/客户收入，逻辑与上面一致
  3. CREDIT_TERM: 贷款的每年还款金额/贷款金额，贷款的还款周期，猜测还款周期短的贷款，用户的短期压力可能会比较大，违约概率高
  4. DAYS_EMPLOYED_PERCENT: 用户工作时间/用户年龄
  5. INCOME_PER_CHILD：客户收入/孩子数量，客户的收入平均到每个孩子身上，同样的收入，如果这个人的家庭很大，孩子很多，那么他的负担可能比较重，违约的可能性可能更高
  6. HAS_HOUSE_INFORMATION : 大家还记得我们在处理缺失值时提到的规律吗， 这里我们根据客户是否有缺失房屋信息设计一个二分类特征，如果未缺失的话是1，缺失的是0
app_train_domain = app_train.copy()
app_test_domain = app_test.copy()

app_train_domain['CREDIT_INCOME_PERCENT'] = app_train_domain['AMT_CREDIT'] / app_train_domain['AMT_INCOME_TOTAL']
app_train_domain['ANNUITY_INCOME_PERCENT'] = app_train_domain['AMT_ANNUITY'] / app_train_domain['AMT_INCOME_TOTAL']
app_train_domain['CREDIT_TERM'] = app_train_domain['AMT_ANNUITY'] / app_train_domain['AMT_CREDIT']
app_train_domain['DAYS_EMPLOYED_PERCENT'] = app_train_domain['DAYS_EMPLOYED'] / app_train_domain['DAYS_BIRTH']
app_train_domain['INCOME_PER_CHILD'] = app_train_domain['AMT_INCOME_TOTAL'] / app_train_domain['CNT_CHILDREN']
app_train_domain['HAS_HOUSE_INFORMATION'] = app_train_domain['COMMONAREA_MEDI'].apply(lambda x:1 if x>0 else 0)

对设计出来的连续性特征查看它们在违约用户和非违约用户中的分布情况，可以发现除CREDIT_TERM这个特征外，其他的特征区分度似乎都不是很明显


  • 再来看一下通过缺失值设计的这个特征，通过下图我们可以看到，缺失房屋信息的用户违约概率要明显高于未缺失用户，这在我们模型的预测中可以算是一个比较有效的特征了


  • 对测试集做同样的处理
    4、建模预测
    最后，利用现有的主数据集先进行一次建模预测，模型的话选择LGB模型，速度快，效果好
# 创建实例
model = lgb.LGBMClassifier(n_estimators=1000, objective = 'binary', 
                           class_weight = 'balanced', learning_rate = 0.05, 
                           reg_alpha = 0.1, reg_lambda = 0.1, 
                           subsample = 0.8, n_jobs = -1, random_state = 50)
        
# 训练模型
model.fit(train_features, train_labels, eval_metric = 'auc',
          eval_set = [(valid_features, valid_labels), (train_features, train_labels)],
          eval_names = ['valid', 'train'], categorical_feature = cat_indices,
          early_stopping_rounds = 100, verbose = 200)
  • 通过lgb自带的函数查看特征的重要性
def plot_feature_importances(df):
    
    # 按照重要性排序
    df = df.sort_values('importance', ascending = False).reset_index()
    
    # 对重要性特征进行标准化
    df['importance_normalized'] = df['importance'] / df['importance'].sum()

    # 作图展示
    plt.figure(figsize = (10, 6))
    ax = plt.subplot()
    
    # 取重要性排名前15名进行展示
    ax.barh(list(reversed(list(df.index[:15]))), 
            df['importance_normalized'].head(15), 
            align = 'center', edgecolor = 'k')
    
    # 设置yticks和labels
    ax.set_yticks(list(reversed(list(df.index[:15]))))
    ax.set_yticklabels(df['feature'].head(15))

    plt.xlabel('Normalized Importance'); plt.title('Feature Importances')
    plt.show()
    
    return df
fi_sorted = plot_feature_importances(fi)

5、利用其他数据集信息
首先是信用局信息，数据集中的每一行代表的是主训练集中的申请人曾经在其他金融机构申请的贷款信息，可以看到数据集中同样有一列是“SK_ID_CURR'，和主训练集中的列一致，可以通过这一列去把辅助训练集和主训练集做left join，但需要注意的一点是，一个SK_ID_CURR可能会对应多个SK_ID_BUREAU，即一个申请人如果在其他金融机构曾经有多条贷款信息的话，这里就会有多条记录，因为模型训练每个申请人在数据集中只能有一条记录，所以说不能直接把辅助训练集去和主训练集join，一般来说需要去计算一些统计特征（groupby操作）

# 读取信用局数据
bureau = pd.read_csv('./bureau.csv')
bureau.head()

#针对每个贷款申请人计算他们在其他金融机构历史上的贷款数量
previous_loan_counts = bureau.groupby('SK_ID_CURR', as_index=False)['SK_ID_BUREAU'].count().rename(columns = {'SK_ID_BUREAU': 'previous_loan_counts'})
previous_loan_counts.head()
  • 然后我们再把计算出来的统计特征和主训练集做left join，可以看到我们的统计特征就出现在了最后一列
app_train = app_train.merge(previous_loan_counts, on = 'SK_ID_CURR', how = 'left')

# 用0填充空值
app_train['previous_loan_counts'] = app_train['previous_loan_counts'].fillna(0)
app_train.head()
  • 在做出新特征后，我们往往还需要检验新特征是否对预测有区分度，不是所有的新特征都是有用的，有些没有用的特征加到数据集里反而会降低预测值，通过查看违约和非违约用户previous_loan_counts的统计属性发现，虽然非违约用户的平均贷款申请数量要略多于违约用户，但差异很小，所以其实很难判断这个特征对预测是否是有用的，我们可以尝试在做一些更多的特征
  • 定义一个查看分布的函数，以后再做出新特征时，我们可以用这个函数快速查看新的特征在违约用户和非违约用户中的分布情况
def kde_target(var_name, df):
    
    # 计算新变量和目标之间的相关性
    corr = df['TARGET'].corr(df[var_name])
    
    # 计算违约和不违约用户的中位数
    avg_repaid = df.ix[df['TARGET'] == 0, var_name].median()
    avg_not_repaid = df.ix[df['TARGET'] == 1, var_name].median()
    
    plt.figure(figsize = (12, 6))
    
    # 作图展示分布情况
    sns.kdeplot(df.ix[df['TARGET'] == 0, var_name], label = 'TARGET == 0')
    sns.kdeplot(df.ix[df['TARGET'] == 1, var_name], label = 'TARGET == 1')
    
    plt.xlabel(var_name); plt.ylabel('Density'); plt.title('%s Distribution' % var_name)
    plt.legend()
kde_target('previous_loan_counts', app_train)


连续型变量特征提取

  • 对于连续型变量，可以采用计算它们的统计值来作为特征
# 定义函数完成特征提取
def agg_numeric(df, group_var, df_name):
    # 移除id变量
    for col in df:
        if col != group_var and 'SK_ID' in col:
            df = df.drop(columns = col)
            
    group_ids = df[group_var]
    numeric_df = df.select_dtypes('number')
    numeric_df[group_var] = group_ids

    # 通过groupby和agg函数计算出统计特征
    agg = numeric_df.groupby(group_var).agg(['count', 'mean', 'max', 'min', 'sum']).reset_index()
    columns = [group_var]
    
    # 列名重命名
    for var in agg.columns.levels[0]:
        if var != group_var:
            for stat in agg.columns.levels[1][:-1]:
                columns.append('%s_%s_%s' % (df_name, var, stat))

    agg.columns = columns
    return agg

bureau_agg_new = agg_numeric(bureau.drop(columns = ['SK_ID_BUREAU']), group_var = 'SK_ID_CURR', df_name = 'bureau')
# 提取的特征加入到主训练集中
app_train = app_train.merge(bureau_agg_new , on = 'SK_ID_CURR', how = 'left')
bureau_agg_new.head()


离散型变量特征提取
对于离散型变量不能计算它们的统计特征，但可以计算离散型特征中每个取值的个数，通过这种方式来获取到一些信息
比如：原始数据是这样的

最终我们希望达到的效果是：

除此之外，我们还可以计算每个用户的每个值的个数在整体中的比例：

# 同样设置一个函数来提取离散变量的特征
def count_categorical(df, group_var, df_name):
 
    # 首先先把数据集中的离散特征变成哑变量
    categorical = pd.get_dummies(df.select_dtypes('object'))
    categorical[group_var] = df[group_var]

    # 通过groupby和agg函数对特征进行统计
    categorical = categorical.groupby(group_var).agg(['sum', 'mean'])
    
    column_names = []
    
    # 重命名列
    for var in categorical.columns.levels[0]:
        for stat in ['count', 'count_norm']:
            column_names.append('%s_%s_%s' % (df_name, var, stat)）
    categorical.columns = column_names
    
    return categorical

bureau_counts = count_categorical(bureau, group_var = 'SK_ID_CURR', df_name = 'bureau')
#做好的特征可以合并到主训练集中
app_train = app_train.merge(bureau_counts , left_on = 'SK_ID_CURR', right_index = True, how = 'left')
app_train.head()


整合所有数据集

  • 重新读取一遍数据集，把数据集还原到初始状态
app_train=pd.read_csv('./application_train.csv')
app_test = pd.read_csv('./application_test.csv')
bureau = pd.read_csv('./bureau.csv')
previous_application = pd.read_csv('./previous_application.csv')
  • 把之前对主训练集做的特征重新加入到数据集
app_train['CREDIT_INCOME_PERCENT'] = app_train['AMT_CREDIT'] / app_train['AMT_INCOME_TOTAL']
app_train['ANNUITY_INCOME_PERCENT'] = app_train['AMT_ANNUITY'] / app_train['AMT_INCOME_TOTAL']
app_train['CREDIT_TERM'] = app_train['AMT_ANNUITY'] / app_train['AMT_CREDIT']
app_train['DAYS_EMPLOYED_PERCENT'] = app_train['DAYS_EMPLOYED'] / app_train['DAYS_BIRTH']
app_train['INCOME_PER_CHILD'] = app_train['AMT_INCOME_TOTAL'] / app_train['CNT_CHILDREN']
app_train['HAS_HOUSE_INFORMATION'] = app_train['COMMONAREA_MEDI'].apply(lambda x:1 if x>0 else 0)


app_test['CREDIT_INCOME_PERCENT'] = app_test['AMT_CREDIT'] / app_test['AMT_INCOME_TOTAL']
app_test['ANNUITY_INCOME_PERCENT'] = app_test['AMT_ANNUITY'] / app_test['AMT_INCOME_TOTAL']
app_test['CREDIT_TERM'] = app_test['AMT_ANNUITY'] / app_test['AMT_CREDIT']
app_test['DAYS_EMPLOYED_PERCENT'] = app_test['DAYS_EMPLOYED'] / app_test['DAYS_BIRTH']
app_test['INCOME_PER_CHILD'] = app_test['AMT_INCOME_TOTAL'] / app_test['CNT_CHILDREN']
app_test['HAS_HOUSE_INFORMATION'] = app_test['COMMONAREA_MEDI'].apply(lambda x:1 if x>0 else 0)
  • 两个函数完成之前对信用局数据中连续变量和离散变量的特征提取
bureau_counts = count_categorical(bureau, group_var = 'SK_ID_CURR', df_name = 'bureau')
bureau_agg_new = agg_numeric(bureau.drop(columns = ['SK_ID_BUREAU']), group_var = 'SK_ID_CURR', df_name = 'bureau')
bureau_agg_new.head()
  • 给训练集和测试集增加信用局相关特征
app_train = app_train.merge(bureau_counts, on = 'SK_ID_CURR', how = 'left')
app_train = app_train.merge(bureau_agg_new, on = 'SK_ID_CURR', how = 'left')

app_test = app_test.merge(bureau_counts, on = 'SK_ID_CURR', how = 'left')
app_test = app_test.merge(bureau_agg_new, on = 'SK_ID_CURR', how = 'left')
  • 对于历史贷款信息数据，做同样的操作
previous_appication_counts = count_categorical(previous_application, group_var = 'SK_ID_CURR', df_name = 'previous_application')
previous_appication_agg_new = agg_numeric(previous_application, group_var = 'SK_ID_CURR', df_name = 'previous_application')

app_train = app_train.merge(previous_appication_counts, on = 'SK_ID_CURR', how = 'left')
app_train = app_train.merge(previous_appication_agg_new, on = 'SK_ID_CURR', how = 'left')
app_test = app_test.merge(previous_appication_counts, on = 'SK_ID_CURR', how = 'left')
app_test = app_test.merge(previous_appication_agg_new, on = 'SK_ID_CURR', how = 'left')

print(app_train.shape)
print(app_test.shape)

特征筛选
在之前的一系列的特征工程中，我们给训练集和测试集增加了很多新的特征，特征也膨胀到了600多列，在最后建模之前，还需要对这些加入的特征再做一次筛选，排除一些具有共线性的特征以提高模型的效果，我们可以计算变量与变量之间的相关系数，来快速移除一些相关性过高的变量，这里可以定义一个阈值是0.8，即移除每一对相关性大于0.8的变量中的其中一个变量

corrs = app_train.corr()

# 设置阈值
threshold = 0.8
above_threshold_vars = {}
for col in corrs:
    above_threshold_vars[col] = list(corrs.index[corrs[col] > threshold])

# Track columns to remove and columns already examined
cols_to_remove = []
cols_seen = []
cols_to_remove_pair = []

# Iterate through columns and correlated columns
for key, value in above_threshold_vars.items():
    # Keep track of columns already examined
    cols_seen.append(key)
    for x in value:
        if x == key:
            next
        else:
            # Only want to remove one in a pair
            if x not in cols_seen:
                cols_to_remove.append(x)
                cols_to_remove_pair.append(key)
            
cols_to_remove = list(set(cols_to_remove))
print('Number of columns to remove: ', len(cols_to_remove))
# 在结果中可以看到，我们一共要移除189列具有高相关性的变量

#训练集和测试集都移除对应的列,把列数降低到426
train_corrs_removed = app_train.drop(columns = cols_to_remove)
test_corrs_removed = app_test.drop(columns = cols_to_remove)

print('Training Corrs Removed Shape: ', train_corrs_removed.shape)
print('Testing Corrs Removed Shape: ', test_corrs_removed.shape)


建模预测
用之前建立好的模型再进行一次建模预测,到此，就完成了这个项目

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