PureBasic Forum
https://www.purebasic.fr/english/

Cubic regression
https://www.purebasic.fr/english/viewtopic.php?f=7&t=73208
Page 1 of 1

Author:  zxretrosoft [ Wed Jul 17, 2019 6:47 am ]
Post subject:  Cubic regression

Hello friends,

I have a question. I need to program a higher regression than quadratic, namely, cubic regression.

To do this, we need to know how the coefficients A, B, C, D are calculated.

I found coefficients A, B for linear:
https://www.easycalculation.com/statistics/learn-regression.php

I found coefficients A, B, C for quadratic:
https://www.easycalculation.com/statistics/learn-quadratic-regression.php

But the coefficients A, B, C, D for cubic regression I can't find anywhere :(

Please, do you have someone tip? How are these coefficients calculated?

Thanks in advance! :idea:

Author:  STARGÅTE [ Wed Jul 17, 2019 10:24 am ]
Post subject:  Re: Cubic regression

The general case is written here:
Polynomial regression - Matrix form and calculation of estimates

I think you can find matrix calculation procedures here in the forum.

Author:  zxretrosoft [ Wed Jul 17, 2019 10:41 am ]
Post subject:  Re: Cubic regression

Yes, thank you.

But,
it's very complicated for me. I can't deduce it. I'm just looking for 4 equations for
A = ...
B = ...
C = ...
D = ...

Nothing more, nothing less.
Exactly what is here (A=... B=... C=...) for quadratic regression:
https://www.easycalculation.com/statistics/learn-quadratic-regression.php

Author:  STARGÅTE [ Wed Jul 17, 2019 11:02 am ]
Post subject:  Re: Cubic regression

Here is a updated version from "alter Mann":
https://www.purebasic.fr/german/viewtopic.php?f=8&t=18033&hilit=Polynomapproximation

You can use Num_ApproxPoly (n.l,m.l,Array x.d(1),Array y.d(1), Array w.d(1),Array c.d(1)) where:
n is the number of points
m is the degree of the polynom
x is the array of x-values
y is the array of y-values
w is an array for weighting (typical 1.0)
c is the resulting array of you coefficients (A,B,C,...)

Code:
; Author: alter Mann , adaptiert von Engeln-Müllges/Uhlig "Numerical Algorithms with C"
; Date: 23. 10. 2008
; OS: Windows
; Demo: Ja
;
; Polynomapproximation (Ausgleichpolynom)
; ( wird leicht instabil bei hohen Polynomgraden und großen Stützpunkt-Werten)


#DBL_EPSILON    =2.2204460492503131e-016     ; kleinster Wert, so dass 1.0+DBL_EPSILON != 1.0

#EPS_QUAD       =(#DBL_EPSILON*#DBL_EPSILON)

EnableExplicit

Procedure.d MathSqrtD (Wert.d)
Protected Wert1.d
If Wert <= 0.0
  Wert1 = 0.0
Else
  Wert1 = Sqr(Wert)
EndIf
ProcedureReturn Wert1
EndProcedure

Macro SQRT(wert)
  MathSqrtD(wert)
EndMacro

Macro SQR2(wert)
  (wert)*(wert)
EndMacro

;************************************************************************
;* ein Polynom P in der Darstellung                                     *
;* P(x)  =  adA[0] + adA[1] * dX + adA[2] * dX^2 + ... + adA[n] * dX^n  *
;* nach dem Hornerschema auswerten                                      *
;*                                                                      *
;* Eingabeparameter:                                                    *
;* =================                                                    *
;* lN: Grad des Polynoms                                                *
;* adA: [0..n]-Vektor mit den Koeffizienten des Polynoms                *
;* dX: Stelle, an der das Polynom auszuwerten ist                       *
;*                                                                      *
;* Funktionswert:                                                       *
;* ==============                                                       *
;* P(x)                                                                 *
;************************************************************************
Procedure.d Num_Horner (lN.l, Array adA.d(1), dX.d) ; Hornerschema zur Polynomauswertung .............*/

  Protected ldErg.d = adA(lN)
  Protected i.l

  For i=lN-1 To 0 Step -1
    ldErg = ldErg * dX + adA(i)
  Next i
  ProcedureReturn ldErg
