SyntaxBomb - Indie Coders

Has anyone written a small program in SmallBasic I could use as a subroutine for fitting a 1st/2nd/3rd degree polynomial of around 10 points?  For example it would output the coefficients a, b &c in ax2+bx+c  It would also be useful if it could calculate the correlation coefficient.

Interesting challenge. I don't have such a program but it seems pretty easy to set up specially if we try not to solve everything at once like for any degree poly but say to make our task easy limit to Quadratic curves because I remember them from HS 50+ years ago!

So, instantly, then I find this on Internet:
https://www.youtube.com/watch?v=jLzkaJk0iZ0

and this:
oops lost link, no matter

So first question about our data, are the x values equally distant from each other? (again looking for easier case to solve)

Hi Stephen, a while ago I wrote a little program for curve fitting. Find below two source codes. The first one is for curve fitting. See the function Fit_Poly3 for some comments. The second source code is a little unit for plotting curves. You only need this for visualization. The whole thing is not  good commented and maybe not straight forward to understand. I hope it is still helpful.

import graphunit as g

Color 15,0
cls

MaxIterations = 30

x1 = seq(-8,8,11)
dim y1(len(x1) - 1)
dim y3(len(x1) - 1)
x2 = seq(-8, 8, 101)
dim y2(len(x2) - 1)
dim c

'Fit_Linear()
'Fit_Exp1()
'Fit_Exp1a()
'Fit_Gauss1()
'Fit_Gauss1a()
'Fit_Poly2()
Fit_Poly3()
'Fit_Sin1()
'Fit_Sin1a()
'Fit_Cos1()
'Fit_Cos1a()
'print y2   



Graph1 = g.GraphData()

'Draw measured data
Graph1.XLim = [min(x1), max(x1)]
Graph1.YLim = [min(y2), max(y2)] * 1.15
Graph1.Position = [50,50]
Graph1.Size = [600, 500]
'GraphOBJ.YLim = [-1.5, 1.5]
Graph1.DrawDots = true
Graph1.LineColor = 2
g.DrawGraph(Graph1, x1, y1)

'Draw fit function
Graph2 = Graph1
Graph2.DrawDots = false
Graph2.LineColor = 3
g.AddGraph(Graph2, x2, y2, [])

' Calc error
'for ii = 0 to len(x1) - 1
'        y3(ii) =  c(0) * cos(c(1) * x1(ii) + c(2))
'next

'ResSTD = ResidualSTD(y1, y3)

'print ResSTD
'print RMSE(y1, y3)

'yDelta = ArrayAdd(y2, ResSTD)
'g.AddGraph(Graph2, x2,yDelta,[])
'yDelta = ArrayAdd(y2, -ResSTD)
'g.AddGraph(Graph2, x2,yDelta,[])

showpage()

'################################################################################


