简介这份资源面向劳动经济学、收入分配与机器学习交叉领域的研究生、学者及政策研究者围绕流动人口劳动收入风险的测算及其收入分配效应展开实证复现。核心方法是用机器学习复原个体收入分布再以方差、偏度和峰度分别刻画收入整体波动性、增长空间与极端收入可能性进而考察风险补偿的异质性及其对收入差距的扩大或缩小作用。资源包共1个PDF文件约756KB内容涵盖数据准备、收入分布复原、风险测算、风险补偿分析、收入分配效应分析与结果可视化等完整环节并附有MATLAB代码及逐段解释便于读者理解分位数回归森林等方法的实现思路与替代方案。已有56人学习适合希望掌握风险度量建模流程、复现实证结论或将其迁移至自身研究场景的读者参考。1. 劳动收入风险测算到底在算什么从流动人口的收入分布说起流动人口的收入问题真正棘手的从来不是均值高低而是分布形态和尾部风险。一个月薪中位数 6000 元的群体可能同时存在大量月入 3000 以下的底层劳动者和少数月入 3 万以上的高技能人才这种右偏厚尾的分布用普通最小二乘回归去拟合得到的只是条件均值完全看不到收入差距的结构性来源。劳动收入风险测算要解决的就是量化这种分布内部的离散程度和不确定性——哪些因素在拉大差距哪些因素在收窄差距不同分位点上各因素的影响方向是否一致。这篇论文复现的核心思路是用分位数回归森林Quantile Regression Forest, QRF替代传统分位数回归在流动人口微观数据上估计不同分位点10%、25%、50%、75%、90%的收入条件分布再基于反事实分布计算基尼系数和收入差距指标最后做多维度风险因素教育、行业、户籍、区域、职业稳定性等的效应分解。适合做劳动经济学实证、收入分配研究、以及想用机器学习方法替代传统计量工具的研究生和从业者。MATLAB 在这个流程里承担数据处理、QRF 训练、分位数预测和基尼系数计算的全链路任务。2. 分位数回归森林为什么比线性分位数回归更适合收入数据原理与选型2.1 收入分布的厚尾和非线性让线性分位数回归力不从心传统分位数回归Koenker Bassett, 1978假设条件分位数是自变量的线性函数通过最小化非对称损失函数求解系数。这个假设在收入数据上经常翻车教育回报率在不同收入水平上差异显著高收入群体教育边际回报更高行业效应在低收入端和高收入端方向可能相反年龄-收入曲线在不同分位点上形状完全不同。线性设定强行把这些异质性压成一条直线估计出来的系数虽然显著但经济含义已经失真。分位数回归森林Meinshausen, 2006的思路完全不同。它基于随机森林框架但不输出条件均值而是输出完整的条件分布函数。具体做法是每棵树的每个叶节点记录落入该节点的所有训练样本的因变量值预测时对新样本找到其落入的叶节点汇总所有树的叶节点样本权重得到经验条件分布再从中提取任意分位点。这个方法不假设任何函数形式自动捕捉非线性和交互效应对厚尾分布尤其友好。选型理由很直接收入数据的分位数处理效应高度异质QRF 能给出每个样本的完整条件分布而非单点估计这是后续做反事实分析和基尼系数分解的前提。MATLAB 没有内置 QRF 函数需要自己实现或调用第三方工具箱但核心逻辑并不复杂。2.2 用 MATLAB 实现 QRF 的核心步骤与参数设置下面是一个可运行的 QRF 实现框架基于 MATLAB 的TreeBagger做改造。关键点在于不直接用TreeBagger的预测输出而是提取每棵树的叶节点样本索引自行构建条件分布。function [qr_models, leaf_samples] trainQRF(X, Y, nTrees, minLeaf, mtry) % trainQRF 训练分位数回归森林 % 输入: % X - 特征矩阵 (n x p) % Y - 响应变量 (n x 1)此处为对数收入 % nTrees - 树的数量建议 500-1000 % minLeaf - 叶节点最小样本数建议 5-20 % mtry - 每次分裂随机选取的特征数建议 p/3 % 输出: % qr_models - 训练好的随机森林模型 % leaf_samples - 每棵树的叶节点样本索引 cell 数组 n size(X, 1); leaf_samples cell(nTrees, 1); % 用 TreeBagger 训练但关闭自带预测 qr_models TreeBagger(nTrees, X, Y, ... Method, regression, ... MinLeafSize, minLeaf, ... NumPredictorsToSample, mtry, ... OOBPrediction, on, ... InBagFraction, 0.632); % 提取每棵树的叶节点信息 for t 1:nTrees tree qr_models.Trees{t}; % 获取叶节点编号 leaf_nodes find(tree.IsBranchNode 0); % 获取每个训练样本落入的叶节点 [~, node_idx] predict(tree, X); leaf_samples{t} node_idx; end end逻辑说明TreeBagger的MinLeafSize控制叶节点最小样本数这个参数直接决定条件分布的平滑程度——太小则过拟合太大则分布估计粗糙。NumPredictorsToSample控制随机性收入数据特征间相关性高时建议取p/3而非默认的p。InBagFraction设为 0.632 是标准自助采样比例。参数说明nTrees取 500 起步1000 通常足够稳定minLeaf在收入数据上建议 10-20因为微观调查样本量通常几千到几万叶节点太小会导致分位数估计方差过大mtry如果特征数在 15-30 之间取 5-10 比较合理。训练完成后预测新样本的条件分布function cond_dist predictQRF(qr_models, leaf_samples, X_new, Y_train) % predictQRF 预测新样本的条件分布 % 输出 cond_dist 为 n_new x n_train 的权重矩阵 n_new size(X_new, 1); n_train length(Y_train); nTrees length(leaf_samples); cond_dist zeros(n_new, n_train); for t 1:nTrees tree qr_models.Trees{t}; % 新样本落入的叶节点 [~, leaf_idx_new] predict(tree, X_new); % 训练样本落入的叶节点 leaf_idx_train leaf_samples{t}; for i 1:n_new % 找到同叶节点的训练样本 same_leaf (leaf_idx_train leaf_idx_new(i)); if sum(same_leaf) 0 cond_dist(i, same_leaf) cond_dist(i, same_leaf) 1; end end end % 归一化为权重 cond_dist cond_dist ./ sum(cond_dist, 2); end这个权重矩阵就是条件分布的经验估计。要取某个分位点对每行按Y_train排序后做加权累积即可。注意Y_train需要是原始尺度如果训练时用了对数收入这里要还原。3. 基尼系数与收入差距的反事实分解从条件分布到分配效应3.1 用条件分布计算反事实基尼系数的完整流程有了 QRF 输出的条件分布就可以做反事实分析。核心逻辑是如果要消除某个因素比如教育差异的影响收入分布会变成什么样具体做法是把所有样本的该特征设为同一值如均值或中位数重新预测条件分布再从中抽样生成反事实收入计算基尼系数。function gini computeGini(income) % computeGini 计算基尼系数 income sort(income); n length(income); cum_income cumsum(income); gini 1 - 2 * sum(cum_income) / (n * sum(income)) 1/n; end function counterfactual_income generateCounterfactual(cond_dist, Y_train, n_samples) % generateCounterfactual 从条件分布中抽样生成反事实收入 n size(cond_dist, 1); counterfactual_income zeros(n, 1); for i 1:n % 按权重抽样 idx randsample(length(Y_train), n_samples, true, cond_dist(i, :)); counterfactual_income(i) mean(Y_train(idx)); end end逻辑说明computeGini用的是标准基尼系数公式的离散形式。generateCounterfactual从每个样本的条件分布中做加权抽样取均值作为该样本的反事实收入。这里n_samples建议取 100-500太小则抽样噪声大太大则计算慢。参数说明反事实分析的关键在于“固定哪个变量”。比如要消除教育的影响就把所有样本的受教育年限设为样本均值其他特征保持不变重新走一遍 QRF 预测。基尼系数的变化量就是该因素对收入差距的贡献。3.2 多维度风险因素的效应分解与结果解读把上述流程包装成循环对每个风险维度做反事实就能得到分解结果。下面是一个完整的分解框架% 假设 X 是特征矩阵Y 是对数收入feature_names 是特征名 risk_dims {education, industry, hukou, region, occupation_stability}; gini_baseline computeGini(exp(Y)); % 基准基尼 for d 1:length(risk_dims) X_cf X; % 将该维度设为均值连续变量或众数分类变量 col_idx find(strcmp(feature_names, risk_dims{d})); if iscontinuous(X(:, col_idx)) X_cf(:, col_idx) mean(X(:, col_idx)); else X_cf(:, col_idx) mode(X(:, col_idx)); end % 重新预测条件分布 cond_dist_cf predictQRF(qr_models, leaf_samples, X_cf, Y); income_cf generateCounterfactual(cond_dist_cf, exp(Y), 200); gini_cf computeGini(income_cf); fprintf(%s 的贡献: %.4f\n, risk_dims{d}, gini_baseline - gini_cf); end结果解读时要注意基尼系数下降幅度越大说明该因素对收入差距的贡献越大。但这里有个常见误解——反事实分析得到的是“关联性贡献”而非“因果效应”。如果要做因果推断需要额外的识别策略如工具变量、双重差分QRF 本身只解决分布估计问题。4. 避坑与排查QRF 做收入分配分析时最容易翻车的五个地方4.1 叶节点样本太少导致分位数估计震荡现象不同随机种子跑出来的基尼系数差异超过 0.02分位数曲线锯齿严重。原因MinLeafSize设得太小比如 1 或 2每个叶节点只有几个样本条件分布估计方差极大。收入数据本身噪声就大叶节点样本少时 QRF 退化成最近邻估计。解决把MinLeafSize调到 10-20同时增加nTrees到 1000。如果样本量本身很小少于 2000考虑先做特征降维或增加正则化。4.2 对数变换后忘记还原导致基尼系数算错现象基尼系数算出来是负数或者大于 1。原因QRF 训练时用了log(income)预测出来的条件分布是对数尺度的直接拿去算基尼系数就错了。基尼系数要求收入是原始尺度。解决在generateCounterfactual里对抽样结果做exp()还原。注意对数收入的均值不等于原始收入均值的对数所以不能先取对数均值再exp必须对每个抽样值单独还原后再取均值。4.3 分类变量当成连续变量处理现象行业、户籍等分类变量在反事实分析中设为“均值”后结果完全不可解释。原因MATLAB 的TreeBagger对分类变量需要显式声明CategoricalPredictors否则会当成连续变量做分裂。行业代码 1-20 被当成数值大小分裂逻辑就错了。解决训练前用categorical()转换分类变量并在TreeBagger中指定CategoricalPredictors参数。反事实时分类变量取众数而非均值。4.4 反事实抽样次数太少导致结果不稳定现象同一个反事实场景跑两次基尼系数差 0.01 以上。原因generateCounterfactual里n_samples设得太小比如 10抽样噪声淹没了真实效应。解决n_samples至少取 100建议 200-500。如果计算资源允许做 10 次重复抽样取均值进一步降低噪声。4.5 特征间高度相关导致分解结果重叠现象教育和职业稳定性的贡献加起来超过基准基尼系数或者某个因素单独看贡献很大但加入其他因素后贡献骤降。原因流动人口数据中教育、职业、行业、收入之间高度相关。反事实分析每次只固定一个变量但其他相关变量还在变导致贡献重叠。解决这是反事实分解的固有局限。可以补充做“序贯分解”——按一定顺序逐个固定变量看边际贡献。或者用 Shapley 值方法做公平分配但计算量会大很多。5. 让 QRF 收入分配分析更稳的几个进阶技巧第一个技巧是分位数交叉验证。不要只看 OOB 误差而是对每个分位点0.1、0.25、0.5、0.75、0.9分别计算预测分位数与实际分位数的偏差。具体做法把样本分成 K 折每折用训练集训练 QRF在验证集上预测各分位点计算实际覆盖率。如果 0.9 分位点的实际覆盖率只有 0.8说明模型在高分位点欠拟合需要增加nTrees或调整MinLeafSize。第二个技巧是变量重要性按分位点分解。标准随机森林的变量重要性是全局的但收入分配研究更关心“哪个因素在高收入端更重要”。做法是对每个分位点计算该分位点预测值对每个特征的偏依赖然后看偏依赖曲线的斜率。MATLAB 没有现成函数需要自己写循环固定其他特征让目标特征在分位数范围内变化观察预测分位点的变化幅度。第三个技巧是基尼系数的 bootstrap 置信区间。反事实分析得到的基尼系数是一个点估计没有不确定性度量。做法是对样本做 200 次 bootstrap 重抽样每次重跑 QRF 和反事实分解得到基尼系数贡献的分布取 2.5% 和 97.5% 分位数作为置信区间。这个计算量很大200 次 QRF 训练但这是让论文结果经得起审稿人质疑的必要步骤。% Bootstrap 置信区间示例 n_boot 200; gini_contrib zeros(n_boot, length(risk_dims)); for b 1:n_boot idx randsample(n, n, true); X_boot X(idx, :); Y_boot Y(idx); [qr_boot, leaf_boot] trainQRF(X_boot, Y_boot, 500, 15, 8); % ... 重复反事实分解流程 gini_contrib(b, :) ...; end % 计算 95% 置信区间 ci_lower prctile(gini_contrib, 2.5, 1); ci_upper prctile(gini_contrib, 97.5, 1);这里n_boot取 200 是底线500 更稳。每次 bootstrap 的nTrees可以降到 300 以节省时间因为 bootstrap 本身的重复已经提供了稳定性。最后一个习惯每次跑完分解先把基准基尼系数和文献里的同类研究对比。如果流动人口样本的基尼系数在 0.35-0.45 之间说明数据和处理基本合理如果低于 0.25 或高于 0.55大概率是数据处理出了问题比如收入没有做通胀调整、极端值没处理、或者样本筛选有偏。这个检查花不了几分钟但能避免后面所有分析建立在错误基础上。希望帮到你。本文还有配套的精品资源点击获取