Я хочу, чтобы двумерный кубический сплайн соответствовал некоторым нерегулярно расположенным данным - то есть функции, которая точно соответствует данным в заданных точках - но также может возвращать значения между ними.
Все, что я могу найти (для нерегулярных разнесенных данных), это scipy.interpolate.SmoothBivariateSpline
.Я не могу понять, как отключить «сглаживание» (независимо от того, какое значение я задаю в параметре s
.
Однако я обнаружил, что могу получить в основном то, что хочу, с помощью scipy.interpolate.griddata
- хотя для этого нужно каждый раз пересчитывать его (т. е. не просто генерировать функцию). Есть ли какая-то разница, принципиально между этими двумя - то есть griddata
делает что-то отличное от «сплайна»? Есть ли способ отключить?сглаживание в SmoothBivariateSpline
или эквивалентной функции, которая не сглаживает?
Ниже приведен скрипт, который я использую для проверки подгонки сплайна к полиному
import numpy as np
from mpl_toolkits.mplot3d import Axes3D
import scipy.optimize
import scipy.interpolate
import matplotlib.pyplot as plt
import numpy.polynomial.polynomial as poly
# Grid and test function
N = 9;
x,y = np.linspace(-1,1, N), np.linspace(-1,1, N)
X,Y = np.meshgrid(x,y)
F = lambda X,Y : X+Y-1*X*Y-(X*Y)**2 -2*X*Y**2 + X**2*Y + 3*np.exp(-((X+1)**2+(Y+1)**2)*5)
Z = F(X,Y)
noise = 0.4
Z *= 1+(np.random.random(Z.shape)*2-1)*noise # noise
# Finer Grid and test function
N2 = 19;
x2,y2 = np.linspace(-1,1, N2), np.linspace(-1,1, N2)
X2,Y2 = np.meshgrid(x2,y2)
Z2 = F(X2,Y2)
# Make data into lists
Xl = X.reshape(X.size)
Yl = Y.reshape(Y.size)
Zl = Z.reshape(Z.size)
# Polynomial fit
# polyval(x,y,p) = p[0,0]+p[0,1]y+p[1,0]x+p[1,1]xy+p[1,2]xy^2 ..., etc
# I use a flat (1D) array for p, so it needs to be reshaped into a 2D array before
# passing to polyval
order = 3
p0 = np.zeros(order**2) # guess parameters (all 0 for now)
f_poly = lambda x,y,p : poly.polyval2d(x,y,p.reshape((order,order))) # Wrapper for our polynomial
errf = lambda p : np.mean((f_poly(Xl,Yl,p.reshape((order,order)))-Zl)**2) # error function to find least square error
sol = scipy.optimize.minimize(errf, p0)
psol = sol['x']
# Spline interpolation
# Bivariate (2D), Smoothed (doesn't fit points *exactly*) cubic (3rd order - i.e. kx=ky=3) spline
spl = scipy.interpolate.SmoothBivariateSpline(Xl, Yl, Zl, kx=3,ky=3)
f_spline = spl.ev
# regular Interpolate
f_interp = lambda x,y : scipy.interpolate.griddata((Xl, Yl), Zl, (x,y), method='cubic')
# Plot
fig = plt.figure(1, figsize=(7,8))
plt.clf()
# poly fit
ax = fig.add_subplot(311, projection='3d')
ax.scatter3D(X2,Y2,Z2,s=3, color='red', label='actual data')
fit = f_poly(X2,Y2, psol)
l = 'order {} poly fit'.format(order)
ax.plot_wireframe(X2,Y2, fit, color='black', label=l)
ax.scatter3D(X,Y,Z, color='blue', label='noisy data')
plt.legend()
print("Average {} error: {}".format(l, np.sqrt(np.mean((fit-Z2)**2))))
# spline fit
ax = fig.add_subplot(312, projection='3d')
ax.scatter3D(X2,Y2,Z2,s=3, color='red', label='actual data')
l = 'smoothed spline'
fit = f_spline(X2,Y2)
ax.plot_wireframe(X2,Y2, fit, color='black', label=l)
ax.scatter3D(X,Y,Z, color='blue', label='noisy data')
plt.legend()
print("Average {} error: {}".format(l, np.sqrt(np.mean((fit-Z2)**2))))
# interp fit
ax = fig.add_subplot(313, projection='3d')
ax.scatter3D(X2,Y2,Z2,s=3, color='red', label='actual data')
l='3rd order interp '
fit=f_interp(X2,Y2)
ax.plot_wireframe(X2,Y2, fit, color='black', label=l)
ax.scatter3D(X,Y,Z, color='blue', label='noisy data')
plt.legend()
print("Average {} error: {}".format(l, np.sqrt(np.mean((fit-Z2)**2))))
plt.show(False)
plt.pause(1)
raw_input('press key to continue') # Change to input() if using python3