我如何操纵样本的 CDF 使其与不同样本的 CDF 相匹配?
How can I manipulate the CDF of a sample such that it matches that of a different sample?
我想使用 CDF 匹配来校正降水的原始模型预测(但该应用程序相当通用)。
假设下面的 CDF B 是观察到的 CDF(我信任的 CDF),我想计算 CDF A 和 B 之间的差异,以便在给定的一天我可以进行降水预报并将其移动到差异介于A和B之间,这样它更能代表B而不是A。
所以对于每个 x 值,我需要获取 A 的 y 值,然后 B 是我需要获取 x 值的相同值,给我 2 个 x 值来计算 a差异。
当然,这只会给我离散的 x 值,其中我知道校正,所以我想我需要做额外的工作来校正介于其他 2 个之间的 x 值。
这是我用来生成示例的 Python 代码:
import numpy.random
import numpy as np
from scipy.interpolate import interp1d
import matplotlib.pyplot as plt
quantiles = [0, 1, 2, 3, 4, 5, 7, 10, 15, 20, 30, 40, 50, 60, 75, 100]
# Generate fake precip data
sample_size = 100000
A = numpy.random.gamma(0.7, scale=50, size=sample_size)
B = numpy.random.gamma(0.5, scale=70, size=sample_size)
ens = (40 - 20) * np.random.random_sample((21)) + 20
# Calculate histograms
A_pdf, edges = np.histogram(A, bins=quantiles)
A_pdf = A_pdf / sample_size
A_cdf = np.cumsum(A_pdf)
B_pdf, edges = np.histogram(B, bins=quantiles)
B_pdf = B_pdf / sample_size
B_cdf = np.cumsum(B_pdf)
# Plot CDFs
plt.figure()
plt.plot(quantiles[1:], A_cdf, 'x-', c='r', lw=3, ms=10, mew=2, label='A')
plt.plot(quantiles[1:], B_cdf, '+-', c='k', lw=3, ms=15, mew=2, label='B')
plt.xticks(quantiles[1:])
plt.legend(loc='upper left')
谢谢大家!
您只需要一个近似 A 的 CDF 的函数,以及一个近似 B 的逆 CDF(或 PPF)的函数。然后您只需计算 qcorrected = PPFB(CDFA(q)).
对于您的示例数据,我们可以简单地将 .cdf
和 .ppf
方法与适当的参数一起用于 scipy.stats.gamma
frozen distributions:
from scipy import stats
distA = stats.gamma(0.7, scale=50)
distB = stats.gamma(0.5, scale=70)
corrected_quantiles = distB.ppf(distA.cdf(quantiles[1:]))
当然,对于真实数据,您不太可能知道真实基础分布的参数。如果您对它们的函数形式有很好的了解,您可以尝试对数据执行最大似然拟合以估计它们:
distA = stats.gamma(*stats.gamma.fit(A))
distB = stats.gamma(*stats.gamma.fit(B))
否则,您可以尝试从您的经验 CDF 中 interpolate/extrapolate,例如使用 scipy.interpolate.InterpolatedUnivariateSpline
:
from scipy.interpolate import InterpolatedUnivariateSpline
# cubic spline interpolation
itp_A_cdf = InterpolatedUnivariateSpline(quantiles[1:], A_cdf, k=3)
# the PPF is the inverse of the CDF, so we simply reverse the order of the
# x & y arguments to InterpolatedUnivariateSpline
itp_B_ppf = InterpolatedUnivariateSpline(B_cdf, quantiles[1:], k=3)
itp_corrected_quantiles = itp_B_ppf(itp_A_cdf(quantiles[1:]))
fig, ax = plt.subplots(1, 1)
ax.hold(True)
ax.plot(quantiles[1:], A_cdf, '-r', lw=3, label='A')
ax.plot(quantiles[1:], B_cdf, '-k', lw=3, label='B')
ax.plot(corrected_quantiles, A_cdf, '--xr', lw=3, ms=10, mew=2, label='exact')
ax.plot(itp_corrected_quantiles, A_cdf, '--+b', lw=3, ms=10, mew=2,
label='interpolated')
ax.legend(loc=5)
我想使用 CDF 匹配来校正降水的原始模型预测(但该应用程序相当通用)。
假设下面的 CDF B 是观察到的 CDF(我信任的 CDF),我想计算 CDF A 和 B 之间的差异,以便在给定的一天我可以进行降水预报并将其移动到差异介于A和B之间,这样它更能代表B而不是A。
所以对于每个 x 值,我需要获取 A 的 y 值,然后 B 是我需要获取 x 值的相同值,给我 2 个 x 值来计算 a差异。
当然,这只会给我离散的 x 值,其中我知道校正,所以我想我需要做额外的工作来校正介于其他 2 个之间的 x 值。
这是我用来生成示例的 Python 代码:
import numpy.random
import numpy as np
from scipy.interpolate import interp1d
import matplotlib.pyplot as plt
quantiles = [0, 1, 2, 3, 4, 5, 7, 10, 15, 20, 30, 40, 50, 60, 75, 100]
# Generate fake precip data
sample_size = 100000
A = numpy.random.gamma(0.7, scale=50, size=sample_size)
B = numpy.random.gamma(0.5, scale=70, size=sample_size)
ens = (40 - 20) * np.random.random_sample((21)) + 20
# Calculate histograms
A_pdf, edges = np.histogram(A, bins=quantiles)
A_pdf = A_pdf / sample_size
A_cdf = np.cumsum(A_pdf)
B_pdf, edges = np.histogram(B, bins=quantiles)
B_pdf = B_pdf / sample_size
B_cdf = np.cumsum(B_pdf)
# Plot CDFs
plt.figure()
plt.plot(quantiles[1:], A_cdf, 'x-', c='r', lw=3, ms=10, mew=2, label='A')
plt.plot(quantiles[1:], B_cdf, '+-', c='k', lw=3, ms=15, mew=2, label='B')
plt.xticks(quantiles[1:])
plt.legend(loc='upper left')
谢谢大家!
您只需要一个近似 A 的 CDF 的函数,以及一个近似 B 的逆 CDF(或 PPF)的函数。然后您只需计算 qcorrected = PPFB(CDFA(q)).
对于您的示例数据,我们可以简单地将 .cdf
和 .ppf
方法与适当的参数一起用于 scipy.stats.gamma
frozen distributions:
from scipy import stats
distA = stats.gamma(0.7, scale=50)
distB = stats.gamma(0.5, scale=70)
corrected_quantiles = distB.ppf(distA.cdf(quantiles[1:]))
当然,对于真实数据,您不太可能知道真实基础分布的参数。如果您对它们的函数形式有很好的了解,您可以尝试对数据执行最大似然拟合以估计它们:
distA = stats.gamma(*stats.gamma.fit(A))
distB = stats.gamma(*stats.gamma.fit(B))
否则,您可以尝试从您的经验 CDF 中 interpolate/extrapolate,例如使用 scipy.interpolate.InterpolatedUnivariateSpline
:
from scipy.interpolate import InterpolatedUnivariateSpline
# cubic spline interpolation
itp_A_cdf = InterpolatedUnivariateSpline(quantiles[1:], A_cdf, k=3)
# the PPF is the inverse of the CDF, so we simply reverse the order of the
# x & y arguments to InterpolatedUnivariateSpline
itp_B_ppf = InterpolatedUnivariateSpline(B_cdf, quantiles[1:], k=3)
itp_corrected_quantiles = itp_B_ppf(itp_A_cdf(quantiles[1:]))
fig, ax = plt.subplots(1, 1)
ax.hold(True)
ax.plot(quantiles[1:], A_cdf, '-r', lw=3, label='A')
ax.plot(quantiles[1:], B_cdf, '-k', lw=3, label='B')
ax.plot(corrected_quantiles, A_cdf, '--xr', lw=3, ms=10, mew=2, label='exact')
ax.plot(itp_corrected_quantiles, A_cdf, '--+b', lw=3, ms=10, mew=2,
label='interpolated')
ax.legend(loc=5)