% This script fits a Gaussian white noise model to each patient's data and
% saves the resulting AIC, AICc, and log likelihood values to a mat file.

% In version 2: restrict to two weeks! Adjust for patients where the two
% weeks of interest doesn't start at start of data set.

clear all;
load PatientData_final.mat;
allpatients = who();

Npats = length(allpatients); % Number of patients
AICvec = zeros(Npats, 1);  % will hold AIC output
AICcvec = zeros(Npats, 1);  % will hold AIC output
loglikevec = zeros(Npats, 1);  % will hold AIC output
for i=1:Npats
    % create strings of command I actually want to use:
    cmd1 = sprintf('%s(:,1);', allpatients{i});
    cmd2 = sprintf('%s(:,10);', allpatients{i});
    % evaluate strings so the pain reports of patient i go into 'painvec'
    % and time values go into 'timevec'
    timevec = eval(cmd1);
    painvec = eval(cmd2);
    % Strip out nan's where data is missing:
    badinds = find(isnan(painvec));
    timevec(badinds)=[];
    painvec(badinds)=[];
    % Exclude everything after 336 hours (two weeks):
    badinds = find(timevec>336);
    tmptime = timevec;
    tmppain = painvec;
    tmptime(badinds) = [];
    tmppain(badinds) = [];
    % Verify that there is more than one report left.  If so, done, but if
    % not, need to start clock at time of _second_ pain report.
    if length(tmppain)>3
        timevec = tmptime;
        painvec = tmppain;
    elseif length(tmppain)==1
        timevec = timevec(2:end);
        painvec = painvec(2:end);
        badinds = find(timevec>timevec(1)+336);
        timevec(badinds) = [];
        painvec(badinds) = [];
    elseif length(tmppain)==2
        timevec = timevec(3:end);
        painvec = painvec(3:end);
        badinds = find(timevec>timevec(1)+336);
        timevec(badinds) = [];
        painvec(badinds) = [];
    elseif length(tmppain)==3
        timevec = timevec(4:end);
        painvec = painvec(4:end);
        badinds = find(timevec>timevec(1)+336);
        timevec(badinds) = [];
        painvec(badinds) = [];
    end
        
    % get mean and std. dev. of pain reports (ignoring nans)
    mu = mean(painvec);
    sigma = std(painvec); 

    % Now the task is to compute the probability of the pain reports 
    % assuming a white noise model with the same mean and std. dev..
    % Start by defining a gaussian with the given mu, sigma:
    gaussfunc = @(x)(exp(-(mu-x).^2./(2*sigma^2))./sqrt(2*pi*sigma^2));
    % The likelihood (prob.) of a given pain report is just this Gaussian evaluated 
    % at the given pain level P.  To get the likelihood of _all_ pain 
    % reports, we multiply all of those likelihood probabilities.  
    % For computational purposes, it's
    % numerically better to compute the log of that (log likehlihood),
    % which means adding up all the log probabilities.
    loggaussfunc = @(x)(-log(2*pi)/2-log(sigma)-(mu-x).^2./(2*sigma^2));
    loglike = sum(loggaussfunc(painvec));
    % The AIC is defined as 2k - 2 log L, where k is the number of
    % parameters and L is the likelihood.  Here k=2 (from mu & sigma):
    k = 2; % number of parameters in model
    n = length(painvec); % number of data points for this patient
    AIC = 2*k - 2*loglike;
    if (n>k-1)
        AICc = AIC + 2*k*(k+1)/(n-k-1);
    else
        AICc = nan;
    end
    fprintf(1, 'AIC for %s: %2.2f (log likelihood %2.2f)\n', ...
        allpatients{i}, AIC, loglike);
    AICvec(i)=AIC;
    AICcvec(i)=AICc;
    loglikevec(i)=loglike;
    
    % For debugging Gaussiand and log Gaussian functions:
    %tmp = (-5:0.1:15);
    %plot(tmp, gaussfunc(tmp), 'r'); hold on;
    %plot(tmp, exp(loggaussfunc(tmp)), ':b');
    %sum(gaussfunc(tmp)*0.1)
end

save('whitenoiseAICvals_twoweeks.mat', 'AICvec', 'AICcvec', 'loglikevec', ...
    'allpatients');