Python 非線形回帰と２次元補間の活用事例

きっかけ

ある紙に書かれたグラフがあるのだが、設計を行うため、これを数値化したい。このため、グラフから座標を読み取り非線形回帰で連続関数に近似する必要が生じた。

グラフでは、(x,y)を与え、zの値を読みとる必要があったが、これを自動化するため、２次元補間を活用することにした。

非線形回帰

プログラムと回帰結果のグラフを以下に示す。
グラフ上の青い点を入力値として、与えた式で再現したものが赤線となる。この回帰はscipy.optimize.leastsqでサクッと行なっている。こういうことが簡単にできることがPythonのいいところである。

import numpy as np
import matplotlib.pyplot as plt
from scipy import optimize


def fit(x,y):
    # curve fitting
    def func(param,x,y):
        aa,bb,cc,dd=param[0],param[1],param[2],param[3]
        return aa*np.sinh((6-x)/bb)+1.0+cc*(6-x)**dd-y
    
    param0=[1,1,1,1]
    result=optimize.leastsq(func,param0,args=(x,y))
    aa,bb,cc,dd=result[0][:]
    print('aa,bb,cc,dd={0:10.3e},{1:10.3e},{2:10.3e},{3:10.3e}'.format(aa,bb,cc,dd))
    xx=np.linspace(x[0],x[-1],111,endpoint=True)
    yy=aa*np.sinh((6-xx)/bb)+1.0+cc*(6-xx)**dd
    
    plt.plot(x,y,'o')
    plt.plot(xx,yy,'-')
    fnameF='_fig_test.jpg'
    plt.savefig(fnameF, dpi=100, bbox_inches="tight", pad_inches=0.1)
    plt.show()
    plt.close()


def main():
    x=np.array([0.500,0.750,1.000,1.250,1.500,2.000,2.500,3.000,3.500,4.000,4.500,5.000,6.000])
    y=np.array([1.250,1.190,1.150,1.125,1.105,1.073,1.052,1.035,1.023,1.015,1.008,1.003,1.000])
    fit(x,y)


if __name__ == '__main__': main()

２次元補間

上記手法などを駆使して、以下のようなグラフが得られたとする。

このグラフからは、横軸$\delta$およびパラメータ$\psi$を指定して縦軸$\mu_{in}$あるいは$\mu_{out}$を読み取る必要がある。
グラフ上の赤点は、最後に示すプログラムでの入力および出力を示している。

ここで、グラフ上の曲線の数値データは既知であるとすると、$(\psi, \delta)$を指定して$\mu$を予測する２次元補間の問題となる。

説明のため、$x=\psi$、$y=\delta$、$z=\mu$とおくと、各曲線は、(x,y,z)の値を有する。ここで２次元補間を行うため、下図に示すように、それぞれの曲線の要素を全て連結して新しい配列(X,Y,Z)を作成する。

この既知の配列(X,Y,Z)=(xx0,yy0,zz0)で構成される曲面上の値z_newを(x_new, y_new)を指定して求めることとなる。
このためには、以下のようにscipy.interporate.griddataを用いる。

 z_new=interpolate.griddata((xx0,yy0),zz0,(x_new,y_new),method='cubic')

上式において、
z_new；求めるべきzの配列
x_new；z値を得るための入力値xの配列
y_new；z値を得るための入力値yの配列
xx0, yy0, zz0；補間のための既知の配列

非線形回帰も含めた一連の作業を実行するプログラムを以下に示す。

import numpy as np
import matplotlib.pyplot as plt
from scipy import interpolate
from scipy import optimize
from PIL import Image


def combine(fnameF):
    def get_concat_h(im1, im2):
        dst = Image.new('RGB', (im1.width + im2.width, im1.height))
        dst.paste(im1, (0, 0))
        dst.paste(im2, (im1.width, 0))
        return dst
    def get_concat_v(im1, im2):
        dst = Image.new('RGB', (im1.width, im1.height + im2.height))
        dst.paste(im1, (0, 0))
        dst.paste(im2, (0, im1.height))
        return dst
    im1 = Image.open('_fig_ifactor_1.jpg')
    im2 = Image.open('_fig_ifactor_2.jpg')
    get_concat_h(im1,im2).save(fnameF)


def barrow(x1,y1,x2,y2):
    col='#333333'
    arst='<->,head_width=2,head_length=5'
    aprop=dict(arrowstyle=arst,connectionstyle='arc3',facecolor=col,edgecolor=col,shrinkA=0,shrinkB=0,lw=1)
    plt.annotate('',
        xy=(x1,y1), xycoords='data',
        xytext=(x2,y2), textcoords='data', fontsize=0, arrowprops=aprop)


def sarrow(x1,y1,x2,y2):
    col='#333333'
    arst='->,head_width=3,head_length=10'
    aprop=dict(arrowstyle=arst,connectionstyle='arc3',facecolor=col,edgecolor=col,shrinkA=0,shrinkB=0,lw=1)
    plt.annotate('',
        xy=(x1,y1), xycoords='data',
        xytext=(x2,y2), textcoords='data', fontsize=0, arrowprops=aprop)


