如何根据经验 cdf 计算并绘制 pdf？

Question

我有两个 numpy 数组，一个是 x 值数组，另一个是 y 值数组，它们一起给出了经验 cdf。例如：

plt.plot(xvalues, yvalues)
plt.show()

我认为数据需要以某种方式平滑才能给出平滑的 pdf。

我想绘制pdf。我怎样才能做到这一点？

原始数据位于： http: //dpaste.com/1HVK5DR。

Answer 1

有两个主要问题：您的数据似乎非常嘈杂，并且间隔不均匀：低端的点采样非常密集，而高端的点采样非常稀疏。这可能会导致数字问题。

因此，首先我建议使用线性插值对数据进行重新采样，以获得等间距的样本：（请注意，相互附加的所有片段都形成一个python 文件的内容。）

import matplotlib.pyplot as plt
import numpy as np
from data import xvalues, yvalues #load data from file


print("#datapoints: {}".format(len(xvalues)))
#don't use every point if your computer is not very fast
xv = np.array(xvalues)[::5]
yv = np.array(yvalues)[::5]

#interpolate to have evenly space data
xi = np.linspace(xv.min(), xv.max(), 400)
yi = np.interp(xi, xv, yv)

然后，为了平滑数据，我建议执行 RBF 回归（=使用“RBF 网络”）。这个想法是拟合形式的曲线

 c(t) = sum a(i) * phi(t - x(i))    #(not part of the program)

其中phi是一些径向基函数。（理论上我们可以使用任何函数。）为了获得非常平滑的结果，我选择了一个非常平滑的函数，即高斯函数：phi(x) = exp( - x^2/sigma^2)其中sigma尚未确定。这些x(i)只是我们可以定义的一些节点。如果我们有一个平滑的函数，我们只需要几个节点。节点的数量还决定了需要完成多少计算。这些a(i)是我们可以优化以获得最佳拟合的系数。在这种情况下，我只使用最小二乘法。

请注意，如果我们可以按照上面的形式编写一个函数，那么计算导数就非常容易了，只需

 c(t) = sum a(i) * phi'(t - x(i))

其中phi'是 的导数phi。#（不是程序的一部分）

关于sigma：通常最好将其选择为我们选择的节点之间步长的倍数。我们选择的越大sigma，得到的函数就越平滑。

#set up rbf network
rbf_nodes = xv[::50][None, :]#use a subset of the x-values as rbf nodes
print("#rbfs: {}".format(rbf_nodes.shape[1]))
#estimate width of kernels:
sigma = 20  #greater = smoother, this is the primary parameter to play with
sigma *= np.max(np.abs(rbf_nodes[0,1:]-rbf_nodes[0,:-1]))


# kernel & derivative
rbf = lambda r:1/(1+(r/sigma)**2)
Drbf = lambda r: -2*r*sigma**2/(sigma**2 + r**2)**2

#compute coefficients of rbf network
r = np.abs(xi[:, None]-rbf_nodes)
A = rbf(r)
coeffs = np.linalg.lstsq(A, yi, rcond=None)[0]
print(coeffs)

#evaluate rbf network
N=1000
xe = np.linspace(xi.min(), xi.max(), N)
Ae = rbf(xe[:, None] - rbf_nodes)
ye = Ae @ coeffs


#evaluate derivative
N=1000
xd = np.linspace(xi.min(), xi.max(), N)
Bd = Drbf(xe[:, None] - rbf_nodes)
yd = Bd @ coeffs


fig,ax = plt.subplots()
ax2 = ax.twinx()

ax.plot(xv, yv, '-')
ax.plot(xi, yi, '-')
ax.plot(xe, ye, ':')

ax2.plot(xd, yd, '-')

fig.savefig('graph.png')
print('done')