Compute or approximate the mode of a normal distribution with skew parameter $\lambda$. How about for Ashour 2010's approximation? I'd like to compute the mode for a normal distribution with skew parameter λ. I assume I can just differentiate, but is there a shortcut or approximation?
I'd like to learn how to do this generally, for my own statistical knowledge. But, more pressingly, I'm a developer and I'm writing a program that uses the approximation in Ashour and Abdul-hameed (2010) of a skewed normal. I need to at least approximate the mode for the distribution therein.
How could I compute or approximate the mode for an Ashour 2010's approximation?
 A: Although you are interested in computing an approximate mode for the skew-normal distribution, it is not actually all that difficult to compute the exact mode.  This might even be easier than computing the model of the approximation in Ashour and Abdel-hameed (2010), so I recommend that you consider using the exact form.  In this answer I will show you how to characterise the exact mode, which you can then obtain numerically via iterative methods.
The skew-normal distribution has a strictly unimodal density function (i.e., it is strictly quasi-concave) so the mode occurs at the unique critical point of the density function.  This means that the mode is characterised by the critical point equation for the density.  To derive this equation, we first note the following derivatives that apply to the standard normal distribution:
$$\frac{d \phi}{dx}(x) = -x \phi(x)
\quad \quad \quad \quad \quad 
\frac{d \Phi}{dx}(\alpha x) = \alpha \phi(\alpha x).$$
Thus, differentiating the density of the skew-normal density function gives:
$$\begin{align}
\frac{df}{dx} (x) 
&= \frac{d}{dx} 2 \phi(x) \Phi(\alpha x) \\[6pt]
&= 2 \Bigg[ \frac{d\phi}{dx} (x) \Phi(\alpha x) + \phi(x) \frac{d\Phi}{dx} (\alpha x) \Bigg] \\[6pt]
&= 2 \Bigg[ -x \phi (x) \Phi(\alpha x) + \alpha \phi(x) \phi(\alpha x) \Bigg] \\[6pt]
&= -2 \phi(x) \Bigg[ x \Phi(\alpha x) - \alpha \phi(\alpha x) \Bigg]. \\[6pt]
\end{align}$$
Consequently, the mode $\hat{x}$ is the unique solution to the (implicit) critical point equation:
$$\hat{x} \Phi(\alpha \hat{x}) = \alpha \phi(\alpha \hat{x}).$$
You can solve this equation numerically to find the mode.  This involves using an appropriate iterative method (e.g., Newton Raphson, bisection, etc.) to solve the critical point equation.
A: Following Ben's derivations, here is a possible implementation in Python:
from scipy.stats import norm
import scipy
from scipy.stats import skewnorm
import numpy as np
import matplotlib.pyplot as plt

def mode_eq(x, alpha):
    '''
    Defines Ben's first order condition
    '''
    return x*norm.cdf(alpha*x) - alpha*norm.pdf(alpha*x)
alpha  = 2
sol = scipy.optimize.fsolve(mode_eq, 0, args=(alpha,), fprime=None, full_output=0, 
                      col_deriv=0, xtol=1.49012e-08, maxfev=0, band=None, 
                           epsfcn=None, factor=100, diag=None)
sol

Giving the solution 0.53075814 (for $\alpha = 2$). This is calculated in around 2.05ms in my machine. Here is a graph of the skew-normal with $\alpha = 2$:

Indeed, the unique mean seems close to $0.5$.
