Python implementation of the Fraser-Suzuki skewed Gaussian function

The Fraser-Suzuki ‘skewed Gaussian’ function (Fraser & Suzuki, Anal. Chem. 1969, 41, 37; Rusch & Lelieur, Anal. Chem. 1973, 45, 1541) is sometimes used to describe the shapes of optical spectra, or chromatographic line shapes. It can be written in the following form.

Here, the function is defined piecewice, since the analytical function in the first line of the definition does not exist for all $x$ , nor for $b = 0$ . Unfortunately, this cannot be fixed. (BTW, it appears, from limited numerical tests, that our piece-wisely defined skewed Gaussian (cases $b \neq 0$ ) is continuous and differentiable in $x$ everywhere.)

Below, I give a Python/numpy implementation of this function. It is ‘vectorised’, i.e. it handles numpy vectors for high-speed computation, but can also handle individual floats as an argument. The skewed Gaussian is parametrised by $A$ (the height of the peak), $x_0$ (the position of the peak), $w$ (proportional to the width of the peak) and $b$ (the ‘skewness’, or asymmetry). It handles the case where the analytical function vanishes (where the argument of the logarithm ceases to be positive). It also handles the case $b =0$ , by switching to the symmetric Gaussian below a certain threshold value of $b$ (where $b$ will considered to be ‘close’ to zero).

This implementation seems to be robust and works for our practical applications, but it may certainly be improved upon. Comments and suggestions are welcome

import numpy as np

def _skew_gauss_vec(x, A, x0, w, b):
	"""vectorised version of skewed Gaussian, called by skew_gauss"""
	ndeps = np.finfo(x.dtype.type).eps
	lim0 = 2.*np.sqrt(ndeps)
	# Through experimentation I found 2*sqrt(machine_epsilon) to be
	# a good safe threshold for switching to the b=0 limit
	# at lower thresholds, numerical rounding errors appear
	if (abs(b) <= lim0): 
		sg = A * np.exp(-4*np.log(2)*(x-x0)**2/w**2)
	else:
		lnterm = 1.0 + ((2*b*(x-x0))/w)
		sg = np.zeros_like(lnterm)
		sg[lnterm>0] =\
			A * np.exp(-np.log(2)*(np.log(lnterm[lnterm>0])/b)**2)
	return sg

def skew_gauss(x, A, x0, w, b):
	"""Fraser-Suzuki skewed Gaussian.

	A: peak height, x0: peak position, 
	w: width, b: skewness"""
	if type(x)==np.ndarray:
		sg = _skew_gauss_vec(x, A, x0, w, b)
	else:
		x = float(x)
		ndeps = np.finfo(type(x)).eps
		lim0 = 2.*np.sqrt(ndeps)
		if (abs(b) <= lim0):
			sg = A * np.exp(-4*np.log(2)*(x-x0)**2/w**2)
		else:
			lnterm = 1.0 + ((2*b*(x-x0))/w)
			if (lnterm>0):
				sg = A * np.exp(-np.log(2)*(np.log(lnterm)/b)**2)
			else:
				sg = 0
	return sg