首页 > 解决方案 > 如何从经验 cdf 计算和绘制 pdf?

问题描述

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

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

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

在此处输入图像描述

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

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

标签: pythonnumpymatplotlibstatisticsprobability

解决方案


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

所以首先我建议使用线性插值对数据进行重新采样以获得等间距的样本:(请注意,所有附加到彼此的片段形成一个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 Network")。这个想法是拟合形式的曲线

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

其中phi是一些径向基函数。(理论上我们可以使用任何函数。)为了得到一个非常平滑的结果,我选择了一个非常平滑的函数,即高斯函数:尚待确定的phi(x) = exp( - x^2/sigma^2)位置。sigmax(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')

在此处输入图像描述


推荐阅读