On an Auxiliary Function for Log-Density Estimation

https://doi.org/10.7892/boris.36904 | downloaded: 31.1.2022

University of Bern

Institute of Mathematical Statistics and Actuarial Science Technical Report 71

On an Auxiliary Function for Log-Density Estimation

Madeleine L. Cule and Lutz D¨umbgen (University of Cambridge and University of Bern)

July 2008, minor revisions in January 2016

Abstract

In this note we provide explicit expressions and expansions for a special functionJwhich appears in nonparametric estimation of log-densities. This function returns the integral of a log-linear function on a simplex of arbitrary dimension. In particular it is used in theR- packageLogCondDEAD by Cule et al. (2007).

Contents

1 Introduction 2

2 The special functionJ(·) 3

2.1 Definition ofJ(·) . . . 3 2.2 A first recursion formula . . . 3 2.3 Another recursion formula . . . 4

3 An expansion forJ(·) 5

4 A recursive implementation ofJ(·)and its partial derivatives 6

5 The special casesd= 1andd= 2 7

5.1 General considerations about a bivariate function . . . 7 5.2 More details for the cased= 1 . . . 8 5.3 The cased= 2 . . . 10 6 Gamma and multivariate beta (Dirichlet) distributions 10

1 Introduction

Suppose one wants to estimate a probability densityfon a certain compact regionC⊂R^d, based on an empirical distributionPˆof a sample fromf. One possibility is to embedCinto a union

S =

m

