Monte Carlo Conditioning on a Sufficient Statistic

Bo H. Lindqvist Department of Mathematical Sciences Norwegian University of Science and Technology N-7491 Trondheim, Norway

Gunnar Taraldsen Department of Physiology and Biomedical Engineering Norwegian University of Science and Technology N-7491 Trondheim, Norway

February 21, 2000

1 Introduction

Let X and T be a random vectors with distribution indexed by a parameter θ ∈ Ω_{Θ}. Suppose further that T is sufficient for θ compared to X. By definition this means that the conditional distribution of X given T does not depend on θ. This property is basic in statistical inference, with applications for example in parameter elimination, construction of optimal estimators and in testing (Lehmann 1997, Lehmann and Casella 1998). In practical cases, however, such conditional distributions may be extremely difficult to derive exactly, so there is a strong need for efficient approximations or Monte Carlo algorithms.

Of course the most straightforward way of computing the conditional distribution of X given T = t would be to fix a parameter value θ_{0} and then compute or simulate the conditional distribution under θ_{0}. In this way the problem is a traditional one. However, such an approach does not make full use of the fact that any value for θ_{0} will give the same conditional distribution. Thus it is conceivable that the possibility of varying the value of θ in an appropriate manner may be more efficient. The key feature of this paper is to present an approach in this direction, where the clue is to equip the parameter space Ω_{Θ} with some suitable σ-finite measure. This idea is from Engen and Lillegård [?], whose main result is essentially our Proposition 1.

The basic assumption is that there is a random vector U taking values in a space Ω_{U} with a known distribution, and functions χ, τ such that (χ(U, θ), τ(U, θ)) and (X, T) have equal distributions for each fixed θ. The idea is then that U can be simulated on a computer and the function (τ, χ) makes it possible to simulate observations of (T, X) for given values of θ.

Now suppose it is given that T = t, and that a sample x_{t} from the conditional distribution of X given T = t is wanted. Since this distribution does not depend on θ, it is reasonable to believe that it can be described in some simple way in terms of the distribution of U, and thus enabling Monte Carlo simulation based on U.

As already mentioned, the basic feature of our general approach is to introduce a σ-finite measure P_{Θ} on Ω_{Θ}. Assuming for a moment that this is a probability measure of a random variable Θ, we assume in addition that U, Θ are stochastically independent. The key result is then Theorem 1 which in the setting of probability measures states that the conditional distribution of X given T = t equals the conditional distribution of χ(U, Θ) given τ(U, Θ) = t. This means that samples x_{t} from the condtional distribution of X given T = t can always be found by sampling (u, θ) from the conditional distribution of (U, Θ) given τ(U, Θ) = t and then computing x_{t} = χ(u, θ). The latter conditional distribution will therefore be of main importance in the paper.

2 Assumptions and Notation

As noted in the introduction, we shall find it convenient to supply the parameter space Ω_{Θ} by a σ-finite measure M_{Θ}. This leads to the need for modifying the definition of conditional expectation to allow for σ-finite distributions. Thus let μ be a σ-finite measure on some measurable space and let Y and Z be measurable functions defined on this space. We shall define what to mean by the conditional expectation μ(Y |Z = z) when the distribution μ_{Z} of Z is σ-finite. Since we will include the case when Y is constant, we do not require Y to be integrable, but require it to be bounded. The conditional expectation μ(Y |Z = z) is defined by requiring

μ(Y ϕ(Z)) = ∫ μ(Y |Z = z)ϕ(z)μ_{Z}(dz)

to hold for all μ_{Z}-integrable ϕ. Here we use the common notation μX to mean the integral of X with respect to μ. Let [Z ∈ A] denote the indicator function of the set (Z ∈ A) = {w|Z(w) ∈ A}. The assumption of σ-finiteness of μ_{Z} ensures that A ↦ μ(Y [Z ∈ A]) is a σ-finite measure. This measure is absolutely continuous with respect to μ_{Z}, and the Radon-Nikodym Theorem [?, p.121, p.123] gives the conditional expectation as its density: μ(Y [Z ∈ dz]) = μ(Y |Z = z)μ_{Z}(dz). This conditional expectation is normalized and has all the other properties of a conventional expectation. If μ is a probability measure it corresponds to the ordinary definition of conditional expectation. In most cases of interest, e.g. when the involved spaces are Borel spaces [?, p.618], the conditional expectation is given by a conditional distribution. We will in the sequel use equality sign between conditional expectations to mean that they are both versions of the conditional expectations and thus equal only a.e.

The general setup of the paper is as follows. Let Ω be a measurable space equipped with a family of probability measures P^{θ} indexed by a parameter θ in Ω_{Θ}. Let further Ω_{T} and Ω_{X} be measurable spaces and let there be given statistics T and X which are measurable functions T : Ω → Ω_{T} and X : Ω → Ω_{X}.

We shall assume that T is sufficient for θ compared to X, which means that the conditional expectation E^{θ}(ϕ(X)|T = t) can be specified independent of θ. We do in particular not demand T to be a function of X when we define sufficiency of T with respect to X. This is in accordance with the general sufficiency concept of Blackwell (1953), which includes the usual concept where it is assumed that T = T (X).

Let M_{U} be a probability measure on a measurable space Ω_{U}, which we shall call the Monte Carlo space. The basic assumption is that there are measurable functions τ and χ such that the distribution of (τ(U, θ), χ(U, θ)) equals the distribution of (T, X) under P^{θ} for each fixed θ. Moreover, let there be given a σ-finite measure M_{Θ} on the parameter space Ω_{Θ}, such that τ(U, Θ) has a σ-finite distribution M_{τ} induced by the product measure M ≡ M_{U} ⊗ M_{Θ}. Functions of U will be called Monte Carlo variables.

3 General results

The following result, which is essentially the main result of Engen and Lillegård (1997), will be generalized in the subsequent theorem. We state it here since it is in a particularly useful form for applications in discrete models, and since it gives much of the flavour of the general method. The proof is elementary and can be compared with the proof given in [?].

Proposition 1 Let ϕ : Ω_{X} → R be bounded and measurable, and define the Monte Carlo variables Z_{t} and Y_{t} by

