如何实现Stata的“nl”包（非线性最小二乘法）到Python

一名用户希望在Python中实现Stata nl（非线性最小二乘法）包。用户尝试使用lmfit和scipy.optimize.leastsq，但都无法正常工作，导致所有参数估计为 NaN。

解决方案

1. 检查输入数据和模型函数是否存在错误

   在使用最小二乘法拟合模型时，需要仔细检查输入数据和模型函数是否存在错误。特别是，检查是否出现了无法计算对数的问题（即x2-b2是否为负）。

2. 使用更鲁棒的优化算法

   lmfit 和 scipy.optimize.leastsq 都可能遇到数值不稳定或收敛缓慢的情况。用户可以选择使用更鲁棒的优化算法，例如 Nelder-Mead 或 Powell 算法。

3. 使用更合适的初始值

   用户给出的初始值 b1=10 和 b2=1990 可能离最优解太远，导致优化算法难以收敛。用户可以尝试使用更合适的初始值，例如通过绘图或使用其他方法来估计参数的初始值。

以下是用 Python 代码实现 Stata 的 nl 函数（非线性最小二乘法）的一个范例：

import numpy as np
from scipy.optimize import least_squares

def objective_function(params, x1, x2, x3, y):
    """
    目标函数（残差平方和）
    """
    b1, b2, b3, b5 = params
    model = x1 + b1 + 0.3*np.log(x2-b2)*b3 - 0.7*x3*b3 + b5*x2
    return model - y

def fit_nls(x1, x2, x3, y, initial_params):
    """
    非线性最小二乘法拟合
    """
    result = least_squares(objective_function, initial_params, args=(x1, x2, x3, y))
    return result

# 数据准备
x1 = np.array([1981, 1982, 1983, 1984, 1985, 1986, 1987, 1988, 1989, 1990,
              1991, 1992, 1993, 1994, 1995, 1996, 1997, 1998, 1999, 2000,
              2001, 2002])
x2 = np.array([15.2824955, 15.2635746, 15.30461597, 15.37490845, 15.41295338,
              15.44680405, 15.48135281, 15.52259159, 15.5565939, 15.57392025,
              15.57197666, 15.60479259, 15.63134575, 15.67116165, 15.69621944,
              15.7329874, 15.77698135, 15.81788635, 15.86141682, 15.89737129,
              15.90485096, 15.92070866])
x3 = np.array([14.56475067, 15.52343941, 16.30871582, 16.76519966, 17.04235458,
              17.25271797, 17.44781876, 17.62217331, 17.71343422, 17.81187439,
              17.89474106, 17.98217583, 18.06685829, 18.16578865, 18.27449799,
              18.38712311, 18.50685883, 18.63579178, 18.76427078, 18.89691544,
              18.99729347, 19.06253433])
y = np.array([2.936807632, 2.908272743, 2.940227509, 3.001846313, 3.030970573,
              3.055702209, 3.081344604, 3.113491058, 3.138068199, 3.144176483,
              3.128887177, 3.14837265, 3.161927223, 3.18959713, 3.202876091,
              3.228042603, 3.260077477, 3.289312363, 3.321393967, 3.34650898,
              3.344522476, 3.351119995])

# 设置初始参数
initial_params = np.array([10, 1990, 5, 12])

# 拟合模型
result = fit_nls(x1, x2, x3, y, initial_params)

# 打印拟合结果
print('拟合结果:')
print('参数估计值:')
print(result.x)
print('残差平方和:')
print(result.fun)

这个示例使用 scipy.optimize.least_squares() 函数进行非线性最小二乘法拟合，并在输出中打印拟合结果。