Image Anal Stereol 2006;25:155-163 Original Research Paper FURTHER RESULTS ON VARIANCES OF LOCAL STEREOLOGICAL ESTIMATORS Z bynEk P awlas1 and E va E. V edel J ensen2 Department of Probability and Mathematical Statistics, Faculty of Mathematics and Physics, Charles University, Sokolovska´ 83, 18675 Prague, Czech Republic, Department of Mathematical Sciences, University of Aarhus, Ny Munkegade, DK-8000 Aarhus C, Denmark e-mail: pawlas@karlin.mff.cuni.cz, eva@imf.au.dk (Accepted October 9, 2006) ABSTRACT In the present paper the statistical properties of local stereological estimators of particle volume are studied. It is shown that the variance of the estimators can be decomposed into the variance due to the local stereological estimation procedure and the variance due to the variability in the particle population. It turns out that these two variance components can be estimated separately, from sectional data. We present further results on the variances that can be used to determine the variance by numerical integration for particular choices of particle shapes. Keywords: local stereology, marked point process, model-based setting, star-shaped set, variance. INTRODUCTION MARKED POINT PROCESSES One of the important unsolved problems in stereology concerns the stereological estimation of particle size distributions without specific assumptions about particle shape. It has been known for some time how to estimate stereologically the mean particle volume for particles of varying shape, cf. Jensen (1998). The resulting distribution of estimated particle volumes has been used as an estimate of the distribution of the true particle volumes. It is clearly important to be able to judge when such a procedure is justified. The particular case of estimating the volume-weighted mean particle volume has recently been treated in Cabo et al. (2003). It is here shown that an estimator based on planar observation is one order of magnitude more efficient than the traditional one based on observation along lines. In the present paper, we concentrate on the ordinary (unweighted) particle volume distribution. It is shown that an estimator of mean particle volume based on planar observations is superior to one based on line observations, especially for elongated particles. The variance of the estimator can be decomposed into the variance due to the stereological estimation procedure and the variance due to the variability in the particle volumes. We will show how to estimate these variance components separately, from sectional data. If the variance due to the stereological estimation procedure is small compared to the variance due to the variability in the particle volumes, the distribution of estimated particle volumes can be regarded as an estimate of the distribution of the true particle volumes. We define the particle model by means of marked point processes. For more details, we refer to Stoyan et al. (1995). Let Ym = {[Xi;Si]} be a marked point process such that Xi is a point in Rn and Si belongs to the space JCd of d-dimensional differentiable manifolds in Rn with finite d-dimensional Hausdorff measure and with the reference point at the origin O. The point Xi then serves as a reference point of Xi + Si, the i-th particle. The marked point process xFm will be assumed to be stationary, i.e., ^m+x = {[Xi+x;Si]} has the same distribution as xFm for every xeR n. Stationarity of Ym implies stationarity of the unmarked point process Y = {Xi}. Denote by A its intensity and assume that 0 < X < oo. Let SO(n,Lr) be the subgroup of SO(n) consisting of rotations keeping an r-dimensional linear subspace Lr fixed. Then, *Fm is said to be invariant under rotations in SO(n,Lr) if BYm = {[BXi;BSi]} and xFm have the same distribution for all rotations BeSO(n,Lr). The intensity measure of the marked point process is defined for A G ßg{R n) and U G 3S{JKd) as Am(AxU) = Eld1A(Xi)1U(Ei), i where 1A(-) stands for the indicator function of the set A and 3ë denotes the Borel c7-algebra. The stationarity 155 Pawlas Z et al: Variances of local stereological estimators of Ym implies the following decomposition Am(AxU)=XV{A)Pm(U), where V = Xn is the volume and Pm is the mark distribution. By S0 we denote a random manifold with distribution Pm. If Ym is invariant under rotations in SO(n,Lr), then BZ0 has the same distribution as S0 for all B e SO(n,Lr). THE LOCAL STEREOLOGICAL ESTIMATORS The local stereological estimators are based on information collected from section planes in Rn through a reference point of the particle. In this section, we present the actual form of local stereological estimators of Hausdorff measure for a generic particle K G Md. Then a section plane of dimension p is a p-dimensional linear subspace (for brevity called p-subspace) of Rn, p = 0,1,...,n. For comprehensive exposition of local stereology, see Jensen (1998). There are various forms of the local estimators, depending on the restriction put on the p-subspace. Denote by J^pnLr the set of p-subspaces containing a fixed r-subspace Lr, 0 < r < p < n. Let And be the d-dimensional Hausdorff measure in Rn and let us use the short notation dxd instead of Xd(dx). Note that the ordinary Lebesgue measure is Ânn = Xn. For K G Md, the local stereological estimator of knd(K), based on a p-subspace Lp G JžfpnL , d - n + p > 0, has the form, cf. Jensen (1998), (5.24), r m (n,d) p,Lr (K,Lp) ||_ un—p On-r Gp-r KnLp xG(Tan[K,x],Lp)-1d xd-n+p (1) where an = 2nnl2/T{n/2) is the surface area of the unit sphere Sn"1 in Rn, nL±_ is the orthogonal projection onto L^r and G(Tan[K,x],Lp) is the so-called G-factor defined in Jensen (1998), Definition 2.9. Note that in the case d = n, G(Tan[K,x],Lp) = 1 for arbitrary Lp. In a design-based setting, Eq. 1 is an unbiased estimator of Xnd(K). Thus, let \in L be the unique probability measure on Jžfpn Lr, invariant under rotations from SO(n,Lr). In what follows, we will write dLnpLr as short for /xpnL (dLp). By an isotropic p-subspace in Rn, containing r the fixed r-subspace Lr, we mean a random p-subspace with constant density with respect to \inpLr . Then, if K satisfies the regularity conditions stated in Jensen (1998), Proposition 5.4, and Lp is an isotropic p-subspace containing Lr, the local estimator m(pn,L,dr)(K,L~ p) is an unbiased estimator of ?nd(K), i.e., Ind (K) = j n m{pn L dr {K,Lp) dLnpLr (2) Example 1. For K G M in R3 there are three local stereological estimators of the volume V(K), m m K,L2) = 2 \\x\\dx2 KL \i'')(K,L1) = 2n \\xfdx1 KL (3,3) ( 2,O (3,3) KC\L2 Kr\Lx (3) (4) m2L(K,L) = n L2) = Tz \\TZLLx\\dx2. (5) KnL2 l The estimators Eq. 3 and Eq. 4 are related by a so-called Rao-Blackwell procedure. We have m(2 3,O,3)(K,L , Li) = E mfO (K, L ~ )\L\ , (6) if L ~ is an isotropic line in L2. For later reference, we also present the local estimator of 1nd{K)2, {n 0 and Vq(yh...,yq) denotes the q-dimensional Hausdorff measure of the parallelepiped spanned by yh...,yq. Assuming that K satisfies the regularity conditions {n,d) from Jensen (1998), Theorem 5.6, m~ (pn,L,dr)(K,L~ p) is an p unbiased estimator of ?nd(K)2 if L~ p is an isotropic p-subspace, containing Lr. Example 2. For d = n = 3, p = 2 and r = 0, the estimator of V(K)2 has the form m ~ (3,3) ( 2,O K,L2)=2? ?2(x,y)dy2dx2 KL KL and ?2(x,y) is twice the area of the triangle with vertices O, x and y. _ 156 Image Anal Stereol 2006;25:155-163 THE VARIANCE OF LOCAL ESTIMATORS We will now return to the model-based case. We let S0 be a generic random particle, distributed according to Pm and denote by Em the expectation with respect to this distribution. Let Lp{0) be a fixed p-subspace in Rn, containing an r-subspace Lr, 0 < r < p < n, d-n + p>0. Below, we give explicit results for the second moment of mpn , L , dr (Z0,Lp{0)). For this purpose, the following proposition is very useful. Proposition 1. Let Pm be invariant under rotations in SO(n,Lr). The estimator mpn, L , dr {^0,Lp{0)) has the same distribution as mpn, L , dr(Z0,Lp), where Lp{0) G Jžfpn , Lr is fixed, Lp is an isotropic p-subspace containing Lr and S0 and Lp are independent. Proof. The result follows from the fact that for any non-negative measurable function h, Emhm (n,d) Sn,L p,Lr -0 ,L p(0) Ln p,Lr hm (n,d) S,L p,Lr \^0 ,L p )) dLp,Lr The left-hand side can be rewritten using the invariance of Pm under rotations in SO(n,Lr) and Jensen (1998), Lemma 8.4, Emh(mpn , L , drf{Z0,Lp 0 = Emhmp,Lr (B E0,Lp 0) = Emh(mpn , L , drf(E0,BLp 0) where B G SO(n,Lr). From invariant measure theory there exists an invariant probability measure an{r) on SO(n,Lr). Note that JžfpnL = {BLp{0) : B G SO(n,Lr)} and f g(BLp 0)anr)(dB)= [ g{Lp)dLnpLr \ , r) p,Lr for any non-negative measurable function g on Jžf pnL . Thus, Emh{m{pn , L , dr [Z0,Lp0 Emh(m{n , L , dr\z0,BLp 0) anr){dB) n,Lr) SO(n,Lr) O(n,Lr) Em h m(pn,L,dr)(?0,BLp(0)) ?(nr)(dB) SO(n,Lr) Em n h m(pn,L,dr)(?0,Lp) dLnp,Lr D When convenient we use the short notation m(S0) for mpn, L , drf {Z0,Lp 0) and m(E0,Lp) for mpn, L , drf(E0,Lp). Using Proposition 1 and Eq. 2, we get Emm(S0) = Em(Z0,Lp) = EmAnd(S0) . (8) Moreover, the relation Eq. 2 for a fixed S0 can be written as E[m(Z0,Lp)\Z0] =Ad (S0), almost surely, and for the variance of m(S0), we get varmm(S0) =varm(S0,Lp) = Em var [m(Z0,Lp) I S0] + varmE [m(Z0,Lp) I S0] = Em var [m(Z0,Lp) | S0] + varm Ad (S0) >varmAnd(S0) . Hence, varmAnd(S0) 0 such that m(K0,Lp{0)) / Xnd(K0) for K0 G A(Lp{0)). For any Lp G Jžfpn , Lr, we can find B G SO(n,Lr) such that Lp =BLp{0}. Since m(K0,Lp{0}) = _ _ 157 Pawlas Z et al: Variances of local stereological estimators m(BKo,BLm), see Jensen (1998), Lemma 8.4, we have m(BK0,Lp) / Xnd(K0) for K0 G A(Lp{0)). Therefore, m(K0,Lp) / Xnd(K0) on the set A = {(K0,Lp) : Ko G A(Lp),Lp G Jžfpn L}, where ALp) = B-A Lpo)). But from the invariance of Pm under rotations in SO(n,Lr) we obtain Pm(A(Lp)) = Pm(ALpo))) > 0 which means that (Pm xupn L)(A)> 0 and this leads us to a contradiction. ' r D If varmm(20) = varmAnd (20), then the Hausdorff measure of the manifold is determined from the local section without error. Such local stereological estimators are exact, i.e., the variance of the estimator is created only by the randomness of particles. The simplest example of a particle with exact local volume estimator is a ball. Proposition 3. Let 20 be an n-dimensional ball in Rn centred at O with probability one. Then var m m n Lr (20 ,1p(0) varmV(20) (10) Proof. If K = bn(O,R) is an n-dimensional ball of radius R with centre at the origin, then for all Lp G Cp n =^n L , I \\nL^x\\n-p dxp kl r Kr\Lp ¦¦¦ (x2r+1 + ---+x2p)!-2Edxp---dx1 {x2+-+x2p) denote the radial function ofK, pK{(o) = max{c:c(oeK}, m G Sn"1 , cf. Gardner (1995), p. 18. Let for co G Sn-1, , , \pK(o)n + pK(-o)n forOeK, PnK{o) = {\\pK(mW-\pK(-coW\ forO^K, be the n-chord function of the set K at O, cf. Gardner (1995), Definition 6.1.1. Furthermore, let Pn,k{L 1 pj = — S pnK{(o)dcop-1 2p Sn-lnLp be the section function, cf. Gardner (1995), Chapter 7. Proposition 5. Let K be a star-shaped set at O. Then m(pn,L,nr)(K,Lp) = nzr 1 pntK{a)\\KLMrpdap-1 2n Sn-lnLp (11) Proof. Using the polar decomposition of Lebesgue measure we obtain / \\nLLx\\n-pdxp KL r Knspm{a} \\x\\n-l\\7iLMrpdxldcop-1 2n S pnK{(o)\\KLMrpd) = «-1/nmin , O , n (So, span{«}) co G Sn"1. Obviously, m(S0) = m(star(E0)), thus Proposition 6 and Proposition 7 can be used for any So if S0 is replaced by star(S0) in the right-hand side of Propositions 6 and 7. Now we turn to the case p > r + 2. Proposition 8. Forp> r+2the second moment of the local estimator is Emmpn , L , drf (So, Lpo))2 = On-rOp-r-l E n-r-lOpr ^ On I tzllx %L\_y Y^\,Y^y\\ p—r dxd dyd COn mSn , O , n (So,span{u)}) = 2 npn,So(ffl), co G Sn"1 . Using the symmetry of S0 (ps0(u>) = ps0(-u>)) and Proposition 1 we obtain the stated result. D Sometimes, it can be useful to have an alternative expression of Emm\n , O , n (S0, L1(0))2. Proposition 7. Let S0 be symmetric and star-shaped set at O. Then for the local estimator with d = n,p=1 andr = 0 we have Emm\n, n\z0,Lm)2 = 2conEm I \\x\\ndxn . Proof. The formula in Proposition 6 can be rewritten as Proof. Under the regularity conditions of (Jensen, 1998), Theorem 5.6, we know that n 2 r r Emm(S0)2 = ^P Em / / / W^\\n~p p-r ^LrZ0nLpz0nLp xG{Tan[Z0,x],Lp)-l\\7iLp\\n-pG{Tan[Z0,y],Lp)-1 Using the generalized Blaschke-Petkantschin formula (Jensen, 1998), Theorem 5.6, with g(x,y) llTrL xll^IlTrL yin -p ?2(?L?r x,?L?r y)n-p weget(d-n + p>0) Emm(S0)2 = «2Em / I nspanju)} On = nco2Em / / ||x||2n 1dx1dLnO = 2conEm I \\x\\ndxn , where in the last step we have used (Jensen, 1998), Proposition 4.1 with g{x) = \\x\\n. D Emmp,Lr (,o,Lp) = 0nri0pr II JL x11n—p II jL y\\n~p xEm.l .1 Z^x^lyy-p dxddyd . So S0 r r Notice that due to the assumed regularity conditions, V2(nLLx, nLi_y) = 0 on a set of (And x And)-measure zero and the integral on the right-hand side is well-defined. The result now follows immediately from V2(x,y) = (||x||2||;y||2-(x,y)2)1/2. D _ _ 159 Pawlas Z et al: Variances of local stereological estimators Note that the second moment of mpn Ld) (So,Lp)) does not depend on the G-factor. For n-dimensional particles the formula for the variance given in Proposition 7 was expressed through integrals over Sn-1 n Lj- in Jensen et al. (1999). Finally, we consider the special case of star-shaped particles and r = 0. Proposition 9. Let S0 be a star-shaped set at O with O G S0 almost surely. Suppose that p>2. Then (n,n) Em mp,O Zo,Lp(o)Y XEm PnßQ(Oh)pnßQ(2}2) ~^ dm?-1 dco^-1 Proof. It is not difficult to see that the result follows immediately from Proposition 8 and the following formulation of polar decomposition of Lebesgue measure f ( x \ k \\x\\ - 1 f ~ 2Sn-1 Sn-1 span{?} 1K{x)g{(o)\\x\\n-ldxld(o n-1 D EXAMPLES In this section we use the results to find explicit expressions of the variance of local stereological estimators for specific particle shapes. THE PLANAR CASE For n = d = 2, p = 1 and r = 0, the variance can easily be determined, using Proposition 6 or Proposition 7 for various shapes of ?0. We give the formulas for varm m(1 2,O,2)(?0,L1(0)) for three particular choices of ?0 centred at O: – rectangle with sides of lengths a, b Eab(a2 + b2)-(Eab)2 ellipse with semiaxes of lengths a, b it1-~2 Eab(a2 + b2)-7i2(Eab)2 equilateral triangle with the side of length a %yf 3(21- 2 8log2) Ea 4_3 (Ea 2 ) 2 TRIAXIAL ELLIPSOIDS We suppose that ?0 is an ellipsoid centred at O and with semiaxes of lengths a, b and c. In R3 there are three local stereological volume estimators, namely Eq. 3, Eq. 4 and Eq. 5. For p = 1 the second moment of m(1 3,O,3)(?0,L1(0)) can be written as (using Proposition 6 and spherical coordinates) 87T 9 o io % cos2 0 cos2

2p S2 S2 xpSo(«1)3pSo(«2)3d«12d«; Note that there is a mistake in Jensen (2000), the constant 4 should be replaced by -X. For fixed S0 the double integral can again be computed numerically. The formula Eq. 14 still holds, the values of k for several choices of ratios a0/b0 and b0/c0 are summarized in Table 1. We have also computed k for _ 160 Image Anal Stereol 2006;25:155-163 an intermediate estimator, usually called the nucleator, cf. Gundersen (1988), m (3'3)f«ni 1,O k-Oj 1 2 mfO (So, span{«i}) m (3,3) 1,O s0, span {^})] isotropic direction in an 2nL2 is orthogonal to «1 where «1 G S2 n L2 is an isotropic plane L2 and u>2 G S1 1L^ is orthogonal to un. Our approach based on numerical integration enables more precise results than those obtained by simulation in Jensen (2000). Note that Vk^1 is the coefficient of error of the local volume estimator for a corresponding ellipsoid with semiaxes a0, b0 and c0. Table 1. The values of k from Eq. 14 for three types of volume estimators and various shapes of ellipsoids. a0/b0 bo/co mi,O(S0) m —iiO(S0) m2,O(S0) 1 1 1 1 2 4 1 1.34440 2.43264 1 1.10854 1.58048 1 1.07945 1.26608 2 2 2 1 2 4 1.51757 2.60921 5.03481 1.16635 1.63012 2.79692 1.11835 1.31307 1.59320 4 4 4 1 2 4 4.42495 8.51728 16.87915 2.45371 4.44465 8.60156 1.59821 2.03585 2.56130 Higher values of k mean higher variance caused by the local stereological estimation. For ball (* = 1) we have an exact estimator with varmm(S0) =V02varp3 =varmA3(S0) . Note that the error is larger for prolate spheroids (b = c) than for the corresponding oblate spheroids (a = b). In view of Eq. 6, it is not surprising that smaller values of error are obtained for the local estimator based on plane sections. In the remainder of this subsection we consider the last local volume estimator Eq. 5. Obviously, it depends on the choice of the fixed line Lx (usually called vertical axis) relative to the ellipsoid. We assume that the vertical axis has the same direction as one of the semiaxes of the ellipsoid (say the one of length c). Then the profile S0 n L2(0) is a planar ellipse with semiaxes of length A and c. Hence, the local estimator has the following form m (3,3) S0,L2(0)' / / A2cr2|cos 3 we do not always have to use numerical integration in order to derive the explicit formula for the variance. For instance, if ?0 is an ellipsoid in R4 with semiaxes a1, a2, a3, a4, then Em m(1 4,O,4)(?0,L1 (0), 7T4Eala2a3a4 [ ^Xa i + 12 ia i=l i=l 8 This result can be derived from Proposition 7, using elliptical coordinates and lengthy but straightforward calculations. satisfy the regularity conditions from Jensen (1998), Theorem 5.6, and let Lp{0) be a fixed p-subspace in Rn, containing Lr, 1 < r+ 1 < p < n, d-n + p > 0. Proposition 8.7 in Jensen (1998) states that if W is a bounded Borel set with positive volume, then 2 A ?1W(Xi (n,d) m~ p,Lr (Xi +?i,Xi +Lp(0))/?(W) is a ratio-unbiased estimator of Em K (So)2. We can also estimate Emmpn L drf (Z0,Lp{0})2 by ESTIMATION OF VARIANCES It is interesting variances to find the estimators of the 2 a2 = varmmfpn L dr (EoLpo)) and o^ varm?nd(?0), separately, from sectional data. First we introduce ratio-unbiased estimators of µ= Em?nd(?0), am = (n,d) Em mp,Lr ^0,Lp(o) and ??2 Em A„ ( L, d2 respectively. Let W be a fixed bounded Borel set in Rn with positive volume. We consider a sample of particles Xi + Si with Xi in the sampling window W. Then the local estimator is determined from the central section (Xi + Si) n (Xi + Lp) for each sampled particle. In order to be in accordance with Eq. 1 we can think that the centred particle Si is sectioned by a fixed p-subspaceLpo) e^L,i.e., m (n,d) (n,d) p n d>(X i + Zi ,X i + L p{0) )=m p n d>(Zi ,L p{0) ) . In what follows we assume that the mark distribution Pm of a stationary marked point process Ym is invariant under rotations in SO(n,Lr). It can be shown that (see Jensen (1998), Proposition 8.5) I1 = ^WiX^mpn L driXi + ZuXi + Lp^/^iW) is a ratio-unbiased estimator of the mean particle Hausdorff measure Em Xnd (S0) = / Xdn (K) Pm (dK) . The proof is based on Campbell’s formula for marked point processes (see Stoyan et al. (1995), Section 4.2) and Eq. 8. Using Eq. 7, we can estimate the second moment in the mark distribution, i.e., EmAnd(S0)2. Let S0 This estimator is again ratio-unbiased, as can be easily seen from Campbell’s formula for marked point processes. The problem is how to estimate M2= ( Emm (n,d) So,L \2 / \2 )) ) = (Em\ (So)] p,Lr ^0 ,L (0) As long as we restrict ourselves to the case of independently marked point process (i.e., the ?i are independent and identically distributed and independent of ?), 1 ^p1i m Xi+si fiY (15) and 22----------fcg_(L)2 ^ (16) ^x-«x-(M) + xj,(W) _ 1 are unbiased estimators of <72 = varmm(So) and o\ = varmAnd(S0), respectively. In the general case of dependent particles we propose following estimators of cr2, o CT — lh(m(Xi + Ei)-m(Xj + Ej))2 (17) or the estimator taking into account edge corrections o CT — ( mjXi+ZiymjXj+Zj) ) 2 l*h V ( (W-Xi)n(W-Xj) ) 2Lh (18) V ( (W-Xi)n(W-Xj) ) where Lh means that the sum is taken over i and j such that Xi G W, Xj G W and ||Xi -Xj|| > h. An appropriate choice of the parameter h > 0 is some trade-off between small value (larger bias) and _ _ _ i 162 Image Anal Stereol 2006;25:155-163 large value (few observations, larger variance). Using Campbell’s formula it is easy to show that o| is in both cases ratio-unbiased if m(Xi + Si) and m(Xj + Sj) are independent whenever ||Xi -Xj\\ > h0. The estimate of of has then the form Gf = a|-a| + c^. (19) In applications, the distribution of estimated particle volumes (or other size parameters) has been used as an estimate of the true particle volume distribution. This procedure is justified if the variance due to the stereological estimation procedure is small compared to the variance due to the variability in the particle population. We can estimate both variances from central sections using Eq. 15 and Eq. 16 or Eq. 17, Eq. 18 and Eq. 19. If the estimates are closed we can expect that the distribution of estimated sizes will be close to the true size distribution. The practical implications of this observation will be investigated elsewhere. ACKNOWLEDGEMENTS This research has been supported by a Marie Curie Fellowship of the European Community Programme Improving the Human Research Potential and the Socioeconomic Knowledge Base under contract number HPMT-CT-2001-00364, the Danish Natural Science Research Council, the Grant Agency of the Czech Republic, project 201/06/0302 and the Czech Ministry of Education, project MSM 0021620839. The work has been presented at the Conference S4G, Prague, June 26–29, 2006. REFERENCES Cabo A, Baddeley A (2003). Estimation of mean particle volume using the set covariance function. Adv Appl Prob 35:27-46. Gardner RJ (1995). Geometric tomography. New York: Cambridge University Press. Gundersen HJG (1988). The nucleator. J Microsc 143:3-45. Jensen EBV (1998). Local stereology. World Scientific, Singapore. Jensen EBV (2000). On the variance of local stereological volume estimators. Image Anal Stereol 19:15-18. Jensen EBV, Petersen L (1999). When are local stereological volume estimators exact? Research report no. 4, Laboratory for Computational Stochastics, University of Aarhus. Stoyan D, Kendall WS, Mecke J (1995). Stochastic geometry and its applications. 2nd ed. New York: John Wiley. 163