Objective Bayesian point and region estimation in location-scale models

Statistics & Operations Research Transactions SORT 31 (1) January-June 2007, 3-44 Statistics & Operations Research Transactions Objective Bayesian point and region estimation in location-scale models c Institut d’Estad´ ıstica de Catalunya [email protected] ISSN: 1696-2281 www.idescat.net/sort Jos´ eM.Bernardo Universitat de Val`encia, Spain Abstract Point and region estimation may both be described as specific decision problems . In point estimation, the action space is the set of possible values of the quantity on interest; in region estimation, the action space is the set of its possible credible regions. Foundations dictate that the solution to these decision problems must depend on both the utility function and the prior distribution. Estimators intended for general use should surely be invariant under one-to-one transformations, and this requires the use of an invariant loss function; moreover, an objective solution requires the use of a prior which does not introduce subjective elements. The combined use of an invariant information-theory based loss function, the intrinsic discrepancy , and an objective prior, the reference prior , produces a general solution to both point and region estimation problems. In this paper, estimation of the two parameters of univariate location-scale models is considered in detail from this point of view, with special attention to the normal model. The solutions found are compared with a range of conventional solutions. MSC: Primary: 62F15, 62C10; secondary: 62F10, 62F15, 62B10 Keywords: Confidence Intervals, Credible Regions, Decision Theory, Intrinsic Discrepancy, Intrinsic Loss, Location-Scale Models, Noninformative Prior, Reference Analysis, Region Estimation, Point Estimation. 1 Introduction Point and region estimation of the parameters of location-scale models have a long, fascinating history which is far from settled. Indeed, the list of contributors to the simpler examples of this class of problems, estimation of the normal mean and estimation of the normal variance, reads like a Who’s Who in 20th century statistics. Address for correspondence: Jos´ e M. Bernardo is Professor of Statistics at the Universitat de Val` encia. Departamento de Estad´ ıstica e I. O., Facultad de Matem` aticas, 46100-Burjassot, Valencia, Spain. [email protected], www.uv.es/bernardo Received: March 2006 4 Objective Bayesian point and region estimation in location-scale models In this paper, an objective Bayesian decision-theoretic solution to both point and region estimation of the parameters of location-scale models is presented, with special attention devoted to the normal model. In marked contrast with most approaches, the solutions found are invariant under one-to-one reparametrization. 1.1 Notation Probability distributions are described through their probability density functions, and no notational distinction is made between a random quantity and the particular values that it may take. Bold italic roman fonts are used for observable random vectors (typically data) and bold italic greek fonts for unobservable random vectors (typically parameters); lower case is used for variables and upper case calligraphic for their dominion sets. The standard mathematical convention of referring to functions, say fx(·) and gx(·)ofx∈X, respectively by f(x)andg(x) is often used. Thus, the conditional probability density of observable data x∈Xgiven ω ω ωis represented by either px(·|ω ω ω) or p(x|ω ω ω), with p(x|ω ω ω)≥0, x∈X,andXp(x|ω ω ω)dx=1, and the posterior density of a non-observable parameter vector θ θ θ∈Θgiven data xis represented by either πθ θ θ(·|x)orπ(θ θ θ|x), with π(θ θ θ|x)≥0andΘπ(θ θ θ|x)dθ θ θ=1. Density functions of specific distributions are denoted by appropriate names. In particular, if xhas a normal distribution with mean μand standard deviation σ, its probability density function will be denoted N(x|μ, σ), and if λhas a gamma distribution with parameters αand β, its probability density function will be denoted Ga(λ|α, β), with E[λ]=α/β,and Var[λ]=α/β2. It is assumed that available data xconsist of one observation from the family F≡{p(x|ω ω ω),x∈X,ω ω ω∈Ω}of probability distributions for x∈X,andthat one is interested in point and region estimation of some function θ θ θ=θ θ θ(ω ω ω)∈Θof the unknown parameter vector ω ω ω. Often, but not necessarily, data consist of a random sample x={x1,...,xn}of some simpler model {q(x|ω ω ω),x∈X,ω ω ω∈Ω},inwhich case, X=Xnand the likelihood function is p(x|ω ω ω)=n j=1q(xj|ω ω ω). Without loss of generality, the original parametric family Fmay be written as F≡{px(·|θ θ θ, λ λ λ),x∈X,θ θ θ∈Θ,λ λ λ∈Λ}(1) in terms of the vector of interest θ θ θ, and a vector λ λ λof nuisance parameters. A point estimator of θ θ θis some function of the data ˜ θ θ θ(x)∈Θsuch that, for each possible set of observed data x,˜ θ θ θ(x) could be regarded as an appropriate proxy for the actual, unknown value of θ θ θ.Ap-credible region of θ θ θis some subset Cp(x,Θ)ofΘwhose posterior probability is p. Within this framework, attention in this paper focuses on problems where data consist of a random sample x={x1,...,xn}from a location-scale model J. M. Bernardo 5 m(x|μ, σ, f), of the form m(x|μ, σ, f)=σ−1f{σ−1(x−μ)},x∈,μ∈,σ>0,(2) where f(.) is some probability density in ,sothat f(y)≥0, f(y)dy =1. Interest lies in either the location parameter μ, the scale parameter σ, or some one-to-one function of these, and the likelihood function is p(x|μ, σ)=n j=1m(xj|μ, σ, f)=σ−nn j=1f{σ−1(xj−μ)}.(3) Standard notation is used for the sample mean and the sample variance, respectively denoted by x=n j=1xj/nand s2=n j=1(xj−x)2/n. Many conventional point estimators of the variance of location-scale models are members of the family of affine invariant estimators, ˜σ2 ν=ns2 ν=1 νn j=1(xj−x)2,ν>0.(4) In particular, with normal data, the MLE of the variance σ2is s2=˜σ2 n, and the unbiased estimator is ˜σ2 n−1. More sophisticated estimators may sometimes be defined in terms of affine estimators; for instance, Stein (1964) and Brown (1968) estimators of the normal variance may respectively be written as ˜σ2 stein =min ˜σ2 n+1,˜σ2 (n+2)/(1+z2),˜σ2 brown =min ˜σ2 n−1,˜σ2 n/(1+z2), where z=x/sis the standardized sample mean. 1.2 Contents Section 2 provides a short review of intrinsic estimation, our approach to both point and region estimation. An information-theory based invariant loss function, the intrinsic discrepancy, is proposed as a reasonable general alternative to the conventional (noninvariant) quadratic loss. As is usually the case in modern literature, point estimation is described as a decision problem where the action space is the set of possible values for the quantity of interest; an intrinsic point estimator is then defined as the Bayes estimator which corresponds to the intrinsic loss and the appropriate reference prior. This provides a general invariant objective Bayes point estimator. Less conventionally, region estimation is also described as a decision problem where, for each p, the action space is the set of possible p-credible regions for the quantity of interest; a p-credible intrinsic region estimator is then defined as the lowest posterior loss p-credible region with respect to the intrinsic loss and the appropriate reference prior. This provides a 6 Objective Bayesian point and region estimation in location-scale models general invariant objective Bayes region estimator which always contains the intrinsic point estimator. In Section 3 location-scale models are analyzed from this point of view. In particular, intrinsic point estimators and intrinsic region estimators are derived for the mean of a normal model, the variance of a normal model, and the scale parameter of a Cauchy model. 2 Intrinsic Estimation 2.1 The intrinsic discrepancy loss function Point estimation of some parameter vector θ θ θ∈Θis customarily described as a decision problem where the action space is the set A={˜ θ θ θ;˜ θ θ θ∈Θ}of possible values of the vector of interest. Foundations dictate (see e.g., Bernardo and Smith, 1994, Ch. 2 and references therein) that to solve this decision problem it is necessary to specify a loss function {˜ θ θ θ, θ θ θ}, such that {˜ θ θ θ, θ θ θ}≥0and{θ θ θ, θ θ θ}=0, which describes, as a function of θ θ θ, the loss suffered from using ˜ θ θ θas a proxy for the unknown value of θ θ θ.Theloss function is context specific, and should be chosen in terms of the anticipated uses of the estimate; however, a number of conventional loss functions have been suggested for those situations where no particular uses are envisaged, as in scientific inference. The simplest of these conventional loss functions (which typically ignore the presence nuisance parameters) is the ubiquitous quadratic loss,{˜ θ θ θ, (θ θ θ, λ λ λ)}=(˜ θ θ θ−θ θ θ)t(˜ θ θ θ−θ θ θ); the corresponding Bayes estimator, if this exists, is the posterior mean,E[θ θ θ|x]. Another common conventional loss function is the zero-one loss,defined as {˜ θ θ θ, (θ θ θ, λ λ λ)}=1, if ˜ θ θ θ does not belong to a -radius neighbourhood of θ θ θ, and zero otherwise; as →0, the corresponding Bayes estimator converges to the posterior mode, Mo[θ θ θ|x]. For details, see, e.g., Bernardo and Smith (1994, p. 257). Example 1 (Normal variance) With the usual objective prior π(μ, σ)=σ−1,the (marginal) reference posterior density of σis the square root inverted gamma π(σ|x)=π(σ|s,n)=n(n−1)/2sn−1 2(n−3)/2Γ[(n−1)/2] σ−ne−1 2ns 2/σ2,n≥2.(5) The quadratic loss in terms of the variance, {˜σ2,σ 2}=(˜σ2−σ2)2, leads to E[σ2|x]= ˜σ2 n−3(which obviously requires n≥3). Similarly, the quadratic loss in terms of the standard deviation, {˜σ, σ}=(˜σ−σ)2, yields E[σ|x]=n 2 Γ[(n−2)/2] Γ[(n−1)/2] s,n≥3.(6) J. M. Bernardo 7 For moderate nvalues, the Stirling approximation of the Gamma functions in (6) produces E[σ|x]2≈˜σ2 n−5/2. Notice that, using the conventional quadratic loss function, the Bayes estimate of σ2is not the square of the Bayes estimate of σ. This lack of invariance is not an special feature of the quadratic loss; on the contrary, this is the case of most conventional loss functions. For instance, the use of the slightly more sophisticated standardized quadratic loss function on the variance, stq(˜σ2,σ 2)=[( ˜σ2/σ2)−1]2(7) yields (if n≥2, for π(σ|x) to be proper) arg min ˜σ2>0∞ 0 [( ˜σ2/σ2)−1]2π(σ|x)dσ=ns 2 n+1=˜σ2 n+1,(8) which is also the minimum risk equivariant estimator (MRIE) of σ2under this loss, while the standardized quadratic loss function it terms of the standard deviation, (˜σ, σ)=[( ˜σ/σ)−1]2yields (again, if n≥2) arg min ˜σ2>0∞ 0 [( ˜σ/σ)−1]2π(σ|x)dσ=n 2 Γ[n/2] Γ[(n+1)/2] s,(9) which is different from (6), and whose square is not (8). Similarly, for the zero-one loss in terms of σ2, the Bayes estimator is the mode of the posterior distribution of σ2, π(σ2|x)=π(σ|x)/(2σ), which is Mo(σ2|x)=˜σ2 n+1, the same as (8), while the Bayes estimator for the zero-one loss in terms of σis Mo(σ|x)=s,theMLEofσ, whose square is obviously not the same as (8). For further information on alternative point estimators of the normal variance, see Brewster and Zidek (1974) and Rukhin (1987). As Example 1 dramatically illustrates, conventional loss functions are typically not invariant under reparametrization. As a consequence, the Bayes estimator φ φ φ∗of a oneto-one transformation φ φ φ=φ φ φ(θ θ θ) of the original parameter θ θ θis not necessarily φ φ φ(θ θ θ∗) and thus, for each loss function, one may produce as many different estimators of the same quantity as alternative parametrizations one is prepared to consider, a less than satisfactory situation. Indeed, scientific applications require this type of invariance. It would certainly be hard to argue that the best estimate of, say the age of the universe is θ∗but that the best estimate of the logarithm of that age is not log(θ∗). Invariant loss functions are required to guarantee invariant estimators. With no nuisance parameters, intrinsic loss functions (Robert, 1996), of the general form (˜ θ θ θ, θ θ θ)={px(.|˜ θ θ θ),px(.|θ θ θ)}shift attention from the discrepancy between the estimate ˜ θ θ θand the true value θ θ θ, to the more relevant discrepancy between the statistical models they label, and they are always invariant under one-to-one reparametrization. The intrinsic discrepancy, introduced by Bernardo and Rueda (2002), is a particular intrinsic loss with specially attractive properties. 8 Objective Bayesian point and region estimation in location-scale models Definition 1 (Intrinsic Discrepancy) The intrinsic discrepancy between two elements px(·|ω ω ω1)and px(·|ω ω ω2)of the parametric family of distributions F={px(·|ω ω ω),x∈ X(ω ω ω),ω ω ω∈Ω},is δx(ω ω ω1,ω ω ω2)=δ{px(.|ω ω ω1),px(.|ω ω ω2)}=min{κx(ω ω ω1|ω ω ω2),κ x(ω ω ω2|ω ω ω1)}, κx(ω ω ωj|ω ω ωi)=X(ω ω ωi) px(x|ω ω ωi)log px(x|ω ω ωi) px(x|ω ω ωj)dx, The intrinsic discrepancy δx{F1,F2}between two subsets F1and F2of Fis the minimum intrinsic discrepancy between its elements, δx(F1,F2)=min ω ω ω1∈F1,ω ω ω2∈F2 δ{px(.|ω ω ω1),px(.|ω ω ω2)} Thus, the intrinsic discrepancy δ(ω ω ω1,ω ω ω2) between two parameter values ω ω ω1and ω ω ω2is the minimum Kullback-Leibler directed logarithmic divergence (Kullback and Leibler, 1951) between the distributions px(.|ω ω ω1)andpx(.|ω ω ω2) which they label. Notice that this is obviously independent of the particular parametrization chosen to describe the distributions. The intrinsic discrepancy is a divergence measure in the class F; indeed, (i) it is symmetric, (ii) it is non-negative and (iii) it is zero if, and only if, px(x|ω ω ω1)= px(x|ω ω ω2) almost everywhere. Notice that in Definition 1 the possible dependence of the sampling space X=X(ω ω ω) on the parameter value ω ω ωis explicitly allowed, so that the intrinsic discrepancy may be used with non-regular models where the support X(ω ω ω1) of, say, px(.|ω ω ω1) may be strictly smaller than the support X(ω ω ω2)ofpx(.|ω ω ω2). The intrinsic discrepancy is also invariant under one-to-one transformations of the random vector x. Moreover, directed logarithmic divergences are additive with respect to conditionally independent observations. Consequently, if x={x1,...,xn}is a random sample from, say qx(·|ω ω ω) so that the probability model is p(x|ω ω ω)=n i=1q(xj|ω ω ω), then the intrinsic discrepancy δx{ω ω ω1,ω ω ω2}between px(·|ω ω ω1)andpx(·|ω ω ω2)issimply nδx{ω ω ω1,ω ω ω2},thatis,ntimes the intrinsic discrepancy between qx(·|ω ω ω1)andqx(·|ω ω ω2). In the context of point estimation, the intrinsic discrepancy leads naturally to the (invariant) intrinsic discrepancy loss δx{˜ θ θ θ, (θ θ θ, λ λ λ)}defined as the intrinsic discrepancy between the assumed model px(·|θ θ θ, λ λ λ) and its closest approximation within the set {px(·|˜ θ θ θ, ˜ λ λ λ),˜ λ λ λ∈Λ}of all models with θ θ θ=˜ θ θ θ. Definition 2 (Intrinsic discrepancy loss) Consider the family of probability distributions F={px(·|θ θ θ, λ λ λ),θ θ θ∈Θ,λ λ λ∈Λ,x∈X(ω ω ω, λ λ λ)}.Theintrinsic discrepancy loss from using ˜ θ θ θas a proxy for θ θ θis δx{˜ θ θ θ, (θ θ θ, λ λ λ)}=inf ˜ λ λ λ∈Λ δx{(˜ θ θ θ, ˜ λ λ λ),(θ θ θ, λ λ λ)}, the intrinsic discrepancy between px(·|θ θ θ, λ λ λ)and the set {px(·|˜ θ θ θ, ˜ λ λ λ),˜ λ λ λ∈Λ}. J. M. Bernardo 9 Notice that the value of δx{˜ θ θ θ, (θ θ θ, λ λ λ)}does not depend on the particular parametrization chosen to describe the problem. Indeed, for any one-to-one reparametrizations φ φ φ=φ φ φ(θ θ θ) and ψ ψ ψ=ψ ψ ψ(λ λ λ), δx{˜ φ φ φ, (φ φ φ, ψ ψ ψ)}=δx{˜ θ θ θ, (θ θ θ, λ λ λ)}(10) so that, as one should surely require, the loss suffered from using ˜ φ φ φ=φ φ φ(˜ θ θ θ) as a proxy for φ φ φ(θ θ θ) is precisely the same as the loss suffered from using ˜ θ θ θas a proxy for θ θ θ, and this is true for any parametrization of the nuisance parameter vector. Under frequently met regularity conditions, the two minimizations required in Definition 2 may be interchanged. This makes analytical derivation of the intrinsic loss considerably simpler. Theorem 1 (Computation of the intrinsic discrepancy loss) Let Fbe a parametric family of probability distributions F={p(x|θ θ θ, λ λ λ),θ θ θ∈Θ,λ λ λ∈Λ,x∈X(θ θ θ, λ λ λ)}, with convex support X(θ θ θ, λ λ λ)for all θ θ θand λ λ λ. Then, δx{˜ θ θ θ, (θ θ θ, λ λ λ)}=inf ˜ λ λ λ∈Λ min κx{˜ θ θ θ, ˜ λ λ λ|θ θ θ, λ λ λ},κ x{θ θ θ, λ λ λ|˜ θ θ θ, ˜ λ λ λ} =min inf ˜ λ λ λ∈Λ κx{˜ θ θ θ, ˜ λ λ λ|θ θ θ, λ λ λ},inf ˜ λ λ λ∈Λ κx{θ θ θ, λ λ λ|˜ θ θ θ, ˜ λ λ λ} Proof. This follows from the fact that, if X(θ θ θ, λ λ λ) is a convex set, then the two directed logarithmic divergences involved in the definition are convex functionals. For details, see Ju´ arez (2004).  Example 2 (Normal variance, continued) Consider x={x1,...,xn}. The directed logarithmic divergence of n i=1N(xj|˜μ, ˜σ) from n i=1N(xj|μ, σ)is κx{˜μ, ˜σ|μ, σ}=n∞ −∞ N(x|μ, σ)logN(x|μ, σ) N(x|˜μ, ˜σ)dx =n 2σ2 ˜σ2−1−log σ2 ˜σ2+(μ−˜μ)2 ˜σ2.(11) This is minimized when ˜μ=μ, to yield inf ˜μ∈Rκx{˜μ, ˜σ|μ, σ}=n 2gσ2 ˜σ2=n 2g(φ), 10 Objective Bayesian point and region estimation in location-scale models where φ=˜σ2/σ2,andg(.)isthelinlog function defined by g(t)=(t−1) −log t,t>0.(12) Notice that g(t)≥0andg(t)=0 if (and only if) t=1. This follows from the fact that g(t) is the absolute distance between log tand its tangent at t=1. Exchanging the roles of (˜μ, ˜σ)and(μ, σ), it is similarly found that κx{μ, σ |˜μ, ˜σ}is also minimized when ˜μ=μ, to yield inf ˜μ∈Rκx{μ, σ |˜μ, ˜σ}=n 2g˜σ2 σ2=n 2g1 φ. Moreover, g(t)<g(1/t) if, and only if, t<1 and hence, using Theorem 1, δx{˜σ, (μ, σ)}=δx{φ}=⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ n 2g(φ)ifφ<1, n 2g(1/φ)ifφ≥1,φ=˜σ2 σ2.(13) Thus, for fixed n, the intrinsic discrepancy loss δx{˜σ, (μ, σ)}only depends on the ratio φ=˜σ2/σ2. The intrinsic discrepancy loss is closely related to the (also invariant) entropy loss, ent{˜σ, σ}=ent{φ}=∞ −∞ N(x|μ, σ)logN(x|μ, σ) N(x|μ, ˜σ)dx =1 2g(φ),(14) 42 0 2 4 1 1 Δstqent logΣ 2Σ2 Figure 1: Intrinsic discrepancy loss (solid line), entropy loss (continuous line), and standardized quadratic loss (dotted line) for point estimation of the normal variance, as a function of ψ=log( ˜σ2/σ2). J. M. Bernardo 11 which Brown (1990) attributes to Stein. Except for the proportionality constant n(which does not affect estimation), the entropy loss (14) is the same as the intrinsic discrepancy loss (13) whenever ˜σ<σ. Indeed, the intrinsic discrepancy loss may be seen as a symmetrized version of the entropy loss. Notice that, for all values of the ratio φ=˜σ2/σ2,δx{φ}=δx{1/φ}; hence, the intrinsic loss equally penalizes overestimation and underestimation. In sharp contrast, both the entropy loss and the often recommended standardized quadratic loss function, which is also a function of the ratio φ, stq(˜σ2,σ 2)=[( ˜σ2/σ2)−1]2=(φ−1)2, clearly underpenalize small estimators, thus yielding estimators of the variance which are too small. This is illustrated in Figure 1, where the functions δx(φ) (for n=1), ent{φ},andstq(φ) are all represented as a function of ϕ=log φ. More conventional loss functions, as the usual quadratic loss, quad(˜σ2,σ 2)=[˜σ2−σ2]2=σ4(φ−1)2, are not even invariant with respect to affine transformations. All this led Stein (1964, p. 156) to write “I find it hard to take the problem of estimating σ2with quadratic loss very seriously”. 2.2 Reference posterior expectation Given data xgenerated by p(x|θ θ θ, λ λ λ), a situation with no prior information about the value of θ θ θis formally described by the reference prior π(λ λ λ|θ θ θ)π(θ θ θ) which corresponds to the model px(·|θ θ θ, λ λ λ)whenθ θ θis the quantity of interest (Bernardo, 1979; Berger and Bernardo, 1992; Bernardo, 2005a). In this case, all available information about θ θ θis encapsulated it its (marginal) reference posterior distribution, π(θ θ θ|x)=Λπ(θ θ θ, λ λ λ|x)dλ λ λ where, by Bayes theorem, π(θ θ θ, λ λ λ|x)∝p(x|θ θ θ, ω ω ω)π(λ λ λ|θ θ θ)π(θ θ θ).If numerical summaries of the information encapsulated in π(θ θ θ|x) are further required in the form of either point or region estimators of θ θ θunder some specified loss function {˜ θ θ θ, (θ θ θ, λ λ λ)}, then the reference posterior expected loss l(˜ θ θ θ|x)=ΘΩ {˜ θ θ θ, (θ θ θ, λ λ λ)}π(θ θ θ, λ λ λ|x)dθ θ θdλ λ λ (from using ˜ θ θ θas a proxy for θ θ θ) has to be evaluated. In view of the arguments given above, attention will focus on the intrinsic discrepancy reference expected loss, or intrinsic expected loss, for short. 18 Objective Bayesian point and region estimation in location-scale models δ{˜μ, (μ, σ)}=δ{θ2}=n 2log[1 +θ2 n],θ=θ(˜μ, μ, σ)=μ−˜μ σ/ √n ,(26) which only depends on the number θof standard deviations which separate ˜μfrom μ. Figure 3 represents the intrinsic loss function (26), as a function of θ, for several values of n. As one might expect, δ{˜μ, (μ, σ)}increases with |θ|. The dependence is essentially quadratic in a neighbourhood of zero, but shows a very reasonable concavity in regions where |θ|is large. Using Definition 3, the intrinsic discrepancy reference expected loss d(˜μ|x) may be written in terms of the reference posterior of θ; indeed, d(˜μ|x)=∞ 0∞ −∞ δ{˜μ, (μ, σ)}π(μ, σ |x)dμdσ =∞ 0 n 2log 1+θ2 nπθ|xdθ. (27) But θ=(˜μ−μ)/(σ/ √n) may be written as a+βwhere, as a function of μand σ, β=(μ−x)/(σ/ √n) has a standard normal reference posterior, and ais the constant a=(x−˜μ)/(σ/ √n). Hence, the conditional posterior distribution of θ2given σis noncentral χ2with one degree of freedom and non centrality parameter a2, π(θ2|x,σ)=χ2(θ2|1,a2),a2=n(x−˜μ)2 σ2. 40 20 0 20 40 2 4 6 8 Θ n2 n3 n10 Figure 3: Intrinsic discrepancy loss for estimation of the normal mean as a function of the number θ=(˜μ−μ)/(σ/ √n)of standard deviations which separates ˜μfrom μ,forn=2,n=3, and n =10. J. M. Bernardo 19 It follows that the intrinsic expected loss d(˜μ|x) only depends on ˜μthrough (x−˜μ)2,and increases with (x−˜μ)2; therefore, the intrinsic estimator of μis ˜μint(x)=arg min ˜μ∈ d(˜μ|x)=arg min ˜μ∈ (x−˜μ)2=x. Moreover, d(˜μ|x) is symmetric around xand, hence, all intrinsic credible regions must be centered at x. In view of (25), this implies that the intrinsic the p-credible regions are just the usual Student-tHPD p-credible intervals Cint p(x)=˜μ;x−qp,ns/√n−1≤˜μ≤x+qp,ns/√n−1,(28) where qp,nis the (p+1)/2 quantile of a standard Student-twith n−1 degrees of freedom. It immediately follows from (28) that Cint pconsist of the set of ˜μvalues such that (x−˜μ)/(s/√n−1) belongs to a probability pcentred interval of a standard Student-t with n−1 degrees of freedom. But, as a function of the data x, the sampling distribution of t(x)=(x−μ)/(s/√n−1) (29) is also a standard Student-twith n−1 degrees of freedom. Hence, for all sample sizes, the expected coverage under sampling of the p-credible intervals (28) is exactly p,and the intrinsic credible regions are exact frequentist confidence intervals. A simple asymptotic approximation to d(˜μ|x), which provides a direct measure in a log-likelihood ratio scale of the expected loss associated to the use of ˜μ, may easily be obtained. Indeed, a variation of the delta method shows that, under appropriate regularity conditions, the expectation of some function y=g(x) of a random quantity xwith mean μxand variance σ2 xmay be approximated by E[g(x)] ≈gμx+σ2 x 2 g(μx) g(μx).(30) On the other hand, the conditional posterior mean of θ2is 1 +a2, and its conditional posterior variance is 2+4a2;butE[σ−2|x]=E[λ|x]=(n−1)/(ns2) (Eq. 24) and hence, the unconditional posterior mean and variance of θ2(˜μ) are, respectively, E[θ2|x]=1+t2,Var[θ2|x]=2+4t2, both functions of the conventional tstatistic (29). Using these in (30) to approximate the posterior expectation of log(1 +θ2/n) required in (27) yields d(˜μ|x)≈n 2log 1+1 n n(1 +t2)+t4 n+t2+1.(31) 20 Objective Bayesian point and region estimation in location-scale models Progressively cruder, but simpler approximations are d(˜μ|x)≈n 2log 1+1 n1+t2≈1 21+t2.(32) Thus, for large n, the intrinsic expected loss d(˜μ|x) is essentially quadratic in the number t=(x−˜μ)/(s/√n−1) of standard deviations which separate xfrom ˜μ. Summarizing, we have thus established Theorem 2 (Intrinsic estimation of the Normal mean) Let xbe a random sample of size n from N(x|μ, σ), with mean and variance x and s2, and let t =√n−1(x−˜μ)/s) be the conventional t statistic. (i) The intrinsic point estimator of μis ˜μint(x)=x. (ii) The unique p-credible intrinsic region for μis the probability centred interval Cint p(x)=x±qp,ns/√n−1, where qp,nis the (p+1)/2quantile of a standard Student-t distribution with n −1 degrees of freedom. For all sample sizes, the frequentist coverage of Cint p(x)is exactly p. (iii) The expected intrinsic loss associated to the use of ˜μas a proxy for μis d(˜μ|x)≈n 2log 1+1 n n(t2+1) +t4 n+t2+1≈n 2log 1+1 n1+t2. As a numerical illustration, a random sample of size n=25 was generated from a standard normal, yielding x=−0.162 and s=0.840. The intrinsic estimator is μ∗= x=−0.162 and the 0.99-intrinsic credible region is the interval [−0.642,0.318]. The exact value of the expected intrinsic loss d(1/3|x), computed from (27) by numerical integration, is 3.768, while (31) and the two approximations in (32) respectively yield 3.781, 3.970 and 4.673. Hence, the observed data may be expected to be about exp(3.768) ≈43 times more likely under the true value of μthat under the closest normal model with μ=1/3, suggesting that the value μ=1/3 is hardly compatible with the observed data. 3.3 Intrinsic estimation of the normal variance It has already been established (Example 2, Eq. 13) that the intrinsic discrepancy loss from using ˜σ2as a proxy for σ2is δx{˜σ2,(μ, σ)}=δx{φ}=⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ n 2g(φ)ifφ<1, n 2g(1/φ)ifφ≥1, (33) J. M. Bernardo 21 where g(t)=(t−1) −log t,andφ=˜σ2/σ2, and that this is also the intrinsic loss δx{˜ ψ, (μ, σ)}from using ˜ ψa a proxy for ψfor any one-to-one function ψ(σ2)ofσ2. Moreover, the reference posterior of φis the gamma distribution of Eq. 18. Hence, the intrinsic estimator of σ2is ˜σ2 int(x)=arg min ˜ σ2>0∞ 0 δx(φ)Ga φn−1 2,ns 2 2˜σ2dφ, where δx(φ) is given by (33). Moreover, it immediately follows from (18) that, as a function of σ, the reference posterior distribution of τ=ns 2/σ2is π(τ|x)=π(τ|n)=χ2(τ|n−1) (34) a central χ2with n−1 degrees of freedom; but φ=cτ/n, with c=˜σ2/s2and, therefore, the expected posterior loss from using ˜σmay further be written as d(˜σ2|s2,n)=d(c|n)=∞ 0 δcτ nχ2(τ|n−1) dτ, c=˜σ2/s2.(35) Thus, the intrinsic estimator of the normal variance is an affine equivariant estimator of the form σ2 int(s,n)=c∗ ns2,c∗ n>0,(36) where c∗ nis the value of cwhich minimizes d(c|n) in (35). The exact value of c∗ nmay be numerically found by one-dimensional numerical integration, followed by numerical optimization. The first row of Table 1 displays the exact values of c∗ nfor several sample sizes. However, good analytical approximations for c∗ nmay be obtained. We first consider a general approximation method. Let ωbe a particular parametrization of the problem, and consider a (variance stabilizing) reference reparametrization Table 1: Exact and alternative approximate values for the intrinsic point estimator of the normal variance σ2 int =c∗ ns2. n2 3 4 5 10 20 50 c∗ n4.982 2.347 1.803 1.569 1.231 1.106 1.041 n n−12 4.000 2.250 1.778 1.563 1.235 1.108 1.041 n n−1e1/(n−1) 5.437 2.473 1.861 1.605 1.242 1.110 1.041 n n−2— 3.000 2.000 1.667 1.250 1.111 1.042 22 Objective Bayesian point and region estimation in location-scale models φ(ω)defined as one with a uniform reference prior. This is given by any solution to the differential equation φ(ω)=π(ω), where π(ω) is the marginal reference prior for ω. Under regularity conditions, the sampling distribution of φ(ˆω), where ˆω=ˆω(x)isthe MLE of ω, and the reference posterior of φ(ω), are both asymptotically normal. Using these approximations, the intrinsic expected loss from using ˜ωis found to be (Bernardo, 2005b, Theo. 4.1) d(˜ω|x)≈n 2σ2 φ+[μφ−φ(˜ω)] 2.(37) where μφand σ2 φare respectively the posterior mean and posterior variance of φ=φ(ω). This is minimized by φ(˜ω)=μφ=E[φ|x]. Hence, in terms of any reference parametrization φ, the intrinsic point estimate is approximately the posterior mean μφ and, by invariance, the intrinsic estimator of any one-to-one function, ψ=ψ(φ)is approximately given by ˜ ψint =ψ(μφ). Thus, ˜ φint(x)≈μφ=E[φ|x],˜ωint(x)≈φ−1{μφ}.(38) Under regularity conditions (see, e.g., Schervish, 1995, Sec. 7.1.3) the delta method may be used to obtain simple approximations to the posterior moments of φis terms of those of ω, namely μφ≈φ{μω}+σ2 ωφ{μω}/2,(39) σ2 φ≈σ2 ω[φ{μω}]2.(40) Substitution into (38) and (37) respectively provide useful approximations to the intrinsic point estimator of ω, and to the expected loss from using ˜ωas a proxy for ω. In the particular case of the normal variance, it is convenient to start from the parametrization in terms of the precision λ=σ−2, whose posterior moments have simple expressions. Since reference priors are consistent under reparametrization, the reference prior for λis π(λ)=π(σ)|∂σ/∂λ|∝λ−1and, therefore, a reference parametrization is φ=φ(λ)=λ π(λ)dλ=λ λ−1dλ=log λ. Notice that the reference prior of φ=log λis indeed uniform, as it is the case for the logarithm of any other power of σ. Using (39) and (40) with the first posterior moments of λ, given in (24), yields ˜ φint(x)≈μφ=E[log λ|x]≈log n−1 ns 2−1 n−1,(41) σ2 φ=Var[log λ|x]≈2 n−1.(42) J. M. Bernardo 23 By invariance, (41) directly provides an approximation to the intrinsic estimator of the variance. This has the form of a modified version of the conventional unbiased estimator ˜σ2 n−1; indeed, since σ2=e−φ, ˜σ2 int(x)=e−˜ φint(x)≈ns 2 n−1e1 n−1=˜σ2 n−1e1 n−1, which, as shown in the third row of Table 1, provides good approximations, even for small values of n. A better analytical approximation to the intrinsic estimator of the normal variance may be obtained making use of the particular features of this example. This is done by separately minimizing the expected value of each of the two functions which enter the definition of the intrinsic discrepancy loss δx{φ}, and using the arithmetic mean of the corresponding results. Indeed, the delta method may be used to approximate both E[g(cτ/n)] and E[g(n/(cτ))] in terms E[τ|n]=n−1 and Var[τ|n]=2(n−1). The approximation to E[g(cτ/n)] is minimized by ˆc∗ 1n=n/(n−1), while the approximation to E[g(n/(cτ))] is minimized by by ˆc∗ 2n=n(n+1)/(n−1)2. As one would expect, their average, ˆc∗ n=ˆc∗ 1n+ˆc∗ 2n 2=n n−12=n n−(2 −n−1),(43) provides a good approximation to the value c∗ nwhich minimizes (35). As shown in the second row of Table 1, the approximation remains good even for small values of n. Combination of (36) and (43) establishes that, for all but very small nvalues, ˜σ2 int(x)=˜σ2 int(s,n)≈n n−12s2.(44) In view of the second expression for ˆc∗ nin (43), a cruder approximation is given by ˜σ2 int ≈˜σ2 n−2. This is larger than the MLE ˆσ2=s2(which divides by nthe sum of squares), and also larger than the conventional unbiased estimate of the variance ˜σ2 n−1(which divides by n−1). Notice that numerical differences between intrinsic and conventional estimators may be large for small values of n. In particular, with only two observations {x1,x2}, the intrinsic estimator of the variance is σ2 int(2,s2)≈5s2=5(x1−x2)2/4; this is 2.5 times larger than the unbiased estimator, (x1−x2)2/2 in this case, which (with good reason) is generally considered to be too small. As shown by (35), the expected intrinsic loss d(˜σ2|n,s2)ofanyaffine equivariant estimator of the variance ˜σ2=kns2, is actually independent of s2and only depends on the sample size n. Moreover, it is easily verified that the expected intrinsic loss d(˜σ2|n,s2) is precisely equal to the frequentist risk associated to the intrinsic discrepancy loss, r(˜σ2 i|n,σ 2)=∞ 0 δ{˜σ2 i(s2),σ 2}p(s2|n,σ 2)ds2. 24 Objective Bayesian point and region estimation in location-scale models Thus, under intrinsic discrepancy loss, the intrinsic estimator ˜σ2 int dominates all affine equivariant estimators. For details, see Bernardo (2006). Region estimation is now considered. As described in Example 4, the intrinsic pcredible region for σis the unique solution Cint p={σ0,σ 1}to the equations system d(σ2 0|x)=d(σ2 1|x),σ1 σ0 π(σ|x)dσ=p. Using (34), this may equivalently be written in terms of τ=ns 2/σ2as d(σ2 0|x)=d(σ2 1|x),ns2/σ2 0 ns2/σ1 χ(τ|n−1) dτ=p.(45) Thus, the unique p-credible intrinsic region for σ2is the interval Cint p(x)=ns 2 χ2 n−1(1 −α),ns 2 χ2 n−1(1 −p−α)(46) where χ2 n−1(q)istheqquantile of a χ2 n−1distribution, and αis the solution to the equation dns 2 χ2 n−1(1 −α)|x=dns 2 χ2 n−1(1 −p−α)|x.(47) By invariance, this provides the intrinsic p-credible region of any one-to-one function of σ2. As a function of the data x, the sampling distribution of ns 2/σ2is also a χ2with n−1 degrees of freedom. Hence, for all sample sizes, the expected coverage under sampling of the p-credible intervals (46) is exactly p. Using (35) to evaluate expected losses, the exact solution to equation (47) may easily be obtained by numerical methods. However, good analytical approximations are may be obtained. Working again in terms of the reference parametrization for this problem, φ= log λ=−2logσ, and using (37), (39) and (40), the expected loss from using ˜ φas a proxy for φis approximately d(˜ φ|x)≈n 22 n−1+˜ φint −˜ φ2.(48) But this is symmetric around ˜ φint =log(˜ λint)=−log(σ2 int) and therefore, to keep those ˜ φpoints with smaller expected loss, any intrinsic credible region for φ=log λmust be (approximately) centered at ˜ φint. Thus, using (42) and (44) this will be of the form J. M. Bernardo 25 Cint p(x,Φ)≈˜ φint ±αpnσφ≈log n−1 n21 s2±γpn 2 n−1(49) where γpn is the solution to the equation ˜ φint+γpn σφ ˜ φint−γpn σφ π(φ|x)dφ=p, or, equivalently since τ=ns 2λ=ns 2eφ, ns2e˜ φint±γpn σφ=(n−1)2 nexp ±γpn 2 n−1, π(τ|x)=χ2(τ|n−1), γpn is the unique solution to the equation Fn−1(n−1)2 ne+γpn √2 n−1−Fn−1(n−1)2 ne−γpn √2 n−1=p,(50) where Fνis the cumulative distribution function of a χ2 νdistribution. A numerical solution to (50) is immediately found with standard statistical software. However, a simple analytical approximation may be derived using the fact that the reference posterior distribution of φ=log λbecomes approximately normal (at a faster rate that any other simple function of σ) as the sample size nincreases. Using this approximation and (49), the p-credible intrinsic region for φis approximated by the interval Cint p(x,Φ)≈log n−1 n21 s2±qp2 n−1(51) where qpis the (p+1)/2 quantile of a standard normal distribution. By invariance, the p-credible intrinsic region for the variance σ2=e−φwill be approximated by s2n n−12e−γnp √2 n−1,s2n n−12e+γnp √2 n−1(52) where γnp is the solution to (50) which, as nincreases, converges to qp,the(p+1)/2 quantile of an standard normal distribution. Summarizing, we have thus established Theorem 3 (Intrinsic estimation of the normal variance) Let xbe a random sample of size n from N(x|μ, σ), with variance s2. 26 Objective Bayesian point and region estimation in location-scale models (i) The intrinsic point estimator ˜σ2 int(x)of σ2is ˜σ2 int(x)=arg min ˜σ>0d(˜σ2|x), d(˜σ2|x)=n 2∞ 0 δ˜σ2τ ns 2χ2(τ|n−1) dτ, δ{θ}=min{g(θ),g(1/θ)},g(θ)=(θ−1) −log θ, ˜σ2 int(x)≈n n−12s2. The intrinsic point estimator ˜σint(x)is the Bayes estimator with respect to the intrinsic discrepancy loss. Besides, it has smaller frequentist risk with respect to this loss than any other affine equivariant estimator. (ii) The unique p-credible intrinsic region Cint p(x)for σ2is the interval Cint p(x)={a(α, x),b(α, p,x)}=ns 2 χ2 n−1(1 −α),ns 2 χ2 n−1(1 −p−α), where χ2 n−1(q)is the q quantile of a χ2 n−1distribution, and αis the solution to the equation d{a(α, x)|x)}=d{b(α, p,x)|x)}. For all sample sizes, the frequentist coverage of Cint p(x)is exactly p. Moreover, Cint p(x)≈s2n n−12e−γnp √2 n−1,e+γnp √2 n−1 where γnp is the solution to the equation Fn−1(n−1)2 ne+γnp √2 n−1−Fn−1(n−1)2 ne−γnp √2 n−1=p. and Fνis the cumulative distribution function of a χ2 νdistribution. As n increases, γnp converges to the (p+1)/2normal quantile. (iii) The expected intrinsic loss associated to the use of ˜σ2is d(˜σ2|s2,n)≈n 22 n−1+log 1 σ2 int −log 1 ˜σ22, with ˜σ2 int(x)≈n2s2/(n−1)2. For the numerical illustration considered in Example 4 (where the sample size was only n=12), the approximation (44) to the intrinsic estimate of σ2yields ˜σ2 int ≈5.104. The approximation (52) to the intrinsic 0.90-credible region yields (2.480,11.507) using the exact solution γnp =1.693 to equation (50), and (2.531,11.272) using the J. M. Bernardo 27 corresponding normal approximation γnp ≈q0.90 =1.645. These approximations may be compared eith the exact values ˜σ2 int =5.090 and Cint 0.90 =(2.481,11.717) numerically found in Example 4. 3.4 Intrinsic estimation of the Cauchy scale parameter We finally consider an example where no analytical expressions are possible. With the use of increasingly complex statistical models, this is fast becoming the norm, rather than the exception, in statistical practice. Let x={x1,...,xn}be a random sample from a Cauchy distribution Ca(x|0,σ), centered at zero with unknown scale parameter σ, so that the likelihood function is p(x|σ)= n  i=1 Ca(xj|0,σ)∝σ−n n  i=11+x2 j σ2−1. The Cauchy distribution does not belong to the exponential family and, therefore, there is no sufficient statistic of finite dimension. There is no analytical expression for the MLE ˆσof the unknown parameter. Fisher information function is n/(2σ2) and, therefore, the posterior distribution of σwill be asymptotically normal, N(σ|ˆσ, √2ˆσ/ √n). Since this is a scale model, the reference prior is π(σ)=σ−1and, using Bayes theorem, the reference posterior is π(σ|x)=σ−1p(x|σ) ∞ 0σ−1p(x|σ)dσ ,(53) which may easily be numerically computed. It may be verified that, provided the data x contain at least two different observations, π(σ|x) has a gamma-like shape with a unique mode. Figure 4 represents the reference posteriors of σwhich correspond to a set of 25 random samples of size n=12, which were all generated from a Cauchy distribution Ca(x|0,2). This may be seen as a graphical representation of the sampling distribution of πσ(·|x),the reference posterior of σ,givenσ=2andn=12. Notice that, although all these posteriors contain indeed the true value σ=2 from which the samples have been simulated (marked in the figure with a solid dot), the variability is very large. The logarithmic divergence κ{σ2|σ1}of Ca(x|0,σ 2) from Ca(x|0,σ 1)is ∞ −∞ Ca(x|0,σ 1)logCa(x|0,σ 1) Ca(x|0,σ 2)dx =log 1 4σ1σ2 +2log(σ1+σ2). Miguel A. G´ omez Villegas Departamento de Estad´ ıstica e Investigaci´ on Operativa Universidad Complutense de Madrid Let me begin by congratulating Professor Bernardo for his excellent job in objective Bayesian analysis. This paper, and the closely related Bernardo (2005), present a unified theory of estimation by point and credible regions based on information ideas he has used previously to define reference priors. The idea originates from the study of both problems as decision problems, where the loss function is the “intrinsic discrepancy” inspired in the Kullblack-Leibler divergence, and defined as the minimum of kx{ θ, λ|θ, λ} and kx{θ, λ| θ, λ}where kx{ θ, λ|θ, λ}=χ(θ,λ) π(x|θ, λ)lnπ(x|θ, λ) π(x| θ, λ)dx An intrinsic point estimator is then defined as the Bayes estimator which corresponds to the intrinsic loss and the appropriate reference prior. A p-credible intrinsic region estimator is defined as the lowest posterior loss p-credible with respect to the intrinsic loss and the appropriate reference prior. Afirst question is: do we need to employ Cint p π(θ|x)dθ≥p with the inequality instead of equality to allow the discrete case? Second, it would be useful to have a better understanding of the proposed approach to applying these ideas to the exponential distribution family instead of location-scale models; this is a family of distributions greater than the other. Professor Bernardo claims that in one-dimensional problems, one may define probability centred credible intervals, and these are invariant under reparametrization. Will it not be necessary to suppose that the transformation is monotonic? Third, on a more philosophical basis, I think that invariance is a compelling argument for point estimations and for credible regions. Indeed both point estimations and credible regions are two answers to the same question: how we can eliminate the uncertainty 36 about θ. Bernardo’s approach permits one to obtain invariance under reparametrization in both problems. Fourth, the chosen examples show the coherence between frequentist inference and Bayesian inference. When intrinsic credible regions that require minimal subjective inputs are employed, exact frequentist confidence regions are obtained, at least in the normal mean and variance. This fact is similar to the one obtained by this discussant in G´ omez-Villegas and Gonz´ alez-P´ erez (2005) and references therein. I wonder if Professor Bernardo has any idea about the essential reasons behind the matching properties between intrinsic credible regions and confidence regions in these cases? Fifth, adopting this approach to credible set construction, I see problems in computations, the posterior intrinsic loss integrated over a large dimensional space. From the point of view of applications, a simple asymptotic approximation to normality should be necessary. In closing, I would like to thank the editor of the journal for giving me the opportunity to discuss this paper. References Bernardo, J. M. (2005). Intrinsic credible regions: an objective Bayesian approach to interval estimation. Test, 14, 317-384. G´ omez-Villegas, M. A. and Gonz´ alez-P´ erez, B. (2005). Bayesian analysis of contingency tables. Communications in Statistics-Theory and Methods, 34, 1743-1754. Dennis V. Lindley [email protected] Two concepts are basic to the ideas of this excellent paper: objectivity and the concept of estimation as a decision problem. In the author’s skilful hands, these lead to reference priors and intrinsic loss functions, and hence, by minimizing expected loss, to estimates which are often superior to the conventional ones. It can be said with some confidence that we have here a solution to the problem Harold Jeffreys first posed around 1939 of providing an objective, coherent method for scientific inference. The development employs several subtle ideas, and considerable mathematical complexity, but one feature that struck me is that the final results are usually fairly simple and look right. An example of this is provided by the loss functions in Figure 3, which have the reasonable convexity property around the true value but, unlike quadratic loss, exhibit sensible concavity at more discrepant values. I would have preferred the loss to have been bounded but, with normal distributions and their thin tails, this scarcely matters. To be bounded may be more important with the fat tails of the Cauchy in Figure 6, in order to avoid paradoxes of the St. Petersburg type. A related point is that although the mathematics can be formidable, at least in the view of some applied statisticians, once it has been done the practitioner can easily use the results in the confidence that the machinery used to produce them is sound. It is comparable to driving a car, without knowing how it was made, but having confidence in the manufacturer. Granted the basic concepts, this is an important paper, but was Jeffreys right to search for objectivity, and was Fisher wrong in dismissing decision concepts from inference? I think Jeffreys was wrong and Fisher was right. At the risk of repeating what I have said before, it seems to me that inference and decision-making are distinct and both are subjective. In other words, the two basic concepts, that provide the foundations of this paper, are suspect. Consider first the fixed likelihood upon which all the arguments in the paper rest. Is it really objective? There are a few cases where substantial evidence for normality exists, but often the normal, or another member of the exponential family, is used merely for mathematical simplicity. With the increased computing power available today, statisticians are less constrained and can use other distributions that appear more realistic, thereby introducing subjectivity. There are some popular data sets that have been repeatedly analysed using different likelihoods. Where is the objectivity there? It is interesting that Bernardo uses one symbol, p, for probabilities of data but another, 38 π, for probabilities of parameters. In reality both pand πreflect beliefs about data and parameters respectively, obey the same rules and do not deserve separate treatments. Inference concerns parameters. (It is more practical to make inference about future data but I do not explore that trail here.) What are these parameters, the θand λof the paper? If our statistical analyses are to be of use in data analysis, θat least ought to relate to something in the real world. Bernardo has only one sentence about this, referring to θas the age of the earth. Putting aside intelligent designers, reputable scientists differ in their views of the age. In other words, their ideas are subjective, so that before relevant data about the age are considered, their different views need to be included. Another relevant fact is that information about the age of the earth does not come from data with normal, or any other objective, likelihood. More conspicuous examples of subjectivity are apparent witch clinical trials, where the different views of drug companies and official bodies are consulted before the trial. This became clear recently when a trial went horribly wrong and experts claimed the probabilities used were, in their opinion, unsound. So often today, θis regarded as nothing more than a symbol, whereas, to be of value, it has to refer to reality and hence influenced by opinions about that reality. These opinions should be incorporated into the analysis, not ignored and replaced by a reference prior, especially when this is improper. All of us have, at some time, expressed an opinion about something without having any intention of basing any action upon that opinion. In statistics, this opinion-forming is inference and means we infer the value of the real θ. In the Bayesian paradigm this is done by means of your probability distribution of θ, given the data and the original information about θ. Whilst it is true that any inference has to be capable of being used as a basis for action, for otherwise what use is it, it is not true that inference has to have immediate actions in mind. In particular, inference does not require a loss function, and certainly not a loss function that ignores reality. In Bayesian terms, there is only one inference, the posterior distribution and, although it may be advantageous to summarize its main features, such approximations scarcely need elaborate techniques, except in the case of many parameters. Inference from data consists in modelling that data in the form of a likelihood depending on parameters, supplying your opinion of the parameters prior to the data, and combining likelihood and prior by Bayes theorem. Finally the nuisance aspect of the parameters is removed by integration. When several people are involved there may be disagreements over likelihood or prior. These may be removed by discussion but, if this fails, the calculations may be repeated under different subjective opinions and the posteriors compared. That science is objective is a myth. Apparent objectivity in science only arises when the data are extensive. This paper explores a field that, in my view, is not in the broad stream of statistics. This is not to deny it great merit, for we now know what that field contains, material of real merit from which all can learn. Mark J. Schervish Carnegie Mellon University, USA I admire Professor Bernardo for his steadfastness and resolution in staying the course of research into reference priors and other so-called objective Bayesian methods. Despite repeated attacks dating back to the discussion of Bernardo (1979) he has continually risen to the challenge of making these methods palatable to practitioners and theoreticians alike. I will not here rehearse all of the criticisms or the support for his work in this area. I refer the interested reader to the various discussions of the papers listed in the reference list to Professor Bernardo’s paper. I will mention just a few problems that I have with the methods as well as what I like about them. To begin with a positive note, I like the idea of having a transformation-equivariant estimation procedure for non-decision-theoretic inference. When one is faced with a decision problem in which a specific loss function is relevant, then one does not care whether one’s inference satisfies an ad hoc criterion such as transformation equivariance. On the other hand, when one merely wishes to report an estimate of some quantity, especially the parameter of a statistical model which most likely is a figment of one’s imagination (model) anyway, then it becomes difficult to explain why the estimate of an equivalent parameter is not the equivalent estimate. Indeed, I believe that the intrinsic discrepancy loss satisfies a slightly stronger invariance than is stated in (10). I believe that one could apply a one-to-one reparameterization of the form φ=φ(θ)and ψ=ψ(λ, θ) and still achieve (10). Of course, a completely general reparameterization would change the meaning of the parameter of interest, and yet the desire for an equivariant estimate would remain. One of the serious concerns with reference priors is their violation of the likelihood principle. The reference priors are different for binomial sampling and negative binomial sampling so that even if the observed data could have come from either sampling scheme, the posterior would depend on the sampling plan. If one were to observe a binomial sample and use the reference prior, and later observe a negative binomial sample, one would get a different inference than if one were to observe the same two samples in the other order. As mentioned earlier, various discussants have described other concerns with the methods advocated in the manuscript, and I will let the reader find them in their original forms. I will add only one other concern that I have, and that is with the use of the description of these methods as “objective”. I suppose that, so long as one agrees with all of the reasons put forth for why such methods should be 40 used, then one will use the methods and they become objective in that sense. But any set of methods could be called objective on those grounds. One of the main strengths of Bayesian methodology is that it forces users to be explicit about the assumptions that they are making. People who think that they are using objective methods are simply borrowing a collection of subjective assumptions and ignoring the fact that choices were made by someone else arriving at those assumptions. When you lay your assumptions out for all to see, you are in a position to evaluate the sensitivity of your inferences to the assumptions. If you hide behind a cloak of objectivity, you may produce the same answer that others produce, but you have lost the ability to see what is the effect of the subjective choices that were made. Rejoinder I am extremely grateful to the three discussants by their thoughtful comments. I will answer them individually. G´omez-Villegas. If the parameter of interest θ θ θis discrete, then we would certainly need to work with regions Csuch that Cπ(θ|θ θ θ)dθ θ θ≥psince, in that special case, not all credible probabilities pwould be attainable. However, point and region estimation are usually done with continuous parameter spaces, and this is indeed the case in the location and scale models considered in this paper. In that situation, the equality may always be obtained. The ideas discussed in the paper may certainly be applied to models in the (generalized) exponential family and it is likely that this would lead to some rather general results. I did not have time and space to do this here, but it is certainly a research line well worth exploring. As Professor G´ omez-Villegas points out, the invariance arguments invoked only refer to monotonic, one-to-one transformations of the parameter. Even though not always explicitly stated, we were indeed always assuming this to be the case. I believe that the exact numerical coincidence between objective credible regions and frequentist confidence interval is the exception, not the rule; when it happens, it is the consequence of the existence of pivotal quantities, so that the reference distribution of the pivot (considered as a function of the parameter) is precisely the same as its sampling distribution (considered as a function of the data). In particular, this coincidence cannot exist if data are discrete, as in the case of binomial or Poisson data. Beyond the particular situations where pivots exist, one may only expect an asymptotic approximation: objective credible regions are typically approximate confidence intervals, the approximation improving with the sample size. Routine application of the methods described in this paper will certainly require either available software producing the exact results (not difficult to write in the standard examples which constitute the vast majority of applications) and/or appropriate analytical approximations. The latter may easily be obtained, as in the examples contained in the paper, by using the normal approximation with the parametrization induced by the appropriate variance-stabilizing transformation, and then making use of the invariance properties of the procedures. Lindley. I am really proud that Professor Lindley may believe that the procedures described provide an objective coherent method for scientific inference in the sense demanded by Jeffreys, and I am very grateful for that comment. 42 It would certainly be better from a foundations viewpoint if the expected loss were bounded, but information measures with continuous parameters are are not bounded (one needs infinite amount of information to know precisely a real number) and yet have all kind of attractive properties. To repeat in print the basics of an argument that Professor Lindley and I have often had in private conversations, (i) I believe, with Jeffreys, that Fisher was wrong in dismissing decision concepts in inference. If, by some reason, you must choose an estimate, then (whether you like it or not) you have a well posed decision problem where the action space is the set of parameter values; then foundations dictate that (to act rationally) you must use a loss function. For instance, in one continuous parameter problems, the median may well be an estimate with good robustness properties, but the fact remains that this would be a good estimate if (and only if ) your loss function is well approximated by a linear, symmetric loss function. (ii) I applaud the use of subjective priors when the problem is simple and small enough for the required probability assessments to be feasible (which is not frequent). But, even in this case, there is no reason while other people should necessarily accept a subjective prior which goes beyond clearly stated assumptions and verifiable (possibly historical) data. There is a clear need for some commonly accepted minimum set of conclusions to be solely derived from assumptions and data, and this is precisely what reference posteriors provide. As their name indicate, they are proposed as a reference, to be compared with subjective posteriors when these are available. This is part of a necessary exercise in sensitivity analysis, by making explicit which parts of the conclusions depend on a particular subjective prior, and which parts are implied by the model assumed and the data obtained. As Professor Lindley points out, although inferential statements are typically used as a basis for action, there are many situations were inferences are to be drawn without any specific action in mind. This is precisely why we suggest the use of an the informationbased loss function. If a particular action is in mind, one should certainly use a context dependent loss function which appropriately describes the decision problem analyzed. It no particular decision problem is in mind, one is bound to use some conventional loss function. We have argued that conventional loss functions (such as the ubiquitous quadratic loss) are often unsatisfactory. Instead, for “pure inference” problems one should try to minimize the information loss due to the use of an estimate of the unknown parameter value; and this, I believe, is appropriately captured by the intrinsic discrepancy loss. Schervish. I am very glad to read that Professor Schervish appreciates the importance of invariant procedures. In teaching, I often start my lectures by stating that any inferential 43 procedure which is not invariant under monotonic transformations of the parameter is suspect, and go on to provide a set of examples of those as “counterexamples” to common statistical procedures. I agree with Professor Schervish on the importance of the likelihood principle, but I believe that the principle is actually compatible with a sensible use of reference distributions. Indeed, a reference posterior encapsulates, by definition, the (minimal) inferential statements you could proclaim about the parameter of a model if your prior was that maximizing the information that data generated from that particular model could possibly provide. If you change the model (even if the new model induces a proportional likelihood function), you change the reference prior. Thus, different reference posteriors corresponding to different sampling schemes with Bernoulli observations provide a collection of conditional answers (one for each sampling scheme one is willing to consider), which may all be part of the sensitivity analysis to changes in the prior mentioned above. Objectivity is indeed an emotionally charged word, and it should be explicitly qualified whenever it is used. No statistical analysis is seriously objective, if only because the choice of both the experiment design and the model used have typically very strong subjective inputs. However, the frequentist paradigm is sold as “objective” just because its conclusions are only conditional on the model assumed and the data obtained, and this objectivity illusion has historically helped frequentist to keep a large share of the statistics market. I claim for the procedures described in this paper the right to use “objective” in precisely the same sense: these are procedures which are only conditional on the assumed model and the observed data. The use of the word “objective” in this precise, limited sense may benefit, I believe, the propagation of the Bayesian paradigm. For a recent discussion of this and related issues see Berger (2006) and ensuing discussion. I fully agree with Professor Schervish on the paramount importance of clearly presenting the assumptions needed for an inferential statement. In the case of reference posteriors this should typically read as a conditional statement of the form: “If available data xhad been generated by model M≡{px(·|ω ω ω),ω ω ω∈Ω}and prior information about θ θ θ(ω ω ω) were minimal with respect to the information about θ θ θ(ω ω ω) that repeated sampling from Mcould possibly provide then, the marginal reference posterior π(θ θ θ|x) encapsulates what could be said about the value of θ θ θ, solely on the basis of that information”. References Berger, J. O. (2006). The case for objective Bayesian analysis. Bayesian Analysis, 1, 385-402 and 457-464 (with discussion).