python - 如何从经验 cdf 计算和绘制 pdf?
问题描述
我有两个 numpy 数组,一个是 x 值数组,另一个是 y 值数组,它们一起给了我经验 cdf。例如:
plt.plot(xvalues, yvalues)
plt.show()
我认为需要以某种方式平滑数据才能提供平滑的 pdf。
我想绘制pdf。我怎样才能做到这一点?
原始数据位于:http ://dpaste.com/1HVK5DR 。
解决方案
有两个主要问题:您的数据似乎很嘈杂,并且间隔不均:低端的点采样非常密集,而高端的点采样非常稀疏。这可能会导致数值问题。
所以首先我建议使用线性插值对数据进行重新采样以获得等间距的样本:(请注意,所有附加到彼此的片段形成一个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)
位置。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')
推荐阅读
- python - 在另一个 python 文件中导入变量不起作用(可能循环导入错误)
- c# - catch 后继续 try 块。在一次尝试中记录所有异常
- google-maps - PrimeFaces gmap (JSF 2.2) 中未显示新标记
- kotlin - IntelliJ 在新的 Ktor 项目中缺少对 slf4j 的引用
- python-3.x - 如何在python中将薪水分类为高/中/低组?
- algorithm - 数据结构与算法最优解说明
- android - 如何在 ViewModel Android 中获取活动上下文?
- django - 使用rest框架的Django自定义url
- php - 在php动态分页中更改记录限制时的问题
- rust - Rust/rocket 将变量传递给端点