%% Experimental Design for a Multivariate Sensitivity Test
% This file is a Monte-Carlo examination, extending the file
% MultivariateTest.m
% In each iteration, a "true" setup is chosen at random
% from the prior distribution,
% The algorithm runs automatically,
% and the Local D-Efficiency is recorded.
%
% This file,too, is a source code for the algorithm found
% for multivariate sequential experimental design for GLM
% using the example from section 4.3 of our paper,
% available at http:/www.math.tau.ac.il/~dms/GLM_Design
%
% This script was written by Hovav Dror, Tel-Aviv University, 2006
%
% If you make changes to the script, or if you implement it
% to a (to be) published problem, please inform us:
% dms@post.tau.ac.il ; hovavdror@gmail.com
%
% further details, including papers describing the algorithm, can be
% found at: http://www.math.tau.ac.il/~dms/GLM_Design


%% Input Parameters

MaxIter=1000; % Number of Monter-Carlo iterations
iter=0; % current iteration

% In this example we assume that all variables are codes to the region [-1,1]
m=4; % Number of variables
p=5; % number of coefficiens

PLD=8; % Recommended number of points in a good Local D-optimal design

% Prior Distribution
betar=[ -3  3
    4  10
    5  11
    -6  0
    -2.5  3.5];
NN=10000;
NI=rand(NN,p);
beta=[];
for j=1:p
    beta(:,j)=betar(j,1)*ones(NN,1) + ( betar(j,2)-betar(j,1) ) * NI(:,j);
end;

family='binomial';
link='logit';
model='linear';

L=48; % Length of the experiment (total number of observations)

% more initialization
N1=500; % N1 - number of vectors to calculate the weighted sum of determinants of information matrices
RealMin=realmin;
clo=clock;
Eff=[];

fprintf('Process began at %g:%g:%1.0f \n',clo(4:6));

%% begin iterations
while iter<MaxIter
    iter=iter+1;
    %% True parameters
    NI=rand(1,p);
    Tbeta=[];
    for j=1:p
        Tbeta(:,j)=betar(j,1) + ( betar(j,2)-betar(j,1) ) * NI(:,j);
    end;
    X0=[]; % This is the matrix with the experimental steps already done
    RResult=[];
    ss=0; % ss=size(X0,1);
    Weights=ones(NN,1)/NN;

for Rem=L:-1:1
    fprintf('%g ',ss+1); % the number of the current observation

    % find the Weighted Median (the "best guess" for the parameter vector)
    betaMed=[];
    for i=1:p
        betai=sortrows([beta(:,i) Weights]);
        betai(:,3)=cumsum(betai(:,2));
        betai(:,3)=abs(betai(:,3)-0.5);
        betaMedi=betai(betai(:,3)==min(betai(:,3)),1);
        betaMed(1,i)=median(betaMedi);% in case there was more than one Mu that met the criterion
    end;

    %% Now, use this "best guess" to choose the next point
    X0t=GLMDAUG(X0,betaMed',m,ss+PLD); % PLD=8
    X0t=[median(X0t) ;X0t ];
    BestSoFar=-inf;
    for i=1:length(X0t)
        ThisOne=0;
        X0ti=[X0; X0t(i,:)];
        if ss<(PLD-1), X0ti=[X0ti ; X0t(2:end,:)]; end;
        F0ti=x2fx(X0ti);

        for j=1:N1
            ThisOnej=max(det(InfoMtrxGLM(F0ti,beta(j,:)',family,link) ) ,RealMin);
            ThisOne=ThisOne + Weights(j)*log( ThisOnej );
        end;
        if ThisOne>BestSoFar
            BestSoFar=ThisOne;
            X0Best=X0t(i,:);
        end;
    end;
    X0=[X0 ; X0Best];
    ss=ss+1;
    F0=x2fx(X0);

    % Choose outcome for the observation from the prior
    Texpc=1/(1+exp(-(Tbeta*F0(end,:)')));
    RResult(end+1)= (rand(1)<Texpc); % Observation

    %% Now modify the Posterior distribution:
    RResultM=ones(NN,1)*RResult;
    ex=exp(-F0*beta');
    bprob=(1 ./ (1+ex))';
    mprob=1-bprob;
    LRi = prod( bprob.^RResultM ,2 ) .* prod( mprob.^(1-RResultM) ,2);
    Weights=LRi/sum(LRi);

end; % for Rem


%% Calculate Efficiency
XO=DGLM(Tbeta,m,size(X0,1),model,family,link,1e-4); % optimal design from True beta
Eff(iter)=(det(InfoMtrxGLM(x2fx(X0),Tbeta,family,link))/det(InfoMtrxGLM(x2fx(XO),Tbeta,family,link)))^(1/p);
MinEff=min(Eff);
Q5Eff=quantile(Eff,0.05);
MedEff=median(Eff);
meanEff=mean(Eff);
if iter>1, stdEff=std(Eff); else stdEff=0; end;
fprintf('\n\n');
fprintf('---------------------------------------------------\n');
fprintf('Iter Last  Min     5%%      Med   Mean  std \n');
fprintf(' %2.0f        %2.2f      %2.2f      %2.2f     %2.2f    %2.2f      %2.2f \n',[iter Eff(iter) MinEff Q5Eff MedEff meanEff stdEff]);
fprintf('---------------------------------------------------\n');


end; % while iter<MaxIter

figure()
hist(Eff);
x=axis;
axis([0 1 x(3:4)]);

ElapsedTime=etime(clock,clo);
fprintf('Total: %2.1f hours, or %2.1f minutes per run \n\n',[ElapsedTime/3600 ElapsedTime/iter/60]);