Z_{t} = ∫ ϕ(χ(U, θ)) · [τ(U, θ) = t] M_{Θ}(dθ), Y_{t} = ∫ [τ(U, θ) = t] M_{Θ}(dθ).

If M_{U}Y_{t} > 0, then

E^{θ}(ϕ(X)|T = t) = \frac{M_{U}Z_{t}}{M_{U}Y_{t}} .

Proof: First we observe that M_{U}Y_{t} is finite since M_{τ} is σ-finite. The composed function ϕ ∘ χ is not more general than χ alone, so in the proof we may assume that ϕ(x) = x. By considering only θ’s such that ¶^{θ}(T = t) ≠ 0 it follows that

E^{θ}(X|T = t) = \frac{E^{θ}(X[T = t])}{¶^{θ}(T = t)} = \frac{∫ E^{θ}(X[T = t]) M_{Θ}(dθ)}{∫ ¶^{θ}(T = t) M_{Θ}(dθ)} = \frac{M(χ(U, Θ)[τ(U, Θ) = t])}{M(τ(U, Θ) = t)}.

The above integration is valid since the sufficiency gives that the fraction is independent of θ. The conclusion follows from the Fubini theorem which allows us to change the order of integration since M_{U} and M_{Θ} are σ-finite.

The above proof shows in particular that the conditional distribution of X given T = t equals the conditional distribution of χ(U, Θ) given τ(U, Θ) = t. The next theorem shows that this fact holds in general, and is in fact basic to our approach.

Theorem 1 If ϕ : Ω_{X} → R is bounded and measurable, then

E^{θ}(ϕ(X)|T = t) = M(ϕ(χ(U, Θ))|τ(U, Θ) = t).  (1)

Proof: As in the proof of Proposition 1 we may ignore the function ϕ. Let ψ_{1}(t) = E^{θ}(X|T = t) and ψ_{2}(t) = M(χ(U, Θ)|τ(U, Θ) = t). We will prove M(ψ_{1}(τ(U, Θ)) ≠ ψ_{2}(τ(U, Θ))) = 0. Let ξ be a M_{τ} -integrable function. The claim follows from

Mξ(τ(U, Θ))ψ_{2}(τ(U, Θ)) = Mξ(τ(U, Θ))χ(U, Θ) = M[M_{U}ξ(τ(U, Θ))χ(U, Θ)]
= M[M_{U}ξ(τ(U, Θ))ψ_{1}(τ(U, Θ))] = Mξ(τ(U, Θ))ψ_{1}(τ(U, Θ)).

The next theorem provides a formula for computing the conditional expectation on the right hand side of (1). A specialization to the Euclidean case is given in a subsequent corollary.

As it turns out, our approach depends heavily on solving the equation τ(u, θ) = t with respect to θ. Introduce therefore the notation Γ(u, t) = {θ : τ(u, θ) = t}. As in Proposition 1 it is seen that the result can be expressed as a quotient of expected values of Monte Carlo variables.

Theorem 2 Suppose the conditional distribution of Θ given τ(u, Θ) = t for each fixed u ∈ Ω_{U} is given as a regular conditional probability distribution ρ(dθ|u, t), and which is such that integrals with respect to ρ(·|u, t) are jointly measurable in u, t. Let μ be a σ-finite measure on Ω_{T} such that τ(u, Θ) for fixed u has density Y_{t}(u) w.r.t. μ. Assume moreover that Y_{t}(u) is jointly measurable in u, t. Then

E^{θ}(ϕ(X)|T = t) = \frac{∫ ∫ ϕ(χ(u, θ))ρ(dθ|u, t)Y_{t}(u)M_{U}(du)}{∫ Y_{t}(u)M_{U}(du)} ≡ \frac{M_{U}Z_{t}}{M_{U}Y_{t}}

where Z_{t} = Y_{t} ∫ ϕ(χ(U, θ))ρ(dθ|U, t)

Proof: Let ν_{u}(dθ, dt) be the measure on Γ(u, ·) ≡ {(θ, t) : τ(u, θ) = t} defined by ∫ ∫ h(θ, t)ν_{u}(dθ, dt) = ∫ ∫ h(θ, t)ρ(dθ|u, t)μ(dt). Then the joint distribution of (Θ, τ(u, Θ)) has the density Y_{t}(u) with respect to the measure ν_{u}. Now for any measurable h : Ω_{U} × Ω_{Θ} × Ω_{T} → R such that (u, θ) ↦ h(u, θ, τ(u, θ)) is M-integrable, we have by Fubini’s theorem,

