dimanche 24 mai 2020

Inconsistent value returned from C++ function vs Python function for skew normal distribution

I have a function to estimate the alpha parameter from the skew normal distribution in C++ and Python. The Python function is written using NumPy and the C++ function uses the STL. My issue is that my C++ implementation is giving me incorrect results. The two functions are essentially identical but the Python version gives me correct results whereas the C++ does not - I have investigated this in some detail and I cannot come to a conclusion as to what's causing the error, any help would be great.

Python Function

import numpy as np
def convert_to_alpha(skew):
    a = np.pi/2 
    skew_ = abs(skew)
    numerator = np.power(skew_, (2/3)) 
    b = (4-np.pi)/2
    b = np.power(b, (2/3))
    denom = numerator + b
    delta = np.sqrt(a * (numerator/denom))
    a = delta/np.sqrt((1-np.power(delta, 2)))
    return a * np.sign(skew)

C++ Function

double convert_to_alpha(double skew)
{
    double pi = 3.141592653589793;
    double a = pi / 2;
    double skew_ = std::abs(skew);
    double numerator = std::pow(skew_, (2 / 3));
    double b = (4 - pi) / 2;
    b = std::pow(b, (2 / 3));
    double denom = numerator + b;
    double delta = std::sqrt(a * (numerator / denom));
    double alpha = delta / std::sqrt((1 - std::pow(delta, 2)));
    if (skew == 0) { return 0; }
    else if (std::signbit(skew) == 1) { return -1 * alpha; } 
    else return alpha; 
}

The Python function returns the values I would expect whereas the C++ function does not, as examples for input 0.99 I'd expect 27.85xxxx or for input 0.5 I'd expect 2.17xxxx which is exactly what I get from the Python implementation, C++ gives me 1.91306.

Also, strangely - regardless of the input, the C++ implementation seems to return 1.91306.

Driver code for C++

#include <cmath>
#include <math.h>
#include <iostream>

int main()
{
    double convert_to_alpha(double skew);
    std::cout << "skew: " << convert_to_alpha(0.99);
    return 0;
}
double convert_to_alpha(double skew)
{
    double pi = 3.141592653589793;
    double a = pi / 2;
    double skew_ = std::abs(skew);
    double numerator = std::pow(skew_, (2 / 3));
    double b = (4 - pi) / 2;
    b = std::pow(b, (2 / 3));
    double denom = numerator + b;
    double delta = std::sqrt(a * (numerator / denom));
    double alpha = delta / std::sqrt((1 - std::pow(delta, 2)));
    if (skew == 0) { return 0; } // if skew is 0 return 0 
    else if (std::signbit(skew) == 1) { return -1 * alpha; } // if skew is negative return -alpha
    else return alpha; // if skew is positive return alpha
}

I'd expect the results to be very similar, definitely not as different as they are currently. I have not encountered an issue like this before so any help figuring out what's causing the inconsistency with the C++ implementation would be very helpful.

Aucun commentaire:

Enregistrer un commentaire