Linear trendline

[Lauderdale]

Hello I am wondering if I can get a linear fit to an array of data and also an r^2 value to determine how linear the data actually is. Or do I just have to buy excel.

[Henko]

Hi, here's a least squares fitting program for polynomiums up to 10 degrees (or higher if you change some dimensions).
The functions are taken from my vector/matrix library "matlib", which is on page 3 of this forum section. The correlation r^2 is not calculated but that should not a problem with the provided functions.

Code: Select all

option base 1 ! randomize
maxx=screen_width() ! maxy=screen_height()
graphics
fill color 0.8,0.8,0.8 ! draw color 0,0,0
fill rect 0,0 to maxx,maxy
degree=2 ! gosub curve_fit  ' set degree of the fit and solve
end

curve_fit:
dim am(20,20),bm(20,20),cm(20,20),av(40),bv(40)
read n                   
for i=1 to n ! read av(i),bv(i) ! next i ' read some point coordinates
points:
data  12,-35,40,-30,20,-25,-10,-15,-25, -10,-10,0,-7
data  10,-15,20,-20,25,-10,30,0,35,20,40,45
vec_uit ("   x ",n,av,10,10,3,1)
vec_uit ("     y ",n,bv,130,10,6,0)
plot(n,av,bv,100,550,600,900,degree)
return

' perform a least squares fit and plot the result
' mode=-1 : plot points and connect them with straight lines
' mode=0  : plot points only
' mode>0  : plot points a best fit polynomium of degree mode
'
def plot (np,xp(),yp(),x1,y1,x2,y2,mode)
dim xc(40),yc(40),coef(11)
xmin=xp(1) ! xmax=xmin ! ymin=yp(1) ! ymax=ymin
for i=2 to np
  if xmin>xp(i) then xmin=xp(i)
  if xmax<xp(i) then xmax=xp(i)
  if ymin>yp(i) then ymin=yp(i)
  if ymax<yp(i) then ymax=yp(i)
  next i
scalex=(x2-x1)/(xmax-xmin) ! scaley=(y2-y1)/(ymax-ymin)
if xmin>0 then
  x_as=0
  for i=1 to np ! xp(i)=xp(i)-xmin ! next i
  end if
if xmax<0 then x_as=1
if xmin<0 and xmax>0 then x_as=-xmin/(xmax-xmin)
if ymin>0 then
  y_as=1
  for i=1 to np ! yp(i)=yp(i)-ymin ! next i
  end if
if ymax<0 then y_as=0
if ymin<0 and ymax>0 then y_as=1+ymin/(ymax-ymin)
xy_grid (x1,y1,x2,y2,x_as,y_as,12,12)
x_as=x1+x_as*(x2-x1)! y_as=y1+y_as*(y2-y1)
fill color 0,0,0
for i=1 to np
  xc(i)=x_as+scalex*xp(i)
  yc(i)=y_as-scaley*yp(i)
  fill circle xc(i)-3,yc(i)-3 to xc(i)+3,yc(i)+3
  next i
if mode<0 then
  draw color 1,0,0 ! draw size 1 ! draw to xc(1),yc(1)
  for i=2 to np ! draw line to xc(i),yc(i) ! next i
  draw color 0,0,0
  endif 
if mode>0 then
  draw color 1,0,0 ! draw size 1
  ls (np,mode,xp,yp,coef)
  draw to x_as+scalex*xmin,y_as-scaley*poly(mode,coef,xmin)
  dx=(xmax-xmin)/48
  for x=xmin+dx to xmax step dx
    draw line to x_as+scalex*x,y_as-scaley*poly(mode,coef,x)
  next x
  draw color 0,0,0
  vec_uit ("   coef'n ",mode+1,coef,270,10,6,4)
end if
end def

' this is a m-degree polynomial least squares fit of a number of
' point in 2 dimensional space
' there are n points, with x- and y-coordinates in vectors x() and y()
' m is the degree of the polynomial (take 1 for straight line fit,
' m=2 for parabola, and so on)
' the coefficients of the best fit polynomial are returned in vector c()
' f.i. for m=2 : y = c(1) + c(2)*x + c(3)*x^2
'
def ls (n,m,x(),y(),c())
dim a1(20,20),a2(20,20),a3(20,20),rl(20)
m=m+1
for i=1 to n
  a1(i,1)=1
  for j=2 to m
    a1(i,j)=a1(i,j-1)*x(i)
  next j
