Zhaitan - 1 year ago 80
Python Question

# Fitting piecewise function: running into NaN problems due to boolean array indexing

So I've been trying to fit to an exponentially modified gaussian function (if interested, 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
``````

In order to prevent overflow problems, one of two equivilent functions must be used depending on whether z is positive or negative (see Alternative forms for computation from previous wikipedia page).

The problem I am having is this: The line
`y += (k2 * np.exp(n2) * sps.erfcx(z)) * (z >= 0)`
is only supposed to add to y when z is positive. But if z is, say, -30,
`sps.erfcx(-30)`
is
`inf`
, and inf * False is NaN. Therefore, instead of leaving y untouched, the resulting y is clustered with NaN. Example:

``````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])
``````

I tried the replacing the line in question with the following:

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

But doing this ran into serious runtime issues. Is there a way to only evaluate
`(k2 * np.exp(n2) * sps.erfcx(z))`
on the condition that
`(z >= 0)`
? Is there some other way to solve this without sacrificing runtime?

Thanks!

EDIT: After Rishi's advice, the following code seems to work much better:

``````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))
``````

How about using numpy.where with something like: `np.where(z >= 0, sps.erfcx(z), sps.erfc(z))`. I'm no numpy expert, so don't know if it's efficient. Looks elegant at least!