用二次多项式拟合函数f（x，y，z）

Question

我正在尝试使用以下二次多项式拟合函数f（x，y，z）：

3d polynomial

某些扭曲的球形表面在三个维度上。问题与固态物理学中有效质量的计算有关。

这里是数据的图片，显示了即使在z方向上的曲率非常低，它的确在所有方向上都呈抛物线下降：3d parabolas

我对与有效质量相对应的系数感兴趣。我有一个xyz坐标数组，它是规则的，并以最大值为中心：

[[ 0.          0.          0.        ]
 [ 0.          0.          0.01282017]
 [ 0.          0.          0.02564034]
 ...
 [-0.05026321 -0.05026321 -0.03846052]
 [-0.05026321 -0.05026321 -0.02564034]
 [-0.05026321 -0.05026321 -0.01282017]]

以及对应的一维标量值数组，每个点一个。围绕此最大值的数据点数范围可以从100到1000。

这是我目前试图用于拟合的代码：

def func(data, mxx, mxy, mxz, myy, myz, mzz):
    x = data[:, 0]
    y = data[:, 1]
    z = data[:, 2]
    return (
        (1 / (2 * mxx)) * (x ** 2)
        + (1 / (1 * mxy)) * (x * y)
        + (1 / (1 * mxz)) * (x * z)
        + (1 / (2 * myy)) * (y ** 2)
        + (1 / (1 * myz)) * (y * z)
        + (1 / (2 * mzz)) * (z ** 2)
    ) + f(0, 0, 0)

energy = data[:, 3]
guess = (mxx, mxy, mxz, myy, myz, mzz)
params, pcov = scipy.optimize.curve_fit(
    func, data, energy, p0=guess, method="trf"
)

其中f（0,0,0）是函数（0，0，0）的值，我使用scipy.interpolate.griddata函数进行检索。

对于这个问题，总的来说，质量应该是负的，并且值在-0.2到-2之间。我正在通过有限差分微分来创建猜测值。

但是，从scipy.interpolate.curve_fit不会得到任何有意义的结果-通常，系数最终以巨大的数字结尾（例如1e9）。我现在完全迷失了。

我在做什么错了（（？

Answer 1

问题之一是您是否适合1/m。尽管从物理学的角度来看这是正确的，但从算法的角度来看却很糟糕。如果拟合算法需要将m的值更改为接近零的符号，则系数会发散。因此，最好拟合mI = 1/m并在以后进行相应的错误处理。在这里，我使用leastsq，它需要对协方差矩阵进行一些额外的计算（因为它返回简化的形式）。我用g()和反质量进行拟合，但是当引入f()并直接拟合m时，您可以立即重现您的问题。

第二点是数据具有偏移量，即，如果x = y = z = 0，则数据为v= -0.0195。这需要引入模型中。

最后，我想说您的数据中已经存在非抛物线行为。

不过，这是它的样子：

import matplotlib.pyplot as plt
import numpy as np
np.set_printoptions(linewidth=300)

from scipy.optimize import leastsq
from scipy.optimize import curve_fit

data = np.loadtxt( "silicon.csv", delimiter=',' )

def f( x, y, z, mxx, mxy, mxz, myy, myz, mzz, offI ):
    out = 1./(2 * mxx) * x * x
    out += 1./( mxy ) * x * y
    out += 1./( mxz ) * x * z
    out += 1./( 2 * myy ) * y * y
    out += 1./( myz ) * y * z
    out += 1./( 2 * mzz ) * z * z
    out += 1./offI
    return out

def g( x, y, z, mxxI, mxyI, mxzI, myyI, myzI, mzzI, off ):
    out = mxxI / 2 * x * x
    out += mxyI  * x * y
    out += mxzI * x * z
    out += myyI / 2 * y * y
    out += myzI  * y * z
    out += mzzI / 2  * z * z
    out += off
    return out


def residuals( params, indata ):
    out = list()
    for x, y, z, v in indata:
        out.append( v - g( x,y, z, *params ) )
    return out

sol, cov, info, msg, ier = leastsq( residuals,  7*[0], args=( data, ), full_output=True)

s_sq = sum( [x**2 for x in residuals( sol, data) ] )/ (len( data ) - len( sol ) )
print "solution"
print sol

masses = [1/x for x in sol]
print "masses:"
print masses

print "covariance matrix:"
covMX = cov * s_sq
print covMX

print "sum of residuals"
print sum( residuals( sol, data) )

### plotting the cuts

fig = plt.figure('cuts')
ax = dict()
for i in range( 1, 10 ):
    ax[i] = fig.add_subplot( 3, 3, i )

dl = np.linspace( -.2, .2, 25)

