Skip to content

csaps.CubicSmoothingSpline provides incorrect derivatives under certain circumstances #48

@shusheer

Description

@shusheer

Certain arrangements of smoothing intensity, weights and xdata, in particular irregularly spaced data where some xvalues are rather close together, causes CubicSmoothingSpline to produce a first derivative that is not the integral of the second derivative. This can be exacerbated by certain weights. In particular, the first derivative can be very non-smooth, even though the second derivative is smooth.

This is important, because a very common use-case for CubicSmoothingSpline (as opposed to say LOWESS smoothing) is to ensure smooth derivatives.

There are at least two possible solutions:
1 - modifying x-values or combining close points, for some definition of close (might be very data-specific).
2 - Sequentially integrating derivatives from the 3rd upwards, to enforce correct 1st and 2nd derivatives (might be quite computationally expensive).
3 - Converting to compressed scipy.CubicSpline form using code discussed at #46 (comment) also avoids this problem, as the spline is effectively resampled at optimally distributed x values.

Code demonstrating the problem (apologies for the large amount of data in the test case - it's some real data from a project I am working on, and I haven't been able to determine exactly which features are causing the problem, though enforcing a minimum x spacing seems to work as shown, and removing weights reduces the issue but doesn't eliminate it).

As you can see, the resulting splines all give quite different estimates of the data, so what we mean by "correct" is up for discussion.

# Problem case giving unsmooth derivatives
# The first derivative of the smoothed function is NOT the integral of the second derivative
# (notice small discontinuities / steps)
# This is particularly acute when certain weight values play a part as well

x = [-288.55699074,-288.28206019,-287.66591435,-287.64053241,-286.65069444,
-286.65068866,-285.52694444,-285.23972222,-284.64372685,-284.44784722,
-283.38444444,-283.38365741,-282.55511574,-282.53914352,-281.63871528,
-280.46857639,-280.46855324,-278.46340278,-278.37118056,-277.58601852,
-277.58574074,-276.54689815,-276.54663194,-275.62396991,-274.54616898,
-274.54615741,-273.55909722,-273.38333333,-272.57217593,-272.5721412,
-269.7065625,-269.51540509,-268.61980324,-268.60930556,-267.57458333,
-267.57457755,-263.58122685,-259.53576389,-258.82375,-258.82333333,
-257.85988426,-257.85986111,-254.25375,-254.25374421,-253.60337963,
-253.60333333,-251.7546875,-251.75461806,-250.94400463,-247.44850694,
-244.78422454,-244.7834838,-243.80513889,-243.80423611,-242.545625,
-242.54511574,-241.58532407,-241.58527778,-240.53894676,-240.53893519,
-239.59284722,-239.59278935,-234.86810185,-234.86706019,-229.50280093,
-229.50270833,-228.67150463,-228.67148148,-227.83733796,-227.63991898,
-225.75958333,-225.75953704,-224.77252315,-224.77251157,-223.98796296,
-223.98793981,-222.53775463,-222.53773148,-221.47532407,-219.68340278,
-219.6830787,-218.77277778,-218.77275463,-216.28890046,-216.28887731,
-215.52909722,-214.40041667,-214.40039352,-212.61553241,-212.61550926,
-211.81768519,-211.47722222,-211.47716435,-208.37016204,-208.36606481,
-207.62565972,-206.73474537,-206.73467593,-205.8547338,-205.83864583,
-204.94078704,-204.85533565,-203.72707176,-203.72493056,-203.72258102,
-203.57979167,-203.57978009,-202.72645833,-202.68732639,-202.54818287,
-201.68228009,-201.51018519,-201.43820602,-200.49405093,-200.49201389,
-199.62027778,-199.50150463,-199.50037037,-199.49247685,-198.74945602,
-198.74351852,-197.64726852,-197.6459375,-197.48726852,-196.49546296,
-196.43032407,-195.52611111,-195.46309028,-194.50974537,-194.50868056,
-193.46177083,-193.46171296,-192.88958333,-192.88957176,-191.64532407,
-191.64497685,-190.5394213,-190.51228009,-189.69407407,-189.69158565,
-188.55555556,-188.43287037,-186.77581019,-186.76670139,-184.90331019,
-184.90259259,-182.46954861,-182.46900463,-181.82518519,-180.79622685,
-180.79261574,-179.48037037,-179.47877315,-178.50547454,-178.50518519,
-177.54893519,-177.5484838,-176.70898148,-176.70897569,-172.47497685,
-171.58761574,-171.58756944,-169.23037037,-169.19318287,-168.54819444,
-165.79189815,-165.79189236,-165.79136574,-164.7587037,-164.75511574,
-161.85717593,-161.85717014,-160.74731481,-160.74730903,-156.4134838,
-156.41269676,-155.81989583,-155.81988426,-154.45269676,-154.43328704,
-153.80079861,-153.80032407,-152.54888889,-152.54819444,-148.93758102,
-147.84275463,-147.84274884,-146.98821759,-145.84561343,-145.84560185,
-140.25047454,-140.25046296,-138.99819444,-138.99818866,-137.92916667,
-137.85449074,-136.48898148,-136.48894676,-124.4325,-122.78548611,
-122.78547454,-121.78831019,-121.7624537,-119.5171875,-119.51571759,
-118.54518519,-118.54511574,-116.81631944,-116.6050463,-115.41570602,
-115.41196759,-114.53983796,-114.51606481,-112.65020833,-112.59722222,
-111.58196759,-111.58173611,-110.28418981,-108.88050926,-108.88037037,
-106.58741898,-106.58740741,-105.82334491,-105.61048611,-102.78446759,
-102.78175926,-101.92847222,-101.92841435,-99.70208333,-99.66606481,
-98.69592593,-98.69592014,-97.68861111,-96.54909722,-96.43082176,
-93.6450463,-93.64504051,-91.82571759,-91.82571181,-83.88231481,
-83.88003472,-82.53726852,-81.8415162,-81.78796296,-78.41912037,
-77.88599537,-77.88598958,-75.8,-75.79999421,-74.91194444,
-74.91193866,-70.52842593,-70.52751157,-66.66107639,-66.49625,
-64.62689815,-64.62497685,-63.84877315,-62.67493056,-62.67486111,
-56.45991898,-56.45849537,-55.80101852,-53.5209375,-53.42854167,
-51.79412037,-51.70625,-50.72030093,-50.72027778,-49.55359954,
-48.56030093,-48.52893519,-47.61050926,-47.61006944,-45.62118056,
-45.61965278,-43.70581019,-43.70578704,-42.61976852,-42.61907407,
-41.67361111,-41.64385417,-32.35585648,-31.6496875,-31.64862269,
-17.80746528,-17.77924769,-16.75157407,-16.75155093,-15.57109954,
-15.57106481,-14.64717593,-13.72175926,-13.67280093,-12.61773148,
-10.63041667,-10.63005787,-9.80740741,-9.80736111,-8.71253472,
-8.70666667,-6.74260417,-6.74258102,-5.34083333,-4.91016204,
-4.90943287,-2.30743056,-2.30638889,-0.42960648,-0.42840278]

y = [70.7572706,70.54032874,70.05606702,70.03624935,69.27432557,
69.27432117,68.43243561,68.22067996,67.78516637,67.64308661,
66.88118275,66.88062559,66.30052443,66.28947425,65.67543639,
64.90642245,64.90640759,63.66449628,63.60904412,63.14039473,
63.14022978,62.52515116,62.52499365,61.97652441,61.32001027,
61.32000308,60.68789985,60.57033966,60.00261623,60.00259097,
57.61865841,57.43806534,56.55456891,56.54383346,55.43956239,
55.43955595,50.38210045,44.53858607,43.50952278,43.50892299,
42.12805664,42.12802358,36.99275737,36.99274908,36.05906926,
36.05900272,33.41368763,33.41358922,32.27359692,27.70135708,
24.84323913,24.84254672,23.98532718,23.98459096,23.0642553,
23.0639276,22.5150554,22.51503233,22.08034126,22.08033742,
21.8389586,21.8389484,23.10692951,23.10771641,29.81761462,
29.81776653,31.21601939,31.21605924,32.68017101,33.03425222,
36.54083419,36.54092348,38.47670459,38.47672766,40.0590836,
40.05913083,43.0771602,43.07720928,45.35765585,49.31879272,
49.31952077,51.38303931,51.38309234,57.16432204,57.16437592,
58.92167454,61.46687347,61.46692457,65.24326265,65.24330938,
66.81838801,67.46953204,67.46964169,72.89983173,72.90644953,
74.08206155,75.44675245,75.44685685,76.74713302,76.77050437,
78.05420313,78.17438359,79.7183499,79.72118421,79.72429379,
79.91227229,79.91228745,80.99326008,81.04109716,81.20997602,
82.219385,82.41175016,82.49142149,83.49504942,83.49713412,
84.36266061,84.47724037,84.47833159,84.4859241,85.18879314,
85.19431565,86.18839601,86.18957211,86.32923784,87.17818962,
87.23249526,87.96788031,88.01785229,88.75344863,88.75424891,
89.51806878,89.51810973,89.91613619,89.91614411,90.73453049,
90.73474985,91.40125092,91.41664377,91.85156726,91.85279549,
92.34222988,92.38575227,92.75900826,92.75984105,92.59792734,
92.59772702,91.20024149,91.19975484,90.56381405,89.28641717,
89.28134229,87.14597333,87.14300865,85.15861765,85.15797276,
82.84855607,82.84738518,80.54827418,80.54825753,66.28762042,
63.00718054,63.00700893,54.31636812,54.18091963,51.84551574,
42.2968528,42.29683372,42.29509707,38.97123725,38.95998918,
30.75854662,30.75853234,28.19489508,28.19488261,21.31146526,
21.31063361,20.72423984,20.72422916,19.66675016,19.65459906,
19.30037669,19.30014055,18.82393793,18.82375279,18.89283954,
19.25883033,19.25883262,19.63216496,20.23324209,20.2332487,
24.12129249,24.12130103,25.03004562,25.03004973,25.76826801,
25.81806455,26.67629851,26.67631891,30.46561539,30.74262981,
30.74263179,30.91672453,30.92128596,31.30818311,31.3084222,
31.4588991,31.45890927,31.69236101,31.71975256,31.87221841,
31.87267866,31.9724062,31.97484795,32.10544092,32.10723373,
32.11864839,32.11864583,32.06417118,31.90691934,31.90689831,
31.39692865,31.3969252,31.14906127,31.07254683,29.76210069,
29.76059886,29.26593574,29.26590077,27.78880075,27.76287183,
27.0440895,27.04408511,26.26627309,25.36699318,25.27332511,
23.11069468,23.11069035,21.80609638,21.80609244,18.29594004,
18.295638,18.20182827,18.21706169,18.22000001,18.82965191,
18.98616378,18.98616554,19.69017271,19.6901748,20.01523736,
20.0152395,21.42764295,21.42784833,21.84871481,21.84950598,
21.78327154,21.78313985,21.72088105,21.59779516,21.59778708,
20.92961501,20.92955839,20.91266214,21.04288478,21.05619122,
21.41830913,21.44484473,21.79485625,21.79486561,22.3360831,
22.90340539,22.92289613,23.5358956,23.53620858,25.13688862,
25.13825808,27.00134805,27.00137215,28.16778118,28.16854978,
29.24170683,29.27634617,42.20507624,43.32162108,43.32331588,
68.1726824,68.22889458,70.29050421,70.29055095,72.69085058,
72.69092166,74.59109794,76.51130697,76.6133282,78.82157314,
83.02355351,83.0243163,84.77591304,84.77601173,87.11188547,
87.12441636,91.32441849,91.32446806,94.32811947,95.25177267,
95.25333675,100.83484899,100.83708178,104.85695,104.85952725]

w = [0.09462211,0.04377652,0.14838558,0.05943898,0.15541495,0.13085324,
0.062704,0.06243595,0.08839641,0.08066833,0.08655245,0.13400375,
0.12432887,0.09684014,0.08021506,0.1265312,0.05262315,0.06449225,
0.0758878,0.08259978,0.0806998,0.18695698,0.08123137,0.07486425,
0.12497271,0.10635125,0.15783934,0.08400718,0.24832489,0.12591398,
0.75945646,0.13893085,0.07963481,0.10749733,0.12091365,0.06358539,
0.01520488,0.03701155,0.16559939,0.21220943,0.08161835,0.16941486,
0.08949668,0.13363251,0.09692045,0.18369284,0.13672927,0.08379181,
0.01858722,0.01133573,0.06237293,0.24676197,0.10859054,0.12052321,
0.13203374,0.14501103,0.14774643,0.09328596,0.04971943,0.15257333,
0.28338405,0.09501884,0.98374631,0.05226721,0.09065792,0.13678137,
0.09321165,0.44103697,0.06999414,0.05989127,0.12692356,0.0915057,
0.11936501,0.10480617,0.1567755,0.07513431,0.07117161,0.06858016,
0.01640047,0.29806503,0.22542,0.12794823,0.21382192,0.80778319,
0.05041481,0.02423374,0.14427077,0.10418953,0.21508381,0.10940715,
0.04888315,0.04528401,0.04778175,0.11278799,0.43211842,0.03932949,
0.11803702,0.07236406,0.09973482,0.11787877,0.10316486,0.03396932,
0.07494178,0.0980055,0.15682156,0.39430807,0.2966579,0.0321651,
0.04839777,0.09632106,0.08455145,0.03593526,0.09416399,0.1499288,
0.20404766,0.05393556,0.05599547,0.08508938,0.03143724,0.03516901,
0.48995939,0.22142592,0.05018331,0.02789187,0.02197713,0.03036933,
0.20079765,0.0914527,0.04390079,0.0632086,0.12006129,0.04469976,
0.10379572,0.06173869,0.04413667,0.03245433,0.1633229,0.0436747,
0.06350709,0.05263706,0.04127021,0.04639606,0.06203611,0.0542154,
0.04355763,0.0831866,0.03767251,0.04846966,0.01901781,0.24288383,
0.05686747,0.0987101,0.11008607,0.0864383,0.17024615,0.23189691,
0.13231802,0.11141714,0.11234526,0.04201772,0.10555506,0.23622534,
0.15991181,0.07504075,0.08703826,0.078157,0.07949321,0.29557024,
0.0676124,0.16409294,0.24457436,0.10823184,0.14073518,0.35402206,
0.13634875,0.08048864,0.11152011,0.0714971,0.11998616,0.07697284,
0.08883651,0.08903601,0.12036493,0.07058305,0.01155317,0.08867653,
0.21117371,0.05630087,0.0809992,0.27670459,0.26026519,0.13855529,
0.15380278,0.23979833,0.20920358,0.1440551,0.13541747,0.1463744,
0.04642054,0.24201843,0.11214505,0.18320577,0.04986542,0.19388245,
0.07125025,0.22815885,0.08084224,0.06743765,0.28559653,0.11051846,
1.71124858,0.10948372,0.1427056,0.14208709,0.21579165,0.14818066,
0.10186744,0.02192575,0.07098707,0.09253821,0.13900954,0.12502331,
0.06113228,0.15379208,0.05778127,0.25567184,0.12133423,0.1361762,
0.06835819,0.17893468,0.07440845,0.05984869,0.02637486,0.06455977,
0.05552673,0.16468979,0.06516525,0.06218418,0.10451724,0.4187673,
0.08585429,0.02827772,0.19387362,0.0435158,0.02646894,0.16485791,
0.12165923,0.11332932,0.06376871,0.1044266,0.059728,0.14428556,
0.10856501,0.16783946,0.12808263,0.06199017,0.112075,0.04430589,
0.17411464,0.09521449,0.05556677,0.10710905,0.16484723,0.09197195,
0.09080481,0.11942555,0.29181033,0.05171834,0.10642797,0.0591445,
0.07264419,0.05472201,0.07241181,0.07863706,0.05298894,0.17034628,
0.13621438,0.12312648,0.20335133,0.10124627,0.08194233,0.08727886,
0.07454792,0.17554491,0.13723729,0.04581896,0.11138872,0.09357495,
0.09733723,0.21097312,0.05863559,0.02409536,0.11577836,0.08534116,
0.02837538,0.09877817,0.03943265,0.09397333,0.16018314,0.09669458,
0.18185679,0.11361261,0.07640622,0.0501189,0.1409333,0.04278513,
0.08166752,0.0946821,0.138107,0.05173292]

xrnd = np.rint(np.array(x)*25.)/25.
xrndv = xrnd[1:]
xrndv[np.flatnonzero(np.diff(xrnd)==0)]+=1./50.
xrndv[np.flatnonzero(np.diff(xrnd)==0)]+=1./100.

xs = np.unique(np.r_[x[0]+np.cumsum(np.repeat(np.diff(x)/4,4)),np.linspace(x[0], x[-1], len(x)*10)])

unweighted = CubicSmoothingSpline(x, y, smooth=5e-06)
weighted = CubicSmoothingSpline(x, y, weights=w, smooth=5e-06)
weightedimprovedx = CubicSmoothingSpline(xrnd, y, weights=w, smooth=5e-06)

def CubicSplineFromCubicSmoothingSpline(smoothing_spline, xrange=None, bc_type='natural', extrapolate=None):
    # If using this to generate a periodic spline from  a csaps spline, data is usually replicated
    # Therefore restrict the xrange to the actual period to be modelled
    x = smoothing_spline.spline.breaks
    y3 = smoothing_spline.spline.coeffs[0,:]*6
    if xrange:
        assert xrange[0] < xrange[1]
        xminidx = max(0,np.flatnonzero(np.r_[xrange[0],x]<=xrange[0])[-1]-1)
        xmaxidx = min(x.size-1,np.flatnonzero(np.r_[x,xrange[1]]>=xrange[1])[0])
        x = x[xminidx:xmaxidx+1].copy()
        x[0]=max(xrange[0],x[0])
        x[-1]=min(xrange[1],x[-1])
        y3 = y3[xminidx:xmaxidx+1]
    idxs_to_use = [0,x.size-1]
    idx_before_zero_crossings = np.where(np.diff(np.sign(y3)))[0]
    idxs_to_use.extend(idx_before_zero_crossings.tolist())
    for aidx,bidx in zip(np.r_[0,idx_before_zero_crossings],np.r_[idx_before_zero_crossings,x.size-1]):
        idx = np.argmax(np.abs(y3[aidx:bidx+1]))+aidx
        idxs_to_use.append(idx)
    x_compressed = np.unique(x[idxs_to_use])
    y_compressed = smoothing_spline(x_compressed)
    if bc_type=='periodic':
        y_compressed[0] = y_compressed[-1] = np.mean(y_compressed[[0,-1]])
    return CubicSpline(x_compressed, y_compressed, bc_type=bc_type, extrapolate=extrapolate)

compressedspline = CubicSplineFromCubicSmoothingSpline(weighted)

plt.plot(xs, unweighted(xs), '-', xs, weighted(xs), '-', xs, weightedimprovedx(xs), '-', xs, compressedspline(xs), '-', x, y, 'o')
plt.show()

for nu in range(1,4):
    fig,(ax1,ax2,ax3,ax4)=plt.subplots(1, 4, figsize=(12,4))
    ax1.plot(xs, unweighted(xs,nu=nu), '-')
    ax1.title.set_text('unweighted nu='+str(nu))
    ax2.plot(xs, weighted(xs,nu=nu), '-')
    ax2.title.set_text('weighted nu='+str(nu))
    ax3.plot(xs, weightedimprovedx(xs,nu=nu), '-')
    ax3.title.set_text('weighted, improved x nu='+str(nu))
    ax4.plot(xs, compressedspline(xs,nu=nu), '-')
    ax4.title.set_text('compressed spline nu='+str(nu))
    plt.show()

Metadata

Metadata

Assignees

No one assigned

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions