英文:
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, math
import numpy as np
import scipy as sci
import matplotlib.pyplot as plt
data = np.array([0.00, 0.00, 11.26399994, 0.00, 0.00, 0.00, 0.00, 0.00, ... ]) # 这里是你的数据
def scigamma(data):
# 这是你的函数实现
return cdf
scicdf = 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, math
import numpy as np
import scipy as sci
import matplotlib.pyplot as plt
data = 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 cdf
scicdf = scigamma(data)
# Method of Moments estimation
mean = np.mean(data)
variance = np.var(data)
alpha = mean**2 / variance
beta = variance / mean
# Generate x-axis values for the curves
x = np.linspace(0, np.max(data), 250)
# Calculate the gamma distribution PDF values
pdf = (x ** (alpha - 1) * np.exp(-x / beta)) / (beta ** alpha * np.math.gamma(alpha))
# Calculate the gamma distribution CDF values
cdf = np.zeros_like(x)
cdf[x > 0] = np.cumsum(pdf[x > 0]) / np.sum(pdf[x > 0])
# Estimate the probability of zero values
num_zeros = np.count_nonzero(data == 0)
zero_probability = np.count_nonzero(data == 0) / len(data)
# Calculate the PDF and CDF values at zero
pdf_zero = zero_probability / (beta ** alpha * np.math.gamma(alpha))
cdf_zero = zero_probability
value = 2.50
percentile = 0.50
index = np.argmax(pdf >= value)
# Calculate the percentile using numerical integration
pct = 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 integration
val = np.trapz(cdf[:index+1], dx=1) + (percentile - cdf[index-1]) * (pdf[index] - pdf[index-1]) / (cdf[index] - cdf[index-1])
# Plot the data histogram
plt.hist(data, bins=30, density=True, alpha=0.5, label='data')
# Plot the gamma distribution CDF curve
plt.plot(x, cdf, 'b', label='Gamma CDF | Custom Fit')
plt.plot(x, scicdf, 'k', label='Gamma CDF | SciPy Fit')
# Set plot labels and legend
plt.xlabel('data')
plt.ylabel('Probability')
plt.legend()
[1]: https://stackoverflow.com/q/75598344/10266106
[2]: https://i.stack.imgur.com/QchWe.png
</details>
# 答案1
**得分**: 1
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 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`.
```python
import numpy as np
import scipy as sci
import matplotlib.pyplot as plt
data = 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`.
```python
import numpy as np
import scipy as sci
import matplotlib.pyplot as plt
data = 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 cdf
scicdf = scigamma(data)
# Method of Moments estimation
mean = np.mean(data)
variance = np.var(data)
alpha = mean**2 / variance
beta = variance / mean
# Generate x-axis values for the curves
# Do not start at 0
x = np.linspace(0.1, np.max(data), 250)
# Calculate the gamma distribution PDF values
pdf = x**(alpha - 1) * np.exp(-x / beta) / (beta**alpha * np.math.gamma(alpha))
# Calculate the gamma distribution CDF values
cdf = 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 integration
value = 2.50
index = np.searchsorted(x, value)
pct = cdf[index]
print(f'x={value} is at percentile {pct:.3f}')
# Calculate the value using numerical integration
pct = 0.50
index = 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.680
percentile 0.5 is x=0.503
Much more directly, though,
param = sci.stats.gamma.fit(data)
# Calculate the percentile
value = 2.50
pct = scipy.stats.gamma.cdf(value, *param)
print(f'x={value} is at percentile {pct:.3f}')
# Calculate the value
pct = 0.50
val = 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 CDF
nonzero_data = data[data.nonzero()]
heaviside = 1 - nonzero_data.size/data.size
cdf = np.full_like(x, fill_value=heaviside)
# second mode of the bimodal distribution assumed still gamma
params = 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 interpolate
f_x_to_percentile = interpolate.interp1d(x, cdf)
x_val = 2.5
print(f"Percentile of value {x_val}: {f_x_to_percentile(x_val)}")
# Percentile of value 2.5: 0.6864894891059873
f_percentile_to_x = interpolate.interp1d(cdf, x)
percentile_val = 0.5
print(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 stats
param = stats.gamma.fit(data)
x_val = 2.5
percentile = stats.gamma(*param).cdf(x_val)
print(f"Percentile of value {x_val}: {percentile}")
# Percentile of value 2.5: 0.7681790728985926
percentile_val = 0.5
value_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 interpolate
f_x_to_percentile = interpolate.interp1d(x, cdf)
x_val = 2.5
print(f"Percentile of value {x_val}: {f_x_to_percentile(x_val)}")
# Percentile of value 2.5: 0.6864894891059873
f_percentile_to_x= interpolate.interp1d(cdf, x)
percentile_val = 0.5
print(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 stats
param = stats.gamma.fit(data)
x_val = 2.5
percentile = stats.gamma(*param).cdf(x_val)
print(f"Percentile of value {x_val}: {percentile}")
# Percentile of value 2.5: 0.7681790728985926
percentile_val = 0.5
value_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
通过集体智慧和协作来改善编程学习和解决问题的方式。致力于成为全球开发者共同参与的知识库，让每个人都能够通过互相帮助和分享经验来进步。
评论