% @(#)LN_BasicStats.m	1.2 02/08/19  

load Pain.dat

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

clf;axis([300 900 -1/2 10.5]);box on
for i=1:12
  text(Engy(Subj==i),Y(Subj==i),num2str(i))
end
title('Subject Pain Rating vs Laser Energy')
fullpage('h')
print

%
%  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('h')
print
%
%  Multiple intercepts, same slope
%


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    
   
%
% Equivalently, could use the F contrast...
%

c2F = ...
   [ 0 0 1 0 0 0 0 0 0 0 0 0 0 0 ;
     0 0 0 1 0 0 0 0 0 0 0 0 0 0 ;
     0 0 0 0 1 0 0 0 0 0 0 0 0 0 ;
     0 0 0 0 0 1 0 0 0 0 0 0 0 0 ;
     0 0 0 0 0 0 1 0 0 0 0 0 0 0 ;
     0 0 0 0 0 0 0 1 0 0 0 0 0 0 ;
     0 0 0 0 0 0 0 0 1 0 0 0 0 0 ;
     0 0 0 0 0 0 0 0 0 1 0 0 0 0 ;
     0 0 0 0 0 0 0 0 0 0 1 0 0 0 ;
     0 0 0 0 0 0 0 0 0 0 0 1 0 0 ;
     0 0 0 0 0 0 0 0 0 0 0 0 1 0 ;
     0 0 0 0 0 0 0 0 0 0 0 0 0 1 ];

% ... but then the math is even worse


%
%  Multiple intercepts, different slopes
%


X3 = [Cons Engy spm_DesMtx(Subj,'-','Subj') ...
      diag(Engy)*spm_DesMtx(Subj,'-','Subj')];
X3d = spm_DesMtx('sca',Cons,'Const',Engy,'Energy');
image(([X3d X3(:,3:3+11) diag(X3d(:,2))*X3(:,3:3+11) ]+1)*32);colormap gray
print
p3 = rank(X3);

B3 = pinv(X3)*Y;
Y3 = X3*B3;


sigs3  = (Y-Y3)'*(Y-Y3)/(n-p3);
  %  sigs2 = 2.6158 
sig3 = sqrt(sigs3)
  %   sig2 = 1.6174 

c3 = [0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ];

c3*B3 
  %  0.0215 

cV3 = c3*pinv(X3)*pinv(X3)'*c3' * sigs3
  %   4.5689e-06 


%
%  T test
%
t3 = c3*B3 / sqrt(cV3)
  %  10.0624 
  

1-spm_Tcdf(t3,n-p3)
   %  0   


  

%
%  Do we need the different slopes?
%

RSS2 = (Y-Y2)'*(Y-Y2);
RSS3 = (Y-Y3)'*(Y-Y3);

F23 = (RSS2-RSS3)/((n-p2)-(n-p3))/ sigs3
   %  2.6217 
   
1-spm_Fcdf(F23,(n-p2)-(n-p3),n-p3)
   %  0.0035