EndProcedure

;*====================================================================*
;*                                                                    *
;*  Num_ChoCec berechnet die Zerlegungsmatrix einer symmetrischen,    *
;*  positiv definiten Matrix. Die Zerlegungsmatrix wird auf mat       *
;*  ueberspeichert.                                                   *
;*                                                                    *
;*====================================================================*
;*                                                                    *
;*   Eingabeparameter:                                                *
;*   ================                                                 *
;*      n        Long n  ( n > 0 )                                    *
;*               Dimension von mat,                                   *
;*               Anzahl der Komponenten des b-Vektors.                *
;*      mat      Double  mat(n,n)                                     *
;*               Matrix des Gleichungssystems.                        *
;*                                                                    *
;*   Ausgabeparameter:                                                *
;*   ================                                                 *
;*      mat      Double  mat(n,n)                                     *
;*               Dekompositionsmatrix, die die Zerlegung von          *
;*               mat enthaelt.                                        *
;*                                                                    *
;*   Rueckgabewert:                                                   *
;*   =============                                                    *
;*      = 0      alles ok                                             *
;*      = 1      n < 1 gewaehlt                                       *
;*      = 2      Matrix ist nicht positiv definit                     *
;*                                                                    *
;*====================================================================*
Procedure.l Num_ChoDec (n.l, Array mat.d(2))

  Protected j.l, k.l, i.l
  Protected sum.d

  If n < 1
    ProcedureReturn 1                         ; n < 1  Fehler !
  EndIf

  If mat(0,0) < #EPS_QUAD
    ProcedureReturn 2                         ; mat ist nicht positiv definit !
  EndIf

  mat(0,0) = SQRT(mat(0,0))
  For j = 1 To n-1 Step 1
    mat(j,0) / mat(0,0)
  Next j

  For i = 1 To n-1 Step 1
    sum = mat(i,i)
    For j = 0 To i-1 Step 1
      sum - SQR2(mat(i,j))
    Next j
    If sum < #EPS_QUAD
      ProcedureReturn 2                       ; nicht positiv definit
    EndIf
    mat(i,i) = SQRT(sum)
    For j = i + 1 To n-1 Step 1
      sum = mat(j,i)
      For k = 0 To i-1 Step 1
        sum - mat(i,k) * mat(j,k)
      Next k
      mat(j,i) = sum / mat(i,i)
    Next j
  Next i

  ProcedureReturn 0
EndProcedure
;*====================================================================*
;*                                                                    *
;*  Num_ChoSol bestimmt die Loesung x des linearen Gleichungssystems  *
;*  B * x = b mit der unteren Dreieckmatrix B wie sie als Ausgabe     *
;*  von chodec bestimmt wird.                                         *
;*                                                                    *
;*====================================================================*
;*                                                                    *
;*   Eingabeparameter:                                                *
;*   ================                                                 *
;*      n        Long n  ( n > 0 )                                    *
;*               Dimension von lmat, Anzahl der Komponenten           *
;*               des b-Vektors u. des Loesungsvektors x.              *
;*      lmat     Double   lmat(n,n)                                   *
;*               untere Dreieckmatrix, wie sie von chodec als Aus-    *
;*               gabe geliefert wird.                                 *
;*      b        Double   b(n)                                        *
;*               Rechte Seite des Gleichungssystems.                  *
;*                                                                    *
;*   Ausgabeparameter:                                                *
;*   ================                                                 *
;*      x        Double   x(n)                                        *
;*               Loesungsvektor des Systems.                          *
;*                                                                    *
;*   Rueckgabewert:                                                   *
;*   =============                                                    *
;*      = 0      alles ok                                             *
;*      = 1      Zerlegungsmatrix unzulaessig oder n < 1              *
;*                                                                    *
;*====================================================================*
Procedure.l Num_ChoSol (n.l, Array lmat.d(2), Array b.d(1), Array x.d(1))

  Protected j.l, k.l
  Protected sum.d

  If n < 1
    ProcedureReturn 1
  EndIf

  If lmat(0,0) = 0.0
    ProcedureReturn 1               ; Unzulaessige Zerlegungsmatrix
  EndIf

  x(0) = b(0) / lmat(0,0)           ; Vorwaertselimination
  For k = 1 To n-1 Step 1
    sum = 0.0
    For j = 0  To k-1 Step 1
      sum + lmat(k,j) * x(j)
    Next j
    If lmat(k,k) = 0.0
      ProcedureReturn 1
    EndIf
    x(k) = (b(k) - sum) / lmat(k,k)
  Next k

  x(n-1) / lmat(n-1,n-1)            ; Rueckwaertselimination
  For k = n - 2 To  0 Step -1
    sum = 0.0
    For j = k + 1 To n-1 Step 1
      sum + lmat(j,k) * x(j)   
    Next j
    x(k) = (x(k) - sum) / lmat(k,k)
  Next k

  ProcedureReturn 0
