Fully monotone interpolation in python

Question

I have some data such as:

I would like to fit a differentiable monotone curve to it. I tried PchipInterpolator class but on a very similar graph it resulted in:

This is not monotone.

How can I fit a monotone curve to data like this?

Here is a sample set of y values for another similar graph:

[0.1109157119023644, 0.20187393816931934, 0.14466318670239758, 0.16535159414166822, 0.05452708697483864, 0.2153046237959556, 0.2200300476272603, 0.21012762463269324, 0.15947100322395022, 0.2819691842129948, 0.15567770052985092, 0.24850595803020692, 0.1329341593280457, 0.15595107081606913, 0.3232021121832229, 0.23707961921686588, 0.2415887076540357, 0.32363506549779797, 0.3584089204036798, 0.29232772580068433, 0.22145994836140775, 0.22797587985241133, 0.2717787840603025, 0.3245255944762287, 0.29301098282789195, 0.32417076823344143, 0.3450906550996232, 0.34272097408024904, 0.3868714875012437, 0.41876692320045755, 0.3544198724867363, 0.33073960954801895, 0.3921033666371904, 0.33349050060172974, 0.3608862044547096, 0.37375822841635425, 0.5396399750708429, 0.4209201143798284, 0.42004773793166883, 0.5217725632679073, 0.5911731474218788, 0.43389609315065386, 0.4287288396176006, 0.43007525393257007, 0.5687062142675405, 0.6030811498722173, 0.5292225577714743, 0.47710974351051355, 0.6182720730381119, 0.6241033581931327, 0.6236788197617511, 0.6643161356364049, 0.5577616524049582, 0.6888440258481371, 0.6867893120660341, 0.6685257606057502, 0.599481675493677, 0.7309075091448749, 0.7644365338580481, 0.6176797601816733, 0.6751467827192018, 0.6452178017908761, 0.6684778262246701, 0.7003380077556168, 0.667035916425416, 0.8434451759113093, 0.8419343615815968, 0.8657695361433773, 0.7392487161484605, 0.8773282098364621, 0.8265679895117846, 0.7246599961191632, 0.7251899061730714, 0.9271640780410231, 0.9180581424305536, 0.8099033021701689, 0.8268585329594615, 0.8519967080830176, 0.8711231413093845, 0.8689802343798663, 0.8299523829217353, 1.0057741699770046, 0.8538130788729608, 0.9662784297225102, 1.023419780920539, 0.913146849759822, 0.9900885996579213, 0.8740638988529978, 0.8900285618419457, 0.9065474574434158, 1.0749522597307315, 1.0319120938258166, 1.0051369663172995, 0.9893558841613622, 1.051384986916457, 1.0327996870915341, 1.0945543972861898, 0.9716604944496021, 1.1490370559566179, 1.1379231481207432, 1.6836433783615088, 1.8162068766097395, 2.072155286917785, 2.0395966998366, 2.191064589600466, 2.1581974932543617, 2.163403843819597, 2.133441151300847, 2.1726053994136922, 2.1157865673629526, 2.2249636455682866, 2.2313062166802147, 2.1731708496472764, 2.315203950110816, 2.1601242661726827, 2.174940281421225, 2.2653635413275945, 2.337227057574145, 2.3645767548381618, 2.3084919291392527, 2.314014515926446, 2.25166717296155, 2.2621157708115778, 2.2644578546265586, 2.313504860292943, 2.398969190357051, 2.309443951779675, 2.278946047410807, 2.4080802287121146, 2.353652872018618, 2.35527529074088, 2.4233001060410784, 2.428767198055608, 2.35677123091093, 2.497135132404064, 2.3978099128437282, 2.3970802609341972, 2.4967434818740024, 2.511209192435555, 2.541001050440798, 2.5760248002036525, 2.5960512284192245, 2.4778408861721037, 2.5757724103530046, 2.631148267999664, 2.538327346218921, 2.4878734713248507, 2.6133797275761066, 2.6282561527857395, 2.6150327104952447, 3.102757164382848, 3.3318503012160905, 3.3907776288198193, 3.6065313558941936, 3.601180295875859, 3.560491539319038, 3.650095006265445, 3.574812155815713, 3.686227315374108, 3.6338261415040867, 3.5661194785086288, 3.5747332336054645, 3.560674343726918, 3.5678550481603635, 3.5342848534390967, 3.4929538312485913, 3.564544653619436, 3.6861775399566126, 3.6390300636595216, 3.6656336332413666, 3.5731185631923945, 3.5965520044069854, 3.537434489989021, 3.5590937423870144, 3.5331656424410083, 3.640652819618705, 3.5971240740252126, 3.641793843012055, 3.6064014089254295, 3.530378938786505, 3.613631139461306, 3.519542268056021, 3.5416251524576, 3.524789618934195, 3.5519951806099512, 3.6435695455293975, 3.6825670484650863, 3.5993379768209217, 3.628367553897596, 3.633290480934276, 3.5772841681579535, 3.602326323397947, 3.518180278272883, 3.531054006706696, 3.5566645495066167, 3.5410992153240985, 3.630762839301216, 3.5924649123201053, 3.646230633817883, 3.568290612034935, 3.638356129262967, 3.566083243271712, 3.6064978645771797, 3.4942864293427633, 3.595438454812999, 3.681726879126678, 3.6501308156903463, 3.5490717955938593, 3.598535359345363, 3.6328331698421654, 3.595159538698094, 3.556715819008055, 3.6292942886764554, 3.6362895697392856, 3.5965220100874093, 3.6103542985016266, 3.5715010140382493, 3.658769915445062, 3.5939686395400416, 3.4974461928859917, 3.5232691556732267, 3.6145687814416614, 3.5682054018341005, 3.648937250575395, 3.4912089018613384, 3.522426560340423, 3.6757968409374637, 3.651348691084845, 3.5395070091675973, 3.5306275536360383, 3.6153498246329883, 3.599762785949876, 3.5351931286962333, 3.6488316987683054, 3.5198301490992963, 3.5696570079786687, 3.561553836008927, 3.5659475947331423, 3.553147100256108, 3.5475591872743664, 3.6097226797553317, 3.6849600324757934, 3.5264731043844413, 3.506658609738451, 3.5535775980874114, 3.5487291053913554, 3.570651383823912, 3.552993371839188, 3.5054297764661846, 3.5723024888238792]