def small_fig(nnn,fsz):
    ds,dp,ell,dd=2,0.75,1,0.75
    xmin,xmax=-1.5,1.5
    ymin,ymax=-1,dd+ell+1
    ddx=0.2
    ddy=ddx/(xmax-xmin)*(ymax-ymin)
    plt.axes((0.88-ddx,0.88-ddy,ddx,ddy))
    plt.xlim([xmin,xmax])
    plt.ylim([ymin,ymax])
    plt.gca().set_aspect('equal',adjustable='box')
    plt.axis('off')
    plt.fill([xmin,xmin,xmax,xmax],[ymin,ymax,ymax,ymin],color='#ffffff')
    x=np.array([xmax,dp/2,dp/2,ds/2,ds/2])
    y=np.array([dd,dd,dd+ell,dd+ell,ymax])
    plt.plot(x,y,'-',color='#000000',lw=2)
    plt.plot(-x,y,'-',color='#000000',lw=2)
    plt.plot([xmin,xmax],[0,0],'-',color='#000000',lw=2)
    dds=0.2
    for i in [1,2,3,4]:
        match i:
            case 1: x1,y1,x2,y2,ss,ddh,ddv,=xmin+0.75,0,xmin+0.75,dd,r'$d$',-dds,0
            case 2: x1,y1,x2,y2,ss,ddh,ddv,=xmin+0.75,dd,xmin+0.75,dd+ell,r'$L$',-dds,0
            case 3: x1,y1,x2,y2,ss,ddh,ddv=-dp/2,dd+ell-0.2,dp/2,dd+ell-0.2,r'$Dp$',0,dds
            case 4: x1,y1,x2,y2,ss,ddh,ddv=-ds/2,ymax-0.5,ds/2,ymax-0.5,r'$Ds$',0,dds
        barrow(x1,y1,x2,y2)
        xs,ys=0.5*(x1+x2)+ddh,0.5*(y1+y2)+ddv
        plt.text(xs,ys,ss,ha='center',va='center',fontsize=fsz)
    plt.text(0,-0.5,r'$\psi=(Dp / Ds)^2$',ha='center',va='center',fontsize=fsz)

    if nnn==1: x1,y1,x2,y2,angA,angB=0,dd+0.2*ell,-dp/2-0.2,dd/2,0,90
    if nnn==2: x1,y1,x2,y2,angA,angB=-dp/2-0.2,dd/2,0,dd+0.2*ell,90,0
    cstyle='angle,angleA={0},angleB={1}'.format(angA,angB)
    aprop=dict(arrowstyle='->',lw=3,connectionstyle=cstyle)
    plt.annotate('', xy=(x1,y1),xytext=(x2,y2),arrowprops=aprop)


def fig(nnn,y0,zz,psi,x_new,y_new,z_new):
    fsz=14
    xmin,xmax,dx= 0, 6, 0.5
    ymin,ymax,dy=1.00, 1.30, 0.02
    iw=8
    ih=10
    plt.figure(figsize=(iw,ih),facecolor='w')
    plt.rcParams['font.size']=fsz
    plt.rcParams['font.family']='sans-serif'

    match nnn:
        case 1: tstr=r'Interference factor for in-flow $\mu$_in'
        case 2: tstr=r'Interference factor for out-flow $\mu$_out'
    plt.xlim([xmin,xmax])
    plt.ylim([ymin,ymax])
    plt.xlabel(r'$\delta=L/D_p$')
    plt.ylabel(tstr)
    plt.xticks(np.arange(xmin,xmax+dx,dx))
    plt.yticks(np.arange(ymin,ymax+dy,dy))
    plt.grid(color='#999999',linestyle='solid')
    #plt.gca().set_aspect('equal',adjustable='box')
    plt.title(tstr,loc='left',fontsize=fsz)
    n,m=zz.shape
    k=15
    ds=1
    for i in range(0,n):
        plt.plot(y0,zz[i,:],'-',color='#000080',lw=1.5,clip_on=False)
        x1,y1=y0[k],zz[i,k]
        x2,y2=x1+ds,y1+0.05*ds*2
        sarrow(x1,y1,x2,y2)
        ss=r'$\psi$={0:.1f}'.format(psi[i])
        if nnn==1 and i==0: ss=r'$\psi\leq${0:.1f}'.format(psi[i])
        plt.text(x2,y2,ss,ha='left',va='center',fontsize=fsz)
    xt,yt,zt=x_new[0],y_new[0],z_new[0] # target
    plt.plot([xmin,yt,yt],[zt,zt,ymin],'--',color='#000000',lw=1)
    plt.plot(yt,zt,'o',ms=6,color='#ff0000',clip_on=False)
    s1='* target'
    s2=r'$\psi$={0:.3f}'.format(xt)
    s3=r'$\delta$={0:.3f}'.format(yt)
    if nnn==1: s4=r'$\mu$_in={0:.3f}'.format(zt)
    if nnn==2: s4=r'$\mu$_out={0:.3f}'.format(zt)
    ss=s1+'\n'+s2+'\n'+s3+'\n'+s4
    xs,ys=xmin+0.75*(xmax-xmin),ymin+0.74*(ymax-ymin)
    dic_box={'fc':'#ffffff','ec':'#000000','boxstyle':'Square','lw':1}
    plt.text(xs,ys,ss,ha='left',va='top',fontsize=fsz,bbox=dic_box)

    small_fig(nnn,fsz)

    #plt.tight_layout()
    fnameF='_fig_ifactor_{0}.jpg'.format(nnn)
    plt.savefig(fnameF, dpi=200, bbox_inches="tight", pad_inches=0.1)
    plt.close()
    #plt.show()