sub CurveFit2(byref x, byref y, FunctionString, byref coeffs, MaxIterations)

   
    'Variables
    local DataLength = len(y) - 1 ' 0 is first element
    local NumCoeffs = len(coeffs) - 1
    local dim r(DataLength,0)
    local dim df_x1(DataLength,NumCoeffs)
    local dim df_x2(DataLength,NumCoeffs)
    local dim J(DataLength,NumCoeffs)
    local dim JT(NumCoeffs,DataLength)
    local r2_old = -1
    local dim r2
    local ii,jj, Temp_A, Temp_B, p

    for jj = 1 to MaxIterations       
       
        select(FunctionString)
            case "gauss1":        'a*exp(-(x-b)^2 /(2*c^2) )
                for ii = 0 to DataLength
                    'Derivative of gaussian function:
                    'see wolfram alpha
                    'i.e. https://www.wolframalpha.com/input/?i=d%2Fdc+a*exp%28-%28x-b%29%5E2%2F%282c%5E2%29%29+%2B+d
                    'Residual vector r
                    r(ii,0) = coeffs(0) * exp(-1*(x(ii) - coeffs(1))^2 / (2*coeffs(2)^2)) - y(ii)
                    'Jakobi-Matrix
                    J(ii,0) = exp(-1*(x(ii) - coeffs(1))^2 / (2*coeffs(2)^2))
                    J(ii,1) = coeffs(0) * (x(ii) - coeffs(1)) / coeffs(2)^2 * exp(-1*(x(ii) - coeffs(1))^2 / (2*coeffs(2)^2))
                    J(ii,2) = coeffs(0)/coeffs(2)^3 * (x(ii)-coeffs(1))^2 *  exp(-1*(x(ii) - coeffs(1))^2 / (2*coeffs(2)^2))
                next
            case "gauss1a":        'a*exp(-(b*x)^2/(2*c^2) ) + d
                for ii = 0 to DataLength
                    r(ii,0) = coeffs(0) * exp(-1*(x(ii) - coeffs(1))^2 / (2*coeffs(2)^2)) + coeffs(3) - y(ii)
                    'Jakobi-Matrix
                    J(ii,0) = exp(-1*(x(ii) - coeffs(1))^2 / (2*coeffs(2)^2))
                    J(ii,1) = coeffs(0) * (x(ii) - coeffs(1)) / coeffs(2)^2 * exp(-1*(x(ii) - coeffs(1))^2 / (2*coeffs(2)^2))
                    J(ii,2) = coeffs(0)/coeffs(2)^3 * (x(ii)-coeffs(1))^2 *  exp(-1*(x(ii) - coeffs(1))^2 / (2*coeffs(2)^2))
                    J(ii,3) = 1
                next
            case "exp1":    'a*exp(b*x)       
                for ii = 0 to DataLength
                    'Residual vector r
                    r(ii,0) = coeffs(0) * exp(coeffs(1) * x(ii)) - y(ii)
                    'Jakobi-Matrix
                    J(ii,0) = exp(coeffs(1) * x(ii))
                    J(ii,1) = coeffs(0) * x(ii) * exp(coeffs(1) * x(ii))
                next
            case "exp1a":    'a*exp(b*x)+c   
                for ii = 0 to DataLength
                    'Residual vector r
                    r(ii,0) = coeffs(0) * exp(coeffs(1) * x(ii)) + coeffs(2) - y(ii)
                    'Jakobi-Matrix
                    J(ii,0) = exp(coeffs(1) * x(ii))
                    J(ii,1) = coeffs(0) * x(ii) * exp(coeffs(1) * x(ii))
                    J(ii,2) = 1
                next
            case "lin":    'a*x + c
                for ii = 0 to DataLength
                    'Residual vector r
                    r(ii,0) = coeffs(0) * x(ii) + coeffs(1) - y(ii)
                    'Jakobi-Matrix
                    J(ii,0) = x(ii)
                    J(ii,1) = 1
                   
                next
            case "poly2":    'a*x^2 + b*x + c
                for ii = 0 to DataLength
                    'Residual vector r
                    r(ii,0) = coeffs(0) * x(ii)^2 + coeffs(1) * x(ii) + coeffs(2) - y(ii)
                    'Jakobi-Matrix
                    J(ii,0) = x(ii)^2
                    J(ii,1) = x(ii)
                    J(ii,2) = 1
                next
            case "poly3":    'a*x^3 + b*x^2 + c*x + d
                for ii = 0 to DataLength
                    'Residual vector r
                    r(ii,0) = coeffs(0) * x(ii)^3 + coeffs(1) * x(ii)^2 + coeffs(2) * x(ii) + coeffs(3) - y(ii)
                    'Jakobi-Matrix
                    J(ii,0) = x(ii)^3
                    J(ii,1) = x(ii)^2
                    J(ii,2) = x(ii)
                    J(ii,3) = 1
                next
            case "sin1":    'a*sin(b*x + c)
                for ii = 0 to DataLength
                    'Residual vector r
                    r(ii,0) = coeffs(0) * sin(coeffs(1) * x(ii) + coeffs(2)) - y(ii)
                    'Jakobi-Matrix
                    J(ii,0) = sin(coeffs(1) * x(ii) + coeffs(2))
                    J(ii,1) = coeffs(0)* x(ii) * cos(coeffs(1) * x(ii) + coeffs(2))
                    J(ii,2) = coeffs(0) * cos(coeffs(1) * x(ii) + coeffs(2))
                next
            case "sin1a":    'a*sin(b*x + c) + d
                for ii = 0 to DataLength
                    'Residual vector r
                    r(ii,0) = coeffs(0) * sin(coeffs(1) * x(ii) + coeffs(2)) + coeffs(3) - y(ii)
                    'Jakobi-Matrix
                    J(ii,0) = sin(coeffs(1) * x(ii) + coeffs(2))
                    J(ii,1) = coeffs(0)* x(ii) * cos(coeffs(1) * x(ii) + coeffs(2))
                    J(ii,2) = coeffs(0) * cos(coeffs(1) * x(ii) + coeffs(2))
                    J(ii,3) = 1
                next
            case "cos1":    'a*cos(b*x + c)
                for ii = 0 to DataLength
                    'Residual vector r
                    r(ii,0) = coeffs(0) * cos(coeffs(1) * x(ii) + coeffs(2)) - y(ii)
                    'Jakobi-Matrix
                    J(ii,0) = cos(coeffs(1) * x(ii) + coeffs(2))
                    J(ii,1) = coeffs(0)* x(ii) * -sin(coeffs(1) * x(ii) + coeffs(2))
                    J(ii,2) = coeffs(0) * -sin(coeffs(1) * x(ii) + coeffs(2))
                next
            case "cos1a":    'a*cos(b*x + c) + d
                for ii = 0 to DataLength
                    'Residual vector r
                    r(ii,0) = coeffs(0) * cos(coeffs(1) * x(ii) + coeffs(2)) + coeffs(3) - y(ii)
                    'Jakobi-Matrix
                    J(ii,0) = cos(coeffs(1) * x(ii) + coeffs(2))
                    J(ii,1) = coeffs(0)* x(ii) * -sin(coeffs(1) * x(ii) + coeffs(2))
                    J(ii,2) = coeffs(0) * -sin(coeffs(1) * x(ii) + coeffs(2))
                    J(ii,3) = 1
                next
        end select
       
   

        '||r||^2 as goodness of fit; smaller is better.
        r2 = 0
        for ii = 0 to DataLength
            r2 = r2 + r(ii,0) * r(ii,0)
        next
       
        'If value of r2 dosen't get smaller anymore, then exit for loop
        if(abs(r2 - r2_old) < 0.1 ) then
                exit sub
        endif
       
        r2_old = r2
       
   
        'Transpose Jakobi matrix
        for ii = 0 to DataLength
            for jj = 0 to NumCoeffs
                JT(jj,ii) = J(ii,jj)
            next               
        next
       
       


        'Solve p: J'J p = -J'r -> Ap = B
        Temp_A = JT * J
        Temp_B = -1 * (JT*r)
        p = LinEqn(Temp_A, Temp_B)

        'New coefficience
        for ii = 0 to NumCoeffs
            coeffs(ii) = coeffs(ii) + p(ii,0)
        next
       
    next

