Cubic regression

Everything else that doesn't fall into one of the other PB categories.
User avatar
zxretrosoft
Enthusiast
Enthusiast
Posts: 169
Joined: Wed May 15, 2013 8:26 am
Location: Czech Republic, Prague
Contact:

Cubic regression

Post by zxretrosoft »

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/statist ... ession.php

I found coefficients A, B, C for quadratic:
https://www.easycalculation.com/statist ... ession.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:
I apologize in advance for bad English
https://zxretrosoft.cz/
User avatar
STARGÅTE
Addict
Addict
Posts: 2067
Joined: Thu Jan 10, 2008 1:30 pm
Location: Germany, Glienicke
Contact:

Re: Cubic regression

Post by STARGÅTE »

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.
PB 6.01 ― Win 10, 21H2 ― Ryzen 9 3900X, 32 GB ― NVIDIA GeForce RTX 3080 ― Vivaldi 6.0 ― www.unionbytes.de
Lizard - Script language for symbolic calculations and moreTypeface - Sprite-based font include/module
User avatar
zxretrosoft
Enthusiast
Enthusiast
Posts: 169
Joined: Wed May 15, 2013 8:26 am
Location: Czech Republic, Prague
Contact:

Re: Cubic regression

Post by zxretrosoft »

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/statist ... ession.php
I apologize in advance for bad English
https://zxretrosoft.cz/
User avatar
STARGÅTE
Addict
Addict
Posts: 2067
Joined: Thu Jan 10, 2008 1:30 pm
Location: Germany, Glienicke
Contact:

Re: Cubic regression

Post by STARGÅTE »

Here is a updated version from "alter Mann":
https://www.purebasic.fr/german/viewtop ... roximation

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: Select all

; 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) 
PB 6.01 ― Win 10, 21H2 ― Ryzen 9 3900X, 32 GB ― NVIDIA GeForce RTX 3080 ― Vivaldi 6.0 ― www.unionbytes.de
Lizard - Script language for symbolic calculations and moreTypeface - Sprite-based font include/module
User avatar
zxretrosoft
Enthusiast
Enthusiast
Posts: 169
Joined: Wed May 15, 2013 8:26 am
Location: Czech Republic, Prague
Contact:

Re: Cubic regression

Post by zxretrosoft »

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.
I apologize in advance for bad English
https://zxretrosoft.cz/
User avatar
NicTheQuick
Addict
Addict
Posts: 1223
Joined: Sun Jun 22, 2003 7:43 pm
Location: Germany, Saarbrücken
Contact:

Re: Cubic regression

Post by NicTheQuick »

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.
The english grammar is freeware, you can use it freely - But it's not Open Source, i.e. you can not change it or publish it in altered way.
User avatar
zxretrosoft
Enthusiast
Enthusiast
Posts: 169
Joined: Wed May 15, 2013 8:26 am
Location: Czech Republic, Prague
Contact:

Re: Cubic regression

Post by zxretrosoft »

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 =
I apologize in advance for bad English
https://zxretrosoft.cz/
User avatar
STARGÅTE
Addict
Addict
Posts: 2067
Joined: Thu Jan 10, 2008 1:30 pm
Location: Germany, Glienicke
Contact:

Re: Cubic regression

Post by STARGÅTE »

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.
How did the calculator calculate it?
Which calculator?

Btw:

Code: Select all

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)
D = 2
C = 4.5
B = -2
A = 0.5
Last edited by STARGÅTE on Wed Jul 17, 2019 12:51 pm, edited 2 times in total.
PB 6.01 ― Win 10, 21H2 ― Ryzen 9 3900X, 32 GB ― NVIDIA GeForce RTX 3080 ― Vivaldi 6.0 ― www.unionbytes.de
Lizard - Script language for symbolic calculations and moreTypeface - Sprite-based font include/module
User avatar
zxretrosoft
Enthusiast
Enthusiast
Posts: 169
Joined: Wed May 15, 2013 8:26 am
Location: Czech Republic, Prague
Contact:

Re: Cubic regression

Post by zxretrosoft »

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.
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.
I apologize in advance for bad English
https://zxretrosoft.cz/
User avatar
STARGÅTE
Addict
Addict
Posts: 2067
Joined: Thu Jan 10, 2008 1:30 pm
Location: Germany, Glienicke
Contact:

Re: Cubic regression

Post by STARGÅTE »

Copy the code (without the DEMO) and add this:

Code: Select all

; 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
A = 0.5
B = -2
C = 4.5
D = 2
PB 6.01 ― Win 10, 21H2 ― Ryzen 9 3900X, 32 GB ― NVIDIA GeForce RTX 3080 ― Vivaldi 6.0 ― www.unionbytes.de
Lizard - Script language for symbolic calculations and moreTypeface - Sprite-based font include/module
User avatar
zxretrosoft
Enthusiast
Enthusiast
Posts: 169
Joined: Wed May 15, 2013 8:26 am
Location: Czech Republic, Prague
Contact:

Re: Cubic regression

Post by zxretrosoft »

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.
I apologize in advance for bad English
https://zxretrosoft.cz/
User avatar
NicTheQuick
Addict
Addict
Posts: 1223
Joined: Sun Jun 22, 2003 7:43 pm
Location: Germany, Saarbrücken
Contact:

Re: Cubic regression

Post by NicTheQuick »

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:
The english grammar is freeware, you can use it freely - But it's not Open Source, i.e. you can not change it or publish it in altered way.
jack
Addict
Addict
Posts: 1336
Joined: Fri Apr 25, 2003 11:10 pm

Re: Cubic regression

Post by jack »

using Maple the symbolic coefficients are

Code: Select all

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: Select all

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);
User avatar
zxretrosoft
Enthusiast
Enthusiast
Posts: 169
Joined: Wed May 15, 2013 8:26 am
Location: Czech Republic, Prague
Contact:

Re: Cubic regression

Post by zxretrosoft »

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) ?
I apologize in advance for bad English
https://zxretrosoft.cz/
User avatar
NicTheQuick
Addict
Addict
Posts: 1223
Joined: Sun Jun 22, 2003 7:43 pm
Location: Germany, Saarbrücken
Contact:

Re: Cubic regression

Post by NicTheQuick »

The generalization is the matrix calculation STARGÅTE posted.
The english grammar is freeware, you can use it freely - But it's not Open Source, i.e. you can not change it or publish it in altered way.
Post Reply