对于分段函数,一次执行多次curve_fit迭代

Question

我正在尝试同时执行Scipy的多次迭代,curve_fit以避免循环,从而提高速度.

这与这个问题非常相似,已经解决了.然而,功能是分段(不连续)的事实使得该解决方案不适用于此.

考虑这个例子:

import numpy as np
from numpy import random as rng
from scipy.optimize import curve_fit
rng.seed(0)
N=20
X=np.logspace(-1,1,N)
Y = np.zeros((4, N))
for i in range(0,4):
    b = i+1
    a = b
    print(a,b)
    Y[i] = (X/b)**(-a) #+ 0.01 * rng.randn(6)
    Y[i, X>b] = 1

这产生了这些数组:

你可以看到哪些是不连续的X==b.我可以通过迭代检索原始值a并b使用curve_fit:

def plaw(r, a, b):
    """ Theoretical power law for the shape of the normalized conditional density """
    import numpy as np
    return np.piecewise(r, [r < b, r >= b], [lambda x: (x/b)**-a, lambda x: 1])


coeffs=[]
for ix in range(Y.shape[0]):
    print(ix)
    c0, pcov = curve_fit(plaw, X, Y[ix])
    coeffs.append(c0)

但这个过程可能非常缓慢依赖的大小X,Y以及循环,所以我想通过努力去加快速度coeffs,而不需要为一个循环.到目前为止,我没有运气.

可能很重要的事情:

• X并且Y只包含正值
• a而且b总是积极的
• 尽管适合于该示例的数据是平滑的(为了简单起见),但实际数据具有噪声

编辑

这是我得到的:

y=np.ma.masked_where(Y<1.01, Y)

lX = np.log(X)
lY = np.log(y)
A = np.vstack([lX, np.ones(len(lX))]).T
m,c=np.linalg.lstsq(A, lY.T)[0]

print('a=',-m)
print('b=',np.exp(-c/m))

但即使没有任何噪音,输出也是:

a= [0.18978965578339158 1.1353633705997466 2.220234483915197 3.3324502660995714]
b= [339.4090881838179 7.95073481873057 6.296592007396107 6.402567167503574]

这比我希望得到的更糟糕.

Answer 1

两个建议：

1. 使用numpy.where（也可能argmin）找到数据变为 1 或略大于 1 时的 X值，并将数据截断到该点——有效地忽略 Y=1 的数据。Y

这可能是这样的：

index_max = numpy.where(y < 1.2)[0][0]
x = y[:index_max]
y = y[:index_max]

2. 使用双对数图中显示的提示，幂律现在在双对数中呈线性。您不需要curve_fit，但可以在vsscipy.stats.linregress上使用。对于您的实际工作，这至少将为后续的配合提供良好的起始值。log(Y)log(Y)

跟进此问题并尝试解决您的问题，您可以尝试以下操作：

import numpy as np 
from scipy.stats import linregress

np.random.seed(0)
npts = 51 
x = np.logspace(-2, 2, npts)
YTHRESH = 1.02

for i in range(5):
    b = i + 1.0 + np.random.normal(scale=0.1)
    a = b + np.random.random()
    y = (x/b)**(-a) + np.random.normal(scale=0.0030, size=npts)
    y[x>b] = 1.0

    # to model exponential decay, first remove the values
    # where y ~= 1 where the data is known to not decay...
    imax = np.where(y < YTHRESH)[0][0]

    # take log of this truncated x and y
    _x = np.log(x[:imax])
    _y = np.log(y[:imax])

    # use linear regression on the log-log data:
    out = linregress(_x, _y)

    # map slope/intercept to scale, exponent
    afit = -out.slope
    bfit = np.exp(out.intercept/afit)

    print(""" === Fit Example {i:3d}
          a  expected {a:4f}, got {afit:4f}
          b  expected {b:4f}, got {bfit:4f}
          """.format(i=i+1, a=a, b=b, afit=afit, bfit=bfit))

希望这足以让您继续前进。