end sub

' Residual standard deviation
func ResidualSTD(byref y, byref yfit)
    local diff, n

    diff = y - yfit
    n = len(y)
   
    return sqr(sumsq(diff) / (n - 2))
end

' root mean squared error
func RMSE(byref y, byref yfit)
    local diff, n

    diff = y - yfit
    n = len(y)
   
    return sqr(sumsq(diff) /n)
end





sub Fit_Exp1()    'a*exp(b*x)
    for ii = 0 to len(x1) - 1
        y1(ii) =  0.7 * exp(0.5 * x1(ii))
    next
    c = [0.5,0.8]
    CurveFit2(x1, y1, "exp1", c, MaxIterations)
   
    for ii = 0 to len(x2) - 1
        y2(ii) = c(0) * exp(c(1)*(x2(ii)))
    next
   
end sub

sub Fit_Exp1a()    'a*exp(b*x)+c   
    for ii = 0 to len(x1) - 1
        y1(ii) =  0.7 * exp(0.5 * x1(ii)) + 1.1
    next
    c = [0.5,0.8,1.6]
    CurveFit2(x1, y1, "exp1a", c, MaxIterations)
   
    for ii = 0 to len(x2) - 1
        y2(ii) = c(0) * exp(c(1)*(x2(ii))) + c(2)
    next
   
end sub

sub Fit_Gauss1()
    for ii = 0 to len(x1) - 1
        y1(ii) =  0.7 * exp(-1*(x1(ii)-1.2)^2/(2*1.8^2))
    next
    c = [0.5,1.1,1.6]
    CurveFit2(x1, y1, "gauss1", c, MaxIterations)
   
    for ii = 0 to len(x2) - 1
        y2(ii) = c(0) * exp(-1*(x2(ii)-c(1))^2/(2*c(2)^2))
    next
end sub

sub Fit_Gauss1a()
    for ii = 0 to len(x1) - 1
        y1(ii) =  0.7 * exp(-1*(x1(ii)-1.2)^2/(2*1.8^2)) + 0.5
    next
    c = [0.5,1.1,1.6, 0]
    CurveFit2(x1, y1, "gauss1a", c, MaxIterations)
   
    for ii = 0 to len(x2) - 1
        y2(ii) = c(0) * exp(-1*(x2(ii)-c(1))^2/(2*c(2)^2)) + c(3)
    next
end sub

