IDL skewed gaussian program incorrect

[M16]

Greetings IDL comrades,

I need to fit some data using a skewed normal/gaussian distribution and I'm not convinced that the IDL program gauss1 works correctly. Very likely it's my ignorance, but it seems to be the IDL routine gauss1 is not correct. I've attached a gif of a few results with the value of skew ranging from 0 to .75. The code I wrote is below. The equation for a skewed gaussian that I find on wiki and other sites is not like what is in the function gauss1.

Am I missing something here?

Thanks very much,

Mark

\;
;   p = [2.2D, 1.4D, 3000.D]
;   x = dindgen(200)*0.1 - 10.
;   y = gauss1(x, p)
;
;   Computes the values of the Gaussian at equispaced intervals
;   (spacing is 0.1).  The gaussian has a mean of 2.2, standard
;   deviation of 1.4, and total area of 3000.
;
; REFERENCES:
;
; MODIFICATION HISTORY:
;   Written, Jul 1998, CM
;   Correct bug in normalization, CM, 01 Nov 1999
;   Optimized for speed, CM, 02 Nov 1999
;
;-

function gauss1, x, p, skew=skew, _EXTRA=extra

p= [2.2D,1.4D,3000.D]

   x = dindgen(200)*0.1 - 10.
   
   
  sz = size(x)
  if sz(sz(0)+1) EQ 5 then smax = 26D else smax = 13.

  if n_elements(p) GE 3 then norm = p(2) else norm = x(0)*0 + 1

  u = ((x-p(0))/(abs(p(1)) > 1e-20))^2
  mask = u LT (smax^2)
  f = norm * mask * exp(-0.5*temporary(u) * mask) / (sqrt(2.D * !dpi)*p(1))
  mask = 0

  if n_elements(skew) GT 0 then $
    f = (1.D + skew * (x-p(0))/p(1))*f

return,f
end




pro trash2,xxx
start
!p.multi=[0,2,2]

for m=0.,5.,.25 do begin
plot,gauss1(skew=m)
if !p.multi[0] eq 0 then subgo

endfor

return
end

[David Grier]

Hi Mark,

You're right that gauss1 is wrong.  This is obvious because the distribution function goes

negative.  The expression for the skew-normal distribution should be

P_a(z) = (1 + erf(a z/sqrt(2))) P(z)

where P(z) is the usual Gaussian distribution and erf(x) is the error function.  Here, a is the

desired skewness, which appears in the argument of the error function.

The expression in gauss1 is what you would get from the first term in the Taylor

expansion of the erf:

P_a(z) = (1 + 2 a z / sqrt(!pi)) P(z),

which would be appropriate for very small values of a.  Note that gauss1 absorbs the

factor of 2/sqrt(!pi) into the definition of the skewness, which also is annoying.

I'd recommend rolling your own skew-normal function using IDL's built-in erf.

All the best,

David