
简介这份资源面向收入分配、劳动经济学与机器学习方向的研究生、学者及政策研究者围绕流动人口劳动收入风险的测算及其分配效应展开。核心方法是用机器学习复原个体收入分布再以方差、偏度和峰度分别刻画收入波动性、增长空间与极端收入可能进而考察风险补偿的异质性及其对收入差距的扩大或缩小作用。资源包共1个PDF文件约756KB内容为论文复现说明与完整MATLAB代码解析涵盖数据准备与预处理、分位数回归森林式的收入分布复原、三阶矩风险测算、风险补偿分析、收入分配效应评估及结果可视化等环节代码含模拟数据生成与逐段注释便于读者理解建模逻辑并迁移到自己的数据。目前已有56人学习适合希望掌握风险度量思路、复现实证流程并获取可运行脚本的读者参考。1. 劳动收入风险测算当分位数回归森林遇上流动人口收入差距流动人口收入分配研究里有个长期被简化处理的问题我们习惯用均值回归看教育、行业、区域对收入的影响但政策制定者真正关心的是尾部——最低收入群体是否被挤压、收入差距是否在代际间固化。传统分位数回归能刻画条件分布但面对流动人口数据中普遍存在的非线性、交互效应和异质性线性设定往往力不从心。分位数回归森林Quantile Regression Forest, QRF把随机森林的集成学习能力和分位数回归的分布视角结合起来不需要预设函数形式就能估计任意分位点上的条件收入分布。这篇实战笔记围绕劳动收入风险测算这个目标把流动人口收入分配效应的完整复现路径拆开从数据构造、QRF建模、基尼系数分解到多维风险因子的边际效应分析。适合有Python或MATLAB基础、正在做收入分配或劳动经济学量化研究的人也适合想把机器学习方法落到社会科学场景的工程师。2. 分位数回归森林为什么比线性分位数回归更适合流动人口数据2.1 流动人口收入数据的三个非线性特征流动人口收入数据有几个绕不开的麻烦。第一收入对数的条件分布随教育年限增加呈现明显的左偏收窄线性分位数回归假设各分位点系数差异只体现在截距上但实际数据里教育回报率在低分位点和高分位点可能差出两倍以上。第二行业和区域变量之间存在强交互建筑业的区域收入差异远大于信息技术业线性模型要显式构造几十个交互项才能捕捉。第三样本自选择问题严重高技能流动人口向特定城市集聚导致协变量分布本身随分位点变化。QRF的处理逻辑是对每一棵树的每个叶节点保留落入该节点的所有训练样本的Y值预测时把目标样本落到各棵树的叶节点汇总所有叶节点里的Y值经验分布再取目标分位数。这样条件分布的形状完全由数据驱动不需要任何函数形式假设。2.2 QRF的核心参数与调参逻辑用Python的quantile_forest库或MATLAB的TreeBagger改造都能实现。关键参数有三个n_estimators控制森林规模一般500到1000棵足够稳定min_samples_leaf决定叶节点最小样本量流动人口数据建议不低于20否则尾部估计方差过大max_features控制分裂时随机选取的特征数通常取总特征数的三分之一到二分之一。from quantile_forest import RandomForestQuantileRegressor import numpy as np # X_train: 特征矩阵, y_train: 对数收入 # 参数说明 # n_estimators800: 森林规模低于500时基尼系数估计波动超过5% # min_samples_leaf25: 叶节点最小样本保证尾部估计稳定 # max_features0.4: 每次分裂随机选40%特征平衡相关性和多样性 qrf RandomForestQuantileRegressor( n_estimators800, min_samples_leaf25, max_features0.4, random_state42 ) qrf.fit(X_train, y_train) # 预测0.1, 0.5, 0.9分位点的条件收入 quantiles [0.1, 0.5, 0.9] y_pred qrf.predict(X_test, quantilesquantiles)这段代码里min_samples_leaf是最需要反复试的参数。设太小0.1分位点的预测值会在不同随机种子下跳变设太大高收入尾部的异质性被抹平。我一般会从15开始试每次加5观察0.1和0.9分位点预测值的标准差变化选标准差下降曲线拐点对应的值。2.3 从条件分布到收入差距测度拿到QRF预测的各分位点条件收入后计算基尼系数有两种路径。一是直接用预测的分位点值构造反事实分布用梯形近似积分算基尼二是用QRF输出的完整条件分布函数对每个样本积分得到期望收入再算基尼。前者计算快但精度受分位点数量影响后者更准但需要保存每棵树的叶节点样本。def gini_from_quantiles(quantile_preds, quantile_levels): 从分位点预测值近似计算基尼系数 quantile_preds: shape (n_samples, n_quantiles) quantile_levels: 如 [0.1, 0.2, ..., 0.9] # 对每个样本用分位点值近似其收入分布 mean_income np.mean(quantile_preds, axis1) # 按均值排序 sorted_idx np.argsort(mean_income) sorted_income mean_income[sorted_idx] n len(sorted_income) # 基尼系数标准公式 cumulative np.cumsum(sorted_income) gini (2 * np.sum((np.arange(1, n1)) * sorted_income) / (n * np.sum(sorted_income)) - (n 1) / n) return gini这里用分位点均值代替真实收入会低估基尼系数约3到5个百分点因为分位点之间的分布信息被丢弃了。如果要做严谨的论文复现建议用QRF的predict方法配合quantiles参数输出更多分位点比如19个或者直接调用quantile_forest的predict_distribution方法。3. 多维度劳动收入风险因子的构造与边际效应测算3.1 风险因子的四个维度与量化方式劳动收入风险不能只看收入波动流动人口的特殊性在于制度性分割和城市融入成本。我一般从四个维度构造风险因子维度具体变量量化方式数据来源就业稳定性合同类型、工作年限、行业波动率合同虚拟变量行业GDP波动流动人口动态监测技能可替代性职业任务 Routine 指数O*NET任务得分映射O*NET数据库城市融入成本房租收入比、通勤时间、社保参与连续变量标准化城市统计年鉴制度性分割户籍限制、子女教育可及性城市落户门槛指数政策文本量化这四个维度不是拍脑袋来的。就业稳定性和技能可替代性决定收入的下行风险城市融入成本和制度性分割决定收入的风险溢价补偿是否充分。QRF的优势在于可以估计每个风险因子在不同分位点上的边际效应而不是只给一个平均效应。3.2 用QRF做条件分位数边际效应具体做法是对每个风险因子构造反事实样本——把该因子分别设为样本的10分位和90分位值其他变量保持不变用训练好的QRF预测两个反事实样本在0.1、0.5、0.9分位点的收入差值就是该因子从低到高的边际效应。def marginal_effect_qrf(qrf, X, feature_idx, q_levels[0.1, 0.5, 0.9]): 计算单个特征从10分位到90分位的边际效应 X: 原始特征矩阵 feature_idx: 要分析的特征列索引 X_low X.copy() X_high X.copy() # 将该特征分别设为10分位和90分位 low_val np.percentile(X[:, feature_idx], 10) high_val np.percentile(X[:, feature_idx], 90) X_low[:, feature_idx] low_val X_high[:, feature_idx] high_val # 预测各分位点收入 pred_low qrf.predict(X_low, quantilesq_levels) pred_high qrf.predict(X_high, quantilesq_levels) # 边际效应 高分位反事实 - 低分位反事实 effects pred_high - pred_low return effects.mean(axis0), effects.std(axis0)这个函数返回每个分位点上的平均边际效应和标准差。标准差很重要——如果某个风险因子在0.1分位点的边际效应标准差很大说明该因子对低收入群体的影响异质性极强政策上需要更精细的 targeting。3.3 收入差距分解QRF Shapley值要回答“哪个风险因子对收入差距贡献最大”可以用Shapley值分解。QRF的预测函数不是可加的但可以用shap库的TreeExplainer对QRF做近似分解。注意quantile_forest的模型对象需要先转成sklearn兼容格式或者直接用shap.Explainer配合自定义预测函数。import shap # 用0.5分位点预测函数构造explainer def median_predict(X): return qrf.predict(X, quantiles[0.5]).flatten() explainer shap.Explainer(median_predict, X_train) shap_values explainer(X_test) # 按风险因子分组汇总Shapley值 risk_groups { 就业稳定性: [0, 1, 2], 技能可替代性: [3, 4], 城市融入成本: [5, 6, 7], 制度性分割: [8, 9] } for group, idx in risk_groups.items(): group_shap np.abs(shap_values.values[:, idx]).mean() print(f{group} 平均|SHAP|: {group_shap:.4f})这里有个坑shap.Explainer默认用interventional特征扰动对QRF这种基于样本的模型tree_path_dependent模式更稳定但计算慢。如果特征间相关性高比如房租收入比和城市GDPShapley值会在相关特征间分摊贡献解释时要小心。4. 避坑与排查QRF做收入分配研究时最容易翻车的五个地方4.1 分位点交叉问题现象预测的0.3分位收入高于0.4分位收入出现分位点交叉。原因QRF对每个分位点独立优化没有单调性约束。当样本量小或叶节点样本少时不同分位点的预测值可能乱序。解决用quantile_forest的monotonic_cst参数强制单调或者后处理时对预测值做排序平滑。我一般会在预测后加一步np.sort沿分位点维度排序简单有效。4.2 外推能力为零现象测试集里某个城市或行业的样本QRF预测值全部接近训练集均值。原因QRF是纯非参数方法对训练集中未出现的特征组合没有外推能力。流动人口数据里城市和行业的组合很多训练集覆盖不全时尾部预测会塌缩。解决要么扩大训练集覆盖要么对城市和行业做目标编码target encoding把高基数类别变量转成连续变量。目标编码要用交叉验证防止泄漏。4.3 基尼系数被低估现象QRF预测收入算出的基尼系数比原始数据低5到8个百分点。原因QRF预测的是条件分位数不是真实收入。分位数之间的分布信息丢失导致预测分布比真实分布更集中。解决用更多分位点至少19个或者用QRF的完整分布预测功能。如果论文要求基尼系数精度建议用predict_distribution输出每个样本的完整条件分布再积分算基尼。4.4 特征重要性误导现象QRF的feature_importances_显示城市GDP重要性最高但Shapley值显示制度性分割贡献最大。原因QRF的默认特征重要性基于分裂时的不纯度减少对高基数类别变量和连续变量有偏好。Shapley值基于预测贡献更可靠。解决论文里报告Shapley值分解结果QRF自带重要性只作为参考。如果一定要用用permutation_importance替代。4.5 样本权重处理现象流动人口数据里某些省份样本量特别大QRF预测被这些省份主导。原因QRF的叶节点样本汇总时没有考虑抽样权重样本量大的组自然主导预测。解决quantile_forest支持sample_weight参数在fit时传入权重。权重可以按省份样本量倒数或按流动人口监测的设计权重设置。5. 用MATLAB复现QRF收入分配分析的替代路径5.1 MATLAB的TreeBagger改造方案如果团队用MATLAB为主可以用TreeBagger加自定义分位数预测函数实现QRF。核心思路是训练TreeBagger回归森林然后用predict方法获取每棵树的叶节点样本索引再手动汇总分位数。% 训练回归森林 % X: 特征矩阵, Y: 对数收入 % MinLeafSize25: 叶节点最小样本对应Python的min_samples_leaf % NumTrees800: 森林规模 rf TreeBagger(800, X, Y, Method, regression, ... MinLeafSize, 25, NumPredictorsToSample, round(size(X,2)*0.4)); % 获取每棵树的叶节点样本 [Y_pred, node_indices] predict(rf, X_test); % 手动计算分位数 % 对每个测试样本收集所有树的叶节点训练样本Y值 n_trees rf.NumTrees; n_test size(X_test, 1); quantile_levels [0.1, 0.5, 0.9]; qrf_preds zeros(n_test, length(quantile_levels)); for i 1:n_test leaf_samples []; for t 1:n_trees % 获取第t棵树中测试样本i落入的叶节点 node node_indices(i, t); % 找到训练集中落入同一叶节点的样本 train_node predict(rf.Trees{t}, X, Trees, t); % 这里需要根据TreeBagger的内部结构提取叶节点样本 % 实际实现时建议用fitrtree逐个训练并保存叶节点样本索引 end qrf_preds(i, :) quantile(leaf_samples, quantile_levels); endMATLAB的TreeBagger不直接暴露叶节点样本索引上面的代码需要改用fitrtree逐棵树训练并手动管理叶节点样本。这条路可行但代码量大适合对MATLAB生态有强依赖的团队。5.2 MATLAB与Python的混合工作流更务实的做法是数据清洗和描述统计在MATLAB里做QRF建模和Shapley分解用Python结果导回MATLAB做可视化。两边用CSV或HDF5交换数据。这样既利用了MATLAB在面板数据处理上的便利又用上了Python的quantile_forest和shap生态。% MATLAB端导出清洗后的数据 writetable(data_clean, income_data.csv); % Python端处理完后导回 % 在MATLAB中读取预测结果 qrf_results readtable(qrf_predictions.csv); % 用MATLAB做分位数回归系数可视化 figure; boxplot(qrf_results.marginal_effect, qrf_results.risk_factor); ylabel(边际效应); title(各风险因子在不同分位点的边际效应分布);这个混合流程的坑在于编码和缺失值处理要统一。MATLAB的readtable默认把空字符串读成NaNPython的pandas读成NaN但类型可能不同。建议导出时统一用-999标记缺失两边都做显式转换。5.3 验证QRF结果稳定性的三个检查不管用Python还是MATLAB跑完QRF后我会做三个检查。第一换三个随机种子重跑看0.1和0.9分位点预测值的相关系数是否都在0.95以上。第二把训练集随机分半分别训练QRF看两半数据在测试集上的基尼系数差异是否小于0.02。第三对关键风险因子用线性分位数回归跑一遍看QRF的边际效应是否在线性模型置信区间内——如果在说明非线性效应不强用线性模型也能交差如果不在QRF的非线性捕捉就是论文的核心贡献点。这三个检查花不了半小时但能避免审稿人问“你的QRF结果稳定吗”时拿不出证据。我吃过这个亏后来每次跑QRF都先把这三个检查跑完再往下做。希望帮到你。本文还有配套的精品资源点击获取