如何 select 从两个拟合优度测试中最佳拟合连续分布？
how to select best fit continuous distribution from two Goodness-to-fit tests?
我调查了问题 Best fit Distribution plots，发现提交的答案使用 Kolmogorov-Smirnov 检验来寻找最佳拟合分布。我还发现有一个 Anderson-Darling 测试也用于获得最佳拟合分布。所以，我有几个问题：
问题一：
如果我有数据并通过 NumPy 直方图传递它，我应该使用什么参数以及我应该将什么输出输入到分布中？
def get_hist(data, data_size):
#### General code:
bins_formulas = ['auto', 'fd', 'scott', 'rice', 'sturges', 'doane', 'sqrt']
# bins = np.histogram_bin_edges(a=data, bins='scott')
# bins = np.histogram_bin_edges(a=data, bins='auto')
bins = np.histogram_bin_edges(a=data, bins='fd')
# print('Bin value = ', bins)
# Obtaining the histogram of data:
# Hist, bin_edges = histogram(a=data, bins=bins, range=np.linspace(start=np.min(data),end=np.max(data),size=data_size), density=True)
# Hist, bin_edges = histogram(a=data, range=np.linspace(np.min(data), np.max(data), data_size), density=True)
# Hist, bin_edges = histogram(a=data, bins=bins, density=True)
# Hist, bin_edges = histogram(a=data, bins=bins, range=(min(data), max(data)), normed=True, density=True)
# Hist, bin_edges = histogram(a=data, density=True)
Hist, bin_edges = histogram(a=data, range=(min(data), max(data)), density=True)
return Hist
问题二：
如果我想结合这两个测试，我该怎么做？哪些参数最适合用于寻找最佳拟合分布？这是我结合这两个测试的尝试。
from statsmodels.stats.diagnostic import anderson_statistic as adtest
def get_best_distribution(data):
dist_names = ['alpha', 'anglit', 'arcsine', 'beta', 'betaprime', 'bradford', 'burr', 'cauchy', 'chi', 'chi2', 'cosine', 'dgamma', 'dweibull', 'erlang', 'expon', 'exponweib', 'exponpow', 'f', 'fatiguelife', 'fisk', 'foldcauchy', 'foldnorm', 'frechet_r', 'frechet_l', 'genlogistic', 'genpareto', 'genexpon', 'genextreme', 'gausshyper', 'gamma', 'gengamma', 'genhalflogistic', 'gilbrat', 'gompertz', 'gumbel_r', 'gumbel_l', 'halfcauchy', 'halflogistic', 'halfnorm', 'hypsecant', 'invgamma', 'invgauss', 'invweibull', 'johnsonsb', 'johnsonsu', 'ksone', 'kstwobign', 'laplace', 'logistic', 'loggamma', 'loglaplace', 'lognorm', 'lomax', 'maxwell', 'mielke', 'moyal', 'nakagami', 'ncx2', 'ncf', 'nct', 'norm', 'pareto', 'pearson3', 'powerlaw', 'powerlognorm', 'powernorm', 'rdist', 'reciprocal', 'rayleigh', 'rice', 'recipinvgauss', 'semicircular', 't', 'triang', 'truncexpon', 'truncnorm', 'tukeylambda', 'uniform', 'vonmises', 'wald', 'weibull_min', 'weibull_max', 'wrapcauchy']
dist_ks_results = []
dist_ad_results = []
params = {}
for dist_name in dist_names:
dist = getattr(st, dist_name)
param = dist.fit(data)
params[dist_name] = param
# Applying the Kolmogorov-Smirnov test
D_ks, p_ks = st.kstest(data, dist_name, args=param)
print("Kolmogorov-Smirnov test Statistics value for " + dist_name + " = " + str(D_ks))
# print("p value for " + dist_name + " = " + str(p_ks))
dist_ks_results.append((dist_name, p_ks))
# Applying the Anderson-Darling test:
D_ad = adtest(x=data, dist=dist, fit=False, params=param)
print("Anderson-Darling test Statistics value for " + dist_name + " = " + str(D_ad))
dist_ad_results.append((dist_name, D_ad))
print(dist_ks_results)
print(dist_ad_results)
for D in range (len(dist_ks_results)):
KS_D = dist_ks_results[D][1]
AD_D = dist_ad_results[D][1]
if KS_D < 0.25 and AD_D < 0.05:
best_ks_D = KS_D
best_ad_D = AD_D
if dist_ks_results[D][1] == best_ks_D:
best_ks_dist = dist_ks_results[D][0]
if dist_ad_results[D][1] == best_ad_D:
best_ad_dist = dist_ad_results[D][0]
print(best_ks_D)
print(best_ad_D)
print(best_ks_dist)
print(best_ad_dist)
print('\n################################ Kolmogorov-Smirnov test parameters #####################################')
print("Best fitting distribution (KS test): " + str(best_ks_dist))
print("Best test Statistics value (KS test): " + str(best_ks_D))
print("Parameters for the best fit (KS test): " + str(params[best_ks_dist])
print('################################################################################\n')
print('################################ Anderson-Darling test parameters #########################################')
print("Best fitting distribution (AD test): " + str(best_ad_dist))
print("Best test Statistics value (AD test): " + str(best_ad_D))
print("Parameters for the best fit (AD test): " + str(params[best_ad_dist]))
print('################################################################################\n')
问题三：
如何获取 Anderson-Darling 检验的 p 值？
问题四：
假设我设法获得了最佳拟合分布，如何根据测试对分布进行排名？喜欢下面的照片。
Goodness-to-fit tests with ranking
编辑 1
我不确定，但是来自 statsmodel 的 normal_ad 是针对任何连续概率分布的一般 Anderson-Darling 检验吗？如果是，我想 select 两个测试通用的分布，如果我按照问题 1 中的相同步骤操作，这是否是正确的方法？另外，如果我想找到最高的 p 值并且在两个测试中都很常见，我如何提取具有 p 值的公共分布名称？
def get_best_distribution(data):
dist_names = ['beta', 'bradford', 'burr', 'cauchy', 'chi', 'chi2', 'erlang', 'expon', 'f', 'fatiguelife', 'fisk', 'gamma', 'genlogistic', 'genpareto', 'invgauss', 'johnsonsb', 'johnsonsu', 'laplace', 'logistic', 'loggamma', 'loglaplace', 'lognorm', 'maxwell', 'mielke', 'norm', 'pareto', 'reciprocal', 'rayleigh', 't', 'triang', 'uniform', 'weibull_min', 'weibull_max']
dist_ks_results = []
dist_ad_results = []
params = {}
for dist_name in dist_names:
dist = getattr(st, dist_name)
param = dist.fit(data)
params[dist_name] = param
# Applying the Kolmogorov-Smirnov test
D_ks, p_ks = st.kstest(data, dist_name, args=param)
print("Kolmogorov-Smirnov test Statistics value for " + dist_name + " = " + str(D_ks))
print("p value (KS test) for " + dist_name + " = " + str(p_ks))
dist_ks_results.append((dist_name, p_ks))
# Applying the Anderson-Darling test:
D_ad, p_ad = adnormtest(x=data, axis=0)
print("Anderson-Darling test Statistics value for " + dist_name + " = " + str(D_ad))
print("p value (AD test) for " + dist_name + " = " + str(p_ad))
dist_ad_results.append((dist_name, p_ad))
# select the best fitted distribution:
best_ks_dist, best_ks_p = (max(dist_ks_results, key=lambda item: item[1]))
best_ad_dist, best_ad_p = (max(dist_ad_results, key=lambda item: item[1]))
print('\n################################ Kolmogorov-Smirnov test parameters #####################################')
print("Best fitting distribution (KS test) :" + str(best_ks_dist))
print("Best p value (KS test) :" + str(best_ks_p))
print("Parameters for the best fit (KS test) :" + str(params[best_ks_dist]))
print('###########################################################################################################\n')
print('################################ Anderson-Darling test parameters #########################################')
print("Best fitting distribution (AD test) :" + str(best_ad_dist))
print("Best p value (AD test) :" + str(best_ad_p))
print("Parameters for the best fit (AD test) :" + str(params[best_ad_dist]))
print('###########################################################################################################\n')
if best_ks_dist == best_ad_dist:
best_common_dist = best_ks_dist
print('##################################### Both test parameters ############################################')
print("Best fitting distribution (Both test) :" + str(best_common_dist))
print("Best p value (KS test) :" + str(best_ks_p))
print("Best p value (AD test) :" + str(best_ad_p))
print("Parameters for the best fit (Both test) :" + str(params[best_common_dist]))
print('###########################################################################################################\n')
return best_common_dist, best_ks_p, params[best_common_dist]
问题 5:
如果我在实施拟合优度检验时出错，请纠正我，获得的 p 值用于检查给定值是否符合上述任何分布。因此，p 值的最大值意味着 p 值位于显着性水平的 %5 以下，因此，例如，Gamma 分布适合数据。我是对的还是没有理解拟合优度测试的主要概念？
第3题很容易用OpenTURNS解决。我通常使用贝叶斯信息准则对分布进行排名，因为它允许对参数较少的分布进行更好的排名。
在下面的示例中，我创建了一个高斯分布并从中生成了一个样本。然后我使用 FittingTest.BIC 函数计算库中 30 个分布的 BIC 分数。然后我使用 np.argsort 函数获取排序后的索引并打印结果。
import openturns as ot
import numpy as np
# Generate a sample
distribution = ot.Normal()
sample = distribution.getSample(100)
tested_factories = ot.DistributionFactory.GetContinuousUniVariateFactories()
nbmax = len(tested_factories)
# Compute BIC scores
bic_scores = []
names = []
for i in range(nbmax):
factory = tested_factories[i]
names.append(factory.getImplementation().getClassName())
try:
fitted_dist, bic = ot.FittingTest.BIC(sample, factory)
except:
bic = np.inf
bic_scores.append(bic)
# Sort the scores
indices = np.argsort(bic_scores)
# Print result
for i in range(nbmax):
factory = tested_factories[i]
name = factory.getImplementation().getClassName()
print(names[indices[i]], ": ", i, bic_scores[indices[i]])
这会产生：
NormalFactory : 0 2.902476153791324
TruncatedNormalFactory : 1 2.9391403094910493
LogisticFactory : 2 2.945101831314491
LogNormalFactory : 3 2.948479498106734
StudentFactory : 4 2.9487326727806438
WeibullMaxFactory : 5 2.9506160993704653
WeibullMinFactory : 6 2.9646030668970464
TriangularFactory : 7 2.9683050343363897
TrapezoidalFactory : 8 2.970676202179786
BetaFactory : 9 3.033244379700322
RayleighFactory : 10 3.0511170157342207
LaplaceFactory : 11 3.0641174552986796
FrechetFactory : 12 3.1472260896504327
UniformFactory : 13 3.1551588725784927
GumbelFactory : 14 3.1928562445001263
HistogramFactory : 15 3.3881831435932748
GammaFactory : 16 3.3925823197940552
ExponentialFactory : 17 3.824030948338899
ArcsineFactory : 18 214.7536151046246
ChiFactory : 19 680.8835152447839
ChiSquareFactory : 20 683.6769102883109
FisherSnedecorFactory : 21 inf
LogUniformFactory : 22 inf
GeneralizedParetoFactory : 23 inf
RiceFactory : 24 inf
DirichletFactory : 25 inf
BurrFactory : 26 inf
InverseNormalFactory : 27 inf
MeixnerDistributionFactory : 28 inf
ParetoFactory : 29 inf
有些分布无法适合此样本。在这些发行版中，我将 BIC 设置为 INF 并将异常包装在 try/except 中。
问题 2. 可以用 NormalityTest.AndersonDarlingNormal class:
来解决
import openturns as ot
distribution = ot.Normal()
sample = distribution.getSample(100)
test_result = ot.NormalityTest.AndersonDarlingNormal(sample)
print(test_result.getPValue())
这会打印：
0.8267360272974381
API 已记录 in the help page of the function, there is an example and the theory is documented here。
我调查了问题 Best fit Distribution plots，发现提交的答案使用 Kolmogorov-Smirnov 检验来寻找最佳拟合分布。我还发现有一个 Anderson-Darling 测试也用于获得最佳拟合分布。所以，我有几个问题：
问题一：
如果我有数据并通过 NumPy 直方图传递它，我应该使用什么参数以及我应该将什么输出输入到分布中？
def get_hist(data, data_size):
#### General code:
bins_formulas = ['auto', 'fd', 'scott', 'rice', 'sturges', 'doane', 'sqrt']
# bins = np.histogram_bin_edges(a=data, bins='scott')
# bins = np.histogram_bin_edges(a=data, bins='auto')
bins = np.histogram_bin_edges(a=data, bins='fd')
# print('Bin value = ', bins)
# Obtaining the histogram of data:
# Hist, bin_edges = histogram(a=data, bins=bins, range=np.linspace(start=np.min(data),end=np.max(data),size=data_size), density=True)
# Hist, bin_edges = histogram(a=data, range=np.linspace(np.min(data), np.max(data), data_size), density=True)
# Hist, bin_edges = histogram(a=data, bins=bins, density=True)
# Hist, bin_edges = histogram(a=data, bins=bins, range=(min(data), max(data)), normed=True, density=True)
# Hist, bin_edges = histogram(a=data, density=True)
Hist, bin_edges = histogram(a=data, range=(min(data), max(data)), density=True)
return Hist
问题二：
如果我想结合这两个测试，我该怎么做？哪些参数最适合用于寻找最佳拟合分布？这是我结合这两个测试的尝试。
from statsmodels.stats.diagnostic import anderson_statistic as adtest
def get_best_distribution(data):
dist_names = ['alpha', 'anglit', 'arcsine', 'beta', 'betaprime', 'bradford', 'burr', 'cauchy', 'chi', 'chi2', 'cosine', 'dgamma', 'dweibull', 'erlang', 'expon', 'exponweib', 'exponpow', 'f', 'fatiguelife', 'fisk', 'foldcauchy', 'foldnorm', 'frechet_r', 'frechet_l', 'genlogistic', 'genpareto', 'genexpon', 'genextreme', 'gausshyper', 'gamma', 'gengamma', 'genhalflogistic', 'gilbrat', 'gompertz', 'gumbel_r', 'gumbel_l', 'halfcauchy', 'halflogistic', 'halfnorm', 'hypsecant', 'invgamma', 'invgauss', 'invweibull', 'johnsonsb', 'johnsonsu', 'ksone', 'kstwobign', 'laplace', 'logistic', 'loggamma', 'loglaplace', 'lognorm', 'lomax', 'maxwell', 'mielke', 'moyal', 'nakagami', 'ncx2', 'ncf', 'nct', 'norm', 'pareto', 'pearson3', 'powerlaw', 'powerlognorm', 'powernorm', 'rdist', 'reciprocal', 'rayleigh', 'rice', 'recipinvgauss', 'semicircular', 't', 'triang', 'truncexpon', 'truncnorm', 'tukeylambda', 'uniform', 'vonmises', 'wald', 'weibull_min', 'weibull_max', 'wrapcauchy']
dist_ks_results = []
dist_ad_results = []
params = {}
for dist_name in dist_names:
dist = getattr(st, dist_name)
param = dist.fit(data)
params[dist_name] = param
# Applying the Kolmogorov-Smirnov test
D_ks, p_ks = st.kstest(data, dist_name, args=param)
print("Kolmogorov-Smirnov test Statistics value for " + dist_name + " = " + str(D_ks))
# print("p value for " + dist_name + " = " + str(p_ks))
dist_ks_results.append((dist_name, p_ks))
# Applying the Anderson-Darling test:
D_ad = adtest(x=data, dist=dist, fit=False, params=param)
print("Anderson-Darling test Statistics value for " + dist_name + " = " + str(D_ad))
dist_ad_results.append((dist_name, D_ad))
print(dist_ks_results)
print(dist_ad_results)
for D in range (len(dist_ks_results)):
KS_D = dist_ks_results[D][1]
AD_D = dist_ad_results[D][1]
if KS_D < 0.25 and AD_D < 0.05:
best_ks_D = KS_D
best_ad_D = AD_D
if dist_ks_results[D][1] == best_ks_D:
best_ks_dist = dist_ks_results[D][0]
if dist_ad_results[D][1] == best_ad_D:
best_ad_dist = dist_ad_results[D][0]
print(best_ks_D)
print(best_ad_D)
print(best_ks_dist)
print(best_ad_dist)
print('\n################################ Kolmogorov-Smirnov test parameters #####################################')
print("Best fitting distribution (KS test): " + str(best_ks_dist))
print("Best test Statistics value (KS test): " + str(best_ks_D))
print("Parameters for the best fit (KS test): " + str(params[best_ks_dist])
print('################################################################################\n')
print('################################ Anderson-Darling test parameters #########################################')
print("Best fitting distribution (AD test): " + str(best_ad_dist))
print("Best test Statistics value (AD test): " + str(best_ad_D))
print("Parameters for the best fit (AD test): " + str(params[best_ad_dist]))
print('################################################################################\n')
问题三：
如何获取 Anderson-Darling 检验的 p 值？
问题四：
假设我设法获得了最佳拟合分布，如何根据测试对分布进行排名？喜欢下面的照片。
Goodness-to-fit tests with ranking
编辑 1
我不确定，但是来自 statsmodel 的 normal_ad 是针对任何连续概率分布的一般 Anderson-Darling 检验吗？如果是，我想 select 两个测试通用的分布，如果我按照问题 1 中的相同步骤操作，这是否是正确的方法？另外，如果我想找到最高的 p 值并且在两个测试中都很常见，我如何提取具有 p 值的公共分布名称？
def get_best_distribution(data):
dist_names = ['beta', 'bradford', 'burr', 'cauchy', 'chi', 'chi2', 'erlang', 'expon', 'f', 'fatiguelife', 'fisk', 'gamma', 'genlogistic', 'genpareto', 'invgauss', 'johnsonsb', 'johnsonsu', 'laplace', 'logistic', 'loggamma', 'loglaplace', 'lognorm', 'maxwell', 'mielke', 'norm', 'pareto', 'reciprocal', 'rayleigh', 't', 'triang', 'uniform', 'weibull_min', 'weibull_max']
dist_ks_results = []
dist_ad_results = []
params = {}
for dist_name in dist_names:
dist = getattr(st, dist_name)
param = dist.fit(data)
params[dist_name] = param
# Applying the Kolmogorov-Smirnov test
D_ks, p_ks = st.kstest(data, dist_name, args=param)
print("Kolmogorov-Smirnov test Statistics value for " + dist_name + " = " + str(D_ks))
print("p value (KS test) for " + dist_name + " = " + str(p_ks))
dist_ks_results.append((dist_name, p_ks))
# Applying the Anderson-Darling test:
D_ad, p_ad = adnormtest(x=data, axis=0)
print("Anderson-Darling test Statistics value for " + dist_name + " = " + str(D_ad))
print("p value (AD test) for " + dist_name + " = " + str(p_ad))
dist_ad_results.append((dist_name, p_ad))
# select the best fitted distribution:
best_ks_dist, best_ks_p = (max(dist_ks_results, key=lambda item: item[1]))
best_ad_dist, best_ad_p = (max(dist_ad_results, key=lambda item: item[1]))
print('\n################################ Kolmogorov-Smirnov test parameters #####################################')
print("Best fitting distribution (KS test) :" + str(best_ks_dist))
print("Best p value (KS test) :" + str(best_ks_p))
print("Parameters for the best fit (KS test) :" + str(params[best_ks_dist]))
print('###########################################################################################################\n')
print('################################ Anderson-Darling test parameters #########################################')
print("Best fitting distribution (AD test) :" + str(best_ad_dist))
print("Best p value (AD test) :" + str(best_ad_p))
print("Parameters for the best fit (AD test) :" + str(params[best_ad_dist]))
print('###########################################################################################################\n')
if best_ks_dist == best_ad_dist:
best_common_dist = best_ks_dist
print('##################################### Both test parameters ############################################')
print("Best fitting distribution (Both test) :" + str(best_common_dist))
print("Best p value (KS test) :" + str(best_ks_p))
print("Best p value (AD test) :" + str(best_ad_p))
print("Parameters for the best fit (Both test) :" + str(params[best_common_dist]))
print('###########################################################################################################\n')
return best_common_dist, best_ks_p, params[best_common_dist]
问题 5:
如果我在实施拟合优度检验时出错，请纠正我，获得的 p 值用于检查给定值是否符合上述任何分布。因此，p 值的最大值意味着 p 值位于显着性水平的 %5 以下，因此，例如，Gamma 分布适合数据。我是对的还是没有理解拟合优度测试的主要概念？
第3题很容易用OpenTURNS解决。我通常使用贝叶斯信息准则对分布进行排名，因为它允许对参数较少的分布进行更好的排名。
在下面的示例中，我创建了一个高斯分布并从中生成了一个样本。然后我使用 FittingTest.BIC 函数计算库中 30 个分布的 BIC 分数。然后我使用 np.argsort 函数获取排序后的索引并打印结果。
import openturns as ot
import numpy as np
# Generate a sample
distribution = ot.Normal()
sample = distribution.getSample(100)
tested_factories = ot.DistributionFactory.GetContinuousUniVariateFactories()
nbmax = len(tested_factories)
# Compute BIC scores
bic_scores = []
names = []
for i in range(nbmax):
factory = tested_factories[i]
names.append(factory.getImplementation().getClassName())
try:
fitted_dist, bic = ot.FittingTest.BIC(sample, factory)
except:
bic = np.inf
bic_scores.append(bic)
# Sort the scores
indices = np.argsort(bic_scores)
# Print result
for i in range(nbmax):
factory = tested_factories[i]
name = factory.getImplementation().getClassName()
print(names[indices[i]], ": ", i, bic_scores[indices[i]])
这会产生：
NormalFactory : 0 2.902476153791324
TruncatedNormalFactory : 1 2.9391403094910493
LogisticFactory : 2 2.945101831314491
LogNormalFactory : 3 2.948479498106734
StudentFactory : 4 2.9487326727806438
WeibullMaxFactory : 5 2.9506160993704653
WeibullMinFactory : 6 2.9646030668970464
TriangularFactory : 7 2.9683050343363897
TrapezoidalFactory : 8 2.970676202179786
BetaFactory : 9 3.033244379700322
RayleighFactory : 10 3.0511170157342207
LaplaceFactory : 11 3.0641174552986796
FrechetFactory : 12 3.1472260896504327
UniformFactory : 13 3.1551588725784927
GumbelFactory : 14 3.1928562445001263
HistogramFactory : 15 3.3881831435932748
GammaFactory : 16 3.3925823197940552
ExponentialFactory : 17 3.824030948338899
ArcsineFactory : 18 214.7536151046246
ChiFactory : 19 680.8835152447839
ChiSquareFactory : 20 683.6769102883109
FisherSnedecorFactory : 21 inf
LogUniformFactory : 22 inf
GeneralizedParetoFactory : 23 inf
RiceFactory : 24 inf
DirichletFactory : 25 inf
BurrFactory : 26 inf
InverseNormalFactory : 27 inf
MeixnerDistributionFactory : 28 inf
ParetoFactory : 29 inf
有些分布无法适合此样本。在这些发行版中，我将 BIC 设置为 INF 并将异常包装在 try/except 中。
问题 2. 可以用 NormalityTest.AndersonDarlingNormal class:
import openturns as ot
distribution = ot.Normal()
sample = distribution.getSample(100)
test_result = ot.NormalityTest.AndersonDarlingNormal(sample)
print(test_result.getPValue())
这会打印：
0.8267360272974381
API 已记录 in the help page of the function, there is an example and the theory is documented here。