The reliability of the AIC method in Cosmological Model Selection

0.35 MB·English

by  Ming Yang Jeremy Tan

Mon.Not.R.Astron.Soc.000,000–000 (0000) Printed26January2013 (MNLATEXstylefilev2.2)
The reliability of the AIC method in Cosmological Model
Selection
2
M.Y.J. Tan,1Rahul Biswas,1,2
1
0 1Department of Physics, University of Illinois at Urbana-Champaign, Urbana, IL 61820, USA
2 2High Energy Physics Division, Argonne National Laboratory, Argonne, IL 60439-4815, USA

n
a
J Released2002XxxxxXX

3
2
ABSTRACT

The Akaike Information Criterion (AIC) has been used as a statistical criterion to
]
O compare the appropriateness of different dark energy candidate models underlying
a particular dataset. Under suitable conditions, AIC is an indirect estimate of the
C
Kullback-LeiblerdivergenceD(T||A)ofacandidatemodelAwithrespecttothetruth
h. T.Thus,adarkenergymodelwithasmallerAIC isrankedasabettermodel,sinceit
p hasasmallerKullback-LeiblerdiscrepancywithT.Inthispaper,weexploretheimpact
- of statistical errors in estimating AIC during model comparison. Using a parametric
o
bootstrap technique, we study the distribution of AIC differences between a set of
r
t candidate models due to different realizations of noise in the data and show that the
s
shape and spread of this distribution can be quite varied. We also study the rate of
a
[ success of the AIC procedure for different values of a threshold parameter popularly
used in the literature. For plausible choices of true dark energy models, our studies
2
suggest that investigating such distributions of AIC differences in addition to the
v
threshold is useful in correctly interpreting comparisons of dark energy models using
5
the AIC technique.
4
1 INTRODUCTION mentary data sets and a definite framework for model
7
5 selection. Several large new surveys such as DES1, BIGBOSS
5. Suppose we wish to choose the ‘best’ model from a 2 LSST 3 EUCLID4 have been planned to study this late
time acceleration by collecting more data (Abbott et al.
0 set of theoretical models (or theories) of a natural phe-
2005; LSST ScienceCollaborations & LSST Project. 2009;
1 nomenon with the aid of the relevant empirical data. How
Schlegel et al. 2009). Of course, even with accumulation
1 can we objectively accomplish this goal? Many statisti-
: cal techniques such as hypothesis testing (Fisher 1925), of more quality data, the importance of analysing the
v model selection process will not diminish, because reliable
Bayesian evidence (Jeffreys 1961), Akaike information cri-
i
X terion(AIC)(Akaike1974),BayesianInformationCriterion discriminating methods can always allow us to exploit the
available data maximally.
r (BIC) (Schwarz 1978), minimum description length (MDL)
a (Barron, Rissanen & Yu1998),etc.havebeendevelopedto
Hence, statistical techniques addressing model selec-
address thequestion of model selection.
tion have been applied to this context (Liddle 2004;
Liddleet al. 2006; Szydl owski & God lowski 2006;
The late time acceleration of expansion of the uni-
Szydl owski, Kurek & Krawiec 2006; Biesiada 2007;
verse has been firmly established (Riess et al. 1998;
Liddle, Mukherjee & Parkinson 2006; Daviset al. 2007;
Garnavich et al. 1998; Perlmutter et al. 1999; Knop et al.
Liddle 2007, 2009; Trotta 2008; Sollerman et al. 2009;
2003; Tonry et al. 2003; Riess et al. 2004; Astier et al.
Biesiada & Piorkowska 2009; Trotta, Kunz& Liddle 2011;
2006; Wood-Vasey et al. 2007; Hicken et al. 2009b;
Biesiada, Malec & Pi´orkowska2011).Manyofthesediscuss
Freedman et al. 2009; Kessler et al. 2009; Schrabback et al.
the use of information criteria like AIC and BIC which are
2010; Blake et al. 2011), but there is no consensus
easy to calculate. A number of works have applied this to
on the physics behind this phenomenon. A number
data.Inparticular,Szydl owski, Kurek & Krawiec(2006)as
of possible explanations such as a small positive cos-
well as Biesiada (2007) used the method of AIC and BIC,
mological constant or vacuum energy, an otherwise
and a compilation of SNIa to compare various late time
unobserved dynamical fields usually called dark energy
acceleration cosmological models.
(Ratra & Peebles 1988; Wetterich 1988; Frieman et al.
1995; Zlatev, Wang & Steinhardt 1999), or a modification
of General Relativity (Dvali,Gabadadze & Porrati 2000;
Deffayet 2001; Carroll et al. 2005) have been proposed
1 http://www.darkenergysurvey.org/
as an explanation. With many models still consistent 2 http://bigboss.lbl.gov/
with current data, it is clear that further progress of 3 http://www.lsst.org/lsst/
the field requires the collection of larger and comple- 4 http://sci.esa.int/euclid
2 M.Y.J. Tan and Rahul Biswas
Davis et al. (2007) and Sollerman et al. (2009) compared another probability distribution. We denote the KL diver-
such models based on the Sloan Digital Sky Survey gence of a candidate model A (which gives the probability
(Kessler et al. 2009) and Equation of State: Supernovae distributionpA)withrespecttothetruthT (withprobabil-
Trace Cosmic Expansion(ESSENCE) Supernova data and itydistributionpT)asD(T||A),wherewehavedefinedT as
high-redshift data (Riess et al. 2004; Wood-Vasey et al. an underlying process consisting of a signal with stochastic
2007) along with a summary of cosmic microwave back- noise; aparticular empirical datum isa single realisation of
ground (CMB) and Baryon Accoustic Oscillations (BAO) this process. If we denote X as the set of all possible out-
data using AIC and BIC. AIC has also been used in com- comes that can be generated by either A or T, and x as an
paring principal component-based models for dark energy element in this set, we can defineD(T||A) as
(Zhao et al.2007).WewishtoexploretheuseofAICinthe
context of selecting the ‘best’ theoretical model describing D(T||A)= dx pT(x)log pT(x).
thelate time acceleration of theuniverse. Zx∈X pA(x)
Inthispaper,wedefineamodelclassasthetotalityofmodel
This paper will focus on the use and validity of apply-
probability distributions with the same parametric form
ing the AIC technique to choose the ‘best’ late time ac-
(but with different parameter values). Within each model
celeration model from SNIa data. For our calculation, we
class, there is a set of parameters (the ‘best’ model) that
use data from the Constitution compilation (Hicken et al.
gives the lowest KL divergence with respect to the truth.
2009b) of the Center for Astrophysics (CFA3) sample
Thus, to choose the ‘best’ model class we must first choose
(Hicken et al.2009a)ESSENCE(Miknaitis et al.2007),Su-
the ‘best’ model (parameter set) from a particular model
pernova Legacy Survey (SNLS) (Astier et al. 2006) and
classastherepresentativeofthemodelclass.Thisisdoneby
‘High-z’ samples (Riess et al. 2007). We will study and re-
themaximumlikelihoodcriterion.Themodelclassselection
finetheAICmethodology in this context.In particular, we
strategy is thus obtained by comparing the KL divergence
wishtostudyitsreliabilitywhensubjectedtoestimatorun-
oftherepresentativemodelsoftheindividualmodelclasses.
certainty. With the AIC technique, the usual approach has
However,thetruthisunknownapriori,soforagivenrepre-
been to use a single number (called AIC) to rank models;
sentative model A, D(T||A) cannot be evaluated directly.
the smaller the AIC the better the model. This concept of
We can solve this problem by computing the AIC value
discreterankingofmodelshasbeenextended(Akaike1983;
(Akaike 1974) which is an asymptotically unbiased estima-
Burnham & Anderson2004)bypayingattentiontothedif-
tor for D(T||A) (up to a fixed offset that is independent of
ferencesbetweenactualAICvalues.However,itisimportant
therepresentativemodels).Sincethefixedoffsetisindepen-
tonotethattheAICvaluesthemselvesareempirically esti-
dent of the models, and hence the choice of model classes,
matedfromthedataandthushavestatisticaluncertainties.
acomparison oftheAICvaluesisausefulsurrogatefor the
Therefore, thereliability oftheestimatesofAICdifferences
strategy of comparing the associated KL divergence of the
isacrucialissue.Whilethereiswork(Shimodaira1997)that
differentcandidatemodelclasses. Wecan computetheAIC
triestoaddressthereliabilityoftheAICtechniquethrough
oftherepresentativemodelsofvariousmodelclasses,which
analytic calculations, we will instead study the reliability
in the asymptotic limit is known (Akaike1974) to be
through numerical simulations.
The paper is organised as follows. In Section 2, we briefly AIC=2k−2log(LML), (1)
review the information theoretic origin of AIC, extending where LML denotes the likelihood evaluated at the set of
this idea further, by reviewing the theoretical implications model parameters that maximize the likelihood L, and k is
oftheAICdifferencesinSection3.InSection4,westudythe thenumberoffreeparametersinthecandidatemodelclass.
estimatoruncertaintyoftheseAICdifferencesinthecontext Whenwemakeafurtherassumptionthatthedistributionof
ofcosmological modelselection. Thisinvolvesacomparison errors follows a Gaussian distribution, this further reduces
of a set of four candidate models by using SNIa, as in the to
work by Szydl owski & God lowski (2006). The reliability of
these AIC estimates are calculated via a bootstrap method AIC=2k+χ2ML, (2)
(Efron & Tibshirani1993),andthiswillbeusedtoevaluate
where χ2 is the usual chi-square evaluated at the maxi-
thevalidityof themodelselection technique.Theresults of ML
mumlikelihoodestimateofthemodelparameters.TheAIC
these simulations and the conclusion will then be given in
valuesaresubsequentlyrankedbythesmallnessoftheirval-
Sections 5 and 6, respectively.
ues; the model class with the smallest value is determined
to bethe‘best’ model class.
It is worth noting that the AIC estimate (Eqn. 2) of the
2 THE AKAIKE INFORMATION CRITERION KL divergence assumes that the number of data points is
(AIC) sufficientlylarge. AICintheform writtenaboveisanunbi-
asedestimatorforlargedatasets.Forsmallerdatasets,the
The Kullback-Leibler (KL) divergence (Cover & Thomas 2k term can be corrected byan additional 2k(2k+1) term to
N−k−1
1991) is a commonly used quantity that measures the dis- approximately correct for the bias due to finiteness of the
crepancy 5 of one probability distribution with respect to dataset, where N is the number of data points in a single
data set. While further studies to obtain a more accurate
expression for this term are possible, for the cases we shall
5 Note that a KL divergence is a non-symmetric measure that consider,thiscorrection isalwaysless than0.06(which will
doesnotobeythetriangleinequality. be seen to be negligible for our purposes) and will only de-
The reliability of the AIC method in Cosmological Model Selection 3
crease in importance when more data is collected. Weshall while any two models with an AIC difference smaller than
therefore ignore this correction altogether in this paper. ∆ are considered of equal rank. Eqn. 3 shows that
threshold
choosing a large enough value of the threshold ∆
Asalreadymentionedexplicitly,theAICprocedureforcom- threshold
implies a high probability that the selected model is truly
paring different model classes actually compares the χ2 for
thebetterone.However, a large valueof ∆ also in-
the best fit representative model from each model class. threshold
creasesthenumberofmodelpairswheretheAICdifferences
Theserepresentativemodelsarederivedfromthemaximum
likelihood estimate of themodel parameters for each of the are in the range −∆threshold < ∆A,B < ∆threshold. Since
this procedure cannot discriminate between such models,
modelclass.Thus,intherestofthispaper,whenwereferto
weshallcallsuchamodelselection resultindeterminate.In
comparing model classes, the calculation only involves the
our convention, we also define the converse of the indeter-
χ2 ofthe(best)representativemodeltakenfrom itsrespec-
tive model class. We sometimes refer to this as comparing minate case as the determinate case (|∆A,B|>∆threshold).
Of course, for a pre-determined choice of ∆ (cor-
models. threshold
responding to a predetermined confidence level), a better
dataset gives a smaller fraction of model pairs which have
indeterminateresults.Asaruleofthumb,auniversalvalue
3 MODEL COMPARISON AND AIC of the threshold ∆threshold = 5, without any regard to the
DIFFERENCES properties of the models under comparison, has been men-
tionedbyLiddle(2007)astheminimumAICdifferencebe-
Since AIC is essentially the measure of discrepancy of a tween two models needed to make a ‘strong’ assertion that
model from the truth, it is intuitively obvious that the one model is better than the other. Such a definition has
smallertheAICdifferencebetweentwomodels,theharderit been used extensively in theliterature.
becomestojudgewhichmodelisbetter;eveniftheAICesti-
mateofthisdifferenceintheKLdivergencecanbeobtained
withoutanestimatorerror,thesmalldifferencewouldmake
it difficult to tell the two probabilistic models apart for a 4 IMPACT OF AIC UNCERTAINTIES IN
small number of observations. Hence, there is a need to as- FINITE DATA SETS
sociateaconfidencelevelfordistinguishingbetweenamodel
A and another model B using the AIC difference between In the preceding section, we have discussed the AIC
them. differences ∆A,B without any regard for the fact that AIC
is a statistical estimate. The associated uncertainty in the
LetP(M)betheprobabilitythatmodelM istrue.Wewish
AICestimate may not be negligible and may be dependent
to relate P(A)/P(B) tothe AICdifference
on the realisation of noise in a particular data set. Thus,
∆A,B ≡AIC(A)−AIC(B)≈2[D(T||A)−D(T||B)], theremustbeastatistical uncertaintyin thevalueof ∆A,B
6 even when estimated from a data set of similar quality
where ≈ implies the asymptotic relation. This asymptotic 7. The ensemble of such observations defines an empirical
relation may not necessarily be realised as the number of
probability distribution of ∆A,B, and the particular value
samples is small. Using a certain set of extra assumptions,
of ∆A,B obtained from the current SN data sets is actu-
Akaike(1983) showed that for a pair of models A and B,
allyasamplevaluedrawnfromthisprobabilitydistribution.
P(A)/P(B)≈exp[−∆A,B/2]. (3)
Ideally, we should be able to study the probability
Since the AIC differences are estimates of the differences
density distribution of ∆A,B under repeated observations
of KL divergences, one can also use the idea of distin-
of results with sample size N: P(∆A,B|EN), where EN
guishing probability distributions to justify Eqn. 3 when
denotes a collection of observation data which individually
the number of samples are sufficiently large (Appendix A).
consists of N observation points and are drawn from
Assuming the truth to be in the set of candidate models,
the underlying truth process. Because producing a large
Burnham & Anderson(2004)extendedEqn.3toobtainthe
subset of EN is impossible, in this paper we instead use a
probability wi of model i by appropriately normalising the bootstrap approach (Efron & Tibshirani 1993) to generate
equation:
‘mock’ empirical data sets and estimate the probability
wi = enxepx[−p[∆−i∆,b/2/]2], (4) distribution of ∆A,B.
j j,b
Perhaps the most frequently used method to produce
wherePthe candidate models are numbered as i = 1,···,n ‘mock’empirical dataisthebootstrap methodproposed by
and b denotes the best model which has the smallest AIC
Efron & Tibshirani(1993).Whenweapplythisapproachto
value among all the candidate models. Either way, Eqn. 3
quantifiestheintuitiveideathatitiseasiertoselectamodel
overanother if theAICdifference is large. 6 ThisisduetoastatisticaluncertaintyintheAICvaluescoming
fromthevariationofχ2 .However,theuncertaintyintheAIC
One can use this idea to modify the AIC methodology de- ML
valuesofthemodelscanbecorrelated,andturnsouttobelarger
scribed above and suppress the probability of obtaining in-
thantheuncertaintyintheAICdifferences.
correct resultsbyintroducingathreshold ∆threshold.Then, 7 In this paper, two Supernovae data sets are said to have the
ratherthanrankingallmodelsaccordingtothesmallnessof same quality when they have the same number of data points,
theirAICvalues,oneadoptstheprocedurewhereamodelA the same set of redshift z values and the same set of error
isrankedtobebetterthanmodelBif∆A,B <−∆threshold, bars(standarddeviation) thatcorrespondstothesetofz values.
4 M.Y.J. Tan and Rahul Biswas
Table1.Differentmodelclassesconsideredinthispaper:WeshowtheevolutionoftheHubblefunctionH(z)withredshiftz,theHubble
constantH0 andotherfreeparametersinthemodels.kisthenumberoffreeparameters.TherespectiveAICvalueswereobtainedfrom
theConstitution compilation.
Model H(z)/H0 FreeParameters k AIC
ΛCDM pΩm(1+z)3+(1−Ωm) Ωm 1 401.35
wCDM qΩm(1+z)3+(1 Ωm)(1+z)3(1+w) Ωm,w 2 403.05
−
CPL qΩm(1+z)3+(1−Ωm)(1+z)3(w0+wa+1)exp[−13+wzaz] Ωm,w0,wa 3 404.66
DGP (pΩm(1+z)3+Ωrc+√Ωrc) Ωm 1 401.13
Ωrc=(1 Ωm)2/4
−
regression models, we need a probability model that speci- distributionofmean0andstandarddeviationσ2:N(0,σ2),
i i
fiesthedistributionofresiduals(eg.Gaussiandistribution). and f(z) is the maximum likelihood estimate of model C
In our case in particular, the distance modulus µ is related (ΛCDM, DGP, etc.) (AppendixC).
to the red shift z, so a regression model has the following
structure µ = f(z)+ǫ, where ǫ is the residual (the error As mentioned above, we wish to simulate the distri-
term).Sincewedonot knowthetruerelation f, wecannot bution of {∆dA,B|d ∈ CN} as a proxy for P(∆A,B|EN). To
obtain the purely empirical distribution of the residuals proceed,wechoose asubset of modelsthat havebeen often
8. Thus, for regression models, the standard bootstrap studied in the literature (Szydl owski & God lowski 2006)
method always involves using a model-dependent probabil- andarelistedinTable.1.Wealsoassumedthattheuniverse
ity distribution of residuals. Daviset al. (2007) extended isspatiallyflatandsetthecurvaturetermΩ intheHubble
k
thismethodtocomparebetweentworegression modelsand function H(z) to zero. Three of these models ΛCDM,
find the standard deviation of their BIC differences; we wCDM, CPL (Chevallier & Polarski 2001; Linder2003) are
will further extend this idea and show that studying the dark energy models with different parametrisations of the
structureofthedistributionofAICdifferencesisimportant. equations of state w(z) = p(z)/ρ(z), where p(z) and ρ(z)
are the pressure and density of dark energy, respectively.
When we wish to check the reliability of a particular Thesemodelsarenested:settingwa =0intheCPLmodel,
model,wecouldusetheprobabilitydistributionofresiduals we obtain the wCDM model; setting w =−1 in the latter
0
based onthemodelitself. Sincewecannothaveanymodel- gives the ΛCDM model. We also use the flat DGP model
free bootstrap data, we must choose a particular model to which is a modified gravity model and cannot be nested in
produce the residual distribution. Therefore, to estimate theprevious classes of models.
∆A,B, we need some model C as a reference probability
model that is used to generate the bootstrap data. Let us To make contact with observational data, we choose
denotetheestimateof∆A,B basedonadatasetdas∆dA,B. candidate models with parameter sets that are ‘best’ for
We wish to produce {∆dA,B|d ∈ CN}, where CN denotes a the Constitution compilation of SNIa data (f(z) being
collection ofparametricbootstrap datagenerated bymodel the maximum likelihood estimation derived model), where
C, which individually consist of N observation points with we use the distance moduli and the error bars in data
the same data quality as our empirical data. Although the (Appendix C). Since we need to find estimates of cos-
needed details are in Appendix B, our approach may be mological parameters for different models by maximising
outlined as follows. likelihoods, we use the results from the more appropriate
MLCS light curve fitter (for RV = 1.7) in (Hicken et al.
Suppose we have a single data set consisting of N ob- 2009b). The details of finding the maximum likelihood and
served results {(z1,µ1,σ1),(z2,µ2,σ2),···,(zN,µN,σN)}, the corresponding AIC value is given in Appendix C; 371
where the zi values denote theith observed redshift, µi the SNIa events were used. The Hubble function, as a function
ith generated distance modulus, and σi the observed error of cosmological parameters in each of these models, along
barsofµi.WecancreateabootstrapsampleCN consisting with the free parameters and AIC values (calculated from
of N observation points based on model C. CN relates the theConstitution compilation) are shown in Table. 1.
same set of coordinates by the relation: µi = f(zi) +ǫi,
i=1,···,N,whereǫi isastochastictermobeyinganormal For the nested models, the parameter values that fit
the data best turn out to be close to the ΛCDM model.
Consequently, one does not gain much in terms of a lower
8 Note that these errorswhich include lightcurve fitting errors, χ2 , while the extra free parameters are penalised to
ML
intrinsic dispersion and peculiar velocity corrections are also re- give higher values of AIC for wCDM and CPL. The DGP
quiredforcalculatingquantities likeχ2 formostmodelselection
model gives the best AIC value, which is only slightly
schemes.Thus,anunderlyingassumptionoftheapplicationofthe
better than the ΛCDM model. It is known that the
AICtechnique asinprevious worksisthat theseerrorestimates
simultaneous useof Cosmic Microwave Background (CMB)
are correct. Since our focus is on the statistical uncertainties in
data from WMAP (Daviset al. 2007; Sollerman et al.
AICafterfollowingotherunderlyingassumptionsusedinthelit-
erature,wealsoassumethattheseerrorestimates arecorrect. 2009) and Large Scale Structure (LSS) data disfavours
The reliability of the AIC method in Cosmological Model Selection 5
the DGP model compared to ΛCDM, since the parameter
X Y DGP ΛCDM wCDM CPL
subspaces that provide the best fits for CMB, LSS and \
SNIa data do not overlap as much as in the ΛCDM model.
We checked that our results are consistent with this, but DGP - 58% 84% 89%
willignoretheCMBandLSSdatatofocusonmethodology.
ΛCDM 59% - 85% 90%
In order to calculate the distribution of ∆d , we adopt
A,B
wCDM 16% 21% - 90%
the best fit model in a model class C as a reference model
to produce 5000 mock data sets of C which is expected
371 CPL 14.0% 16% 17% -
to be similar in quality to the Constitution compilation
of Supernova data. All the model classes in Table. 1 are
successively chosen as the reference model C. Following
Table 2. Percentage of cases where the AIC method with a
our definition of ‘similar quality’, each simulated data set
threshold ∆threshold = 0 selects the correct model over other
has the same set of redshift values and error bars as the
candidate modelsconsidered. Wedefined thedifference andper-
Constitution compilation, while their apparent magnitudes centageusingthefollowingconvention.Ifwedefinethebootstrap
are those expected from a noisy realisation of the reference data as being produced by the reference model X (rows) and
model. Thebasic steps are: we are comparing it against model Y (columns), the difference
is defined as the AIC value of X minus the AIC value of Y and
(i) We produce a mock data set consisting of 5000 real- is denoted by the symbol ∆d . The table counts the percent-
X,Y
aisbaotivoen.s of d ∈ C371 for a reference model C as outlined age of negative ∆dX,Y,d ∈ XN. Note that a value greater than
50%indicates thatthecorrectmodelischosenamajorityofthe
time. We use this cut-off of 50% or more as a simple definition
(ii) Candidate models A and B are fitted to d ∈ C
371 forstatisticalself-consistency.
by maximising the likelihood and the AIC values of these
models A and B are computed through Eqn.2.
(iii) Thus, for each element d ∈ C we can make the
371 disqualified underthecurrent numberof data points in the
AIC difference ∆d . We study {∆d | d ∈ C } for the
A,B A,B 371 sample. In another example, we compared the DGP and
reference model C byplotting a histogram.
ΛCDM models in the table and noticed ∆d < 0,
ΛCDM,DGP
Weshould notethattheprobability distributionof ∆d is 59% of the time when ΛCDM model is the reference (i.e.,
A,B
due to errors introduced by the stochastic noise term, but P(DGP/ΛCDM)=0.594), and42%ofthetimewhenDGP
ignores the effect of the uncertainties in the cosmological model is the reference (i.e., 1−P(ΛCDM/DGP) = 0.422).
parameters of thereference model C. ThismeansthatifthereferencemodelwaseitherΛCDMor
DGP and we apply AIC to compare between them using a
zero ∆ , AIC will only slightly favour the reference.
threshold
The AIC technique is only statistically self-consistent (for
5 BOOTSTRAP STUDY RESULTS a zero ∆threshold) when applied to compare DGP and
ΛCDM while we cannot use AIC to self-consistently study
We first consider the issue of statistical self-consistency in the other models under the current level of observation
the following sense: When the reference model is C, does quality. However, a test that gives the right answer 3 out
the bootstrap AIC method outlined above choose model C of 5 times is unreliable, since we can only do a single
as a better model than the rest? The distribution of the empirical test from our actual data. ΛCDM and DGP
values of ∆dC,A for d ∈ CN can tell us about the statistical cannot be distinguished significantly using either reference
self-consistency.Westartwiththecasewhen∆ =0. models. Thus, looking at both examples, we must conclude
threshold
If P(A/C) ≡ P({d ∈ CN|∆dC,A < 0}) is larger than some that there are insufficient data points to tell reliably the
predefined proportion, which must not be less than 1/2, models apart using AIC when ∆threshold = 0. Another
we may say model C is better than A when C is the trend that results from the insufficiency of data points is
reference. The P(A/C) results for this case is summarised that the AIC procedure tends to favour models with a
inTable.2.TheP-values9forthetableareextremelysmall. smaller number of free parameters. The trend persists even
when we later increase ∆ . This seems to indicate
threshold
From Table. 2, we obtained ∆d < 0, 89% that the addition of extra free parameters does not signifi-
DGP,CPL
of the time, when DGP model is the reference (i.e., cantlyimprovetheχ2fitforthenumberofdatapointsused.
P(CPL/DGP) = 0.892), and 83% of the time, when CPL
model is the reference (i.e., 1−P(DGP/CPL)) = 0.836). We next study the behaviour when the threshold pa-
This means that, when we use a threshold ∆threshold = 0, rameter is increased, ∆threshold = 2 and 5, corresponding
the DGP model is significantly favoured over the CPL to choices made in the literature to moderate and strong
model even if the CPL model is the reference. This means evidence. Unlike the previous ∆threshold = 0 case, we
that the AIC method is statistically inconsistent when it have to consider the effect of indeterminate cases, which
is applied to the CPL model, and would be automatically wedefinedinsection3asthecasewhen|∆A,B|<∆threshold.
In order to study the reliability of the AIC technique
9 ThehighestP-valuefoundwas3×10−8for∆DGP,ΛCDM.Most atdifferentvalueofthethresholdparameter∆threshold >0,
oftheP-valuesfoundwereatthelevelofmachineprecision. we analyse the probability of the selected model being
6 M.Y.J. Tan and Rahul Biswas
incorrect for different values of ∆ . To do so, we of free parameters than the reference model, fdet tends
threshold false
definethefollowing: to gently increase before steeply decreasing at a certain
value of ∆ . From Fig. 1, we can see that this
threshold
f = Numberof samples with |∆A,B|<∆threshol(d5) sharp decline happens at roughly twice the difference in
ind
Numberof samples the number of free parameters between the candidate and
fall = Numberof samples with ∆A,B >∆threshold(6) reference models. We also study these quantities for the
false Numberof samples case where the reference model has less parameters than
fdet = Numberof samples with ∆A,B >∆threshold(7) a candidate model in Fig. 2. In this case, ffdaetlse actually
false Numberof samples(1−find) increases rapidly at the point where ∆threshold is equal
f is the fraction of cases where the AIC procedure has to the difference in the number of free parameters in the
ind
an indeterminate result for a given value of the threshold two models. This implies that if the model underlying the
∆ , so a high value of f reflects the inadequacy empirical data was similar to the reference models studied,
threshold ind
of the data to discriminate between the pair of models in it is improbable that the data set would provide an AIC
question with a certain level of confidence for a relevant difference for the considered models greater than a large
∆ . Using A as the reference model, fall is the threshold ∆ (eg. 5) as shown by f . However,
threshold false threshold ind
fraction among all cases, where the AIC procedure results if this dataset did yield a AIC difference larger than a
inanincorrectmodelselection.fdet isthefraction among predetermined∆ ,itdoesnotnecessarily meanthat
false threshold
determinate cases (|∆A,B| > ∆threshold) where the AIC the AIC selected model has a high probability of being the
procedure results in an incorrect selection, and reflects the true underlying model. This is because fdet for a given
false
ratio of correct to incorrect model selections. Our results ∆ depends on the model pairs being considered,
threshold
are summarized in Tables. 3 and 4. In each table, the suggesting that even having AIC differences larger than a
rows correspond to different reference models X, while the specified threshold ∆threshold does not guarantee reliability
columns list candidate models Y that were compared with of theAIC selection process.
thereference model.
For each pair of reference model X and candidate In order to intuitively understand what leads to these
model Y, we show the fractions (find,ffaalllse,ffdaetlse). First examples, it is instructive to consider the shapes of the
we reconsider the ∆threshold = 0 case in Table. 2. By distribution of AICdifferences when comparing theΛCDM
definition, the proportion of ∆threshold =0 cases where the and CPL models or the ΛCDM and DGP models. We note
AICmethodisnotstatisticallyself-consistent(1−P(A/B)) that a comparison of nested models will always involve
is equivalent to both ffaalllse and ffdaetlse. We note that the strongly asymmetric, exponential-like distribution of AIC
values of ffalse are large, indicating a high failure rate and differences, while the comparison of non-nested models
an unsatisfactory procedure. could result in almost symmetric distributions. This is
because the χ2 for the general case cannot be smaller
ML
We expect these failures to be suppressed when we thanthemorespecificcase.Forexample,ΛCDMisaspecial
choose higher values of the threshold ∆threshold. When we case of the more general wCDM and CPL models. Hence,
increase thethreshold∆threshold to2and then5,wenotice when we fit these models against any data or simulated
theexpectedsuppressionofffdaetlse.However,ffdaetlse doesnot data(regardlessofthereferencemodelusedtogenerateit),
decrease as dramatically as ffaalllse. thewCDM and CPL χ2ML valuescannot belarger than the
ΛCDM χ2 value and are usually smaller. This results in
ML
We study the behaviour of ffdaetlse, ffaalllse and find in a sharp edge along a maximum value of ∆ΛCDM,CPL and
Fig. 1 for different values of ∆threshold as well as different an exponential-like distribution of AIC differences between
choices of candidate and reference models; using bootstrap thesemodels.
simulations. The fdet values become dominated by noise
false
as find increases, since the calculation is made from an We show examples of such a Gaussian-like distribu-
ever decreasing number of determinate bootstrap cases. tion in Fig. 5 and exponential-like distribution in Fig. 4.
Hence, the values of ffdaetlse at large values of find should As an aside, we note that the different reference models
be ignored. Nevertheless, we can still study the regime for considered do not make much difference in the shape of
smaller find. It should also be noted that for the case of these distributions in the figures. Since AIC = χ2ML +2k,
∆threshold =5, the proportion of indeterminate results was the edge of these exponential-like histograms is shifted by
high (Table. 4). For example, when CPL is the reference twice the difference between their number of parameters.
model, the proportion of ∆dCPL,ΛCDM and ∆dCPL,DGP Hence for the comparison between the CPL model and
between ±5are both approximately 98%. ΛCDM model in Fig. 1 and Fig. 2, where the difference
between the χ2 is almost zero, the distribution of AIC
ML
We also note, as expected, that find increases asymp- differences ∆ΛCDM,CPL has a sharp edge at approximately
totically to one as we increase ∆threshold. Increasing 4. This implies that any AIC difference greater than a
∆threshold monotonically suppresses ffaalllse, but not nec- ∆threshold value of approximately 4 will exclusively select
essarily fdet . fdet decreases with increasing ∆ theCPLmodel,irrespectiveofwhetherthereferencemodel
false false threshold
for the cases where the DGP model was wrongly picked used was a CPL model or a ΛCDM model. At lower values
over the reference ΛCDM model. However, for the other of ∆threshold, the AIC selection procedure tends to select
cases, when the candidate model has a smaller number the model with the lower numbers of parameters since the
χ2 values are approximately the same. On the other
ML
The reliability of the AIC method in Cosmological Model Selection 7
X Y DGP ΛCDM wCDM CPL
\
DGP - (97,1,26)% (86,5,32)% (25,4,5)%
ΛCDM (97,1,24)% - (94,4,66)% (24,4,5)%
wCDM (87,9,71)% (90,1,11)% - (96,3,79)%
CPL (31,64,92)% (32,61,90)% (94,0,0)% -
Table 3. Failure Rate for ∆threshold = 2 : The values in the parentheses show find, the percentage of total cases where the AIC
procedurehasanindeterminateresult;fall ,thefractionoftotalcaseswheretheAICprocedureresultsinanincorrectmodelselection;
false
andfdet ,thepercentageofdeterminatecaseswheretheAICwronglyselectsacandidate modelYoverthereferencemodelX.
false
X Y DGP ΛCDM wCDM CPL
\
DGP - (100,0,“-”)% (99,1,100)% (99,1,100)%
ΛCDM (100,0,“-”)% - (99,1,100)% (99,1,100)%
wCDM (100,0,0)% (99,1,0)% - (99,1,0)%
CPL (99,0,0)% (98,0,0)% (99,0,0)% -
Table 4. Failure Rate for ∆threshold = 5 : The values in the parentheses show find, the percentage of total cases where the AIC
procedurehasanindeterminateresult;fall ,thefractionoftotalcaseswheretheAICprocedureresultsinanincorrectmodelselection;
false
andfdet ,thepercentageofdeterminatecaseswheretheAICwronglyselectsacandidate modelYoverthereferencemodelX.
false
hand,thereisnosimilar constraintrelatingtheχ2 values 68.3% of the differences. We chose 68.3%, because 0.683 is
ML
oftheΛCDMandDGP(atthebestfitvalue);consequently approximatelytheprobability offindinganoutcomewithin
the histogram turns out to be Gaussian-like. They have astandarddeviationofthemeanofaGaussiandistribution.
the same number of free parameters so this distribution is
roughly centred about zero (the difference in χ2 between As an illustration, we look at statistical uncertainty
ML
these models). Obviously, if the quality of data was much of the distribution of ∆d for the two separate
ΛCDM,DGP
better in terms of the error bars on each observation, or cases where the DGP and ΛCDM models are the refer-
having a larger number of observations, the differences in ences. The distribution of the differences in both cases are
theχ2 termswouldbemuchlarger forthesamechoiceof Gaussian-like. We notice that the standard deviation of
ML
referencemodels used.Insuchcases, themodelcomparison ∆d is 0.89 with DGP model as the reference and
ΛCDM,DGP
would be purely data-driven, deriving its discriminatory 0.83withΛCDMmodelasthereference.TheAICdifference
power from the ‘fit’ term of AIC. This agrees with our between the DGP model and ΛCDM model ∆ΛCDM,DGP
intuitive idea that a better dataset should be able to to in our original AIC analysis is 0.22 and smaller than the
resolve models better. error bars. This result means that any subsequent analysis
based on the value of ∆ΛCDM,DGP observed is unreliable.
The shapes of the distribution of AIC differences are also Theerrorbarscouldalsobesignificantevenif∆ΛCDM,DGP
was larger than the error bars as we would have to modify
important when we study the spread of the statistical
any subsequent AIC difference analysis to include this un-
uncertainty of the distribution of the AIC differences since
certainty. When we look at ∆d (exponential-like
thestatistical spreadofthedistributioncannotbespecified ΛCDM,CPL
distribution), the error bar region ranges from −4.00 to
unless we know the shapes beforehand. Due to the struc-
−2.23 for the case when the reference is the ΛCDM model
tural differences between the two kinds of distributions, we
and ranges from −4.00 to −1.42 when the reference is the
must definethe‘error’ bars according totheshape in order
to make any useful comments about the uncertainty of the CPL model. The ∆ΛCDM,CPL value calculated from the
empirical data was found to be −3.31. However, the error
estimate. For the case of the Gaussian-like distribution,
bar region range of approximately 2 would make the value
we define the error bars to be the standard deviation of
of −3.31 less certain. Instead of quantifying the odds ratio
the statistical distribution of differences. For the case of
given by Eqn. 3 as having a value of 0.19, we now make a
the exponential-like distribution, we define an error bar
statement about its uncertainty by saying that the odds
region as the range that begins at the sharp edge of the
ratio can be a value between 0.14 and 0.48. These are
distributionandstopsatthepointwheretherangecontains
8 M.Y.J. Tan and Rahul Biswas
1.0 1.0
f
ind
0.8 Ref=LCDM 0.8 Ref=LCDM
ns0.6 Cand=DGP ns0.6 Cand=CPL ffadleste
o o
i i
t t
rac0.4 ffadleste rac0.4 find
f f
0.2 0.2
f all f all
false false
0.0 0.0
0 1 2 3 4 5 0 1 2 3 4 5
D D
threshold threshold
1.0
f det Figure2. Probabilityofthecandidatemodel(Cand)CPLbeing
false f
ind selectedoverthereferencemodel(Ref)ΛCDMfordifferentvalues
0.8 Ref=wCDM of the threshold. In this case, fdet actually increases sharply
false
at about twice the difference in free parameters between these
models,showingthatinthiscaseAICtends toincorrectlyselect
ns0.6 Cand=LCDM themodelwithahighernumberofparameters.
io f all
t false
c
ra0.4 just two examples in which one can carry out an analysis
f
to determine the reliability of the model likelihood ratio
0.2 P(A)/P(B)that is calculated in Eqn.3.
We note that the statistical uncertainty of the AIC differ-
0.0
ences obtained above is smaller than the ∆ value
0 1 2 3 4 5 threshold
of 5 already mentioned above (Liddle 2007; Davis et al.
D
threshold 2007; Sollerman et al. 2009). However, it is still significant
1.0 enough and needsto be considered in our analysis.
f det
0.8 false Ref=CPL
f all
s0.6 false Cand=LCDM
n
6 CONCLUSION
o
i
t
c
a0.4 AIC has been widely used as a technique for model selec-
r
f f tion. Most commonly, this has been applied by computing
ind
0.2 the AIC values for each candidate model through Eqn. 2
and selecting the model with the smallest AIC value as
the best model. The issue of considering the magnitudes
0.0
of AIC differences between the models to indicate the
0 1 2 3 4 5
relative plausibility or confidence in the models has also
D been addressed by Akaikeand elaborated by Burnham and
threshold
Anderson through the use of Akaikeweights Eqn. 4. In the
Figure 1. Probabilities of different candidate models (Cand) field of cosmology, AIC has been used in selecting models
beingselectedoverdifferentreferencemodels(Ref)byusingAIC underlying the late time acceleration of the universe using
fordifferentvaluesof∆thresholdforthecaseof(a)(top)thecan- data. There have also been suggestions of a rule of thumb
didatemodelDGPbeingpickedoverthereferencemodelΛCDM thatstatesanAICdifferenceof5ormorewouldgivestrong
(b) (middle) the candidate model ΛCDM being picked over the
evidence for the model with the smallest AIC value. This
reference model wCDM, and (c) (bottom) the candidate model
approximately corresponds to a ratio of model likelihoods
ΛCDM being picked over the reference model CPL. The solid
of 12 or more.
blue curve shows the number of incorrect results as a fraction
fdet of the cases where the procedure returns a determinate
false
result.Theblackdashedcurveshowsthenumberofincorrectre- In this paper, we propose a method for calculating
sultsasafractionofthetotalnumberofsimulations.Thefraction the ratio of the model likelihood between models A and
ffdaetlseisextremelynoisyandshouldbeignoredwhenthefraction B based on their AIC differences ∆A,B (Appendix A),
of determinate cases find (red, dotted) is large. The plots show using the idea of probabilistic model distinguishability by
that increasing ∆threshold always decreases ffaalllse, ffdaetlse does Balasubramanian (1997). This is an alternative method
not necessarily decrease. For a comparison of the models used, of arriving at the odds ratio of Eqn. 3 assuming that the
this shows that AIC tends to incorrectly select models with a
lowernumberofparameters.Inthecomparisonofnestedmodels,
The reliability of the AIC method in Cosmological Model Selection 9
0.7 0.7
LCDM LCDM
0.35 0.35
0 0
-5 -3 -1 1 3 5 -5 -3 -1 1 3 5
0.7 0.7
wCDM wCDM
0.35 0.35
0 0
-5 -3 -1 1 3 5 -5 -3 -1 1 3 5
0.7 0.7
DGP DGP
0.35 0.35
0 0
-5 -3 -1 1 3 5 -5 -3 -1 1 3 5
0.7 0.7
CPL CPL
0.35 0.35
0 0
-5 -3 -1 1 3 5 -5 -3 -1 1 3 5
Figure 4. Probability distributions of the AIC differences be-
Figure 3. Probability distributions of the AIC differences be- tween ΛCDM and CPL (∆d ) for different reference
ΛCDM,CPL
tween ΛCDM and DGP (∆d ) for different reference model C. C is used to generate the respective bootstrap sam-
ΛCDM,DGP
modelC.C isusedtogeneratetherespectivebootstrapsamples ples and is written in the right upper corner of the figures. The
and is written in the right upper corner of the figures. The hor- horizontalaxisindicatesthe∆d valuewhilethevertical
ΛCDM,CPL
izontal axis indicates the ∆d value while the vertical axis indicates their relative frequency. If the process underlying
ΛCDM,DGP
axis indicates their relative frequency. If the process underlying ourobservationswasreallythebestfitmodelofclassC,thenthe
ourobservationswasreallythebestfitmodelofclassC,thenthe valuesoftheAICdifferences underdifferentrealisationsofnoise
vwaoluuledshoafvtehtehAeIhCistdoigffrearmendceisstruinbduetriodniff{∆erdΛenCtDrMea,lDisGaPtio|nds∈ofCn3o7i1s}e wshoouwldn.hTavheetehxephoinsetnotgiraalm-likdeisdtrisibtruibtiuotnio{n∆sdΛoCbsDeMrv,eCdPaLre|ddu∈eCto37th1}e
shown. Vertical lines show the respective AIC differences that factthatΛCDMandCPLarenestedmodels.Verticallinesshow
werederivedfromtheobserveddata. the respective AIC differences that were derived from the ob-
serveddata.
10 M.Y.J. Tan and Rahul Biswas
AIC difference between the candidate models is a perfect exponential-like distribution with a sharp edge and one
unbiased estimator of the difference between the models’ sided tail. Thus, in order touse AICreliably, one must pay
KLdivergenceswithrespecttothetruth.Arelatedanalysis properattention to thestatistical variation of ∆A,B.
using the AIC differences in accordance with the Akaike
weights also derives the same equation. The analysis of In this paper, we made a number of assumptions to
∆A,B was extended further by investigating the statistical study the AIC technique. All calculations in this paper
uncertainty of this estimate. Our focus was not necessarily were only for an assumed reference. Since the empirical
on the ‘best’ cosmological theory. Thus, we did not use data does not give us EN, there is no way to know the
the most exhaustive data sets, nor explore in detail the actual distribution of ∆A,B. However, we should note that
systematics associated with the surveysconsidered. AIC is a model comparison technique that assumes that
one of the model classes contains the reference C. By
To this end, we studied the distribution of the differ- restricting C to the set of candidate models, we can at
ences of AIC estimates given a certain quality of data least look for statistical self-consistency in that assump-
P(∆A,B|EN). Since we do not know the exact process un- tion. It should be emphasised again that the reference
derlyingtheempirical data(EN in Section 4),weapproach models used were the best fit models and did not take
this problem by studying surrogate processes, where the into account the statistical uncertainty of the individual
reference or generating model is assumed to be one of four model parameters. That can be taken into account by
(best fit) candidate models for the late time acceleration of sampling the distribution of parameters. As mentioned
the universe; following the approach used by Davis et al. before, the whole point of this simulation is to highlight
(2007). These models were listed in Table. 1 with the best the statistical distribution of the AIC differences under
fit free parameters equal to the maximum likelihood values different reference models and look for statistical inconsis-
of each of these models. This was obtained by fitting 371 tencies under each of the assumptions. Another issue that
SNIaextracted from theConstitution compilation. should be noted is that the exact variance in the data is
unknown and that the error bars in the data may not be
Our simulations have demonstrated that, given the reflective of the true error bounds. While we use these in
data used, there was insufficient data to reliably use our simulations, we note that the correctness of these error
AIC to tell all the models apart; in agreement with the estimateswasanassumptionofpreviousAICcomputations.
general consensus (Daviset al. 2007). For the case of
∆ =0, the failure rate of the technique was shown In summary, the reliability of the AIC estimator is an
threshold
to be particularly unsatisfactory. We also studied the important issue that should be taken into consideration
reliabilityoftheAICtechniquewhen∆ isincreased when using the AIC technique to select models. It should
threshold
to 2 and 5. Increasing ∆ results in increasing the alsobenotedthatsuchconsiderationsarenotjustrestricted
threshold
number of cases where we cannot make a conclusion based to AIC but any technique that relies on the maximum
on the AIC procedure. This was demonstrated by f , likelihood estimators. This should be borne in mind when
ind
whichcalculates thefraction ofcaseswhenthedifferencein applyingthe techniquestoany statistical analysis.
AIC values between two models is less than ∆ . We
threshold
also studied fdet , the proportion of cases where the AIC
false
procedure using a threshold ∆ gives an incorrect
threshold
resultasafractionofcaseswherewecanmakeaconclusion 7 ACKNOWLEDGEMENTS
(i.e. |∆A,B| > ∆threshold). We showed that ffdaetlse does not
necessarily decrease in the same universal way with an We would like to thank Yoshi Oono for checking the
increasing ∆threshold. Therefore, even when AIC chooses a manuscriptandprovidingmanyinsightsaboutthenatureof
model class (with a high level of ∆threshold), the result is AIC and the bootstrap technique. M.Y.J. Tan was funded
unreliableforatleast somemodelswithin thatmodelclass. in part by Yoshi Oono. During the course of this work, R.
The demonstrated examples would perhaps not arise if the Biswas was partially supported by NSF AST 07-08849 and
datawasgoodenoughthatthedifferencesinχ2MLwaslarge. NSF AST 09-08693 ARRA and then by Argonne National
Laboratory,operated byUChicago Argonne,LLC. Workat
We also calculated the respective statistical uncertainty Argonnewas supported byU.S.Departmentof Energy,Of-
(∼ 1 σ error bars) of ∆dA,B and showed it was even larger ficeofScienceundercontractDE-AC02-06CH11357.Finally,
than the observed ∆A,B between some of the models. This we would like to express gratitude to the anonymous re-
gives us a way to gauge the adequacy in the number of viewerwhoseinputsubstantially improvedthetoneand di-
data points since the statistical uncertainty would become rection of this paper.
smaller than the observed differences when there is a
sufficient amount of data. As an important example, we
considered the ΛCDM and DGP models since they were
shown to have the two lowest AIC values in Table. 1. It APPENDIX A: CONFUSION PROBABILITY
was shown that the statistical uncertainty of ∆dΛCDM,DGP AND MODEL LIKELIHOOD
was larger than the observed ∆ΛCDM,DGP, making it
difficult to determine the better model between the two. This is an almost identical repeat of Balasubramanian
From our simulation, we also showed that the shapes of (1997) explanation of error probabilities which is framed in
the distribution of the AIC differences can be quite varied, the language of hypothesis testing. It is reproduced for the
ranging from a symmetric Gaussian-like distribution to an convenience of the reader. Suppose {x1,x2,...,xN} ∈ XN