英文:
Finding Percentiles and Values From Calculated Gamma Distribution
问题
Here is the translated text from your provided content:
背景
我正在使用Numpy（ndarray）计算一系列最佳拟合伽马曲线的工作，这是一个关于这个问题起源的先前问题，可以在这里找到。
以前使用过Scipy（scipy.gamma.stats），但是该库并不针对多维数组进行优化，因此编写了一个最基本的函数来满足目标。我已经成功地拟合了一个曲线（尽管不像Scipy那样干净），该曲线如下所示。
当前问题
我想要获取给定值的百分位数，反之亦然，沿着计算得到的伽马分布。但是，我得到的拟合曲线上的期望值与预期值不符。例如，提供50th百分位数会得到一个值为4.471的值，这与下面显示的曲线拟合不符。可以进行哪些修改或整体更改以从提供的数据中得出百分位数和值？
图表
代码
import sys, os, mathimport numpy as npimport scipy as sciimport matplotlib.pyplot as pltdata = np.array([0.00, 0.00, 11.26399994, 0.00, 0.00, 0.00, 0.00, 0.00, ... ]) # 这里是你的数据def scigamma(data):# 这是你的函数实现return cdfscicdf = scigamma(data)# 这是更多的代码，包括伽马分布拟合和绘图等。
请注意，这是您提供的内容的翻译，其中包括代码部分。如果您需要任何进一步的帮助或有其他问题，请随时提问。<details><summary>英文:</summary>**Background**I am working on computing a series of best-fit gamma curves for a 2-D dataset in Numpy (ndarray), a prior question for the genesis of this can be found [here][1].Scipy was previously utilized (scipy.gamma.stats), however this library is not optimized for multi-dimensional arrays and a barebones function was written to meet the objective. I've successfully fit (albeit not as cleanly as Scipy) a curve to the dataset, which is provided below.**Current Issue**I want to obtain the percentile of a given value and vice versa along the calculated gamma distribution. However, I'm not obtaining expected values off the fitted curve. For example, providing the 50th percentile yields a value of 4.471, which does not match up with the curve fit shown below. What modifications or wholesale alterations can be made to yield both percentiles and values from supplied data?**Graph**[![Gamma Fitting Output][2]][2]**Code**import sys, os, mathimport numpy as npimport scipy as sciimport matplotlib.pyplot as pltdata = np.array([0.00, 0.00, 11.26399994, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 17.06399918, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 8.33279991, 0.00, 7.54879951, 0.00, 0.00, 0.00, 4.58799982, 7.9776001, 0.00, 0.00, 0.00, 0.00, 11.45040035, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 18.73279953, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 8.94559956, 0.00, 7.73040009, 0.00, 0.00, 0.00, 5.03599977, 8.62639999, 0.00, 0.00, 0.00, 0.00, 11.11680031, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 14.37839985, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 8.16479969, 0.00, 7.30719948, 0.00, 0.00, 0.00, 3.41039991, 7.17280006, 0.00, 0.00, 0.00, 0.00, 10.0099199963, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 13.97839928, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 7.6855998, 0.00, 6.86559963, 0.00, 0.00, 0.00, 3.21600008, 7.93599987, 0.00, 0.00, 0.00, 0.00, 11.55999947, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 18.76399994, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 10.0033039951, 0.00, 8.10639954, 0.00, 0.00, 0.00, 4.76480007, 6.87679958, 0.00, 0.00, 0.00, 0.00, 11.42239952, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 19.42639732, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 10.0052400017, 0.00, 8.2567997, 0.00, 0.00, 0.00, 5.08239985, 7.9776001, 0.00, 0.00, 0.00, 0.00, 10.0099839973, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 11.5855999, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 7.88399982, 0.00, 5.96799994, 0.00, 0.00, 0.00, 3.07679987, 7.81360006, 0.00, 0.00, 0.00, 0.00, 11.51119995, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 20.0030959892, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 10.0050879955, 0.00, 8.20479965, 0.00, 0.00, 0.00, 5.51599979, 9.02879906, 0.00, 0.00])def scigamma(data):param = sci.stats.gamma.fit(data)x = np.linspace(0, np.max(data), 250)cdf = sci.stats.gamma.cdf(x, *param)value = np.round((sci.stats.gamma.cdf(0.10, *param) * 100), 2)percentile = np.round((sci.stats.gamma.ppf(50.00, *param) * 100), 2)return cdfscicdf = scigamma(data)# Method of Moments estimationmean = np.mean(data)variance = np.var(data)alpha = mean**2 / variancebeta = variance / mean# Generate x-axis values for the curvesx = np.linspace(0, np.max(data), 250)# Calculate the gamma distribution PDF valuespdf = (x ** (alpha - 1) * np.exp(-x / beta)) / (beta ** alpha * np.math.gamma(alpha))# Calculate the gamma distribution CDF valuescdf = np.zeros_like(x)cdf[x > 0] = np.cumsum(pdf[x > 0]) / np.sum(pdf[x > 0])# Estimate the probability of zero valuesnum_zeros = np.count_nonzero(data == 0)zero_probability = np.count_nonzero(data == 0) / len(data)# Calculate the PDF and CDF values at zeropdf_zero = zero_probability / (beta ** alpha * np.math.gamma(alpha))cdf_zero = zero_probabilityvalue = 2.50percentile = 0.50index = np.argmax(pdf >= value)# Calculate the percentile using numerical integrationpct = np.trapz(pdf[:index+1], dx=1) + (value - pdf[index]) * (cdf[index] - cdf[index-1]) / (pdf[index-1] - pdf[index])index = np.argmax(cdf >= percentile)# Calculate the value using numerical integrationval = np.trapz(cdf[:index+1], dx=1) + (percentile - cdf[index-1]) * (pdf[index] - pdf[index-1]) / (cdf[index] - cdf[index-1])# Plot the data histogramplt.hist(data, bins=30, density=True, alpha=0.5, label='data')# Plot the gamma distribution CDF curveplt.plot(x, cdf, 'b', label='Gamma CDF | Custom Fit')plt.plot(x, scicdf, 'k', label='Gamma CDF | SciPy Fit')# Set plot labels and legendplt.xlabel('data')plt.ylabel('Probability')plt.legend()[1]: https://stackoverflow.com/q/75598344/10266106[2]: https://i.stack.imgur.com/QchWe.png</details># 答案1**得分**: 1I'm not clear on why you don't trust `sci.stats`. If it's because of performance, have you actually profiled? It's worth noting that neither your fit nor scipy's fit are particularly good, because your data are strongly skewed toward 0.When initializing your second `x`, don't start it at 0 because your calculations are not well-defined for that point, which (a) dumps warnings on the console and (b) requires a bunch of indexing gymnastics later.I'm also not clear on why you're attempting `trapz`. Why integrate when you already have a CDF? Just use the CDF; and also don't `argmax`; use `searchsorted`.Finally, your plotting is slightly a disaster. For apples-to-apples comparison it's very important that you pass `cumulative` to `hist`.```pythonimport numpy as npimport scipy as sciimport matplotlib.pyplot as pltdata = np.array([0. , 0. , 11.26399994, 0. , 0. ,0. , 0. , 0. , 0. , 0. ,17.06399918, 0. , 0. , 0. , 0. ,0. , 0. , 8.33279991, 0. , 7.54879951,0. , 0. , 0. , 4.58799982, 7.9776001 ,0. , 0. , 0. , 0. , 11.45040035,0. , 0. , 0. , 0. , 0. ,0. , 0. , 18.73279953, 0. , 0. ,0. , 0. , 0. , 0. , 8.94559956,0. , 7.73040009, 0. , 0. , 0. ,5.03599977, 8.62639999, 0. , 0. , 0. ,0. , 11.11680031, 0. , 0. , 0. ,0. , 0. , 0. , 0. , 14.37839985,0. , 0. , 0. , 0. , 0. ,0. , 8.16479969, 0. , 7.30719948, 0. ,0. , 0. , 3.41039991, 7.17280006, 0. ,0. , 0. , 0. , 10.00992 , 0. ,0. , 0. , 0. , 0. , 0. ,0. , 13.97839928, 0. , 0. , 0. ,0. , 0. , 0. , 7.6855998 , 0. ,6.86559963, 0. , 0. , 0. , 3.21600008,7.93599987, 0. , 0. , 0. , 0. ,11.55999947, 0. , 0. , 0. , 0. ,0. , 0. , 0. , 18.76399994, 0. ,0. , 0. , 0. , 0. , 0. ,10.003304 , 0. , 8.10639954, 0. , 0. ,0. , 4.76480007, 6.87679958, 0. , 0. ,0. , 0. , 11.42239952, 0. , 0. ,0. , 0. , 0. , 0. , 0. ,19.42639732, 0. , 0. , 0. , 0. ,0. , 0. , 10.00524 , 0. , 8.2567997 ,0. , 0. , 0. , 5.08239985, 7.9776001 ,0. , 0. , 0. , 0. , 10.009984 ,0. , 0. , 0. , 0. , 0. ,0. , 0. , 11.5855999 , 0. , 0. ,0. , 0. , 0. , 0. , 7.88399982,0. , 5.96799994, 0. , 0. , 0. ,3.07679987, 7.81360006, 0. , 0. , 0. ,0. , 11.51119995, 0. , 0. , 0. ,0. , 0. , 0. , 0. , 20.00309599,0. , 0. , 0. , 0. , 0. ,0. , 10.005088 , 0. , 8.20479965, 0. ,<details><summary>英文:</summary>I'm not clear on why you don't trust `sci.stats`. If it's because of performance, have you actually profiled? It's worth noting that neither your fit nor scipy's fit are particularly good, because your data are strongly skewed toward 0.When initialising your second `x`, don't start it at 0 because your calculations are not well-defined for that point, which (a) dumps warnings on the console and (b) requires a bunch of indexing gymnastics later.I'm also not clear on why you're attempting `trapz`. Why integrate when you already have a CDF? Just use the CDF; and also don't `argmax`; use `searchsorted`.Finally, your plotting is slightly a disaster. For apples-to-apples comparison it's very important that you pass `cumulative` to `hist`.```pythonimport numpy as npimport scipy as sciimport matplotlib.pyplot as pltdata = np.array([0. , 0. , 11.26399994, 0. , 0. ,0. , 0. , 0. , 0. , 0. ,17.06399918, 0. , 0. , 0. , 0. ,0. , 0. , 8.33279991, 0. , 7.54879951,0. , 0. , 0. , 4.58799982, 7.9776001 ,0. , 0. , 0. , 0. , 11.45040035,0. , 0. , 0. , 0. , 0. ,0. , 0. , 18.73279953, 0. , 0. ,0. , 0. , 0. , 0. , 8.94559956,0. , 7.73040009, 0. , 0. , 0. ,5.03599977, 8.62639999, 0. , 0. , 0. ,0. , 11.11680031, 0. , 0. , 0. ,0. , 0. , 0. , 0. , 14.37839985,0. , 0. , 0. , 0. , 0. ,0. , 8.16479969, 0. , 7.30719948, 0. ,0. , 0. , 3.41039991, 7.17280006, 0. ,0. , 0. , 0. , 10.00992 , 0. ,0. , 0. , 0. , 0. , 0. ,0. , 13.97839928, 0. , 0. , 0. ,0. , 0. , 0. , 7.6855998 , 0. ,6.86559963, 0. , 0. , 0. , 3.21600008,7.93599987, 0. , 0. , 0. , 0. ,11.55999947, 0. , 0. , 0. , 0. ,0. , 0. , 0. , 18.76399994, 0. ,0. , 0. , 0. , 0. , 0. ,10.003304 , 0. , 8.10639954, 0. , 0. ,0. , 4.76480007, 6.87679958, 0. , 0. ,0. , 0. , 11.42239952, 0. , 0. ,0. , 0. , 0. , 0. , 0. ,19.42639732, 0. , 0. , 0. , 0. ,0. , 0. , 10.00524 , 0. , 8.2567997 ,0. , 0. , 0. , 5.08239985, 7.9776001 ,0. , 0. , 0. , 0. , 10.009984 ,0. , 0. , 0. , 0. , 0. ,0. , 0. , 11.5855999 , 0. , 0. ,0. , 0. , 0. , 0. , 7.88399982,0. , 5.96799994, 0. , 0. , 0. ,3.07679987, 7.81360006, 0. , 0. , 0. ,0. , 11.51119995, 0. , 0. , 0. ,0. , 0. , 0. , 0. , 20.00309599,0. , 0. , 0. , 0. , 0. ,0. , 10.005088 , 0. , 8.20479965, 0. ,0. , 0. , 5.51599979, 9.02879906, 0. ,0. ])def scigamma(data: np.ndarray) -> np.ndarray:param = sci.stats.gamma.fit(data)x = np.linspace(0, np.max(data), 250)cdf = sci.stats.gamma.cdf(x, *param)return cdfscicdf = scigamma(data)# Method of Moments estimationmean = np.mean(data)variance = np.var(data)alpha = mean**2 / variancebeta = variance / mean# Generate x-axis values for the curves# Do not start at 0x = np.linspace(0.1, np.max(data), 250)# Calculate the gamma distribution PDF valuespdf = x**(alpha - 1) * np.exp(-x / beta) / (beta**alpha * np.math.gamma(alpha))# Calculate the gamma distribution CDF valuescdf = np.cumsum(pdf) / np.sum(pdf)def plot() -> None:plt.hist(data, bins=50, density=True, cumulative=True, alpha=0.5, label='data')plt.plot(x, cdf, 'b', label='Gamma CDF | Custom Fit')plt.plot(x, scicdf, 'k', label='Gamma CDF | SciPy Fit')plt.xlabel('data')plt.ylabel('Probability')plt.legend()plt.show()def check_percentiles() -> None:# Calculate the percentile using numerical integrationvalue = 2.50index = np.searchsorted(x, value)pct = cdf[index]print(f'x={value} is at percentile {pct:.3f}')# Calculate the value using numerical integrationpct = 0.50index = np.searchsorted(cdf, pct)val = cdf[index]print(f'percentile {pct} is x={val:.3f}')check_percentiles()plot()
x=2.5 is at percentile 0.680percentile 0.5 is x=0.503
Much more directly, though,
param = sci.stats.gamma.fit(data)# Calculate the percentilevalue = 2.50pct = scipy.stats.gamma.cdf(value, *param)print(f'x={value} is at percentile {pct:.3f}')# Calculate the valuepct = 0.50val = scipy.stats.gamma.ppf(pct, *param)print(f'percentile {pct} is x={val:.3f}')
It's also worth noting that if you want scipy to produce a fit similar to yours, then pass method='MM' to fit():
As in the comments, a single gamma will not fit your data well. Consider instead a simple bimodal where the first CDF is a Heaviside step and the second CDF is still assumed gamma:
x = np.linspace(start=0, stop=20, num=250)# first mode of the bimodal distribution produces a heaviside step in the CDFnonzero_data = data[data.nonzero()]heaviside = 1 - nonzero_data.size/data.sizecdf = np.full_like(x, fill_value=heaviside)# second mode of the bimodal distribution assumed still gammaparams = scipy.stats.gamma.fit(nonzero_data)print('Gamma params:', params)cdf += scipy.stats.gamma.cdf(x, *params)*(1 - heaviside)plt.hist(data, bins=50, density=True, cumulative=True, alpha=0.5, label='data')plt.plot(x, cdf, label='fit')plt.xlabel('data')plt.ylabel('Probability')plt.legend(loc='right')plt.show()
Gamma params: (3.9359070026702017, 1.3587246910542263, 2.0440735483494272)
答案2
得分: 0
使用以下变量：
x = np.linspace(0, np.max(data), 250)cdf[x > 0] = np.cumsum(pdf[x > 0]) / np.sum(pdf[x > 0])
你可以使用 scipy 的 interpolate 来拟合函数，以获取从数值到百分位数和反之亦然的结果：
from scipy import interpolatef_x_to_percentile = interpolate.interp1d(x, cdf)x_val = 2.5print(f"Percentile of value {x_val}: {f_x_to_percentile(x_val)}")# Percentile of value 2.5: 0.6864894891059873f_percentile_to_x = interpolate.interp1d(cdf, x)percentile_val = 0.5print(f"Value at percentile {percentile_val*100}: {f_percentile_to_x(percentile_val)}")# Value at percentile 50.0: 1.015695412520532
更新:
如果你决定继续使用 scipy 来拟合你的分布，你可以直接使用 param 对象中的拟合参数，使用 .cdf() 和 .ppf() 方法获取结果：
from scipy import statsparam = stats.gamma.fit(data)x_val = 2.5percentile = stats.gamma(*param).cdf(x_val)print(f"Percentile of value {x_val}: {percentile}")# Percentile of value 2.5: 0.7681790728985926percentile_val = 0.5value_at_percentile = stats.gamma(*param).ppf(0.5)print(f"Value at percentile {percentile_val*100}: {value_at_percentile}")# Value at percentile 50.0: 0.7536517552749796
英文:
Using the variables below
x = np.linspace(0, np.max(data), 250)cdf[x > 0] = np.cumsum(pdf[x > 0]) / np.sum(pdf[x > 0])
you could fit a curve using scipy's interpolate to fit functions to obtain percentile from values and vice versa:
from scipy import interpolatef_x_to_percentile = interpolate.interp1d(x, cdf)x_val = 2.5print(f"Percentile of value {x_val}: {f_x_to_percentile(x_val)}")# Percentile of value 2.5: 0.6864894891059873f_percentile_to_x= interpolate.interp1d(cdf, x)percentile_val = 0.5print(f"Value at percentile {percentile_val*100}: {f_percentile_to_x(percentile_val)}")# Value at percentile 50.0: 1.015695412520532
Update:
If you decide to continue using scipy to fit your distribution, you can directly use the fitted parameters in your param object to obtain the results using the .cdf() and .ppf() methods:
from scipy import statsparam = stats.gamma.fit(data)x_val = 2.5percentile = stats.gamma(*param).cdf(x_val)print(f"Percentile of value {x_val}: {percentile}")# Percentile of value 2.5: 0.7681790728985926percentile_val = 0.5value_at_percentile = stats.gamma(*param).ppf(0.5)print(f"Value at percentile {percentile_val*100}: {value_at_percentile}")# Value at percentile 50.0: 0.7536517552749796
通过集体智慧和协作来改善编程学习和解决问题的方式。致力于成为全球开发者共同参与的知识库，让每个人都能够通过互相帮助和分享经验来进步。
评论