Image Anal Stereol 2013;32:117-125 Original Research Paper LEVY-BASED ERROR PREDICTION IN CIRCULAR SYSTEMATIC SAMPLING Kristjana Yr Jonsdottirb and Eva B. Vedel Jensen Department of Mathematics, Aarhus University, Ny Munkegade 118, 8000 Aarhus C, Denmark e-mail: kyj@imf.au.dk, eva@imf.au.dk (Received October 8, 2012; revised May 14, 2013; accepted May 15, 2013) ABSTRACT In the present paper, L6vy-based error prediction in circular systematic sampling is developed. A modelbased statistical setting as in Hobolth and Jensen (2002) is used, but the assumption that the measurement function is Gaussian is relaxed. The measurement function is represented as a periodic stationary stochastic process X obtained by a kernel smoothing of a L6vy basis. The process X may have an arbitrary covariance function. The distribution of the error predictor, based on measurements in n systematic directions is derived. Statistical inference is developed for the model parameters in the case where the covariance function follows the celebrated p-order covariance model. Keywords: Fourier series, L6vy basis, planar particles, stationary stochastic processes, stereology, systematic sampling. INTRODUCTION A long-standing problem in stereology is variance estimation in systematic sampling. One class of problems involves estimation of an integral of the form f 2n Q = x(e ) de Jo (1) where x(d) is an integrable function on [0,2n), called the measurement function. The estimator typically considered is based on circular systematic sampling and takes the form - 2n n-1 ^ 2n i. Qn = — V x(0 + —) , n > 1 , 1-1 1-1 i=0 where 0 is uniformly distributed in [0,2n). For instance, if Y is a bounded convex planar set containing the origin O, examples of Eq. 1 are x(e ) 1 r(0)2 if Q = area of Y, h(0) if Q = boundary length of Y, where r(d) and h(d) are the radial function and the support function of Y in direction d, respectively. The geometric identity (Eq. 1) is in these cases a consequence of polar decomposition in the plane and an identity for mean width (Schneider, 1993, Eq. 5.3.12). If instead Y is a bounded convex spatial set containing O, the volume and the surface area of Y may be estimated by a two-step procedure which involves circular systematic sampling in a section through O and the use of the cubed radial function or the squared support function (Gundersen, 1988; Cruz-Orive, 2005). Yet another example is volume estimation by the so-called vertical rotator (Jensen and Gundersen, 1993). In Gual-Arnau and Cruz-Orive (2000), a design-based procedure of approximating the variance of Qn, based on modelling the covariogram of x(d) by a polynomial model, is developed. Hobolth and Jensen (2002) consider a model-based procedure, where the measurement function is assumed to be a realization of a periodic stationary stochastic Gaussian process X = {X(0) : 0 e [0,2n)}. It is shown in Hobolth and Jensen (2002) that the covariogram model considered in the paper by Gual-Arnau and Cruz-Orive (2000) is a special case of a p-order covariance model for the stochastic process X in the model-based set-up. The p-order covariance model is given by c(e ) = cov(x (e ), x (o)) = h + V ^s cos(se ), s=2 V1 = a + ß (s2p - 22p) , s > 2 , (2) where the model parameters satisfy p > 1/2, a, p > 0. In Hobolth et al. (2003), this parametric covariance function has been used in the modelling of the radial function of a random star-shaped planar particle. In this case, p determines the smoothness of the particles boundary while a and p determine the 'global' and the 'local' shape of the particle, respectively. (Note that in Eq. 2, Ai is set to zero which ensures that the reference point of the particle is approximately the centre of mass.) The model-based counterpart of the design-based methodology provided in Gual-Arnau and Cruz-Orive (2000) was further developed in Jónsdóttir et al. (2006), where the general form of the p-order covariance model (Eq. 2) was used to obtain a more accurate approximation of the prediction error E(Qn - Q)2. In Hobolth and Jensen (2002) and Jónsdóttir et al. (2006), the process X is assumed to be Gaussian. Motivated by the fact that powers of the radial function and the support function are used in practice, we will in this paper consider non-Gaussian models, obtained as a kernel smoothing of a so-called Lévy basis. As we will show, it is possible under the Lévy-based model to derive the distribution of the error predictor Qn — Q which may be markedly non-Gaussian for the moderate sizes of n used in practice. Lévy-based modelling has been popular in recent years, e.g. in the modelling of turbulent flows, spatio-temporal growth, spatial point processes and random fields (Barndorff-Nielsen and Schmiegel, 2004; Jónsdóttir et al., 2008; Hellmund et al., 2008; Jónsdóttir et al., 2013). More specifically, we will consider stochastic processes of the form f 2n X(6)= m +/ k(6 — 0)Z(d0) , 6 e [0,2n) , Jg where m determines the mean of the process, Z is a homogeneous and factorizable Lévy basis on [0,2n) and k is a deterministic kernel function. In principle, any covariance model, including the p-order covariance model, can be induced under this modelling framework, by assuming a specific form of the kernel function (see the next section). Under the p-order model, it is easy to control the local and global fluctuations of the stochastic process X. The Lévy-based models with p-order covariance thus constitute a flexible and tractable model class. In particular, this model class has more structure than the non-Gaussian models considered in Hobolth et al. (2003) and this allows us to derive distributional results. The composition of the remaining part of the paper is as follows. First, a theoretical background for stationary periodic processes with period 2n, based on kernel smoothing of a Lévy basis, is given. Then, estimation of E(Qn — Q)2 under the general Lévy-based model is discussed. The distribution of the error predictor Qn — Q under the Lévy-based model is derived, and it is shown how to estimate this distribution. An example of random particles simulated from a Lévy-based model is given together with the distribution of the n-point area estimator of these particles. Finally, a discussion is provided. Some technical derivations are deferred to an appendix. LÉVY-BASED STOCHASTIC PROCESSES ON THE CIRCLE This section provides an overview of stationary periodic processes on [0,2n) based on integration with respect to a Lévy basis. For further details on the general theory on Lévy bases, in particular, the integration with respect to a Lévy basis, the reader is referred to Barndorff-Nielsen and Schmiegel (2004) and Hellmund et al. (2008). Let X = {X(6) : 6 e [0,2n)} be a 2n periodic stationary stochastic process on [0,2n), given by f 2n X(6) = m + / k(6 — 0)Z(d0) , 6 e [0,2n) , (3) JG where m determines the mean of the process, Z is a homogeneous and factorizable Lévy basis on [0,2n) and k is an even periodic kernel function with period 2n and a Fourier representation k(6) = + £ ^scos(s6) . (4) s= 1 The model (Eq. 3) is the continuous analogue of the following discrete model X(6)= m + £k(6 — 0)Z(0) , 0 where the sum is over an equally spaced set of angles and the random variables Z(0) are independent and identically distributed, the common distribution being infinitely divisible. The integral in Eq. 3 is formally defined as a limit in probability (Rajput and Rosinski, 1989). A spatio-temporal version of Eq. 3 has previously been considered in Jónsdóttir et al. (2008). The Lévy basis Z is extended by Z(A + 2nm) = Z(A) for all m e Z and all Borel sets A e B([0,2n)). A Lévy basis has the property that Z(A1),...,Z(An) are independent when A1,..., An e B([0,2n)) are disjoint and Z(A) is infinitely divisible for any A e B([0,2n)). The assumption of homogeneity implies that all the finite-dimensional distributions of Z are translation invariant. If Z is Gaussian, the integral in Eq. 3 exists if k is L2-integrable with respect to the Lebesgue measure on [0,2n). When Z is a so-called Lévy jump basis (e.g., Gamma or inverse Gaussian Lévy basis), the integral exists if k is integrable with respect to the Lebesgue measure on [0,2n) and if /R|r|V(dr) < where V is the Lévy measure associated with Z. These results follow from Hellmund et al. (2008, Lemma 1). An important entity associated with a Levy basis is its spot variable 7' which is infinitely divisible. Without loss of generality we will in what follows assume that the spot variable 7' is centered, 7' = W — E(W), where W is an infinitely divisible random variable. The following theorem characterizes the distribution of the stochastic variable X(6). Theorem 1. The stochastic variable X(6) has the cumulant generating function kx (t ) = log E(etX (6)) = t(M — 2rc&E(W))+/ kw(tk($))d$ , Jo (5) where KW (t) is the cumulant generating function ofW. Moreover, the derivatives of KX (t) are given by (when they exist) KX(t)= M — 2rc§)E(W)+ / KW(tk($))k($) d$ , o K{X](t) = rK^)(tk($))k($)rd$ , r > 2 , (6) ) (r) where K'^)(t) denotes the r'th derivative ofKW(t). Proof. The result is obtained by using that the cumulant generating function of the integral f • 7 = /02n f (6)7(d6) of a function f with respect to a homogeneous and factorizable Levy basis 7 is given by f 2n Kf.z(t) = J kz(tf (6)) d6 , (7) where 7' is the spot variable associated with 7, cf. Hellmund et al. (2008, Eq. 10). □ the underlying spot variable, e.g., when W is Gamma distributed, X is called a Gamma Levy process, irrespectively of the choice of the kernel function. We will typically assume that k2(w) = 1, i.e., the skewness and kurtosis of W are equal to the third and fourth cumulant, respectively. In Table 1, we give the cumulant generating function, third and fourth cumulants of W for the three distributions mentioned above. Note that as W has unit variance, the Gamma and inverse Gaussian Levy bases are only determined by a single parameter n > 0. The Levy measures V of the Gamma and inverse Gaussian Levy bases satisfy the condition /R |r|V(dr) < ~ for the existence of the integral in Eq. 3, cf. Hellmund et al. (2008, Example 3). In principle, any covariance function can be modelled within this set-up. This can be seen, using the theorem below. Theorem 2. The stochastic process X has a mean value M and a covariance function c(6) = A0 + £ Ascos(s6) , 6 e [0,2n), s=1 where A0 = 2n§?K2(W) , As = ^s2K2(W) , s > 1 . (8) Proof. Using that c(6 )= Cov(X (6), X (0)) = K2(W) k(6 — $)k—) d0 , 0 we easily obtain Eq. 8. □ It follows that the cumulants of the stochastic variable X(6 ) are given by K1(X (6 )) = M, f 2n Kr(X(6)) = Kr(W) k($)rd$ , r > 2 , 0 where Kr (W) denotes the r-th cumulant of the stochastic variable W. Possible choices of the distribution of W are the Gaussian, Gamma and inverse Gaussian distributions. When the kernel function is proportional to an indicator function, k(6) = c1A(6) for A e B([0,2n)), the marginal distribution of X(6) will be of the same type as that of W. Otherwise, the marginal distribution will not be as simple, but the process X will inherent the name of the distribution of Using Theorem 2, we can construct the candidate kernel k that induces a given covariance function c. For instance, if c follows the p-order model (Eq. 2), then k is of the form specified in Eq. 4 with Ao 5o = 2nK2(W) ' §1 = 0, ^KK2(W)[a + ß(s2p -22p)] ' s> 2 . This choice of kernel will for a Gaussian Levy basis give a well-defined integral in Eq. 3 if p > 1/2, while for a Gamma or an inverse Gaussian basis the integral is well-defined if p > 1. 1 § s Table 1. Examples of the distribution of W, together with the corresponding cumulant generating function, and the third and fourth cumulants. distribution Gaussian Gamma Inverse Gaussian w n (0, i) r(n, yn)_ig(n3, n ) Kw(t) 12/2 -niog(i -t/yn) n4(i - Vi -2t/n2) K3(w ) 0 2/yn 3/n2 K4(w )_o_6/n_15/n4 ESTIMATING E(Qn - Q)2 UNDER THE LEVY-BASED MODEL In Hobolth and Jensen (2002), the focus was on the prediction error E(Qn — Q)2. If the covariance function of X has the following Fourier expansion c(6) = Xo + £ As cos(s6) , 6 e [0,2n) , s= 1 then it was shown in Hobolth and Jensen (2002) that E(Qn — Q)2 = £ Xnk , (9) k=1 where Xnk is the Fourier coefficient of order n • k in the Fourier expansion of the covariance function of X. Note that in Hobolth and Jensen (2002), circular systematic sampling on [0,1) instead of [0,2n) is considered, so Eq. 9 represents an adjusted version of Hobolth and Jensen (2002, Eq. 8). In Hobolth et al. (2003), a procedure for estimating the prediction error under a Gaussian p-order model was developed, based on a Fourier expansion of X X(6)=A0 + £ (As cos(s6)+ Bs sin(s6)) , s=1 6 e [0,2n) , where =d means equality in distribution and 1 r2n A0 = — X(6) d6 , 2n J0 1 r2n As = -/ X(6) cos(s6) d6 , n J0 1 c2n Bs = -/ X(6) sin(s6)d6 . n J0 When X is a periodic stationary Gaussian process, the Fourier coefficients of X become independent and normally distributed, As ~ Bs ~ N(0, As). As suggested in Hobolth et al. (2003), the parameters a and p in the p-order model can then be estimated using maximum likelihood estimation based on the first S Fourier coefficients, r , m A 1 ( (a2 + b2) A Lo,S(a,p) = n^w—^exp -^h—À , , v s=2 2nXs(a,p) V 2Xs(a,p)J (10) where Xs(a,p), s = 2,...,S, satisfy (2) and as and bs, s = 2,..., S, denote discretized Fourier coefficients of X. In this section, we will show that this procedure can be used under the general Lévy-based model. The following theorem gives the distribution of the Fourier coefficients and their relations under the general Lévy-based model. In the Gaussian case, the theorem can be found e.g., in Dufour and Roy (1976). Theorem 3. The stationary Lévy-based stochastic process X can be written in terms of its Fourier coefficients as X(6)= Ao + £ (Ascos(s6)+ Bs sin(s6)) , s=1 6 e [0,2n) , where A0 = u + Ç0Z([0,2n )), As = Çs cos(s^ )Z(d0 ) J0 f 2n Bs = Çs sin(s0)Z(d0) . (11) J0 Moreover, the Fourier coefficients are pairwise uncorrelated and the Fourier coefficients of order s have the same distribution which is characterized by the cumulant generating function KAs(t) = KBs(t) = Kjj (tÇs), where Ku (t)= Kw(t cos(6 )) d6 . J0 Proof. Writing the kernel function in terms of its Fourier representation and then calculate the Fourier coefficients of X gives Eq. 11. The Fourier coefficients are uncorrelated as for all r, s > 1, Cov(As,Br) = K2(W) cos(s^) sin(r0) d^ = 0 , J0 and 1-2% Cov(As,Ar) = k2(W) cos(s0)cos(r0) d0 = 0 , Jo Cov(Bs,Br) = K2(WW sin(s0) sin(r^)d^ = 0 , Jo for all r, s > 1, r = s. The cumulant generating function of As is given by KAs (t)=/ Kz/ (t ^ cos(s^)) d0 Jo 1 2ns s Jo Kz/ (t ^ cos(0 )) d0 = Ku (&t ) and a similar argument shows that KBs (t) = KU (£st). □ The cumulant generating function of As and Bs yield simple expressions for their cumulants, which are given by K1 (Ao) = M , Kr(Ao) = 2n^oKr(W) , r > 2 , and (r- 1)" Kr (As) = Kr (Bs) = (—Kr (W )1(r even) , r!! for s > 1 and r > 1. Here and in the following, we let for a positive integer n, n!! 2■ 4••• n , if neven, 1-3---n , if nuneven. This means that As and Bs have mean, variance, skewness and kurtosis of the following form: K1(As)= 0 Y1(As)= 0 K2(As) = n^s2K2 (W ) 3 Y2(As) = 4nY2(W) , where /2(W) is the kurtosis of W. Moreover, the normalized Fourier coefficients of order s = 1,2,..., obtained by As = As/^s and BBs = will all have the same distribution characterized by the cumulant generating function KU (t). From Theorems 2 and 3, it follows that for a non-Gaussian Lévy basis Z the Fourier coefficients will be uncorrelated with variance Xs (a, p). Furthermore, the distribution of the Fourier coefficients in the non-Gaussian and Gaussian model only differs in even cumulants of order four and higher. Therefore, Eq. 10 can be regarded as a pseudo-likelihood function for (a, p) also in the non-Gaussian case. Fig. 1 shows the small difference in the saddlepoint densities of the normalized Fourier coefficients and the Gaussian density for different values of n in the case of a Gamma Levy basis. Furthermore, a simulation study indicated that the estimates of (a, p) are robust against deviations from Gaussianity in the underlying distribution, but the mean square error of the estimates increases somewhat as n decreases. " \» " V ' \ h h Fig. 1. The saddlepoint density of the normalized Fourier coefficients for n = 2 (red line), n = 4 (green line) and n = 16 (blue line). The density in a Gaussian model is shown for comparison (black line). THE DISTRIBUTION OF Qn - Q UNDER THE LEVY-BASED MODEL In the previous section, we have seen that the method developed in Hobolth et al. (2003) for estimating E(Qn — Q)2 based on a Gaussian process is robust against departures from the distributional assumption. In this section, we will derive the distribution of the error predictor Qn — Q which may be markedly non-Gaussian. Theorem 4. Under the Levy-based model, the error predictor is distributed as ^ f2n Q„ — Q - 2n / )Z(d0) , (12) J0 where ) = £ cos(sn0) . s=1 2 4 6 The distribution of Qn — Q is characterized by its cumulant generating function K; r2n (t) = Kw(2tnkn())d . Jo 1\ (As 0 0 K2(As) = 7Z^2 = Xs(cc,I3) , 9 K(AS) = 2nrj Fig. 3 shows the corresponding saddlepoint densities of the estimated area of the particles for two different values of n. □ Qn-Q Proof. Recall that 271V*1 , 27n\ Qn = — £x(0 + —), (n> 1), n to n where 0 is uniformly distributed in |().171 ¡n). Without loss of generality we can assume that 0 = 0 as the distribution of Qn does not depend on 0. The result given in Eq. 12 is obtained by observing that the mean of the kernel functions is given by 1 2ni °° - y" k(--) = + y" ^nC0s(sn) . " j~0 " s=l The expression for the cumulant generating function of Qn — Q is a consequence of Eq. 7. □ As the cumulant generating function of Qn — Q has a simple form, its cumulants are easily available. In particular, it enables us to obtain a saddlepoint approximation of its density. An alternative is to use Theorem 4 for simulating the distribution of Qn — Q. Example. Let us consider a Levy-based model (Eq. 3) for X with a Gamma Levy basis Z and k chosen such that the covariance function of X follows a p-order model. Under this model the Fourier coefficients /\iV, and Bs, s > 1, of X(0) have mean, variance, skewness and kurtosis, V =2 V 17 = 16 Fig. 2. Realizations of particles obtained by assuming that the squared radial function is given by a Gamma Levy process with a p-order covariance function. The values of p, a and ¡5 were p = 2, log a = 6 and log/3 = — 3. Each column corresponds to realizations for a fixed value ofr\ f i] 2.4. 16/ This model may, for instance, be used to model the squared radial function of random star-shaped planar particles containing the origin. Fig. 2 shows examples of particles simulated from such a model, using different values of rj. The value of p, a and ¡5 was p = 2, log a = 6 and log /3 = —3. We used r\ = 2,4,16; a low value of r\ corresponds to an underlying distribution with a heavy tail. The value of r\ controls the frequency and size of the irregularities of the boundary of the particles. Small values of r\ will produce particles with few large fluctuations on the particle boundary and less smaller fluctuations. Higher values of r\ will produce particles with more frequently occurring moderate fluctuations across the boundary. Fig. 3. The saddlepoint density of the n-point area estimate for n 5 (stippled) and n 10 (full line). The different colours represent densities for the three particles considered: r\ = 2 (red lines), r\ = 4 (green lines) and r\ = 16 (blue lines). Densities for a Gaussian model is shown for comparison (black lines). In applications, it is needed to estimate the parameter n of the underlying Levy basis Z, i.e., the parameter of the distribution of W. For a given kernel function k, we have a simple expression for the cumulant generating function of X and its derivatives (when they exist), cf. Theorem 1. We suggest estimating the parameter n determining the Levy basis by considering the saddlepoint approximation of the density of X (0). For more details on saddlepoint approximations, cf. Jensen (1995). The first order saddlepoint approximation of the density is given by f(x) = 1 V2n KX (t(x)) eKx (t(x))-t(x)x where t (x) is the solution to the saddlepoint equation KX (t)= x . Note that the saddlepoint equation is non-linear when Z is non-Gaussian, but can be solved numerically, using a Newton method for a given kernel function k. A better approximation of the density of X (6) is obtained by multiplying the density with the correction factor c(t(x)) = 1 + KX4)(t(x)) 8KX (?(x))2 5 24 K(x\Kx))\ '(t(x))3/V ' Given an estimation of the kernel function k we can establish a pseudo-likelihood function based on the observations of the stochastic process X given by L(n ) = n f(*(fli)). i=1 Here, f(x(0i)) is calculated using the approximated kernel function S ks(0) = £ Is cos(s0) , s=2 where |s = \JAs(a,p)/n is obtained using the estimates of (a, p). Note that we need to normalize the densities f(xi) for each value of n, when maximizing the likelihood function L(n). If the estimated likelihood function is an increasing function of n , this suggests that the underlying Levy basis is Gaussian. The cumulant generating function and its derivatives have simple analytic expressions, when the underlying Levy basis is a Gamma basis or Inverse Gaussian basis. These expressions can be derived by combining Theorem 1 and Table 1. For a Gamma Levy process the cumulant generating function of X (0) and its derivatives are given by 2n Kx(t)= t(u - ) - n Jo log^1 r kX (t ) = u - 2nWñ+n Jo KXr)(t ) = (r - 1)!n Jo ) o vn - tk($) 2n k(0 )r tk(0 ) d0 , d0 o (vn-tk(0))r d0, r> 2 ' For an inverse Gaussian Lévy process the cumulant generating function of X (0 ) and its derivatives are given by Kx (t ) = t (u - 2) + n4 /0 2n 1 - ^ ) kx (t ) = u - 2n^on2 + n3/ —¡= Jo vn K^r)(t ) = (2r - 3)!!n 3 /q d0, 2 - 2tk(0) k(0 )r -t d^ , r > 2 . 0 (n2 - 2tk(0))r-2 In both cases, the saddlepoint approximation of the density of X (0) are easily obtainable using numerical integration and a Newton algorithm for finding the saddlepoint. Note that the saddlepoint approximation of the density function for an arbitrary n can be written in terms of the saddlepoint approximation of the density, the cumulant generating derivatives and saddlepoint solution for n = 1. CONCLUSION AND PERSPECTIVES We have developed a Levy-based error prediction in circular systematic sampling. In contrast to previous model-based methods, we consider a flexible class of non-Gaussian measurement functions based on kernel smoothing of a homogeneous Levy basis. In particular, we have derived the distribution of the error predictor in circular systematic sampling. The modelling framework allows us to consider in principle any given covariance structure of the measurement function, in particular the popular p-order covariance model which enables controlling the local and global fluctuations of the measurement function. Relation to generalized p-order models. Note that as the Levy-based process X is strictly stationary, it can be shown that X has a polar expansion of the form X(0) = M + V2 £ s/Cs cos(s(0 - Ds)) , s=0 where the random variables 1 o o 2n Cs = 2 (A2 + b2)= AsZs , Ds - U[0,—] , s > 1 , are independent, cf. the Appendix, and E(Zs) = 1. Here, the variable Zs can be expressed as 1 r2n r2n 2nK2(W )Jo Jo cos(s(e - 0))Z(de)Z(d0) . This shows that the Levy-based models are closely related to the generalized p-order models proposed in Hobolth et al. (2003), but the Levy-based models have more structure that allows for derivation of distributional results. Modelling particles in 2D and 3D. The model (Eq. 3) can be used directly to model the shape of featureless two-dimensional particles by assuming that a particle Y is a stochastic deformation of a template particle Y0. If r0 is a radial function of a template particle, we let the radial function of Y be of the form R(0) = r0(d) + X(0), where X is a zero mean Levy-based stochastic process. The strength of this technique is two-folded. Firstly, the global and local fluctuations of the Levy-based stochastic process are controlled by the variance of the Fourier coefficients which are determined by the kernel function. Secondly, the underlying Levy basis determines the frequency and size of the irregularities of the process. Finally, the methodology presented can be extended to model the shape of three-dimensional featureless particles, by considering Levy-based processes on the unit sphere S2. Hansen et al. (2011) consider three-dimensional Levy particles using different covariance models. The focus is here on the Hausdorff dimension of the boundary of particles obtained using a Gaussian basis. Approximations of densities. The saddlepoint approximation was applied here to obtain an approximation of the density of X (0) in the Levy-based stochastic model. As the cumulant generating function is easily obtainable for stochastic convolutions of the type X (Ç ) = - Ç )Z(dn ), Ç e Rn where f : Rn ^ R and Z is a homogeneous Levy basis on Rn, the saddlepoint approximation of the density is an attractive tool for studying Levy-based convolution models in general. For the stochastic processes considered here, other types of approximations of densities can also be considered based on approximating As and Bs by differences of variables from the same family of distributions as W. As an example one could consider approximating the density of the Fourier coefficients by a type II McKay distribution (Holm and Alouini, 2004), when the underlying Levy basis is a Gamma Levy basis. Choice of kernel function. Finally, it should be emphasized that assuming that the kernel function is even does not affect the flexibility of the induced covariance model. When k is not necessarily even, œ k(6) = 1. It is easily seen that œ X(6 + h) = Ao + £(As(h) cos(s6)+ Bs(h) sin(s6)) , s=1 where As(h)\ _ ( cos(sh) sin(shA (As\ ^r(As Bs(h) V- sin(sh) cos(sh) / \B: V Bs and V is a rotation matrix. As for each h e [0,2n), {X(0): d e [0,2n)} and {X(d + h): d e [0,2n)} have the same distribution, As(h)\ _V(As Bs(h) Bs s s We now have that E . /As/v^+B „(As /y^JTB _ VE Es ._{bs/Va2Tb2) ~ V U/VAJTB^J _ VEs ' and consequently Es is uniformly distributed on the unit circle, i.e., arctan(Bs/As) is uniformly distributed on [0,2n) and Ds is uniformly distributed on [0,2n/s). Now consider the conditional distribution of \/2£¡ _ VAJTBJ given Es. As A JÏ£s | Es _ es ~ A JÏ£s | VEs _ es, the conditional distribution does not depend on es and hence y/2£s is independent of Es and Ds. ACKNOWLEDGEMENTS The authors want to express their gratitude to Jan Pedersen for fruitful discussions. This research was supported by MINDLab and Centre for Stochastic Geometry and Advanced Bioimaging, supported by the Villum Foundation. REFERENCES Barndorff-Nielsen OE, Schmiegel J (2004). Lévy based tempo-spatial modeling; with applications to turbulence. Uspekhi Mat Nauk 159:63-90. Cruz-Orive LM (2005). A new stereological principle for test lines in three-dimensional space. J Microsc 219:1828. Dufour JM, Roy R (1976). On spectral estimation for a homogeneous random process on the circle. Stoch Proc Appl 4:107-20. Gual-Arnau X, Cruz-Orive LM (2000). Systematic sampling on the circle and the sphere. Adv Appl Prob SGSA 32:628-47. Gundersen HJG (1988). The nucleator. J Microsc 151:3-21. Hansen LV, Thorarinsdottir TL, Gneiting T (2011). L6vy particles: Modelling and simulating star-shaped random sets. Research report 4, CSGB, Department of Mathematical Sciences, Aarhus University. Submitted. Hellmund G, Prokesovä M, Jensen EBV (2008). L6vy-based Cox point processes. Adv Appl Prob SGSA 40:603-29. Hobolth A, Jensen EBV (2002). A note on design-based versus model-based variance estimation in stereology. Adv Appl Prob SGSA 34:484-90. Hobolth A, Pedersen J, Jensen EBV (2003). A continuous parametric shape model. Ann Inst Statist Math 55:22742. Holm H, Alouini MS (2004). Sum and difference of two squared correlated Nakagami variates in connection with the McKay distribution. IEEE Trans Commun 52:1367-76. Jensen EBV, Gundersen HJG (1993). The rotator. J Microsc 170:35-44. Jensen JL (1995). Saddlepoint Approximations. Oxford: Clarendon Press. Jönsdöttir KY, Hoffmann LM, Hobolth A, Jensen EBV (2006). On variance estimation in circular systematic sampling. J Microsc 222:212-6. Jönsdöttir KY, R0nn-Nielsen A, Mouridsen K, Jensen EBV (2013). L6vy based modelling in brain imaging. Scand J Stat, in press. Jönsdöttir KY, Schmiegel J, Jensen EBV (2008). L6vy-based growth models. Bernoulli 14:62-90. Rajput R, Rosinski J (1989). Spectral representation of infinitely divisible processes. Prob Theory Relat Fields 82:451-87. Schneider R (1993). Convex Bodies: The Brunn-Minkowski Theory. Cambridge: Cambridge University Press.