平滑数据并找到最大值

Question

我有一个数据集，里面有两个变量，x和y。我想找出在什么样的x值下，y的值达到最大。现在我做法是直接找出那个能让我得到最大y的x值。但这样做不太好，因为我的数据有点杂乱，所以我想先对数据进行平滑处理，然后再找最大值。

到目前为止，我尝试用R语言的npreg（核回归）这个功能来平滑我的数据，使用的是np这个包，得到了下面这个曲线：

但我不太确定怎么找到最大值。

我希望能在Python中解决以下问题：

1) 对数据进行平滑处理（不一定要用核回归）

2) 使用平滑后的数据找到y值最大的x值

x   y
-20 0.006561733
-19 -4.48E-08
-18 -4.48E-08
-17 -4.48E-08
-16 0.003281305
-15 0.00164063
-14 0.003280565
-13 0.003282537
-12 -4.48E-08
-11 0.003281286
-10 0.004921239
-9  0.00491897
-8  -1.52E-06
-7  0.004925867
-6  -1.27E-06
-5  0.009839438
-4  0.001643726
-3  -4.48E-08
-2  2.09E-06
-1  -0.001640027
0   0.006559627
1   0.001636958
2   2.36E-06
3   0.003281469
4   0.011481469
5   0.004922279
6   0.018044207
7   0.011483134
8   0.014765087
9   0.008201379
10  0.00492497
11  0.006560482
12  0.009844796
13  0.011483199
14  0.008202129
15  0.001641621
16  0.004921645
17  0.006563377
18  0.006561068
19  0.008201004

Answer 1

我不太确定你主要想解决什么问题？是想要更好的平滑效果，还是找到最小值，或者是想全部用Python来做？如果你在R语言中已经取得了不错的进展，为什么还要转到Python呢？我发现R语言里的内置函数supsmu通常能很好地进行非参数平滑处理。这是我在R中处理这个问题的方法。

smooth <- do.call(supsmu, data)
min.idx <- which.min(smooth$y)
min.point <- c(smooth$x[min.idx], smooth$y[min.idx])

Answer 2

你可以试试使用平滑样条函数：

import numpy as np
from scipy import interpolate
x = range(-20,20)
y = [0.006561733, -4.48e-08, -4.48e-08, -4.48e-08, 0.003281305, 0.00164063, 0.003280565, 0.003282537, -4.48e-08, 0.003281286, 0.004921239, 0.00491897, -1.52e-06, 0.004925867, -1.27e-06, 0.009839438, 0.001643726, -4.48e-08, 2.09e-06, -0.001640027, 0.006559627, 0.001636958, 2.36e-06, 0.003281469, 0.011481469, 0.004922279, 0.018044207, 0.011483134, 0.014765087, 0.008201379, 0.00492497, 0.006560482, 0.009844796, 0.011483199, 0.008202129, 0.001641621, 0.004921645, 0.006563377, 0.006561068, 0.008201004]

tck = interpolate.splrep(x,y) # pass in s= some value to change smoothing: 
                              # higher = smoother, s=0 for no smoothing

xnew = np.arange(-20, 20, 0.1)
ynew = interpolate.splev(xnew,tck,der=0)

现在，xnew 和 ynew 里包含了更细致的拟合结果，你可以用下面的代码找到最大值：

max_index = np.argmax(ynew)
max_value = ynew[max_index]
max_x = xnew[max_index]

抱歉，我没法测试这个；我现在用的电脑上没有安装scipy等库……不过这应该能给你一些思路。

参考资料：http://docs.scipy.org/doc/scipy/reference/tutorial/interpolate.html

Answer 3

我会对数据使用高斯滤波器来进行平滑处理：

# first, make a function to linearly interpolate the data
f = scipy.interpolate.interp1d(x,y)

# resample with 1000 samples
xx = np.linspace(-20,19, 1000)

# compute the function on this finer interval
yy = f(xx)

# make a gaussian window
window = scipy.signal.gaussian(200, 60)

# convolve the arrays
smoothed = scipy.signal.convolve(yy, window/window.sum(), mode='same')

# get the maximum
xx[np.argmax(smoothed)]

这是处理后平滑的结果：

最大值出现在6.93。

在scipy.signal中还有很多其他的窗口函数和滤波选项。想了解更多，可以查看文档。