Python中Tobit回归报错ValueError: Input must be 1- or 2-d的排查

我用Python实现基于极大似然估计的Tobit回归模型，代码分为数据定义、估计器构建、模型估计三部分，但在提取参数标准误时触发ValueError: Input must be 1- or 2-d。以下是完整代码、报错堆栈及数据模拟代码：

核心实现代码 #

import numpy as np
from scipy.optimize import minimize


# 定义因变量和自变量
X = data.iloc[:, 1:]  
y = data.iloc[:, 0]

# 为自变量添加常数项列
X = np.c_[np.ones(X.shape[0]), X]

# 定义Tobit模型的似然函数
def likelihood(params, y, X, lower, upper):
    beta = params[:-1]
    sigma = params[-1]
    mu = X @ beta
    prob = (1 / (sigma * np.sqrt(2 * np.pi)) * np.exp(-0.5 * ((y - mu) / sigma)**2))
    prob[y < lower] = 0
    prob[y > upper] = 0
    return -np.log(prob).sum()

# 设置参数初始值和截尾上下界
params_init = np.random.normal(size=X.shape[1] + 1)
bounds = [(None, None) for i in range(X.shape[1])] + [(1e-10, None)]

# 执行MLE估计
res = minimize(likelihood, params_init, args=(y, X, 0, 100), bounds=bounds, method='L-BFGS-B')

# 提取估计参数和标准误
params = res.x
stderr = np.sqrt(np.diag(res.hess_inv))  # 报错行

# 打印结果
print(f'Coefficients: {params[:-1]}')
print(f'Standard Errors: {stderr[:-1]}')
print(f'Sigma: {params[-1]:.4f}')

报错信息 #

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-245-5f39f416cc07> in <module>
     31 # Extract the estimated parameters and their standard errors
     32 params = res.x
---> 33 stderr = np.sqrt(np.diag(res.hess_inv))
     34 
     35 # Print the results

/opt/anaconda3/lib/python3.8/site-packages/numpy/core/overrides.py in diag(*args, **kwargs)

/opt/anaconda3/lib/python3.8/site-packages/numpy/lib/twodim_base.py in diag(v, k)
    307         return diagonal(v, k)
    308     else:
---> 309         raise ValueError("Input must be 1- or 2-d.")
    310 
    311 

ValueError: Input must be 1- or 2-d.

数据模拟代码 #

import pandas as pd
import numpy as np

data = pd.DataFrame()

# 生成残疾与非残疾人群的面试概率
interview_prob_disabled = np.random.normal(38.63, 28.72, 619)
interview_prob_enabled = np.random.normal(44.27, 28.19, 542)
interview_prob = np.append(interview_prob_disabled, interview_prob_enabled)

# 修正变量范围并取整
interview_prob = np.clip(interview_prob, 0, 100)
interview_prob = np.round(interview_prob)

# 添加面试概率变量
data['Interview Probabilities'] = interview_prob

# 添加其他变量
data['Age'] = np.random.randint(18, 65, size=len(interview_prob))
data['Gender'] = np.random.choice(['Male', 'Female'], size=len(interview_prob))
data['Employment Status'] = np.random.choice(['Employed', 'Unemployed', 'Retired'], size=len(interview_prob))
data['Education Level'] = np.random.choice(['High School', 'College', 'Vocational', 'Graduate School'], size=len(interview_prob))

# 添加残疾状态虚拟变量
data['Disability Status'] = np.append(np.repeat('Disabled', 619), np.repeat('Non-disabled', 542))

# 转换分类变量为数值
data['Gender'] = data['Gender'].map({'Male': 0, 'Female': 1})
data['Employment Status'] = data['Employment Status'].map({'Employed': 0, 'Unemployed': 1})
data['Education Level'] = data['Education Level'].map({'High School': 0, 'College': 1, 'Vocational': 2, 'Graduate School': 3})
data['Disability Status'] = data['Disability Status'].map({'Disabled': 1, 'Non-disabled': 0})

print(data)

错误原因 #

你使用的L-BFGS-B优化方法返回的res.hess_inv不是常规的NumPy数组，而是**scipy.sparse.linalg.LinearOperator对象**——这是L-BFGS-B用来近似逆Hessian矩阵的内存高效型结构，并非二维数组。而np.diag()只能处理1维或2维数组，因此触发类型不匹配的错误。

解决方案 #

有两种可靠的方法提取标准误：

方法1：将LinearOperator转换为密集矩阵 #

调用LinearOperator的todense()方法将其转为二维数组，再提取对角线元素：

# 修正后的标准误提取代码
stderr = np.sqrt(np.diag(res.hess_inv.todense()))

方法2：通过矩阵乘法获取对角线元素 #

利用逆Hessian矩阵与单位矩阵的点积，直接提取对角线值（避免转换为密集矩阵，内存效率更高）：

# 修正后的标准误提取代码
n_params = len(params)
hess_inv_diag = res.hess_inv.dot(np.eye(n_params))[np.arange(n_params), np.arange(n_params)]
stderr = np.sqrt(hess_inv_diag)

额外注意：似然函数的正确性 #

虽然不是当前报错的原因，但你的Tobit似然函数存在逻辑缺陷：对于截尾样本（y <= lower或y >= upper），直接将概率设为0会导致np.log(0)出现无穷大值，正确的做法应该使用累积分布函数（CDF）计算截尾概率：

from scipy.stats import norm

def likelihood(params, y, X, lower, upper):
    beta = params[:-1]
    sigma = params[-1]
    mu = X @ beta
    
    # 非截尾样本的概率密度
    pdf = norm.pdf(y, loc=mu, scale=sigma)
    # 左截尾样本的累积概率（y <= lower）
    left_cdf = norm.cdf(lower, loc=mu, scale=sigma)
    # 右截尾样本的生存概率（y >= upper）
    right_sf = norm.sf(upper, loc=mu, scale=sigma)
    
    # 构建似然值
    prob = np.where(y < lower, left_cdf, 
                    np.where(y > upper, right_sf, pdf))
    return -np.log(prob).sum()

这个修正能保证似然函数在截尾样本处的数值稳定性，避免优化过程中出现异常。

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