#### xx
xdata = [ [ x, v ] for x,y,z,v in data if ( abs(y)<1e-3 and abs(z) < 1e-3 ) ]
vl = np.fromiter( ( f( x, 0, 0, *masses ) for x in dl ), np.float )
ax[1].plot( *zip(*sorted( xdata ) ), ls='', marker='o')
ax[1].plot( dl, vl )

#### xy
xydata = [ [ x, v ] for x, y, z, v in data if ( abs( x - y )<1e-2 and abs(z) < 1e-3 ) ]
vl = np.fromiter( ( f( xy, xy, 0, *masses ) for xy in dl ), np.float )
ax[2].plot( *zip(*sorted( xydata ) ), ls='', marker='o')
ax[2].plot( dl, vl )

#### xz
xzdata = [ [ x, v ] for x, y, z, v in data if ( abs( x - z )<1e-2 and abs(y) < 1e-3 ) ]
vl = np.fromiter( ( f( xz, 0, xz, *masses ) for xz in dl ), np.float )
ax[3].plot( *zip(*sorted( xzdata ) ), ls='', marker='o')
ax[3].plot( dl, vl )

#### yy
ydata = [ [ y, v ] for x, y, z, v in data if ( abs(x)<1e-3 and abs(z) < 1e-3 ) ]
vl = np.fromiter( ( f( 0, y, 0, *masses ) for y in dl ), np.float )
ax[5].plot( *zip(*sorted( ydata ) ), ls='', marker='o' )
ax[5].plot( dl, vl )

#### yz
yzdata = [ [ y, v ] for x, y, z, v in data if ( abs( y - z )<1e-2 and abs(x) < 1e-3 ) ]
vl = np.fromiter( ( f( 0, yz, yz, *masses ) for yz in dl ), np.float )
ax[6].plot( *zip(*sorted( yzdata ) ), ls='', marker='o')
ax[6].plot( dl, vl )

#### zz
zdata = [ [ z, v ] for x, y, z, v in data if ( abs(x)<1e-3 and abs(y) < 1e-3 ) ]
vl = np.fromiter( ( f( 0, 0, z, *masses ) for z in dl ), np.float )
ax[9].plot( *zip(*sorted( zdata ) ), ls='', marker='o' )
ax[9].plot( dl, vl )

#### some diag
ddata = [ [ z, v ] for x, y, z, v in data if ( abs(x - y)<1e-3 and abs(x - z) < 1e-3 ) ]
vl = np.fromiter( ( f( d, d, d, *masses ) for d in dl ), np.float )
ax[7].plot( *zip(*sorted( ddata ) ), ls='', marker='o' )
ax[7].plot( dl, vl )

#### some other diag
ddata = [ [ z, v ] for x, y, z, v in data if ( abs(x - y)<1e-3 and abs(x + z) < 1e-3 ) ]
vl = np.fromiter( ( f( d, d, -d, *masses ) for d in dl ), np.float )
ax[8].plot( *zip(*sorted( ddata ) ), ls='', marker='o' )
ax[8].plot( dl, vl )

plt.show()

这将提供以下输出：

solution
[-1.46528595  0.25090717  0.25090717 -1.46528595  0.25090717 -1.46528595 -0.01993436]
masses:
[-0.6824606499739905, 3.985537743156507, 3.9855376943660676, -0.6824606473928339, 3.9855377322848344, -0.6824606467055248, -50.16463861555409]
covariance matrix:
[
[ 4.76417852e-03 -1.46907683e-12 -8.57639600e-12 -2.21281938e-12 -2.38444957e-12  8.42981521e-12 -2.70034183e-05]
[-1.46907683e-12  9.17104397e-04 -7.10573582e-13  1.32125214e-11  7.44553140e-12  1.29909935e-11 -1.11259046e-13]
[-8.57639600e-12 -7.10573582e-13  9.17104389e-04 -8.60004172e-12 -6.14797647e-12  8.27070243e-12  3.11127064e-14]
[-2.21281914e-12  1.32125214e-11 -8.60004172e-12  4.76417860e-03 -4.20477032e-12  9.20893224e-12 -2.70034186e-05]
[-2.38444957e-12  7.44553140e-12 -6.14797647e-12 -4.20477032e-12  9.17104395e-04  1.50963408e-11 -7.28889534e-14]
[ 8.42981530e-12  1.29909935e-11  8.27070243e-12  9.20893175e-12  1.50963408e-11  4.76417849e-03 -2.70034182e-05]
[-2.70034183e-05 -1.11259046e-13  3.11127064e-14 -2.70034186e-05 -7.28889534e-14 -2.70034182e-05  5.77019926e-07]
]
sum of residuals
4.352727352163743e-09

...在这里，如果一个不在主轴之一上，则某些一维切割显示出与抛物线行为的明显偏离。