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