Ntrials=10;
Nsub=12;
N=Ntrials*Nsub;

pop_effect=10;
sigma_b_2=1;
ymin=pop_effect-3*sqrt(sigma_b_2); 
ymax=pop_effect+3*sqrt(sigma_b_2);

% Within subject variance
lambda=rand(Nsub,1);
%lambda(5)=10;

% Generate design matrices
X1=kron(eye(Nsub),ones(Ntrials,1));
X2=ones(Nsub,1);

% Generate first-level covariance components
C1=zeros(N,N);
for n=1:Nsub,
    q=zeros(Ntrials*Nsub,1);
    start=(n-1)*Ntrials+1;
    stop=start+Ntrials-1;
    q(start:stop)=1;
    Q{n}=diag(q);
    C1=C1+lambda(n)*Q{n};
end

% Second level
lambda(n+1)=sigma_b_2;
C2=lambda(n+1)*eye(Nsub);
Q{Nsub+1}=X1*eye(Nsub)*X1';

% Error covariance for collapsed model
C=C1+X1*C2*X1';

% Generate and plot data

sub=randn(Nsub,1)*sqrt(lambda(Nsub+1))+pop_effect;
y=[];
figure
for n=1:Nsub,
    trials=randn(Ntrials,1)*sqrt(lambda(n))+sub(n);
    subplot(1,Nsub,n);
    plot(ones(Ntrials,1),trials,'x');
    axis([0 2 ymin ymax]);
    set(gca,'XTickLabel',[]);
    if n>1
        set(gca,'YTickLabel',[]);
    end
    xlabel(int2str(n));
    y=[y;trials];
end


% Plot design matrix and covariance components
figure
subplot(2,2,1);
imagesc(X1);
colormap gray
title('X_1');

subplot(2,2,2);
imagesc(C1);
colormap gray
title('C_1');

subplot(2,2,3);
imagesc(X2);
colormap gray
title('X_2');

subplot(2,2,4);
imagesc(C2);
colormap gray
title('C_2');

% ReML estimation
X=X1*X2;
YY=y*y';
[RC,Rh,RPh]=spm_reml(YY,X,Q,1);

% Display results
figure
true_var=lambda(1:Nsub);
est_var=Rh(1:Nsub);
mv=max([true_var; est_var]);
plot(true_var,est_var,'rx');
hold on
plot([0 mv],[0 mv],'-');
title('Error variances');
xlabel('True');
ylabel('ReML estimate');