I wrote some python code to estimate a parameter using the Maximum Likelihood. I'm using the Newton-Raphson Method to solve the problem. However, I need to convert it to C and integrate it with the rest of the Software.

I am not familiar with C , how can I convert the following bock in Python to C :

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
    #Iterate until we find the solution within the desired epsilon
    while abs(fdc>epsilon or i<max_iter):
        i = i 1
        x1 = xhat - (fdc/sdc)
        xhat = x1
        fdc = first_derivative(x)
        print('The ML estimate of xhat is', xhat)
    return xhat

raphson_newton(xhat)

Given the toy example above, xhat should be around -0.5343967677954681.

I have tried the following but it's not converging to the same value. Not sure where I am getting it wrong.

#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;
}

CodePudding user response：

There are several issues in your C code:

• In first() and second() you need to iterate over elements of the vector y (or x as it is named in your Python code and thus also in my code below).
• In newton_raphson(), you are changing max_iter and there is no check if the result is already within the tolerance. Generally, the code can be made to resemble the Python code better.
• The iteration is started with the value 1 instead of start and the global variable xhat is not used.

There is still room for improvement, but the following should work:

#include <cmath>
#include <iostream>
#include <vector>

#include <cmath>

unsigned max_iter = 100;

std::vector<double> x = {-1.94, 0.59, -5.98, -0.08, -0.77};
double start = -0.77;

//Derivative of the function

double first(double xhat)
{
    double tfd = 0.0;
    for(auto &xi: x) tfd  = (xi - xhat) / (1   std::pow(xi - xhat, 2));
    double fd = 2 * tfd;
    return fd;
}

// Second derivative of the function

double second(double xhat)
{
    double tsd = 0.0;
    for(auto &xi: x) tsd  = (std::pow(xi - xhat, 2) - 1) / std::pow(1   std::pow(xi - xhat, 2), 2);
    double sd = 2 * tsd;
    return sd;
}

double newton_raphson(double xhat)
{
    double fdc = first(xhat);
    double sdc = second(xhat);
    double tolerance = 0.001;
    unsigned i = 0;
    // Iterate until we find a root within the desired tolerance
    while(i < max_iter && std::abs(fdc) > tolerance)
    {
        i  ;
        xhat -= fdc/sdc;
        fdc = first(xhat);
    }

    return xhat;
}

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