python中的完全单调插值 [英] Fully monotone interpolation in python

问题描述

我有一些数据，例如:

我想为它拟合一条可微的单调曲线.我尝试了

这不是单调的.

如何为这样的数据拟合单调曲线?

这是另一个类似图形的 y 值样本集:

<预> <代码> [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]

解决方案

这里有一个基本上是 5 行 python 的单调曲线拟合器，带有 numpy和来自 scipy.signal 的低通滤波器:

#!/usr/bin/env python"""https://stackoverflow.com/questions/56551114/fully-monotone-interpolation-in-python"""# 也可以看看# https://en.wikipedia.org/wiki/Monotone-spline 又名 I-spline# https://scikit-learn.org/stable/modules/isotonic.html# 丹尼斯 2020 年 3 月 2 日从 __future__ 导入师，print_function将 numpy 导入为 np从 scipy 导入信号作为 sig从 matplotlib 导入 pyplot 作为 plt将 seaborn 作为 sns 导入def butter_filtfilt( x, Wn=0.5, axis=0 ):""" 黄油 ( 2, Wn ), filtfilt轴 0 每列，-1 每行"""b, a = sig.butter( N=2, Wn=Wn )return sig.filtfilt( b, a, x, axis=axis, method="gust" ) # 两次，向前向后定义整数( x ):返回 x.round().astype(int)def minavmax( x ):返回最小 av 最大 %.3g %.3g %.3g" % (x.min(), x.mean(), x.max() )def pvec( x ):n = len(x)//25 * 25返回 "%s \n%s \n" % (minavmax( x ),ints( x[ - n : ]) .reshape( -1, 25 ))#………………………………………………………………………………………………......................................def monofit( y, Wn=0.1 ):"""单调递增曲线拟合"""y = np.asarray(y).squeeze()打印(\n{ monofit:y %d %s Wn %.3g" % (len(y), minavmax( y ), Wn ))ygrad = np.gradient( y )打印(grad y:"，pvec(ygrad))# 低通滤波器  -gradsmooth = butter_filtfilt( ygrad, Wn=Wn )打印(梯度平滑:"，pvec(梯度平滑))ge0 = np.fmax(梯度平滑，0)ymono = np.cumsum( ge0 ) # 整合，对前几个敏感ymono += (y - ymono).mean()错误 = y - ymono打印(y - ymono:"，pvec(错误))errstr = "平均值 |y - monofit|: %.2g" % np.abs( err ).mean()打印( errstr )打印(}\n")返回 ymono, err, errstr#………………………………………………………………………………………………......................................如果 __name__ == "__main__":导入系统np.set_printoptions( threshold=20, edgeitems=15, linewidth=120,formatter = dict( float = lambda x: "%.2g" % x )) # 浮点数组 %.2g打印(80 *=")thispy = sys.argv[0]infile = sys.argv[1] 如果 len(sys.argv) >1 \否则so-mono.txt"Wn = 0.1 # ?params = "%s %s Wn %g" % (thispy, infile, Wn)打印(参数)y = np.loadtxt( infile ) * 100打印(y:"，y)ymono, err, errstr = monofit( y, Wn=Wn )如果 1:sns.set_style("whitegrid")图, ax = plt.subplots( figsize=[10, 5] )plt.subplots_adjust( left=.05, right=.99, bottom=.05, top=.90 )fig.suptitle(简单的单调曲线拟合:np.gradient | 低通滤波器 | 剪辑 <0 | 整合 \n"+ errstr, multialignment="left" )ax.plot(ymono，颜色=橙色")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")打印(写作"，png )plt.savefig( png )plt.show()

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]

解决方案

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()

这篇关于python中的完全单调插值的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！