next up previous index
Next: Generating Normal(0,1) Random Variables Up: Octave Programs Previous: Octave Programs   Index

Click for printer friendely version of this HowTo


Univariate General Linear Models

function beta=general_linear(x,y,c)

# y is a column (nx1) vector of observations
# x is the design matrix (nxm)
# c is the contrast matrix
# NOTE: it is assumed that c is set up so that theta = 0

  n = max(size(x)); # number of observations
  p = min(size(x)); # number of parameters
  num_tests = min(size(c)); # number of tests

  beta = inv(x'*x)*x'*y # estimate the parameters

  # test model
  numerator = beta'*c'*inv((c*inv(x'*x)*c'))*c*beta / num_tests;
  denominator = (y'*y - y'*x*beta)/(n-p); # this is also known as MSE
  f_test = numerator / denominator
  p_value = 1 - f_cdf(f_test, num_tests, n-p) # calculate p-value

  if (num_tests == 1)
    t_test = sqrt(f_test)
  end

  # generate graphical output
  MIN_F_DIST = 0;
  MAX_F_DIST = 20;
  NUM_STEPS = 100;

  graph_x = linspace(MIN_F_DIST, MAX_F_DIST, NUM_STEPS);
  graph_y = f_pdf(graph_x, num_tests, n-p);

  description = sprintf("g;F(%d, %d)-Dist;", num_tests, n-p);
  plot(graph_x, graph_y, description, f_test, 0, "r*;Your Data;")

endfunction



Frank Starmer 2004-05-19
>