sub Fit_Linear()
    for ii = 0 to len(x1) - 1
        y1(ii) =  0.7 * x1(ii) + 1.2
    next
    c = [0.5, 1.1]
    CurveFit2(x1, y1, "lin", c, MaxIterations)
   
    for ii = 0 to len(x2) - 1
        y2(ii) = c(0) * x2(ii) + c(1)
    next
end sub

sub Fit_Poly2()
       
    for ii = 0 to len(x1) - 1
        y1(ii) =  0.1 * x1(ii)^2 + 0.2 * x1(ii) - 0.8
    next
    c = [0.5, 1.1, 0.3]
   
    CurveFit2(x1, y1, "poly2", c, MaxIterations)

    for ii = 0 to len(x2) - 1
        y2(ii) = c(0) * x2(ii)^2 + c(1) * x2(ii) + c(2)
    next
end sub

sub Fit_Poly3()
   
    ' Create data points; this could be for example a measurment   
    for ii = 0 to len(x1) - 1
        y1(ii) =  0.1 * x1(ii)^3 + 0.2 * x1(ii)^2 + 0.1 * x1(ii) - 0.8
    next
   
    ' Define starting values for fit, should be not too far away
    c = [0.5, 1.1, 0.3, 0.5]
   
    ' Do the fit
    CurveFit2(x1, y1, "poly3", c, MaxIterations)

    ' Create a curve using the fit results for plotting
    for ii = 0 to len(x2) - 1
        y2(ii) =  c(0) * x2(ii)^3 + c(1) * x2(ii)^2 + c(2) * x2(ii) + c(3)
    next
end sub

sub Fit_Sin1()
       
    for ii = 0 to len(x1) - 1
        y1(ii) =  0.7 * sin(0.5*x1(ii) + 0.7)
    next
    c = [0.4, 0.5, 0.1]

    CurveFit2(x1, y1, "sin1", c, MaxIterations)

    for ii = 0 to len(x2) - 1
        y2(ii) =  c(0) * sin(c(1) * x2(ii) + c(2))
    next
end sub

sub Fit_Sin1a()
       
    for ii = 0 to len(x1) - 1
        y1(ii) =  0.7 * sin(0.5*x1(ii) + 0.7) + 0.2
    next
    c = [0.4, 0.5, 0.1, 0.1]

    CurveFit2(x1, y1, "sin1a", c, MaxIterations)

    for ii = 0 to len(x2) - 1
        y2(ii) =  c(0) * sin(c(1) * x2(ii) + c(2)) + c(3)
    next
end sub

sub Fit_Cos1()
       
    for ii = 0 to len(x1) - 1
        y1(ii) =  (0.7+rnd/5 )* cos((0.5+rand/5)*x1(ii) + (0.7+rnd/5))
    next
    c = [0.4, 0.5, 0.1]

    CurveFit2(x1, y1, "cos1", c, MaxIterations)

    for ii = 0 to len(x2) - 1
        y2(ii) =  c(0) * cos(c(1) * x2(ii) + c(2))
    next
end sub

sub Fit_Cos1a()
       
    for ii = 0 to len(x1) - 1
        y1(ii) =  0.7 * cos(0.5*x1(ii) + 0.7) + 0.2
    next
    c = [0.4, 0.5, 0.1, 0.1]

    CurveFit2(x1, y1, "cos1a", c, MaxIterations)

    for ii = 0 to len(x2) - 1
        y2(ii) =  c(0) * cos(c(1) * x2(ii) + c(2)) + c(3)
    next
end sub    

func ArrayAdd(a, c)

    local l, b
    l = ubound(a)
    dim b(l)
       
    for ii = 0 to l
        b[ii] =  a[ii] + c
    next
   
    return b

end



Unit graphunit

export GraphData
export DrawGraph
export AddGraph

func GraphData()

    local TempG
   
    TempG.Position = [0,0]
    TempG.Size = [xmax,ymax]
    TempG.XLim = []
    TempG.YLim = []
    TempG.DrawAxis = true
    TempG.DrawDots = false
    TempG.LineColor = 1
    TempG.AxisColor = rgb(255,255,255)
    TempG.DrawLabels = true
    TempG.DrawTicks = true
    TempG.TickSeparationX = 80
    TempG.TickSeparationY = 80
   
    TempG.InternalAddGraph = false
    TempG.InternalSizePlotArea = [0,0,0,0]
      
    GraphData = TempG
end