Answer 1

Here's a monotone curve fitter in essentially 5 lines of python, with numpy and a lowpass filter from scipy.signal:

#!/usr/bin/env python
"""https://stackoverflow.com/questions/56551114/fully-monotone-interpolation-in-python """
# see also
# https://en.wikipedia.org/wiki/Monotone-spline aka I-spline
# https://scikit-learn.org/stable/modules/isotonic.html
# denis 2 March 2020

from __future__ import division, print_function
import numpy as np
from scipy import signal as sig

from matplotlib import pyplot as plt
import seaborn as sns


def butter_filtfilt( x, Wn=0.5, axis=0 ):
    """ butter( 2, Wn ), filtfilt
        axis 0 each col, -1 each row
    """
    b, a = sig.butter( N=2, Wn=Wn )
    return sig.filtfilt( b, a, x, axis=axis, method="gust" )  # twice, forward backward

def ints( x ):
    return x.round().astype(int)

def minavmax( x ):
    return "min av max %.3g %.3g %.3g" % (
            x.min(), x.mean(), x.max() )

def pvec( x ):
    n = len(x) // 25 * 25
    return "%s \n%s \n" % (
        minavmax( x ),
        ints( x[ - n : ]) .reshape( -1, 25 ))

#...............................................................................
def monofit( y, Wn=0.1 ):
    """ monotone-increasing curve fit """
    y = np.asarray(y).squeeze()
    print( "\n{ monofit: y %d %s  Wn %.3g " % (
        len(y), minavmax( y ), Wn ))
    ygrad = np.gradient( y )
    print( "grad y:", pvec( ygrad ))

        # lowpass filter --
    gradsmooth = butter_filtfilt( ygrad, Wn=Wn )
    print( "gradsmooth:", pvec( gradsmooth ))

    ge0 = np.fmax( gradsmooth, 0 )

    ymono = np.cumsum( ge0 )  # integrate, sensitive to first few
    ymono += (y - ymono).mean()

    err = y - ymono
    print( "y - ymono:", pvec( err ))
    errstr = "average |y - monofit|: %.2g" % np.abs( err ).mean()
    print( errstr )
    print( "} \n" )

    return ymono, err, errstr

#...............................................................................
if __name__ == "__main__":
    import sys

    np.set_printoptions( threshold=20, edgeitems=15, linewidth=120,
            formatter = dict( float = lambda x: "%.2g" % x ))  # float arrays %.2g
    print( 80 * "=" )

    thispy = sys.argv[0]
    infile = sys.argv[1] if len(sys.argv) > 1 \
        else "so-mono.txt"
    Wn = 0.1  # ?
    params = "%s %s  Wn %g " % (thispy, infile, Wn)
    print( params )

    y = np.loadtxt( infile ) * 100
    print( "y:", y )

    ymono, err, errstr = monofit( y, Wn=Wn )

    if 1:
        sns.set_style("whitegrid")
        fig, ax = plt.subplots( figsize=[10, 5] )
        plt.subplots_adjust( left=.05, right=.99, bottom=.05, top=.90 )
        fig.suptitle(
"Easy monotone curve fit: np.gradient | lowpass filter | clip < 0 | integrate \n"
            + errstr, multialignment="left" )

        ax.plot( ymono, color="orangered" )
        j = np.where( ymono < y )[0]
        xax = np.arange( len(y) )
        plt.vlines( xax[j], ymono[j], y[j], color="blue", lw=1 )
        j = np.where( ymono > y )[0]
        plt.vlines( xax[j], y[j], ymono[j], color="blue", lw=1 )

        png = thispy.replace( ".py", ".png" )
        print( "writing", png )
        plt.savefig( png )
        plt.show()