home *** CD-ROM | disk | FTP | other *** search
/ Shareware Overload / ShartewareOverload.cdr / progm / grafx.zip / SFIT.FOR < prev    next >
Text File  |  1989-03-21  |  2KB  |  77 lines
  1. c    program main
  2. c    call gset                            !  enter graphic mode
  3. c    call sfit
  4. c    call aset                            !  return to alpha mode
  5. c    end    
  6. c    -----------------------------------------------------------------------
  7.     subroutine sfit
  8. c
  9. c    ...fit a straight line to data
  10. c
  11. c    (C) Copyright 1988, 1989 by Jim Farrell      All Rights Reserved.
  12. c
  13.     parameter (n=100)
  14.     real x(n),y(n)
  15. c
  16.     fn(x)=30.0+2.0/3.0*x                !  linear function
  17. c
  18.     ikr=1.0
  19.     xmax=real(n)/2.0
  20.     xmin=-xmax
  21.     ymin=fn(xmin)
  22.     ymax=fn(xmax)
  23.     dx=(xmax-xmin)/(n-1)
  24.     dy=(ymax-ymin)/5.0
  25.     call putstr(10,20,' DEMONSTRATION MULTIPLE PLOT OVERLAYS')
  26.     call putstr(12,20,'       FOR A LEAST SQUARES FIT')
  27.     call pause('WHEN READY ')        !  display message and pause
  28.     call gcls
  29.     call mulplt0(xmin,ymin-dy,xmax,ymax+dy)    !  set max/mins
  30.     do 200 i=1,n                        !  compute line from function
  31.         x(i)=xmin+dx*real(i-1)
  32.         y(i)=fn(x(i))
  33. 200    continue
  34.     call mulplt(x,y,n,4)                !  plot as connected lines
  35.     call pause('REVIEW PLOT OF PURE LINE')
  36.     do 300 i=1,n                        !  corrupt data with noise & bias
  37.         x(i)=xmin+(xmax-xmin)*rand(ikr)
  38.         y(i)=fn(x(i))+dy*sign(rand(ikr),rand(ikr)-0.3)
  39. 300    continue
  40.     call mulplt(x,y,n,0)                !  plot as points
  41.     call pause('REVIEW PLOT OF CORRUPTED DATA')
  42.     call fits(x,y,n,a,b)                !  fit slope & intercept
  43.     do 400 i=1,n                        !  compute line based on fit
  44.         x(i)=xmin+dx*real(i-1)
  45.         y(i)=a+b*x(i)
  46. 400    continue
  47.     call mulplt(x,y,n,4)                !  plot as connected lines
  48.     call pause('REVIEW PLOT OF COMPUTED LINE')
  49.     return
  50.     end
  51. c    -----------------------------------------------------------------------
  52.     subroutine fits(x,y,ndata,a,b)
  53. c
  54. c    ...fit a straight line to data (i.e. determine slope and intercept)
  55. c
  56.     real x(*),y(*)
  57. c
  58.     sx=0.0
  59.     sy=0.0
  60.     do 130 i=1,ndata
  61.         sx=sx+x(i)
  62.         sy=sy+y(i)
  63. 130    continue
  64.     sxoss=sx/real(ndata)
  65.     st2=0.0
  66.     b=0.0
  67.     do 140 i=1,ndata
  68.         t=x(i)-sxoss
  69.         st2=st2+t**2
  70.         b=b+t*y(i)
  71. 140    continue
  72.     b=b/st2
  73.     ss=real(ndata)
  74.     a=(sy-sx*b)/ss
  75.     return
  76.     end
  77.