# Python - Differentiating Cubic Spline numerically or analytically

I'm trying to take data in numpy arrays, interpolate them and retrieve this function to differentiate/integrate at different stages.

Here is how I got here, data included:

import numpy as np

x = np.array([  50.,   55.,   60.,   65.,   70.,   75.,   80.,   85.,   90.,
95.,  100.,  105.,  110.,  115.,  120.,  125.,  130.,  135.,
140.,  145.,  150.,  155.,  160.,  165.,  170.,  173.,  174.,
175.,  176.,  177.,  178.,  179.,  180.,  181.,  182.,  183.,
184.,  185.,  186.,  187.,  188.,  189.,  190.,  191.,  192.,
193.,  194.,  195.,  196.,  197.,  198.,  199.,  200.,  201.,
202.,  203.,  204.,  205.,  206.,  207.,  208.,  209.,  210.,
211.,  212.,  213.,  214.,  215.,  216.,  217.,  218.,  219.,
220.,  221.,  222.,  223.,  224.,  225.,  226.,  227.,  228.,
229.,  230.,  231.,  232.,  233.,  234.,  235.,  236.,  237.,
238.,  239.,  240.,  245.,  250.,  255.,  260.,  265.,  270.,
275.,  280.,  285.,  290.,  295.,  300.,  305.,  310.,  315.,
320.,  325.])

y = np.array([  1.00000000e-02,   1.00000000e-02,   1.00000000e-02,
1.00000000e-02,   1.00000000e-02,   1.00000000e-02,
1.00000000e-02,   1.50000000e-02,   2.00000000e-02,
2.50000000e-02,   3.50000000e-02,   4.00000000e-02,
5.00000000e-02,   6.00000000e-02,   8.00000000e-02,
1.05000000e-01,   1.30000000e-01,   1.70000000e-01,
2.20000000e-01,   2.75000000e-01,   3.50000000e-01,
4.50000000e-01,   5.75000000e-01,   7.30000000e-01,
9.40000000e-01,   1.08500000e+00,   1.14000000e+00,
1.19500000e+00,   1.25500000e+00,   1.32000000e+00,
1.38500000e+00,   1.45000000e+00,   1.52500000e+00,
1.60000000e+00,   1.68000000e+00,   1.76500000e+00,
1.84500000e+00,   1.94000000e+00,   2.03500000e+00,
2.13500000e+00,   2.24000000e+00,   2.35000000e+00,
2.47500000e+00,   2.59000000e+00,   2.72500000e+00,
2.85000000e+00,   2.99000000e+00,   3.15000000e+00,
3.29500000e+00,   3.46500000e+00,   3.63000000e+00,
3.81500000e+00,   4.00500000e+00,   4.20500000e+00,
4.41500000e+00,   4.62500000e+00,   4.86500000e+00,
5.10500000e+00,   5.36000000e+00,   5.62000000e+00,
5.89500000e+00,   6.18000000e+00,   6.49000000e+00,
6.80500000e+00,   7.14000000e+00,   7.49500000e+00,
7.87000000e+00,   8.26000000e+00,   8.68500000e+00,
9.12500000e+00,   9.58500000e+00,   1.00800000e+01,
1.06000000e+01,   1.11500000e+01,   1.17300000e+01,
1.23400000e+01,   1.29850000e+01,   1.36650000e+01,
1.43700000e+01,   1.51150000e+01,   1.58900000e+01,
1.67000000e+01,   1.75450000e+01,   1.84200000e+01,
1.93150000e+01,   2.02200000e+01,   2.11600000e+01,
2.20800000e+01,   2.30300000e+01,   2.39950000e+01,
2.49600000e+01,   2.59400000e+01,   2.69150000e+01,
3.18400000e+01,   3.68050000e+01,   4.17550000e+01,
4.67450000e+01,   5.17150000e+01,   5.66950000e+01,
6.16750000e+01,   6.66600000e+01,   7.16350000e+01,
7.66150000e+01,   8.15900000e+01,   8.65700000e+01,
9.15500000e+01,   9.65300000e+01,   1.01520000e+02,
1.06490000e+02,   1.11470000e+02])


With some assistance, I got to this point.

import matplotlib.pyplot as plt
from scipy.interpolate import UnivariateSpline

y_spl = UnivariateSpline(x,y,s=0,k=4)
x_range = np.linspace(x[0],x[-1],1000)

plt.plot(x,y,'ro',label = 'data')
plt.plot(x_range,y_spl(x_range))


This seems like a reasonably good fit. I want to take this polynomial and retrieve its results symbolically, because my numerical derivatives appear unnatural, even for this spline.

y_spl_2d = y_spl.derivative(n=2)
plt.plot(x_range,y_spl_2d(x_range))


I think even visually, this second derivative just doesn't cut it. The initial function appears pretty "smooth" in a math-sense.

What are my options from here? Can I retrieve the initial function symbolically and perform the derivative this way (analytically?). What other numerical procedures can I implement to smooth my second derivative? It doesn't appear to me at any point that my second derivative should be negative (and in a modeling sense with what I'm physically modeling, it should never be negative).

Can I retrieve the polynomial from the fit, or can I make numerical modifications to smooth my second derivative?

Edit: the goal is to retrieve a probability density function. I would also like to restrict this final derivative to be strictly positive and integrate to 1, but I think my problems are bigger than this right now. Those are secondary issues.

Here is the first derivative, which looks like a very fine-looking CDF, so I'm becoming more and more sure that my problems are numerical. Should I interpolate this result, as well, and differentiate? I'm sure these problems are very common - what are my next steps.

y_spl_2d = y_spl.derivative(n=1)
plt.plot(x_range,y_spl_2d(x_range), 'ro')


## 1 Answer

The scipy spline interpolation routine can create a smoothed spline that doesn't exactly interpolate the given points but which trades off smoothness against how closely it interpolates noisy points. You could turn up the smoothing to get a more stable result.