如何将pygam模型转换为SciPy BSpline

Question

我想在两个世界中取长补短: Scipy和pygam。具体地说，我想使用自定义的pygam惩罚函数惩罚我的数据的三阶导数，然后我想使用结果作为Scipy BSpline。后者允许容易地计算数据的n阶导数。当然，gam.coef_给了我模型的系数，但是结点是什么呢？如果使用gam.terms.get_coef_indices()-1来接近，但结果实际上应该与gam.predict(x)相同，即

在ev-br的有用评论之后，我们现在知道我们可以通过pygam.utils.get_knot_edges获得结边。因此，一段有希望的代码将是

plt.figure()
res0=pygam.LinearGAM(terms=pygam.terms.s(0),n_splines=len(x),penalties=MyPenalty,lam=0.3,fit_intercept=False, basis='ps').fit(x,y,1/np.sqrt(w))
plt.plot(x,y,'x',zorder=10, label='data')
knot_edges=utils.gen_edge_knots(x,dtype='numerical')
knots=np.linspace(knot_edges[0],knot_edges[-1],len(res0.coef_))
plt.plot(x,res0.predict(x), label='PyGAM fit')
plt.plot(x,scipy.interpolate.BSpline(knots,res0.coef_,3)(x), label='BSpline from PyGAM parameters')
   
    
plt.legend()

遗憾的是，输出结果并不令人满意：

​

Answer 1

根据https://github.com/dswah/pyGAM/issues/265的说法，pyfam适合统一的结。可以通过pygam.utils.gen_edge_knots提取第一个和最后一个结的位置。

Answer 2

不幸的是，B样条方法没有标准化(我在过去的40年里开发了五种不同的B样条拟合算法，没有一种算法与其他算法兼容！)，仔细检查pygam和scipy的结果表明，不幸的是不兼容。首先，节点要么被移动和缩放(scipy splprep)，要么没有(pygam)，第二，由于不同的计算方法，B样条坐标和系数是不兼容的。尽管它们之间的转换可能是可能的，但这可能会被证明是几乎难以处理的。

另一点是，如果你想使用平滑，splprep提供了这一点(同意，它可能不像pygam方法那样容易和优雅)。

这里是一种方法(使用数值微分)，与一组树木生长数据的spldev导数相比，该方法可能适用于pygam。

import numpy as np
import scipy
from scipy.interpolate import interp1d, splprep, splev
from scipy.sparse import csr_matrix
import matplotlib.pyplot as plt
import pandas as pd
import math

from pygam import LinearGAM, s, f, utils
import pygam as pygam

X = [1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24]

Y = [5,5,5.045,5.234,5.369,5.558,5.693,
     5.819,5.864,5.882,5.882,5.882,5.882,
     5.882,5.906,6.074,6.398,6.674,6.842,
     6.998,7.07,7.082,7.082,7.082
    ]


plt.figure()
# Plot data
plt.plot(X,Y,'x',zorder=10, label='data')

# First, pygam
GAM=pygam.LinearGAM(terms=pygam.terms.s(0),n_splines=len(X),
                     lam=0.3,fit_intercept=False,
                     basis='ps').fit(X,Y)

knot_edges=utils.gen_edge_knots(X,dtype='numerical')
knots=np.linspace(knot_edges[0],knot_edges[-1],len(GAM.coef_))
print('GAM Knots')
print(knots)
coeffs = GAM.coef_
print('GAM Coefficients')
print(coeffs)
GAMY = GAM.predict(X)
plt.plot(X,GAMY, label='PyGAM fit')
print('GAMY fit')
print(GAMY)

# Second, splprep
tck, u = splprep([X,Y],s=0) # force interpolation by s=0
# where t = knots (scaled 0 to 1)
#       c = coefficients
#       k = degree of spline (default = 3)
#       u = An array of the values of the parameter.
new_points = splev(u,tck, der=0)
plt.plot(new_points[0], new_points[1], '.', color='k', label='splprep fit')
print('tck')
print(tck)
print('u')
print(u)

#Ynew = scipy.interpolate.BSpline(knots,res0.coef_,3)
#plt.plot(X,Ynew)
         #(X,label='BSpline from PyGAM parameters'))
   
plt.legend()
plt.show()

# Derivatives, splev
dev_points_x, dev_points_y = splev(u, tck, der=1) # first derivative
plt.plot(u*23+1, dev_points_y/23, '--ro', label='splev derivative')
#u and y must be scaled and shifted

# Deriv, GAM (by numerical difference method)
XX = GAM.generate_X_grid(term=0, n=200)
#now do derivative(s)
Yprime = np.zeros(200) ## Y prime
GAMY = GAM.predict(XX) # fitted Y on XX
for i in range(1,len(Yprime)-1):
    Yprime[i] = (GAMY[i+1] - GAMY[i-1])/(XX[i+1] - XX[i-1])
# compute central difference on XX[i]
# now fix up Yp1[0] and Yp1[end]
#Yp1[0] = Yp1[1] # just substitution
Yprime[0] =  (GAMY[1] - GAMY[0])/(XX[1] - XX[0])  # forward difference
#Yp1[len(Yp1)-1] = Yp1[len(Yp1) - 2] # just substitution
last = len(Yprime)-1
Yprime[last] = (GAMY[last] - GAMY[last-1])/(XX[last] - XX[last-1]) # backward diff


plt.plot(XX,Yprime,'-b',label='GAM numerical deriv' )

plt.legend()
plt.show()

​

​

​

​