在Python中,如何计算给定分布示例列表的值的概率?

2024-03-29 09:23:33 发布

您现在位置:Python中文网/ 问答频道 /正文

不确定这是否属于统计数据,但我正在尝试使用Python来实现这一点。我基本上只有一个整数列表:

data = [300,244,543,1011,300,125,300 ... ]

我想知道给定这个数据值出现的概率。 我使用matplotlib绘制了数据的柱状图,并得到了这些:

enter image description here

enter image description here

在第一个图中,数字表示序列中的字符数。在第二个图表中,它是以毫秒为单位测量的时间量。最小值大于零,但不一定有最大值。这些图表是用上百万个例子创建的,但我不确定我是否可以对分布做任何其他假设。我想知道一个新值的概率,因为我有几百万个值的例子。在第一张图中,我有几百万个不同长度的序列。例如,想知道200长度的概率。

我知道对于一个连续的分布,任何精确点的概率都应该是零,但是给定一系列新的值,我需要能够说出每个值的可能性有多大。我已经浏览了一些numpy/scipy概率密度函数,但是我不确定在运行scipy.stats.norm.pdf(data)之类的函数之后,应该从哪个函数中选择,或者如何查询新的值。似乎不同的概率密度函数对数据的拟合不同。鉴于直方图的形状,我不知道如何决定使用哪一种。


Tags: 数据函数列表datamatplotlib图表绘制序列
3条回答

由于您似乎没有考虑特定的分布,但您可能有很多数据样本,我建议使用非参数密度估计方法。您描述的数据类型之一(时间单位为毫秒)显然是连续的,用于连续随机变量概率密度函数(PDF)的非参数估计的一种方法是您已经提到的直方图。但是,正如您将在下面看到的,Kernel Density Estimation (KDE)可能更好。您描述的第二类数据(序列中的字符数)是离散类型的。在这里,核密度估计也很有用,可以看作是一种平滑技术,适用于没有足够数量的样本来处理离散变量的所有值的情况。

估计密度

下面的示例演示如何首先从2个高斯分布的混合中生成数据样本,然后应用核密度估计来查找概率密度函数:

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.mlab as mlab
from sklearn.neighbors import KernelDensity

# Generate random samples from a mixture of 2 Gaussians
# with modes at 5 and 10
data = np.concatenate((5 + np.random.randn(10, 1),
                       10 + np.random.randn(30, 1)))

# Plot the true distribution
x = np.linspace(0, 16, 1000)[:, np.newaxis]
norm_vals = mlab.normpdf(x, 5, 1) * 0.25 + mlab.normpdf(x, 10, 1) * 0.75
plt.plot(x, norm_vals)

# Plot the data using a normalized histogram
plt.hist(data, 50, normed=True)

# Do kernel density estimation
kd = KernelDensity(kernel='gaussian', bandwidth=0.75).fit(data)

# Plot the estimated densty
kd_vals = np.exp(kd.score_samples(x))
plt.plot(x, kd_vals)

# Show the plots
plt.show()

这将生成以下绘图,其中真实分布显示为蓝色,直方图显示为绿色,使用KDE估计的PDF显示为红色:

Plot

如您所见,在这种情况下,由直方图近似的PDF不是很有用,而KDE提供了更好的估计。然而,对于大量的数据样本和适当的仓位大小选择,直方图也可能产生一个很好的估计。

在KDE情况下可以调整的参数是内核带宽。可以将内核看作是估计PDF的构建块,Scikit Learn中提供了几个内核函数:gaussian、tophat、epanechnikov、exponential、linear、cosine。改变带宽允许您调整偏差-方差权衡。较大的带宽将导致增加的偏差,这是好的,如果你有较少的数据样本。带宽越小,方差越大(估计中包含的样本越少),但如果有更多的样本可用,则估计效果越好。

计算概率

对于PDF,概率是通过计算一系列值上的积分得到的。正如您所注意到的,这将导致特定值的概率为0。

Scikit Learn似乎没有用于计算概率的内置函数。然而,很容易估计PDF在一个范围内的积分。我们可以在这个范围内对PDF进行多次求值,并将得到的值乘以每个求值点之间的步长。在下面的示例中,通过步骤step获得N样本。

# Get probability for range of values
start = 5  # Start of the range
end = 6    # End of the range
N = 100    # Number of evaluation points 
step = (end - start) / (N - 1)  # Step size
x = np.linspace(start, end, N)[:, np.newaxis]  # Generate values in the range
kd_vals = np.exp(kd.score_samples(x))  # Get PDF values for each x
probability = np.sum(kd_vals * step)  # Approximate the integral of the PDF
print(probability)

请注意,kd.score_samples生成数据样本的对数似然性。因此,需要np.exp来获得可能性。

同样的计算也可以使用内置的SciPy积分方法进行,这样会得到更精确的结果:

from scipy.integrate import quad
probability = quad(lambda x: np.exp(kd.score_samples(x)), start, end)[0]

例如,对于一次运行,第一个方法将概率计算为0.0859024655305,而第二个方法生成0.0850974209996139

好的,我提供这个作为起点,但是估计密度是一个非常广泛的话题。对于涉及序列中字符数量的情况,我们可以从一个直接的频率学家的角度,使用经验概率对其进行建模。在这里,概率本质上是百分率概念的推广。在我们的模型中,样本空间是离散的,并且都是正整数。好吧,然后你只需计算发生的次数,除以事件总数,就可以得到概率的估计值。在任何观测值为零的地方,我们对概率的估计为零。

>>> samples = [1,1,2,3,2,2,7,8,3,4,1,1,2,6,5,4,8,9,4,3]
>>> from collections import Counter
>>> counts = Counter(samples)
>>> counts
Counter({1: 4, 2: 4, 3: 3, 4: 3, 8: 2, 5: 1, 6: 1, 7: 1, 9: 1})
>>> total = sum(counts.values())
>>> total
20
>>> probability_mass = {k:v/total for k,v in counts.items()}
>>> probability_mass
{1: 0.2, 2: 0.2, 3: 0.15, 4: 0.15, 5: 0.05, 6: 0.05, 7: 0.05, 8: 0.1, 9: 0.05}
>>> probability_mass.get(2,0)
0.2
>>> probability_mass.get(12,0)
0

现在,对于计时数据,将其建模为连续分布更为自然。不要使用参数化方法假设数据具有某种分布,然后将该分布与数据匹配,而应采用非参数化方法。一个简单的方法是使用kernel density estimate。您可以简单地将其视为平滑直方图的一种方法,以提供连续的概率密度函数。有几个可用的库。对于单变量数据,最直接的方法可能是scipy:

>>> import scipy.stats
>>> kde = scipy.stats.gaussian_kde(samples)
>>> kde.pdf(2)
array([ 0.15086911])

要获取某个时间间隔内的观测概率:

>>> kde.integrate_box_1d(1,2)
0.13855869478828692

这是一个可能的解决方案。计算原始列表中每个值的出现次数。给定值的未来概率是其过去发生率,即过去发生率除以原始列表的长度。在Python中非常简单:

x是给定的值列表

from collections import Counter
c = Counter(x)

def probability(a):
    # returns the probability of a given number a
    return float(c[a]) / len(x)

相关问题 更多 >