www.neuron.yale.edu

Code: Select all

/*
Fits y = mx + b to data by linear regression.
Expects x and y data in Vectors, returns m and b. 
$o1  x data
$o2  y data
$&3  m
$&4  b
*/

proc linfit() { local num, xmean, ymean, sxy, s1sq  localobj xvec, yvec
  num = $o1.size()
  if ($o1.size() < 2) {
    print "Data must contain two or more points."
    quit()
  }

  xmean = $o1.mean
  ymean = $o2.mean

  xvec = $o1.c.sub(xmean)
  yvec = $o2.c.sub(ymean)

  sxy = xvec.c.mul(yvec).sum/(num-1)
  s1sq = xvec.sumsq()/(num-1)

  $&3 = sxy/s1sq
  $&4 = -$&3*xmean + ymean
}

Code: Select all

load_file("nrngui.hoc")
load_file("linfit.hoc")

// generate test data
NUMPTS = 10
objref xdat, ydat
xdat = new Vector(NUMPTS)
ydat = new Vector(NUMPTS)
objref r
r = new Random()
r.uniform(-1, 1)
xdat.setrand(r)
ydat.setrand(r)

/*
rearrange data so x values are monotonically increasing
$o1  xdata
$o2  ydata
*/

proc sortxy() { local ii  localobj xvec, yvec, inx
  inx = $o1.sortindex()
  xvec = new Vector()
  yvec = new Vector()
  xvec.index($o1, inx)
  yvec.index($o2, inx)
  $o1 = xvec
  $o2 = yvec
}

sortxy(xdat, ydat)

// plot data, letting graph axes extend slightly beyond range
objref g
g = new Graph(0)
g.size(-1.1,1.1, -1.1,1.1)
g.view(-1.1, -1.1, 2.2, 2.2, 212, 109, 300.48, 200.32)
ydat.mark(g, xdat, "O", 6)

// calculate slope and intercept
m = b = 0
linfit(xdat, ydat, &m, &b)
print "slope ", m, " y intercept ", b

// plot regression line

func rline() {
  return m*$1 + b
}

objref yfit
yfit = xdat.c.apply("rline")
yfit.mark(g, xdat, "+", 12, 2, 1)

g.color(2)
g.brush(2)
g.beginline()
g.line(-1.1, rline(-1.1))
g.line(1.1, rline(1.1))
g.flush()