Anonymous user
Polynomial regression: Difference between revisions
no edit summary
m (→{{header|Phix}}: added syntax colouring the hard way) |
No edit summary |
||
Line 685:
=={{header|FreeBASIC}}==
<lang FreeBASIC>
#include "crt.bi" 'for rounding only
Type vector
End Type
Type matrix
Dim As Double element(Any,Any)
Declare Function inverse() As matrix
Declare Function transpose() As matrix
private:
Declare Function GaussJordan(As vector) As vector
End Type
'mult operators
Operator *(m1 As matrix,m2 As matrix) As matrix
Dim rows As Integer=Ubound(m1.element,1)
Dim columns As Integer=Ubound(m2.element,2)
If Ubound(m1.element,2)<>Ubound(m2.element,1) Then
Print "Can't do"
Exit Operator
End If
Dim As matrix ans
Redim ans.element(rows,columns)
Dim rxc As Double
For r As Integer=1 To rows
For c As Integer=1 To columns
rxc=0
For k As Integer = 1 To Ubound(m1.element,2)
rxc=rxc+m1.element(r,k)*m2.element(k,c)
Next k
ans.element(r,c)=rxc
Next c
Next r
Operator= ans
End Operator
Operator *(m1 As matrix,m2 As vector) As vector
Dim rows As Integer=Ubound(m1.element,1)
Dim columns As Integer=Ubound(m2.element,2)
If Ubound(m1.element,2)<>Ubound(m2.element) Then
Print "Can't do"
Exit Operator
End If
Dim As vector ans
Redim ans.element(rows)
Dim rxc As Double
For r As Integer=1 To rows
rxc=0
For k As Integer = 1 To Ubound(m1.element,2)
rxc=rxc+m1.element(r,k)*m2.element(k)
Next k
ans.element(r)=rxc
Next r
Operator= ans
End Operator
Function matrix.transpose() As matrix
Dim As matrix b
Redim b.element(1 To Ubound(this.element,2),1 To Ubound(this.element,1))
For i As Long=1 To Ubound(this.element,1)
For j As Long=1 To Ubound(this.element,2)
b.element(j,i)=this.element(i,j)
Next
Next
Return b
End Function
Function matrix.GaussJordan(rhs As vector) As vector
Dim As Integer n=Ubound(rhs.element)
Dim As vector ans=rhs,r=rhs
Dim As matrix b=This
#macro pivot(num)
For p1 As Integer = num To n - 1
For p2 As Integer = p1 + 1 To n
If Abs(b.element(p1,num))<Abs(b.element(p2,num)) Then
Swap r.element(p1),r.element(p2)
For g As Integer=1 To n
Swap b.element(p1,g),b.element(p2,g)
Next g
End If
Line 708 ⟶ 773:
#endmacro
For k As Integer=1 To n-1
pivot(k)
For row As Integer =k To n-1
If b.element(row+1,k)=0 Then Exit For
Var f=b.element(k,k)/b.element(row+1,k)
r.element(row+1)=r.element(row+1)*f-r.element(k)
For g As Integer=1 To n
b.element((row+1),g)=b.element((row+1),g)*f-b.element(k,g)
Next g
Next row
Next k
'back substitute
For z As Integer=n To 1 Step -1
ans.element(z)=r.element(z)/b.element(z,z)
For j As Integer = n To z+1 Step -1
ans.element(z)=ans.element(z)-(b.element(z,j)*ans.element(j)/b.element(z,z))
Next j
End Function
Function matrix.inverse() As matrix
Var ub1=Ubound(this.element,1),ub2=Ubound(this.element,2)
Dim As matrix
Dim As vector
Redim temp.element(1 To ub1):Redim null_.element(1 To ub1)
Redim ans.element(1 To ub1,1 To ub2)
For a As
ans.element(b,a)=temp.element(b)
Next b
Next a
Return ans
End Function
'vandermode of x
Function vandermonde(x_values() As Double,w As Long) As matrix
Dim As matrix mat
Var n=Ubound(x_values)
Redim mat.element(1 To n,1 To w)
For a As Integer=1 To n
For b As Integer=1 To w
mat.element(a,b)=x_values(a)^(b-1)
Next b
Next a
Return mat
End Function
'main preocedure
Sub regress(x_values() As Double,y_values() As Double,ans() As Double,n As Long)
Redim ans(1 To Ubound(x_values))
Dim As matrix m1= vandermonde(x_values(),n)
Dim As matrix T=m1.transpose
Dim As vector y
Redim y.element(1 To Ubound(ans))
For n As Long=1 To Ubound(y_values)
y.element(n)=y_values(n)
Next n
Dim As vector result=(((T*m1).inverse)*T)*y
Redim Preserve ans(1 To n)
For n As Long=1 To Ubound(ans)
ans(n)=result.element(n)
Next n
End Sub
'Evaluate a polynomial at x
Function polyeval(Coefficients() As Double,Byval x As Double) As Double
Dim As Double acc
For i As Long=Ubound(Coefficients) To Lbound(Coefficients) Step -1
acc=acc*x+Coefficients(i)
Next i
Return acc
End Function
Function CRound(Byval x As Double,Byval precision As Integer=30) As String
If precision>30 Then precision=30
Dim As zstring * 40 z:Var s="%." &str(Abs(precision)) &"f"
sprintf(z,s,x)
If Val(z) Then Return Rtrim(Rtrim(z,"0"),".")Else Return "0"
End Function
Function show(a() As Double,places as long=10) As String
Dim As String s,g
For n As Long=Lbound(a) To Ubound(a)
If n<3 Then g="" Else g="^"+Str(n-1)
if val(cround(a(n),places))<>0 then
s+= Iif(Sgn(a(n))>=0,"+","")+cround(a(n),places)+ Iif(n=Lbound(a),"","*x"+g)+" "
end if
Next n
Return s
End Function
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 ans()
regress(x(),y(),ans(),3)
print show(ans())
sleep
</lang>
{{out}}
<pre>
+1 +2*x +3*x^2
</pre>
=={{header|GAP}}==
| |||