% sfran.m
y = [ 22.6 20.5 20.8; 22.6 21.2 20.5; ...
      17.3 16.2 16.6; 21.4 23.7 23.2; ...
      20.9 22.2 22.6; 14.5 10.5 12.3; ...
      20.8 19.1 21.3; 17.4 18.6 18.6; ...
      25.1 24.8 24.9; 14.9 16.3 16.6];
f = @(theta) -sfran_loglike(theta(1),theta(2),theta(3),y);
yhat = mean(y,2);
theta0 = [mean(yhat) var(yhat) mean(var(y,0,2)) ];
thetahat = fminsearch(f,theta0)