EndProcedure

;*====================================================================*
;*                                                                    *
;*  Die Funktion cholesky dient zur Loesung eines linearen            *
;*  Gleichungssystems:  mat * x = b                                   *
;*  mit der n x n Koeffizientenmatrix mat und der rechten Seite b     *
;*  nach dem Cholesky-Verfahren.                                      *
;*                                                                    *
;*  Die Eingabematrix muss symmetrisch und positiv definit sein,      *
;*  d.h. fuer alle n Vektoren y != 0 muss gelten:                     *
;*                                           T                        *
;*                                          y  * mat * y > 0.         *
;*                                                                    *
;*  cholesky arbeitet nur auf der unteren Dreieckmatrix incl. der     *
;*  Diagonalen, so dass gegenueber dem Gaussverfahren eine erhebliche *
;*  Reduktion des Speicherplatzes erreicht wird; es genuegt daher     *
;*  diesen Teil der Matrix an cholesky zu uebergeben.                 *
;*                                                                    *
;*====================================================================*
;*                                                                    *
;*   Anwendung:                                                       *
;*   =========                                                        *
;*                                                                    *
;*      Lineare Gleichungssysteme mit n x n Koeffizientenmatrix,      *
;*      die symmetrisch und positiv definit ist.                      *
;*                                                                    *
;*====================================================================*
;*                                                                    *
;*   Steuerparameter:                                                 *
;*   ===============                                                  *
;*      mod      Long mod;                                            *
;*               Aufrufart von cholesky:                              *
;*       = 0     Bestimmung der Zerlegungsmatrix und Berechnung       *
;*               der Loesung des Gleichungssystems.                   *
;*       = 1     Nur Berechnung der Zerlegungsmatrix; wird auf        *
;*               mat ueberspeichert.                                  *
;*       = 2     Nur Loesung des Gleichungssystems; zuvor muss je-    *
;*               doch die Zerlegungsmatrix bestimmt sein. Diese       *
;*               Aufrufart wird verwendet, falls bei gleicher         *
;*               Matrix lediglich die rechte Seite des Systems vari-  *
;*               iert, z. B. zur Berechnung der Inversen.             *
;*                                                                    *
;*   Eingabeparameter:                                                *
;*   ================                                                 *
;*      n        Long n;  ( n > 0 )                                   *
;*               Dimension von mat, Anzahl der Komponenten            *
;*               des b-Vektors u. des Loesungsvektors x.              *
;*      mat      Double  mat(n,n)                                     *
;*                 mod = 0, 1: Matrix des Gleichungssystems.          *
;*                 mod = 2   : Zerlegungsmatrix.                      *
;*      b        Double   b(n)         ( bei mod = 0, 2 )             *
;*               Rechte Seite des Gleichungssystems.                  *
;*                                                                    *
;*   Ausgabeparameter:                                                *
;*   ================                                                 *
;*      mat      Double mat(n,n)       ( bei mod = 0, 1 )             *
;*               Dekompositionsmatrix, die die Zerlegung von          *
;*               mat enthaelt.                                        *
;*      x        Double  x(n);         ( bei mod = 0, 2 )             *
;*               Loesungsvektor des Systems.                          *
;*                                                                    *
;*   Rueckgabewert:                                                   *
;*   =============                                                    *
;*      = 0      alles ok                                             *
;*      = 1      n < 1 gewaehlt oder falsche Eingabeparameter         *
;*      = 2      Matrix ist nicht positiv definit                     *
;*      = 3      Falsche Aufrufart                                    *
;*                                                                    *
;*====================================================================*
;*                                                                    *
;*   Benutzte Funktionen:                                             *
;*   ===================                                              *
;*      Num_ChoDec()  : Bestimmt die Dekomposition                    *
;*      Num_ChoSol()  : Loest das lineare Gleichungssystem            *
;*                                                                    *
;*====================================================================*
Procedure Num_Cholesky  (mod.l,n.l,Array mat.d(2), Array b.d(1), Array x.d(1))

  Protected rc.l

  If n < 1 
    ProcedureReturn 1
  EndIf

  If mod = 0  ; Zerlegung bestimmen und Gleichungssystem loesen
    rc = Num_ChoDec (n, mat())
    If rc = 0
      ProcedureReturn Num_ChoSol (n, mat(), b(), x())
    Else
      ProcedureReturn rc
    EndIf
  ElseIf mod = 1 ; Nur Zerlegung berechnen
    ProcedureReturn Num_ChoDec (n, mat())
  ElseIf mod = 2 ; Nur Loesung bestimmen
    ProcedureReturn Num_ChoSol (n, mat(), b(), x())
  EndIf 
  ProcedureReturn 3 ; Falsche Aufrufart
