function [ la,le,ci,lr,pvalue ] = baumgartner2014_likelistat( p,tang,rang,target,response,varargin )
%BAUMGARTNER2014_LIKELISTAT Likelihood statistics for evaluation of model performance
% Usage: [la,le,ci,lr,pvalue] = baumgartner2014_likelistat(p,tang,rang,target,response,varargin)
%
% Input parameters:
% p : pdf matrix
%
% tang : polar angles of possible target angles
%
% rang : polar angles of possible response angles
%
% target : target polar angles of localization test
%
% response : response polar angles of localization test
%
% varargin : Use 'normalize' for normalization of likelihoods.
% 1 corresponds to unitary likelihood. This is the default.
% Use 'original' according to Langendijk et al. (2002).
%
% Output parameters:
% la : actual likelihood
%
% le : expected likelihood
%
% ci : 99% confidence interval for expected likelihood
%
% lr : reference likelihoods
%
%
% BAUMGARTNER2014_LIKELISTAT calculates the likelihood statistics for
% the evaluation of the model performance.
%
% The reference likelihoods are given in the following order:
%
% '1st dim' unimodal (1 gaussian dist.: std=17 deg, mu=0)
%
% '2nd dim' bimodal (2 gaussians: mu1=0, mu2=180)
%
% '3rd dim' trimodal (mu1=0, mu2=90, mu3=180)
%
% '4th dim' unitary
%
% See also: baumgartner2014
%
% Url: http://amtoolbox.org/amt-1.5.0/doc/modelstages/baumgartner2014_likelistat.php
% #StatusDoc: Perfect
% #StatusCode: Perfect
% #Verification: Verified
% #Requirements: SOFA CircStat M-SIGNAL M-Stats O-Statistics
% #Author: Robert Baumgartner (2014), Acoustics Research Institute, Vienna, Austria
% This file is licensed unter the GNU General Public License (GPL) either
% version 3 of the license, or any later version as published by the Free Software
% Foundation. Details of the GPLv3 can be found in the AMT directory "licences" and
% at <https://www.gnu.org/licenses/gpl-3.0.html>.
% You can redistribute this file and/or modify it under the terms of the GPLv3.
% This file is distributed without any warranty; without even the implied warranty
% of merchantability or fitness for a particular purpose.
definput.flags.normlikeli = {'normalize','original'};
definput.flags.plot = {'no_plot','plot'};
definput.flags.jackknife = {'','jackknife'};
definput.keyvals.runs = 100;
[flags,kv]=ltfatarghelper({'runs'},definput,varargin);
runs = kv.runs;
if iscell(p)
pa = cell(length(p),1); % predicted probabilities at actual positions
pe = cell(length(p),1); % predicted probabilities at expected positions
for ll = 1:length(p)
% Exclude target-response pairs not defined in prediction matrix
tvalid = target{ll} < max(tang{ll}) & target{ll} > min(tang{ll});
rvalid = response{ll} < max(rang{ll}) & response{ll} > min(rang{ll});
idvalid = tvalid & rvalid;
targets = target{ll}(idvalid);
responses = response{ll}(idvalid);
nt(ll) = sum(idvalid);
pa{ll} = interp2(tang{ll},rang{ll},p{ll},targets,responses);
pa{ll} = pa{ll}';
% Define target angle boundaries necessary for computation of L_expect
tangbound = tang{ll}(:)+0.5*diff([tang{ll}(1)-diff(tang{ll}(1:2));tang{ll}(:)]);
% Indices of target positions
post=zeros(size(targets));
for ii = 1:length(post)
post(ii) = find(tangbound>=targets(ii),1);
end
% Generate response angles
pe{ll} = zeros(runs,length(post));
for ind=1:runs
posr=zeros(size(post));
for jj = 1:length(post)
posr(jj) = local_discreteinvrnd(p{ll}(:,post(jj)),1);
pe{ll}(ind,jj)=p{ll}(posr(jj),post(jj));
end
end
end
Nt = sum(nt);
% Actual likelihood
la=-2*sum(log([pa{:}]+eps));
% Expected likelihood
lex = -2*sum(log([pe{:}]+eps),2);
le=mean(lex);
err=2.58*std(lex); % standard deviation
% err=2.626*std(lex); % 99%-quantile of Studen-t dist.
ci=[le-err; le+err]; % 99% confidence interval
% ci = [min(lex) max(lex)];
% Probability that L_actual occurs as L_expected
[muhat,sigmahat] = normfit(lex);
if la > muhat
pvalue = 2*(1-normcdf(la,muhat,sigmahat));
else
pvalue = 2*normcdf(la,muhat,sigmahat);
end
% Reference Likelihoods ------------------------------------------------
s = 17; % std in degrees
idmed = ceil(length(p)/2);
% unimodal
pdfuni = exp(-0.5 * (rang{idmed}./s).^2); % normpdf
pdfuni = pdfuni /sum(pdfuni);
iduni = local_discreteinvrnd(pdfuni,Nt,runs);
lr(1,1) = mean(-2*sum(log(pdfuni(iduni))));
% unitary
lr(4,1) = -2*Nt*log(1/length(rang{1}));
% normalization of likelihoods, 1 corresponds to unitary
if flags.do_normalize
la = la/lr(4);
le = le/lr(4);
ci = ci/lr(4);
lr = lr/lr(4);
end
else
% target = round(target); % !!!!!!!!!!!!!!!!
tangbound = tang(:)+0.5*diff([tang(1)-diff(tang(1:2));tang(:)]); % !!!!!!!!!
% Actual Likelihood ----------------------------------------------------
% Exclude undefined target angles
tvalid = target < max(tang) & target > min(tang);
rvalid = response < max(rang) & response > min(rang);
idvalid = tvalid & rvalid;
target = target(idvalid);
response = response(idvalid);
nt = sum(idvalid);
pa = interp2(tang,rang,p,target,response);
% Evaluate actual likelihood
la=-2*sum(log(pa+eps)); % actual likelihood
% Jackknife actual likelihood
if flags.do_jackknife
jackdiv = 2;
njack = nt/jackdiv; % resample
lax = zeros(runs,1);
idjack = zeros(njack,runs);
for ind=1:runs
idjack(:,ind) = randi(nt,njack,1);
pax = interp2(tang,rang,p,target(idjack(:,ind)),response(idjack(:,ind)));
lax(ind)=-2*sum(log(pax+eps));
end
end
% Expected Likelihood --------------------------------------------------
post=zeros(size(target)); % indices of target positions
for ii = 1:nt
post(ii) = find(tangbound>=target(ii),1);
end
lex=zeros(runs,1);
for ind=1:runs
if flags.do_jackknife
postsel = post(idjack(:,ind));
else
postsel = post;
end
posr=zeros(size(postsel));
pe=posr;
for jj = 1:length(postsel)
posr(jj) = local_discreteinvrnd(p(:,postsel(jj)),1);
pe(jj)=p(posr(jj),postsel(jj));
end
lex(ind)=-2*sum(log(pe+eps));
end
le=mean(lex); % expected likelihood
err=2.58*std(lex); % 99%-quantile of standard normal dist.; std(lex) not standard error
% err=2.626*std(lex);% 99%-quantile of Studen-t dist.; std(lex) as empirical standard deviation
ci=[le-err le+err]; % 99% confidence interval
% Bhattacharyya distance
if flags.do_jackknife
x = min([lax;lex]) : max([lax;lex]);
plax = hist(lax,x);
plax = plax/sum(plax);
plex = hist(lex,x);
plex = plax/sum(plex);
Bdist = -log(sum(sqrt(plax .* plex)));
else
Bdist = nan;
end
% Reference Likelihoods ------------------------------------------------
lr = zeros(4,1);
s = 17; % std in degrees
% unimodal
pdfuni = normpdf(rang,0,s);
pdfuni = pdfuni /sum(pdfuni);
iduni = local_discreteinvrnd(pdfuni,nt,runs);
lr(1) = mean(-2*sum(log(pdfuni(iduni))));
% bimodal
pdfbi = 0.5*normpdf(rang,0,s) + 0.5*normpdf(rang,180,s);
pdfbi = pdfbi /sum(pdfbi);
idbi = local_discreteinvrnd(pdfbi,nt,runs);
lr(2) = mean(-2*sum(log(pdfbi(idbi))));
% trimodal
pdftri = 1/3*normpdf(rang,0,s) + 1/3*normpdf(rang,90,s) + 1/3*normpdf(rang,180,s);
pdftri = pdftri /sum(pdftri);
idtri = local_discreteinvrnd(pdftri,nt,runs);
lr(3) = mean(-2*sum(log(pdftri(idtri))));
% unitary
lr(4) = -2*nt*log(1/length(rang));
% normalization of likelihoods, 1 corresponds to unitary
if flags.do_normalize
la = la/lr(4);
le = le/lr(4);
ci = ci/lr(4);
lr = lr/lr(4);
end
end
if flags.do_plot
plot_baumgartner2014_likelistat(la,le,ci,lr)
end
end
function [ X ] = local_discreteinvrnd(p,n,m)
% DISCRETEINVRND implements an inversion method for a discrete distribution
% with probability mass vector p for n trials
% Usage: [ X ] = discreteinvrnd(p,n)
%
% AUTHOR : Robert Baumgartner
if ~exist('m','var')
m=1;
end
p = p/sum(p); % ensure probability mass vector
c = cumsum(p);
t = max(c)*rand(n,m); % rand returns ]0,1]
X = zeros(n,m);
for jj = 1:m
for ii = 1:n
X(ii,jj) = find(c >= t(ii,jj) ,1);
end
end
end