拟合分段函数：因布尔数组索引而遇到NaN问题

Question

所以我一直在尝试拟合一个指数修正的高斯函数（如果有兴趣，https://en.wikipedia.org/wiki/Exponentially_modified_Gaussian_distribution）

import numpy as np
import scipy.optimize as sio
import scipy.special as sps

def exp_gaussian(x, h, u, c, t):
    z = 1.0/sqrt(2.0) * (c/t - (x-u)/c)   #not important
    k1 = k2 = h * c / t * sqrt(pi / 2)    #not important
    n1 = 1/2 * (c / t)**2 - (x-u)/t       #not important
    n2 = -1 / 2 * ((x - u) / c)**2        #not important
    y = np.zeros(len(x))
    y += (k1 * np.exp(n1) * sps.erfc(z)) * (z < 0)
    y += (k2 * np.exp(n2) * sps.erfcx(z)) * (z >= 0)
    return y

为了防止溢出问题，必须使用两个等价函数中的一个，这取决于z是正的还是负的（请参阅上一页wikipedia中的可选计算形式）。你知道吗

我遇到的问题是：当z为正时，线y += (k2 * np.exp(n2) * sps.erfcx(z)) * (z >= 0)应该只加到y上。但是如果z是，比如说，-30，sps.erfcx(-30)是inf，inf*False是NaN。因此，结果y不是保持y不变，而是与NaN聚集在一起。示例：

x = np.linspace(400, 1000, 1001)
y = exp_gaussian(x, 100, 400, 10, 5)
y
array([ 84.27384586,  86.04516723,  87.57518493, ...,          nan,
                nan,          nan])

我尝试用以下内容替换有问题的行：

y += numpy.nan_to_num((k2 * np.exp(n2) * sps.erfcx(z)) * (z >= 0))

但是这样做会遇到严重的运行时问题。有没有办法只在(z >= 0)的条件下计算(k2 * np.exp(n2) * sps.erfcx(z))？有没有其他方法可以在不牺牲运行时间的情况下解决这个问题？你知道吗

谢谢！你知道吗

编辑：根据Rishi的建议，下面的代码似乎工作得更好：

def exp_gaussian(x, h, u, c, t):
    z = 1.0/sqrt(2.0) * (c/t - (x-u)/c)
    k1 = k2 = h * c / t * sqrt(pi / 2)
    n1 = 1/2 * (c / t)**2 - (x-u)/t
    n2 = -1 / 2 * ((x - u) / c)**2
    return = np.where(z >= 0, k2 * np.exp(n2) * sps.erfcx(z), k1 * np.exp(n1) * sps.erfc(z))

Answer 1

您可以做的一件事是创建一个掩码并重用它，这样就不需要对它进行两次评估。另一个想法是在结尾只使用nan\u to\u num一次

mask = (z<0)
y += (k1 * np.exp(n1) * sps.erfc(z)) * (mask)
y += (k2 * np.exp(n2) * sps.erfcx(z)) * (~mask)
y = numpy.nan_yo_num(y)

试试看这是否有帮助。。。你知道吗

Answer 2

把numpy.where和这样的东西一起使用怎么样：np.where(z >= 0, sps.erfcx(z), sps.erfc(z))。我可不是什么numpy专家，所以不知道这是否有效。至少看起来很优雅！你知道吗