sub DrawGraph(Byref G, Byref XData , Byref YData)

    local Temp
    local EndIndex = (len(XData)) - 1
    local dim Poly
    local MaxValueX = max(XData)
    local MaxValueY = max(YData)
    local ii
    local SizePlotArea
    local XLim, YLim
  
    XLim = G.XLim
    YLim = G.YLim
   
    if(empty(XLim)) then
        XLim = [XData(0), XData(len(XData)-1)]
    endif
   
    if(empty(YLim)) then
        YLim = [min(YData), max(YData)]
    endif
   
    'Color G.AxisColor
   
    if(!G.InternalAddGraph)
        if(G.DrawLabels) then
            OffsetY = textheight("Tq") * 1.1
            OffsetX = textwidth("-1.23") * 1.1
            SizePlotArea = [G.Size(0) - OffsetX, G.Size(1) - OffsetY]
        else
            OffsetY = 0
            OffsetX = 0
            SizePlotArea = [G.Size(0) - OffsetX, G.Size(1) - OffsetY]
        endif
        G.InternalSizePlotArea = SizePlotArea
    else
        SizePlotArea = G.InternalSizePlotArea      
    endif
   
    'Calculate points
    for ii = 0 to EndIndex
       
        if(XData(ii) >= XLim(0) AND XData(ii) <= XLim(1)) then
            Poly << SizePlotArea(0) / ( XLim(1) - XLim(0) ) * (XData(ii) - XLim(0))
           
            if(YData(ii) < YLim(0)) then
                Poly << SizePlotArea(1)
            elseif(YData(ii) > YLim(1)) then
                Poly << 0
            else
                Poly << SizePlotArea(1) - SizePlotArea(1) / ( YLim(1) - YLim(0) ) * (YData(ii) - YLim(0))
            endif
           
        endif
   
    next

    if(G.DrawDots) then
       
        for ii = 0 to len(Poly) - 1 Step 2
            circle G.Position(0) + Poly(ii) + OffsetX, G.Position(1) + Poly(ii+1), 3, 1, G.LineColor
        next
       
    else
        DrawPoly Poly(), G.Position(0) + OffsetX, G.Position(1), 1, G.LineColor
    endif
   
   
    if(G.DrawAxis) then
        rect G.Position(0) + OffsetX, G.Position(1), G.Position(0)+ SizePlotArea(0) + OffsetX, G.Position(1) + SizePlotArea(1) , G.AxisColor
           
       
        if(G.DrawLabels or G.DrawTicks) then
           
            'NumberOfTicksX = floor(SizePlotArea(0) / G.TickSeparationX)
            'NumberOfTicksY = floor(SizePlotArea(1) / G.TickSeparationY)
            'G.TickSeparationX = SizePlotArea(0) / NumberOfTicksX
            'G.TickSeparationY = SizePlotArea(1) / NumberOfTicksY
                                   
            NumberOfTicksX =(SizePlotArea(0) / G.TickSeparationX)
            NumberOfTicksY =(SizePlotArea(1) / G.TickSeparationY)
                               
            'X-Labels
            posy = G.Position(1) + SizePlotArea(1)
            for ii = 0 to NumberOfTicksX * G.TickSeparationX step G.TickSeparationX
                posx = G.Position(0) + ii + OffsetX
                if(G.DrawTicks) then
                    line posx, posy, posx, posy-5, G.AxisColor
                endif
                if(G.DrawLabels) then
                    at posx,posy+4: print  round((XData(EndIndex) - XData(0))/NumberOfTicksX * ii/G.TickSeparationX + XData(0), 2)
                endif
            next
           
            'Y-Labels
            posx = G.Position(0) + OffsetX
            for ii = 0 to NumberOfTicksY * G.TickSeparationY step G.TickSeparationY
                posy = G.Position(1) + SizePlotArea(1) - ii
                if(G.DrawTicks) then
                    line posx, posy, posx+5, posy, G.AxisColor
                endif
                if(G.DrawLabels) then
                    at posx - OffsetX,posy - OffsetY/3: print  round((max(YData) - min(YData))/NumberOfTicksY * ii/G.TickSeparationY + min(YData), 2)
                                       
                endif
            next
        endif
    endif
   
end

sub AddGraph(Byref G, Byref XData, Byref YData, LineColor)
   
    local TempG
   
    if(LineColor == [])
        G.LineColor++
        if(G.LineColor > 15) then G.LineColor = 1
    endif
   
    TempG = G
    TempG.DrawAxis = false
    TempG.DrawTicks = false
    TempG.InternalAddGraph = true
   
    DrawGraph(TempG, XData, YData)
