使用 MSE 将二次函数拟合到数据

Fit quadratic function to data using MSE

所以我的想法是(从神经网络人员那里借来的)如果我有数据集 D,我可以通过首先计算误差相对于参数 (a, b,和 c),然后做一个最小化错误的小更新。我的问题是,下面的代码实际上并不符合曲线。对于线性的东西,类似的方法可行,但二次方似乎由于某种原因失败了。你能看出我做错了什么吗(假设或只是实现错误)

编辑:问题不够具体:以下代码不能很好地处理数据偏差。出于某种原因,它会以某种方式更新 a 和 b 参数,从而使 c 落在后面。这种方法类似于机器人学(使用雅可比矩阵寻找路径)和神经网络(基于错误寻找参数)所以它不是不合理的算法,现在的问题是,为什么这个具体实现没有产生我期望的结果。

在下面的Python代码中,我使用math作为m,MSE是一个计算两个数组之间的均方误差的函数。除此之外,代码是独立的

代码:

def quadraticRegression(data, dErr):
    a = 1 #Starting values
    b = 1
    c = 1

    a_momentum = 0 #Momentum to counter steady state error
    b_momentum = 0
    c_momentum = 0

    estimate = [a*x**2 + b*x + c for x in range(len(data))] #Estimate curve
    error = MSE(data, estimate) #Get errors 'n stuff
    errorOld = 0

    lr = 0.0000000001 #learning rate
    while abs(error - errorOld) > dErr:
        #Fit a (dE/da)
        deda = sum([ 2*x**2 * (a*x**2 + b*x + c - data[x]) for x in range(len(data)) ])/len(data)
        correction = deda*lr
        a_momentum = (a_momentum)*0.99 + correction*0.1 #0.99 is to slow down momentum when correction speed changes
        a = a - correction - a_momentum
    
        #fit b (dE/db)
        dedb = sum([ 2*x*(a*x**2 + b*x + c - data[x]) for x in range(len(data))])/len(data)
        correction = dedb*lr
        b_momentum = (b_momentum)*0.99 + correction*0.1
        b = b - correction - b_momentum

        #fit c (dE/dc)
        dedc = sum([ 2*(a*x**2 + b*x + c - data[x]) for x in range(len(data))])/len(data)
        correction = dedc*lr
        c_momentum = (c_momentum)*0.99 + correction*0.1
        c = c - correction - c_momentum

        #Update model and find errors
        estimate = [a*x**2 +b*x + c for x in range(len(data))]
        errorOld = error
        print(error)
        error = MSE(data, estimate)

    return a, b, c, error

对我来说,您的代码看起来完全正确!至少算法是正确的。我已将您的代码更改为使用 numpy 而不是纯 Python 进行快速计算。我还配置了一些参数,比如改变动量和学习率,也实现了 MSE.

然后我用matplotlib画了剧情动画。最后在动画上看起来你的回归实际上试图使曲线适合数据。虽然在最后一次迭代中拟合 sin(x) for x in [0; 2 * pi] 看起来像线性近似,但仍然尽可能接近数据点以求二次曲线。但是对于 sin(x) for x in [0; pi] 它看起来像是理想的近似值(它从大约 12-th 迭代开始拟合)。

i-th帧动画只是用dErr = 0.7 ** (i + 15)做回归。

我的动画对于 运行 脚本来说有点慢,但是如果你像这样添加参数 save python script.py save 它会 render/save 到 line.gif情节绘图动画。如果您 运行 没有参数的脚本,它将在您的 PC 屏幕上实时显示动画 plotting/fitting。

完整代码在图形之后,代码需要通过 运行ning 安装一些 python 模块一次 python -m pip install numpy matplotlib

接下来是 sin(x) for x in (0, pi):

接下来是 sin(x) for x in (0, 2 * pi):

接下来是 abs(x) for x in (-1, 1):

# Needs: python -m pip install numpy matplotlib
import math, sys
import numpy as np, matplotlib.pyplot as plt, matplotlib.animation as animation
from matplotlib.animation import FuncAnimation


x_range = (0., math.pi, 0.1) # (xmin, xmax, xstep)
y_range = (-0.2, 1.2) # (ymin, ymax)
num_iterations = 50

def f(x):
    return np.sin(x)

def derr(iteration):
    return 0.7 ** (iteration + 15)
    
    
def MSE(a, b):
    return (np.abs(np.array(a) - np.array(b)) ** 2).mean()

def quadraticRegression(*, x, data, dErr):
    x, data = np.array(x), np.array(data)
    assert x.size == data.size, (x.size, data.size)

    a = 1 #Starting values
    b = 1
    c = 1

    a_momentum = 0.1 #Momentum to counter steady state error
    b_momentum = 0.1
    c_momentum = 0.1

    estimate = a*x**2 + b*x + c #Estimate curve
    error = MSE(data, estimate) #Get errors 'n stuff
    errorOld = 0.

    lr = 10. ** -4 #learning rate
    while abs(error - errorOld) > dErr:
        #Fit a (dE/da)
        deda = np.sum(2*x**2 * (a*x**2 + b*x + c - data))/len(data)
        correction = deda*lr
        a_momentum = (a_momentum)*0.99 + correction*0.1 #0.99 is to slow down momentum when correction speed changes
        a = a - correction - a_momentum
    
        #fit b (dE/db)
        dedb = np.sum(2*x*(a*x**2 + b*x + c - data))/len(data)
        correction = dedb*lr
        b_momentum = (b_momentum)*0.99 + correction*0.1
        b = b - correction - b_momentum

        #fit c (dE/dc)
        dedc = np.sum(2*(a*x**2 + b*x + c - data))/len(data)
        correction = dedc*lr
        c_momentum = (c_momentum)*0.99 + correction*0.1
        c = c - correction - c_momentum

        #Update model and find errors
        estimate = a*x**2 +b*x + c
        errorOld = error
        #print(error)
        error = MSE(data, estimate)

    return a, b, c, error


        
fig, ax = plt.subplots()
fig.set_tight_layout(True)

x = np.arange(x_range[0], x_range[1], x_range[2])
#ax.scatter(x, x + np.random.normal(0, 3.0, len(x)))
line0, line1 = None, None

do_save = len(sys.argv) > 1 and sys.argv[1] == 'save'

def g(x, derr):
    a, b, c, error = quadraticRegression(x = x, data = f(x), dErr = derr)
    return a * x ** 2 + b * x + c
    
def dummy(x):
    return np.ones_like(x, dtype = np.float64) * 100.

def update(i):
    global line0, line1
    
    de = derr(i)
    
    if line0 is None:
        assert line1 is None
        line0, = ax.plot(x, f(x), 'r-', linewidth=2)
        line1, = ax.plot(x, g(x, de), 'r-', linewidth=2, color = 'blue')
        ax.set_ylim(y_range[0], y_range[1])
        
    if do_save:
        sys.stdout.write(str(i) + ' ')
        sys.stdout.flush()

    label = 'iter {0} derr {1}'.format(i, round(de, math.ceil(-math.log(de) / math.log(10)) + 2))
    line1.set_ydata(g(x, de))
    ax.set_xlabel(label)
    return line1, ax

if __name__ == '__main__':
    anim = FuncAnimation(fig, update, frames = np.arange(0, num_iterations), interval = 200)
    if do_save:
        anim.save('line.gif', dpi = 200, writer = 'imagemagick')
    else:
        plt.show()