PureBasic Forum

 It is currently Thu Jun 04, 2020 11:49 pm

 All times are UTC + 1 hour

 Page 1 of 1 [ 15 posts ]
 Print view Previous topic | Next topic
Author Message
 Post subject: Cubic regressionPosted: Wed Jul 17, 2019 6:47 am
 Enthusiast

Joined: Wed May 15, 2013 8:26 am
Posts: 165
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:

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?

_________________

Top

 Post subject: Re: Cubic regressionPosted: Wed Jul 17, 2019 10:24 am

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

_________________

Top

 Post subject: Re: Cubic regressionPosted: Wed Jul 17, 2019 10:41 am
 Enthusiast

Joined: Wed May 15, 2013 8:26 am
Posts: 165
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:

_________________

Top

 Post subject: Re: Cubic regressionPosted: Wed Jul 17, 2019 11:02 am

Joined: Thu Jan 10, 2008 1:30 pm
Posts: 1273
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

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                                     *
;* 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 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

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
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
#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 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)

StartDrawing(ImageOutput(#Image1))
Box(0,0,#MaxX,#MaxY,RGB(50,50,50))
StopDrawing()

Repeat
Event = WaitWindowEvent()

Select Event

If Aktion = #Option1 And EventType() = #PB_EventType_LeftClick And Anzahl < #MaxAnz
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()

EndIf

Case #Button
Event = #PB_Event_CloseWindow

Case #Option1
If Aktion <> #Option1
Anzahl = 0
StartDrawing(ImageOutput(#Image1))
Box(0,0,#MaxX,#MaxY,RGB(50,50,50))
StopDrawing()
Aktion = #Option1
EndIf

Case #Option2
If Aktion <> #Option2
If Anzahl >= 2
Repeat

If i > 0
MessageRequester("Error!","Error (" + Str(i) + ") during the calculation!")
StartDrawing(ImageOutput(#Image1))
Box(0,0,#MaxX,#MaxY,RGB(50,50,50))
StopDrawing()
Aktion = #Option1
Anzahl = 0
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
If i>XA
LineXY(i-1,YA,i,YE,RGB(255,255,255))
EndIf
YA = YE
Next i
StopDrawing()
Aktion = #Option2
EndIf
Else
MessageRequester("Error!","To less points")
EndIf
EndIf
EndSelect
EndSelect
Until Event = #PB_Event_CloseWindow
EndIf
CloseWindow(#Window)

_________________

Top

 Post subject: Re: Cubic regressionPosted: Wed Jul 17, 2019 12:11 pm
 Enthusiast

Joined: Wed May 15, 2013 8:26 am
Posts: 165
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.

_________________

Top

 Post subject: Re: Cubic regressionPosted: Wed Jul 17, 2019 12:32 pm
 Enthusiast

Joined: Sun Jun 22, 2003 7:43 pm
Posts: 523
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.

_________________
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

 Post subject: Re: Cubic regressionPosted: Wed Jul 17, 2019 12:40 pm
 Enthusiast

Joined: Wed May 15, 2013 8:26 am
Posts: 165
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 =

_________________

Top

 Post subject: Re: Cubic regressionPosted: Wed Jul 17, 2019 12:44 pm

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

_________________

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

Top

 Post subject: Re: Cubic regressionPosted: Wed Jul 17, 2019 12:47 pm
 Enthusiast

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

_________________

Top

 Post subject: Re: Cubic regressionPosted: Wed Jul 17, 2019 1:03 pm

Joined: Thu Jan 10, 2008 1:30 pm
Posts: 1273
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

_________________

Top

 Post subject: Re: Cubic regressionPosted: Wed Jul 17, 2019 1:12 pm
 Enthusiast

Joined: Wed May 15, 2013 8:26 am
Posts: 165
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!

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.

_________________

Top

 Post subject: Re: Cubic regressionPosted: Wed Jul 17, 2019 2:52 pm
 Enthusiast

Joined: Sun Jun 22, 2003 7:43 pm
Posts: 523
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.

_________________
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

 Post subject: Re: Cubic regressionPosted: Wed Jul 17, 2019 4:58 pm

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

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

 Post subject: Re: Cubic regressionPosted: Wed Jul 17, 2019 9:23 pm
 Enthusiast

Joined: Wed May 15, 2013 8:26 am
Posts: 165
Location: Czech Republic, Prague
Yes jack! This is exactly it!

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) ?

_________________

Top

 Post subject: Re: Cubic regressionPosted: Wed Jul 17, 2019 9:55 pm
 Enthusiast

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

Top

 Display posts from previous: All posts1 day7 days2 weeks1 month3 months6 months1 year Sort by AuthorPost timeSubject AscendingDescending
 Page 1 of 1 [ 15 posts ]

 All times are UTC + 1 hour

Who is online

Users browsing this forum: No registered users and 5 guests

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

Search for:
 Jump to:  Select a forum ------------------ PureBasic    Coding Questions    Game Programming    3D Programming    Assembly Programming    The PureBasic Editor    The PureBasic Form Designer    General Discussion    Feature Requests and Wishlists    Tricks 'n' Tips Bug Reports    Bugs - Windows    Bugs - Linux    Bugs - Mac OSX    Bugs - IDE    Bugs - Documentation OS Specific    AmigaOS    Linux    Windows    Mac OSX Miscellaneous    Announcement    Off Topic Showcase    Applications - Feedback and Discussion    PureFORM & JaPBe    TailBite