Anonymous user
Polynomial regression: Difference between revisions
no edit summary
No edit summary |
|||
Line 470:
3.0000
</pre>
=={{header|freebasic}}==
<lang freebasic>
Sub GaussJordan(matrix() As Double,rhs() As Double,ans() As Double)
Dim As Integer n=Ubound(matrix,1)
Redim ans(0):Redim ans(1 To n)
Dim As Double b(1 To n,1 To n),r(1 To n)
For c As Integer=1 To n 'take copies
r(c)=rhs(c)
For d As Integer=1 To n
b(c,d)=matrix(c,d)
Next d
Next c
#macro pivot(num)
For p1 As Integer = num To n - 1
For p2 As Integer = p1 + 1 To n
If Abs(b(p1,num))<Abs(b(p2,num)) Then
Swap r(p1),r(p2)
For g As Integer=1 To n
Swap b(p1,g),b(p2,g)
Next g
End If
Next p2
Next p1
#endmacro
For k As Integer=1 To n-1
pivot(k) 'full pivoting
For row As Integer =k To n-1
If b(row+1,k)=0 Then Exit For
Var f=b(k,k)/b(row+1,k)
r(row+1)=r(row+1)*f-r(k)
For g As Integer=1 To n
b((row+1),g)=b((row+1),g)*f-b(k,g)
Next g
Next row
Next k
'back substitute
For z As Integer=n To 1 Step -1
ans(z)=r(z)/b(z,z)
For j As Integer = n To z+1 Step -1
ans(z)=ans(z)-(b(z,j)*ans(j)/b(z,z))
Next j
Next z
End Sub
'Interpolate through points.
Sub Interpolate(x_values() As Double,y_values() As Double,p() As Double)
Var n=Ubound(x_values)
Redim p(0):Redim p(1 To n)
Dim As Double matrix(1 To n,1 To n),rhs(1 To n)
For a As Integer=1 To n
rhs(a)=y_values(a)
For b As Integer=1 To n
matrix(a,b)=x_values(a)^(b-1)
Next b
Next a
'Solve the linear equations
GaussJordan(matrix(),rhs(),p())
End Sub
'Evaluate a polynomial at x by Horner method
Function polyeval(Coefficients() As Double,byval x As Double) As Double
Dim As Double acc
For i As Integer=Ubound(Coefficients) To Lbound(Coefficients) Step -1
acc=acc*x+Coefficients(i)
Next i
Return acc
End Function
'======================== SET UP THE POINTS ===============
Dim As Double x(1 To ...)={0,1,2,3,4,5,6,7,8,9,10}
Dim As Double y(1 To ...)={1,6,17,34,57,86,121,162,209,262,321}
Redim As Double Poly(0)
'Get the polynomial Poly()
Interpolate(x(),y(),Poly())
'print coefficients to console
print "Polynomial Coefficients:"
print
For z As Integer=1 To Ubound(Poly)
If z=1 Then
Print "constant term ";tab(20);Poly(z)
Else
Print tab(8); "x^";z-1;" = ";tab(20);Poly(z)
End If
Next z
sleep
</lang>
<pre>
Polynomial Coefficients:
constant term 1
x^ 1 = 2
x^ 2 = 3
x^ 3 = 0
x^ 4 = 0
x^ 5 = 0
x^ 6 = 0
x^ 7 = 0
x^ 8 = 0
x^ 9 = 0
x^ 10 = 0
</pre>
=={{header|GAP}}==
| |||