I have a data set of points that I need to fit to a normal distribution, where the points approximate the curve of the distribution's second integral. Examples of such curves are given below in the first image, and an example of the data set is given in the second image.

How would one approach curve fitting here? I was thinking of implementing non-linear regression (Gauss-Newton), but that method assumes that the partial derivatives of the target function are known. Is there any clean way to approximate the residuals and jacobian matrix for this problem? Any existing software libraries, or references of interest?

Thank you

second integral of normal distribution


  • 2
    $\begingroup$ Am I understaning you correctly that you want to fit the function $y = \Psi(z)$ where $z=(x-\mu)/\sigma$ and $\Psi(z) = \int^z \Phi(t) dt$ where $\Phi(t)$ is the CDF of a normal distribution? If so you can derive the Jacobian analytically. $\endgroup$
    – Dave
    Aug 31, 2021 at 22:19
  • 2
    $\begingroup$ @Monty It's not clear what difficulty you anticipate with computation of the partial derivatives. $\endgroup$
    – Glen_b
    Sep 1, 2021 at 0:38

2 Answers 2


Let $\Psi(z) = \int_{- \infty}^z \Phi(t) \> dt$. Note $\Phi$ is approximately constant as $\vert z \vert \to \infty$. This means that that $\Psi$ is going to be linear in these same regions.

Because of this behaviour, I think it is reasonable to use a natural spline to approximate $\Psi$. A Natural spline is linear in the tails of the data (just like our target function). The question becomes: where the knots should be placed? I'm sure the knots could be placed intelligently leveraging properties of $\Phi$ and the Gaussian density, but I'll just let the rms library pick them for me today. If you were able to intelligently select knot locations, then you could pass them to rcs using the parms argument.

Using R...


# Generate data
x = seq(-8, 8, 0.25)
f = pnorm(x)
Psi = cumtrapz(x, f)

train_ix = abs(x)<=5
test_ix = abs(x)>5

xtrain = x[train_ix]
xtest = x[test_ix]

ytrain = Psi[train_ix]
ytest = Psi[test_ix]

model = ols(ytrain ~ rcs(xtrain))

plot(x, predict(model, newdata=list(xtrain=x)), col='red', type='l')
points(xtrain, ytrain)
points(xtest, ytest, col='blue', pch=2)

Which produces the following fit

enter image description here

The triangles here are data which were not used to fit the model. Because we have used a natural spline (which is linear in its tails) we do fairly well. Shown below is the absolute error between model and integral (pay attention to those $x$ which are larger than 5 in absolute value, those are the test points)

enter image description here

Because $\Phi$ is not actually linear in the tails (although nearly) our model suffers slightly as evidenced by the growth in absolute error. The relative error is quite good in the right tail, but appears to explode in the left tail since $\Psi$ is quite small in those regions.

As for a model formula, here it is

enter image description here

Apologies for the awful formatting, I the latex function from rms doesn't seem to play nice with typesetting here.


Analytical derivations for applicability to Gauss-Newton are simple after all. I mistakenly thought the derivatives were hard.

For residuals, the cumulative normal distribution function is equal to the Gauss Error function by scaling and translation. The antiderivative is exactly $z\ \text{erf}\ z + \frac{e^{-z^2}}{\sqrt{\pi}}$ and can be derived by integration by parts. Software libraries usually include the normal CDF (and its inverse) as a built-in.

The Jacobian is computable by the chain rule and the fundamental theorem of calculus using again the normal CDF. Set $z = (x - \mu) / \sigma$.

$$\frac{d\Psi(z)}{d\mu} = \frac{d\Psi(z)}{dz} \frac{dz}{d\mu} = - \frac{\Phi(z)}{\sigma}$$ $$\frac{d\Psi(z)}{d\sigma} = \frac{d\Psi(z)}{dz} \frac{dz}{d\sigma} = - \frac{z\ \Phi(z)}{\sigma}$$

Convergence is a non-issue for this class of functions.

  • $\begingroup$ Isn't $dz/d\mu = -1/\sigma$?. Also, your $\Phi(x)$ should be $\Phi(z)$'s. $\endgroup$
    – Dave
    Sep 1, 2021 at 13:47
  • $\begingroup$ I got $d\Psi/d\mu = -\Phi(z)/\sigma$ and $d\Psi/d\mu = -(z/\sigma) \Phi(z)$ (the latter just re-expresses what you've got). $\endgroup$
    – Dave
    Sep 1, 2021 at 13:48
  • $\begingroup$ @Dave You are correct. As you can tell, my calculus is slightly rusty ;). $\endgroup$ Sep 1, 2021 at 15:23

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.