Python转C++：牛顿-拉夫森法参数估计代码不收敛求助

我用Python实现了基于极大似然估计的参数求解，采用牛顿-拉夫森法完成计算，现在需要将代码转换为C以集成到现有软件中，但我不熟悉C，自己写的C代码无法收敛到正确结果（正确结果约为 -0.5343967677954681），请帮忙排查问题并给出正确的C实现。

原Python实现代码 #

import numpy as np

x = np.array([-1.94, 0.59, -5.98, -0.08, -0.77])
start = np.median(x)
xhat = start
max_iter = 20
epsilon = 0.001

def first_derivative(xhat):
    fd = 2 * sum((x - xhat) / (1 + (x - xhat)**2))
    return fd

def second_derivative(xhat):
    sd = 2 * sum((((x - xhat)**2) - 1) / ((1 + (x - xhat)**2)**2))
    return sd

def raphson_newton(xhat):
    fdc = first_derivative(xhat)
    sdc = second_derivative(xhat)
    xhat = start
    i = 0
    # 迭代直到满足精度要求或达到最大迭代次数
    while abs(fdc) > epsilon and i < max_iter:
        i += 1
        x1 = xhat - (fdc / sdc)
        xhat = x1
        fdc = first_derivative(xhat)
        sdc = second_derivative(xhat)
        print('当前ML估计值xhat为', xhat)
    return xhat

result = raphson_newton(xhat)
print('最终ML估计值xhat为', result)

注：原Python代码中循环条件存在小错误，已修正为abs(fdc) > epsilon and i < max_iter，确保逻辑正确。

我尝试的C++代码（无法收敛） #

#include <cmath>
#include <iostream>
#include <vector>
using namespace std;

#include <cmath>

double max_iter = 100;
double start = -0.77;
double xhat = start;

vector<double> y = {-1.94, 0.59, -5.98, -0.08, -0.77};

//Derivative of the function

double first(double y)
{
    double tfd = (y - xhat) / (1 + pow(y - xhat, 2));
    double fd = 2 * tfd;
    return fd;
}

// Second derivative of the function

double second(double y)
{
    double tsd = (pow(y - xhat, 2) - 1) / pow(1 + pow(y - xhat, 2), 2);
    double sd = 2 * tsd;
    return sd;
}

double newton_raphson(double xhat)
{
    double tolerance = 0.001;
    double x1;
    int i = 0;
    // Iterate until we find a root within the desired tolerance
    do
    {
        double x1 = xhat - first(xhat) / second(xhat);
        xhat = x1;
        max_iter= i++;
    } while ( i < max_iter);

    return double (xhat);
}

int main()
{
    double xhat = newton_raphson(1);
    cout << "xhat: " << xhat << endl;
    return 0;
}

问题分析 #

你的C++代码存在多处逻辑错误，导致无法收敛：

• 导数计算错误：原Python代码中的一阶、二阶导数是对整个数组的项求和，而你的C++代码只计算单个元素的项，没有完成求和操作。
• 全局变量滥用：全局变量xhat导致计算时上下文混乱，迭代过程中无法正确传递当前的估计值。
• 迭代逻辑错误：循环中错误修改max_iter变量导致迭代次数失控；每次迭代没有重新计算当前xhat对应的一阶、二阶导数；循环条件未检查导数精度要求，只限制了迭代次数。
• 初始值错误：main函数中传递初始值为1，而正确的初始值应该是数组的中位数（-0.77）。

修正后的C++实现 #

#include <cmath>
#include <iostream>
#include <vector>
#include <algorithm> // 用于计算中位数

using namespace std;

// 计算一阶导数：对整个数组求和
double first_derivative(double xhat, const vector<double>& data) {
    double sum = 0.0;
    for (double y : data) {
        double diff = y - xhat;
        sum += diff / (1 + pow(diff, 2));
    }
    return 2 * sum;
}

// 计算二阶导数：对整个数组求和
double second_derivative(double xhat, const vector<double>& data) {
    double sum = 0.0;
    for (double y : data) {
        double diff = y - xhat;
        double diff_sq = pow(diff, 2);
        double denom = pow(1 + diff_sq, 2);
        sum += (diff_sq - 1) / denom;
    }
    return 2 * sum;
}

double newton_raphson(double initial_xhat, const vector<double>& data) {
    const int max_iter = 20;
    const double epsilon = 0.001;
    double xhat = initial_xhat;
    int i = 0;
    double fdc = first_derivative(xhat, data);
    
    while (abs(fdc) > epsilon && i < max_iter) {
        i++;
        double sdc = second_derivative(xhat, data);
        // 防止分母为0的情况
        if (abs(sdc) < 1e-8) {
            cout << "二阶导数接近0，无法继续迭代" << endl;
            break;
        }
        double x1 = xhat - fdc / sdc;
        xhat = x1;
        fdc = first_derivative(xhat, data);
        cout << "当前ML估计值xhat为: " << xhat << endl;
    }
    
    return xhat;
}

// 计算数组的中位数
double calculate_median(vector<double> data) {
    sort(data.begin(), data.end());
    int n = data.size();
    if (n % 2 == 1) {
        return data[n / 2];
    } else {
        return (data[n/2 - 1] + data[n/2]) / 2.0;
    }
}

int main() {
    vector<double> data = {-1.94, 0.59, -5.98, -0.08, -0.77};
    double initial_xhat = calculate_median(data);
    double result = newton_raphson(initial_xhat, data);
    cout << "最终ML估计值xhat为: " << result << endl;
    return 0;
}

代码说明 #

• 移除全局变量：将数据数组和估计值作为函数参数传递，避免上下文混乱。
• 实现求和逻辑：一阶、二阶导数函数遍历整个数组，计算所有项的和，与Python逻辑一致。
• 修正迭代逻辑：同时检查导数绝对值是否小于精度阈值和是否达到最大迭代次数；每次迭代重新计算当前xhat对应的一阶、二阶导数；添加分母为0的判断，防止程序崩溃。
• 自动计算中位数：新增calculate_median函数，自动计算数组的中位数作为初始值，无需手动指定。

运行修正后的代码，会输出与Python一致的结果：-0.534397（保留六位小数）。

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