def ifactor(psi0,delta0):
    def fit(nnn,y,z):
        # curve fitting
        def func(param,x,y):
            aa,bb,cc,dd=param[0],param[1],param[2],param[3]
            return aa*np.sinh((6-x)/bb)+1.0+cc*(6-x)**dd-y
        param0=[1,1,1,1]
        result=optimize.leastsq(func,param0,args=(y,z))
        aa,bb,cc,dd=result[0][:]
        print('aa,bb,cc,dd={0:10.3e},{1:10.3e},{2:10.3e},{3:10.3e}'.format(aa,bb,cc,dd))
        yy=np.linspace(y[0],y[-1],111,endpoint=True)
        zz=aa*np.sinh((6-yy)/bb)+1.0+cc*(6-yy)**dd
        plt.plot(y,z,'o')
        plt.plot(yy,zz,'-')
        fnameF='_fig{0}.jpg'.format(nnn)
        plt.savefig(fnameF, dpi=100, bbox_inches="tight", pad_inches=0.1)
        plt.close()
        return yy,zz

    def inp_in():
        # input data for in-flow
        psi=np.array([0.4,0.5,0.6,0.7,0.8,0.9,1.0])
        zref=np.array([1.250,1.230,1.193,1.154,1.103,1.055,1.000])
        y=np.array([0.500,0.750,1.000,1.250,1.500,2.000,2.500,3.000,3.500,4.000,4.500,5.000,6.000])
        z=np.array([1.250,1.190,1.150,1.125,1.105,1.073,1.052,1.035,1.023,1.015,1.008,1.003,1.000])
        rr=(zref[:]-1)/(zref[0]-1)
        return psi,rr,y,z

    def inp_out():
        # input data for out-flow
        psi=np.array([0.0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0])
        zref=np.array([1.248,1.220,1.198,1.180,1.158,1.137,1.115,1.092,1.068,1.040,1.000])
        rr=(zref[:]-1)/(zref[0]-1)
        y=np.array([0.500,0.750,1.000,1.250,1.500,2.000,2.500,3.000,3.500,4.000,4.500,5.000,6.000])
        z=np.array([1.248,1.176,1.135,1.113,1.093,1.066,1.043,1.028,1.017,1.010,1.005,1.003,1.000])
        rr=(zref[:]-1)/(zref[0]-1)
        return psi,rr,y,z

    # Calculation of interference factors
    for nnn in [1,2]:
        if nnn==1 and psi0<0.4: psi1=0.4
        else: psi1=psi0
        x_new=np.array([psi1])
        y_new=np.array([delta0])
        # define the curves
        match nnn:
            case 1: psi,rr,y,z=inp_in() # in-flow
            case 2: psi,rr,y,z=inp_out() # out-flow
        y0,z0=fit(nnn,y,z) # curve fitting
        n,m=len(rr),len(y0)
        zz=np.zeros((n,m),dtype=np.float64)
        for i in range(0,n):
            zz[i,:]=(z0[:]-1.0)*rr[i]+1.0
        # Estimate z_new using interpolation (griddata)
        x0=np.zeros(len(y0),dtype=np.float64)
        xx0=x0+psi[0]
        yy0=y0
        zz0=zz[0,:]
        for i in range(1,n):
            xx0=np.concatenate([xx0,x0+psi[i]])
            yy0=np.concatenate([yy0,y0])
            zz0=np.concatenate([zz0,zz[i,:]])
        xxx,yyy=np.meshgrid(xx0,yy0)
        z_new=interpolate.griddata((xx0,yy0),zz0,(x_new,y_new),method='cubic')
        fig(nnn,y0,zz,psi,x_new,y_new,z_new)
        
        match nnn:
            case 1: mu_i=z_new[0]
            case 2: mu_o=z_new[0]
    return mu_i,mu_o


def main():
    psi0,delta0=0.75,2.70
    mu_i,mu_o=ifactor(psi0,delta0)
    print('mu_i,mu_o={0:6.3f},{1:6.3f}'.format(mu_i,mu_o))
    fnameF='fig_ifactor.jpg'
    combine(fnameF)
    
    

#---------------
# Execute
#---------------
if __name__ == '__main__': main()

以　上