end

Wow J7M looks very experienced even professional. 

For me though, I never did import export with sb. I did try units and I think you need to include those files with your code if you want someone else to use that code (as well as for importing).

Me I want something simpler, plug in n amount of points, out dumps my coefficients for Quad with some sort or reliability factor. 

as I wrote, you don't need the plotting unit it is only for visualization. Here a cleaner version. Only the fitting part.  As an example the fit for a polygon of 3. power.

' Fit polygon of 3. power: y = ax^3 + bx^2 + cx + d

' Create data points; this could be for example a measurement
' as an example a poly3 with the following coefficients is used
' a = 0.1, b = 0.2, c = 0.1, d = -0.8

x = [-5,-4,-3,-2,-1,0,1,2,3,4,5]
y = [-8.8,-4.4,-2,-1,-0.8,-0.8,-0.4,1,4,9.2,17.2]

' Define starting values for fit, should be not too far away from the expected coefficients
c = [0.5, 1.1, 0.3, 0.5]

' Do the fit, c will hold the fitted coefficients
CurveFit(x, y, "poly3", c, 30)

print "Fitted coefficienc are: "; c


'################################################################################


sub CurveFit(byref x, byref y, FunctionString, byref coeffs, MaxIterations)

'Variables
local DataLength = len(y) - 1 ' 0 is first element
local NumCoeffs = len(coeffs) - 1
local dim r(DataLength,0)
local dim df_x1(DataLength,NumCoeffs)
local dim df_x2(DataLength,NumCoeffs)
local dim J(DataLength,NumCoeffs)
local dim JT(NumCoeffs,DataLength)
local r2_old = -1
local dim r2
local ii,jj, Temp_A, Temp_B, p

for jj = 1 to MaxIterations

select(FunctionString)
case "gauss1": 'a*exp(-(x-b)^2 /(2*c^2) )
for ii = 0 to DataLength
'Derivative of gaussian function:
'see wolfram alpha
'i.e. https://www.wolframalpha.com/input/?i=d%2Fdc+a*exp%28-%28x-b%29%5E2%2F%282c%5E2%29%29+%2B+d
'Residual vector r
r(ii,0) = coeffs(0) * exp(-1*(x(ii) - coeffs(1))^2 / (2*coeffs(2)^2)) - y(ii)
'Jakobi-Matrix
J(ii,0) = exp(-1*(x(ii) - coeffs(1))^2 / (2*coeffs(2)^2))
J(ii,1) = coeffs(0) * (x(ii) - coeffs(1)) / coeffs(2)^2 * exp(-1*(x(ii) - coeffs(1))^2 / (2*coeffs(2)^2))
J(ii,2) = coeffs(0)/coeffs(2)^3 * (x(ii)-coeffs(1))^2 *  exp(-1*(x(ii) - coeffs(1))^2 / (2*coeffs(2)^2))
next
case "gauss1a": 'a*exp(-(b*x)^2/(2*c^2) ) + d
for ii = 0 to DataLength
r(ii,0) = coeffs(0) * exp(-1*(x(ii) - coeffs(1))^2 / (2*coeffs(2)^2)) + coeffs(3) - y(ii)
'Jakobi-Matrix
J(ii,0) = exp(-1*(x(ii) - coeffs(1))^2 / (2*coeffs(2)^2))
J(ii,1) = coeffs(0) * (x(ii) - coeffs(1)) / coeffs(2)^2 * exp(-1*(x(ii) - coeffs(1))^2 / (2*coeffs(2)^2))
J(ii,2) = coeffs(0)/coeffs(2)^3 * (x(ii)-coeffs(1))^2 *  exp(-1*(x(ii) - coeffs(1))^2 / (2*coeffs(2)^2))
J(ii,3) = 1
next
case "exp1": 'a*exp(b*x)
for ii = 0 to DataLength
'Residual vector r
r(ii,0) = coeffs(0) * exp(coeffs(1) * x(ii)) - y(ii)
'Jakobi-Matrix
J(ii,0) = exp(coeffs(1) * x(ii))
J(ii,1) = coeffs(0) * x(ii) * exp(coeffs(1) * x(ii))
next
case "exp1a": 'a*exp(b*x)+c
for ii = 0 to DataLength
'Residual vector r
r(ii,0) = coeffs(0) * exp(coeffs(1) * x(ii)) + coeffs(2) - y(ii)
'Jakobi-Matrix
J(ii,0) = exp(coeffs(1) * x(ii))
J(ii,1) = coeffs(0) * x(ii) * exp(coeffs(1) * x(ii))
J(ii,2) = 1
next
case "lin": 'a*x + c
for ii = 0 to DataLength
'Residual vector r
r(ii,0) = coeffs(0) * x(ii) + coeffs(1) - y(ii)
'Jakobi-Matrix
J(ii,0) = x(ii)
J(ii,1) = 1

