找出曲线内的最小值

2024-03-28 12:32:32 发布

您现在位置:Python中文网/ 问答频道 /正文

总结: 我有一个函数,我想通过曲线拟合,这是分段的,但它的片段没有一个干净的分析形式。我已经通过包含一个有点简陋的for循环让这个过程正常工作了,但是它运行得非常慢:对于N=10000的相对数据集来说,运行时间是1-2分钟。我希望就如何(a)使用numpy广播来加速操作(无for循环)或(b)做一些完全不同的事情来获得相同类型的结果,但要快得多的建议。你知道吗

我使用的函数z=f(x,y; params)y中逐段单调递增,直到它在y_0点达到最大值f,然后饱和并变为常数。我的问题是断点y_0不是解析的,因此需要一些优化。这将相对容易,除了我在xy中都有真实的数据,并且断点是拟合参数c的函数。所有数据xyz都有检测噪声。你知道吗

示例问题: 下面的函数已被更改,以便更容易说明我要处理的问题。是的,我知道它是可以解析解的,但我真正的问题是不是。

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)是通过取fWRT 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那么多次。你知道吗


Tags: oftheto函数forlenvaluenp
1条回答
网友
1楼 · 发布于 2024-03-28 12:32:32

有几件事可以优化。一个问题是数据结构。如果我正确地理解了代码,您可以查找所有xmax。但是,您创建了这样的结构:相同的值被重复了一次又一次。因此,在这里您将浪费大量的计算工作。你知道吗

我不确定f的评估在现实中有多困难,但我假设它的成本并不比优化要高很多。所以在我的解决方案中,我只是计算整个数组,寻找最大值,然后改变后面的值。你知道吗

我想我的代码也可以优化,但现在看来:

import numpy as np
from scipy.optimize import leastsq
import matplotlib.pyplot as plt


def function( yArray, x=1, c=1 ):
    out = np.fromiter( ( y * ( x - c * y ) for y in yArray ), np.float )
    pos = np.argmax( out )
    outMax = out[ pos ]
    return np.fromiter( ( x if i < pos else outMax for i, x in enumerate( out ) ), np.float )


def residuals( param, xArray, yList, zList ):
    c = param
    zListTheory = [ function( yArray, x=X, c=c ) for yArray, X in zip( yList, xArray ) ]
    diffList = [ zArray - zArrayTheory for zArray, zArrayTheory in zip( zList, zListTheory ) ]
    out = [ item  for diffArray in diffList for item in diffArray ]
    return out


'''  Create and plot 'real' data   '''
delta = 0.01
xArray = np.array( [ 0.5, 1., 1.5, 2. ] )  #keep this as parameter
#create sweep from 1 to the current value of x, with spacing=delta
yList = [ np.arange( 0, X, delta ) for X in xArray ] ## as list of arrays
zList = [ function( yArray, x=X, c=0.75 ) for yArray, X in zip( yList, xArray ) ]


fig = plt.figure( 1, ( 12, 8 ) )
ax = fig.add_subplot( 1, 1, 1 )
for y,z in zip( yList, zList ):
    ax.plot( y, z )

#introduce random instrumentation noise
yRList =[ yArray + np.random.normal( 0, 0.02, size=len( yArray ) ) for yArray in yList ]
zRList =[ zArray + np.random.normal( 0, 0.02, size=len( zArray ) ) for zArray in zList ]

ax.set_prop_cycle( None )
for y,z in zip( yRList, zRList ):
    ax.plot( y, z, marker='o', markersize=2, ls='' )

sol, cov, info, msg, ier = leastsq( residuals, x0=.9, args=( xArray, yRList, zRList ), full_output=True )
print "c_est = {:.3g} ".format( sol[0] )
plt.show()

提供

>> c_est = 0.752

原始图形和噪声数据 Theory and Data

相关问题 更多 >