I tweaked your code to use the Intel Decimal Floating-Point Math Library viewtopic.php?t=29203
I hope you guys don't mind, if this post is a bother then feel free to delete it
Code: Select all
#include once "libbid.bi"
type range
as BID.Decimal128 _max,_min,sd
as long _maxi,_mini
end type
function mean(a() as BID.Decimal128,R as range) as BID.Decimal128
'R=type(-1e10,1e10,0,0,0)
R._min="1e-10"
R._max="-1e-10"
R.sd=0
R._maxi=0
R._mini=0
dim as long lb=lbound(a),ub=ubound(a)
dim as BID.Decimal128 acc,m
for n as long=lb to ub
acc+=a(n)
if R._max<a(n) then R._max=a(n):R._maxi=n
if R._min>a(n) then R._min=a(n):R._mini=n
next
m=acc/(ub-lb+1)
acc=0
for n as long=lb to ub
acc+=(a(n)-m)*(a(n)-m)
next
R.sd =sqr(acc/(ub-lb))
return m
end function
Type vector
Dim As BID.Decimal128 element(Any)
End Type
Type matrix
Dim As BID.Decimal128 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 BID.Decimal128
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 BID.Decimal128
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
Next p2
Next p1
#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
Next z
Function = ans
End Function
Function matrix.inverse() As matrix
Var ub1=Ubound(this.element,1),ub2=Ubound(this.element,2)
Dim As matrix ans
Dim As vector temp,null_
Redim temp.element(1 To ub1):Redim null_.element(1 To ub1)
Redim ans.element(1 To ub1,1 To ub2)
For a As Integer=1 To ub1
temp=null_
temp.element(a)=1
temp=GaussJordan(temp)
For b As Integer=1 To ub1
ans.element(b,a)=temp.element(b)
Next b
Next a
Return ans
End Function
'vandermode of x
Function vandermonde(x_values() As BID.Decimal128,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 BID.Decimal128,y_values() As BID.Decimal128,ans() As BID.Decimal128,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 BID.Decimal128,Byval x As BID.Decimal128) As BID.Decimal128
Dim As BID.Decimal128 acc
For i As Long=Ubound(Coefficients) To Lbound(Coefficients) Step -1
acc=acc*x+Coefficients(i)
Next i
Return acc
End Function
Dim As BID.Decimal128 sample()
Dim As Long f = Freefile
Dim As Double d
Open "Dinosaur.txt" For Input As #f
Do Until EOF(f)
Redim Preserve sample(Ubound(sample) + 1)
Input #f, d
sample(Ubound(sample))=d
Loop
Close #f
Screenres 1070,500
width 1070\8,500\16
dim as BID.Decimal128 x(),y()
dim as long counter
For I As Long = 0 To Ubound(sample)
counter+=1
redim preserve x(1 to counter)
redim preserve y(1 to counter)
x(counter)=i
y(counter)=sample(i)
If I = 0 Then
Pset(I+200, sample(I) - 26500)
Else
Line -(I+200, sample(I) - 26500)
End If
Next I
Redim As BID.Decimal128 ans()
redim as BID.Decimal128 result(0 to Ubound(sample))
regress(x(),y(),ans(),9)
For x As Single=0 To Ubound(sample)
result(x)=polyeval(ans(),x)
Var f=polyeval(ans(),x)-26500
If x=0 Then Pset(x+200,f),5 Else Line-(x+200,f),2
Next
dim as range R
dim as BID.Decimal128 m
var mn=mean(sample(),R)
print "Mean ";mn;" Standard dev. ";R.sd;" Min ";R._min;" Max ";R._max;" difference ";R._max-R._min
print
mn=mean(result(),R)
color 2
print "Mean ";mn;" Standard dev. ";R.sd;" Min ";R._min;" Max ";R._max;" difference ";R._max-R._min
Sleep
question
could this problem be considered as noise in data?
if so, then perhaps an algorithm to reduce noise could be used