EndProcedure


;************************************************************************
;* die Koeffizienten eines Ausgleichspolynoms n-ten Grades nach der     *
;* diskreten Fehlerquadratmethode von Gauss berechnen.                  *
;* Das Verfahren beruht auf den Gaussschen Normalgleichungen, die mit   *
;* dem Choleskyverfahren geloest werden: Die Koeffizienten c[j] des     *
;* Vektors                                                              *
;*                                                                      *
;*        ( w[0] * (c[0] + c[1] * x[0]^1 + ... + c[n] * x[0]^n) )       *
;*        ( w[1] * (c[0] + c[1] * x[1]^1 + ... + c[n] * x[1]^n) )       *
;*        ( ...                            ...                  )       *
;*        ( w[m] * (c[0] + c[1] * x[m]^1 + ... + c[n] * x[m]^n) )       *
;*                                                                      *
;* sollen so berechnet werden, dass er der rechten Seite                *
;*                                                                      *
;*                        ( w[0] * y[0] )                               *
;*                        ( w[1] * y[1] )                               *
;*                        ( ...         )                               *
;*                        ( w[m] * y[m] )                               *
;*                                                                      *
;* in der Euklidischen Norm moeglichst nahe kommt.                      *
;* c ergibt sich dann als Loesung des linearen Gleichungssystems        *
;*                        a * c  =  b              (Normalgleichungen)  *
;* mit                                                                  *
;*      a[i][j] = w[0] * x[0]^(j+i)    + ... + w[m] * x[m]^(j+i),       *
;*      b[i]    = w[0] * y[0] * x[0]^i + ... + w[m] * y[m] * x[m]^i     *
;* fuer i,j=0(1)n.                                                      *
;*                                                                      *
;* Eingabeparameter:                                                    *
;* =================                                                    *
;* m: Nummer der letzten Stuetzstelle                                   *
;* n: Hoechstgrad des algebraischen Ausgleichspolynoms                  *
;* x: [0..m]-Vektor mit den Stuetzstellen                               *
;* y: [0..m]-Vektor mit den Funktionswerten zu den Stuetzstellen        *
;* w: [0..m]-Vektor mit den Gewichten                                   *
;*                                                                      *
;* Ausgabeparameter:                                                    *
;* =================                                                    *
;* c: [0..n]-Vektor mit den Koeffizienten des Ausgleichspolynoms        *
;*                                                                      *
;* Funktionswert:                                                       *
;* ==============                                                       *
;* 0: kein Fehler                                                       *
;* 1: Fehler in den Eingabedaten: n zu klein oder m zu klein            *
;* 2: Fehler in choly()                                                 *
;*                                                                      *
;************************************************************************
Procedure Num_ApproxPoly (n.l,m.l,Array x.d(1),Array y.d(1), Array w.d(1),Array c.d(1))

  Dim a.d(n,n)         ; die [0..n,0..n]-Matrix des Gleichungssystems
  Dim b.d(n)           ; [0..n]-Vektor mit der rechten Seite des Gleichungssystems
  Protected summand.d ; Hilfsvariable zur Berechnung eines Matrixelements
  Protected  i.l, j.l ; Laufvariablen


  If n < 0 Or m < n           ; n oder m zu klein?
    ProcedureReturn 1
  EndIf

  ; die erste Spalte, die letzte Zeile der Matrix und die rechte Seite des Gleichungssystems berechnen
  For i = 0 To n  Step 1
    a(i,0) = 0.0
    a(n,i) = 0.0
    b(i)   = 0.0
  Next i
                                   
  For j = 0 To m Step 1         

    summand = w(j)
    For i = 0 To n-1 Step 1
      a(i,0) + summand
      b(i)   + summand * y(j)
      summand * x(j)
    Next i
    a(n,0) + summand
    b(n)   + summand * y(j)
    For i = 1 To n Step 1
      summand * x(j)
      a(n,i) + summand
    Next i
  Next j

  ;die Matrix vervollstaendigen, indem immer wieder Zeile i + 1 versetzt in Zeile i kopiert wird
  For i = n - 1 To  0 Step -1
    For j = 0 To n-1 Step 1   
      a(i,j + 1) = a(i + 1,j)
    Next j
  Next i

  i = Num_Cholesky(0, n + 1, a(), b(), c());        ; das Gleichungssystem loesen

  If i <> 0              ; Fehler beim Loesen des Gleichungsystems?
    ProcedureReturn 2
  EndIf
  ProcedureReturn 0
