It is currently Wed Sep 18, 2019 9:26 pm

All times are UTC + 1 hour




Post new topic Reply to topic  [ 15 posts ] 
Author Message
 Post subject: Cubic regression
PostPosted: Wed Jul 17, 2019 6:47 am 
Offline
Enthusiast
Enthusiast
User avatar

Joined: Wed May 15, 2013 8:26 am
Posts: 150
Location: Czech Republic, Prague
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:

_________________
I apologize in advance for bad English


Top
 Profile  
Reply with quote  
 Post subject: Re: Cubic regression
PostPosted: Wed Jul 17, 2019 10:24 am 
Offline
Addict
Addict
User avatar

Joined: Thu Jan 10, 2008 1:30 pm
Posts: 1213
Location: Germany, Glienicke
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.

_________________
ImageImage


Top
 Profile  
Reply with quote  
 Post subject: Re: Cubic regression
PostPosted: Wed Jul 17, 2019 10:41 am 
Offline
Enthusiast
Enthusiast
User avatar

Joined: Wed May 15, 2013 8:26 am
Posts: 150
Location: Czech Republic, Prague
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

_________________
I apologize in advance for bad English


Top
 Profile  
Reply with quote  
 Post subject: Re: Cubic regression
PostPosted: Wed Jul 17, 2019 11:02 am 
Offline
Addict
Addict
User avatar

Joined: Thu Jan 10, 2008 1:30 pm
Posts: 1213
Location: Germany, Glienicke
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)

_________________
ImageImage


Top
 Profile  
Reply with quote  
 Post subject: Re: Cubic regression
PostPosted: Wed Jul 17, 2019 12:11 pm 
Offline
Enthusiast
Enthusiast
User avatar

Joined: Wed May 15, 2013 8:26 am
Posts: 150
Location: Czech Republic, Prague
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


Top
 Profile  
Reply with quote  
 Post subject: Re: Cubic regression
PostPosted: Wed Jul 17, 2019 12:32 pm 
Offline
Enthusiast
Enthusiast
User avatar

Joined: Sun Jun 22, 2003 7:43 pm
Posts: 426
Location: Germany, Saarbrücken
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.

_________________
Electronics, Crazy & Interesting Stuff, all that with text, image and sound? Click here!

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.


Top
 Profile  
Reply with quote  
 Post subject: Re: Cubic regression
PostPosted: Wed Jul 17, 2019 12:40 pm 
Offline
Enthusiast
Enthusiast
User avatar

Joined: Wed May 15, 2013 8:26 am
Posts: 150
Location: Czech Republic, Prague
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


Top
 Profile  
Reply with quote  
 Post subject: Re: Cubic regression
PostPosted: Wed Jul 17, 2019 12:44 pm 
Offline
Addict
Addict
User avatar

Joined: Thu Jan 10, 2008 1:30 pm
Posts: 1213
Location: Germany, Glienicke
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

_________________
ImageImage


Last edited by STARGÅTE on Wed Jul 17, 2019 12:51 pm, edited 2 times in total.

Top
 Profile  
Reply with quote  
 Post subject: Re: Cubic regression
PostPosted: Wed Jul 17, 2019 12:47 pm 
Offline
Enthusiast
Enthusiast
User avatar

Joined: Wed May 15, 2013 8:26 am
Posts: 150
Location: Czech Republic, Prague
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.

_________________
I apologize in advance for bad English


Top
 Profile  
Reply with quote  
 Post subject: Re: Cubic regression
PostPosted: Wed Jul 17, 2019 1:03 pm 
Offline
Addict
Addict
User avatar

Joined: Thu Jan 10, 2008 1:30 pm
Posts: 1213
Location: Germany, Glienicke
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

_________________
ImageImage


Top
 Profile  
Reply with quote  
 Post subject: Re: Cubic regression
PostPosted: Wed Jul 17, 2019 1:12 pm 
Offline
Enthusiast
Enthusiast
User avatar

Joined: Wed May 15, 2013 8:26 am
Posts: 150
Location: Czech Republic, Prague
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


Top
 Profile  
Reply with quote  
 Post subject: Re: Cubic regression
PostPosted: Wed Jul 17, 2019 2:52 pm 
Offline
Enthusiast
Enthusiast
User avatar

Joined: Sun Jun 22, 2003 7:43 pm
Posts: 426
Location: Germany, Saarbrücken
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:

_________________
Electronics, Crazy & Interesting Stuff, all that with text, image and sound? Click here!

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.


Top
 Profile  
Reply with quote  
 Post subject: Re: Cubic regression
PostPosted: Wed Jul 17, 2019 4:58 pm 
Offline
Addict
Addict

Joined: Fri Apr 25, 2003 11:10 pm
Posts: 1199
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);


Top
 Profile  
Reply with quote  
 Post subject: Re: Cubic regression
PostPosted: Wed Jul 17, 2019 9:23 pm 
Offline
Enthusiast
Enthusiast
User avatar

Joined: Wed May 15, 2013 8:26 am
Posts: 150
Location: Czech Republic, Prague
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


Top
 Profile  
Reply with quote  
 Post subject: Re: Cubic regression
PostPosted: Wed Jul 17, 2019 9:55 pm 
Offline
Enthusiast
Enthusiast
User avatar

Joined: Sun Jun 22, 2003 7:43 pm
Posts: 426
Location: Germany, Saarbrücken
The generalization is the matrix calculation STARGÅTE posted.

_________________
Electronics, Crazy & Interesting Stuff, all that with text, image and sound? Click here!

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.


Top
 Profile  
Reply with quote  
Display posts from previous:  Sort by  
Post new topic Reply to topic  [ 15 posts ] 

All times are UTC + 1 hour


Who is online

Users browsing this forum: No registered users and 12 guests


You cannot post new topics in this forum
You cannot reply to topics in this forum
You cannot edit your posts in this forum
You cannot delete your posts in this forum

Search for:
Jump to:  

 


Powered by phpBB © 2008 phpBB Group
subSilver+ theme by Canver Software, sponsor Sanal Modifiye