load Pain.dat

n    = size(Pain,1);
Y    = Pain(:,1);   % Subjective rating
Engy = Pain(:,2);
Cons = ones(n,1);


plot(Engy,Y,'x');set(gca,'ylim',[-1/2 10.5])
title('Subject Pain Rating vs Laser Energy')
fullpage('v')
print

delete(gcf)

%
%  Simple correlation model
%
%

image((spm_DesMtx('sca',Cons,'Const',Engy,'Energy')+1)*32);colormap gray
print

X1 = [Cons Engy];
p1 = 2

B1 = inv(X1'*X1)*X1'*Y;
Y1 = X1*B1;

sigs1  = (Y-Y1)'*(Y-Y1)/(n-p1);
  %  sigs1 =  3.5277  
sig1 = sqrt(sigs1)
  %   sig1 =  1.8782 

c1 = [0 1];

c1*B1 
  %  0.0080

cV1 = c1*inv(X1'*X1)*c1' * sigs1
  % 1.4875e-06      


%
%  T test
%
t1 = c1*B1 / sqrt(cV1)
  %  6.5489 
  
  
% P(T > 6.5489)  = 1 - P( T <=  6.5489)  = F( 6.5489 )

1-spm_Tcdf(t1,n-p1)
   %  1.3966e-10   


plot(Engy,Y,'x');set(gca,'ylim',[-1/2 10.5])
line([min(Engy) max(Engy)],B1(1)+[min(Engy) max(Engy)]*B1(2))
title('Subject Pain Rating vs Laser Energy')
fullpage('v')
print


%
%  Multiple intercepts, same slope
%

Subj = Pain(:,3);

X2 = [Cons Engy spm_DesMtx(Subj,'-','Subj')];
X2d = spm_DesMtx('sca',Cons,'Const',Engy,'Energy');
image(([X2d X2(:,3:end)]+1)*32);colormap gray
print
p2 = rank(X2);

B2 = inv(X2'*X2)*X2'*Y;
% Get message:  Warning: Matrix is close to singular or badly scaled.    
% because X2 is not uniquely determined
% Use pinv instead

B2 = pinv(X2)*Y;
Y2 = X2*B2;


sigs2  = (Y-Y2)'*(Y-Y2)/(n-p2);
  %  sigs2 =  2.7906  
sig2 = sqrt(sigs2)
  %   sig2 =  1.6705 

c2 = [0 1 0 0 0 0 0 0 0 0 0 0 0 0 ];

c2*B2 
  %  0.0176 

cV2 = c2*pinv(X2)*pinv(X2)'*c2' * sigs2
  %  2.4210e-06   


%
%  T test
%
t2 = c2*B2 / sqrt(cV2)
  %  11.2850   
  

1-spm_Tcdf(t2,n-p2)
   %  0   


  

%
%  Do we need the extra intercepts?
%

RSS1 = (Y-Y1)'*(Y-Y1);
RSS2 = (Y-Y2)'*(Y-Y2);

F12 = (RSS1-RSS2)/((n-p1)-(n-p2))/ sigs2
   %  7.6755  
   
1-spm_Fcdf(F12,(n-p1)-(n-p2),n-p2)
   %  1.6425e-11