next i
mat_trans (n,m,a1,a2)
mat_mul (m,n,a2,a1,a3)
mat_vec(m,n,a2,y,rl)
mat_inv (m,a3,a1)
mat_vec(m,n,a1,rl,c)
end def

def vec_uit (n$,n,v(),x,y,l,d)
ys=y-30
if n$<>"" then
  draw text n$ at x,y
  ys=ys+30
end if
for i=1 to n
  a$=n2a$(i,3,0) & " " & n2a$(v(i),l,d)
  draw text a$ at x-12,ys+25*i
next i
end def

def mat_trans (n,m,mat(,),matt(,))
for i=1 to n
  for j=1 to m
    matt(j,i)=mat(i,j)
  next j
next i
end def

def mat_mul (n,m,mata(,),matb(,),matc(,))
for i=1 to n
  for j=1 to n
    tot=0
    for k=1 to m
      tot=tot+mata(i,k)*matb(k,j)
    next k
  matc(i,j)=tot
  next j
next i
end def

def mat_vec (n,m,mat(,),vin(),vuit())
for i=1 to n
  tot=0
  for j=1 to m
    tot=tot+mat( i,j)*vin(j)
  next j
  vuit(i)=tot
next i
end def

def mat_inv (nvar,a(,),ainv(,))
dim w(25,50)                      
for i=1 to nvar                 
  for j=1 to nvar ! w(i,j)=a(i,j) ! w(i,j+nvar)=0  ! next j
  w(i,i+nvar)=1
  next i
for piv=1 to nvar
  fac=w(piv,piv)
  for j=piv to piv+nvar ! w(piv,j)=w(piv,j)/fac ! next j
  for i=1 to nvar
    if i<>piv then
      fac=w(i,piv)
      for j=piv to piv+nvar ! w(i,j)=w(i,j)-fac*w(piv,j) ! next j
      endif
    next i
  next piv
for i=1 to nvar
  for j=1 to nvar ! ainv(i,j)=w(i,j+nvar) ! next j
  next i
end def

def xy_grid (xtop,ytop,xbot,ybot,alfa,beta,xpix,ypix)
if xpix=0 or ypix=0 then
  fill rect xtop-6,ytop-6 to xbot+6,ybot+6
  return
end if
xc=xtop+alfa*(xbot-xtop)
yc=ytop+beta*(ybot-ytop)
graphics lock
draw color 0,0,1
draw size 2
draw line xtop,yc to xbot,yc
draw line xc,ytop to xc,ybot
draw color 0,0,0
draw size 3
draw rect xtop-4,ytop-4 to xbot+4,ybot+4
draw size 1
draw alpha 0.3
for x=xc+xpix to xbot step xpix
  draw line x,ytop to x,ybot
next x
for x=xc-xpix to xtop step -xpix
  draw line x,ytop to x,ybot
next x
for y=yc+ypix to ybot step ypix
  draw line xtop,y to xbot,y
next y
for y=yc-ypix to ytop step -ypix
  draw line xtop,y to xbot,y
next y
draw alpha 1
graphics unlock
end def

def poly (n,coef(),x)
res=coef(n+1)
for i=n to 1 step-1
  res=res*x+coef(i)
next i
poly=res
end def

def n2a$(num,lang,dec)
b$="               "
fac=10^dec
num$=floor(fac*num)/fac
tot=lang-len(num$)
if tot<1 then tot=1
rem test ("", tot)
a$=substr$(b$,1,tot) & num$
n2a$=a$
end def

def stat (n,a())
min=a(1)
max=a(1)
som=a(1)
for i=2 to n
  if min>a(i) then min=a(i)
  if max<a(i) then max=a(i)
  som=som+a(i)
next i
a(n+1)=min
a(n+2)=max
a(n+3)=som
a(n+4)=som/n
end def

def test (in$,in)
fill rect 30,900 to 400,940
a$="test " & in$ & in
draw text a$ at 30,900
button "cont" title "doorgaan ?" at 440,890
loopje:
if button_pressed("cont")=0 then goto loopje
button "cont" delete
fill rect 30,900 to 400,940
end def


