总结:
我有一个函数,我想通过曲线拟合,这是分段的,但它的片段没有一个干净的分析形式。我已经通过包含一个有点简陋的for
循环让这个过程正常工作了,但是它运行得非常慢:对于N=10000的相对数据集来说,运行时间是1-2分钟。我希望就如何(a)使用numpy广播来加速操作(无for
循环)或(b)做一些完全不同的事情来获得相同类型的结果,但要快得多的建议。你知道吗
我使用的函数z=f(x,y; params)
在y
中逐段单调递增,直到它在y_0
点达到最大值f
,然后饱和并变为常数。我的问题是断点y_0
不是解析的,因此需要一些优化。这将相对容易,除了我在x
和y
中都有真实的数据,并且断点是拟合参数c
的函数。所有数据x
、y
和z
都有检测噪声。你知道吗
示例问题: 下面的函数已被更改,以便更容易说明我要处理的问题。是的,我知道它是可以解析解的,但我真正的问题是不是。
f(x,y; c) =
y*(x-c*y), for y <= x/(2*c)
x**2/(4*c), for y > x/(2*c)
断点y_0 = x/(2*c)
是通过取f
WRT y
的导数并求最大值得到的。最大值f_max=x**2/(4*c)
是通过将y_0
放回f
找到的。问题是断点既依赖于x
-值,也依赖于拟合参数c
,因此我无法在内部循环之外计算断点。你知道吗
代码 我已经将点数减少到~500点,以允许代码在合理的时间内运行。我的真实数据超过了10000点。你知道吗
import numpy as np
from scipy.optimize import curve_fit,fminbound
import matplotlib.pyplot as plt
def function((x,y),c=1):
fxn_xy = lambda x,y: y*(x-c*y)
y_max = np.zeros(len(x)) #create array in which to put max values of y
fxn_max = np.zeros(len(x)) #array to put the results
'''This loop is the first part I'd like to optimize, since it requires
O(1/delta) time for each individual pass through the fitting function'''
for i,X in enumerate(x):
'''X represents specific value of x to solve for'''
fxn_y = lambda y: fxn_xy(X,y)
#reduce function to just one variable (y)
#by inputting given X value for the loop
max_y = fminbound(lambda Y: -fxn_y(Y), 0, X, full_output=True)
y_max[i] = max_y[0]
fxn_max[i] = -max_y[1]
return np.where(y<=y_max,
fxn_xy(x,y),
fxn_max
)
''' Create and plot 'real' data '''
delta = 0.01
xs = [0.5,1.,1.5,2.] #num copies of each X value
y = []
#create repeated x for each xs value. +1 used to match size of y, below
x = np.hstack([X]*(int(X//delta+1)) for X in xs)
#create sweep from 1 to the current value of x, with spacing=delta
y = np.hstack([np.arange(0, X, delta) for X in xs])
z = function((x,y),c=0.75)
#introduce random instrumentation noise to x,y,z
x += np.random.normal(0,0.005,size=len(x))
y += np.random.normal(0,0.005,size=len(y))
z += np.random.normal(0,0.05,size=len(z))
fig = plt.figure(1,(12,8))
axis1 = fig.add_subplot(121)
axis2 = fig.add_subplot(122)
axis1.scatter(y,x,s=1)
#axis1.plot(x)
#axis1.plot(z)
axis1.set_ylabel("x value")
axis1.set_xlabel("y value")
axis2.scatter(y,z, s=1)
axis2.set_xlabel("y value")
axis2.set_ylabel("z(x,y)")
'''Curve Fitting process to find optimal c'''
popt, pcov = curve_fit(function, (x,y),z,bounds=(0,2))
axis2.scatter(y, function((x,y),*popt), s=0.5, c='r')
print "c_est = {:.3g} ".format(popt[0])
结果如下所示,带有“真实”x、y、z值(蓝色)和拟合值(红色)。
注:我的直觉是找出如何广播x
变量s.t。我可以在fminbound中使用它。但这可能是幼稚的。思想?你知道吗
谢谢大家!你知道吗
编辑:澄清一下,x
-值在这样的组中并不总是固定的,而是可以随着y
-值保持稳定而进行扫描。不幸的是,这就是为什么我需要一些方法来处理x
那么多次。你知道吗
有几件事可以优化。一个问题是数据结构。如果我正确地理解了代码,您可以查找所有
x
的max
。但是,您创建了这样的结构:相同的值被重复了一次又一次。因此,在这里您将浪费大量的计算工作。你知道吗我不确定
f
的评估在现实中有多困难,但我假设它的成本并不比优化要高很多。所以在我的解决方案中,我只是计算整个数组,寻找最大值,然后改变后面的值。你知道吗我想我的代码也可以优化,但现在看来:
提供
原始图形和噪声数据
相关问题 更多 >
编程相关推荐