Paul Lambert
  • Home
  • Publications
  • Software & Tutorials
  • Talks
  • Courses
  • Interactive graphs

Derivatives and integrals of spline functions.

One advantage of using spline functions is that we get analytical expressions for the derivative and integral of the spline function.

I will illustrate this using B-splines. I first generate some data and plot it.

. clear

. set obs 1000
Number of observations (_N) was 0, now 1,000.

. range x 0 5

. gen y = 2 + sin(x) + rnormal(0,0.1)

. scatter y x, color(%30) ylabel(0(0.5)3.5, format(%3.1f))

We use B-splines to fit a non-linear function to this data and then obtain the derivative and integral of this function.

If our spline function estimate f(x), we can generate the derivative, f′(x) by obtaining the derivative of each spline variable w.r.t. x.

In gensplines this is done by including the dgen() option.

. gensplines x, type(bs) df(4) gen(_bs) dgen(_dbs)

. global knots `r(knots)'

. regress y _bs*

      Source |       SS           df       MS      Number of obs   =     1,000
-------------+----------------------------------   F(4, 995)       =  12637.11
       Model |  501.188667         4  125.297167   Prob > F        =    0.0000
    Residual |  9.86543954       995  .009915015   R-squared       =    0.9807
-------------+----------------------------------   Adj R-squared   =    0.9806
       Total |  511.054106       999  .511565672   Root MSE        =    .09957

------------------------------------------------------------------------------
           y | Coefficient  Std. err.      t    P>|t|     [95% conf. interval]
-------------+----------------------------------------------------------------
        _bs1 |   1.179382   .0307769    38.32   0.000     1.118987    1.239777
        _bs2 |   1.422978   .0239687    59.37   0.000     1.375943    1.470013
        _bs3 |  -1.309876   .0273988   -47.81   0.000    -1.363642    -1.25611
        _bs4 |  -.8256583   .0198255   -41.65   0.000     -.864563   -.7867537
       _cons |    1.92711   .0150208   128.30   0.000     1.897634    1.956586
------------------------------------------------------------------------------

. predict mu
(option xb assumed; fitted values)

Below I add the fitted regression line to the scatter plot

. twoway (scatter y x, color(%30))         ///
>        (line mu x, sort),                ///
>        ylabel(0(0.5)3.5, format(%3.1f))  ///
>        legend(off)

To generate the dervative of the fitted function we can multipy the coefficients by the variables created by the dgen() option, the derivative of each spline variable w.r.t. x..

. gen double deriv1 = _b[_bs1]*_dbs1 + _b[_bs2]*_dbs2 + _b[_bs3]*_dbs3 + _b[_bs4]*_dbs4      

As a comparison I calculate the derivatives numerically using the dydx command.

. dydx mu x, gen(deriv2) double

. compare deriv1 deriv2

                                        ---------- Difference ----------
                            Count       Minimum      Average     Maximum
------------------------------------------------------------------------
deriv1<deriv2                 514     -2.23e-06    -6.57e-09   -1.11e-16
deriv1=deriv2                   1
deriv1>deriv2                 485      1.39e-17     2.21e-08    8.32e-06
                       ----------
Jointly defined              1000     -2.23e-06     7.35e-09    8.32e-06
                       ----------
Total                        1000

We can plot the function and can see that the derivative is zero at the turning points of the original function.

. line deriv2 x

In order to calculate the integral of the function we need to generate the integral of each of the spline variables. This can be done by using the type(ibs) option of gensplines. If we use df(4) as above, then the knots will be in exactly the same location. Alternatively, we could have saved the knots when we generated the original spline variables and used the allknots() option.

. gensplines x, type(ibs) df(4) gen(_ibs) 

We need to integrate the intercept and then incude the integrated spline functions and multiply each by the associated coefficient.

. gen double integ1 = _b[_cons]*x + _b[_bs1]*_ibs1 + _b[_bs2]*_ibs2 + _b[_bs3]*_ibs3 + _b[_bs
> 4]*_ibs4      

We can compare to the numerically integrated function.

. integ mu x, gen(integ2) double

number of points = 1000

integral         = 10.73512

. compare integ1 integ2

                                        ---------- Difference ----------
                            Count       Minimum      Average     Maximum
------------------------------------------------------------------------
integ1<integ2                   1     -1.38e-11    -1.38e-11   -1.38e-11
integ1=integ2                   1
integ1>integ2                 998      1.05e-12     2.82e-12    8.18e-12
                       ----------
Jointly defined              1000     -1.38e-11     2.80e-12    8.18e-12
                       ----------
Total                        1000

And plot the integral at each value of x.

. line integ1 x

If we are just interested in the integral at one value of x, we just pass this value to gensplines. For example, the integral at x=5 is,

. gensplines 5, type(ibs) allknots(${knots}) gen(_k) 

. lincom _b[_cons]*5 + _b[_bs1]*_k1 + _b[_bs2]*_k2 + _b[_bs3]*_k3 + _b[_bs4]*_k4

 ( 1)  1.25*_bs1 + 1.25*_bs2 + 1.25*_bs3 + .625*_bs4 + 5*_cons = 0

------------------------------------------------------------------------------
           y | Coefficient  Std. err.      t    P>|t|     [95% conf. interval]
-------------+----------------------------------------------------------------
         (1) |   10.73512   .0157442   681.85   0.000     10.70422    10.76602
------------------------------------------------------------------------------

Note we also get a confidence interval for the integral as I used lincom.

Adding the integral of the intercept can be avoided by generating an extra spline variable and then using the nonconstant option when fitting the model. Below I do this and then obtain the integral in a frame. By naming the integrated B-spline variables the same as the spline variables, I can use the predict command to obtain the integrated function.

. drop _*

. gensplines x, type(bs) df(4) gen(_bs) intercept

. regress y _bs*, nocons

      Source |       SS           df       MS      Number of obs   =     1,000
-------------+----------------------------------   F(5, 995)       >  99999.00
       Model |  5108.18252         5   1021.6365   Prob > F        =    0.0000
    Residual |  9.86543954       995  .009915015   R-squared       =    0.9981
-------------+----------------------------------   Adj R-squared   =    0.9981
       Total |  5118.04796     1,000  5.11804796   Root MSE        =    .09957

------------------------------------------------------------------------------
           y | Coefficient  Std. err.      t    P>|t|     [95% conf. interval]
-------------+----------------------------------------------------------------
        _bs1 |    1.92711   .0150208   128.30   0.000     1.897634    1.956586
        _bs2 |   3.106492   .0191994   161.80   0.000     3.068816    3.144168
        _bs3 |   3.350088   .0261145   128.28   0.000     3.298842    3.401334
        _bs4 |   .6172346   .0191994    32.15   0.000     .5795586    .6549105
        _bs5 |   1.101452   .0150208    73.33   0.000     1.071976    1.130928
------------------------------------------------------------------------------

.  
. cap frame drop integral

. frame put y x, into(integral)

. frame integral {
.   gensplines x, type(ibs) df(4) gen(_bs) intercept
.   predict mu,
(option xb assumed; fitted values)
.   line mu x, ytitle(Integral of fitted values)
. }