M(h(U, Θ, τ(U, Θ))
= ∫ ∫ h(u, θ, τ(u, θ))M_{Θ}(dθ)M_{U}(du)
= ∫ ∫ h(u, θ, t)Y_{t}(u)dν_{u}M_{U}(du)
= ∫ ∫ ∫ h(u, θ, t)Y_{t}(u)ρ(dθ|u, t)M_{U}(du)μ(dt)  (2)

In particular, if h(u, θ, t) is taken as a function of t alone, it follows that the marginal distribution of τ(U, Θ) w.r.t. μ is

M_{τ} (dt) = ∫ Y_{t}(u)M_{U}(du)μ(dt)  (3)

To prove the theorem we need by the definition of conditional expectation to show that for any M_{τ} -integrable function ψ we have

M(ϕ(χ(U, Θ))ψ(τ(U, Θ))) = ∫ \frac{∫ ∫ ϕ(χ(u, θ))Y_{t}(u)ρ(dθ|u, t)M_{U}(du)}{∫ Y_{t}(u)M_{U}(du)}ψ(t)M_{τ} (dt)

By (3) the right hand side equals ∫ ∫ ∫ ϕ(χ(u, θ))Y_{t}(u)ρ(dθ|u, t)ψ(t)M_{U}(du)μ(dt) and the result therefore follows directly from (2).

The following corollary specializes to the case where Ω_{Θ} and Ω_{T} are subsets of the same Euclidean space R^{k}. Let det ∂_{θ}τ(u, θ) be the determinant of the matrix of partial derivatives of τ(u, θ) for fixed u.

Corollary 1 Let Ω_{Θ} and Ω_{T} be open subsets of R^{k} with positive Lebesgue-measure, and assume M_{Θ}(θ) = f(θ)dθ where dθ means Lebesgue measure on Ω_{Θ}. Let furthermore θ ↦ τ(u, θ) be differentiable for each u, and assume that Γ(u, t) is a finite set for each fixed u, t.

Then

E^{θ}(ϕ(X)|T = t) = \frac{∫ ∑_{θ̃∈Γ(u,t)} ϕ(χ(u, θ̃))|\frac{f(θ)}{det ∂_{θ}τ(u,θ)}|_{θ=θ̃}M_{U}(du)}{∫ ∑_{θ̃∈Γ(u,t)} |\frac{f(θ)}{det ∂_{θ}τ(u,θ)}|_{θ=θ̃}M_{U}(du)}

Proof: The result follows from the preceding theorem by noting that the density of τ(u, Θ) with respect to Lebesgue measure is

∑_{θ̃∈Γ(u,t)} |\frac{f(θ)}{det ∂_{θ}τ(u, θ)}|_{θ=θ̃}

while the conditional distribution of Θ given τ(u, Θ) = t is concentrated on Γ(u, t) and given by

\frac{|\frac{f(θ)}{det ∂_{θ}τ(u,θ)}|_{θ=θ^{∗}}}{∑_{θ̃∈Γ(u,t)} |\frac{f(θ)}{det ∂_{θ}τ(u,θ)}|_{θ=θ̃}}; θ^{∗} ∈ Γ(u, t)

Remark 1: As remarked earlier, conditional expectations given T = t are to be regarded as measuarable functions of t, with equality between two expressions meaning equality a.e. In applications the measurable sets in Ω_{T} are typically Borel sets of some topology and often there exist versions of the conditional expectations which are continuous functions of t. In these cases equality a.e. will imply equality everywhere.

Remark 2: The essential feature of Theorem 2 is to evaluate the right hand side of (1) by double expectation by conditioning on U. Formally,

M(ϕ(χ(U, Θ))|τ(U, Θ) = t) = E_{U|τ} M(ϕ(χ(U, Θ))|U, τ(U, Θ) = t)  (4)

where the outer expectation is with respect to the conditional distribution of U given τ(U, Θ) = t. With notation as in Theorem 2, this conditional distribution is given as

\frac{Y_{t}(u)M_{U}(du)}{∫ Y_{t}(u)M_{U}(du)}

while the inner expectation in (4) is ∫ ϕ(χ(u, θ))ρ(dθ|u, t).

Of course we may well change the roles of Θ and U in (4), thus conditioning on (Θ, τ(U, Θ)) instead of (U, τ(U, Θ)). We would then recognize the inner expectation as simply the conditional expectation of ϕ(X) given T = t. Since this by sufficiency does not depend on θ, the outer expectation becomes trivial and we are back to the left hand side of (1).

This indicates that the essential feature of Theorem 2 is that computations of conditional expectations given τ(U, θ) = t for fixed θ are changed into the (often much simpler) problem of computing conditional expectations given τ(u, Θ) = t for fixed u. It may in fact be the case that the equation τ(u, θ) = t has a single solution for θ (see next section).

Remark 3: The following derivation for the case where both U and Θ are discretely distributed, illustrates the main idea of the theorem. Let for simplicity ϕ be the identity function.

E(X|T = t) = E(χ(U, Θ)|τ(U, Θ) = t)

∑_{u} ∑_{θ} χ(u, θ)P (U = u, Θ = θ|τ(U, Θ) = t)

= ∑_{u} ∑_{θ} χ(u, θ)\frac{P (U = u, Θ = θ, τ(U, Θ) = t)}{P (τ(U, Θ) = t)}

= \frac{∑_{u} ∑_{θ} χ(u, θ)P (Θ = θ, τ(U, Θ) = t|U = u)P (U = u)}{∑_{u} P (τ(U, Θ) = t|U = u)P (U = u)}

= \frac{∑_{u} ∑_{θ} χ(u, θ)P (Θ = θ, τ(u, Θ) = t)P (U = u)}{∑_{u} P (τ(u, Θ) = t)P (U = u)}

= \frac{∑_{u}\{∑_{θ} χ(u, θ)P (Θ = θ|τ(u, Θ) = t)\}P (τ(u, Θ) = t)P (U = u)}{∑_{u} P (τ(u, Θ) = t)P (U = u)}

where the last formula corresponds to the expression in Theorem 2. Alternative expressions can now be found which correspond to the discussion in the previous remark:

= ∑_{u}{∑_{θ} χ(u, θ)P (Θ = θ|τ(u, Θ) = t)}P (U = u|τ(U, Θ) = t)
= ∑_{u} E{χ(u, Θ)|τ(u, Θ) = t}P (U = u|τ(U, Θ) = t)

Remark 4: Proposition 1 follows directly from Theorem 2 by noting that ρ(dθ|u, t) is now the measure on Ω_{Θ} with density at θ given by

\frac{[τ(u, θ) = t]M_{Θ}(dθ)}{Y_{t}(u)}

4 Unique solution for θ

Assumption 1 (Unique θ.) Assume that there is for each fixed pair (u, t) a unique solution for θ of the equation τ(u, θ) = t.

The assumption implies that each Γ(u, t) has a single element, which we shall denote by θ̂(u, t). It follows that the conditional distribution ρ(dθ|u, t) is concentrated at θ̂(u, t). The Monte Carlo variable X_{t} = χ(U, θ̂(U, t)) will play an important role in the following.

As a consequence of the results of the previous section we get the following result valid under Assumption 1:

Theorem 3 Make Assumption 1. If ϕ : Ω_{X} → R is bounded and measurable, then

E^{θ}(ϕ(X)|T = t) = M(ϕ(X_{t})|τ(U, Θ) = t).  (5)

Let μ be a σ-finite measure on Ω_{T} such that τ(u, Θ) for fixed u has density Y_{t}(u) w.r.t. μ. Assume moreover that Y_{t}(u) is jointly measurable in u, t. Then

E^{θ}(ϕ(X)|T = t) = \frac{M_{U}ϕ(X_{t})Y_{t}}{M_{U}Y_{t}}  (6)

Suppose that Ω_{Θ} and Ω_{T} are open subsets of R^{k} with positive Lebesgue-measure, and that M_{Θ}(θ) = f(θ)dθ where dθ means Lebesgue measure on Ω_{Θ}. If θ ↦ τ(u, θ) is differentiable for each u, then

Y_{t}(u) = |\frac{f(θ)}{det ∂_{θ}τ(u, θ)}|_{θ=θ̂(u,t)}  (7)

Recall that the simple algorithm given in the introduction, if valid, would imply that X_{t} under Assumption 1 has the same distribution as X given T = t. As already noted in the introduction, this is not the case in general, however. If it is the case, then we shall say for short that “X_{t} has the correct distribution”.

By (5) in Theorem 2 it follows that X_{t} has the correct distribution if and only if X_{t} and τ(U, Θ) are independent (extending the concept of independence in an obvious way to the case of σ-finite distributions). The formula (6) in Theorem 2 implies, on the other hand, that X_{t} has the correct distribution whenever X_{t} is stochastically independent of Y_{t}. This happens in particular if we can find a measure M_{Θ} such that Y_{t}(u) is constant in u. Under the assumptions leading to (7) this corresponds to finding a function f(θ) such that Y_{t}(u) becomes constant in u. This is the case provided

| det ∂_{θ}τ(u, θ)|_{θ=θ̂(u,t)} = J(θ̂(u, t), t)  (8)

for some function J. Suppose namely that we are interested in the conditional distribution of X given T = t_{0} for some specific value t_{0}. Then we may choose

f(θ) = J(θ, t_{0})  (9)

which leads to a function Y_{t}(u) satisfying Y_{t_{0}}(u) = 1. Then under continuity conditions (see Remark 1 of the previous section) it follows from (6) that X_{t} has the correct distribution for all t.

Moreover, if we under condition (10) choose instead f(θ) = ψ(θ)(J(θ, t_{0})) which gives Y_{t_{0}}(u) = ψ(θ̂(u, t_{0})). Then (6) implies

M_{U}(ϕ(X_{t_{0}})ψ(θ̂(U, t_{0}))) = M_{U}(ϕ(X_{t_{0}}))M_{U}(ψ(θ̂(U, t_{0})))

for all bounded and measurable ϕ and ψ, which finally implies that X_{t} and θ̂(U, t) are independent for all t when (10) holds.

Note that the continuity conditions can be avoided under the condition

| det ∂_{θ}τ(u, θ)|_{θ=θ̂(u,t)} = a(t)J(θ̂(u, t))  (10)

in which case we can choose f(θ) = J(θ) and get Y_{t}(u) = a(t). Now (6) directly gives E^{θ}(ϕ(X)|T = t) = M_{U}ϕ(X_{t}).

We close the section by showing that X_{t} has the correct distribution when a certain pivotal property holds.

Assumption 2 (The pivotal case.) Assume that there exist measurable functions r : Ω_{U} → Ω_{V}, τ̃ : Ω_{V} × Ω_{Θ} → Ω_{T} with τ(u, θ) = τ̃(r(u), θ) such that the equation τ̃(v, θ) = t has a unique solution v = ṽ(θ, t) for all (θ, t).

In words this means that τ(u, θ) depends on u only through r(u), and that for a given θ this r(u) is uniquely determined by t. Moreover, if Assumption 2 holds, then ṽ(θ, T) is a pivotal quantity in the classical meaning. To see this, note the identity r(U) = ṽ(θ, τ(U, θ)), which has a distribution which is not a function of θ. The theorem thus follows since t = τ̃(v, θ) is the unique solution of v = ṽ(θ, t).

Theorem 4 Make Assumptions 1 and 2. Then X_{t} has the correct distribution. Furthermore X_{t} is independent of r(U) and hence of θ̂(U, t) and each τ(U, θ).

Proof: We will use the general rule E(X|Y = y) = E(X|h(Y) = h(y)), valid whenever h is one-one and measurable, repeatedly below. The idea is to replace the condition τ(U, θ) = t by the condition θ̂(U, t) = θ and thereby prove that the equality _{?} = in the introduction actually is an equality under the given assumptions. We have τ(U, θ) = τ̃(r(U), θ). The condition τ(U, θ) = t can therefore be replaced by r(U) = ṽ(θ, t). This can in turn be replaced by θ̃(r(U), t) = θ. Now θ̂(U, t) = θ̃(r(U), t) and this finishes the proof. The independence result of the theorem follows since r(U) and τ(U, θ) generate the same σ-field.

We finally note that under the assumptions of Theorem 4 and the additional assumption of the Jacobi case, condition (10) is always satisfied. Moreover, it is possible to give an explicit formula for the density f(θ) which makes the Y_{t}(u) of (7) identically equal to 1.

First it follows from Assumption 2 that r(u) = ṽ(θ̂(u, t), t). Hence τ(u, θ) = τ̃(ṽ(θ̂(u, t), t), θ), so det∂_{θ}τ(u, θ) = det∂_{θ}τ̃(ṽ(θ̂(u, t), t), θ). By substituting θ̂(u, t) for θ it is seen that (10) holds, and the desired f(θ) is given by

f(θ) = |det∂_{θ}τ̃(v, θ)|_{v=ṽ(θ,t)}  (11)

5 Monte Carlo Conditioning

The previous sections are devoted to deriving expressions for conditional expectations E^{θ}(ϕ(X)|T = t) in terms of expectations of functions of U. The idea is that these expectations can be evaluated by Monte Carlo simulation. Basically, we can estimate M_{U}h(U) by (1/N) ∑_{i=1}^{N} h(u_{i}) where u_{1}, …, u_{N} is a computer generated sequence corresponding to N independent realizations of U. Standard theory on Monte Carlo simulation contains methods for improving on this naive approach of estimating M_{U}h(U). We shall, however, not pursue this matter here.

In Section ... we investigated conditions under which X_{t}(u) from Assumption 1 directly gives samples from the conditional distribution X given T = t, thus not having to explicitly define a measure M_{Θ} on Ω_{Θ}. In general, however, we need to specify a measure M_{Θ} and the choice may influence the effectiveness of the algorithm. It should be stressed that M_{Θ} in a practical application with given t may well be chosen to depend on the value of t. Equation (10) gives an example of how t may enter the choice of distribution M_{Θ}.

Another instance where one wants to choose M_{Θ} according to t is discussed by Engen and Lillegård (1997). In discrete models the sets Γ(u, t) are typically intervals (see e.g. Example xxx). The endpoints may however be difficult to find explicitly. Engen and Lillegård (1997) suggest to let M_{Θ} be counting measure on a grid of points with center at θ^{∗}(t), the maximum likelihood estimate for θ when T = t.

A closer study of the formulae of Theorem 2? and 3? indicate that M_{Θ} should preferably be chosen to minimize the variation in u of the “weights” Y_{t}(u). As an example, consider the formula in Theorem 2old where the weights are given by

Y_{t}(u) = |\frac{f(θ)}{det ∂_{θ}τ(u, θ)}_{θ=θ̂(u,t)}  (12)

Here the density f(θ) should be chosen so that f(θ̂(u, t)) as a function of u is roughly proportional to |det∂_{θ}τ(u, θ)|_{θ=θ̂(u,t)}. In the pivotal case it follows by Theorem ... that we may choose f such that (12) is constant in u. For example, in the case where X is a sample from N(μ, σ) it will be shown in Example 6 that such an f is given by f(μ, σ) = 1/σ, which is the standard noninformative (improper) prior for this case. This motivates the use of vague prior distributions from Bayesian theory as natural choices for distributions M_{Θ}, for example Jeffreys’ priors (Jeffreys, 1946).

As a simple motivation, consider the (pivotal) case where θ is one-dimensional and τ(u, θ) = r(u)θ. This means that T/θ is a pivotal quantity. Assume that M_{U}r(U) = 1. In this case ∂_{θ}τ(u, θ) = r(u) and θ̂(u, t) = t/r(u). Hence Y_{t}(u) in (12) becomes constant in u by choosing f(θ) = 1/θ. Assuming that T is the maximum likelihood estimator of θ, then it can be seen that 1/θ is in fact the Jeffreys’ prior here. In connection with the use of Jeffreys’ priors, it is also interesting to note that the terms f(θ)/|det∂_{θ}τ(u, θ)| in Theorems 2 and 3? remain unchanged under one-to-one reparametrizations of θ if f is chosen by Jeffreys’ rule.

For a chosen measure M_{Θ}, we are able to compute conditional expectations E^{θ}(ϕ(X)|T = t) by simulation for given fundtions ϕ. In some cases it may be of interest, however, to simulate direct samples from the conditional distribtuion of X given T = t. This was the intention of Algorithm 1 in the introduction. In Section ... we showed how this can be done in certain special cases through X_{t}(u), but in general it is not straightforward to obtain such samples.

By Theorem 1 we need to sample (u, θ) from the conditional distribution of (U, Θ) given τ(U, Θ) = t and then compute x_{t} = χ(u, θ). Theorem 2 implies that the appropriate distribution from which to sample (u, θ) is

\frac{Y_{t}(u)ρ(dθ|u, t)M_{U}(du)}{∫ Y_{t}(u)ρ(dθ|u, t)M_{U}(du)}  (13)

In practice it may be easier to start by sampling (u, θ) as follows: (i) Draw u from the distribution M_{U}. (ii) Given U = u, draw θ from the conditional distribution of Θ given τ(u, Θ) = t. This means that (u, θ) is sampled from the distribution

Y_{t}(u)ρ(dθ|u, t)M_{U}(du)  (14)

It thus remains to modify samples coming from (14) to become samples from (13), for example rejection sampling. This requires, however, that we are able to bound the function Y_{t}(u) from above. This may be impossible or impractical in many cases. In such cases we are led to using for example Markov Chain Monte Carlo methods (Tierney 1994) or the SIR-algorithm (Rubin, 1988).

6 Examples

We now proceed by giving four examples with the aim of illustrating the theoretical results, as well as indicating their relevance in statistical inference problems.

[Normal distribution] This example illustrates in a well known setting the basic ideas of the paper. Moreover, it gives an example of a case where the algorithm of the introduction holds.

Let X = (X_{1}, …, X_{n}) be a vector of i.i.d. N(μ, σ) variables and let T = (X̄, S) be the sufficient statistic. Here X̄ = ∑_{i=1}^{n} X_{i}/n and S^{2} = ∑_{i=1}^{n}(X_{i} − X̄)^{2}/(n − 1). For simulation in this model we use a vector U = (U_{1}, …, U_{n}) of i.i.d. N(0, 1) variables and simulate X and T by, respectively,

χ(U, θ) = (μ + σU_{1}, …, μ + σU_{n})
τ(U, θ) = (μ + σŪ, σS_{U})

where θ = (μ, σ).

Now Assumption 1 holds with

θ̂(U, t) = (x̄ − \frac{Ū}{S_{U}} s, \frac{1}{S_{U}} s)

where t = (x̄, s), and Ū, S_{U} are similar to X̄, S. Also, Assumption 2 holds with r(U) = (Ū, S_{U}) and hence Theorem 4 implies that samples from the conditional distribution of X given T = t are given by

X_{t} = χ(U, θ̂(U, t)) = (x̄ + \frac{U_{1} − Ū}{S_{U}} s, …, x̄ + \frac{U_{n} − Ū}{S_{U}} s)

Thus X_{t} has a multivariate t-distribution.

The above formulas show that X_{t} and τ(U, θ) are stochastically independent (use for example Basu’s theorem, Lehmann and Casella, 1998, p. 42), likewise that X_{t} and θ̂(U, t) are independent. Thus from the discussion following Theorem ?? we can conclude that X_{t} has the correct distribution even without using Theorem 4.

Further, det∂_{t}θ̂(U, t) = 1/S_{U} so (??) holds with J(θ, t) = σ/s and we may hence choose f(μ, σ) = 1/σ for the density defining M_{Θ}. Accidentally, this is Jeffreys’ suggested joint prior for (μ, σ) when μ and σ are given independent priors (referanse).

[Nonhomogenous Poisson process] This example provides a counterexample to the algorithm of the introduction in the case where θ̂(U, t) is uniquely determined.

Consider a nonhomogenous Poisson process on the time interval [0, 1], with intensity function λ(t) = exp(δ + θt). Let N be the number of events, observed at times X_{1}, X_{2}, …, X_{N}. It is well known that given N = n, the X_{1}, X_{2}, …, X_{n} are distributed as the order statistic of n i.i.d. observations from the distribution function on [0, 1] given by

F (x, θ) = ∫_{0}^{x} λ(t)/ ∫_{0}^{1} λ(t)dt = \{ \begin{matrix}\frac{e^{θx}−1}{e^{θ}−1}, & θ ≠ 0 \\ x, & θ = 0 \end{matrix}  (15)

This distribution function depends analytically on the parameter θ ∈ R. The statistic T = ∑_{i=1}^{n} X_{i} is sufficient compared to X = (X_{1}, …, X_{n}).

Let U = (U_{1}, U_{2}, …, U_{n}) be a vector of n i.i.d. random variables uniformly distributed on [0, 1]. Inversion based on the distribution function (15) gives (X, T) ∼ (χ(U, θ), τ(U, θ)) where χ(u, θ) = (η(u_{1}, θ), …, η(u_{n}, θ)), τ(u, θ) = ∑_{i=1}^{n} η(u_{i}, θ), and

η(v, θ) = \{ \begin{matrix}\frac{log(1+(e^{θ}−1)v)}{θ} & θ ≠ 0 \\ v & θ = 0 \end{matrix}

The function η(v, θ) is analytic and strictly increasing in θ. This follows from the equation F (η(v, θ), θ) = v which determined η together with the observation that the distribution function F (x, θ) is strictly decreasing in θ. Consequently Assumption 1 is fulfilled and θ̂(u, t) is uniquely determined by τ(u, θ̂(u, t)) = t.

Consider for a moment the case when n = 2. It is then easy to check that the conditional distribution of X_{1} given T = t is uniform on (0, t) for 0 < t ≤ 1, and uniform on (t − 1, 1) for 1 < t < 2. Still we can prove that X_{t,1} = η(U_{1}, θ̂(U, t)) is not uniformly distributed for any 0 < t < 2. This can be demonstrated in general by numerical calculation, but we treat only the special case t = 1. Direct calculation gives

η(U_{1}, θ̂(U, t)) = \frac{log \frac{1−U_{1}}{U_{2}}}{log \frac{1−U_{1}}{U_{2}} + log \frac{1−U_{2}}{U_{1}}}

which is easily verified to be non-uniformly distributed and thus establishes the announced counter example.

Return now to the case of general n and suppose we want to compute E(ϕ(X)|T = t). By Theorem 2 we need to compute | det ∂_{θ}τ(u, θ)|_{θ=θ̂(u,t)}, and we get

∂_{θ}τ(u, θ)|_{θ=θ̂(u,t)} = \frac{e^{θ̂(u,t)}}{θ̂(u, t)} ∑_{i=1}^{n} \frac{u_{i}}{1 + (e^{θ̂(u,t)} − 1)u_{i}} − \frac{t}{θ̂(u, t)}

In principle we can use Theorem 2 with any choice of f(θ). As suggested in subsection 3.1 we should aim at flattening out the Y_{t}(u). A straightforward computation gives the square root of Fisherinformation to be (1/θ^{2} − exp(θ)/(exp(θ) − 1)^{2})^{1/2} which is the Jeffreys’ prior here and thus a possible choice for f.

We finally note that since θ = 0 corresponds to the X_{i} being uniformly distributed on [0, 1], we consider in fact the purely probabilistic problem of finding the conditional distribution of n i.i.d. uniform random variables given their sum.

[Logistic regression] In this example we have discrete observations and a continuous parameter. In this case we will typically have interval solutions for θ of the equation τ(u, θ) = t.

Let X = (X_{1}, …, X_{n}) be independent random variabels such that X_{i} is Bernoulli-distributed with

P (X_{i} = 1) = p_{i}(θ, δ) = \frac{e^{θ+δz_{i}}}{1 + e^{θ+δz_{i}}}

where z_{1}, …, z_{n} are given real covariates and −∞ < δ, θ < ∞ are parameters. Suppose δ is the parameter of interest. Since T = ∑_{i=1}^{n} X_{i} is sufficient for θ for fixed δ, we can eliminate the nuisance parameter θ by conditioning on T.

Thus, consider the value of δ as fixed at δ_{0} in the following. We shall see how the method of Proposition ?? can be used to find the conditional distribution of X given T = t. For example, this distribution is of interest when testing the null hypothesis H_{0} : δ = δ_{0}.

Let U = (U_{1}, …, U_{n}) be i.i.d. and uniform on [0, 1]. Then X under θ can be simulated as χ(U, θ) = (η_{1}(U_{1}, θ), η_{2}(U_{2}, θ), …, η_{n}(U_{n}, θ)) where

η_{i}(u_{i}, θ) = [u_{i} < p_{i}(θ, δ_{0})],  (16)

and T is therefore simulated by τ(U, θ) = ∑_{i=1}^{n} η_{i}(u_{i}, θ).

From(16) follows that

η_{i}(u_{i}, θ) = 1 ↔ θ > logit(u_{i}) − δ_{0}z_{i} ≡ ψ_{i}(u)

where logit(u) = log(u/(1 − u)). By ordering the ψ_{i}(u) as ψ_{(1)}(u) < ψ_{(2)}(u) < ⋯ ψ_{(n)}(u) it follows that τ(u, θ) = t if and only if ψ_{(t)}(u) < θ < ψ_{(t+1)}(u). Thus the set of θ satisfying τ(u, θ) is an interval and

Y_{t}(u) = M_{Θ}(ψ_{(t)}(u), ψ_{(t+1)}(u))  (17)

where M_{Θ} is a measure on the real numbers.

Since each η_{i}(u_{i}, θ) is constant (either 0 or 1) on the interval (ψ_{(t)}(u), ψ_{(t+1)}(u)) it follows that Assumption 1 holds and that we can use Proposition ?? for computing conditinal expectations. Also, we can use the SIR algorithm discussed in the previous section to obtain approximate samples from the conditional distribution.

Given the measure M_{Θ} the computation of weights Y_{t} is straightforward since we have the explicit solution (17). For choice of measure, it is by (17) mathematically convenient to use Lebesgue measure. Otherwise one may compute Jeffreys’ prior distribution or an approximation of it.

For the special case δ_{0} = 0 the conditional distribution of T is of course well known. We shall show that the above derived algorithm simplifies in this case, since we can show that X_{t}(u) and Y_{t}(u) are independent. Indeed, now ψ_{i}(u) = logit(u_{i}) so

X_{t}(u) = ([u_{1} ≤ u_{(t)}], [u_{2} ≤ u_{(t)}], …, [u_{n} ≤ u_{(t)}])

which is clearly stochastically independent of (u_{(1)}, …, u_{(n)}). Thus samples from the conditional distribution are given as X_{t}(u), and it is seen that our algorithm simplifies to a standard one.

It is interesting to note that the algorithm derived in this example also solves the following nonstatistical problem. Suppose that X = (X_{1}, …, X_{n}) are independent random variabels where X_{i} is Bernoulli-distributed with P (X_{i} = 1) = p_{i}, i = 1, …, n, and that the conditional distribution of X given T = ∑_{i=1}^{n} X_{i} is wanted. Then with z_{i} = logit(p_{i}) for i = 1, …, n and δ_{0} = 1 it is seen that this corresponds to θ = 0 in the above statistical setup. The suggested algorithm can therefore be used to solve the present problem as well.

[Censored exponential varables]

In this example we show how Theorem 1 (or 2?) can be applied in a case where the equation τ(u, θ) = t may have no solution for θ for given u and t.

Let n units with potential lifetimes Y_{1}, Y_{2}, …, Y_{n} be observed from time 0, but assume that the observation of the ith unit is censored at a given time c_{i} > 0. This means that we observe only X_{i} = min(Y_{i}, c_{i}), i = 1, …, n. Suppose Y_{1}, …, Y_{n} are i.i.d. exponentially distributed with expected value θ. Then the likelihood of X_{1}, …, X_{n} can be written (1/θ^{R}) exp(−S/θ) where R = ∑_{i}[X_{i} < c_{i}] is the number of noncensored observations and S = ∑_{i} X_{i} is the sum of all observations. Here T = (R, S) is sufficient compared to X = (X_{1}, …, X_{n}).

It should be remarked that the potential censoring times c_{i} are assumed known also for the units where X_{i} < c_{i}. For example this is the case if n machines are observed from possibly different starting points in time, and until a common end point. Let c_{1}, …, c_{n} be fixed, known numbers in the following.

In order to simulate X for a given θ we let U = (U_{1}, …, U_{n}) be a vector of n i.i.d. exponential random variables with expected value 1 and simulate X by means of χ(U, θ) = (η_{1}(U_{1}, θ), …, η_{n}(U_{n}, θ)) where

η_{i}(U_{i}, θ) = min(θU_{i}, c_{i}), i = 1, …, n

Thus T is simulated by τ(U, θ) = (γ(U, θ), ψ(U, θ)) where γ(U, θ) = ∑_{i}[θU_{i} < c_{i}] and ψ(U, θ) = ∑_{i} η_{i}(U_{i}, θ).

Let now t = (r, s) be fixed with 0 < r ≤ n, 0 < s < ∑_{i} c_{i} and consider the equations γ(u, θ) = r, ψ(u, θ) = s for a given u. Since ψ(u, θ) is stricly increasing in θ, from 0 to ∑_{i} c_{i}, there is a unique θ̃ which satisfies ψ(u, θ̃) = s. However, this θ̃ may not solve γ(u, θ) = r, even though the solution for θ of the latter equation is an interval. In the cases where indeed γ(u, θ̃) = r, let I(u) = {i_{1}, …, i_{r}} be the set of indices i for which θ̃u_{i} < c_{i}. With this notation we can then express the solution θ̃ as

θ̂(u, (r, s)) = \frac{s − ∑_{i∉I(u)} c_{i}}{∑_{i∈I(u)} u_{i}}  (18)

We shall use Theorem 1 for computation of conditional expectations E^{θ}(ϕ(X)|T = t). First we need to equip Ω_{Θ} = R^{+} with a measure. Let this measure be given by a density f(θ) with respect to Lebesgue measure. We then need to find the density Y_{t}(u) of τ(u, Θ) with respect to the product measure of counting measure and Lebesgue measure. This is found from

g(u, (r, s))dr ds = M_{Θ}{θ : γ(u, θ) = r, s ≤ ψ(u, θ) ≤ s + ds}
= \{ \begin{matrix}0 & if no solution of τ(u, θ) = (r, s) \\ f(θ̂(u, (r, s)))ds/ ∑_{i∈I(u)} u_{I} & otherwise \end{matrix}

Furthermore, the W_{t}(u) in Theorem ?? is given by

W_{t}(u) = \{ \begin{matrix}0 & if no solution of τ(u, θ) = (r, s) \\ X_{t}(u) ≡ χ(u, θ̂(u, (r, s)) & otherwise \end{matrix}

The choice f(θ) = 1/θ yields in particular the simple weight function

Y_{t}(u) = (s − ∑_{i∉I(u)} c_{i})^{−1}

An important special case is when the c_{i} are all equal. In this case we can sample directly from the conditional distribution of X given T = t by using X_{t}(u) = χ(u, θ̂(u, (r, s)) whenever θ̂(u, (r, s)) is defined.

It seems clear that the method of this example may be ineffective if n is large, since the probability of drawing a u for which a θ̂ exists may then be small.

[Multiple solutions for θ]

Consider the distribution on R with density function

f(x, θ) = \{ \begin{matrix}\frac{θx^{θ−1}}{a^{θ}+b^{θ}} & if 0 < x < a \\ \frac{θ(x−a)^{θ−1}}{a^{θ}+b^{θ}} & if a < x < a + b \end{matrix}

where a, b > 0 are given constants and θ > 0 is the unknown parameter. It is straightforward to show that if X = (X_{1}, …, X_{n}) is a sample from this distribution, then

T = ∑_{X_{i}<a} log X_{i} + ∑_{X_{i}>a} log(X_{i} − a)

is sufficient for θ.

The given density, and more general versions of it, appear naturally when we study several independent nonhomogeneous Poisson processes with a common intensity function λ(t) = αθt^{θ−1} (the so called power law process, see e.g. Crowder et al. 198x), but observed over time intervals with different lengths. As an example, consider two processes, observed over time intervals (0, a) and (0, b), respectively. Suppose we glue together the processes into one process observed on (0, a + b). Then, if we condition on the total number n of events, the event times in the glued process are distributed as the order statistic of n i.i.d. observations from the density f(x, θ).

For given values of θ we can simulate X from the density f(x, θ) by ordinary inversion using U with uniform distribution on (0, 1). In fact we then simulate X_{i} by

η(u_{i}, θ) = \{ \begin{matrix}u_{i}(a^{θ} + b^{θ}))^{1/θ} & if u_{i} < a^{θ}/(a^{θ} + b^{θ}) \\ u_{i}(a^{θ} + b^{θ}) − a^{θ})^{1/θ} + a & if u_{i} > a^{θ}/(a^{θ} + b^{θ}) \end{matrix}

By the above definition of T it is from this also clear how to define τ(u, θ) as an appropriate sum of terms of the form η(u_{i}, θ).

The point of this example is to show that the equation τ(u, θ) = t may have a varying number of solutions for θ, depending on the values of u and t. As an illustration we have in Figure xxx plotted the function τ(u, θ) as a function of θ when n = 2, a = 3, b = 1 and u = (0.5, 0.9). It is seen that there are two solutions for θ if t < 0.2, one solution if t > 1.9 and no solution if 0.2 < t < 1.9. As an example, for t = 0 we get the solutions θ_{1} = 0.925, θ_{2} = 1.801. Thus the formula of “Corollary” is needed for computation of conditional expectations.

7 Concluding remarks

Si noe om hva vi har gjort....

Andre ting (skal ikke nummereres nødvendigvis i den endelige versjon):

1. The results obtained in this paper are believed to be of greatest value in small sample cases where limit results based on asymptotic normality may not be trustworthy. Note that bootstrapping may be a reasonable “intermediate” alternative in many cases. For example, if θ is a nuisance parameter under the null hypothesis H_{0} of a testing problem, the approach of the present paper amounts to sampling conditionally on a sufficient statistic under H_{0}, while bootstrapping means to sample from an estimated distribution under H_{0}. The latter approach is thus not exact, but may work well at least for for moderate samples (Efron and Tibshirani, 1993).

2. In particular we have shown that the simple algorithm of the introduction holds if θ̂ is unique and if the model in addition has a certain pivotal structure. Remark that, as noted by Mood et al. (1974), it is always possible to find a pivotal quantity when X_{1}, …, X_{n} is a sample from a continuous distribution F (x, θ); simply by taking −2 ∑_{i=1}^{n} log F (X_{i}, θ)) which is chi-square distributed with 2n degrees of freedom. However, as stated in Theorem ?? our pivotal condition is stronger since it requires that the sufficient statistic can be recovered from the pivotal quantity. For example, this would not be the case in Example 6.

3. In the introduction we mentioned the possibility of constructing statistical problems from probabilistic ones. This can typically be done by introducing exponential families with the correct sufficient statistic. The main problem in practice is, however, to define a simulation algorithm based on a suitable U as in this paper. This problem is discussed for example in Bølviken and Skovlund (1996).

4. The results of Section ?? may be used in theoretical derivations of conditional distributions in particular cases. As an example, in the normal case of Example 6 it follows easily from our setup that the conditional distribution of X given X̄, S is a multivariate t-distribution. Also, Theorem ?? may be applied to find pivotal quantitites, and also to prove independence results.

5. Consider now a Bayesian setting where T is an observed sufficient statistic. Suppose there is given a prior distribtuion for the parameter θ, which we identify with M_{Θ}. In the notation of the present paper, the posterior distribution of Θ given the observation T = t equals the conditional distribution of Θ given τ(U, Θ) = t. The latter distribtuion is therefore the marginal distribtuion of Θ in the distribution given by Theorem 2new. For example, under the uniqueness conditions of Theorem ..., the posterior expectation of a function h(Θ) would be

E(h(Θ)|τ(U, Θ) = t)\frac{∫ h(θ̂(u, t))|\frac{f(θ)}{det ∂_{θ}τ(u,θ)}|_{θ=θ̂(u,t)}M_{U}(du)}{∫ |\frac{f(θ)}{det ∂_{θ}τ(u,θ)}|_{θ=θ̂(u,t)}M_{U}(du)}

where f(θ) is interpreted as the prior density for θ.

In the special case of normal distributions as given in Example 6 we get the posterior expectation of a function h(μ, σ) to be given by

\frac{E(h(x̄ − \frac{Ū}{S_{U}} s, \frac{1}{S_{U}} s)\frac{1}{S_{U}} f(x̄ − \frac{Ū}{S_{U}} s, \frac{1}{S_{U}} s))}{E(\frac{1}{S_{U}} f(x̄ − \frac{Ū}{S_{U}} s, \frac{1}{S_{U}} s))}

This expression exactly corresponds to evaluating the posterior expectation by ordinary importance sampling, drawing (μ, σ) from the distribution of θ̂(U, t) ≡ (x̄− \frac{Ū}{S_{U}} s, \frac{1}{S_{U}} s) which turns out to equal Fisher’s fiducial distribution for this case (see e.g. Welsh, 1996).

The formula can be easily evaluated by simulation by noting that Ū and S_{U} are independent, Ū ∼ N(0, 1/n) and (n − 1)S_{U}^{2} ∼ χ_{n−1}^{2}.

Bo Lindqvist Department of mathematical sciences Norwegian University of Science and Technology N-7034 TRONDHEIM

NORWAY Gunnar Taraldsen Department of Physiology and Biomedical Engineering Norwegian University of Science and Technology Center of Medicine and Technology N-7005 TRONDHEIM

NORWAY