Mata面板数据线性规划问题：个体循环与目标函数设定报错

需求回顾 #

基于有序probit模型得到的个体教育类别概率，通过线性规划为每个个体分配低（z）、中（s）、高（v）教育类别，需满足：

• 每个个体仅分配一个类别（Xz,Xs,Xv为0/1，且和为1）
• 人群中各教育类别的总人数符合给定的边际约束（∑Xz=Qz，∑Xs=Qs，∑Xv=Qv）
• 目标是最大化所有个体分配类别的概率总和

原代码运行时出现错误：invalid dimension of coefficients The coefficients must be a rowvector with length at least 1.，核心问题在于模型构建逻辑和代码实现的多处错误。

关键错误分析 #

1. 系数维度不匹配：
   原代码中coefs = (Pz, Ps, Pv)直接使用了整个n×3的概率矩阵，而每个线性规划实例需要对应变量的1维行向量。即使按个体循环，也应该取当前个体的概率行向量P[i,.]。

2. 全局约束处理逻辑错误：
   原代码试图在每个个体的LP中加入全局总和约束，这完全不合理——单个个体的LP只能控制自身的变量，无法影响其他个体的选择，全局边际约束必须在包含所有个体变量的统一LP模型中定义。

3. 矩阵维度倒置：
   X = J(m, n, .)将类别数作为行数、个体数作为列数，与需求（每行对应一个个体）不符，应改为X = J(n, m, .)。

4. LP对象未重置：
   循环中重复使用同一个lp对象，未清空之前的约束设置，导致约束累积引发维度错误。

修正后的代码实现 #

正确的做法是构建一个包含所有个体变量的LP模型，将每个个体的三个类别选择作为变量，统一设置目标函数和约束：

clear
input float(Pz Ps Pv Qz Qs Qv) 
 .0226813 .8045667    .172752 3 5 2
.09679342 .8531785  .05002809 3 5 2
.06026781 .8577851   .0819471 3 5 2
.15252922 .8199771 .027493654 3 5 2
.24786794 .7403268 .011805267 3 5 2
 .0226813 .8045667    .172752 3 5 2
.09538278 .8537293   .0508879 3 5 2
.05995392 .8576999  .08234616 3 5 2
.08492604 .8570921  .05798185 3 5 2
.09858042 .8524512  .04896835 3 5 2
end

* 提取Q的标量值（所有行Q相同，取第一行）
local Qz = Qz[1]
local Qs = Qs[1]
local Qv = Qv[1]

mata
putmata Pz Ps Pv, replace
P = (Pz, Ps, Pv)
n = rows(P)          // 个体数
m = cols(P)          // 类别数
total_vars = n*m     // 总变量数：每个个体3个变量

// 初始化LP对象
lp = LinearProgram()
lp.setMaxOrMin("max")

// 设置目标函数系数：每个个体的三个概率依次排列
coefs = J(1, total_vars, .)
for (i=1; i<=n; i++) {
    coefs[1, (i-1)*m + 1 :: i*m] = P[i,.]
}
lp.setCoefficients(coefs)

// 设置约束1：每个个体的三个类别变量和为1（共n个约束）
eq_lhs = J(n, total_vars, 0)
eq_rhs = J(n, 1, 1)
for (i=1; i<=n; i++) {
    eq_lhs[i, (i-1)*m + 1 :: i*m] = 1
}
lp.setEquality(eq_lhs, eq_rhs)

// 设置约束2：全局边际分布约束（共3个约束）
margin_lhs = J(m, total_vars, 0)
margin_rhs = (`Qz', `Qs', `Qv')'
for (k=1; k<=m; k++) {
    margin_lhs[k, k :: m :: total_vars] = 1
}
lp.setEquality(margin_lhs, margin_rhs)

// 设置变量边界：所有变量0<=X<=1
lower = J(1, total_vars, 0)
upper = J(1, total_vars, 1)
lp.setBounds(lower, upper)

// 求解LP
lp.optimize()

// 提取结果并整理为n×3的X矩阵
sol = lp.parameters()
X = J(n, m, .)
for (i=1; i<=n; i++) {
    X[i,.] = sol[(i-1)*m + 1 :: i*m]
}

// 由于LP得到的是连续解，需转换为0/1整数解（这里采用四舍五入，可根据需求调整）
X = round(X)

// 验证约束是否满足
sum_X = colsum(X)
disp("实际边际分布：")
disp(sum_X)

// 导出到Stata
st_matrix("X_matrix", X)
end

* 查看结果
matrix list X_matrix

补充说明 #

• 上述代码中，线性规划得到的是连续解，通过round(X)转换为0/1整数解，若需要严格满足整数约束，可考虑使用Stata的整数规划命令（如intprog）或Mata的整数规划扩展。
• 验证步骤可确认实际边际分布是否与给定的Qz/Qs/Qv一致，若存在偏差，可调整整数化方式（如贪心算法选择概率最高的类别，同时调整满足边际约束）。

内容的提问来源于stack exchange，提问作者Katarina V