[

j=1

S_j

of simplicesS_j ⊂ R^d with pairwise disjoint interior. By a simplex inR^d we mean the convex hull ofd+ 1points. Then we consider the familyG=G(S₁, . . . , Sm)of all continuous functions ψ:S→Rwhich are linear on each simplexS_j. Now

ψˆ := arg max

ψ∈G

Z

S

ψ dPˆ− Z

S

exp(ψ(x))dx

(1) defines a maximum likelihood estimatorfˆ:= exp( ˆψ)of a probability density onS, based onP.ˆ For existence and uniqueness of this estimator see, for instance, Cule et al. (2008).

To computeψˆexplicitly, note thatψ ∈ G is uniquely determined by its values at the corners (extremal points) of all simplicesS_j, andR

ψ dPˆ is a linear function of these values. The second integral in (1) may be represented as follows: LetSjbe the convex hull ofx0j,x1j, . . . ,xdj ∈R^d, and setyij :=ψ(xij). Then

Z

S

exp(ψ(x))dx =

m

X

i=1

Z

Si

exp(ψ(x))dx =

m

X

i=1

D_j·J(y_0j, y_1j, . . . , y_dj),

where

D_j := det

x_1j−x_0j,x_2j−x_0j, . . . ,x_dj−x_0j , whileJ(·)is an auxiliary function defined and analyzed subsequently.

2 The special function J(·)

2.1 Definition ofJ(·) Ford∈Nlet

T_d :=

n

u∈(0,1)^d:

d

X

i=1

ui <1 o

. Then fory₀, y₁, . . . , y_d∈Rwe define

J(y0, y1, . . . , y_d) :=

Z

T_d

exp

(1−u+)y0+

d

X

i=1

uiyi

du

withu+:=Pd i=1ui.

Standard considerations in connection with beta- and gamma-distributions as described in Sec- tion 6 reveal the following alternative representation:

J(y0, y1, . . . , yd) := 1 d!Eexp

X^d

i=0

Biyi

withBi = Bd,i := Ei Pd

s=0Esand stochastically independent, standard exponential random variablesE₀, E₁, . . . , E_d. This representation shows clearly that J(·) is symmetric in its argu- ments.

An often useful identity is

J(y0, y1, . . . , y_d) = exp(y∗)J(y0−y∗, y1−y∗, . . . , y_d−y∗) for anyy∗ ∈R. (2) 2.2 A first recursion formula

Ford= 1one can computeJ(y₀, y₁)explicitly:

J(y₀, y₁) = Z 1

0

exp (1−u)y₀+uy₁ du =







exp(y₁)−exp(y₀) y1−y0

ify0 6=y1, exp(y0) ify0 =y1. Ford≥2one may use the following recursion formula:

J(y₀, y₁, . . . , y_d) =











J(y₁, y₂, . . . , y_d)−J(y₀, y₂, . . . , y_d)

y₁−y₀ ify₀ 6=y₁,

∂

∂y1

J(y₁, y₂, . . . , y_d) ify₀ =y₁.

(3)

SinceJ(y0, y1, . . . , yd)is continuous iny0, y1, . . . , yd, it suffices to verify (3) in case ofy0 6=y1. We may identifyT_d with the set

(v,u) : u ∈ T_d−1, v ∈ (0,1−u₊) . Then it follows from Fubini’s theorem that

J(y₀, y₁, . . . , y_d)

= Z

Td−1

Z 1−u+

0

exp

(1−u+−v)y0+vy1+

d

X

i=2

ui−1yi

dv du

= Z

T_d−1

exp (1−u+−v)y0+vy1+Pd

i=2ui−1yi

y₁−y₀

1−u₊ v=0

du

= Z

T_d−1

exp (1−u+)y1+Pd

i=2ui−1yi

−exp (1−u+)y0+Pd

i=2ui−1yi

y₁−y₀ du

= J(y1, y2, . . . , yd)−J(y0, y2, . . . , yd) y1−y0

.

2.3 Another recursion formula

It is well-known that for any integer0≤j < d, Ei

Pj s=0E_s

!j

i=0

, B :=

Pj i=0Ei

Pd s=0Es

, Ei

Pd

s=j+1Es

!d

i=j+1

are stochastically independent withB ∼Beta(j+ 1, d−j); see also Section 6. Hence we end up with the following recursive identity:

J(y0, y1, . . . , yd)

= j!(d−j−1)!

d! E J(By₀, . . . , By_j)J((1−B)y_j+1, . . . ,(1−B)y_d)

= Z 1

0

u^j(1−u)^d−j−1J(uy0, . . . , uyj)J((1−u)yj+1, . . . ,(1−u)yd)du with

J(r) := exp(r).

Here we used the well-known identity Z

(1−u)<sup>`</sup>u<sup>m</sup>du = `!m!

(`+m+ 1)! for integers`, m≥0. (4) Plugging inj=d−1into the previous recursive equation leads to

J(y₀, y₁, . . . , y_d) = Z 1

0

u^d−1J(uy₀, . . . , uyd−1) exp((1−u)y_d)du. (5)

3 An expansion for J(·)

Withy¯:= (d+ 1)⁻¹Pd

i=0yi andzi:=yi−y¯one may write

J(y0, y1, . . . , y_d) = exp(¯y)J(z0, z1, . . . , z_d) by virtue of (2). Note thatz+:=Pd

i=0zi= 0. Asz:= (zi)^d_i=0→0, d!J(z0, z1, . . . , zd)

= 1 +

d

X

i=0

IE(B_i)z_i+1 2

d

X

i,j=0

IE(B_iB_j)z_iz_j+1 6

d

X

i,j,k=0

IE(B_iB_jB_k)z_iz_jz_k+O(kzk⁴).

It follows from Lemma 6.1 that IE

Y^d

i=0

B_i^ki

=

d

Y

i=0

ki! .

[d+k+]k+ for integersk0, k1, . . . , kd≥0.

In particular,

IE(B0) = 1 d+ 1, IE(B₀²) = 2

[d+ 2]₂, IE(B₀B₁) = 1 [d+ 2]₂, IE(B₀³) = 6

[d+ 3]₃, IE(B₀²B1) = 2

[d+ 3]₃, IE(B0B1B2) = 1 [d+ 3]₃. Consequently,P_d

i=0IE(B_i)z_i = IE(B₀)z₊ = 0, [d+ 2]₂

d

X

i,j=0

IE(B_iB_j)z_iz_j =

d

X

i,j=0

1_[i=j]·2 + 1_[i6=j]

z_iz_j

=

d

X

i,j=0

1_[i=j]+ 1 z_iz_j

=

d

X

i=0

z_i²+z₊²

=

d

X

i=0

z_i²,

and

[d+ 3]₃

d

X

i,j,k=0

IE(B_iB_jB_k)z_iz_jz_k

=

d

X

i,j,k=0

1_[i=j=k]·6 + 1[#{i,j,k}=2]·2 + 1[#{i,j,k}=3]

zizjzk

=

d

X

i,j,k=0

1_[i=j=k]·5 + 1[#{i,j,k}=2]+ 1 z_iz_jz_k

= 5

d

X

i=0

z_i³+ 3

d

X

s,t=0

1_[s6=t]z_s²z_t+z₊³

= 5

d

X

i=0

z_i³+ 3

d

X

s=0

z_s²z₊−3

d

X

s=0

z_s³+z₊³

= 2

d

X

i=0

z_i³.

Consequently,

J(y0, y1, . . . , yd) = exp(¯y) 1

d!+ 1

2(d+ 2)!

d

X

i=0

z_i²+ 1 3(d+ 3)!

d

X

i=0

z_i³+O kzk⁴ . (6)

4 A recursive implementation of J(·) and its partial derivatives

By means of (3) and the Taylor expansion (6) one can implement the functionJ(·)in a recursive fashion. In what follows we use the abbreviation

ya:b =

((y_a, . . . , y_b) ifa≤b () ifa > b

To computeJ(y_0:d)we assume without loss of generality thaty₀≤y₁ ≤ · · · ≤y_d. It follows from (3) and symmetry ofJ(·)that

J(y_0:d) = J(y_1:d)−J(y0:d−1) yd−y0

ify06=yd. This formula is okay numerically ifyd−y0is not too small. Otherwise one should use (6). This leads to the the pseudo code in Table 1.

To avoid messy formulae, one can express partial derivatives ofJ(·)in terms of higher order versions ofJ(·)by means of the recursion (3). For instance,

∂J(y_0:d)

∂y0

= lim

→0

J(y₀+, y_1:d)−J(y₀, y_1:d)

= lim

→0 J(y₀, y₀+, y_1:d)

= J(y₀, y₀, y_1:d).

AlgorithmJ ←J(y, d, ) ify_d−y₀< then

¯ y←Pd

i=0y_i/(d+ 1) z₂ ←Pd

i=0(y_i−y)¯ ²/2 z₃ ←Pd

i=0(y_i−y)¯ ³/3

J ←exp(¯y) 1/d! +z₂/(d+ 2)! +z₃/(d+ 3)!

else

J ← J(y_1:d, d−1, )−J(y0:d−1, d−1, )

/(y_d−y₀) end if.

Table 1: Pseudo-code forJ(y)with ordered input vectory.

Similarly,

∂²J(y_0:d)

∂y₀² = lim

→0

J(y₀+, y_1:d)−J(y₀, y_1:d)

−J(y₀, y_1:d)−J(y₀−, y_1:d)

= 2 lim

→0

J(y0, y0+, y_1:d)−J(y0, y0−, y_1:d) 2

= 2 lim

→0 J(y0, y0−, y0+, y_1:d)

= 2J(y0, y0, y0, y_1:d), while

∂²J(y_0:d)

∂y0∂y1

= lim

→0

J(y₀+, y₁+, y_2:d)−J(y₀, y₁+, y_2:d)

− J(y₀+, y₁, y_2:d)−J(y₀, y₁, y_2:d)

= lim

→0

J(y0, y0+, y1+, y2:d)−J(y0, y0+, y1, y2:d)

= lim

→0 J(y0, y0+, y1, y1+, y2:d)

= J(y0, y0, y1, y1, y2:d).

5 The special cases d = 1 and d = 2

For small dimensiondit may be worthwhile to work with non-recursive implementations of the functionJ(·). Here we collect and extend some results of D¨umbgen et al. (2007).

5.1 General considerations about a bivariate function

In view of (3) we consider an arbitrary functionf :R→Rwhich is infinitely often differentiable.

Then

h(r, s) :=







f(s)−f(r)

s−r ifs6=r f⁰(r) ifs=r

defines a smooth and symmetric functionh:R² →Rsuch that h(r, s) = f⁰(r) +f⁰⁰(r)

2 (s−r) +O (s−r)²

ass→r.

Its first partial derivatives of order one and two are given by

∂h(r, s)

∂r =











f(s)−f(r)−f⁰(r)(s−r)

(s−r)² ifs6=r,

f⁰⁰(r)

2 + f⁰⁰⁰(r)

6 (s−r) +O (s−r)²

ass→r,

∂²h(r, s)

∂r² =











2 f(s)−f(r)−f⁰(r)(s−r)

−(s−r)²f⁰⁰(r)

(s−r)³ ifs6=r,

f⁰⁰⁰(r)

3 +f⁰⁰⁰⁰(r)

12 (s−r) +O (s−r)²

ass→r,

∂²h(r, s)

∂r∂s =











(s−r) f⁰(r) +f⁰(s)

−2 f(s)−f(r)

(s−r)³ ifs6=r, f⁰⁰⁰(r)

6 +f⁰⁰⁰⁰(r)

12 (s−r) +O (s−r)²

ass→r.

The other partial derivatives of order one and two follow via symmetry considerations.

5.2 More details for the cased = 1

Recall that

J(r, s) = Z 1

0

exp (1−u)r+us du =







exp(s)−exp(r)

s−r ifr 6=s, exp(r) ifr =s.

This is just the function introduced by D¨umbgen, H¨usler and Rufibach (2007). Let us recall some properties and formulae for the corresponding partial derivatives

J_a,b(r, s) := ∂^a+b

∂r^a∂s^b J(r, s) = Z ₁

0

(1−u)^au^bexp((1−u)r+us)du.

Note first that

Ja,b(r, s) = Jb,a(s, r) = exp(r)Ja,b(0, s−r).

Thus it suffices to derive formulae for(r, s) = (0, y)andb≤a. It follows from (4) that Ja,0(0, y) =

Z 1 0

(1−u)^a

∞

X

k=0

u^k k!y^kdu

=

∞

X

k=0

1 k!

Z 1 0

(1−u)^au^kdu·y^k

=

∞

X

k=0

a!

(k+a+ 1)!y^k

= a!

y^a+1

exp(y)−

a

X

`=0

y^`

`!

.

In particular,

J_1,0(0, y) = exp(y)−1−y y²

= 1 2+ y

6 +y² 24 + y³

120+O(y⁴) (y→0), J_2,0(0, y) = 2(exp(y)−1−y−y²/2)

y³

= 1 3+ y

12+ y² 60 + y³

360+O(y⁴) (y→0), J3,0(0, y) = 6(exp(y)−1−y−y²/2−y³/6)

y⁴

= 1 4+ y

20+ y² 120+ y³

840+O(y⁴) (y→0), J4,0(0, y) = 24(exp(y)−1−y−y²/2−y³/6−y⁴/24)

y⁵

= 1 5+ y

30+ y²

210+ y³

1680+O(y⁴) (y→0).

Another general observation is that J_a,b(r, s) =

Z 1 0

(1−u)^a(1−(1−u))^bexp((1−u)r+us)du

=

b

X

i=0

b i

(−1)ⁱJ_a+i,0(r, s).

In particular,

J_a,1(r, s) = J_a,0(r, s)−J_a+1,0(r, s),

J_a,2(r, s) = J_a,0(r, s)−2J_a+1,0(r, s) +J_a+2,0(r, s).

On the other hand,

J_a,b(0, y) =

∞

X

k=0

y^k k!

Z 1 0

(1−u)^au^k+bdu

=

∞

X

k=0

a![k+b]_b (k+a+b+ 1)!y^k with[r]0 := 1and[r]m :=Qm−1

i=0 (r−i)for integersm >0. In particular, J1,1(0, y) = exp(y)(y−2) + 2 +y

y³

= 1 6 + y

12 +y² 40 + y³

180+O(y⁴) (y →0).

5.3 The cased = 2

Our recursion formula (3) yields

J(r, s, t) =







J(s, t)−J(r, t)

s−r ifr 6=s, J10(r, t) ifr =s.

Because ofJ’s symmetry we may rewrite this in terms of the order statisticsy₍₀₎ ≤y₍₁₎ ≤y₍₂₎of (y_i)²_i=0as

J(r, s, t) =











J(y₍₁₎, y₍₂₎)−J(y₍₀₎, y₍₁₎)

y₍₂₎−y₍₀₎ ify₍₀₎ < y₍₂₎, exp(y₍₀₎)

2 ify₍₀₎ =y₍₂₎.

For fixed third argumentt, this functionJ(r, s, t)corresponds toh(r, s)in Section 5.1 with f(x) :=J(x, t). Thus

∂J(r, s, t)

∂r =











J(s, t)−J(r, t)−J1,0(r, t)(s−r)

(s−r)² ifr 6=s,

J_2,0(r, t)

2 +J_3,0(r, t)(s−r)

6 +O (s−r)²

ass→r.

Moreover,

∂²J(r, s, t)

∂r² =











2 J(s, t)−J(r, t)−J_1,0(r, t)(s−r)

−(s−r)²J_2,0

(s−r)³ ifr6=s,

J3,0(r, t)

3 +J4,0(r, t)(s−r)

12 +O (s−r)²

ass→r,

∂²J(r, s, t)

∂r∂s =











J1,0(r, t) +J1,0(s, t)

(s−r)−2 J(s, t)−J(r, t)

(s−r)³ ifr 6=s,

J_3,0(r, t)

6 +J_4,0(r, t)(s−r)

12 +O (s−r)²

ass→r.

6 Gamma and multivariate beta (Dirichlet) distributions

LetG₀, G₁, . . . , G_mbe stochastically independent random variables withG_i ∼Gamma(a_i)for certain parametersai >0. That means, for any Borel setA⊂(0,∞),

IP(Gi ∈A) = Z

A

Γ(ai)⁻¹y^ai−1exp(−y)dy.

Now we definea+:=Pm

i=0ai,G+:=Pm

i=0Gi and

B˜ := (Gi/G+)^m_i=0, B := (Gi/G+)^m_i=1.

Note that B˜ is contained in the unit simplex in R^m+1, while B is contained in the open set T_m =

u ∈ (0,1)^m : u+ < 1 withu+ := Pm

i=1ui. We also define u0 := 1−u+ for any u∈ T_m.

Lemma 6.1. The random vectorB and the random variableG+ are stochastically independent.

Moreover,

G₊ ∼ Gamma(a₊) whileBis distributed according to the Lebesgue density

f(u) := Γ(a₊) Qm

i=0Γ(ai)

m

Y

i=0

u^a_iⁱ⁻¹

onT_m. For arbitrary numbersk₀, k₁, . . . , k_m ≥0andk₊ :=Pm i=0k_i, IE

Y^m

i=0

B^k_iⁱ

= Γ(a+) Γ(a₊+k₊)

m

Y

i=0

Γ(ai+ki) Γ(a_i) . As a by-product of this lemma we obtain the following formula:

Corollary 6.2. For arbitrary numbersa₀, a₁, . . . , a_m >0, Z

Tm

m

Y

i=0

u^a_iⁱ⁻¹du = Γ(a+)⁻¹

m

Y

i=0

Γ(ai).

Proof of Lemma 6.1. Note that G = (G_i)^m_i=0 my be written as Ξ(G₊,B) with the bijective mappingΞ : (0,∞)× T_m →(0,∞)^m+1,

Ξ(s,u) := (su_i)^m_i=0. Note also that

detDΞ(s,u) = det

















u0 −s −s · · · −s u₁ s 0 · · · 0 u2 0 s . .. ... ... ... ... ... 0 um 0 · · · 0 s

















= det

















1 0 0 · · · 0 u₁ s 0 · · · 0 u2 0 s . .. ...

... ... ... ... 0 um 0 · · · 0 s

















= s^m.

Thus the distribution of(G+,B)has a Lebesgue densityhon(0,∞)× T_mwhich is given by h(s,u) =

m

Y

i=0

Γ(ai)⁻¹Ξ(s,u)^a_iⁱ⁻¹exp(−Ξ(s,u)i)

·

detDΞ(s,u)

=

m

Y

i=0

Γ(a_i)⁻¹(su_i)^ai−1exp(−su_i)

·s^m

= s^a+−1exp(−s)

m

Y

i=0

Γ(ai)⁻¹u^a_iⁱ⁻¹

= Γ(a+)⁻¹s^a+−1exp(−s)·f(u).

Since this is the density ofGamma(a₊)atstimesf(u), we see thatG₊andBare stochastically independent, whereG+ has distributionGamma(a+), and thatf is indeed a probability density onT_mdescribing the distribution ofB.

The fact thatf integrates to one overT_mentails Corollary 6.2. But then we can conclude that IEY^m

i=0

B_i^k(i)

= Z

Tm

m

Y

i=0

u^a_i^i+ki−1du.Z

Tm

m

Y

i=0

u^a_iⁱ⁻¹du

= Γ(a+) Γ(a₊+k₊)

m

Y

i=0

Γ(ai+ki) Γ(a_i) .

References

[1] M.L. CULE, R.B. GRAMACY, and R.J. SAMWORTH (2007). LogConcDEAD, An R package for log-concave density estimation in arbitrary dimensions. Available from http://cran.r-project.org/.

[2] M.L. CULE, R.J. SAMWORTHand M.I. STEWART(2008). Maximum likelihood estimation of a multidimensional log-concave density. Preprint.

[3] L. D ¨UMBGEN, A. H ¨USLERand K. RUFIBACH(2007). Active set and EM algorithms for log- concave densities based on complete and censored data. Technical report 61, IMSV, University of Bern (arXiv:0707.4643).