next
case "poly2": 'a*x^2 + b*x + c
for ii = 0 to DataLength
'Residual vector r
r(ii,0) = coeffs(0) * x(ii)^2 + coeffs(1) * x(ii) + coeffs(2) - y(ii)
'Jakobi-Matrix
J(ii,0) = x(ii)^2
J(ii,1) = x(ii)
J(ii,2) = 1
next
case "poly3": 'a*x^3 + b*x^2 + c*x + d
for ii = 0 to DataLength
'Residual vector r
r(ii,0) = coeffs(0) * x(ii)^3 + coeffs(1) * x(ii)^2 + coeffs(2) * x(ii) + coeffs(3) - y(ii)
'Jakobi-Matrix
J(ii,0) = x(ii)^3
J(ii,1) = x(ii)^2
J(ii,2) = x(ii)
J(ii,3) = 1
next
case "sin1": 'a*sin(b*x + c)
for ii = 0 to DataLength
'Residual vector r
r(ii,0) = coeffs(0) * sin(coeffs(1) * x(ii) + coeffs(2)) - y(ii)
'Jakobi-Matrix
J(ii,0) = sin(coeffs(1) * x(ii) + coeffs(2))
J(ii,1) = coeffs(0)* x(ii) * cos(coeffs(1) * x(ii) + coeffs(2))
J(ii,2) = coeffs(0) * cos(coeffs(1) * x(ii) + coeffs(2))
next
            case "sin1a": 'a*sin(b*x + c) + d
for ii = 0 to DataLength
'Residual vector r
r(ii,0) = coeffs(0) * sin(coeffs(1) * x(ii) + coeffs(2)) + coeffs(3) - y(ii)
'Jakobi-Matrix
J(ii,0) = sin(coeffs(1) * x(ii) + coeffs(2))
J(ii,1) = coeffs(0)* x(ii) * cos(coeffs(1) * x(ii) + coeffs(2))
J(ii,2) = coeffs(0) * cos(coeffs(1) * x(ii) + coeffs(2))
                    J(ii,3) = 1
next
            case "cos1": 'a*cos(b*x + c)
for ii = 0 to DataLength
'Residual vector r
r(ii,0) = coeffs(0) * cos(coeffs(1) * x(ii) + coeffs(2)) - y(ii)
'Jakobi-Matrix
J(ii,0) = cos(coeffs(1) * x(ii) + coeffs(2))
J(ii,1) = coeffs(0)* x(ii) * -sin(coeffs(1) * x(ii) + coeffs(2))
J(ii,2) = coeffs(0) * -sin(coeffs(1) * x(ii) + coeffs(2))
next
            case "cos1a": 'a*cos(b*x + c) + d
for ii = 0 to DataLength
'Residual vector r
r(ii,0) = coeffs(0) * cos(coeffs(1) * x(ii) + coeffs(2)) + coeffs(3) - y(ii)
'Jakobi-Matrix
J(ii,0) = cos(coeffs(1) * x(ii) + coeffs(2))
J(ii,1) = coeffs(0)* x(ii) * -sin(coeffs(1) * x(ii) + coeffs(2))
J(ii,2) = coeffs(0) * -sin(coeffs(1) * x(ii) + coeffs(2))
                    J(ii,3) = 1
next
end select


'||r||^2 as goodness of fit; smaller is better.
r2 = 0
for ii = 0 to DataLength
r2 = r2 + r(ii,0) * r(ii,0)
next

'If value of r2 dosen't get smaller anymore, then exit for loop
if(abs(r2 - r2_old) < 0.1 ) then
exit sub
endif

r2_old = r2


'Transpose Jakobi matrix
for ii = 0 to DataLength
for jj = 0 to NumCoeffs
JT(jj,ii) = J(ii,jj)
next
next

'Solve p: J'J p = -J'r -> Ap = B
Temp_A = JT * J
Temp_B = -1 * (JT*r)
p = LinEqn(Temp_A, Temp_B)

'New coefficience
for ii = 0 to NumCoeffs
coeffs(ii) = coeffs(ii) + p(ii,0)
next

next

end sub