def box(xlu,ylu,width,height,xm,ym,scale)
' xlu and ylu = left upper corner in pixel coordinates
' width and height = size of graphics box in pixels
' xm and ym = relative position of graph origin
' scale = number of pixels per graph unit
'
org_x=xlu+xm*width ! org_y=ylu+(1-ym)*height
xrl=xlu+width ! yrl=ylu+height
refresh off
draw size 4
draw rect xlu-4,ylu-4 to xrl+4,yrl+4
draw size 1 ! draw alpha .2
for x=org_x to xrl step scale ! draw line x,ylu to x,yrl ! next x
for x=org_x to xlu step -scale ! draw line x,ylu to x,yrl ! next x
for y=org_y to yrl step scale ! draw line xlu,y to xrl,y ! next y
for y=org_y to ylu step -scale ! draw line xlu,y to xrl,y ! next y
draw color 0,0,1 ! draw size 2 ! draw alpha 1
draw line xlu,org_y to xrl,org_y
draw line org_x,ylu to org_x,yrl
refresh
end def

def draw_it(xs,ys,xe,ye,clip)
' all 4 coordinates in graph units
'
x1=box.org_x+box.scale*xs ! y1=box.org_y-box.scale*ys
x2=box.org_x+box.scale*xe ! y2=box.org_y-box.scale*ye
if not clip then draw
dx=x2-x1 ! dy=y2-y1
if dx<>0 then p=dy/dx else p=sign(dy)*1000
q=y1-x1*p
if x1<box.xlu then            ' clip left side of the box
  if x2<box.xlu then return
  x1=box.xlu ! y1=p*x1+q
  end if
if x2<box.xlu then ! x2=box.xlu ! y2=p*x2+q ! end if
if x2>box.xrl then             ' clip right side
  if x1>box.xrl then return
  x2=box.xrl ! y2=p*x2+q
  end if
if x1>box.xrl then ! x1=box.xrl ! y1=p*x1+q ! end if
if y1<box.ylu then             ' clip upper side
  if y2<box.ylu then return
  y1=box.ylu ! x1=(y1-q)/p 
  end if
if y2<box.ylu then ! y2=box.ylu ! x2=(y2-q)/p ! end if
if y2>box.yrl then             ' clip lower side
  if y1>box.yrl then return
  y2=box.yrl ! x2=(y2-q)/p
  end if
if y1>box.yrl then ! y1=box.yrl ! x1=(y1-q)/p ! end if
draw:
draw line x1,y1 to x2,y2
end def

[Mr. Kibernetik]

Nice program!

[Lauderdale]

Thanks I appreciate it

[Lauderdale]

Ok so I put this into my app, and all I get are a bunch of numbers on my screen. Do I have to edit the x and y values within the program? What do coefficients 1, 2, and 3 mean? Is that for a, b and c of y = ax^2+bx+c?

[Henko]

Ok so I put this into my app, and all I get are a bunch of numbers on my screen. Do I have to edit the x and y values within the program? What do coefficients 1, 2, and 3 mean? Is that for a, b and c of y = ax^2+bx+c?

I'll adapt the program a little. Don't be surprised or upset that there's very little documentation, because it is a program for my own use and not a commercial program

Indeed the coordinations of the points are presently in a read/data construct. Ill change that. Moreover i will restrict the polynomial degree to 1, a straight line and remove the graphics and other plot functions that are not used in this program. So you will have a to-the-point program for lineair regression. Try to program then the r^2 correlation coefficient yourself, i'll jump in if it doesn't work out.

But all this tomorrow, it is sleeping time over here.

[Henko]

Hi, Lauderdale. Here's your simple regression program.
The coordinates of the point are now stuffed in a separate, editable datafile, named "xypoints". This file is created at the very first time you run the program and you can edit it using the SmartBasic editor. Presently it will contain nine test points. If you modify the file, don't forget to put the number of points at the beginning of the file.
Your r^2 factor has yet to be calculated.

Code: Select all

' lineair regression 2 
' fits a straight line through a number of points in R2
' (look into "matlib" for nonlineair and 3D fitting)
' coded by Henko, feb/2015
'
nmax=100
dim x(nmax),y(nmax), ycalc(nmax), err(nmax)
graphics ! graphics clear ! draw color 0,0,0 ! option base 1
button "start" text "start calculation" at 250,30 size 200,50
once_only()

loop1:                            ' start and load the data file
  if button_pressed("start") then
  load_points() ! button "start" hide ! goto cont1
  end if
  goto loop1