EndProcedure
;=======================================================================================
;- Demo

#Image1 = 1
#Window = 2
#ImageGad = 3
#Button = 4
#Option1 = 5
#Option2 = 6
#MaxX = 800
#MaxY = 700
#MaxAnz = 1000

Dim X.d(#MaxAnz)
Dim Y.d(#MaxAnz)
Dim W.d(#MaxAnz)
Dim C.d(#MaxAnz)
Define i.l
Define Anzahl.l = 0
Define Grad.l = 0
Define Event.l
Define Aktion.l = #Option1
Define Input$
Define PosX.l
Define PosY.l
Define Xmin.d
Define Xmax.d
Define XA.l
Define XE.l
Define YA.l
Define YE.l

CreateImage(#Image1,#MaxX,#MaxY)

If OpenWindow(#Window,0,0,#MaxX+10,#MaxY+55,"Polynomial regression",#PB_Window_ScreenCentered)
   ImageGadget (#ImageGad,  5,10, #MaxX,#MaxY,ImageID(#Image1))
   ButtonGadget(#Button,5,#MaxY+20,100,25, "Quit")
   OptionGadget(#Option1,110,#MaxY+15,100,20,"Add Points")
   OptionGadget(#Option2,110,#MaxY+35,100,20,"Polynomial regression")
   
  StartDrawing(ImageOutput(#Image1))
  Box(0,0,#MaxX,#MaxY,RGB(50,50,50))
  StopDrawing()
  SetGadgetState(#ImageGad,ImageID(#Image1))

  SetGadgetState(#Option1,1) 

  Repeat
    Event = WaitWindowEvent()
     
    Select Event

      Case #PB_Event_Gadget
       
        Select EventGadget()
         
          Case #ImageGad
            If Aktion = #Option1 And EventType() = #PB_EventType_LeftClick And Anzahl < #MaxAnz
              PosX = WindowMouseX(#Window) - GadgetX(#ImageGad)
              PosY = WindowMouseY(#Window) - GadgetY(#ImageGad)
              X(Anzahl) = PosX
              Y(Anzahl) = PosY
              W(Anzahl) = 1.0
              Anzahl + 1

              StartDrawing(ImageOutput(#Image1))
              Box(0,0,#MaxX,#MaxY,RGB(50,50,50))
              For i = 0 To Anzahl-1   
                Box   (Int(X(i)) - 2,Int(Y(i)) - 2,4,4,RGB(250,250,0))
              Next
              StopDrawing()
              SetGadgetState(#ImageGad,ImageID(#Image1))
                           
            EndIf
         
          Case #Button
            Event = #PB_Event_CloseWindow
             
          Case #Option1
            If Aktion <> #Option1
              Anzahl = 0
              Grad = 0               
              StartDrawing(ImageOutput(#Image1))
              Box(0,0,#MaxX,#MaxY,RGB(50,50,50))
              StopDrawing()
              SetGadgetState(#ImageGad,ImageID(#Image1))               
              Aktion = #Option1
            EndIf
               
          Case #Option2
            If Aktion <> #Option2
              If Anzahl >= 2               
                Repeat
                  Grad = Anzahl-1
                  Input$ = InputRequester("Degree of polynomial", "(min.:1/max."+Str(Grad)+"):", Str(Grad))
                  Grad = Val(Input$)
                Until Grad>0 And Grad < Anzahl
                 
                i = Num_ApproxPoly(Grad,Anzahl-1,X(),Y(),W(),C())
                If i > 0
                  MessageRequester("Error!","Error (" + Str(i) + ") during the calculation!")
                  StartDrawing(ImageOutput(#Image1))
                  Box(0,0,#MaxX,#MaxY,RGB(50,50,50))
                  StopDrawing()
                  SetGadgetState(#ImageGad,ImageID(#Image1))
                  Aktion = #Option1
                  Anzahl = 0
                  SetGadgetState(#Option1,1)       
                Else
                  StartDrawing(ImageOutput(#Image1))
                  Box(0,0,#MaxX,#MaxY,RGB(50,50,50))
                  Xmin = #MaxX
                  Xmax = 0.0
                  For i = 0 To Anzahl-1   
                    Box   (Int(X(i)) - 2,Int(Y(i)) - 2,4,4,RGB(250,250,0))
                    If X(i) < Xmin
                      Xmin = X(i)
                    EndIf
                    If X(i) > Xmax
                      Xmax = X(i)
                    EndIf
                  Next
                  XA = Int(Xmin)
                  XE = Int(Xmax)
                  For i=XA To XE Step 1
                    YE = Num_Horner(Grad,C(),i)
                    If i>XA
                      LineXY(i-1,YA,i,YE,RGB(255,255,255))
                    EndIf
                    YA = YE
                  Next i
                  StopDrawing()
                  SetGadgetState(#ImageGad,ImageID(#Image1))
                  Aktion = #Option2
                EndIf
              Else
                MessageRequester("Error!","To less points")
                SetGadgetState(#Option1,1)
              EndIf
            EndIf                       
        EndSelect
    EndSelect
  Until Event = #PB_Event_CloseWindow
EndIf
CloseWindow(#Window)

Author:  zxretrosoft [ Wed Jul 17, 2019 12:11 pm ]
Post subject:  Re: Cubic regression

Yes, this is an eventuality.

But it is necessary to solve it numerically.

x=1 y=5
x=2 y=7
x=3 y=11
x=4 y=20

Result (by the calculator):
A=0.5
B=-2
C=4.5
D=1.999999999

How did the calculator calculate it?
There is a formula for ABCD.

Author:  NicTheQuick [ Wed Jul 17, 2019 12:32 pm ]
Post subject:  Re: Cubic regression

zxretrosoft wrote:
But it is necessary to solve it numerically.

Was this a question? Then no. You only have to solve it numerically if you have less input values than output values.

Author:  zxretrosoft [ Wed Jul 17, 2019 12:40 pm ]
Post subject:  Re: Cubic regression

Of course.

The answer is very very simple.

The quadratic regression formulas for ABC are:
a = { [ Σ x2 y * Σ xx ] - [Σ xy * Σ xx2 ] } / { [ Σ xx * Σ x2x 2] - [Σ xx2 ]2 }
b = { [ Σ xy * Σ x2x2 ] - [Σ x2y * Σ xx2 ] } / { [ Σ xx * Σ x2x 2] - [Σ xx2 ]2 }
c = [ Σ y / n ] - { b * [ Σ x / n ] } - { a * [ Σ x 2 / n ] }

The cubic regression formulas for ABCD are:
a =
b =
c =
d =

Author:  STARGÅTE [ Wed Jul 17, 2019 12:44 pm ]
Post subject:  Re: Cubic regression

Have you read the article?
Have you seen my code?
Of course you can use a very very long formula for A,B,C,D or you can just solve/calculate it with the matrix algorithm. The result of A, B, C and D is in the Array c.d()
Your formula for degree 2 or 3 is nothing else like this matrix formula, just rewritten.

Quote:
How did the calculator calculate it?

Which calculator?

Btw:
Code:
x(0) = 1 : y(0) = 5
x(1) = 2 : y(1) = 7
x(2) = 3 : y(2) = 11
x(3) = 4 : y(3) = 20

w(0) = 1
w(1) = 1
w(2) = 1
w(3) = 1

Num_ApproxPoly (3, 3, x(), y(), w(), c())

Debug "D = "+c(0)
Debug "C = "+c(1)
Debug "B = "+c(2)
Debug "A = "+c(3)

Quote:
D = 2
C = 4.5
B = -2
A = 0.5

Author:  zxretrosoft [ Wed Jul 17, 2019 12:47 pm ]
Post subject:  Re: Cubic regression

STARGÅTE wrote:
Have you read the article?
Have you seen my code?
Of course you can use a very very long formula for A,B,C,D or you can just solve/calculate it with the matrix algorithm. The result of A, B, C and D is in the Array c.d()
Your formula for degree 2 or 3 is nothing else like this matrix formula, just rewritten.

Quote:
How did the calculator calculate it?

Which calculator?


Yes, thank you, I'm studying the code :wink: It's difficult for me because I can't speak English. It takes me to understand something.
In the code, where can I enter those XY values?

P.S. Calculator Casio fx9860.

Author:  STARGÅTE [ Wed Jul 17, 2019 1:03 pm ]
Post subject:  Re: Cubic regression

Copy the code (without the DEMO) and add this:
Code:
; Enter here all points x1,y1 ; x1,y2 ; ... :
Define PointList.s = "1,5 ; 2,7 ; 3,11 ; 4,20"

Define Points.i = CountString(PointList, ";")
Define Dim X.d(Points)
Define Dim Y.d(Points)
Define Dim W.d(Points)
Define Dim C.d(3) ; Degree of polynom
Define I

For I = 0 To Points
   X(I) = ValD(StringField(StringField(PointList, I+1, ";"), 1, ","))
   Y(I) = ValD(StringField(StringField(PointList, I+1, ";"), 2, ","))
   W(I) = 1.0
Next

If Num_ApproxPoly(Points, 3, X(), Y(), W(), C()) = 0
   Debug "A = "+C(3)
   Debug "B = "+C(2)
   Debug "C = "+C(1)
   Debug "D = "+C(0)
Else
   Debug "Error"
EndIf
Quote:
A = 0.5
B = -2
C = 4.5
D = 2

Author:  zxretrosoft [ Wed Jul 17, 2019 1:12 pm ]
Post subject:  Re: Cubic regression

Yes it works right! Your code is great and it's the solution I'm looking for. I'll still use the code. Thank you so much STARGÅTE! 8)

But I still search for 4 ABCD equations that I can't find. I don't care if they're long, but I'd like to know them.

Author:  NicTheQuick [ Wed Jul 17, 2019 2:52 pm ]
Post subject:  Re: Cubic regression

zxretrosoft wrote:
But I still search for 4 ABCD equations that I can't find. I don't care if they're long, but I'd like to know them.

I don't think this will help you in any way. Just understand that you better should use matrix calculations with more than 3 values. :wink:

Author:  jack [ Wed Jul 17, 2019 4:58 pm ]
Post subject:  Re: Cubic regression

using Maple the symbolic coefficients are
Code:
A = y1/(x1^3-x1^2*x2-x1^2*x3-x1^2*x4+x1*x2*x3+x1*x2*x4+x1*x3*x4-x2*x3*x4)-y2/(x1*x2^2-x1*x2*x3-x1*x2*x4+x1*x3*x4-x2^3+x2^2*x3+x2^2*x4-x2*x3*x4)+y3/(x1*x2*x3-x1*x2*x4-x1*x3^2+x1*x3*x4-x2*x3^2+x2*x3*x4+x3^3-x3^2*x4)-y4/(x1*x2*x3-x1*x2*x4-x1*x3*x4+x1*x4^2-x2*x3*x4+x2*x4^2+x3*x4^2-x4^3)

B = -(x2+x3+x4)*y1/(x1^3-x1^2*x2-x1^2*x3-x1^2*x4+x1*x2*x3+x1*x2*x4+x1*x3*x4-x2*x3*x4)+(x1+x3+x4)*y2/(x1*x2^2-x1*x2*x3-x1*x2*x4+x1*x3*x4-x2^3+x2^2*x3+x2^2*x4-x2*x3*x4)-(x1+x2+x4)*y3/(x1*x2*x3-x1*x2*x4-x1*x3^2+x1*x3*x4-x2*x3^2+x2*x3*x4+x3^3-x3^2*x4)+(x1+x2+x3)*y4/(x1*x2*x3-x1*x2*x4-x1*x3*x4+x1*x4^2-x2*x3*x4+x2*x4^2+x3*x4^2-x4^3)

C = (x2*x3+x2*x4+x3*x4)*y1/(x1^3-x1^2*x2-x1^2*x3-x1^2*x4+x1*x2*x3+x1*x2*x4+x1*x3*x4-x2*x3*x4)-(x1*x3+x1*x4+x3*x4)*y2/(x1*x2^2-x1*x2*x3-x1*x2*x4+x1*x3*x4-x2^3+x2^2*x3+x2^2*x4-x2*x3*x4)+(x1*x2+x1*x4+x2*x4)*y3/(x1*x2*x3-x1*x2*x4-x1*x3^2+x1*x3*x4-x2*x3^2+x2*x3*x4+x3^3-x3^2*x4)-(x1*x2+x1*x3+x2*x3)*y4/(x1*x2*x3-x1*x2*x4-x1*x3*x4+x1*x4^2-x2*x3*x4+x2*x4^2+x3*x4^2-x4^3)

D = -x2*x3*x4*y1/(x1^3-x1^2*x2-x1^2*x3-x1^2*x4+x1*x2*x3+x1*x2*x4+x1*x3*x4-x2*x3*x4)+x4*x3*x1*y2/(x1*x2^2-x1*x2*x3-x1*x2*x4+x1*x3*x4-x2^3+x2^2*x3+x2^2*x4-x2*x3*x4)-x1*x2*x4*y3/(x1*x2*x3-x1*x2*x4-x1*x3^2+x1*x3*x4-x2*x3^2+x2*x3*x4+x3^3-x3^2*x4)+x1*x2*x3*y4/(x1*x2*x3-x1*x2*x4-x1*x3*x4+x1*x4^2-x2*x3*x4+x2*x4^2+x3*x4^2-x4^3)

you could probably rewrite them using sums and therefore making the expressions shorter and more aesthetically pleasing, I leave that up to you
one could achieve the same result using the free Maxima CAS https://wxmaxima-developers.github.io/wxmaxima/
maxima code follows
Code:
load ("eigen")$

m:matrix([1,x1,x1^2,x1^3],[1,x2,x2^2,x2^3], [1,x3,x3^2,x3^3], [1,x4,x4^2,x4^3]);
y:matrix([y1],[y2], [y3], [y4]);
b:invert_by_lu(m).y$

A:ratsimp(b[4][1]);
B:ratsimp(b[3][1]);
C:ratsimp(b[2][1]);
D:ratsimp(b[1][1]);

subst([x1=1, x2=2,x3=3,x4=4,y1=5,y2=7,y3=11,y4=20],A);
subst([x1=1, x2=2,x3=3,x4=4,y1=5,y2=7,y3=11,y4=20],B);
subst([x1=1, x2=2,x3=3,x4=4,y1=5,y2=7,y3=11,y4=20],C);
subst([x1=1, x2=2,x3=3,x4=4,y1=5,y2=7,y3=11,y4=20],D);

Author:  zxretrosoft [ Wed Jul 17, 2019 9:23 pm ]
Post subject:  Re: Cubic regression

Yes jack! This is exactly it! 8)

But it seems to me that it's just for 4 input values. Can we generalize it anyway? (i.e. for any number of values) ?

Author:  NicTheQuick [ Wed Jul 17, 2019 9:55 pm ]
Post subject:  Re: Cubic regression

The generalization is the matrix calculation STARGÅTE posted.

Page 1 of 1 All times are UTC + 1 hour
Powered by phpBB © 2000, 2002, 2005, 2007 phpBB Group
http://www.phpbb.com/