Support Forum / Форум поддержки

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

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

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: