Analysis of spherical monofractal and multifractal random fields

Analysis of Spherical Monofractal and Multifractal Random Fields with Cosmological Applications Nikolai N. Leonenkoa , Ravindi Nanayakkarab , Andriy Olenko∗,b a School arXiv:2004.14522v1 [math.PR] 30 Apr 2020 b Department of Mathematics, Cardiff University, Senghennydd Road, Cardiff CF24 4AG, UK of Mathematics and Statistics, La Trobe University, Bundoora, VIC 3086, Australia Abstract The Rényi function plays an important role in the analysis of multifractal random fields. For random fields on the sphere, there are three models in the literature where the Rényi function is known explicitly. The theoretical part of the article presents multifractal random fields on the sphere and develops specific models where the Rényi function can be computed explicitly. Then these results are applied to the Cosmic Microwave Background Radiation data collected from the Planck mission. The main statistical model used to describe these data in the literature is isotropic Gaussian fields. We present numerical multifractality studies and methodology based on simulating random fields, computing the Rényi function and the multifractal spectrum for different scenarios and actual CMB data. Key words: Rényi Function, Random Fields, Multifractality, Monofractality, Cosmic Microwave Background Radiation 1. Introduction The concept of multifractality initially emerged in the context of physics. B. Mandelbrot showed the significance of scaling relations in turbulence modelling. Subsequently this concept developed to mathematical models and examining their fine scale characteristics [12]. A multifractal pattern is a type of a fractal pattern that scales with multiple scaling rules in contrast to monofractals that have only scaling rule. A fractal dimension explores the change in characteristics with respect to the change in the scale used. In general, a multifractal scheme is a fractal scheme where its dynamics cannot be explained by a single fractal dimension. More details and references can be found in [14]. Multifractal structures are typical in nature. Multifractal models have been extensively used in the fields of geophysics, genomics, image modelling, finance, hydrodynamic turbulence, meteorology, internet traffic, etc., see references in [3]. Multifractal behaviour has been discovered in stochastic processes as ∗ Corresponding author Email addresses: LeonenkoN@Cardiff.ac.uk (Nikolai N. Leonenko), D.Nanayakkara@latrobe.edu.au (Ravindi Nanayakkara), A.Olenko@latrobe.edu.au (Andriy Olenko) Preprint submitted to Elsevier November 19, 2021 well, see [11, 12, 13]. Multifractal products of stochastic processes have been investigated by [25] with applications of time series in economics and teletraffic. New teletraffic models have been explored and random multifractal measure constructions by considering the stationarity of the processes’ increments were proposed. This methodology has been first introduced in the groundwork [18] on multiplicative chaos and T -martingales of positive type. [16] has shown the multifractal nature of specific Lévy processes and demonstrated that the multifractal spectra of such processes depicts a linear pattern rather than a concave pattern which has been noticed in the actual teletraffic data. [29] has studied the multifractal properties such as scaling exponents of the structure function and the Rényi function of random cascade measures under various conditions. Multifractal analysis has been an important technique in the examination of singular measures and the multifractal spectrum of the random measures based on self-similar processes [8]. The main methods to construct random multifractal structures are based on stochastic processes, branching processes and binomial cascades [34]. The Rényi function plays an important role in the analysis of multifractal random fields. There are several scenarios where the Rényi function was computed for the one-dimensional case and time-series. However, there are very few multidimensional models where it is given in an explicit form. [23] computed the Rényi function for three classes of multifractal random fields on the sphere. It showed some major schemes with regard to the Rényi function which reveal the multifractality of random fields that are homogeneous and isotropic. The Cosmic Microwave Background(CMB) is the radiation from the universe since 380,000 years from the Big Bang. This elongated time period is very short compared to the age of the universe which is of 14 billion years. The CMB is an electromagnetic radiation residue from it’s earliest stage. The CMB depicts variations which corresponds to different regions and represents the roots for all future formation including the solar system, stars and galaxies. At the beginning, the universe was very hot and dense and formation of atoms was impossible. The atoms were split as electrons and protons. That time the universe constituted of a plasma or ionised gas. Then the universe started to expand and cool down. Thus, it had been possible for the atoms to reconcile. This phenomenon is known as “Epoch of combination” and since that time photons have been able to move freely escaping from the opaque of the early universe. The first light which eliminated from this process is termed as the cosmic microwave background, see [37]. In 2009, the European Space Agency launched the mission Planck to study the CMB. The frequency range captured by the Planck is much wider and its resolution is higher than that of the previous space mission WMAP. The CMB’s slight variations were measured with a high precision, see [33]. One of the aims of the mission Planck was to verify the standard model of cosmology using this achieved greater resolution and to find out fluctuations from the specified standard model of cosmology. According to the standard model of the CMB, the Universe is homogeneous and isotropic. This means that almost every part of the universe has very similar properties and that they do not differ based on the direction of the space. However, various research argue that it’s not the case [15, 20, 26, 28, 31, 36]. The motivation of this paper is to develop 2 several multifractal models and the corresponding statistical methodology and use them and other existing models to study whether the Cosmic Microwave Background Radiation data has a multifractal behaviour. The first aim of this paper is to present three known multifractal models for random fields defined on the sphere and suggest several simpler models for which the Rényi function can be explicitly computed. For all these models, we provide plots that illustrate typical multifractal behaviour and their dependence on the scaling parameter. Secondly, we demonstrate the direct probability approach that can be employed to check whether assumptions on models’ parameters guarantee the form of the Rényi function. This approach is less general than the one that is based on martingales for q ∈ [1, 2], see [23, 25]. The advantage is that the proposed methodology is simple and can also be used for q > 2. Finally, we discuss the methodology of computing the Rényi functions and provide various numerical studies of the actual CMB data. The proposed models and methodology can find various applications to other spherical data. The plan of the paper is as follows. Section 2 provides main notations and definitions related to the theory of random fields. Section 3 introduces spherical random fields. Section 4 gives results related to the theory of multifractality and the Rényi function. Section 5 describes the direct probability approach to get the conditions on the limit random measure µ. Section 6 provides results for three known models of the Rényi function for spherical random fields of the exponential type. Section 7 proposes new models based on power transformations of Gaussian fields. Section 8 presents numerical studies including computing and fitting the empirical Rényi functions for CMB data from different sky windows and models. Finally, the conclusions and some new problems are given in Section 9. All numerical studies were conducted by using Maple 2019.0 and R 3.6.3 software, in particular, the R packages ‘rcosmo’ [9, 10] and ‘RandomFields’ [35]. A reproducible version of the code in this paper is available in the folder “Research materials” from the website https://sites.google.com/site/olenkoandriy/. 2. Main notations and definitions This section presents background materials in the random fields theory and multifractal analysis methodology. Most of the material included in this and next two sections are based on [12, 21, 22, 24, 25, 27]. Let S ⊂ Rn be a multidimensional set, k · k denotes the Euclidean distance in Rn and sn−1 (1) = {u ∈ d Rn : kuk = 1}. The notation | · | is used for the Lebesgue measure on Rn . {·} = {·} will stand for the equality of finite dimensional distributions. Definition 2.1 A random field is a function ξ(ω, x) : Ω × S → Rm such that ξ(ω, x) is a random vector for each x ∈ S. For simplicity it will also be denoted by ξ(x), x ∈ S. When n = 1 ξ(x) is a random process. When S ⊆ Rn , n > 1, then ξ(x) is termed as a random field. It is called a vector random field for m > 1. In this paper, we mainly concentrate on scalar random fields ξ(x), x ∈ S, n > 1 and m = 1. 3 If {ξ(x1 ), ..., ξ(xN ), x1 , ..., xN ∈ S} is a set of random variables belonging to a Gaussian system for each N ≥ 1, then ξ(x), x ∈ S, is called Gaussian. We assume that all random variables ξ(x) are defined on the same probability space (Ω, A, P ). Random fields are indexed by elements of subsets of Rn . Definition 2.2 A second order random field is a random function ξ : S → L2 (Ω, A, P ), S ⊂ Rn . In other words, the random variable ξ(x), x ∈ S, satisfies E|ξ(x)|2 < +∞. Thus, a second order random field over S is a family {ξ(x), x ∈ S} of square integrable random variables. Definition 2.3 A real valued random field ξ(x), x ∈ Rn , satisfying E[ξ 2 (x)] < ∞ is homogeneous (in the wide sense) if its mathematical expectation m(x) = E[ξ(x)] and covariance function B(x, y) = cov(ξ(x), ξ(y)) are invariant with respect to the Abelian group G = (Rn , +) of shifts in Rn , that is m(x) = m(x + τ ), B(x, y) = B(x + τ, y + τ ), for any x, y, τ ∈ Rn . That is, for homogeneous random fields E[ξ(x)] = const, and the covariance function B(x, y) = B(x − y) depends only on the difference x − y. It will be also assumed that E[ξ(x)] = 0 without loss of generality. n The covariance function B(x − y) of a homogeneous random field is a non-negative definite kernel on R × Rn , that is, for any r ≥ 1, x(j) ∈ Rn , zj ∈ C, j = 1, ..., r, r X i,j=1 B(x(i) − x(j) )zi z¯j ≥ 0. If the covariance function B(x) is continuous at x = 0, then the field is mean-square continuous for each x ∈ Rn and vice versa. Definition 2.4 The random field ξ(x) is isotropic on Rn , if Eξ(x)2 < ∞, and its first and second order moments are invariant with respect to the group of rotations on the sphere, i.e. Eξ(x) = Eξ(gx), Eξ(x)ξ(y) = Eξ(gx)ξ(gy), for every g ∈ SO(n). Definition 2.5 A real valued random field ξ(x), x ∈ Rn , satisfying E[ξ 2 (x)] < ∞ is homogeneous and isotropic if its mathematical expectation m(x) = E[ξ(x)] = const, and the covariance function B(x, y) = cov(ξ(x), ξ(y)) = B(kx − yk) depends only on the Euclidean distance ρxy = kx − yk between x and y. It means that its mathematical expectation m(x) and covariance function B(x, y) are invariant with respect to shifts, rotations and reflections in Rn : m(x) = m(x + τ ), m(x) = m(gx), B(x, y) = B(x + τ, y + τ ), B(x, y) = B(gx, gy), for any x, y ∈ Rn , g ∈ SO(n), where SO(n) is the group of rotations in Rn . 4 Definition 2.6 A stochastic process {X(t), t ≥ 0} is self-similar if for any non-random a > 0, there exists d non-random b > 0 such that {X(at)} = {bX(t)}. For self-similar, continuous at 0 and non-trivial X(t), the constant b must be equal aH , a > 0, where d H ≥ 0. Thus, {X(at)} = {aH X(t)}. The constant H is known as the Hurst parameter. We will call the process {X(t), t ≥ 0} H − ss(self-similar) or H − sssi(self-similar stationary increments) if its increments are stationary. The concept of multifractal processes was motivated by establishing the following scaling rule of selfsimilar processes. d Definition 2.7 A stochastic process X(t) is multifractal if {X(ct)} = {M (c)X(t)}, where M (c) is a random variable independent of X(t) for every c > 0 and the distribution of M (c) does not depend on t. The process is self-similar if M (c) is non-random for every c > 0 and M (c) = cH . The scaling factor d M (c) satisfies {M (ab)} = {M1 (a)M2 (b)} for every selection of constants a and b and random M1 and M2 that are independent copies of M . This establishes the characteristic of the deterministic factor H − ss processes (ab)H = aH bH . Another definition of multifractality is Definition 2.8 A stochastic process X(t) is multifractal if there exist non-random functions c(q) and τ (q) such that for all t, s ∈ T , q ∈ Q, E|X(t) − X(s)|q = c(q)|t − s|τ (q) , where T and Q are intervals on the real line with positive length and 0 ∈ T . The function τ (q) is known as the scaling function. The interval Q may include negative values. Instead of the increments of the process, the definition can also be established on the moments of the process. i.e. E|X(t)|q = c(q)tτ (q) . Above definitions coincide if the increments are stationary. If {X(t)} is H − sssi, then it holds that τ (q) = Hq. 3. Spherical random fields This section introduces some basic notations of the theory of random fields on a sphere. The sphere is a simplest case of a manifold in Rn . For simplicity, we consider only the case n = 3. Let us denote the 3-dimensional unit ball as B 3 = {x ∈ R3 : kxk ≤ 1}. The spherical surface in R3 with a given radius r > 0 is s2 (r) = {x ∈ R3 : kxk = r}, with the corresponding Lebesgue measure on the sphere σr (du) = σr (dθ · dϕ) = r2 sin θdθdϕ, (θ, ϕ) ∈ s2 (1). The Kronecker delta is a function defined as: δij =   0,  1, if i 6= j, if i = j. 5 For ν > − 21 , we use the Bessel function of the first kind of order ν Jν (z) = ∞ X (−1)m m=0  z 2m+ν 2 [m!Γ(m + ν + 1)]−1 , z > 0. A spherical random field T = {T (r, θ, ϕ) : 0 ≤ θ ≤ π, 0 ≤ ϕ ≤ 2π, r > 0} is a random function, which is defined on the sphere s2 (r). We deal with a spherical real-valued mean-square continuous random field T with the constant mean M and finite second order moments. In this case the covariance function B(cos θ) is a continuous function on [0, π). Definition 3.1 A real valued random field T (x), x ∈ s2 (1), with E[T (x)] < ∞ and E[T (x)] = 0 is isotropic if E[T (x1 )T (x2 )] = B(cos θ), x1 , x2 ∈ s2 (1), depends only on the geodesic (angular) distance cos θ between x1 and x2 . An isotropic spherical random field on s2 (r) can be expanded in a Laplace series in the mean-square sense. T (r, θ, ϕ) = M + l ∞ X X Ylm (θ, ϕ)am l (r), (1) l=0 m=−l where {Ylm (θ, ϕ)} represents the spherical harmonics defined as m Ylm (θ, ϕ) = cm l exp (imϕ)Pl (cos θ), l = 0, 1, ..., m = 0, ±1, ..., ±l, with cm l = (−1) m  2l + 1 (l − m)! 4π (l + m)! 1/2 , and the Legendre polynomials Plm (cos θ) having degree l and order m, see [23]. The notation T̃ (x) = T (r, θ, ϕ), x ∈ R3 , will be used to highlight the random field’s dependence on Euclidean coordinates. The random coefficients of the Laplace series can be computed as the mean-square stochastic integrals via the inversion arguments as Z π Z 2π Z m (θ, ϕ)r 2 sin θdθdϕ = (r) = T (r, θ, ϕ)Y am l l 0 0 where u = x kxk ∈ s2 (1), r = kxk. ′ s2 (1) T̃ (ru)Ylm (u)σ1 (du), (2) ′ The covariance functions E(T (r, θ, ϕ)T (r, θ , ϕ )) of the isotropic random fields depend only on the an′ ′ gular distance θ = θP Q between the points P = (r, θ, ϕ) and Q = (r, θ , ϕ ). For spherical isotropic random fields it possesses ′ ′ ′ l m m Eam l (r)al′ (r) = δl δm Cl (r), 2 E|am l (r)| = Cl (r), m = 0, ±1, ..., ±l. (3) The angular power spectrum of the isotropic random field T (r, θ, ϕ) is defined as the functional series {C1 (r), C2 (r), ..., Cl (r), ...}. 6 From (1), (2) and (3) we obtain ′ ′ Cov(T (r, θ, ϕ), T (r, θ , ϕ )) = ∞ 1 X (2l + 1)Cl (r)Pl (cos θ), 4π l=1 where P∞ l=1 (2l + 1)Cl (r) < ∞ for every r > 0. If T (r, θ, ϕ) is an isotropic Gaussian field defined on the sphere s2 (r), then the coefficients am l (r) are independent Gaussian random variables that are complex-valued with ′ ′ ′ m l m Eam l (r)al′ (r) = δm δl Cl (r). Eam l (r) = 0, The homogeneous and isotropic random field T (x), x ∈ R3 , which is restricted to the sphere s2 (r) has ′ ′ ′ l m m Eam l (r)al′ (s) = δl δm Cl (r, s), where Cl (r, s) = 2π 2 Z ∞ Jl+ 12 (µr)Jl+ 12 (µs) (µr)1/2 (µs)1/2 0 G(dµ), l = 1, 2, ..., and G(·) is a finite measure defined on the Borel sets of [0, ∞) satisfying Z ∞ 2 σ = V ar{T̃ (0)} = G(dµ) < ∞. 0 4. Rényi function and multifractal spectrum This section introduces basic notations and concepts of the multifractal theory and Rényi functions. Consider a random field Λ(x, ω), x ∈ R3 , ω ∈ Ω, that is measurable, homogeneous and isotropic (HIRF) on the 3-dimensional Euclidean space R3 . It will be called the mother field. For simplicity it will be denoted as Λ(x) = Λ(x, ω). Condition 1: Let a random field Λ(x), x ∈ R3 , satisfy Λ(x) = 1, Λ(x) > 0, 2 ρΛ (kx − yk), Cov(Λ(x), Λ(y)) = RΛ (kx − yk) = σΛ 2 where ρ(0) = 1 and σΛ < ∞. Let Λ(i) (x), x ∈ R3 , i = 0, 1, 2, ..., be a sequence of independent copies of the random field Λ(·). We consider the re-scaling of Λ(·) defined as Λ(i) (bi x), where b > 1 is a constant called a scaling factor and bi x is the product of a vector x by a scalar bi . A finite-product field on B 3 is defined by Λk (x) = k Y Λ(i) (bi x), k = 1, 2, .... i=0 Then one can introduce the random measure µk (·) on the Borel σ-algebra B of a unit ball B 3 by Z Λk (y)dy, A ∈ B, k = 0, 1, 2, .... µk (A) = y∈A 7 d We denote by µk − → µ, k → ∞, the weak convergence of the measures µk to some measure µ. It means that for all continuous functions g(y), y ∈ B 3 , it holds Z Z g(y)µk (dy) → g(y)µ(dy), k → ∞. B3 B3 Definition 4.1 The Rényi function of a random measure µ is a non-random function defined by P (m) log2 E l µ(Bl )q , T (q) = lim inf (m) m→∞ log2 |Bl | (m) where {Bl , l = 0, 1, ..., 2m − 1, m = 1, 2, ..., } denotes the mesh formed by the mth level dyadic decomposi- tion of the unit ball B 3 . The main result on the form of the Rényi function is given by the following theorem. Theorem 4.1 [23] Suppose that Condition 1 holds. (i) Assume that the correlation function ρΛ (kx−yk) = ρ(r) of the field Λ(·) satisfies the following condition |ρΛ (r)| ≤ Ce−γr , r > 0, (4) for some positive constants C and γ. Then, for the scaling factor b > k → ∞, on B 3 . p d 3 2 , the measures µ − 1 + σΛ k → µ, (ii) If for some range q ∈ Q = [q− , q+ ], both E q Λ(0) < ∞ and Eµq (B 3 ) < ∞, then the Rényi function T (q) of µ is given by T (q) = q − 1 − 1 logb EΛq (x), 3 q ∈ Q. Similarly, for spherical random fields on s2 (1), one can introduce an analogous approach. Condition 2: Let the random field Λ̃(x), x ∈ s2 (1), satisfy 2 V arΛ̃(x) = σΛ < ∞, E Λ̃(x) = 1, ′ ′ Cov(Λ(θ, ϕ), Λ(θ , ϕ )) = Λ(x) > 0, ∞ 1 X (2l + 1)Cl Pl (cos θ), 4π l=1 ∞ X l=1 (2l + 1)Cl < ∞. Let Λ̃(i) (x), x ∈ s2 (1), i = 0, 1, 2, ..., be a sequence of independent copies of the field Λ̃(·). Consider Λ̃(i) (bi x), where b > 1 is a scaling factor, bi x := (1, bi × θ, bi × ϕ) ∈ s2 (1), and the modulus algebra is used π 2π to compute the products. Define the finite product fields on s2 (1) by Λ̃k (x) = k Y Λ̃(i) (bi x), k = 1, 2, .... i=0 Let us introduce the random measure µk (·) on the Borel σ-algebra B of s2 (1) as Z Λk (y)dy, k = 0, 1, 2, ..., A ∈ B. µk (A) = A 8 (5) d We denote by µk − → µ, k → ∞, the weak convergence of the measures µk to some non-degenerate measure µ. It means that for all continuous functions g(y), y ∈ s2 (1), Z Z g(y)µ(dy), g(y)µk (dy) → k → ∞. s2 (1) s2 (1) The Rényi function of the random measure µ defined on s2 (1) is determined as T (q) = lim inf log2 E (m) l (m) q µ(Sl ) (m) log2 |Sl | m→∞ where {Sl P , (6) , l = 0, 1, ..., 2m − 1} is the mesh constructed by mth level dyadic decomposition of the spherical surface of s2 (1). Theorem 4.2 [23] Suppose that Condition 2 holds and the isotropic random field Λ̃(·) is the restriction to the sphere s2 (1) of the HIRF Λ(x), x ∈ R3 , with the correlation function ρΛ (kx − yk) = ρ(r). Under the same assumptions as in Theorem 4.1, the Rényi function T (q) of the limit measure µ on s2 (1) is given by T (q) = q − 1 − 1 logb EΛq (t), 2 q ∈ Q. The multifractal or singularity spectrum is defined via the Legendre transform as f (h) = inf (hq − T (q)). q (7) and is used to describe local fractal dimensions of random fields. 5. Conditions on measure µ The random measure µ in the previous section was defined as a weak limit of the measures µk . Therefore, it would be difficult to check moment conditions on µ as its probability distribution is not explicitly known. 2 This section demonstrates how to prove sufficient conditions on the scaling factor b and the variance σΛ that guarantee Eµq (B 3 ) < ∞. The general method to obtain such conditions for the range q ∈ [1, 2] uses martingale L2 convergence, see, for example, [25]. In this section we use the direct probability approach, which is more elementary. From the weak convergence of the measures µk to µ, it follows that Eµqk (B 3 ) → Eµq (B 3 ), k → ∞. By the Lyapunov’s inequality Eµqk (B 3 ) ≤ (Eµ2k (B 3 ))q/2 , for q ∈ [1, 2]. 2 Therefore, to guarantee Eµq (B 3 ) < +∞, q ∈ [1, 2], it is sufficient to provide such b and σΛ that sup Eµ2k (B 3 ) < +∞. k∈N 9 By (4), the non-negativity of Λk (y) and independence of Λ(i) it holds Z Z Z Z 2 3 E[Λk (y)Λk (ỹ)]dỹdy Λk (y)Λk (ỹ)dỹdy = Eµk (B ) = E B3 Z = = = Z Z B3 B3 ≤ Z Z Z B3 Z E B3 Z Λ(i) (ybi )Λ(i) (ỹbi )dỹdy = i=0 k  Y B 3 i=0 k Y k Y B3 B3 i (i) B3 B3 i EΛ(i) (ybi )Λ(i) (ỹbi )dỹdy i=0 (i) (Cov(Λ(ybi ), Λ(ỹbi )) + 1)dỹdy = k Y B3 k Y  E(Λ (yb ) − 1)(Λ (ỹb ) − 1) + EΛ (yb ) + EΛ (ỹb ) − 1 dỹdy (i) B 3 i=0 Z Z B3 (1 + B 3 i=0 i 2 σΛ Ce−γky−ỹkb )dỹdy Z ≤ Z Z B3 B3 Z k Y ∞ Y B 3 i=0 ≤ Z B3 Z ∞ Y 2 eσΛ Ce (i) i
2 1 + σΛ ρΛ (ky − ỹkbi ) dỹdy B 3 i=0 From the inequality 1 + a ≤ ea , it follows that Eµ2k (B 3 ) i i 2 (1 + σΛ Ce−γky−ỹkb )dỹdy. −γky−ỹkbi dỹdy. B 3 i=0 Introducing the new variables z = y, z̃ = y − ỹ, one obtains Eµ2k (B 3 ) ≤ Z dz B3 Z ∞ Y 2 eσΛ Ce −γkz̃kbi dz̃, B 3 −B 3 i=0 where B 3 − B 3 = {z̃ : z̃ = y − ỹ, y, ỹ ∈ B 3 }. Hence, by using the spherical change of variables, Eµ2k (B 3 ) 2 As eσΛ Ce −rbi 3 ≤ |B | Z diam(B 3 ) r ∞ Y 2 0 i i=0 |B 3 | dr = 3 . γ Z γdiam(B 3 ) r2 0 ∞ Y 2 eσΛ Ce −rbi dτ. i=0 is a decreasing function of r, for n(r) = max(0, −[logb (r)]), r > 0, we obtain ∞ Y e 2 Ce−rb σΛ n(r)−1 i ≤ i=0 ≤e Notice that e 2 σΛ Ce−γrb ∞ Y 2 Cn(r) σΛ ∞ Y e Y e 2 Ce−rb σΛ i=0 2 Ce−rb σΛ 2 −bi 2 = e σΛ C i+n(r) P∞ ≤e i=0 e−b i=0 =e 2C σΛ e P∞ i=0 e−(b−1) i 2 eσΛ Ce −rbi i=n(r) i=0 eσΛ Ce ∞ Y i =e 10 i 2 Cn(r) σΛ ∞ Y 2 eσΛ Ce −bi . i=0 2 ≤ e σΛ C P∞ i=0 2C σΛ 1 e 1−e−(b−1) e−(1+(b−1) < +∞. i) Therefore, 2 Eµ2k (B 3 ) σΛ C |B 3 | ≤ 3 e e(1−e−(b−1) ) γ 3 |B | = 3 e γ The integral is finite if 2 − 2 σΛ C ln (b) Z 2C σΛ e(1−e−(b−1) ) > −1, i.e. b > e Z γdiam(B 3 ) 2 z 2 eσΛ Cn(z) dz 0 γdiam(B 3 ) σ2  2 z max 1, z 0 2C σΛ 3 Λ − ln (b)  dz. . Thus, we proved the following result. Theorem 5.1 Let the mother field Λ(x) > 0, x ∈ R3 , satisfy the conditions EΛ(x) = 1, 2 V arΛ(x) = σΛ < +∞, 2 Cov(Λ(x), Λ(y)) = σΛ ρΛ (kx − yk), |ρΛ (τ )| ≤ Ce−γτ , τ > 0, p σ2 C 2 , e Λ3 ). and the scaling factor b > max( 3 1 + σΛ d Then the measures µk − → µ, k → ∞, and Eµq (B 3 ) < +∞, for q ∈ [1, 2]. Remark 5.1 The approach developed in this section can be used to obtain conditions on the mother field that guarantee Eµqk (B 3 ) < +∞, for q in the range [1, Q], where Q > 2. For example, using the Lyapunov’s inequality for q ∈ [1, 4] Eµqk (B 3 ) ≤ (Eµ4k (B 3 ))q/4 , the conditions on b and σλ2 that guarantee Eµ4k (B 3 ) < +∞ are also sufficient for Eµq (B 3 ) < +∞, q ∈ [1, 4]. Then, it follows from Eµ4k (B 3 ) = Z B3 Z B3 Z B3 Z k Y B 3 i=0  E 4 Y j=1  Λ(i) (yi bi ) that one needs some assumptions on the fourth order moments E( Q4 j=1 4 Y dyi , j=1 Λ(yj bi )) or cumulants of the mother field Λ(·). We will provide an example of such conditions in Section 7 for Model 4. 6. Rényi functions of exponential models For the random fields on the sphere, there are three models where the Rényi function is known explicitly, see [23]. These models were obtained for exponential type spherical random fields. Model 1 Let the random field Λ(x) be given as   1 2 Λ(x) = exp Y (x) − σY , 2 where Y (x), x ∈ R3 , is a zero-mean Gaussian, measurable, separable random field with the covariance function σY2 ρY (r), ρY (0) = 1. The following result provides the conditions and the explicit form of the Rényi function for Model 1. 11 Theorem 6.1 [23] Let for Model 1 the correlation function satisfy 0 < |ρY (r)| ≤ Ce−γr , r > 0, for some positive C and γ and b > exp{ 2 σY 3 }. If Y (x), x ∈ s2 (1), is a spherical isotropic random field that is a restriction of Y (x), x ∈ R3 , on the sphere  s2 (1), then the random measures (5) generated by the spherical fields Λ(x) = exp Y (x) − 12 σY2 , x ∈ s2 (1), converge weakly to the random measure µ. The corresponding Rényi function is  2    σY σ2 − q2 − 1, q ∈ [1, 2]. T (q) = q 1 + Y 4 ln b 4 ln b (8) Model 2 Let the random field Λ(x) be of the form  1 cZ = − ln 1 − λ Λ(x) = exp {Z(x) − cZ } , β , where Z(x), x ∈ R3 , is a gamma-correlated HIRF with the correlation function ρZ (r). The field Z(x) has the marginal density p(u) = λβ β−1 −λu u e , Γ(β) (9) u, λ, β ∈ (0, +∞), and the bivariate density β−1    √  u1 · u2 · α (u1 u2 /α) 2 u1 + u2 p0 (u1 , u2 ; α) = Iβ−1 2 , exp − Γ(β)(1 − α) 1−α 1−α (10) where α ∈ [0, 1], λ, β, and γ are constant parameters, Iv (·) is the modified Bessel function of the first
2k+v P∞ −1 (k!Γ(k + v + 1)) . kind given by Iv (z) = k=0 z2 Then the covariance function of the mother random field is ! e−2cz ρΛ (r) =
β − 1 1 − λ2 + λ22 (1 − ρZ (τ )) e−2cz
β − 1 1 − λ1 !−1 . The following result gives the Rényi function and the corresponding conditions for Model 2. Theorem 6.2 [23] Suppose that for Model 2 the parameter λ > 2 and the correlation function satisfies 0 < |ρZ (r)| ≤ Ce−γr , r > 0, for some positive constants C and γ. Then, for the parameters (β, λ) from the set Lβ,λ  = {(β, λ) : b > 1 + 1 λ2 1− 2 λ  β2 , λ > 2, β > 0}, d the measures µk − → µ, k → ∞. The Rényi function of µ is given by       β q 1 β + logb 1 − − 1. T (q) = q 1 − logb 1 − 2 λ 2 λ where q ∈ Q = {0 < q < λ} ∩ [1, 2] ∩ Lβ,λ . 12 (11) Model 3 Let the mother random field be Λ(x) = exp {U (x) − cU } , x ∈ R3 , where U (x) = −Z −1 (x), and Z(x), x ∈ R3 , is a gamma-correlated HIRF with the densities given by (9) and (10) and the correlation function ρZ (r). Theorem 6.3 [23] Suppose that for Model 3 the correlation function satisfies 0 < |ρZ (r)| ≤ Ce−γr , r > 0, for some positive constants C and γ. Then, for any (β, λ) ∈ Lβ,λ and b > d  1 2 β √ Γ(β)2 2 −1 Kβ (2 2λ) √ λβ/2 [Kβ (2 λ)]2 the measures µk − → µ when k → ∞. The Rényi function of measure µ is   β/2    p  cU  1 1 2λ T (q) = q 1 + − logb q β/2 Kβ (2 qλ) − 1 + logb , 2 ln b 2 2 Γ(β) (12) where q ∈ Q = [1, 2] ∩ Lβ,λ , Kλ (x) is the modified Bessel function of the third kind and √ ! 2λβ/2 Kβ (2 λ) . cU = ln Γ(β) Let α(q) denote the qth order singularity exponent and be defined by α(q) = d T (q). dq (13) Then the multifractal spectrum defined by (7) can be expressed as a function of α by f (α(q)) = q · α(q) − T (q). (14) For Model 1 it is easy to see from (8) that α(q) = 1 + f (α(q)) = 1 − σY2 σY2 − q, 4 ln(b) 2 ln(b) σY2 q2 , 4 ln(b) q ∈ [1, 2]. (15) By (11) we obtain for Model 2 α(q) = 1 − f (α(q)) = 1 +   β β 1 + logb 1 − , 2 λ 2 ln(b)(q − λ) β 2 ′    q q − logb 1 − . ln(b)(q − λ) λ For Model 3 it follows from (12) and Kβ (q) = − 12 (Kβ−1 (q) + Kβ+1 (q)), see 9.6.26 in [1], that √ √ √ cU λ(Kβ−1 (2 qλ) + Kβ+1 (2 qλ)) β √ √ α(q) = 1 + , − + 2 ln(b) 4 ln(b)q 2 ln(b)Kβ (2 qλ) q 13 (16)  β/2  p β β 1 2λ − logb + logb (q β/2 Kβ (2 qλ)) 2 Γ(β) 4 ln(b) 2 √ √ √ qλ(Kβ−1 (2 qλ) + Kβ+1 (2 qλ)) √ + . 2 ln(b)Kβ (2 qλ) f (α(q)) = 1 + (17) Summarising the above results we obtain Theorem 6.4 Let the corresponding conditions of Theorems 6.1, 6.2 and 6.3 are satisfied for Models 1, 2 and 3. Then the multifractal spectra of these models are given by (15), (16) and (17) respectively. The plots shown in Figure 1 illustrate behaviours of the Rényi functions and multifractal spectra for Models 1, 2 and 3. For these numerical examples, we used the following values of the parameters: b = 2, σY = 1, λ = 3, and β = 2. Notice that these values of b, λ and β satisfy the conditions in Lβ,λ . We also selected (0, 3) as the range of q values. It is slightly wider than the range [1, 2] in the theorems and allows better visualisation of T (q) and f (α), see Section 8.1 on the way to check its validity. a) Rényi function of Model 1 b) Rényi function of Model 2 c) Rényi function of Model 3 d) Spectrum of Model 1 e) Spectrum of Model 2 f) Spectrum of Model 3 Figure 1: Examples of Rényi functions and multifractal spectra for Models 1, 2 and 3. Figure 1 shows that the Rényi functions of the Models 1 and 2 have parabolic shapes while the Rényi 14 function of the Model 3 is closer to a linear shape on the interval (0, 3). Also, comparing the plots for Models 1, 2 and 3, we can see that the Rényi functions of Model 1 and 2 exhibit a concave down increasing and decreasing behaviour within q ∈ (0, 3), whereas for Model 3 it increases. The multifractal spectra of Models 1, 2 and 3 show a concave down increasing behaviour within q ∈ (0, 3). 7. Models based on power transformations of Gaussian fields In the previous section, we considered three models based on an exponential transformation of Gaussian or gamma-correlated HIRF. This section propose few much simpler scenarios where conditions of the theorems from Section 4 are satisfied. First, note that the condition Λ(x) > 0 in [23] can be relaxed to Λ(x) > 0 almost sure. Model 4 Let Λ(x) = Y 2 (x), where Y (x), x ∈ R3 , is a zero-mean unit variance Gaussian HIRF with a covariance function ρY (τ ), τ ≥ 0. Then it follows that 2 σΛ = V arΛ(x) = E(Y 4 (x)) − 1 = 2, EΛ(x) = E(Y 2 (x)) = ρY (0) = 1, Cov(Λ(x), Λ(y)) = E(Y 2 (x) − 1)(Y 2 (y) − 1) = 2ρ2Y (kx − yk). To compute the covariance, we used the property E(Hk (Y (x))Hl (Y (y))) = δkl k!ρkY (kx − yk), x, y ∈ R3 , (18) where Hk (u), k ≥ 0, u ∈ R, are the Hermite polynomials, see [32]. For k = 2, the Hermite polynomial of order 2 is H2 (u) = u2 − 1. Thus, Model 4 satisfies Conditions 1 and 2. Note that the condition |ρΛ (r)| ≤ Ce−γr , r > 0, γ > 0, is equivalent to ′ ′ ′ |ρY (r)| ≤ C e−γ r , r > 0, γ > 0. (19) So, if (19) is satisfied, then one can apply Theorems 4.1 and 4.2 and the Rényi function of the limit measure equals to T (q) = q − 1 − 1 logb EY 2q (x). 2 Note that for p > −1 and Z ∼ N (µ, σ 2 ) E|Z − µ|p = σ p Therefore, we obtain the following result. 15 2p/2 Γ( p+1 ) √ 2 . π (20) Theorem 7.1 Suppose that for Model 4, the correlation function of Y (x) satisfies |ρY (r)| ≤ Ce−γr , r > p 2 0, γ > 0, and b > max( 3 1 + σλ2 , eσΛ C/3 ). d → µ, k → ∞, and the corresponding Rényi function is Then the measures µk − 1 T (q) = q − 1 − logb 2  2q Γ(q + 12 ) √ π  , q ∈ [1, 2]. (21) Example Now we demonstrate how the approach developed in Section 5 can be used to obtain the moment conditions for q ∈ [1, 4] in the case of Model 4. By Remark 5.1, it is enough to check that   Z Z Z Z Y 4 4 k Y Y Y 2 (yj bi ) sup Eµ4k (B 3 ) = sup E dyj < +∞. k∈N k∈N B3 B3 B3 B 3 i=0 Notice, that by Wick’s theorem   4 Y X E Y 2 (yj bi ) = j=1 Y j=1 j=1 Cov(Y (yj bi ), Y (yj̃ bi )), (22) p∈P42 (j,j̃)∈p where the sum is over all parings p of {1, 1, 2, 2, 3, 3, 4, 4}, which are distinct ways of partitioning {1, 1, 2, 2, 3, 3, 4, 4} into pairs (i, j). The product in (22) is over all pairs contained in p, see [17]. Notice that for the pairing p∗ = {(1, 1), (2, 2), (3, 3), (4, 4)}. Y Cov(Y (yj bi ), Y (yj̃ bi )) = 4 Y EY 2 (yj bi ) = 1. j=1 (j,j̃)∈p∗ In all other cases of pairing there is at least one pair (j, j̃) such that j 6= j̃. Therefore, the expectation Q4 E( j=1 Y 2 (yj bi )) equals X Y 1+ Cov(Y (yj bi ), Y (yj̃ bi )). p∈P42 (j,j̃)∈p p6=p∗ As, 1 + a < ea , it can be estimated by    X Y   Cov(Y (yj bi ), Y (yj̃ bi )) . exp    2 p∈P4 (j,j̃)∈p p6=p∗ As at least for one pairing (j, j̃) ∈ p 6= p∗ it holds that j 6= j̃, then one can use the upper bound i |Cov(Y (yj bi ), Y (yj̃ bi ))| ≤ σY2 Ce−γkyj −yj̃ kb , and the approach from Section 5. 16 Namely, 4 Y E 2 i Y (yj b ) j=1 ! ≤e P p∈p2 4 Q (j,j̃)∈p 2 Ce σY −γkyj −y kbi j̃ 2 ≤ e(max(σΛ C,1)) 4 P 1≤j≤j̃≤4 e −γkyj −y kbi j̃ . Hence, sup k∈N ≤  ≤ Z B3 Z Z B3 B3 Z Z Z B3 B3 B3 Z Z Z B3 B3 Z Z B3 Z k Y B 3 i=0 ∞ Y  E 4 Y j=1 2 e(max(σΛ C,1))  Y 2 (yj bi ) 4P 1≤j≤j̃≤4 4 Y dyj j=1 −γkyj −y kbi j̃ e B 3 i=0 ∞ Y e P 2 6(max(σΛ C,1))4 1≤j≤j̃≤4 −γkyj −y kbi j̃ e B 3 i=0 4 Y j=1 where the last inequality follows from the generalized Hölder’s inequality K Y k=1 where PK k=1 fk 1 ≤ K Y k=1  dyj  , kfk kpk , pk −1 = 1. In our case K = 6 is the number of different j and j̃ satisfying 1 ≤ j ≤ j̃ ≤ 4. Finally, similar to the proof in Section 5, from (21) we obtain the condition b>e σ(max(σΛ C,1))4 3 . Now we show that the assumption Λ(x) > 0 almost surely is indeed not restrictive and it is easy to construct a modification of Model 4 with Λ(x) > 0. Model 4’ Let Λ(x) = Y 2 (x) · (1 − ε) + ε, ε ∈ (0, 1), where Y (x), x ∈ Rn , is a zero-mean unit variance Gaussian HIRF with a covariance function ρY (τ ), τ ≥ 0. It is easy to see that EΛ(x) = (1 − ε)EY 2 (x) + ǫ = 1, 2 σΛ = V arΛ(x) = 2(1 − ε)2 < +∞, Cov(Λ(x), Λ(y)) = 2(1 − ε)2 ρ2Y (kx − yk). ′ Hence, Model 4’ satisfies Conditions 1 and 2 and |ρ(r)| ≤ Ce−γr , r > 0, γ > 0, if |ρY (r)| ≤ C ′ e−γ r , r > ′ 0, γ > 0. Therefore, we obtain 1 logb E((1 − ε)Y 2 (x) + ε)q 2 q  1 ε q . = q − 1 − logb (1 − ε) − logb E Y 2 (x) + 2 2 1−ε Tε (q) = q − 1 − and 1 lim Tε (q) = q − 1 − logb E(Y 2q (x)), 2 17 ε→0 which coincides with (21). The next model generalizes Model 4 to an arbitrary even power of a Gaussian random field. Model 5 Let Λ(x) = Y 2k (x), where Y (x), x ∈ R3 , is a zero-mean Gaussian HIRF with the variance  √ − k1 π σ 2 = 2k Γ(k+ and a covariance function ρY (r), r ≥ 0. 1 ) 2 Then it follows from (20) that EΛ(x) = EY 2k (x) = σ 2k 2 σΛ = V arΛ(x) = E(Y 4k ) − 1 =  √ π k 2 Γ(k + 12 ) 2 2k Γ(k + 12 ) √ = 1, π 22k Γ(2k + 12 ) −1= π √ πΓ(2k + 12 ) − 1 < +∞. Γ2 (k + 12 ) To compute the covariance function we use (18) and the following Hermite expansion [1, page 775] z 2k = (2k)! k X H2k−2i (z) . i i!(2k − 2i)! 2 i=0 Therefore, Cov(Λ(x), Λ(y)) = E(Y 2k (x) − 1)(Y 2k (y) − 1) = π 22k Γ2 (k + 12 ) E(Ỹ 2k (x)Ỹ 2k (y)) − 1 k = ((2k)!)2 X E[H2k−2i (Ỹ (x))H2k−2i (Ỹ (y))] π −1 1 2k 2 22i (i!)2 ((2k − 2i)!)2 2 Γ (k + 2 ) i=0 k = ((2k)!)2 where Ỹ (x) = Y (x)/ function  √ π 2k Γ(k+ 12 ) X ρ̃2k−2i (kx − yk) π − 1, 22k Γ2 (k + 12 ) i=0 22i (i!)2 (2k − 2i)! 1/2k (23) is a zero-mean unit variance Gaussian HIRF with the covariance 2k Γ(k + 21 ) √ π ρ̃(kx − yk) = !1/k ρY (kx − yk). Notice, that for i = k in (23) by the Legendre duplication formula ((2k)!)2 π Γ2 (2k + 1)π (2k)2 Γ2 (2k)π = = 1. = 1 1 24k k 2 22−4k πΓ2 (2k) 22k Γ2 (k + 2 )22k (k!)2 24k Γ2 (k + 2 )k 2 Γ2 (k) Hence, k−1 2i ((2k)!)2 π X (2k Γ(k + 12 ))2+ k Cov(Λ(x), Λ(y)) = 2k 2 ρ̃2k−2i (kx − yk). 2 Γ (k + 12 ) i=0 22i (i!)2 (2k − 2i)!π 1/2k ′ ′ ′ Therefore, if |ρ̃(r)| ≤ C e−γ r , r > 0, γ > 0, then the covariance function of Model 5 satisfies the condition |ρΛ (r)| ≤ Ce−γr , r > 0, γ > 0, and the Rényi function equals 1 1 T (q) = q − 1 − logb EY 2kq (x) = q − 1 − logb 2 2 Finally, we obtain 18 2kq Γ(kq + 12 ) √ π ! . (24) −γr Theorem 7.2 Suppose , r > 0,  that for Model  5 the correlation function of Y (x) satisfies |ρY (r)| ≤ Ce 2C p σΛ 3 2 1 + σΛ , e 3 γ > 0, and b > max . d Then the measures µk − → µ, k → ∞, and the Rényi function is given by (24) for q ∈ [1, 2]. The following model shows how vector-valued random fields can be used to construct mother fields. Model 6 Let Λ(x) = 2 kY (x), where Y (x) ∼ χ2 (k), and the HIRF field Y (x), x ∈ R3 , has a covariance function ρY (r), r ≥ 0. By properties of the chi-square distribution, it follows that EΛ(x) = 2 EY (x) = 1, k V arΛ(x) = Cov(Λ(x), Λ(y)) = 4 2 V arY (x) = < +∞, k2 k 4 ρY (kx − yk). k2 Notice that if Y (x) = 21 (Z12 (x) + ... + Zk2 (x)), x ∈ R3 , where Zi (x), i = 1, ..., k, are independent zero-mean unit variance components of k-dimensional vector Gaussian HIRF with a covariance function ρZ (r), r ≥ 0 of each component, then Cov(Λ(x), Λ(y)) = 2 4 k 2 · ρZ (kx − yk) = ρ2Z (kx − yk). 2 k 2 k Therefore, Model 6 satisfies Conditions 1 and 2 and |ρΛ (r)| ≤ Ce−γr , r > 0, γ > 0, if |ρY (r)| ≤ C ′ e−γ ′ ′ r or ′ |ρZ (r)| ≤ C ′ e−γ r , r ≥ 0, γ > 0. The corresponding Rényi function is 1 T (q) = q − 1 − logb 2    q   1 1 2 2 q − 1 − logb EY (x) = q 1 − logb k 2 k 2 q Γ(q 2 + k2 ) Γ( k2 ) ! . (25) Summarising, one gets −γr Theorem 7.3 Suppose that , function in Model 6 satisfies the inequality |ρY (r)| ≤ Ce  the correlation 2C p σΛ 3 2 . r > 0, γ > 0, and b > max 1 + σΛ , e 3 d → µ, k → ∞, and for q ∈ [1, 2] the Rényi function is given by (25). Then the measures µk − ′ Note that Γ (x) = ψ(x)Γ(x), where ψ(x) is the digamma function defined by  Z ∞  −t e e−xt ψ(x) = − dt. t 1 − e−t 0 Then, it follows from (13), (14) and (21) that for Model 4 ψ(q + 12 ) 1 logb 2 − , 2 2 ln 2   Γ(q + 21 ) qψ(q + 21 ) 1 √ f (α(q)) = 1 + logb . − 2 2 ln 2 π α(q) = 1 − 19 (26) Analogously, for Model 5 one gets from (24) kψ(kq + 21 ) k logb 2 − , 2 2 ln 2   kqψ(kq + 12 ) Γ(kq + 21 ) 1 √ − . f (α(q)) = 1 + logb 2 2 ln 2 π α(q) = 1 − Finally, it follows from (25) that for Model 6   ψ(q + k2 ) 1 2 1 − logb 2 − , α(q) = 1 − logb 2 k 2 2 ln 2 ! Γ(q + k2 ) qψ(q + k2 ) 1 f (α(q)) = 1 + logb . − k 2 2 ln 2 Γ( 2 ) (27) (28) Summarising, we obtain Theorem 7.4 Let the corresponding conditions of Theorems Theorems 7.1, 7.2 and 7.3 are satisfied for Models 4, 5 and 6. Then the multifractal spectra of these models are given by (26), (27) and (28) respectively. The plots in Figure 2 demonstrate the behaviours of the Rényi functions and multifractal spectra of Models 4, 5 and 6. To plot Figure 2 we used similar settings and coordinate ranges as in Figure 1. The following values of parameters were chosen to produce the plots: b = 2, k = 2 and σY = 1. Figure 2 shows that similar to Models 1 and 2, Models 4, 5 and 6 exhibit a parabolic-type behaviour. The spread of T (q) values is wider for Model 5 when compared to Models 4 and 6 within q ∈ (0, 3). The Rényi functions of Models 4, 5 and 6 manifest a concave down increasing and decreasing behaviour within q ∈ (0, 3). Similar to Models 1, 2 and 3, the multifractal spectra of Models 4, 5 and 6 show a concave down increasing behaviour within q ∈ (0, 3). Finally, we illustrate the impact of parameter b on the Rényi function using Models 1, 2, 3, 4, 5 and 6. For Model 1, σY = 1 was selected. Then, for Models 2 and 3, the parameters λ = 3 and β = 2 were used. The parameter k = 2 was chosen for Models 5 and 6. Figure 3 suggests that the Rényi functions for all the models exhibit a similar pattern. It suggests that the Rényi functions for Models 1, 3 and 4 are more concave than for other models for small values of b. The dispersion of T (q) values on the interval (0, 3) is smallest for Model 3 and largest for Model 5. Concavities of the Rényi functions become smaller when b increases and the other functions have almost linear behaviors for q ∈ (0, 3) for large values of b. 8. Numerical studies 8.1. Simulation methodology There are numerous models for which explicit expressions for the Rényi function in terms of elementary functions or even series are not available. Also, in the majority of cases obtaining explicit mathematical formulae for Rényi functions is a difficult problem and rigorous results were derived only for some ranges of the parameter q. For example, for all models in [23] and this paper, T (q) was derived for q ∈ [1, 2] only. 20 a) Rényi function of Model 4 b) Rényi function of Model 5 c) Rényi function of Model 6 d) Spectrum of Model 4 e) Spectrum of Model 5 f) Spectrum of Model 6 Figure 2: Examples of Rényi functions and multifractal spectra for Models 4, 5 and 6. These results are also likely to be true for wider ranges of q, but showing it requires new proof strategies, see Section 5 and the example in Section 7. For such difficult cases, random field simulations can be used to obtain realizations of random fields from theoretical models and compute empirical Rényi functions. These empirical Rényi functions can be used as a substitute of exact mathematical functions for verifications whether real data are consistent with considered theoretical models. The R package RandomFields provides a wide range of simulation techniques and algorithms for random fields, see [35] for more detail. The function RFSimulate simulates Gaussian random fields for a given covariance or variogram model and parameters defined by the arguments RMnugget and MTrend. To compute a sample Rényi function the ratio log2 E P l (m) q ) µ(Sl (m) log2 |Sl | from (6) for large values of m was used. This ratio was replaced by an empirical estimator that employs the HEALPix structure. For the HEALPix (m) resolution n = 1024 with 12582912 pixels (in the following denoted by i), 196608 sets Sl with 64 pixels per set were used to estimate the ratio. As CMB measurements M (i) can take negative values they were transformed to non-negative ones by subtracting their minimum value: M̃ (i)=M (i) − min (M (i)). Then, 21 a) Rényi functions for Model 1 b) Rényi functions for Model 2 c) Rényi functions for Model 3 d) Rényi functions for Model 4 e) Rényi functions for Model 5 f) Rényi functions for Model 6 Figure 3: Dependence of the Rényi function on the parameter b. (m) the terms E(µ(Sl q P P M̃ (i) . Finally, the empirical Rényi ))q were estimated by µ̂l q = (m) M̃ (i) / i∈S l function was computed as T̂ (q) = q l µ̂l ) . (m) log2 |Sl | log2 ( P The empirical multifractal spectrum was estimated by fˆ(α) = q · α̂(q) − T̂ (q), where α̂(q) = P l ((µ̂l q / P µ̂l q ) · ln(µ̂l )) (m) log2 |Sl | . For the selected large number of pixels T̂ (q) and fˆ(q) provide reliable estimations and can be used for wider intervals of q values than [1, 2]. In this paper we considered intervals (0.5, 3) and (−10, 10) when it was required. Figure 4(a) shows a realization of a multifractal random field in a large spherical window. The field was obtained from a Gaussian mother random field Y (x) with the exponential covariance model and its variance 22 equals 2. This covariance function has an exponential form and obviously satisfies the inequality in (4). As −0.5 0.2 0.0 0.4 0.5 T(q) f(α) 0.6 1.0 0.8 1.5 1.0 an approximation of the limit field a finite product field Λ40 (x) with b = 3 was used. 0.5 1.0 1.5 2.0 2.5 3.0 0.8 q a) Realization of a multifractal b) Sample Rényi function with the random field fitted Log-Normal model 1.0 1.2 1.4 1.6 1.8 α c) f (α) versus α Figure 4: Analysis of a simulated multifractal random field First 40 realizations Yi (bi x), i = 1, ..., 40, were simulated on a sphere by the package RandomFields. Then the finite-product field Λ40 (x) was computed by transforming the simulated values according to the P40 formula exp( i=1 Yi (bi x) − 40). The dot plot of the empirical Rényi function is shown in Figure 4(b). The solid straight line is used as a reference to see departures from the fitted Model 1. It is clear that the empirical Rényi function and the theoretical one from (8) are very close on an interval that is wider than [1, 2]. Figure 4(c) shows the spread of the multifractal spectrum. The simulation studies suggest that the theoretical results from previous sections also hold for intervals wider than in the theorems. 8.2. Computing the Rényi function for CMB data In this section, empirical Rényi functions were calculated for real cosmological data obtained from the NASA/IPAC Infrared Science Archive [30]. Figure 5 gives examples of sky windows CMB data from which were used to get empirical Rényi functions in the following examples. Extensive numerical studies were conducted for different windows in various sky locations. As in all cases we obtained rather similar results, we restrict our presentation only to few typical examples. The R package rcosmo was used for computations and visualizations, see [9, 10] for more details. For small windows the function fRen was slightly modified to change the support of the measure µ from the whole sky to the selected window. First the Rényi function was computed for the whole sky. The obtained sample Rényi function is shown in Figure 6(a) by dots. The straight line in the Figure 6(a) was drawn to assess departures of the sample Rényi function from a linear behaviour. The difference of the sample Rényi function and the linear function 23 a) Large window b) Medium window c) Small and very small windows Figure 5: Different sky windows of CMB data that connects the points (1,0) and (2,1) is shown in Figure 6(b). It is clear that the departure from a linear behaviour is not substantial. Figure 6(c) shows the function α(q) and Figure 6(d) plots the function f (α(q)) versus α(q). As it was discussed in Section 8.1 to compute α(q) and f (α(q)), we used the formula for Rényi functions for the range (−10, 10) as simulation studies and analysis in [12] suggest the same analytical form of the Rényi function as for the range [1, 2]. All these plots confirm only very small multifractality of the CMB data. Similar results were also obtained for different sky windows, see, for example, Figures 7(a), 7(b), 7(d) and 7(e). The Rényi functions, multifractal spectra, similar analysis and plots were produced for different window sizes of the CMB unit sphere. Large, medium, small and very small window sizes with areas 1.231, 0.4056, 0.0596 and 0.0017 were selected, see Figure 5. The Rényi function was computed for small windows located at different places of the sky sphere such as near the pole, near the equator and other places of the sphere. Although different window sizes of the sphere were investigated, there’s not that much of evidence to suggest that we have substantial multifractality. The ranges of y scale in Figure 6(b), 7(b) and 7(e) suggest that this multifractality is very small. These results and the variations of the values of T̂ (q) between windows suggest that collecting data at very fine scales and further tests of hypothesis are required. We plan to develop tests of hypothesis about Rényi functions in future publications and use them for new high resolution data that will be available from the next generation CMB experiments CMB-S4 [19]. As the obtained plots are rather similar we present only two of them for large and small windows in Figure 7. Then, all models from Sections 6 and 7 were used to fit the empirical Rényi function. For the log-normal model we present the results for all windows. For other models, only results for CMB data in a large window are given. Similar results were also obtained for other windows. To fit models to empirical Rényi functions several methods were employed. For the log-normal model the simple linear regression approach was used 24 1.015 0.00012 0.00010 1.0 0.00008 1.010 0.8 1.005 α(q) Difference 0.00006 0.6 T(q) 0.00000 0.995 0.00002 1.000 0.00004 0.4 0.2 0.0 1.0 1.2 1.4 1.6 1.8 2.0 1.0 1.2 1.4 q 1.6 1.8 2.0 −10 −5 0 a) Sample Rényi function versus b) Difference of sample Rényi linear function function and linear function 10 c) α(q) versus q 0e+00 Difference 0.0 0.92 −2e−06 0.2 0.94 0.4 −1e−06 T(q) f(α) 0.96 0.6 0.98 0.8 1e−06 1.0 1.00 5 q q 0.995 1.000 1.005 1.010 1.015 1.0 α d) f (α) versus α 1.2 1.4 1.6 1.8 2.0 1.0 1.2 1.4 q 1.6 1.8 2.0 q e) Sample Rényi function with the f) Difference between sample Rényi fitted Log-Normal model function and the fitted Log-Normal model Figure 6: Whole sky data analysis whereas for the other models the non-linear regression approach was applied. As the Rényi function of the log-normal model is specified by (8), substituting a = 2 σY 4 ln b , results in the form T (q) = a(−q 2 + q) + q − 1. Then the R function “lm” was used for a simple linear regression fit with the intercept 0 to T (q) − q + 1. The values of the parameter a and the root mean square error for deviations of Model 1 from the empirical Rényi function are given in Table 1. Figure 6(e) demonstrates the fit of the log-normal model(shown in the red colour) to the empirical Rényi function. As this plot is rather similar for all other models and windows we present only the plots of residuals in Figures 6(f), 7(c), 7(f) and Figure 8. As the estimated value of a is close to zero, the fit of Model 1 gives an almost degenerated case, when 25 1e−06 0.00012 1.00 0.00002 −2e−06 0.00004 −1e−06 Difference 0.00006 Difference 0.00008 0e+00 0.00010 0.98 0.96 f(α) 0.94 0.00000 0.92 0.995 1.000 1.005 1.010 1.015 1.0 1.2 1.4 α 1.6 1.8 2.0 1.0 1.4 1.6 1.8 2.0 q b) Difference with linear c) Difference with Model 1 for large window function for large window window 0.0e+00 Difference 0.00010 −1.5e−06 0.00000 0.92 −1.0e−06 0.00005 0.96 Difference 0.98 5.0e−07 0.00015 1.00 1.0e−06 a) f (α) versus α for large 0.94 f(α) 1.2 q 0.990 0.995 1.000 1.005 1.010 1.015 α −5.0e−07 0.990 1.0 1.2 1.4 1.6 1.8 2.0 1.0 1.2 1.4 q 1.6 1.8 2.0 q d) f (α) versus α for small e) Difference with linear f) Difference with Model 1 for small window function for small window window Figure 7: Analysis of large and small sky windows data either σY2 is very small or b is very large, which is consistent with the plot in Figure 3. The results in Table 1 also confirm that multifractality is very small as for all observation windows a is almost zero and αmax − αmin is very small. Next, for the log-gamma model specified by (11) we used the reparameterisation A = 2 β ln(b), B = λ−1 and considered the non-linear model T (q) − q + 1 = A−1 (ln(1 − Bx) − x ln(1 − B)). The command “nlsLM” from the R package minpack.lm with appropriate initial values was used to fit the model to the sample values of T̂ (q) − q + 1. The values of estimated parameters were  = 0.029407 and B̂ = 0.005469 with RM SE = 1.7198 × 10−6 . The corresponding values of b, λ and β satisfy the assumptions of Theorem 6.2. 1 For the log-negative-inverse-gamma model given by (12) the reparameterisation A = 2 ln(b) , B = β, √ C = λ and application of “nlsLM” resulted in  = 0.254719, B̂ = 5.695755, Ĉ = 0.386207 and RM SE = 26 Observation window [αmin , αmax ] αmax - αmin Whole Sky [0.9916, 1.0165] RMSE a 0.024917 0.000513 1.3602 · 10−6 Large [0.9908, 1.0167] 0.025846 0.000555 1.3590 · 10−6 Medium [0.9893, 1.0159] 0.026620 0.000629 1.1033 · 10−6 Small [0.9867, 1.0170] 0.030219 0.000745 7.9095 · 10−7 Very Small [0.9842, 1.0543] 0.070150 0.001500 1.3949 · 10−5 Difference −4e−06 −3e−06 −2e−10 −2e−05 −2e−06 0e+00 0e+00 Difference 2e−10 2e−05 4e−10 1e−06 0e+00 −1e−06 Difference 4e−05 2e−06 Table 1: Analysis of different sky windows data with Model 1 1.0 1.2 1.4 1.6 1.8 2.0 1.0 1.2 1.4 q 1.6 1.8 2.0 1.0 1.2 1.4 q 1.8 2.0 b) Difference with Model 3 c) Difference with Model 4 Difference 0e+00 1e−05 −2e−05 −2e−04 −1e−05 −1e−04 0e+00 Difference 1e−04 2e−05 2e−04 3e−05 3e−04 a) Difference with Model 2 1.6 q 1.0 1.2 1.4 1.6 1.8 2.0 1.0 q 1.2 1.4 1.6 1.8 2.0 q d) Difference with Model 5 e) Difference with Model 6 Figure 8: Differences between the sample Rényi function and the fitted models 2.6201 × 10−10 . Note, that the obtained Ĉ corresponds to λ that is outside of Lβ,λ in Theorem 6.3. For parameters satisfying the conditions of Theorem 6.3 RMSE is substantially larger. As the results in Theorem 6.3 might be also true for other parameters(see the discussion in Section 8.1) we used the obtained 27 value of Ĉ. 1 Model 4 was fit by using a linear regression model with the parameter A = − 2 log 2 b. The estimated parameter  was −0.000667 and RM SE = 2.56533 × 10−5 . For Models 5 and 6 the non-linear regression approach and the R function “nls” were used. For Model 5 the estimates were found as  = −0.000762, k̂ ≈ 1 and RM SE = 1.6347 × 10−5 . Finally, for Model 6 the estimated parameters were  = 0.000269, k̂ = 1 and RM SE = 1.8393 × 10−4 . Figures 7(c) and 8 demonstrate that departures of the fitted models from the empirical Rényi function are very small, but have different patters. The numerical studies suggested that Models 2 and 3 are more flexible than the other models. However, to fit these models one has to very carefully choose initial values of the parameters for the nls estimation. Different initial values can lead to different results which can be a potential issue for data which are similarly to CMB show minor multifractality. Also, nls method’s rates of convergence for Models 2 and 3 are very slow. Models 1, 4, 5 and 6 have less parameters and are less flexible than Models 2 and 3. However, in many cases they give a reasonable fit very quickly, are robust to the choice of initial values and more computationally efficient. All models gave a good fit to the empirical Rényi functions. The analysis in this section suggests no significant or very small multifractality for the currently available resolution of CMB measurements. 9. Conclusion This paper investigates the multifractal behaviour of spherical random fields and some applications to cosmological data from the mission Planck. The aim of this paper is to introduce several multifractal models for random fields on a sphere and to propose simpler models where the Rényi function can be computed explicitly. All Rényi functions for the specified models exhibit either parabolic or approximately linear behaviours. We present the Rényi function computations for different CMB sky windows located at different places of the sphere. Finally we fit the specified models to actual CMB data. All models fit to the data. The analysis suggests that there may exist a very minor multifractality of the data. Some related problems and extensions of the current research that would be interesting for future studies: • Develop statistical tests for different types of Rényi functions; • Prove that the theoretical results and the formulae for the Rényi functions are also valid for the values of q outside the interval [1, 2], see [6, 7, 12]; • Study other models based on vector random fields (similar to Model 6), where the Rényi functions can be computed explicitly; • Investigate changes of the Rényi functions depending on evolutions of random fields driven by SPDEs on the sphere, see [2, 4, 5]; 28 • Apply the developed models and methodology to new high-resolution CMB data from future CMB-S4 surveys [19]. Acknowledgments This research was partially supported under the Australian Research Council’s Discovery Projects funding scheme (project number DP160101366). References [1] Abramowitz, M., & Stegun, I. A. (1948). Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables. Dover Publications, New York. [2] Anh, V. V., Broadbridge, P., Olenko, A., & Wang, Y. G. (2018). On approximation for fractional stochastic partial differential equations on the sphere. Stoch Environ Res Risk Assess, 32(9), 2585–2603. [3] Anh, V. V., Leonenko, N. N., & Shieh, N.-R. (2008). Multifractality of products of geometric Ornstein-Uhlenbeck-type processes. Adv Appl Probab, 40(4), 1129–56. [4] Broadbridge, P., Kolesnik, A. D., Leonenko, N., & Olenko, A. (2019). Random spherical hyperbolic diffusion. J Stat Phys, 177(5), 889–916. [5] Broadbridge, P., Kolesnik, A. D., Leonenko, N., Olenko, A., & Omari, D. (2020). Spherically restricted random hyperbolic diffusion. Entropy, 22(2), 217. [6] Denisov, D., & Leonenko, N. (2016). Limit theorems for multifractal products of geometric stationary processes. Bernoulli, 22(4), 2579–2608. [7] Denisov, D., & Leonenko, N. (2016). Multifractal scenarios for products of geometric Levy-based stationary models. Stoch Anal Appl , 34(4), 610–643. [8] Falconer, K. J. (1994). The multifractal spectrum of statistically self-similar measures. J Theor Probab, 7(3), 681–702. [9] Fryer, D., Li, M., & Olenko, A. (2020). rcosmo: R Package for Analysis of Spherical, HEALPix and Cosmological Data. arXiv:1907.05648 will appear in The R Journal. [10] Fryer, D., Olenko, A., Li, M., & Wang, Y. (2019). rcosmo: Cosmic Microwave Background Data Analysis. URL: https://CRAN.R-project.org/package=rcosmo R package version 1.1.0. [11] Grahovac, D. (2020). Multifractal processes: Definition, properties and new examples. Chaos Solitons Fractals, 134 , 109735. [12] Grahovac, D., & Leonenko, N. (2018). Bounds on the support of the multifractal spectrum of stochastic processes. Fractals, 26(4), 1850055. [13] Grahovac, D., & Leonenko, N. N. (2014). Detecting multifractal stochastic processes under heavy-tailed effects. Chaos Solitons Fractals, 65 , 78–89. [14] Harte, D. (2001). Multifractals: Theory and Applications. Chapman and Hall/CRC, Boca Raton. [15] Hill, J. C. (2018). Foreground biases on primordial non-Gaussianity measurements from the CMB temperature bispectrum: Implications for Planck and beyond. Phys Rev D, 98(8), 083542. [16] Jaffard, S. (1999). The multifractal nature of Lévy processes. Probab Theory Relat Fields, 114(2). [17] Janson, S. (1997). Gaussian Hilbert Spaces. Cambridge University Press, Cambridge. [18] Kahane, J.-P. (1987). Positive martingales and random measures. Chinese Ann Math B , 8(1), 1–12. [19] Kevork, A., Graeme, A., Peter, A., Zeeshan, A., Steven W., A., David, A. et al. (2019). CMB-S4 Science Case, Reference Design, and Project Plan. arXiv:1907.04473. 29 [20] Kogut, A., Banday, A., Bennett, C., Górski, K., Hinshaw, G., Smoot, G. et al. (1996). Tests for non-Gaussian statistics in the DMR four-year sky maps. Astrophys J Lett, 464(1), L29–L33. [21] Lang, A., & Schwab, C. (2015). Isotropic Gaussian random fields on the sphere: Regularity, fast simulation and stochastic partial differential equations. Ann Appl Probab, 25(6), 3047–94. [22] Leonenko, N. (1999). Limit Theorems for Random Fields with Singular Spectrum. Springer, Dordrecht. [23] Leonenko, N., & Shieh, N.-R. (2013). Rényi function for multifractal random fields. Fractals, 21(2), 1350009. [24] Malyarenko, A. (2012). Invariant Random Fields on Spaces with a Group Action. Springer-Verlag, Berlin. [25] Mannersalo, P., Norros, I., & Riedi, R. H. (2002). Multifractal products of stochastic processes: construction and some basic properties. Adv Appl Probab, 34(4), 888–903. [26] Marinucci, D. (2004). Testing for non-Gaussianity on cosmic microwave background radiation: A review. Stat Sci, 19(2), 294–307. [27] Marinucci, D., & Peccati, G. (2011). Random Fields on the Sphere: Representation, Limit Theorems and Cosmological Applications. Cambridge University Press, New York. [28] Minkov, M., Pinkwart, M., & Schupp, P. (2019). Entropy methods for CMB analysis of anisotropy and non-Gaussianity. Phys Rev D, 99(10), 103501. [29] Molchan, G. (1996). Scaling exponents and multifractal dimensions for independent random cascades. Commun Math Phys, 179(3), 681–702. [30] NASA/IPAC infrared science archive (2019). https://irsa.ipac.caltech.edu/data/Planck/release_2/all-sky-maps/ maps/component-maps/cmb/. Accessed 24 September 2019. [31] Novikov, D., Schmalzing, J., & Mukhanov, V. (2000). On non-Gaussianity in the Cosmic Microwave Background. Astron Astrophys, 364(1). [32] Peccati, G., & Taqqu, M. S. (2011). Wiener Chaos: Moments, Cumulants and Diagrams: A Survey with Computer Implementation. Springer-Verlag, Milan. [33] Planck and the cosmic microwave background (2020). https://www.esa.int/Science_Exploration/Space_Science/Pla nck/Planck_and_the_cosmic_microwave_background. Accessed 4 April 2020. [34] Riedi, R. H. (2002). Multifractal processes. In P. Doukhan, G. Oppenheim, & M. Taqqu (Eds.), Theory and Applications of Long-Range Dependence (pp. 625–716). BirkhÃďuser, Basel. [35] Schlather, M., Malinowski, A., Oesting, M., Boecker, D., Strokorb, K., Engelke, S. et al. (2019). RandomFields: Simulation and Analysis of Random Fields. URL: https://cran.r-project.org/package=RandomFields R package version 3.3.6. [36] Starck, J.-L., Aghanim, N., & Forni, O. (2004). Detection and discrimination of cosmological non-Gaussian signatures by multi-scale methods. Astron Astrophys, 416(1), 9–17. [37] The Cosmic Microwave Background (2020). http://planck.cf.ac.uk/science/cmb. Accessed 4 April 2020. 30