%% Experimental Design for a Multivariate Sensitivity Test
% This file 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
%
% At each step, the algorithm will suggest a location for the 
% next observations, such as: Recommendation:   -1.0000    0.4564   -1.0000 1.0000
% You can either press 'enter' to accept, or enter a different location
% such as:  "enter"=accept, or enter value:    [0.2 0.25 0 0]
%
% 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

% 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=40; % 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;
X0=[]; % This is the matrix with the experimental steps already done
RResult=[];
ss=0;

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

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;

      % Recommend a value for the next observations
      fprintf('Recommendation:');
      disp(X0Best);
    %text1=sprintf('Recommendation: %g  ; "enter"=accept, or enter value: ',X0Best);
    XChosen=input(' "enter"=accept, or enter value:'); % input where to take the next observations
    if isempty(XChosen), XChosen=X0Best; end;

    RResultTemp=[];
    X0=[X0 ; XChosen];
    ss=ss+1;
    F0=x2fx(X0);

    while isempty(RResultTemp) || ~isnumeric(RResultTemp) || (RResultTemp~=0 & RResultTemp~=1)
        RResultTemp=input('Enter the result (0 /1): ');  % Input the result of the observation (0=no response, 1=response)
    end;
    RResult(end+1)=RResultTemp;

    %% 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






