import scipy.stats as ss
import numpy as np
#import statsmodels.api as sm
### Tests on residuals
[docs]def normal_Kolmogorov_Smirnov(sample):
"""The moon illumination expressed as a percentage.
:param astropy sun: the sun ephemeris
:param astropy moon: the moon ephemeris
:return: a numpy array like indicated the moon illumination.
:rtype: array_like
"""
mu, sigma = ss.norm.fit(sample)
#use mu sigma for anomaly, 0,1 for rescaling???
KS_stat, KS_pvalue = ss.kstest(sample, 'norm', args=(0, 1))
# the sample is likely Gaussian-like if KS_stat (~ maximum distance between sample and theoritical distribution) -> 0
# the null hypothesis can not be rejected ( i.e the distribution of sample come from a Gaussian) if KS_pvalue -> 1
KS_judgement = 0
if KS_pvalue > 0.01:
KS_judgement = 1
if KS_pvalue > 0.05:
KS_judgement = 2
return KS_stat, KS_pvalue, KS_judgement
[docs]def normal_Anderson_Darling(sample):
"""Compute a Anderson-Darling test on the sample versus a normal distribution with mu = 0, sigma = 1
:param array_like sample: the sample you want to check the "Gaussianity"
:returns: the Anderson-Darling statistic, the Anderson-Darling critical values associated to the significance
level of 15 % and the Anderson-Darling judgement
:rtype: float, array_like, array_like
"""
AD_stat, AD_critical_values, AD_significance_levels = ss.anderson(sample)
# the sample is likely Gaussian-like if AD_stat (~ maximum distance between sample and theoritical distribution) -> 0
# the null hypothesis can not be rejected ( i.e the distribution of sample come from a Gaussian) if AD_pvalue -> 1
AD_judgement = 0
if AD_stat < 2*AD_critical_values[-1]:
AD_judgement = 1
if AD_stat < AD_critical_values[-1]:
AD_judgement = 2
return AD_stat, AD_critical_values[-1], AD_judgement
[docs]def normal_Shapiro_Wilk(sample):
"""Compute a Shapiro-Wilk test on the sample versus a normal distribution with mu = 0, sigma = 1
:param array_like sample: the sample you want to check the "Gaussianity"
:returns: the Shapiro-Wilk statistic and its related p_value
:rtype: float, float
"""
SW_stat, SW_pvalue = ss.shapiro(sample)
# the null hypothesis can not be rejected ( i.e the distribution of sample come from a Gaussian) if SW_stat -> 1
# the null hypothesis can not be rejected ( i.e the distribution of sample come from a Gaussian) if SW_pvalue -> 1
# Judegement made on the STATISTIC because 'W test statistic is accurate but the p-value may not be" (see scipy doc)
SW_judgement = 0
if SW_pvalue > 0.01:
SW_judgement = 1
if SW_pvalue > 0.05:
SW_judgement = 2
return SW_stat, SW_pvalue, SW_judgement
### Statistics fit quality metrics
[docs]def normalized_chi2(chi2, n_data, n_parameters) :
"""Compute the chi^2/dof
:param float chi2: the chi^2
:param int n_data: the number of data_points
:param int n_parameters: the number of model parameters
:returns: the chi^2/dof and the chi2dof_judgement
:rtype: float
"""
chi2_sur_dof = chi2/(n_data-n_parameters)
chi2dof_judgement = 0
if chi2_sur_dof < 2 :
chi2dof_judgement = 2
return chi2_sur_dof,chi2dof_judgement