cont1:
sx=0 ! sx2=0 ! sy=0 ! sxy=0       ' perform the lineair fit
for i=1 to np
  sx+=x(i) ! sx2+=x(i)*x(i) ! sy+=y(i) ! sxy+=x(i)*y(i)
  next i
det=np*sx2-sx*sx
yo=(sx2*sy-sx*sxy)/det ! dir=(np*sxy-sx*sy)/det

for i=1 to np                      ' prepare and print output
  ycalc(i)=yo+dir*x(i) ! err(i)=y(i)-ycalc(i)
  next i
t$="fitting line : y = " & yo & " + " & dir & " * " & "x"
draw text t$ at 100,20
vec_out(" x-coords",np,x,10,60,6,1)
vec_out(" y-coords",np,y,160,60,6,1)
vec_out("y-calculated",np,ycalc,310,60,8,3)
vec_out("   errors",np,err,480,60,8,3)
draw text "Lauderdale r^2 factor = " & "?" at 200,60+30*np
end

def load_points()
  filename$="xypoints"
  file filename$ input .np
  for i=1 to .np ! file filename$ input .x(i),.y(i) ! next i
end def

' print a vector at position x,y with number length l and precision d
'
def vec_out (n$,n,v(),x,y,l,d)
ys=y-30
if n$<>"" then ! draw text n$ at x,y ! ys=ys+30 ! end if
for i=1 to n
  a$=n2a$(i,3,0) & " " & n2a$(v(i),l,d)
  draw text a$ at x-12,ys+25*i
next i
end def

' formatting floats for graphics mode
'
def n2a$(num,lang,dec)
b$="               "
fac=10^dec ! num$=floor(fac*num)/fac ! tot=lang-len(num$)
if tot<1 then tot=1
a$=substr$(b$,1,tot) & num$ ! return a$
end def

def once_only()
  f$="xypoints" ! if file_exists(f$) then return
  np=9 ! file f$ print np
  for i=1 to np ! read xt,yt ! file f$ print xt,yt ! next i
  data -2,1,2,4,3,3.5,7,6,10,6,14,7,15,11,18,11,21,9
end def

[Lauderdale]

Code: Select all

  print ""
print ""
print "   R squared calculator"
print ".     by 'Lauderdale'"
print ""
varr = 1
button "c" text "touch to continue" at 10,450 size 300,50
while varr = 1
if button_pressed ("c")=1 then varr = 2
end while
button "c" delete
refresh
print ""
print "First you will input the data points (x,y).  After, you input the coefficients of your linear regression equation y = ax + b.  You will enter a and b in the appropriate boxes."
varr = 2
button "d" text "touch to start" at 10,450 size 300,50
while varr = 2
if button_pressed ("d")=1 then varr = 3
end while
button "d" delete

input "# of data points ":data
totaly = 0
if data < 5 then
print ""
print ""
print ""
print "Sorry but you need nore data points"
goto end
endif
dim points(data,2)
for k = 0 to data-1
j = k+1
input "x"&j&"":xx
points(k,0) = xx
input "y"&j&"":yy
points(k,1) = yy
totaly = totaly + yy
next k
input "slope coefficient a":a
input "y-intercept b":b

avrg = totaly/j
ssa = 0
ssb = 0
for m = 0 to data-1
ssa = (points(m,1)-avrg)^2 + ssa
ssb = ((a*points(m,0)+b)-points(m,1))^2 + ssb
next m


if ssa = 0 then 
print "only space goats can divide by zero."
else
rsq = 1 - (ssb/ssa)
print ""
print ""
print ""
print "r^2 = "&rsq&""
endif
end:  

Here is the r^2 program, I'll give yours a try In a little bit

[Henko]

if ssa = 0 then
print "only space goats can divide by zero."
else
rsq = 1 - (ssb/ssa)

If ssa=0 then rsq=1 else rsq=1-ssb/ssa

Another point: i prefer to see the result in graphical presentation:
Imagine 2 "clouds" of points. One scattered around a straight line, and the other scattered about a parabolic section, but both having the same r^2 value of .85 (example). Looking at the r^2 values, the only conclusion is a "rather good" correlation in both cases.
But looking at the graphic representations, you can see the "rather good" correlation, but in the second case also that a much better correlation could be obtained by a 2nd degree fit.

[Lauderdale]

Ok thanks for the help, I will look into higher degree fits as well.