从密度分布中抽样随机值
Sample Random values from a density distribution
大家好，
我正在尝试从 0 到 1 之间的随机值进行采样，权重由上面的数据提供。我使用 scipy.stats.gaussian_kde 及其 .resample(n) 方法找到了该问题的部分解决方案。我的主要问题是，因为我的大部分数据非常接近于 0，所以重新采样 returns 一堆负数会搞乱我以后的计算。
有没有办法限制我的重采样全部大于零，而不改变样本space？我考虑过只取所有东西的绝对值去除底片，但我不知道这是否反映了分布权重。
澄清一下，我重新采样的每个值 (n) 都对应于我代码中的一个特定变量，所以我不能只删除小于零的数字。
# Here is a little sample dataset if you need something to work this out!
import numpy as np
data = np.array([0.147, 0.066, 0.017, 0.011, 0.040, 0.087, 0.024, 0.127, 0.071, 0.127,
0.027, 0.008, 0.067, 0.032, 0.247, 0.028, 0.122, 0.304, 0.074, 0.119])
# Thank you!
您可以使用支持度不包含负数的分布。例如，从指数分布中抽样可能适用于您提供的示例数组：
import numpy as np
from scipy.stats import expon
import matplotlib.pyplot as plt
data = np.array([0.147, 0.066, 0.017, 0.011, 0.040, 0.087, 0.024, 0.127, 0.071, 0.127, 0.027, 0.008, 0.067, 0.032, 0.247, 0.028, 0.122, 0.304, 0.074, 0.119])
# fit exponential model using data
loc, scale = expon.fit(data)
# plot histogram and model
fig, ax = plt.subplots()
ax.hist(data, density = True)
x = np.linspace(0.01, 1, 200)
ax.plot(x, expon.pdf(x, loc, scale), 'k-')
plt.show()
# sample from your modelled distribution using your fitted loc and scale parameters
sample = expon.rvs(loc, scale)
您可以使用拒绝抽样或逆变换抽样直接对数据进行抽样。
该算法有 4 个步骤：
- 计算数据的累积分布函数 CDF。
- 从区间 [0,1] 的均匀分布中抽取一个随机数
u。
- 求 CDF 的倒数
InverseCDF。
- 计算
X = InverseCDF(u) 其中 X 将分配为 CDF(X)。
下面是一个 Python 示例代码，使用 NumPy 来说明算法：
import numpy as np
import numpy.random as ra
from matplotlib import pyplot as plt
'''
Step 1: Compute cumulative distribution function
'''
data = [[1, 2, 3, 4, 5, 6],[100, 200, 400, 310, 130, 50]]
event_value = np.array(data[0])
event_frequency = np.array(data[1])
prob = event_frequency/float(sum(event_frequency))
cum_prob = np.cumsum(prob)
print(prob)
print(cum_prob)
'''
Step 2: Sample uniform random numbers
'''
N = 10000
U = ra.uniform(0, 1, N)
'''
Step 3: Generate Samples from Inverse CDF:
'''
sample_x = [int(event_value[np.argwhere(cum_prob == min(cum_prob[(cum_prob - u) > 0]))]) for u in U]
'''
Step 4: Comparison
'''
sample_x = (np.array(sample_x)-1).astype(int)
times = np.arange(0,6,1)
lc = np.bincount(sample_x, minlength=len(times))
plot1, = plt.plot(lc/float(sum(lc)), 'r--', label='Sampled events')
plot2, = plt.plot(prob,'g',label='Original data')
plt.xlabel('Event Value')
plt.ylabel('Probability')
plt.legend(handles=[plot1,plot2])
plt.show()
输出：
要完成 Ben Devries 的回答，您有多种选择来处理这种情况。您面临的是一个 cut-off 为零的分布（让我们指出这可能不是这种情况，但了解您的数据，您似乎确定这是不可能的，这没关系！）。 Gaussian KDE 处理得不是很好，因为它们通常是在真实 space.
上定义的
明智的选择是从 KDE 切换到参数估计。这意味着您假设概率密度的一种形式（基于数据的形式和您对它们来源的了解）并尝试调整概率密度参数（例如许多分布的位置和比例）以使密度适合数据.在您的情况下，分布看起来很像指数分布。
如果你没有头绪，你可以尝试坚持使用 KDE（顺便说一句 non-parametric 估计的一种形式），并使用随机变量转换来管理 cut-off：你尝试对你的数据应用一个函数，这样它就可以很容易地被 KDE 拟合，并在之后很容易地返回到原始分布。
在那里，对数转换似乎是应用的完美函数。 A link 用于详细信息：https://thirdorderscientist.org/homoclinic-orbit/2013/10/24/kernel-density-estimation-for-random-variables-with-bounded-support-mdash-the-transformation-trick
在下一个脚本中，我尝试绘制提到的各种方法（直接 KDE，如 Ben DeVries 所说，用指数拟合），以及转换变量上的 KDE（我也为此重新采样变量，并绘制直方图):
import numpy as np
from scipy.stats import expon, gaussian_kde
import matplotlib.pyplot as plt
data = np.array([0.147, 0.066, 0.017, 0.011, 0.040, 0.087, 0.024, 0.127, 0.071, 0.127,
0.027, 0.008, 0.067, 0.032, 0.247, 0.028, 0.122, 0.304, 0.074, 0.119])
# gaussian KDE
kde = gaussian_kde(data)
x_kde = np.linspace(-0.2, 0.6, 10000)
y_kde = kde.evaluate(x_k)
# fit an exponential distribution
loc, scale = expon.fit(data)
x_exp = np.linspace(0.01, 0.6, 10000)
y_exp = expon.pdf(x_exp, loc, scale)
# use random variable transformation
log_data = np.log(data)
kde_log = gaussian_kde(log_data)
x_kl = np.linspace(-5, -0.5, 10000) # don't hesitate to plot kde_log over x_kl !
x_klt = np.exp(x_kl) # turn back to x (y = log(x) <-> x = exp(y))
y_klt = [kde_log.evaluate(x) / np.exp(x) for x in x_kl] # turn back to fX(x) = fY(log(x)) / x
fig, ax = plt.subplots(1, 1, dpi=120)
ax.hist(data, label="data", density=True)
ax.plot(x_kde, y_kde, label="kde")
ax.plot(x_exp, y_exp, label="exp")
ax.plot(x_klt, y_klt, label="kde_transformed")
#let's try to resample using KDE.resample(n)
x_resample = kde_log.resample(1000)[0]
y = [np.exp(x) for x in x_resample]
ax.hist(y, density=True, bins=100, histtype='step', #plot using 'step' to see something
label="resampled_from_transformed")
ax.legend()
剧情：
如您所见：
- 指数与转换后的 KDE 非常匹配。这似乎是一个安全的选择。
- 转换后的 KDE 没有像 0 附近的指数拟合那样表现出尖锐的 cut-off；这是我认为的典型技术（https://stats.stackexchange.com/questions/403532/can-kernel-density-estimation-estimate-an-exponential-distribution 了解更多详情）
大家好，
我正在尝试从 0 到 1 之间的随机值进行采样，权重由上面的数据提供。我使用 scipy.stats.gaussian_kde 及其 .resample(n) 方法找到了该问题的部分解决方案。我的主要问题是，因为我的大部分数据非常接近于 0，所以重新采样 returns 一堆负数会搞乱我以后的计算。
有没有办法限制我的重采样全部大于零，而不改变样本space？我考虑过只取所有东西的绝对值去除底片，但我不知道这是否反映了分布权重。
澄清一下，我重新采样的每个值 (n) 都对应于我代码中的一个特定变量，所以我不能只删除小于零的数字。
# Here is a little sample dataset if you need something to work this out!
import numpy as np
data = np.array([0.147, 0.066, 0.017, 0.011, 0.040, 0.087, 0.024, 0.127, 0.071, 0.127,
0.027, 0.008, 0.067, 0.032, 0.247, 0.028, 0.122, 0.304, 0.074, 0.119])
# Thank you!
您可以使用支持度不包含负数的分布。例如，从指数分布中抽样可能适用于您提供的示例数组：
import numpy as np
from scipy.stats import expon
import matplotlib.pyplot as plt
data = np.array([0.147, 0.066, 0.017, 0.011, 0.040, 0.087, 0.024, 0.127, 0.071, 0.127, 0.027, 0.008, 0.067, 0.032, 0.247, 0.028, 0.122, 0.304, 0.074, 0.119])
# fit exponential model using data
loc, scale = expon.fit(data)
# plot histogram and model
fig, ax = plt.subplots()
ax.hist(data, density = True)
x = np.linspace(0.01, 1, 200)
ax.plot(x, expon.pdf(x, loc, scale), 'k-')
plt.show()
# sample from your modelled distribution using your fitted loc and scale parameters
sample = expon.rvs(loc, scale)
您可以使用拒绝抽样或逆变换抽样直接对数据进行抽样。 该算法有 4 个步骤：
- 计算数据的累积分布函数 CDF。
- 从区间 [0,1] 的均匀分布中抽取一个随机数
u。 - 求 CDF 的倒数
InverseCDF。 - 计算
X = InverseCDF(u)其中 X 将分配为CDF(X)。
下面是一个 Python 示例代码，使用 NumPy 来说明算法：
import numpy as np
import numpy.random as ra
from matplotlib import pyplot as plt
'''
Step 1: Compute cumulative distribution function
'''
data = [[1, 2, 3, 4, 5, 6],[100, 200, 400, 310, 130, 50]]
event_value = np.array(data[0])
event_frequency = np.array(data[1])
prob = event_frequency/float(sum(event_frequency))
cum_prob = np.cumsum(prob)
print(prob)
print(cum_prob)
'''
Step 2: Sample uniform random numbers
'''
N = 10000
U = ra.uniform(0, 1, N)
'''
Step 3: Generate Samples from Inverse CDF:
'''
sample_x = [int(event_value[np.argwhere(cum_prob == min(cum_prob[(cum_prob - u) > 0]))]) for u in U]
'''
Step 4: Comparison
'''
sample_x = (np.array(sample_x)-1).astype(int)
times = np.arange(0,6,1)
lc = np.bincount(sample_x, minlength=len(times))
plot1, = plt.plot(lc/float(sum(lc)), 'r--', label='Sampled events')
plot2, = plt.plot(prob,'g',label='Original data')
plt.xlabel('Event Value')
plt.ylabel('Probability')
plt.legend(handles=[plot1,plot2])
plt.show()
输出：
要完成 Ben Devries 的回答，您有多种选择来处理这种情况。您面临的是一个 cut-off 为零的分布（让我们指出这可能不是这种情况，但了解您的数据，您似乎确定这是不可能的，这没关系！）。 Gaussian KDE 处理得不是很好，因为它们通常是在真实 space.
上定义的明智的选择是从 KDE 切换到参数估计。这意味着您假设概率密度的一种形式（基于数据的形式和您对它们来源的了解）并尝试调整概率密度参数（例如许多分布的位置和比例）以使密度适合数据.在您的情况下，分布看起来很像指数分布。
如果你没有头绪，你可以尝试坚持使用 KDE（顺便说一句 non-parametric 估计的一种形式），并使用随机变量转换来管理 cut-off：你尝试对你的数据应用一个函数，这样它就可以很容易地被 KDE 拟合，并在之后很容易地返回到原始分布。
在那里，对数转换似乎是应用的完美函数。 A link 用于详细信息：https://thirdorderscientist.org/homoclinic-orbit/2013/10/24/kernel-density-estimation-for-random-variables-with-bounded-support-mdash-the-transformation-trick
在下一个脚本中，我尝试绘制提到的各种方法（直接 KDE，如 Ben DeVries 所说，用指数拟合），以及转换变量上的 KDE（我也为此重新采样变量，并绘制直方图):
import numpy as np
from scipy.stats import expon, gaussian_kde
import matplotlib.pyplot as plt
data = np.array([0.147, 0.066, 0.017, 0.011, 0.040, 0.087, 0.024, 0.127, 0.071, 0.127,
0.027, 0.008, 0.067, 0.032, 0.247, 0.028, 0.122, 0.304, 0.074, 0.119])
# gaussian KDE
kde = gaussian_kde(data)
x_kde = np.linspace(-0.2, 0.6, 10000)
y_kde = kde.evaluate(x_k)
# fit an exponential distribution
loc, scale = expon.fit(data)
x_exp = np.linspace(0.01, 0.6, 10000)
y_exp = expon.pdf(x_exp, loc, scale)
# use random variable transformation
log_data = np.log(data)
kde_log = gaussian_kde(log_data)
x_kl = np.linspace(-5, -0.5, 10000) # don't hesitate to plot kde_log over x_kl !
x_klt = np.exp(x_kl) # turn back to x (y = log(x) <-> x = exp(y))
y_klt = [kde_log.evaluate(x) / np.exp(x) for x in x_kl] # turn back to fX(x) = fY(log(x)) / x
fig, ax = plt.subplots(1, 1, dpi=120)
ax.hist(data, label="data", density=True)
ax.plot(x_kde, y_kde, label="kde")
ax.plot(x_exp, y_exp, label="exp")
ax.plot(x_klt, y_klt, label="kde_transformed")
#let's try to resample using KDE.resample(n)
x_resample = kde_log.resample(1000)[0]
y = [np.exp(x) for x in x_resample]
ax.hist(y, density=True, bins=100, histtype='step', #plot using 'step' to see something
label="resampled_from_transformed")
ax.legend()
剧情：
如您所见：
- 指数与转换后的 KDE 非常匹配。这似乎是一个安全的选择。
- 转换后的 KDE 没有像 0 附近的指数拟合那样表现出尖锐的 cut-off；这是我认为的典型技术（https://stats.stackexchange.com/questions/403532/can-kernel-density-estimation-estimate-an-exponential-distribution 了解更多详情）