' Residual standard deviation
func ResidualSTD(byref y, byref yfit)
    local diff, n

    diff = y - yfit
    n = len(y)
   
    return sqr(sumsq(diff) / (n - 2))
end

' root mean squared error
func RMSE(byref y, byref yfit)
    local diff, n

    diff = y - yfit
    n = len(y)
   
    return sqr(sumsq(diff) /n)
end

Is Stephen still with us?

Me, I'd still want to plug in n amount points (actually very few needed for quadratic) be told the points do not suggest a Quadratic curve or be told the coefficients to Quadratic with some reliability factor.

I think we need some Guassian elimination code that I think I did see in Library examples from way back.
Might need to sort on x's first just to be sure they are in order, then do difference calcs twice then ready for Guass to solve for 3 unknowns for standard form of Quadratic if the 2nd differences are constant or close enough (if not then not Quadratic goodbye!).

Ah my memory is still serving me, from sb Library of old code:
[pre][color=#000000]
'''gaussj.bas
''version 1.0
''by adolfo leon sepulveda
'' 20/08/2004
''solve a nxn equation system with gauss jordan method
''It''s equivatent to LinEqn function of SmallBASIC
''example:
'' 0x+1y=1
'' 1x+1y=2
''
''Input:
'' Nro Equat: 2
''  a(1,1) : ? 0
''  a(1,2) : ? 1
''  b(1) : ? 1
''  a(2,1) : ? 1
''  a(2,1) : ? 1
''  b(2) : ? 2
'' result:
'' x=1
'' y=1
'' invert matrix:
'' -1 1
''  1 0
n=0
dim a(50,50)
dim x(50)
main
End ''program
sub main()
 readata a,x,n
 gaussj a,x,n
 Printinversa a,n
 Printsolution x,n
end
sub  readata( byref a, byref x, byref n)
print "gauss jordan"
print "nro equat:  ";
input n
print "coef y term indep "
for i=0 to n-1
  for j=0 to n
    if j<>n then
      print "a(";i+1;", ";j+1;") : ";
    else
        print "b(";i+1;") : ";
    endif
    input a(i,j)
  next
  x(i)=a(i,n)
next
end
sub gaussj(byref a, byref b,n)
local mayor,tmp,pivinv
dim indco(50),indfi(50),piv(50)
local i,j,k,t,h,col,fil
for j = 0 to n-1
  piv(j)=0
next j
for i = 0 to n-1
  mayor = 0.0
  for j = 0 to n-1
    if piv(j) <> 1 then
      for k = 0 to n-1
          if piv(k) = 0 then
          if abs(a(j,k)) >= mayor then
            mayor = abs(a(j,k))
            fil = j
            col = k
          endif
          elseif piv(k) > 1 then
            print "error: singular matrix"
            Exit Sub
          endif
      next
      endif
  next
  piv(col) = piv(col)+1
  if fil <> col then
    for t = 0 to n-1
      tmp = a(fil,t)
        a(fil,t)  = a(col,t)
        a(col,t) = tmp
    next
    tmp = b(fil)
    b(fil) = b(col)
    b(col) = tmp
  endif
  indfi(i) = fil
  indco(i) = col
  if a(col,col) = 0 then
    print "error: singular matrix"
    Exit Sub
  endif
  pivinv= 1.0/a(col,col)
  a(col,col) = 1.0
  for t=0 to n-1
    a(col,t) = a(col,t)*pivinv
  next
  b(col) = b(col)*pivinv
  for h = 0 to n-1
  if h <> col then
    tmp = a(h,col)
      a(h,col) = 0.0
      for t = 0 to n-1
        a(h,t) = a(h,t) - a(col,t)*tmp
      next
      b(h) = b(h) - b(col)*tmp
  endif
  next
next
for t = n-1 to 0
  if indfi(t) <> indco(t) then
    for k = 0 to n - 1
      tmp = a(k,indfi(t)) = a(k,indco(t))
      a(k,indco(t)) = tmp
    next
  endif
next
end
sub Printsolution(x,n)
 print "solucion"
for j=0 to n-1
  print "x";j+1;"=";x(j)
next
end
sub Printinversa(a,n)
print "matriz inversa"
for i=0 to n-1
  for j=0 to n-1
    print using "##0.00"; a(i,j);
  next
  print
next
end[/color][/pre]


 PS these last 2 posts from me refer to method I linked above in my first reply to Stephen.
They walk you through to point where you have 3 equations and 3 unknowns from which you apply Guassj.bas after some mods of course setting up the matrix for Guass.