# Linear Models with Response Functions Based on the Laplace Distribution: Statistical Formulae and Their Application to Epigenomics

Authors: C. Z. W. Hassell Sweatman; G. C. Wake; A. B. Pleasants; C. A. McLean; A. M. Sheppard
Journal: ISRN Probability and Statistics
 (2013)

---

## Abstract

The statistical application considered here arose in epigenomics, linking the DNA methylation proportions measured at specific genomic sites to characteristics such as phenotype or birth order. It was found that the distribution of errors in the proportions of chemical modification (methylation) on DNA, measured at CpG sites, may be successfully modelled by a Laplace distribution which is perturbed by a Hermite polynomial. We use a linear model with such a response function. Hence, the response function is known, or assumed well estimated, but fails to be differentiable in the classical sense due to the modulus function. Our problem was to estimate coefficients for the linear model and the corresponding covariance matrix and to compare models with varying numbers of coefficients. The linear model coefficients may be found using the (derivative-free) simplex method, as in quantile regression. However, this theory does not yield a simple expression for the covariance matrix of the coefficients of the linear model. Assuming response functions which are𝒞2 except where the modulus function attains zero, we derive simple formulae for the covariance matrix and a log-likelihood ratio statistic, using generalized calculus. These original formulae enable a generalized analysis of variance and further model comparisons.


---

## Body

## 1. Introduction and Motivation

This work arose in a biological context, in epigenomics, namely, the modelling of the distribution of errors in the proportions of chemical modification (methylation) on DNA, measured at specific genomic sites (CpG sites). It was observed that this error distribution may be suitably modelled by a truncated Laplace distribution perturbed by a Hermite polynomial.This error distribution was first noticed in Sequenom measurements but has wider application. A survey of data generated by measurements on the Infinium, Illumina, Affymetrix, and MeDIP2 machines showed similar characteristics to that of the Sequenom, where such an amended Laplace distribution was required to properly describe the probability density function. It is thought that variation in the scattering angle of light in the measurement processes common to all of these platforms is responsible for the frequencies in the tails of the measurement distributions not conforming to a simple Laplace density and requiring our proposed amendment. Without the amendment the Laplace density gives tail probabilities for the deviations that are too high, potentially leading to an incorrect failure to reject a null hypothesis. Because the observed frequency distribution of epigenomic and gene expression measurements appears to be a common feature of molecular biology, it is important that the process of estimation and inference under the amended Laplace probability density be studied. The paper reports results from a study of estimation and inference under the amended Laplace density.We extend the theory of linear models as given in [1] to deal with response variables with distributions more general than the exponential family. We consider the Laplace distribution, and amendments thereof, with probability density functions which have abrupt changes in gradient due to the presence of the modulus function. The theory in this paper corresponds to least absolute error (LAE) (or least absolute deviation (LAD)) regression [2–4], also called median regression [5], when the response function is the Laplace distribution without modification. We focus on coefficient estimation for a linear model with a response variable distribution assumed to be a truncated and/or perturbed version of the Laplace distribution, estimating the standard errors of these coefficients and understanding the asymptotic theory.The theory of generalized linear models as described in [1] covers the case of distributions from the exponential family. These distributions have probability density functions which are twice continuously differentiable (𝒞2), everywhere on their support. The usual expressions for the standard errors of the model coefficients for the generalized linear models in [1] are derived using Taylor series and assume distributions with probability density functions which are 𝒞2 everywhere on their support. They cannot be applied to our model due to the presence of the modulus function.The modified Laplace probability density functions considered here have a sharp peak at the maximum. Maximum likelihood estimation (MLE) of coefficients may be done by non-gradient methods, such as the simplex method. However, the usual classical expressions for the standard errors of the coefficients, the information matrix, and the log-likelihood ratio statistic do not apply due to lack of differentiability. We derive expressions for generalized versions of these quantities using generalized functions. Consequently, we show that our MLE is asymptotically normal.The method we present to estimate these statistics could in principle be applied to other probability density functions exhibiting abrupt changes in gradient. Response function parameters are assumed known or previously estimated. The theory is applied to find the standard errors for coefficients of a linear model, assuming the response function has a truncated Laplace distribution with added kurtosis, due to perturbation by a Hermite polynomial. To illustrate the application, we show how birth order can be linked to methylation status at two CpG sites in the promotor of the H19 gene.


## 2. The Model

### 2.1. The Expectation Is Linear

Let(1)y=(y1,y2,…,yn)T∈ℝn
be a vector of response variables; let
(2)X=Xn,m=(1x12⋯x1m1x22⋯x2m⋮1xn2⋯xnm)
be an n×m matrix of explanatory variables (real-valued). The subscripts denote the dimensions and will be omitted when these are assumed fixed (in Sections 2 and 3). Let
(3)β=(β1,β2,…,βm)T∈ℝm
be a vector of coefficients for our linear model and assume that
(4)E(y)=Xβ.
Then each component of the deviation (or error) vector
(5)z=y-Xβ
has expectation zero. The explanatory variables may be continuous or discrete. We assume n≥m. Let rX denote the rank of X; then rX≤m. In practice, we usually have n, the number of data points, much larger than m, the number of coefficients. Our goal is estimating the components of β by ML principles, and determining their standard errors, given a set of response variables y, explanatory variables X, and a response variable distribution based on the Laplace distribution as described below. In terms of generalized linear models, the link function is assumed to be the identity.


### 2.2. The Distribution of the Deviations—A Modified Laplace Distribution

Example 1.
Letf:ℝ→ℝ be defined by
(6)f(z;p)=(pQ(p))e(-p|z|),
where p>0 is a real-valued parameter and Q(p) is a real-valued normalizing function defined so that
(7)∫-∞∞f(z;p)dz=1.
Then f is the probability density function for the Laplace distribution, with scale parameter p, centred at the origin, and with unbounded support. It is not differentiable at the origin in the classical sense.
The method of MLE for the response function (6) corresponds to least absolute error (LAE) regression [2–4]. However, the theory of LAE regression is not sufficiently general for our epigenomic modelling problem. We next describe the more general response functions we require.Example 2.
Now consider the case of bounded support. For finiteB>0, define f:[-B,B]→ℝ by
(8)f(z;p)=(pQ(p;B))e(-p|z|),
where p>0 is a real-valued scale parameter and Q(p;B) is a real-valued normalizing function defined so that
(9)∫-BB‍f(z;p)dz=1.
Then f is the probability density function for the truncated Laplace distribution with scale parameter p, centred at the origin, and with bounded support [-B,B].Example 3.
More generally, consider perturbations of the truncated Laplace probability density function of the following form. Let(10)f(z;p,q)=(pQ(p,g,q;B))e(-p|z|)g(|z|;q),
where real-valued g(z;q) is equal to the constant map
(11)g1(z)=1
plus a perturbation, q is a vector of parameters for g, and parameter vector p=(p,q)∈ℝr, r≥1. If r=1, g has no parameters. Here, Q(p,g,q;B) is a real-valued normalizing function defined so that
(12)∫-BB‍f(z;p)dz=1.
We assume that there exists some ϵ>0 such that g(z;q) is 𝒞3 in z on (-ϵ,B+ϵ) and that g(z;q)>0 on [0,B], for fixed parameter vector q. Note that, as a consequence of using the modulus function, g(|z|;q) will not be differentiable with respect to z, at z=0, in general. As in Example 1, f is not differentiable at the origin due to the use of the modulus function.Example 4.
We could allow unbounded support if∫-∞∞‍e(-p|z|)g(|z|;q)dz is finite.Example 5.
Now consider our motivating example, a truncated Laplace distribution with bounded support[-1,1], perturbed by adding kurtosis. Such a distribution is used to model the deviations in the proportions of methylation measured at gene promoter CpG sites. Specifically, to fit with observations, kurtosis is added to a Laplace probability density function with bounded support by adding a third-order Hermite polynomial to give an amended version.
Consider(13)f(z;p)=(pQ(p,g2,q;B))e(-p|z|)g2(|z|;q).
Here, p=(p,q), q≥0 is small, B=1,
(14)g2(z;q)=1+qH3(z),
and H3(z)=z3-3z is the third-order Hermite polynomial. Solving for Q yields
(15)f(z;p,q)=p4e(-p|z|)[1+qH3(|z|)]2[(p3-3qp2+6q)-e-p(p3(1-2q)+6pq+6q)].Example 6.
The functionsg3(z;q)=1-qz and g4(z;q)=e-qz2, assuming small positive q and bounded support, could be used in (10) to model distributions similar to the Laplace but with thinner tails.
We restrict to symmetric distributions satisfyingf(z)=f(-z).


## 2.1. The Expectation Is Linear

Let(1)y=(y1,y2,…,yn)T∈ℝn
be a vector of response variables; let
(2)X=Xn,m=(1x12⋯x1m1x22⋯x2m⋮1xn2⋯xnm)
be an n×m matrix of explanatory variables (real-valued). The subscripts denote the dimensions and will be omitted when these are assumed fixed (in Sections 2 and 3). Let
(3)β=(β1,β2,…,βm)T∈ℝm
be a vector of coefficients for our linear model and assume that
(4)E(y)=Xβ.
Then each component of the deviation (or error) vector
(5)z=y-Xβ
has expectation zero. The explanatory variables may be continuous or discrete. We assume n≥m. Let rX denote the rank of X; then rX≤m. In practice, we usually have n, the number of data points, much larger than m, the number of coefficients. Our goal is estimating the components of β by ML principles, and determining their standard errors, given a set of response variables y, explanatory variables X, and a response variable distribution based on the Laplace distribution as described below. In terms of generalized linear models, the link function is assumed to be the identity.


## 2.2. The Distribution of the Deviations—A Modified Laplace Distribution

Example 1.
Letf:ℝ→ℝ be defined by
(6)f(z;p)=(pQ(p))e(-p|z|),
where p>0 is a real-valued parameter and Q(p) is a real-valued normalizing function defined so that
(7)∫-∞∞f(z;p)dz=1.
Then f is the probability density function for the Laplace distribution, with scale parameter p, centred at the origin, and with unbounded support. It is not differentiable at the origin in the classical sense.
The method of MLE for the response function (6) corresponds to least absolute error (LAE) regression [2–4]. However, the theory of LAE regression is not sufficiently general for our epigenomic modelling problem. We next describe the more general response functions we require.Example 2.
Now consider the case of bounded support. For finiteB>0, define f:[-B,B]→ℝ by
(8)f(z;p)=(pQ(p;B))e(-p|z|),
where p>0 is a real-valued scale parameter and Q(p;B) is a real-valued normalizing function defined so that
(9)∫-BB‍f(z;p)dz=1.
Then f is the probability density function for the truncated Laplace distribution with scale parameter p, centred at the origin, and with bounded support [-B,B].Example 3.
More generally, consider perturbations of the truncated Laplace probability density function of the following form. Let(10)f(z;p,q)=(pQ(p,g,q;B))e(-p|z|)g(|z|;q),
where real-valued g(z;q) is equal to the constant map
(11)g1(z)=1
plus a perturbation, q is a vector of parameters for g, and parameter vector p=(p,q)∈ℝr, r≥1. If r=1, g has no parameters. Here, Q(p,g,q;B) is a real-valued normalizing function defined so that
(12)∫-BB‍f(z;p)dz=1.
We assume that there exists some ϵ>0 such that g(z;q) is 𝒞3 in z on (-ϵ,B+ϵ) and that g(z;q)>0 on [0,B], for fixed parameter vector q. Note that, as a consequence of using the modulus function, g(|z|;q) will not be differentiable with respect to z, at z=0, in general. As in Example 1, f is not differentiable at the origin due to the use of the modulus function.Example 4.
We could allow unbounded support if∫-∞∞‍e(-p|z|)g(|z|;q)dz is finite.Example 5.
Now consider our motivating example, a truncated Laplace distribution with bounded support[-1,1], perturbed by adding kurtosis. Such a distribution is used to model the deviations in the proportions of methylation measured at gene promoter CpG sites. Specifically, to fit with observations, kurtosis is added to a Laplace probability density function with bounded support by adding a third-order Hermite polynomial to give an amended version.
Consider(13)f(z;p)=(pQ(p,g2,q;B))e(-p|z|)g2(|z|;q).
Here, p=(p,q), q≥0 is small, B=1,
(14)g2(z;q)=1+qH3(z),
and H3(z)=z3-3z is the third-order Hermite polynomial. Solving for Q yields
(15)f(z;p,q)=p4e(-p|z|)[1+qH3(|z|)]2[(p3-3qp2+6q)-e-p(p3(1-2q)+6pq+6q)].Example 6.
The functionsg3(z;q)=1-qz and g4(z;q)=e-qz2, assuming small positive q and bounded support, could be used in (10) to model distributions similar to the Laplace but with thinner tails.
We restrict to symmetric distributions satisfyingf(z)=f(-z).


## 3. Maximum Likelihood Estimation

### 3.1. The Log-Likelihood Function

Let(16)f(z,p)
be a probability density function with parameter vector p, as described in Section 2. Let z1,…,zn be a sequence of independent and identically distributed deviations with joint probability density function
(17)f(z1,…,zn;p)=f(z1(β),…,zn(β);p)=Πi=1nf(zi(β);p).
This is also the likelihood function
(18)Lz(z;p)=Lβ(β;p;X,y)=L(z(β);p)=f(z(β);p),
which may be regarded as a function of z or β; here, the subscript reflects our point of view. We use the log-likelihood function l in the estimation of β, where, using various notation,
(19)lz(z;p)=loge(Lz(z;p))=lβ(β;p;X,y)=l(z(β);p).
Substituting measured values of yi and known inputs xij into lβ, we obtain a function of β and p. The parameters p are assumed known, but if not, may be estimated separately. In our biological application, they are estimated by MLE and are assumed fixed for a particular measuring process. Hence, we have a function of β, the coefficients of our linear model.Our aim is to find a maximum likelihood estimator (MLE) denotedβ^n∈ℝm, that is, some point at which l attains its maximum value. The subscript n corresponds to the number of deviations. Now l is a continuous function. If B is finite, it has compact support in ℝm. Since a continuous function on a compact set attains its supremum, the existence of a MLE for lβ is guaranteed. Even if B=∞, we may consider truncations with finite bounds Bk=k, k=1,2,…. Since l is maximized when the zi are small, truncating f to [-k,k] for k large enough will not affect the set of points at which l attains its maximum. We show in Section 3.3 that a MLE β^n is not necessarily unique.Case 1.
Ifg(z;q)=g1(z)=1 and so f is the Laplace probability density function with parameter p as in (6) or (8), then, for β∈ℝm such that every zi(β) is in the support of f,
(20)lβ(β;p;X,y)=lβ,n,m(β;p;Xn,m,y)=-nlog(Q)+nlog(p)+∑i=1n(-p|zi(β)|)=-nlog(Q)+nlog(p)+∑i=1n(-p|yi-(Xβ)i|),
where Q=2 if the support of f is ℝ. If the support of f is [-B,B] then
(21)Q=Q(p,g1;B)=2(1-e-pB).
The theory of LAE regression (corresponding to MLE using Laplace distributions without modification as response functions) may be found in various texts, for example, [3]. Here, it is proved that there exists a MLE β^n corresponding to at least rX zero errors. We are concerned with the extension of these ideas to the case of perturbed and truncated Laplace response functions. For the truncated Laplace distribution, we prove that there exists a MLE β^n corresponding to at least rX zero errors. Consider the following more general case.Case 2.
Iff is a perturbed Laplace probability density function with perturbing function g(z;q) and bounded support as in (10), then for β∈ℝm such that |zi(β)|≤B, i=1,2,…,n,
(22)lβ(β;p;X,y)=lβ,n,m(β;p;Xn,m,y)=-nlog(Q(p,g,q;B))+nlog(p)+∑i=1n(-p|zi(β)|)+∑i=1nlog(g(|zi(β)|;q))=-nlog(Q(p,g,q;B))+nlog(p)+∑i=1n(-p|yi-(Xβ)i|)+∑i=1nlog(g(|yi-(Xβ)i|;q)).
Note since g(zi;q) is strictly positive on [0,B], so is f(zi;p) and so log(f(zi;p)) is well-defined on [-B,B], i=1,2,…,n.
In Section3.3, we show that if the perturbing function g is such that logg(z;q) is convex and non-increasing on [0,B], there exists a MLE corresponding to at least rX data points. We give an upper bound on |dlogg(z;q)/dz| on [0,B] which, if not exceeded, ensures that there exists a MLE corresponding to at least one data point. We apply these results to our motivating example, the Laplace distribution with added kurtosis, described below.Example 7.
Ifg(z;q)=g2(z;q)=1+qH3(z) and so f is a Laplace probability density function with scale parameter p, with added kurtosis and bounded support as in (13), then for β∈ℝm such that |zi(β)|≤1, i=1,2,…,n,
(23)lβ(β;p,q;X,y)=lβ,n,m(β;p;Xn,m,y)=-nlog(Q(p,g2,q;1))+nlog(p)+∑i=1n(-p|zi(β)|)+∑i=1nlog(1+qH3(|zi(β)|))=-nlog(Q(p,g2,q;1))+nlog(p)+∑i=1n(-p|yi-(Xβ)i|)+∑i=1nlog(1+qH3(|yi-(Xβ)i|)).
Now lβ is not differentiable in the classical sense with respect to the linear model coefficients βj when any zi=0. Hence, we cannot assume that lβ is differentiable at a MLE. This paper addresses this issue firstly by proposing a non-gradient method of coefficient estimation and secondly (in Section 4) by using generalized functions to calculate statistical estimates including estimates of standard errors. In Section 5, we discuss the asymptotic theory of our MLE.


### 3.2. Coefficient Estimation Dealing with Abrupt Changes in Gradient

AlthoughLβ and lβ are continuous functions of β, their first derivatives are not. Consider the geometry of the coefficient space ℝm, where β∈ℝm. For each index i, since the response function distribution is defined in terms of absolute values, we can find a hyperplane Hi0 in ℝm on which L=Lβ and l=lβ are not differentiable, defined by setting zi=0. Let xiT=(xi1,xi2,…,xim), the ith row of X, it is never the zero vector since xi1=1, i=1,2,…,n. Choose w∈ℝm so that xiTw=0. Then
(24)β=w+yixixiTxi
yields zi=0. Let Hi0 be the set of all such β. For example, for m=2, for each error term there is a line in ℝ2 on which Lβ and lβ have a sharp ridge. By inspection of the geometry, we would expect the values of β which maximize lβ to be either on the union of the hyperplanes or very close to intersections of the hyperplanes Hi0, i=1,2,…,n. Imagine searching in β-space near the hyperplanes Hi0. Even if Lβ has a local maximum near but not on the union of the hyperplanes, it would be difficult to use a method based on the gradient of either Lβ or lβ, since the gradient changes sharply whenever we cross one of the Hi0. The simplex method of coefficient estimation, which does not require any partial derivatives, suits this geometry.


### 3.3. The Maximum Likelihood Estimator Corresponds to a Data Point

#### 3.3.1. Convex and Non-Increasing Perturbations of the Laplace Probability Density Function

Consider the probability density function(25)f(z;p,q)=(pQ(p,g,q;B))e(-p|z|)g(|z|;q),
with support [-B,B], for some finite B>0 as described in Section 2.2 (Example 3). Recall we assume that there exists some ϵ>0 such that g(z;q) is 𝒞3 in z on (-ϵ,B+ϵ) and that g(z;q)>0 on [0,B], for fixed parameter q. Let
(26)ΩB={z∈ℝn:|zi|≤B,i=1,…,n}=[-B,B]n.
We consider the log-likelihood function
(27)lz(z;p):ΩB⟶ℝ
(conditional on p) with its domain restricted to
(28)𝒜={z(β):β∈ℝm}∩ΩB={y-Xβ:β∈ℝm}∩ΩB,
that is, constrained to 𝒜 or equivalently, consider the log-likelihood function
(29)lβ(β;p;X,y):ℬ⟶ℝ,
constrained to ℬ, where ℬ⊂ℝm is defined as {β∈ℝm:z(β)∈ΩB}. Note that 𝒜=𝒜(X,y,B) and similarly ℬ=ℬ(X,y,B). Also, if the function Z is defined by
(30)Z:ℝm⟶ℝn,β⟼z(β)=y-Xβ,
then Z(ℬ)=𝒜.Lemma 8.
Iflog(g):(-ϵ,B+ϵ)→ℝ is convex and non-increasing (i.e., dlogg(z;q)/dz≤0 on [0,B], then there exists a maximum of lβ:ℬ→ℝ corresponding to at least rX data points. That is, there exists β^n∈ℬ⊂ℝm such that the constrained lβ attains its maximum at β^n, and there exists at least rX indices ij∈{1,2,…,n} such that zij(β^n)=0, j=1,2,…,rX.Corollary 9.
Letf(z;p) be the Laplace probability density function (8) with support [-B,B] and parameter p. Then there exists a maximum of lβ:ℬ→ℝ corresponding to at least rX data points.Proof of Corollary9.
Letg be the constant map g1(z)=1; then log(g(z))=0, and hence log(g(z)) is convex and non-increasing on [0,B]. Corollary 9 follows from Lemma 8.Note that Corollary9 could be proved directly by linear programming theory. Linear programming has been applied to the problem of minimising the sum of absolute errors in various applications (see [6], a survey paper, and also [7]). The convex analysis results we require are in Appendix A.Proof of Lemma8.
To begin, assume thatX has full rank m, recall n≥m and that f has bounded support. Then 𝒜 is a compact convex subset of an m-dimensional affine subspace of ℝn. Since the mapping Z is linear and has full rank m, the inverse image ℬ=Z-1(𝒜) is a compact convex subset of ℝm. We partition ℬ⊂ℝm, which is the support of Lβ, into a finite collection of compact convex sets, so that, on each subset, lβ is convex.
LetHiδ be the hyperplane in ℝm corresponding to the error term zi(β)=δ. Then Hi-B and HiB are the hyperplanes in ℝm corresponding to errors zi=-B and zi=B, respectively. It follows that the log-likelihood function
(31)lβ(β;p;X,y)=-nlog(Q(p,g,q;B))+nlog(p)+∑i=1n(-p|zi(β)|)+∑i=1nlog(g(|zi(β)|;q))=-nlog(Q(p,g,q;B))+nlog(p)+∑i=1n(-p|yi-(Xβ)i|)+∑i=1nlog(g(|yi-(Xβ)i|;q))
is a convex function in between the the hyperplanes Hi-B, Hi0, and HiB, i=1,…,n. These 3n hyperplanes divide the domain ℬ in ℝm into at most 2n open sets bounded by (but not intersecting) the hyperplanes. Each such open set (and hence its closure) may be labelled by a set of n signs. For any β∈ℝm such that 0<|zi(β)|<B, i=1,…,n; (sgn(z1(β)),sgn(z2(β)),…,sgn(zn(β))) labels the open set containing β.
Next, considerℝn as the union of its orthants, which we denote 𝒪k, k=1,…,2n. We assume the orthants are closed sets. For example, the non-negative orthant is {z∈ℝn:zi≥0,i=1,…,n}. In the interior of any 𝒪k, the sign of zi does not change, i=1,…,n. Relabel the open subsets ℬk=ℬk(X,y,B), where k∈{1,2,…,2n}, so that Z(ℬk)⊂int(𝒪k). Let 𝒜k=Z(ℬk)=𝒜k(X,y,B).
Now,ℬ=⋃kcl(ℬk), where cl(ℬk) denotes the closure of the set. Since cl(ℬk) is bounded by hyperplanes, it is convex. It is closed and bounded and hence compact. Choose k∈{1,…,2n}, such that ℬk is non-empty. Since continuous functions are bounded on compact sets, the supremum of lβ, when restricted to cl(ℬk), must be attained at one or more points in cl(ℬk). By Corollary A.3, the supremum (in our case the maximum) of lβ on cl(ℬk) is attained on the whole set or on a union of faces of dimension less than m or at a vertex. Since there are a finite number of sets to consider, the maximum of lβ must occur at a vertex but might occur, for example, on the whole of a set or on a union of faces. This is important to consider when using search algorithms such as the simplex method; as repeated application with different starting points may give a set of solutions which, for example, lie on a line segment. Note the following points.(i)
Assuming thatrX=m, any vertex of the set in ℝm at which lβ attains its maximum must correspond to m data points, possibly more (degeneracy). This is due to the fact that in ℝn, the gradient ∇lz(z;p) points in the direction of the boundary of the corresponding orthant and away from the boundary of ΩB.(ii)
A MLE is not necessarily unique.(iii)
IfrX<m, then we may apply the same reasoning to a subspace of ℝm of dimension rX on which the error mapping has full rank rX.(iv)
Since at the MLE the absolute values of the deviations|zi(β^n)| will all be small, this proof for finite B may be extended to B=∞.Lemma8 is useful but we need to know what happens for more general perturbing functions g. First, we consider the Laplace distribution without perturbation.


#### 3.3.2. The Truncated Laplace Probability Density Function

Forn=1, let
(32)f(z;p)=(pQ(p;1))e(-p|z|)
be the Laplace probability density function (8) with scale parameter p and with support [-1,1]. Then
(33)Q(p;1)=2(1-e-p),f(z;p)=p2(1-e-p)e(-p|z|),
and
(34)log(f(z;p))=log(p2(1-e-p))-p|z|.
The latter is a piecewise affine function in z. It has a maximum value of log(p/(2(1-e-p))) when z=0. More generally, for n>1 and independent deviations (error terms) z1,z2,…,zn,
(35)lz(z;p)=∑i=1nlog(f(zi;p))=∑i=1n(-p|zi|)+nlog(p2(1-e-p)).
Hence, the log-likelihood function is, up to a constant term, a piecewise linear function in the error terms, with a maximum attained when all the errors are zero. However, we must restrict our domain to 𝒜⊂ℝn, or equivalently to ℬ⊂ℝm. The log-likelihood function
(36)lβ(β;p)=-p∑i=1n|(y-Xβ)i|+nlog(p2(1-e-p))
is a piecewise linear function, up to a constant term, in between the hyperplanes Hi-1, Hi0, and Hi+1, i=1,…,n. Let Xj denote the jth column of X, j=1,2…,m. Let
(37)Sp{X}=Sp{X1,X2,…,Xm}⊂ℝn,
denote the span of the columns of X.Nowlz:cl(𝒜k)↦ℝ has a critical point at z(β)∈𝒜k if and only if the gradient ∇lz=(∂l/∂z1,…,∂l/∂zn)T (evaluated at z(β)) is orthogonal to Sp{X}. The columns of X are tangent vectors to 𝒜 at this point. The gradient ∇lz=-p(sgn(z1(β)),sgn(z2(β)),…,sgn(zn(β)))T is constant in the interior of any orthant. If we travel along a straight line path in any orthant, lz either always increases, always decreases, or remains constant. Hence, we will not find an isolated local maximum or minimum in 𝒜k, an open set, for the constrained lz.We need to be aware of the case where𝒜k lies in or very close to a level set of lz. We might need to test for this. This happens when the sign vector (sgn(z1),…,sgn(zn))T is orthogonal to Sp{X}, or nearly so. In the former case, the constrained l is constant on cl(𝒜k). In the latter case, the constrained l will differ very little around the maximum on cl(𝒜k). If this is the case for all the 𝒜k, then the ML values for the coefficients βj will not be sharply defined (will have large variance).


#### 3.3.3. More General Perturbations of the Laplace Probability Density Function

The question is, given a nontrivial perturbing functiong(z;q), does the log-likelihood function attain its maximum at a data point? We have given conditions on g in Lemma 8 which are sufficient to ensure the maximum is attained at a data point. We give a more general criterion in Lemma 10.Assume thatlog(g) is non-linear in any orthant. Otherwise, we can write g in the form of a scaled Laplace distribution and apply Lemma 8. Then l is the sum of an affine function and a non-linear function in any orthant. This affine function is
(38)lp,B=-nlog(Q(p,g1;B))+nlog(p)+∑i=1n(-p|zi|),
the log-likelihood function corresponding to the Laplace distribution. The non-linear function is
(39)lnlin=-nlog(Q(p,g,q;B))+nlog(Q(p,g1;B))+∑i=1nlogg(|zi|;q).
Then l=lp,B+lnlin, and so ∇l=∇lp,B+∇lnlin, where
(40)∇lp,B=-p(sgn(z1),sgn(z2),…,sgn(zn))T,∇lnlin=(dlog(g(z1;q))dz1,dlog(g(z2;q))dz2,…,dlog(g(zn;q))dzn)T.It is possible that there exist orthants𝒪k in which the set 𝒜k(X,y,B) is orthogonal to the gradient ∇lp,B. It is possible that lp,B attains its global maximum (with respect to 𝒜(X,y,B)) on the whole of such an 𝒜k(X,y,B). Hence, it is important to consider the behaviour of lnlin in such orthants.Lemma 10.
Assume bounded support and let(41)γ=sup0≤z≤B{|dg(z;q)/dzg(z;q)|}=sup0≤z≤B{|dlog(g(z;q))dz|},
where sup denotes supremum. Then, if γ<p, the supremum or maximum of l is attained at a data point. In the special case that Sp{X} is orthogonal to (-sgn(z1),-sgn(z2),…,-sgn(zn))T in any orthant, it may be that the supremum is also attained elsewhere.Proof of Lemma10.
Choose𝒜k(X,y,B)=𝒜(X,y,B)∩𝒪k, where k∈{1,2…,2n}, such that the set 𝒜k(X,y,B) is non-empty, choose z∈𝒜k(X,y,B) (an open set relative to 𝒜), and let w(k,X) be the projection of
(42)∇lp,Bp=(-sgn(z1),-sgn(z2),…,-sgn(zn))T
onto Sp{X}. Consider the case w(k,X) is non-zero. Then the affine function lp,B is not constant on 𝒜k(X,y,B). Given any v in 𝒜k(X,y,B), there exists some ϵ>0 and  a map
(43)hp,B:(-ϵ,ϵ)⟶ℝt⟼lp,B(v+tw(k,X))
which is increasing at t=0, since
(44)(dhp,Bdt)t=0=∇lp,BT(d(v+tw(k,X))dt)t=0=∇lp,BTw(k,X)=(w(k,X))Tw(k,X).
Now consider
(45)h:(-ϵ,ϵ)⟶ℝt⟼l(v+tw(k,X)).
Since g(z;q) is continuous and strictly positive on [0,B], γ is finite. We show that if γ<p, dh/dt>0 at t=0, that is, at v. Specifically,
(46)dhdt=-p∑i=1n(sgn(vi+twi(k,X))wi(k,X))+∑i=1n(dg(|zi|;q)/d|zi|)h(t)g(|vi+twi(k,X)|;q)×(sgn(vi+twi(k,X)))wi(k,X)≥-p∑i=1n(sgn(vi+twi(k,X))wi(k,X))+γ∑i=1n(sgn(vi+twi(k,X))wi(k,X))=(1-γp)(w(k,X))Tw(k,X)>0
if γ<p and w(k,X) is non-zero. Then the supremum of l, constrained to cl(𝒜k(X,y,B)), must be attained on the boundary of 𝒜k(X,y,B), relative to 𝒜, and this must be at a data point due to the direction of ∇l.
In the limiting case, wherew(k,X)=0, but we still have γ<p, we find that the supremum of l, constrained to cl(𝒜k(X,y,B)) (relative to 𝒜), must be attained at a data point but may be attained at other points in 𝒜k(X,y,B) as well. Assume this is not the case, that is, there exists no data point at which l attains its supremum when constrained to cl(𝒜k(X,y,B)). We find a contradiction as follows. Assume there exists some v∈𝒜k(X,y,B), an open set (relative to 𝒜), such that
(47)lz(v)-sup{lz(z):z∈bd(𝒪k)∩cl(𝒜k(X,y,B))}=τ>0,
where bd denotes boundary. Imagine perturbing the columns of X slightly (and continuously) to obtain X′ so that lβ, a continuous function, conditional on X, changes by at most τ/3, that is,
(48)|lβ(β;p;X,y)-lβ(β;p;X′,y)|<τ3
for all β∈ℬk(X,y,B)∩ℬk(X′,y,B), and so that now w(k,X′)≠0. We also require the perturbation small enough that
(49)|sup{lz(z;p):z∈bd(𝒪k)∩cl(𝒜k(X,y,B))}(𝒜k(X′,y,B))-sup{lz(z;p):z∈bd(𝒪k)∩cl(𝒜k(X′,y,B))}|<τ3.
Now v=y-Xβv for some unique βv. Let v′=y-X′βv. If the perturbation of X is small enough, then βv∈ℬk(X′,y,B) and so v′∈𝒜k(X′,y,B).
Now,|lz(v′)-lz(v)|<τ/3. Hence,
(50)lz(v′)-sup{lz(z):z∈bd(𝒪k)∩cl(𝒜k(X′,y,B))}>τ3>0,
which is not possible since w(k,X′)≠0. Hence, we have a contradiction. Now, for completeness, assume there exists some r∈bd(𝒜k(X,y,B)), r not a data point (that is some ri=±B) such that
(51)lz(r)-sup{lz(z):z∈bd(𝒪k)∩cl(𝒜k(X,y,B))}=τ>0.
Then there exists v∈𝒜k(X,y,B), such that
(52)lz(r)-lz(v)<τ6.
Hence,
(53)lz(v)-sup{lz(z):z∈bd(𝒪k)∩cl(𝒜k(X,y,B))}>5τ6>0.
The contradiction follows as above. Hence, we have proved the lemma.


#### 3.3.4. An Amended Laplace Probability Density Function with Added Kurtosis

We may apply Lemma10 to show that, for the amended Laplace probability density function with added kurtosis (Example 5) and for realistic values of p and q, the maximum of the log-likelihood function must be attained at a data point. Here, g(z;q)=g2(z;q)=1+q(z3-3z)=1+qH3(z) on [0,1], that is, g2(|z|;q)=1+q(H3(|z|)) on [-1,1]. In n dimensions,
(54)l(β;p;X,y)=-nlog(Q(p,g2,q;1))+nlog(p)+∑i=1n(-p|yi-(Xβ)i|)+∑i=1nlog(1+qH3(|yi-(Xβ)i|))=-nlog(Q(p,g2,q;1))+nlog(p)+∑i=1n(-p|yi-(Xβ)i|)+∑i=1nlog(1+q(|yi-(Xβ)i|3-3|yi-(Xβ)i||yi-(Xβ)i|3)q(|yi-(Xβ)i|3).
Consider, for z∈[0,1],
(55)dlog(g2(z;q))dz=d(log(1+q(z3-3z)))dz=3q(z2-1)1+q(z3-3z),d2log(g2(z;q))dz2=d2(log(1+q(z3-3z)))dz2=(3q)(2z-q(3+z4))(1+q(z3-3z))2.
When 0<q<0.5, d2log(g2(z;q))/dz2 is negative when z=0 and positive when z=1, and so the continuous monotonic function log(g2) has a point of inflection, where d2log(g2(z;q))/dz2=0 in (0,1), at say ω(q). Hence, log(g2) and l are both concave and convex on [0,1]. When q is small, the concavity occurs close to the origin and the functions are convex on most of (0,1). Moreover, dlog(g2(z;q))/dz is always negative on (0,1), and its limit as z→∞ is also negative, so we will not get isolated local maxima in one dimension. The function log(g2(z;q)) and its first and second derivatives are plotted in Figures 1, 2, and 3, respectively; setting q=0.025, a typical value. Reading from Figure 2, the upper bound γ is approximately 0.075, when q=0.025. Since typically p≥3, the criterion (γ<p) for Lemma 10 is satisfied.Figure 1
The non-linear part of the log-likelihood function,n=1; log(g(z;q)), q=0.025.Figure 2
The first derivative oflog(g(z;q)), q=0.025.Figure 3
The second derivative oflog(g(z;q)), q=0.025.


#### 3.3.5. Non-Increasing Perturbations Both Concave and Convex

In certain situations we might find the criteria for both Lemmas8 and 10 are not satisfied but that logg(z;q) is non-increasing on [0,B] and convex on [w(q),B], where w(q) is small. Since a sum of convex functions is convex (see [8]), l will be convex wherever logg(|zi|;q) is convex for i=1,2,…,n. For example, for g=g2 and B=1, we can prove a partial result as follows. By restricting to the domain
(56)ℬ(X,y,[ω(q),1])=⋂i=1n{β∈ℝm:ω(q)≤|zi(β)|≤1}⊂ℬ(X,y,1)
or equivalently by substituting for each set 𝒜k(X,y,1) the subset
(57){z∈∩𝒜k(X,y,1):ω(q)<|zi|<1,i=1,…,n}⊂𝒜k(X,y,1),
and applying convex function theory as for Lemma 8 we can prove that there exists β1∈ℬ(X,y,[ω(q),1]) at which l attains its maximum, and there exists at least one index i such that zi(β1)=ω(q). Hence, there exists some β2∈ℬ(X,y,1)=ℬ(X,y,[0,1]) at which l attains its maximum, and there exists at least one index i such that zi(β2)≤ω(q). In other words, there exists some point at which l has a maximum at which at least one of the errors is small. It makes sense to search for the maxima of l near or at vertices.Forg=g2, we may apply Lemma 10 and so do not need this partial result, but for a perturbation similar in shape to g2, non-increasing everywhere on [0,1], with logg(z;q) convex everywhere on (ω(g;q),1], concave everywhere on [0,ω(g;q)), for some small positive ω(g;q), and with steep slope at the point of inflection (z=ω(g;q)), such analysis would be useful.


## 3.1. The Log-Likelihood Function

Let(16)f(z,p)
be a probability density function with parameter vector p, as described in Section 2. Let z1,…,zn be a sequence of independent and identically distributed deviations with joint probability density function
(17)f(z1,…,zn;p)=f(z1(β),…,zn(β);p)=Πi=1nf(zi(β);p).
This is also the likelihood function
(18)Lz(z;p)=Lβ(β;p;X,y)=L(z(β);p)=f(z(β);p),
which may be regarded as a function of z or β; here, the subscript reflects our point of view. We use the log-likelihood function l in the estimation of β, where, using various notation,
(19)lz(z;p)=loge(Lz(z;p))=lβ(β;p;X,y)=l(z(β);p).
Substituting measured values of yi and known inputs xij into lβ, we obtain a function of β and p. The parameters p are assumed known, but if not, may be estimated separately. In our biological application, they are estimated by MLE and are assumed fixed for a particular measuring process. Hence, we have a function of β, the coefficients of our linear model.Our aim is to find a maximum likelihood estimator (MLE) denotedβ^n∈ℝm, that is, some point at which l attains its maximum value. The subscript n corresponds to the number of deviations. Now l is a continuous function. If B is finite, it has compact support in ℝm. Since a continuous function on a compact set attains its supremum, the existence of a MLE for lβ is guaranteed. Even if B=∞, we may consider truncations with finite bounds Bk=k, k=1,2,…. Since l is maximized when the zi are small, truncating f to [-k,k] for k large enough will not affect the set of points at which l attains its maximum. We show in Section 3.3 that a MLE β^n is not necessarily unique.Case 1.
Ifg(z;q)=g1(z)=1 and so f is the Laplace probability density function with parameter p as in (6) or (8), then, for β∈ℝm such that every zi(β) is in the support of f,
(20)lβ(β;p;X,y)=lβ,n,m(β;p;Xn,m,y)=-nlog(Q)+nlog(p)+∑i=1n(-p|zi(β)|)=-nlog(Q)+nlog(p)+∑i=1n(-p|yi-(Xβ)i|),
where Q=2 if the support of f is ℝ. If the support of f is [-B,B] then
(21)Q=Q(p,g1;B)=2(1-e-pB).
The theory of LAE regression (corresponding to MLE using Laplace distributions without modification as response functions) may be found in various texts, for example, [3]. Here, it is proved that there exists a MLE β^n corresponding to at least rX zero errors. We are concerned with the extension of these ideas to the case of perturbed and truncated Laplace response functions. For the truncated Laplace distribution, we prove that there exists a MLE β^n corresponding to at least rX zero errors. Consider the following more general case.Case 2.
Iff is a perturbed Laplace probability density function with perturbing function g(z;q) and bounded support as in (10), then for β∈ℝm such that |zi(β)|≤B, i=1,2,…,n,
(22)lβ(β;p;X,y)=lβ,n,m(β;p;Xn,m,y)=-nlog(Q(p,g,q;B))+nlog(p)+∑i=1n(-p|zi(β)|)+∑i=1nlog(g(|zi(β)|;q))=-nlog(Q(p,g,q;B))+nlog(p)+∑i=1n(-p|yi-(Xβ)i|)+∑i=1nlog(g(|yi-(Xβ)i|;q)).
Note since g(zi;q) is strictly positive on [0,B], so is f(zi;p) and so log(f(zi;p)) is well-defined on [-B,B], i=1,2,…,n.
In Section3.3, we show that if the perturbing function g is such that logg(z;q) is convex and non-increasing on [0,B], there exists a MLE corresponding to at least rX data points. We give an upper bound on |dlogg(z;q)/dz| on [0,B] which, if not exceeded, ensures that there exists a MLE corresponding to at least one data point. We apply these results to our motivating example, the Laplace distribution with added kurtosis, described below.Example 7.
Ifg(z;q)=g2(z;q)=1+qH3(z) and so f is a Laplace probability density function with scale parameter p, with added kurtosis and bounded support as in (13), then for β∈ℝm such that |zi(β)|≤1, i=1,2,…,n,
(23)lβ(β;p,q;X,y)=lβ,n,m(β;p;Xn,m,y)=-nlog(Q(p,g2,q;1))+nlog(p)+∑i=1n(-p|zi(β)|)+∑i=1nlog(1+qH3(|zi(β)|))=-nlog(Q(p,g2,q;1))+nlog(p)+∑i=1n(-p|yi-(Xβ)i|)+∑i=1nlog(1+qH3(|yi-(Xβ)i|)).
Now lβ is not differentiable in the classical sense with respect to the linear model coefficients βj when any zi=0. Hence, we cannot assume that lβ is differentiable at a MLE. This paper addresses this issue firstly by proposing a non-gradient method of coefficient estimation and secondly (in Section 4) by using generalized functions to calculate statistical estimates including estimates of standard errors. In Section 5, we discuss the asymptotic theory of our MLE.


## 3.2. Coefficient Estimation Dealing with Abrupt Changes in Gradient

AlthoughLβ and lβ are continuous functions of β, their first derivatives are not. Consider the geometry of the coefficient space ℝm, where β∈ℝm. For each index i, since the response function distribution is defined in terms of absolute values, we can find a hyperplane Hi0 in ℝm on which L=Lβ and l=lβ are not differentiable, defined by setting zi=0. Let xiT=(xi1,xi2,…,xim), the ith row of X, it is never the zero vector since xi1=1, i=1,2,…,n. Choose w∈ℝm so that xiTw=0. Then
(24)β=w+yixixiTxi
yields zi=0. Let Hi0 be the set of all such β. For example, for m=2, for each error term there is a line in ℝ2 on which Lβ and lβ have a sharp ridge. By inspection of the geometry, we would expect the values of β which maximize lβ to be either on the union of the hyperplanes or very close to intersections of the hyperplanes Hi0, i=1,2,…,n. Imagine searching in β-space near the hyperplanes Hi0. Even if Lβ has a local maximum near but not on the union of the hyperplanes, it would be difficult to use a method based on the gradient of either Lβ or lβ, since the gradient changes sharply whenever we cross one of the Hi0. The simplex method of coefficient estimation, which does not require any partial derivatives, suits this geometry.


## 3.3. The Maximum Likelihood Estimator Corresponds to a Data Point

### 3.3.1. Convex and Non-Increasing Perturbations of the Laplace Probability Density Function

Consider the probability density function(25)f(z;p,q)=(pQ(p,g,q;B))e(-p|z|)g(|z|;q),
with support [-B,B], for some finite B>0 as described in Section 2.2 (Example 3). Recall we assume that there exists some ϵ>0 such that g(z;q) is 𝒞3 in z on (-ϵ,B+ϵ) and that g(z;q)>0 on [0,B], for fixed parameter q. Let
(26)ΩB={z∈ℝn:|zi|≤B,i=1,…,n}=[-B,B]n.
We consider the log-likelihood function
(27)lz(z;p):ΩB⟶ℝ
(conditional on p) with its domain restricted to
(28)𝒜={z(β):β∈ℝm}∩ΩB={y-Xβ:β∈ℝm}∩ΩB,
that is, constrained to 𝒜 or equivalently, consider the log-likelihood function
(29)lβ(β;p;X,y):ℬ⟶ℝ,
constrained to ℬ, where ℬ⊂ℝm is defined as {β∈ℝm:z(β)∈ΩB}. Note that 𝒜=𝒜(X,y,B) and similarly ℬ=ℬ(X,y,B). Also, if the function Z is defined by
(30)Z:ℝm⟶ℝn,β⟼z(β)=y-Xβ,
then Z(ℬ)=𝒜.Lemma 8.
Iflog(g):(-ϵ,B+ϵ)→ℝ is convex and non-increasing (i.e., dlogg(z;q)/dz≤0 on [0,B], then there exists a maximum of lβ:ℬ→ℝ corresponding to at least rX data points. That is, there exists β^n∈ℬ⊂ℝm such that the constrained lβ attains its maximum at β^n, and there exists at least rX indices ij∈{1,2,…,n} such that zij(β^n)=0, j=1,2,…,rX.Corollary 9.
Letf(z;p) be the Laplace probability density function (8) with support [-B,B] and parameter p. Then there exists a maximum of lβ:ℬ→ℝ corresponding to at least rX data points.Proof of Corollary9.
Letg be the constant map g1(z)=1; then log(g(z))=0, and hence log(g(z)) is convex and non-increasing on [0,B]. Corollary 9 follows from Lemma 8.Note that Corollary9 could be proved directly by linear programming theory. Linear programming has been applied to the problem of minimising the sum of absolute errors in various applications (see [6], a survey paper, and also [7]). The convex analysis results we require are in Appendix A.Proof of Lemma8.
To begin, assume thatX has full rank m, recall n≥m and that f has bounded support. Then 𝒜 is a compact convex subset of an m-dimensional affine subspace of ℝn. Since the mapping Z is linear and has full rank m, the inverse image ℬ=Z-1(𝒜) is a compact convex subset of ℝm. We partition ℬ⊂ℝm, which is the support of Lβ, into a finite collection of compact convex sets, so that, on each subset, lβ is convex.
LetHiδ be the hyperplane in ℝm corresponding to the error term zi(β)=δ. Then Hi-B and HiB are the hyperplanes in ℝm corresponding to errors zi=-B and zi=B, respectively. It follows that the log-likelihood function
(31)lβ(β;p;X,y)=-nlog(Q(p,g,q;B))+nlog(p)+∑i=1n(-p|zi(β)|)+∑i=1nlog(g(|zi(β)|;q))=-nlog(Q(p,g,q;B))+nlog(p)+∑i=1n(-p|yi-(Xβ)i|)+∑i=1nlog(g(|yi-(Xβ)i|;q))
is a convex function in between the the hyperplanes Hi-B, Hi0, and HiB, i=1,…,n. These 3n hyperplanes divide the domain ℬ in ℝm into at most 2n open sets bounded by (but not intersecting) the hyperplanes. Each such open set (and hence its closure) may be labelled by a set of n signs. For any β∈ℝm such that 0<|zi(β)|<B, i=1,…,n; (sgn(z1(β)),sgn(z2(β)),…,sgn(zn(β))) labels the open set containing β.
Next, considerℝn as the union of its orthants, which we denote 𝒪k, k=1,…,2n. We assume the orthants are closed sets. For example, the non-negative orthant is {z∈ℝn:zi≥0,i=1,…,n}. In the interior of any 𝒪k, the sign of zi does not change, i=1,…,n. Relabel the open subsets ℬk=ℬk(X,y,B), where k∈{1,2,…,2n}, so that Z(ℬk)⊂int(𝒪k). Let 𝒜k=Z(ℬk)=𝒜k(X,y,B).
Now,ℬ=⋃kcl(ℬk), where cl(ℬk) denotes the closure of the set. Since cl(ℬk) is bounded by hyperplanes, it is convex. It is closed and bounded and hence compact. Choose k∈{1,…,2n}, such that ℬk is non-empty. Since continuous functions are bounded on compact sets, the supremum of lβ, when restricted to cl(ℬk), must be attained at one or more points in cl(ℬk). By Corollary A.3, the supremum (in our case the maximum) of lβ on cl(ℬk) is attained on the whole set or on a union of faces of dimension less than m or at a vertex. Since there are a finite number of sets to consider, the maximum of lβ must occur at a vertex but might occur, for example, on the whole of a set or on a union of faces. This is important to consider when using search algorithms such as the simplex method; as repeated application with different starting points may give a set of solutions which, for example, lie on a line segment. Note the following points.(i)
Assuming thatrX=m, any vertex of the set in ℝm at which lβ attains its maximum must correspond to m data points, possibly more (degeneracy). This is due to the fact that in ℝn, the gradient ∇lz(z;p) points in the direction of the boundary of the corresponding orthant and away from the boundary of ΩB.(ii)
A MLE is not necessarily unique.(iii)
IfrX<m, then we may apply the same reasoning to a subspace of ℝm of dimension rX on which the error mapping has full rank rX.(iv)
Since at the MLE the absolute values of the deviations|zi(β^n)| will all be small, this proof for finite B may be extended to B=∞.Lemma8 is useful but we need to know what happens for more general perturbing functions g. First, we consider the Laplace distribution without perturbation.


### 3.3.2. The Truncated Laplace Probability Density Function

Forn=1, let
(32)f(z;p)=(pQ(p;1))e(-p|z|)
be the Laplace probability density function (8) with scale parameter p and with support [-1,1]. Then
(33)Q(p;1)=2(1-e-p),f(z;p)=p2(1-e-p)e(-p|z|),
and
(34)log(f(z;p))=log(p2(1-e-p))-p|z|.
The latter is a piecewise affine function in z. It has a maximum value of log(p/(2(1-e-p))) when z=0. More generally, for n>1 and independent deviations (error terms) z1,z2,…,zn,
(35)lz(z;p)=∑i=1nlog(f(zi;p))=∑i=1n(-p|zi|)+nlog(p2(1-e-p)).
Hence, the log-likelihood function is, up to a constant term, a piecewise linear function in the error terms, with a maximum attained when all the errors are zero. However, we must restrict our domain to 𝒜⊂ℝn, or equivalently to ℬ⊂ℝm. The log-likelihood function
(36)lβ(β;p)=-p∑i=1n|(y-Xβ)i|+nlog(p2(1-e-p))
is a piecewise linear function, up to a constant term, in between the hyperplanes Hi-1, Hi0, and Hi+1, i=1,…,n. Let Xj denote the jth column of X, j=1,2…,m. Let
(37)Sp{X}=Sp{X1,X2,…,Xm}⊂ℝn,
denote the span of the columns of X.Nowlz:cl(𝒜k)↦ℝ has a critical point at z(β)∈𝒜k if and only if the gradient ∇lz=(∂l/∂z1,…,∂l/∂zn)T (evaluated at z(β)) is orthogonal to Sp{X}. The columns of X are tangent vectors to 𝒜 at this point. The gradient ∇lz=-p(sgn(z1(β)),sgn(z2(β)),…,sgn(zn(β)))T is constant in the interior of any orthant. If we travel along a straight line path in any orthant, lz either always increases, always decreases, or remains constant. Hence, we will not find an isolated local maximum or minimum in 𝒜k, an open set, for the constrained lz.We need to be aware of the case where𝒜k lies in or very close to a level set of lz. We might need to test for this. This happens when the sign vector (sgn(z1),…,sgn(zn))T is orthogonal to Sp{X}, or nearly so. In the former case, the constrained l is constant on cl(𝒜k). In the latter case, the constrained l will differ very little around the maximum on cl(𝒜k). If this is the case for all the 𝒜k, then the ML values for the coefficients βj will not be sharply defined (will have large variance).


### 3.3.3. More General Perturbations of the Laplace Probability Density Function

The question is, given a nontrivial perturbing functiong(z;q), does the log-likelihood function attain its maximum at a data point? We have given conditions on g in Lemma 8 which are sufficient to ensure the maximum is attained at a data point. We give a more general criterion in Lemma 10.Assume thatlog(g) is non-linear in any orthant. Otherwise, we can write g in the form of a scaled Laplace distribution and apply Lemma 8. Then l is the sum of an affine function and a non-linear function in any orthant. This affine function is
(38)lp,B=-nlog(Q(p,g1;B))+nlog(p)+∑i=1n(-p|zi|),
the log-likelihood function corresponding to the Laplace distribution. The non-linear function is
(39)lnlin=-nlog(Q(p,g,q;B))+nlog(Q(p,g1;B))+∑i=1nlogg(|zi|;q).
Then l=lp,B+lnlin, and so ∇l=∇lp,B+∇lnlin, where
(40)∇lp,B=-p(sgn(z1),sgn(z2),…,sgn(zn))T,∇lnlin=(dlog(g(z1;q))dz1,dlog(g(z2;q))dz2,…,dlog(g(zn;q))dzn)T.It is possible that there exist orthants𝒪k in which the set 𝒜k(X,y,B) is orthogonal to the gradient ∇lp,B. It is possible that lp,B attains its global maximum (with respect to 𝒜(X,y,B)) on the whole of such an 𝒜k(X,y,B). Hence, it is important to consider the behaviour of lnlin in such orthants.Lemma 10.
Assume bounded support and let(41)γ=sup0≤z≤B{|dg(z;q)/dzg(z;q)|}=sup0≤z≤B{|dlog(g(z;q))dz|},
where sup denotes supremum. Then, if γ<p, the supremum or maximum of l is attained at a data point. In the special case that Sp{X} is orthogonal to (-sgn(z1),-sgn(z2),…,-sgn(zn))T in any orthant, it may be that the supremum is also attained elsewhere.Proof of Lemma10.
Choose𝒜k(X,y,B)=𝒜(X,y,B)∩𝒪k, where k∈{1,2…,2n}, such that the set 𝒜k(X,y,B) is non-empty, choose z∈𝒜k(X,y,B) (an open set relative to 𝒜), and let w(k,X) be the projection of
(42)∇lp,Bp=(-sgn(z1),-sgn(z2),…,-sgn(zn))T
onto Sp{X}. Consider the case w(k,X) is non-zero. Then the affine function lp,B is not constant on 𝒜k(X,y,B). Given any v in 𝒜k(X,y,B), there exists some ϵ>0 and  a map
(43)hp,B:(-ϵ,ϵ)⟶ℝt⟼lp,B(v+tw(k,X))
which is increasing at t=0, since
(44)(dhp,Bdt)t=0=∇lp,BT(d(v+tw(k,X))dt)t=0=∇lp,BTw(k,X)=(w(k,X))Tw(k,X).
Now consider
(45)h:(-ϵ,ϵ)⟶ℝt⟼l(v+tw(k,X)).
Since g(z;q) is continuous and strictly positive on [0,B], γ is finite. We show that if γ<p, dh/dt>0 at t=0, that is, at v. Specifically,
(46)dhdt=-p∑i=1n(sgn(vi+twi(k,X))wi(k,X))+∑i=1n(dg(|zi|;q)/d|zi|)h(t)g(|vi+twi(k,X)|;q)×(sgn(vi+twi(k,X)))wi(k,X)≥-p∑i=1n(sgn(vi+twi(k,X))wi(k,X))+γ∑i=1n(sgn(vi+twi(k,X))wi(k,X))=(1-γp)(w(k,X))Tw(k,X)>0
if γ<p and w(k,X) is non-zero. Then the supremum of l, constrained to cl(𝒜k(X,y,B)), must be attained on the boundary of 𝒜k(X,y,B), relative to 𝒜, and this must be at a data point due to the direction of ∇l.
In the limiting case, wherew(k,X)=0, but we still have γ<p, we find that the supremum of l, constrained to cl(𝒜k(X,y,B)) (relative to 𝒜), must be attained at a data point but may be attained at other points in 𝒜k(X,y,B) as well. Assume this is not the case, that is, there exists no data point at which l attains its supremum when constrained to cl(𝒜k(X,y,B)). We find a contradiction as follows. Assume there exists some v∈𝒜k(X,y,B), an open set (relative to 𝒜), such that
(47)lz(v)-sup{lz(z):z∈bd(𝒪k)∩cl(𝒜k(X,y,B))}=τ>0,
where bd denotes boundary. Imagine perturbing the columns of X slightly (and continuously) to obtain X′ so that lβ, a continuous function, conditional on X, changes by at most τ/3, that is,
(48)|lβ(β;p;X,y)-lβ(β;p;X′,y)|<τ3
for all β∈ℬk(X,y,B)∩ℬk(X′,y,B), and so that now w(k,X′)≠0. We also require the perturbation small enough that
(49)|sup{lz(z;p):z∈bd(𝒪k)∩cl(𝒜k(X,y,B))}(𝒜k(X′,y,B))-sup{lz(z;p):z∈bd(𝒪k)∩cl(𝒜k(X′,y,B))}|<τ3.
Now v=y-Xβv for some unique βv. Let v′=y-X′βv. If the perturbation of X is small enough, then βv∈ℬk(X′,y,B) and so v′∈𝒜k(X′,y,B).
Now,|lz(v′)-lz(v)|<τ/3. Hence,
(50)lz(v′)-sup{lz(z):z∈bd(𝒪k)∩cl(𝒜k(X′,y,B))}>τ3>0,
which is not possible since w(k,X′)≠0. Hence, we have a contradiction. Now, for completeness, assume there exists some r∈bd(𝒜k(X,y,B)), r not a data point (that is some ri=±B) such that
(51)lz(r)-sup{lz(z):z∈bd(𝒪k)∩cl(𝒜k(X,y,B))}=τ>0.
Then there exists v∈𝒜k(X,y,B), such that
(52)lz(r)-lz(v)<τ6.
Hence,
(53)lz(v)-sup{lz(z):z∈bd(𝒪k)∩cl(𝒜k(X,y,B))}>5τ6>0.
The contradiction follows as above. Hence, we have proved the lemma.


### 3.3.4. An Amended Laplace Probability Density Function with Added Kurtosis

We may apply Lemma10 to show that, for the amended Laplace probability density function with added kurtosis (Example 5) and for realistic values of p and q, the maximum of the log-likelihood function must be attained at a data point. Here, g(z;q)=g2(z;q)=1+q(z3-3z)=1+qH3(z) on [0,1], that is, g2(|z|;q)=1+q(H3(|z|)) on [-1,1]. In n dimensions,
(54)l(β;p;X,y)=-nlog(Q(p,g2,q;1))+nlog(p)+∑i=1n(-p|yi-(Xβ)i|)+∑i=1nlog(1+qH3(|yi-(Xβ)i|))=-nlog(Q(p,g2,q;1))+nlog(p)+∑i=1n(-p|yi-(Xβ)i|)+∑i=1nlog(1+q(|yi-(Xβ)i|3-3|yi-(Xβ)i||yi-(Xβ)i|3)q(|yi-(Xβ)i|3).
Consider, for z∈[0,1],
(55)dlog(g2(z;q))dz=d(log(1+q(z3-3z)))dz=3q(z2-1)1+q(z3-3z),d2log(g2(z;q))dz2=d2(log(1+q(z3-3z)))dz2=(3q)(2z-q(3+z4))(1+q(z3-3z))2.
When 0<q<0.5, d2log(g2(z;q))/dz2 is negative when z=0 and positive when z=1, and so the continuous monotonic function log(g2) has a point of inflection, where d2log(g2(z;q))/dz2=0 in (0,1), at say ω(q). Hence, log(g2) and l are both concave and convex on [0,1]. When q is small, the concavity occurs close to the origin and the functions are convex on most of (0,1). Moreover, dlog(g2(z;q))/dz is always negative on (0,1), and its limit as z→∞ is also negative, so we will not get isolated local maxima in one dimension. The function log(g2(z;q)) and its first and second derivatives are plotted in Figures 1, 2, and 3, respectively; setting q=0.025, a typical value. Reading from Figure 2, the upper bound γ is approximately 0.075, when q=0.025. Since typically p≥3, the criterion (γ<p) for Lemma 10 is satisfied.Figure 1
The non-linear part of the log-likelihood function,n=1; log(g(z;q)), q=0.025.Figure 2
The first derivative oflog(g(z;q)), q=0.025.Figure 3
The second derivative oflog(g(z;q)), q=0.025.


### 3.3.5. Non-Increasing Perturbations Both Concave and Convex

In certain situations we might find the criteria for both Lemmas8 and 10 are not satisfied but that logg(z;q) is non-increasing on [0,B] and convex on [w(q),B], where w(q) is small. Since a sum of convex functions is convex (see [8]), l will be convex wherever logg(|zi|;q) is convex for i=1,2,…,n. For example, for g=g2 and B=1, we can prove a partial result as follows. By restricting to the domain
(56)ℬ(X,y,[ω(q),1])=⋂i=1n{β∈ℝm:ω(q)≤|zi(β)|≤1}⊂ℬ(X,y,1)
or equivalently by substituting for each set 𝒜k(X,y,1) the subset
(57){z∈∩𝒜k(X,y,1):ω(q)<|zi|<1,i=1,…,n}⊂𝒜k(X,y,1),
and applying convex function theory as for Lemma 8 we can prove that there exists β1∈ℬ(X,y,[ω(q),1]) at which l attains its maximum, and there exists at least one index i such that zi(β1)=ω(q). Hence, there exists some β2∈ℬ(X,y,1)=ℬ(X,y,[0,1]) at which l attains its maximum, and there exists at least one index i such that zi(β2)≤ω(q). In other words, there exists some point at which l has a maximum at which at least one of the errors is small. It makes sense to search for the maxima of l near or at vertices.Forg=g2, we may apply Lemma 10 and so do not need this partial result, but for a perturbation similar in shape to g2, non-increasing everywhere on [0,1], with logg(z;q) convex everywhere on (ω(g;q),1], concave everywhere on [0,ω(g;q)), for some small positive ω(g;q), and with steep slope at the point of inflection (z=ω(g;q)), such analysis would be useful.


## 3.3.1. Convex and Non-Increasing Perturbations of the Laplace Probability Density Function

Consider the probability density function(25)f(z;p,q)=(pQ(p,g,q;B))e(-p|z|)g(|z|;q),
with support [-B,B], for some finite B>0 as described in Section 2.2 (Example 3). Recall we assume that there exists some ϵ>0 such that g(z;q) is 𝒞3 in z on (-ϵ,B+ϵ) and that g(z;q)>0 on [0,B], for fixed parameter q. Let
(26)ΩB={z∈ℝn:|zi|≤B,i=1,…,n}=[-B,B]n.
We consider the log-likelihood function
(27)lz(z;p):ΩB⟶ℝ
(conditional on p) with its domain restricted to
(28)𝒜={z(β):β∈ℝm}∩ΩB={y-Xβ:β∈ℝm}∩ΩB,
that is, constrained to 𝒜 or equivalently, consider the log-likelihood function
(29)lβ(β;p;X,y):ℬ⟶ℝ,
constrained to ℬ, where ℬ⊂ℝm is defined as {β∈ℝm:z(β)∈ΩB}. Note that 𝒜=𝒜(X,y,B) and similarly ℬ=ℬ(X,y,B). Also, if the function Z is defined by
(30)Z:ℝm⟶ℝn,β⟼z(β)=y-Xβ,
then Z(ℬ)=𝒜.Lemma 8.
Iflog(g):(-ϵ,B+ϵ)→ℝ is convex and non-increasing (i.e., dlogg(z;q)/dz≤0 on [0,B], then there exists a maximum of lβ:ℬ→ℝ corresponding to at least rX data points. That is, there exists β^n∈ℬ⊂ℝm such that the constrained lβ attains its maximum at β^n, and there exists at least rX indices ij∈{1,2,…,n} such that zij(β^n)=0, j=1,2,…,rX.Corollary 9.
Letf(z;p) be the Laplace probability density function (8) with support [-B,B] and parameter p. Then there exists a maximum of lβ:ℬ→ℝ corresponding to at least rX data points.Proof of Corollary9.
Letg be the constant map g1(z)=1; then log(g(z))=0, and hence log(g(z)) is convex and non-increasing on [0,B]. Corollary 9 follows from Lemma 8.Note that Corollary9 could be proved directly by linear programming theory. Linear programming has been applied to the problem of minimising the sum of absolute errors in various applications (see [6], a survey paper, and also [7]). The convex analysis results we require are in Appendix A.Proof of Lemma8.
To begin, assume thatX has full rank m, recall n≥m and that f has bounded support. Then 𝒜 is a compact convex subset of an m-dimensional affine subspace of ℝn. Since the mapping Z is linear and has full rank m, the inverse image ℬ=Z-1(𝒜) is a compact convex subset of ℝm. We partition ℬ⊂ℝm, which is the support of Lβ, into a finite collection of compact convex sets, so that, on each subset, lβ is convex.
LetHiδ be the hyperplane in ℝm corresponding to the error term zi(β)=δ. Then Hi-B and HiB are the hyperplanes in ℝm corresponding to errors zi=-B and zi=B, respectively. It follows that the log-likelihood function
(31)lβ(β;p;X,y)=-nlog(Q(p,g,q;B))+nlog(p)+∑i=1n(-p|zi(β)|)+∑i=1nlog(g(|zi(β)|;q))=-nlog(Q(p,g,q;B))+nlog(p)+∑i=1n(-p|yi-(Xβ)i|)+∑i=1nlog(g(|yi-(Xβ)i|;q))
is a convex function in between the the hyperplanes Hi-B, Hi0, and HiB, i=1,…,n. These 3n hyperplanes divide the domain ℬ in ℝm into at most 2n open sets bounded by (but not intersecting) the hyperplanes. Each such open set (and hence its closure) may be labelled by a set of n signs. For any β∈ℝm such that 0<|zi(β)|<B, i=1,…,n; (sgn(z1(β)),sgn(z2(β)),…,sgn(zn(β))) labels the open set containing β.
Next, considerℝn as the union of its orthants, which we denote 𝒪k, k=1,…,2n. We assume the orthants are closed sets. For example, the non-negative orthant is {z∈ℝn:zi≥0,i=1,…,n}. In the interior of any 𝒪k, the sign of zi does not change, i=1,…,n. Relabel the open subsets ℬk=ℬk(X,y,B), where k∈{1,2,…,2n}, so that Z(ℬk)⊂int(𝒪k). Let 𝒜k=Z(ℬk)=𝒜k(X,y,B).
Now,ℬ=⋃kcl(ℬk), where cl(ℬk) denotes the closure of the set. Since cl(ℬk) is bounded by hyperplanes, it is convex. It is closed and bounded and hence compact. Choose k∈{1,…,2n}, such that ℬk is non-empty. Since continuous functions are bounded on compact sets, the supremum of lβ, when restricted to cl(ℬk), must be attained at one or more points in cl(ℬk). By Corollary A.3, the supremum (in our case the maximum) of lβ on cl(ℬk) is attained on the whole set or on a union of faces of dimension less than m or at a vertex. Since there are a finite number of sets to consider, the maximum of lβ must occur at a vertex but might occur, for example, on the whole of a set or on a union of faces. This is important to consider when using search algorithms such as the simplex method; as repeated application with different starting points may give a set of solutions which, for example, lie on a line segment. Note the following points.(i)
Assuming thatrX=m, any vertex of the set in ℝm at which lβ attains its maximum must correspond to m data points, possibly more (degeneracy). This is due to the fact that in ℝn, the gradient ∇lz(z;p) points in the direction of the boundary of the corresponding orthant and away from the boundary of ΩB.(ii)
A MLE is not necessarily unique.(iii)
IfrX<m, then we may apply the same reasoning to a subspace of ℝm of dimension rX on which the error mapping has full rank rX.(iv)
Since at the MLE the absolute values of the deviations|zi(β^n)| will all be small, this proof for finite B may be extended to B=∞.Lemma8 is useful but we need to know what happens for more general perturbing functions g. First, we consider the Laplace distribution without perturbation.


## 3.3.2. The Truncated Laplace Probability Density Function

Forn=1, let
(32)f(z;p)=(pQ(p;1))e(-p|z|)
be the Laplace probability density function (8) with scale parameter p and with support [-1,1]. Then
(33)Q(p;1)=2(1-e-p),f(z;p)=p2(1-e-p)e(-p|z|),
and
(34)log(f(z;p))=log(p2(1-e-p))-p|z|.
The latter is a piecewise affine function in z. It has a maximum value of log(p/(2(1-e-p))) when z=0. More generally, for n>1 and independent deviations (error terms) z1,z2,…,zn,
(35)lz(z;p)=∑i=1nlog(f(zi;p))=∑i=1n(-p|zi|)+nlog(p2(1-e-p)).
Hence, the log-likelihood function is, up to a constant term, a piecewise linear function in the error terms, with a maximum attained when all the errors are zero. However, we must restrict our domain to 𝒜⊂ℝn, or equivalently to ℬ⊂ℝm. The log-likelihood function
(36)lβ(β;p)=-p∑i=1n|(y-Xβ)i|+nlog(p2(1-e-p))
is a piecewise linear function, up to a constant term, in between the hyperplanes Hi-1, Hi0, and Hi+1, i=1,…,n. Let Xj denote the jth column of X, j=1,2…,m. Let
(37)Sp{X}=Sp{X1,X2,…,Xm}⊂ℝn,
denote the span of the columns of X.Nowlz:cl(𝒜k)↦ℝ has a critical point at z(β)∈𝒜k if and only if the gradient ∇lz=(∂l/∂z1,…,∂l/∂zn)T (evaluated at z(β)) is orthogonal to Sp{X}. The columns of X are tangent vectors to 𝒜 at this point. The gradient ∇lz=-p(sgn(z1(β)),sgn(z2(β)),…,sgn(zn(β)))T is constant in the interior of any orthant. If we travel along a straight line path in any orthant, lz either always increases, always decreases, or remains constant. Hence, we will not find an isolated local maximum or minimum in 𝒜k, an open set, for the constrained lz.We need to be aware of the case where𝒜k lies in or very close to a level set of lz. We might need to test for this. This happens when the sign vector (sgn(z1),…,sgn(zn))T is orthogonal to Sp{X}, or nearly so. In the former case, the constrained l is constant on cl(𝒜k). In the latter case, the constrained l will differ very little around the maximum on cl(𝒜k). If this is the case for all the 𝒜k, then the ML values for the coefficients βj will not be sharply defined (will have large variance).


## 3.3.3. More General Perturbations of the Laplace Probability Density Function

The question is, given a nontrivial perturbing functiong(z;q), does the log-likelihood function attain its maximum at a data point? We have given conditions on g in Lemma 8 which are sufficient to ensure the maximum is attained at a data point. We give a more general criterion in Lemma 10.Assume thatlog(g) is non-linear in any orthant. Otherwise, we can write g in the form of a scaled Laplace distribution and apply Lemma 8. Then l is the sum of an affine function and a non-linear function in any orthant. This affine function is
(38)lp,B=-nlog(Q(p,g1;B))+nlog(p)+∑i=1n(-p|zi|),
the log-likelihood function corresponding to the Laplace distribution. The non-linear function is
(39)lnlin=-nlog(Q(p,g,q;B))+nlog(Q(p,g1;B))+∑i=1nlogg(|zi|;q).
Then l=lp,B+lnlin, and so ∇l=∇lp,B+∇lnlin, where
(40)∇lp,B=-p(sgn(z1),sgn(z2),…,sgn(zn))T,∇lnlin=(dlog(g(z1;q))dz1,dlog(g(z2;q))dz2,…,dlog(g(zn;q))dzn)T.It is possible that there exist orthants𝒪k in which the set 𝒜k(X,y,B) is orthogonal to the gradient ∇lp,B. It is possible that lp,B attains its global maximum (with respect to 𝒜(X,y,B)) on the whole of such an 𝒜k(X,y,B). Hence, it is important to consider the behaviour of lnlin in such orthants.Lemma 10.
Assume bounded support and let(41)γ=sup0≤z≤B{|dg(z;q)/dzg(z;q)|}=sup0≤z≤B{|dlog(g(z;q))dz|},
where sup denotes supremum. Then, if γ<p, the supremum or maximum of l is attained at a data point. In the special case that Sp{X} is orthogonal to (-sgn(z1),-sgn(z2),…,-sgn(zn))T in any orthant, it may be that the supremum is also attained elsewhere.Proof of Lemma10.
Choose𝒜k(X,y,B)=𝒜(X,y,B)∩𝒪k, where k∈{1,2…,2n}, such that the set 𝒜k(X,y,B) is non-empty, choose z∈𝒜k(X,y,B) (an open set relative to 𝒜), and let w(k,X) be the projection of
(42)∇lp,Bp=(-sgn(z1),-sgn(z2),…,-sgn(zn))T
onto Sp{X}. Consider the case w(k,X) is non-zero. Then the affine function lp,B is not constant on 𝒜k(X,y,B). Given any v in 𝒜k(X,y,B), there exists some ϵ>0 and  a map
(43)hp,B:(-ϵ,ϵ)⟶ℝt⟼lp,B(v+tw(k,X))
which is increasing at t=0, since
(44)(dhp,Bdt)t=0=∇lp,BT(d(v+tw(k,X))dt)t=0=∇lp,BTw(k,X)=(w(k,X))Tw(k,X).
Now consider
(45)h:(-ϵ,ϵ)⟶ℝt⟼l(v+tw(k,X)).
Since g(z;q) is continuous and strictly positive on [0,B], γ is finite. We show that if γ<p, dh/dt>0 at t=0, that is, at v. Specifically,
(46)dhdt=-p∑i=1n(sgn(vi+twi(k,X))wi(k,X))+∑i=1n(dg(|zi|;q)/d|zi|)h(t)g(|vi+twi(k,X)|;q)×(sgn(vi+twi(k,X)))wi(k,X)≥-p∑i=1n(sgn(vi+twi(k,X))wi(k,X))+γ∑i=1n(sgn(vi+twi(k,X))wi(k,X))=(1-γp)(w(k,X))Tw(k,X)>0
if γ<p and w(k,X) is non-zero. Then the supremum of l, constrained to cl(𝒜k(X,y,B)), must be attained on the boundary of 𝒜k(X,y,B), relative to 𝒜, and this must be at a data point due to the direction of ∇l.
In the limiting case, wherew(k,X)=0, but we still have γ<p, we find that the supremum of l, constrained to cl(𝒜k(X,y,B)) (relative to 𝒜), must be attained at a data point but may be attained at other points in 𝒜k(X,y,B) as well. Assume this is not the case, that is, there exists no data point at which l attains its supremum when constrained to cl(𝒜k(X,y,B)). We find a contradiction as follows. Assume there exists some v∈𝒜k(X,y,B), an open set (relative to 𝒜), such that
(47)lz(v)-sup{lz(z):z∈bd(𝒪k)∩cl(𝒜k(X,y,B))}=τ>0,
where bd denotes boundary. Imagine perturbing the columns of X slightly (and continuously) to obtain X′ so that lβ, a continuous function, conditional on X, changes by at most τ/3, that is,
(48)|lβ(β;p;X,y)-lβ(β;p;X′,y)|<τ3
for all β∈ℬk(X,y,B)∩ℬk(X′,y,B), and so that now w(k,X′)≠0. We also require the perturbation small enough that
(49)|sup{lz(z;p):z∈bd(𝒪k)∩cl(𝒜k(X,y,B))}(𝒜k(X′,y,B))-sup{lz(z;p):z∈bd(𝒪k)∩cl(𝒜k(X′,y,B))}|<τ3.
Now v=y-Xβv for some unique βv. Let v′=y-X′βv. If the perturbation of X is small enough, then βv∈ℬk(X′,y,B) and so v′∈𝒜k(X′,y,B).
Now,|lz(v′)-lz(v)|<τ/3. Hence,
(50)lz(v′)-sup{lz(z):z∈bd(𝒪k)∩cl(𝒜k(X′,y,B))}>τ3>0,
which is not possible since w(k,X′)≠0. Hence, we have a contradiction. Now, for completeness, assume there exists some r∈bd(𝒜k(X,y,B)), r not a data point (that is some ri=±B) such that
(51)lz(r)-sup{lz(z):z∈bd(𝒪k)∩cl(𝒜k(X,y,B))}=τ>0.
Then there exists v∈𝒜k(X,y,B), such that
(52)lz(r)-lz(v)<τ6.
Hence,
(53)lz(v)-sup{lz(z):z∈bd(𝒪k)∩cl(𝒜k(X,y,B))}>5τ6>0.
The contradiction follows as above. Hence, we have proved the lemma.


## 3.3.4. An Amended Laplace Probability Density Function with Added Kurtosis

We may apply Lemma10 to show that, for the amended Laplace probability density function with added kurtosis (Example 5) and for realistic values of p and q, the maximum of the log-likelihood function must be attained at a data point. Here, g(z;q)=g2(z;q)=1+q(z3-3z)=1+qH3(z) on [0,1], that is, g2(|z|;q)=1+q(H3(|z|)) on [-1,1]. In n dimensions,
(54)l(β;p;X,y)=-nlog(Q(p,g2,q;1))+nlog(p)+∑i=1n(-p|yi-(Xβ)i|)+∑i=1nlog(1+qH3(|yi-(Xβ)i|))=-nlog(Q(p,g2,q;1))+nlog(p)+∑i=1n(-p|yi-(Xβ)i|)+∑i=1nlog(1+q(|yi-(Xβ)i|3-3|yi-(Xβ)i||yi-(Xβ)i|3)q(|yi-(Xβ)i|3).
Consider, for z∈[0,1],
(55)dlog(g2(z;q))dz=d(log(1+q(z3-3z)))dz=3q(z2-1)1+q(z3-3z),d2log(g2(z;q))dz2=d2(log(1+q(z3-3z)))dz2=(3q)(2z-q(3+z4))(1+q(z3-3z))2.
When 0<q<0.5, d2log(g2(z;q))/dz2 is negative when z=0 and positive when z=1, and so the continuous monotonic function log(g2) has a point of inflection, where d2log(g2(z;q))/dz2=0 in (0,1), at say ω(q). Hence, log(g2) and l are both concave and convex on [0,1]. When q is small, the concavity occurs close to the origin and the functions are convex on most of (0,1). Moreover, dlog(g2(z;q))/dz is always negative on (0,1), and its limit as z→∞ is also negative, so we will not get isolated local maxima in one dimension. The function log(g2(z;q)) and its first and second derivatives are plotted in Figures 1, 2, and 3, respectively; setting q=0.025, a typical value. Reading from Figure 2, the upper bound γ is approximately 0.075, when q=0.025. Since typically p≥3, the criterion (γ<p) for Lemma 10 is satisfied.Figure 1
The non-linear part of the log-likelihood function,n=1; log(g(z;q)), q=0.025.Figure 2
The first derivative oflog(g(z;q)), q=0.025.Figure 3
The second derivative oflog(g(z;q)), q=0.025.


## 3.3.5. Non-Increasing Perturbations Both Concave and Convex

In certain situations we might find the criteria for both Lemmas8 and 10 are not satisfied but that logg(z;q) is non-increasing on [0,B] and convex on [w(q),B], where w(q) is small. Since a sum of convex functions is convex (see [8]), l will be convex wherever logg(|zi|;q) is convex for i=1,2,…,n. For example, for g=g2 and B=1, we can prove a partial result as follows. By restricting to the domain
(56)ℬ(X,y,[ω(q),1])=⋂i=1n{β∈ℝm:ω(q)≤|zi(β)|≤1}⊂ℬ(X,y,1)
or equivalently by substituting for each set 𝒜k(X,y,1) the subset
(57){z∈∩𝒜k(X,y,1):ω(q)<|zi|<1,i=1,…,n}⊂𝒜k(X,y,1),
and applying convex function theory as for Lemma 8 we can prove that there exists β1∈ℬ(X,y,[ω(q),1]) at which l attains its maximum, and there exists at least one index i such that zi(β1)=ω(q). Hence, there exists some β2∈ℬ(X,y,1)=ℬ(X,y,[0,1]) at which l attains its maximum, and there exists at least one index i such that zi(β2)≤ω(q). In other words, there exists some point at which l has a maximum at which at least one of the errors is small. It makes sense to search for the maxima of l near or at vertices.Forg=g2, we may apply Lemma 10 and so do not need this partial result, but for a perturbation similar in shape to g2, non-increasing everywhere on [0,1], with logg(z;q) convex everywhere on (ω(g;q),1], concave everywhere on [0,ω(g;q)), for some small positive ω(g;q), and with steep slope at the point of inflection (z=ω(g;q)), such analysis would be useful.


## 4. Statistics for Linear Model Coefficients Assuming Perturbed and Truncated Laplace Response Functions

### 4.1. Dealing with Abrupt Changes in Gradient

The inclusion of the modulus (absolute value) function in the Laplace probability density function (6) (and variations thereof) is the cause of abrupt changes in the gradient of the log-likelihood function lβ (see Section 4.2). This section is devoted to dealing with the problems which are associated with these abrupt changes, encountered when deriving statistical formulae, for example, for standard errors.The fact thatlβ is not differentiable in the classical sense at a local maximum means that the assumptions made in the derivation of the usual classical formulae for the information matrix, the expected value of the Hessian of the log-likelihood function and the variance-covariance matrix for the model coefficients βj, j=1,2,…,m, are not met. For 𝒞2 probability density functions (and 𝒞2 log-likelihood functions), these formulae are derived using Taylor series. We find alternative expressions for these quantities assuming the truncated and/or perturbed Laplace response functions (as defined in Section 2) which are 𝒞3 where the modulus function is non-zero. In Section 5 these expressions will be used to prove the asymptotic convergence of our MLE to a random variable with a normal distribution.


### 4.2. Differentiation in a Generalized Sense

The following generalized functions are required to determine the first and second partial derivatives of the log-likelihood functionlβ, with respect to the coefficients βj. These derivatives are needed for the calculation of the standard errors. We require
(58)sgn(z)={1,z>0,0,z=0,-1,z<0
and δ(z) which is the delta function, that is, δ(z)=0 except at z=0, and ∫-∞∞‍δ(z)dz=1. These expressions and the modulus function are connected by
(59)d|z|dz=sgn(z),(60)dsgn(z)dz=2δ(z),
where the differentiation is taken in the generalized sense (see [9, 10]). Hence, for zi∈[-B,B], the generalized derivative
(61)dg(|zi|;q)dzi=sgn(zi)dg(|zi|;q)d|zi|.
Also, the derivative of the delta function may be defined via integration by parts, assuming h:ℝ→ℝ is 𝒞1, we have
(62)∫-∞∞‍δ′(t)h(t)dt=-∫-∞∞‍δ(t)h′(t)dt=-h′(0).In Section5, we investigate the behaviour of our model as n→∞, and so use subscripts to clarify the variables under consideration (z or β) and/or the dimension of the space(s) under consideration. Using (10) and (59), it follows that
(63)∂l∂zi=∂lz,n∂zi=-psgn(zi)+ddzilog(g(|zi|;q))=-psgn(zi)+sgn(zi)g(|zi|;q)(dg(|zi|;q)d|zi|).
Since z=y-Xn,mβ (see (5)),
(64)∂l∂βj=∂lβ,n,m∂βj=p∑i=1nxijsgn(zi)-∑i=1nxijsgn(zi)g(|zi|;q)dg(|zi|;q)d|zi|.
Letting ∇lβ,n,m denote the gradient of lβ,n,m and letting ∇lz,n denote the gradient of lz,n, we have
(65)∇lβ,n,m=-Xn,mT∇lz,n.
In addition, using (60) and omitting the dependence of g upon its parameters for brevity,
(66)∂2lz,n∂zi2=-2pδ(zi)+2δ(zi)1g(|zi|)(dg(|zi|)d|zi|)+sgn(zi)ddzi(1g(|zi|)dg(|zi|)d|zi|)=-2pδ(zi)+2δ(zi)1g(|zi|)(dg(|zi|)d|zi|)+(sgn(zi))2dd|zi|(1g(|zi|)dg(|zi|)d|zi|)=-2pδ(zi)+2δ(zi)1g(|zi|)(dg(|zi|)d|zi|)+(sgn(zi))2×[1g(|zi|)d2g(|zi|)d|zi|2-(1g(|zi|)2)(dg(|zi|)d|zi|)2]=-2pδ(zi)+2δ(zi)1g(|zi|)(dg(|zi|)d|zi|)+(sgn(zi))2(g(|zi|))2[g(|zi|)d2g(|zi|)d|zi|2-(dg(|zi|)d|zi|)2],
and, if i≠j,
(67)∂2lz,n∂zi∂zj=0.Letℋz,n denote the generalized Hessian of lz,n, where
(68)ℋz,n=(∂2l∂z12∂2l∂z1∂z2⋯∂2l∂z1∂zn∂2l∂z2∂z1∂2l∂z22⋯∂2l∂z2∂zn⋮∂2l∂zn∂z1∂2l∂zn∂z2⋯∂2l∂zn2),
and let E(ℋz,n) denote its expected value. Then ℋz,n and E(ℋz,n) are diagonal matrices. Since the diagonal elements of E(ℋz,n) are all equal (see Section 4.6), this matrix is a multiple of the identity. We have
(69)E(ℋz,n)=E(∂2lz,n∂z12)In,
where In denotes the n×n identity matrix. Let ℋβ,n,m denote the generalized Hessian of lβ,n,m, where
(70)ℋβ,n,m=(∂2l∂β12∂2l∂β1∂β2⋯∂2l∂β1∂βm∂2l∂β2∂β1∂2l∂β22⋯∂2l∂β2∂βm⋮∂2l∂βm∂β1∂2l∂βm∂β2⋯∂2l∂βm2),
and let E(ℋβ,n,m) denote its expected value. Then
(71)ℋβ,n,m=Xn,mTℋz,nXn,m
and
(72)E(ℋβ,n,m)=Xn,mTE(ℋz,n)Xn,m=E(∂2lz,n∂z12)Xn,mTXn,m.
If any zi=0, then the ith diagonal element of ℋz,n is infinite. In this case ℋβ,n,m has infinite components. We prove in Section 4.6 that E(∂2lz,n/∂z12) is finite.


### 4.3. A Classical Relation to Be Generalized

Let𝒥β,n,m denote the m×m Fisher information matrix, where
(73)(𝒥β,n,m)jk=E[(∂lβ,n,m∂βj-E(∂lβ,n,m∂βj))(∂lβ,n,m∂βk-E(∂lβ,n,m∂βk))].
The components of our MLE β^n depend on the errors zi and have a distribution whose variance-covariance matrix is denoted 𝒱β,n,m, where
(74)(𝒱β,n,m)jk=E[(βj-E(βj))(βk-E(βk))],j=1,2,…,m, and k=1,2,…,m. If the log-likelihood function l was sufficiently smooth around the region of interest (i.e., around its maximum value), then Taylor series expansions could be used to derive a relationship between E(ℋβ,n,m), 𝒥β,n,m, and 𝒱β,n,m, namely,
(75)𝒱β,n,m=𝒥β,n,m-1=(-E(ℋβ,n,m))-1,
(see [1]). However, our l is not sufficiently smooth, and so we cannot make use of this relationship without further justification. In general, (75) does not hold, assuming a truncated (and possibly perturbed) Laplace distribution.In Sections4.5 and 4.6 we calculate the expected values of the first and second partial derivatives of the log-likelihood function lβ,n,m using generalized functions; this enables us to derive a generalized Taylor series expansion for the log-likelihood function about a maximum even when the maximum is, for example, on a ridge or at a vertex. Also, this enables us to derive an expression for the generalized variance-covariance matrix for the MLEs of the model coefficients and an expression for the generalized log-likelihood ratio statistic. These formulae differ from the standard formulae for the case of smooth log-likelihood functions, although their form is similar. Specifically, in our case, we prove that 𝒥β,n,m is a multiple of E(ℋβ,n,m), but that the multiple is not −1, rather, it is a negative real number that depends on p, the perturbation g, its parameters q, and the bound B. We prove that our generalized 𝒱β,n,m is a multiple of 𝒥β,n,m-1, where the multiple is a positive real number depending on p, g, q, and B. We assume independent error distributions.


### 4.4. The Mean and Variance of the Partial Derivatives of the Log-Likelihood Function

The mean and the variance of∂l/∂zi are required in the calculation of 𝒥β,n,m. Recall f(z1,…,zn;p)=Πi=1nf(zi;p) is the joint probability density function for the independent deviations (errors), and that
(76)l=lz,n=∑i=1nlog(f(zi;p)),
so
(77)∂l∂zi=∂lz,n∂zi=∂log(f(zi;p))∂zi.
Let
(78)μi(p,g,q;B)=E(∂lz,n∂zi)=∫ΩB‍(∂log(f(zi;p))∂zi)(f(z;p))dz1⋯dzn=∫-BB‍(∂log(f(zi;p))∂zi)(f(zi;p))dzi,i=1,…,n; then μ(p,g,q;B)=μi(p,g,q;B) is independent of index i. Since f(zi;p)=f(-zi;p), i=1,…,n, L and l are symmetric about the origin, and so
(79)μ=μ(p,g,q;B)=0
for any choice of p, g, q, and B. Let
(80)νi(p,g,q;B)=var(∂lz,n∂zi)=∫ΩB‍(∂log(f(zi;p))∂zi)2(f(z;p))dz1⋯dzn=∫-BB‍(∂log(f(zi;p))∂zi)2(f(zi;p))dzi,
where var denotes variance. Using (63) for ∂lz,n/∂zi,
(81)νi(p,g,q;B)=2∫0B‍(F(zi;p))2f(zi;p)dzi,
where, for z∈(0,B],
(82)F(z;p)=F(z;p,g,q)=-p+1g(z)dg(z;q)dz
and F(0;p)=0. Since νi(p,g,q;B) is independent of index i, we omit the subscript. Hence, the information matrix 𝒥z,n, defined with respect to lz,n, is
(83)𝒥z,n=νIn.Ifg(z;q)=g1(z)=1, then F(z;p)=F(z;p,g1)=-p, and so
(84)ν(p,g1;B)=p2.
For the nontrivial perturbing function g2, for fixed p and q, one can show that ν depends on B by direct calculation.


### 4.5. The Information Matrix

We calculate the information matrix𝒥β,n,m (conditional on p, g, q, and B). We are trying to quantify the steepness of the slope of lβ,n,m around a maximum, in the directions represented by the coefficients βj. If lβ,n,m is very flat in one direction, then the model coefficient representing that direction is not well defined (will have large variance). When calculating 𝒥β,n,m, we are taking into account the behaviour of the gradient of lβ,n,m on a whole neighbourhood of the MLE (how it differs from the expected value) and discontinuities on sets of measure zero can be accommodated. Recall (64), for j=1,2,…,m,
(85)∂lβ,n,m∂βj=p∑i=1nxijsgn(zi)-∑i=1nxijsgn(zi)g(|zi|;q)(dg(|zi|;q)d|zi|).
Hence, omitting some subscripts on l for brevity, we obtain
(86)E(∂lβ,n,m∂βj)=∫ΩB‍∂l∂βj(z1,…,zn)f(z1,…,zn;p)dz1⋯dzn=∫ΩB‍(∑i=1n(∂l∂zi∂zi∂βj))×f(z1;p)⋯f(zn;p)dz1⋯dzn=∑i=1n(-xij)∫-BB‍(∂l∂zi)f(zi;p)dzi=∑i=1n(-xij)E(∂l∂zi)=∑i=1n(-xij)μ=0.
Since the expectations of the partial derivatives of lβ,n,m are zero, the diagonal elements of 𝒥β,n,m are an indication of the steepness of the gradient around the maximum likelihood estimate. Now,
(87)(𝒥β,n,m)jk=E[(∂l∂βj-E(∂l∂βj))(∂l∂βk-E(∂l∂βk))]=E[(∂l∂βj)(∂l∂βk)]=∫ΩB‍(∂l∂βj)(∂l∂βk)f(z;p)dz1⋯dzn=∫ΩB‍(∑i=1n(∂l∂zi∂zi∂βj))×(∑t=1n(∂l∂zt∂zt∂βk))f(z;p)dz1⋯dzn=∫ΩB‍(∑i=1n(-xij∂l∂zi))×(∑t=1n(-xtk∂l∂zt))×f(z1;p)⋯f(zn;p)dz1⋯dzn=∫ΩB‍(∑i=1n(xij∂l∂zi)(xik∂l∂zi))×f(z1;p)⋯f(zn;p)dz1⋯dzn,
since the cross-terms indexed by i≠t equal zero by symmetry, that is,
(88)∫zi=-B,zt=-Bzi=B,zt=B‍(xij∂l∂zi)(xtk∂l∂zt)f(zi;p)×f(zt;p)dzidzt=xijxtkμ2=0.
So
(89)(𝒥β,n,m)jk=∫ΩB‍(∑i=1n(xijxik(∂l∂zi)2))×f(z1;p)⋯f(zn;p)dz1⋯dzn=∑i=1n(xijxik)∫-BB‍(∂l∂zi)2f(zi;p)dzi=ν(p,g,q;B)∑i=1n(xijxik).
Hence,
(90)𝒥β,n,m=Xn,mT𝒥z,nXn,m=ν(p,g,q;B)Xn,mTXn,m.Ifg(z;q)=g1(z)=1, then 𝒥β,n,m=p2Xn,mTXn,m.


### 4.6. The Expected Value of the Generalized Hessian

In order to calculate the expected value of the generalized Hessian we require(91)E(ℋβ,n,m)=Xn,mTE(ℋz,n)Xn,m=E(∂2lz,n∂z12)Xn,mTXn,m,
and
(92)∂2lz,n∂zi2=-2pδ(zi)+2δ(zi)1g(|zi|;q)(dg(|zi|;q)d|zi|)+(sgn(zi))2(g(|zi|;q))2×[g(|zi|;q)d2g(|zi|;q)d|zi|2-(dg(|zi|;q)d|zi|)2].
Let
(93)ζi(p,g,q;B)=E(∂2lz,n∂zi2)=∫ΩB‍∂2l∂zi2f(z;p)dz1⋯zn=∫-BB‍∂2l∂zi2f(zi;p)dzi=-2pf(0;p)+2f(0;p)1g(0;q)×(dg(|zi|;q)d|zi|)zi=0+∫-BB‍1g(|zi|;q)d2g(|zi|;q)d|zi|2f(zi;p)dzi-∫-BB‍1(g(|zi|;q))2(dg(|zi|;q)d|zi|)2×f(zi;p)dzi,
since (dg(|zi|;q)/d|zi|)zi=0=(dg(u;q)/du)u=0, where u=|zi|. The quantity ζi(p,g,q;B) does not depend on i, so we omit the index. Hence,
(94)E(ℋβ,n,m)=ζ(p,g,q;B)Xn,mTXn,m.
If g(z;q)=g1(z)=1, then
(95)ζ(p,g1;B)=-2p2Q(p,g1;B)=-p2(1-e-pB),
and so although the generalised Hessian ℋβ,n,m has infinite elements, its expected value is negative definite. By continuity, if g(z;q) is close enough to the constant map g1, the expected value E(ℋβ,n,m) will still be negative definite. Note that if g has negative slope at the origin, the peak of lz,n at the origin becomes sharper, compared to that for the case g=g1. If g has positive slope at the origin, we see the opposite effect.


### 4.7. The Generalized Variance-Covariance Matrix for the Model Coefficients

We use a generalized Taylor series expansion (in the coefficientsβj) to approximate lβ,n,m by a negative definite quadratic function about a local maximum. Although we know that, for finite n, this approximation is not exact, we show in Section 5 that we would expect it to become more accurate as n→∞. Assuming Xn,m and y∈ℝn are given, conditional on p and q, we could write a Taylor series expansion for lβ,n,m about a ML estimator β^=β^n∈ℬ as follows:
(96)lβ,n,m(β)≈lβ,n,m(β^n)+(β-β^n)T(∇lβ)β^n+(12)(β-β^n)Tℋβ^n(β-β^n)+⋯,
if f and hence Lβ,n,m and lβ,n,m were sufficiently smooth. Now our probability density function f and hence lβ,n,m are not sufficiently smooth, but we can replace the second derivative of lβ,n,m by its expected value using generalized functions. Let Δβ=β-β^n. This yields the approximation
(97)lβ,n,m(β)≈lβ,n,m(β^n)+(12)(Δβ)TE(ℋβ,n,m)Δβ.
Assuming E(ℋβ,n,m) is nonsingular, which is true when g=g1, we ignore higher order terms. Equation (97) provides an indication of the behaviour of lβ,n,m about a maximum since, for example, if g=g1, then E(ℋβ,n,m) is negative definite.Next we consider the score function∇lβ,n,m and use a Taylor series approximation incorporating generalized functions (about a local maximum β^n) to derive a relationship between the expected value of the generalized Hessian E(ℋβ,n,m), the information matrix 𝒥β,n,m, and the generalized variance-covariance matrix 𝒱β,n,m. This approximation is
(98)(∇lβ,n,m)β≈(∇lβ,n,,m)β^n+E(ℋβ,n,m)Δβ=E(ℋβ,n,m)Δβ.
Since E(∇lβ,n,m)=0 and E(ℋβ,n,m) has full rank, E(β^n)=β. We multiply each side by its own transpose and take expected values to obtain
(99)E[(∇lβ,n,m)β(∇lβ,n,m)βT]=(E(ℋβ,n,m))E[Δβ(Δβ)T](E(ℋβ,n,m))T.
Hence,
(100)𝒥β,n,m=(E(ℋβ,n,m))𝒱β,n,m(E(ℋβ,n,m))T=(E(ℋβ,n,m))𝒱β,n,m(E(ℋβ,n,m)),
so
(101)𝒱β,n,m=𝒱β,n,m(p,g,q,Xn,m;B)=(E(ℋβ,n,m))-1𝒥β,n,m(E(ℋβ,n,m))-1=ν(p,g,q;B)(ζ(p,g,q;B))2(Xn,mTXn,m)-1.
Here, we require Xn,m to have full rank and that ζ(p,g,q;B)≠0, which is certainly true when g=g1.


### 4.8. Generalized Statistical Expressions and Relations

We have shown that the expected value of the generalized HessianE(ℋβ,n,m) and the information matrix 𝒥β,n,m are multiples of Xn,mTXn,m. Specifically,
(102)E(ℋβ,n,m)=Xn,mTE(ℋz,n)Xn,m=E(∂2l(z1;p,g,q)∂z12)Xn,mTXn,m=ζ(p,g,q;B)Xn,mTXn,m,𝒥β,n,m=E((∂l(z1;p,g,q)∂z1)2)Xn,mTXn,m=ν(p,g,q;B)Xn,mTXn,m
(see (94) and (90)). Hence, if ν(p,g,q;B)≠0 (true if g=g1), then
(103)E(ℋβ,n,m)=E(∂2l(z1;p,g,q)/∂z12)E((∂l(z1;p,g,q)/∂z1)2)𝒥β,n,m=ζ(p,g,q;B)ν(p,g,q;B)𝒥β,n,m.
In addition, assuming that the scalar ζ(p,g,q;B) is also non-zero (true if g=g1) and that Xn,m has full rank m, we have proved the following relations:
(104)(𝒱β,n,m(p,g,q,Xn,m;B))-1=(ζ(p,g,q;B))2ν(p,g,q;B)Xn,mTXn,m(105)=ζ(p,g,q;B)ν(p,g,q;B)E(ℋβ,n,m)(106)=(ζ(p,g,q;B))2(ν(p,g,q;B))2𝒥β,n,m,
to be used in the derivation of the generalized log-likelihood ratio statistic (Section 4.11).


### 4.9. Statistical Relations for the Laplace Distribution

Ifg(z,q)=g1(z)=1, then using (21), (84), and (95):
(107)Q(p,g1;B)=2(1-e-pB),ν(p,g1;B)=p2,ζ(p,g1;B)=-2p2Q(p,g1;B)=-p2(1-e-pB).
Hence,
(108)𝒥β,n,m=p2Xn,mTXn,m,E(ℋβ,n,m)=(-p2)(1-e-pB)Xn,mTXn,m,
so
(109)𝒱β,n,m=(1-e-pB)2p2(Xn,mTXn,m)-1.
Here,
(110)𝒥β,n,m=-(1-e-pB)E(ℋβ,n,m),(111)𝒱β,n,m=(1-e-pB)2𝒥β,n,m-1,
where p>0, B>0, and Xn,m has full rank m. Hence, (75), derived for 𝒞2 distributions, does not hold for the Laplace distribution with bounded support. However, (75) describes the limiting behaviour, as B→∞.


### 4.10. Statistical Relations for a Laplace Distribution with Added Kurtosis

Now consider our motivating example, a Laplace distribution with bounded support[-1,1] amended by adding kurtosis. In this case g(z,q)=g2(z;q)=1+qH3(z) (see (14)), B=1,
(112)Q(p,q)=2[(p3-3qp2+6q)-e-p(p3(1-2q)+6pq+6q)]p3,f(z;p,q)=p4e(-p|z|)[1+qH3(|z|)]2[(p3-3qp2+6q)-e-p(p3(1-2q)+6pq+6q)].
Also, using (82) and (81), we obtain
(113)F(z;p,q)=F(z;p,g2,q)=-p+3q(z2-1)1+q(z3-3z),ν(p,g2,q;1)=2∫01‍(F(z;p,q))2f(z;p,q)dz,
and using (93), we obtain
(114)ζ(p,g2,q;1)=-2pf(0;p,q)+2f(0;p,q)×1g(0;q)(dg(|z|;q)d|z|)z=0+∫-11‍1g(|z|;q)d2g(|z|;q)d|z|2f(z;p,q)dz-∫-11‍1(g(|z|;q))2(dg(|z|;q)d|z|)2×f(z;p,q)dz=-2p2Q(p,q)+-6pqQ(p,q)+2∫01‍1g(z;q)d2g(z;q)dz2f(z;p,q)dz-2∫01‍1(g(z;q))2(dg(z;q)dz)2×f(z;p,q)dz=-2p2-6pqQ(p,q)+2∫01‍6qz(1+q(z3-3z))f(z;p,q)dz-2∫01‍1(1+q(z3-3z)2)×(3q(z2-1))2f(z;p,q)dz.
For typical values p=5.254 and q=0.025, numerical integration gives ν=28.3561 and ζ=-28.4957, to four decimal places. In this case,
(115)𝒥β,n,m=28.3561Xn,mTXn,m,E(ℋβ,n,m)=-28.4957Xn,mTXn,m,
and so
(116)𝒱β,n,m=28.3561(-28.4957)2(Xn,mTXn,m)-1.
Here
(117)𝒥β,n,m=28.3561-28.4957E(ℋβ,n,m),𝒱β,n,m=(28.3561)2(-28.4957)2𝒥β,n,m-1.
In this example, (75) could be used as an approximation but does not exactly describe the relationship between E(ℋβ,n,m), 𝒥β,n,m, and 𝒱β,n,m.


### 4.11. The Generalized Log-Likelihood Ratio Statistic

The log-likelihood ratio statistic enables us to assess the adequacy of a model. It enables us to compare a model withM coefficients (parameters) with a model of interest which differs only in that it has fewer coefficients, say P, with M>P≥1, see [1], for example. We wish to compare a linear model with m coefficients, the βj, for j=1,2,…,m with a lesser linear model with fewer coefficients. The aim is to decide whether or not the excluded coefficients are useful. Our comparison is conditional on the parameters p and q.LetLρ,n,M(ρ1,…,ρM;p,q;Xn,M,y) denote the likelihood function for the model with M parameters, and let Lψ,n,P(ψ1,…,ψP;p,q;Xn,P,y) denote the likelihood function for the lesser model. Let ρ^ (or ρ^n) be the maximum likelihood estimator of ρ=(ρ1,…,ρM)T, and let ψ^ (or ψ^n) be the maximum likelihood estimator of ψ=(ψ1,…,ψP)T. Then the likelihood ratio
(118)λ=Lρ,n,M(ρ^;p,q;Xn,M,y)Lψ,n,P(ψ^;p,q;Xn,P,y)
is a ratio of two probabilities and will be greater than one since the model with M parameters provides the more complete description of the model. In our application, M=m, ρ=β, and usually (but not necessarily) P=m-1. We show that a multiple (conditional on p and q) of
(119)log(λ)=log(Lρ,n,M(ρ^;p,q;Xn,M,y))-log(Lψ,n,P(ψ^;p,q;Xn,P,y))
has a chi-squared distribution as follows.The derivation of the log-likelihood ratio statistic (for smooth functions) may be found in the textbooks. For example, for generalized linear models (see [1]) where the log-likelihood function is smooth in a neighbourhood of its maximum, D=2log(λ) is distributed approximately as χ2(M-P,δD). Here δD is a noncentrality parameter, a positive constant which will be near zero if the lesser model fits the data almost as well as the model with more coefficients. Consider our situation in which the log-likelihood function is continuous at a maximum, but where this maximum occurs at a vertex or possibly on a ridge in (β,l)-space, and generalized calculus is required to consider Taylor series expansions. Then, as in (97), around ρ^, replacing the Hessian of log(LM)=lM by its expected value, we obtain the following approximation:
(120)log(Lρ,n,M(E(ρ^);p,q;Xn,M,y))-log(Lρ,n,M(ρ^;p,q;Xn,M,y))=log(Lρ,n,M(E(ρ^)))-log(Lρ,n,M(ρ^))≈12(E(ρ^)-ρ^)TE(ℋρ,n,M)(E(ρ^)-ρ^).So, by (105) assuming that the scalar factors ν(p,g,q;B) and ζ(p,g,q;B) are both non-zero and that Xn,M has full rank
(121)2E((∂2l(z1;p,g,q))/∂z12)E(((∂l(z1;p,g,q))/∂z1)2)(lρ,n,M(E(ρ^))-lρ,n,M(ρ^))≈(E(ρ^)-ρ^)Tζ(p,g,q;B)ν(p,g,q;B)E(ℋρ,n,M)(E(ρ^)-ρ^)=(E(ρ^)-ρ^)T𝒱ρ,n,M-1(E(ρ^)-ρ^),
which has the distribution χ2(M) if the MLE has a normal distribution. We show in Section 5 that the distribution of the MLE is asymptotically normal. Hence, if n is large enough, this will be a good approximation. Similarly
(122)2E((∂2l(z1;p,g,q))/∂z12)E(((∂l(z1;p,g,q))/∂z1)2)(lψ,n,P(E(ψ^))-lψ,n,P(ψ^))≈(E(ψ^)-ψ^)Tζ(p,g,q;B)ν(p,g,q;B)E(ℋψ,n,P)(E(ψ^)-ψ^)=(E(ψ^)-ψ^)T𝒱ψ,n,P-1(E(ψ^)-ψ^),
which has the distribution χ2(P), approximately. Noting that
(123)ζ(p,g,q;B)ν(p,g,q;B)=-1(1-e-pB)<0,
when g=g1, let
(124)Dgen=-2ζ(p,g,q;B)ν(p,g,q;B)log(λ).
Then Dgen, the log-likelihood ratio statistic calculated with generalized functions, a positive number, may be expressed as
(125)Dgen=-2ζ(p,g,q;B)ν(p,g,q;B)×(lρ,n,M(ρ^;p,q;Xn,M,y)-lψ,n,P(ψ^;p,q;Xn,P,y))=-2ζ(p,g,q;B)ν(p,g,q;B)×(lρ,n,M(ρ^;p,q;Xn,M,y)-lρ,n,M(E(ρ^);p,q;Xn,M,y))+2ζ(p,g,q;B)ν(p,g,q;B)×(lψ,n,P(ψ^;p,q;Xn,P,y)-lψ,n,P(E(ψ^);p,q;Xn,P,y))-2ζ(p,g,q;B)ν(p,g,q;B)×(lρ,n,M(E(ρ^);p,q;Xn,M,y)-lψ,n,P(E(ψ^);p,q;Xn,P,y)).
Hence, Dgen is the sum of three terms; the first (positive) has the distribution χ2(M) (approximately). The second (negative) has the distribution χ2(P) (approximately). The third is a positive constant (say δgen) that depends on p and q. Hence, Dgen is distributed approximately as χ2(M-P,δgen). If the lesser model gives a good description of the data, then δgen will be small. The generalized log-likelihood ratio statistic Dgen is easily calculated and is hence a potentially useful statistic for assessing our linear model for which the log-likelihood function is not 𝒞2. Note if g=g1, Dgen=(2log(λ))/(1-e-pB), and so Dgen→2log(λ), as B→∞.


## 4.1. Dealing with Abrupt Changes in Gradient

The inclusion of the modulus (absolute value) function in the Laplace probability density function (6) (and variations thereof) is the cause of abrupt changes in the gradient of the log-likelihood function lβ (see Section 4.2). This section is devoted to dealing with the problems which are associated with these abrupt changes, encountered when deriving statistical formulae, for example, for standard errors.The fact thatlβ is not differentiable in the classical sense at a local maximum means that the assumptions made in the derivation of the usual classical formulae for the information matrix, the expected value of the Hessian of the log-likelihood function and the variance-covariance matrix for the model coefficients βj, j=1,2,…,m, are not met. For 𝒞2 probability density functions (and 𝒞2 log-likelihood functions), these formulae are derived using Taylor series. We find alternative expressions for these quantities assuming the truncated and/or perturbed Laplace response functions (as defined in Section 2) which are 𝒞3 where the modulus function is non-zero. In Section 5 these expressions will be used to prove the asymptotic convergence of our MLE to a random variable with a normal distribution.


## 4.2. Differentiation in a Generalized Sense

The following generalized functions are required to determine the first and second partial derivatives of the log-likelihood functionlβ, with respect to the coefficients βj. These derivatives are needed for the calculation of the standard errors. We require
(58)sgn(z)={1,z>0,0,z=0,-1,z<0
and δ(z) which is the delta function, that is, δ(z)=0 except at z=0, and ∫-∞∞‍δ(z)dz=1. These expressions and the modulus function are connected by
(59)d|z|dz=sgn(z),(60)dsgn(z)dz=2δ(z),
where the differentiation is taken in the generalized sense (see [9, 10]). Hence, for zi∈[-B,B], the generalized derivative
(61)dg(|zi|;q)dzi=sgn(zi)dg(|zi|;q)d|zi|.
Also, the derivative of the delta function may be defined via integration by parts, assuming h:ℝ→ℝ is 𝒞1, we have
(62)∫-∞∞‍δ′(t)h(t)dt=-∫-∞∞‍δ(t)h′(t)dt=-h′(0).In Section5, we investigate the behaviour of our model as n→∞, and so use subscripts to clarify the variables under consideration (z or β) and/or the dimension of the space(s) under consideration. Using (10) and (59), it follows that
(63)∂l∂zi=∂lz,n∂zi=-psgn(zi)+ddzilog(g(|zi|;q))=-psgn(zi)+sgn(zi)g(|zi|;q)(dg(|zi|;q)d|zi|).
Since z=y-Xn,mβ (see (5)),
(64)∂l∂βj=∂lβ,n,m∂βj=p∑i=1nxijsgn(zi)-∑i=1nxijsgn(zi)g(|zi|;q)dg(|zi|;q)d|zi|.
Letting ∇lβ,n,m denote the gradient of lβ,n,m and letting ∇lz,n denote the gradient of lz,n, we have
(65)∇lβ,n,m=-Xn,mT∇lz,n.
In addition, using (60) and omitting the dependence of g upon its parameters for brevity,
(66)∂2lz,n∂zi2=-2pδ(zi)+2δ(zi)1g(|zi|)(dg(|zi|)d|zi|)+sgn(zi)ddzi(1g(|zi|)dg(|zi|)d|zi|)=-2pδ(zi)+2δ(zi)1g(|zi|)(dg(|zi|)d|zi|)+(sgn(zi))2dd|zi|(1g(|zi|)dg(|zi|)d|zi|)=-2pδ(zi)+2δ(zi)1g(|zi|)(dg(|zi|)d|zi|)+(sgn(zi))2×[1g(|zi|)d2g(|zi|)d|zi|2-(1g(|zi|)2)(dg(|zi|)d|zi|)2]=-2pδ(zi)+2δ(zi)1g(|zi|)(dg(|zi|)d|zi|)+(sgn(zi))2(g(|zi|))2[g(|zi|)d2g(|zi|)d|zi|2-(dg(|zi|)d|zi|)2],
and, if i≠j,
(67)∂2lz,n∂zi∂zj=0.Letℋz,n denote the generalized Hessian of lz,n, where
(68)ℋz,n=(∂2l∂z12∂2l∂z1∂z2⋯∂2l∂z1∂zn∂2l∂z2∂z1∂2l∂z22⋯∂2l∂z2∂zn⋮∂2l∂zn∂z1∂2l∂zn∂z2⋯∂2l∂zn2),
and let E(ℋz,n) denote its expected value. Then ℋz,n and E(ℋz,n) are diagonal matrices. Since the diagonal elements of E(ℋz,n) are all equal (see Section 4.6), this matrix is a multiple of the identity. We have
(69)E(ℋz,n)=E(∂2lz,n∂z12)In,
where In denotes the n×n identity matrix. Let ℋβ,n,m denote the generalized Hessian of lβ,n,m, where
(70)ℋβ,n,m=(∂2l∂β12∂2l∂β1∂β2⋯∂2l∂β1∂βm∂2l∂β2∂β1∂2l∂β22⋯∂2l∂β2∂βm⋮∂2l∂βm∂β1∂2l∂βm∂β2⋯∂2l∂βm2),
and let E(ℋβ,n,m) denote its expected value. Then
(71)ℋβ,n,m=Xn,mTℋz,nXn,m
and
(72)E(ℋβ,n,m)=Xn,mTE(ℋz,n)Xn,m=E(∂2lz,n∂z12)Xn,mTXn,m.
If any zi=0, then the ith diagonal element of ℋz,n is infinite. In this case ℋβ,n,m has infinite components. We prove in Section 4.6 that E(∂2lz,n/∂z12) is finite.


## 4.3. A Classical Relation to Be Generalized

Let𝒥β,n,m denote the m×m Fisher information matrix, where
(73)(𝒥β,n,m)jk=E[(∂lβ,n,m∂βj-E(∂lβ,n,m∂βj))(∂lβ,n,m∂βk-E(∂lβ,n,m∂βk))].
The components of our MLE β^n depend on the errors zi and have a distribution whose variance-covariance matrix is denoted 𝒱β,n,m, where
(74)(𝒱β,n,m)jk=E[(βj-E(βj))(βk-E(βk))],j=1,2,…,m, and k=1,2,…,m. If the log-likelihood function l was sufficiently smooth around the region of interest (i.e., around its maximum value), then Taylor series expansions could be used to derive a relationship between E(ℋβ,n,m), 𝒥β,n,m, and 𝒱β,n,m, namely,
(75)𝒱β,n,m=𝒥β,n,m-1=(-E(ℋβ,n,m))-1,
(see [1]). However, our l is not sufficiently smooth, and so we cannot make use of this relationship without further justification. In general, (75) does not hold, assuming a truncated (and possibly perturbed) Laplace distribution.In Sections4.5 and 4.6 we calculate the expected values of the first and second partial derivatives of the log-likelihood function lβ,n,m using generalized functions; this enables us to derive a generalized Taylor series expansion for the log-likelihood function about a maximum even when the maximum is, for example, on a ridge or at a vertex. Also, this enables us to derive an expression for the generalized variance-covariance matrix for the MLEs of the model coefficients and an expression for the generalized log-likelihood ratio statistic. These formulae differ from the standard formulae for the case of smooth log-likelihood functions, although their form is similar. Specifically, in our case, we prove that 𝒥β,n,m is a multiple of E(ℋβ,n,m), but that the multiple is not −1, rather, it is a negative real number that depends on p, the perturbation g, its parameters q, and the bound B. We prove that our generalized 𝒱β,n,m is a multiple of 𝒥β,n,m-1, where the multiple is a positive real number depending on p, g, q, and B. We assume independent error distributions.


## 4.4. The Mean and Variance of the Partial Derivatives of the Log-Likelihood Function

The mean and the variance of∂l/∂zi are required in the calculation of 𝒥β,n,m. Recall f(z1,…,zn;p)=Πi=1nf(zi;p) is the joint probability density function for the independent deviations (errors), and that
(76)l=lz,n=∑i=1nlog(f(zi;p)),
so
(77)∂l∂zi=∂lz,n∂zi=∂log(f(zi;p))∂zi.
Let
(78)μi(p,g,q;B)=E(∂lz,n∂zi)=∫ΩB‍(∂log(f(zi;p))∂zi)(f(z;p))dz1⋯dzn=∫-BB‍(∂log(f(zi;p))∂zi)(f(zi;p))dzi,i=1,…,n; then μ(p,g,q;B)=μi(p,g,q;B) is independent of index i. Since f(zi;p)=f(-zi;p), i=1,…,n, L and l are symmetric about the origin, and so
(79)μ=μ(p,g,q;B)=0
for any choice of p, g, q, and B. Let
(80)νi(p,g,q;B)=var(∂lz,n∂zi)=∫ΩB‍(∂log(f(zi;p))∂zi)2(f(z;p))dz1⋯dzn=∫-BB‍(∂log(f(zi;p))∂zi)2(f(zi;p))dzi,
where var denotes variance. Using (63) for ∂lz,n/∂zi,
(81)νi(p,g,q;B)=2∫0B‍(F(zi;p))2f(zi;p)dzi,
where, for z∈(0,B],
(82)F(z;p)=F(z;p,g,q)=-p+1g(z)dg(z;q)dz
and F(0;p)=0. Since νi(p,g,q;B) is independent of index i, we omit the subscript. Hence, the information matrix 𝒥z,n, defined with respect to lz,n, is
(83)𝒥z,n=νIn.Ifg(z;q)=g1(z)=1, then F(z;p)=F(z;p,g1)=-p, and so
(84)ν(p,g1;B)=p2.
For the nontrivial perturbing function g2, for fixed p and q, one can show that ν depends on B by direct calculation.


## 4.5. The Information Matrix

We calculate the information matrix𝒥β,n,m (conditional on p, g, q, and B). We are trying to quantify the steepness of the slope of lβ,n,m around a maximum, in the directions represented by the coefficients βj. If lβ,n,m is very flat in one direction, then the model coefficient representing that direction is not well defined (will have large variance). When calculating 𝒥β,n,m, we are taking into account the behaviour of the gradient of lβ,n,m on a whole neighbourhood of the MLE (how it differs from the expected value) and discontinuities on sets of measure zero can be accommodated. Recall (64), for j=1,2,…,m,
(85)∂lβ,n,m∂βj=p∑i=1nxijsgn(zi)-∑i=1nxijsgn(zi)g(|zi|;q)(dg(|zi|;q)d|zi|).
Hence, omitting some subscripts on l for brevity, we obtain
(86)E(∂lβ,n,m∂βj)=∫ΩB‍∂l∂βj(z1,…,zn)f(z1,…,zn;p)dz1⋯dzn=∫ΩB‍(∑i=1n(∂l∂zi∂zi∂βj))×f(z1;p)⋯f(zn;p)dz1⋯dzn=∑i=1n(-xij)∫-BB‍(∂l∂zi)f(zi;p)dzi=∑i=1n(-xij)E(∂l∂zi)=∑i=1n(-xij)μ=0.
Since the expectations of the partial derivatives of lβ,n,m are zero, the diagonal elements of 𝒥β,n,m are an indication of the steepness of the gradient around the maximum likelihood estimate. Now,
(87)(𝒥β,n,m)jk=E[(∂l∂βj-E(∂l∂βj))(∂l∂βk-E(∂l∂βk))]=E[(∂l∂βj)(∂l∂βk)]=∫ΩB‍(∂l∂βj)(∂l∂βk)f(z;p)dz1⋯dzn=∫ΩB‍(∑i=1n(∂l∂zi∂zi∂βj))×(∑t=1n(∂l∂zt∂zt∂βk))f(z;p)dz1⋯dzn=∫ΩB‍(∑i=1n(-xij∂l∂zi))×(∑t=1n(-xtk∂l∂zt))×f(z1;p)⋯f(zn;p)dz1⋯dzn=∫ΩB‍(∑i=1n(xij∂l∂zi)(xik∂l∂zi))×f(z1;p)⋯f(zn;p)dz1⋯dzn,
since the cross-terms indexed by i≠t equal zero by symmetry, that is,
(88)∫zi=-B,zt=-Bzi=B,zt=B‍(xij∂l∂zi)(xtk∂l∂zt)f(zi;p)×f(zt;p)dzidzt=xijxtkμ2=0.
So
(89)(𝒥β,n,m)jk=∫ΩB‍(∑i=1n(xijxik(∂l∂zi)2))×f(z1;p)⋯f(zn;p)dz1⋯dzn=∑i=1n(xijxik)∫-BB‍(∂l∂zi)2f(zi;p)dzi=ν(p,g,q;B)∑i=1n(xijxik).
Hence,
(90)𝒥β,n,m=Xn,mT𝒥z,nXn,m=ν(p,g,q;B)Xn,mTXn,m.Ifg(z;q)=g1(z)=1, then 𝒥β,n,m=p2Xn,mTXn,m.


## 4.6. The Expected Value of the Generalized Hessian

In order to calculate the expected value of the generalized Hessian we require(91)E(ℋβ,n,m)=Xn,mTE(ℋz,n)Xn,m=E(∂2lz,n∂z12)Xn,mTXn,m,
and
(92)∂2lz,n∂zi2=-2pδ(zi)+2δ(zi)1g(|zi|;q)(dg(|zi|;q)d|zi|)+(sgn(zi))2(g(|zi|;q))2×[g(|zi|;q)d2g(|zi|;q)d|zi|2-(dg(|zi|;q)d|zi|)2].
Let
(93)ζi(p,g,q;B)=E(∂2lz,n∂zi2)=∫ΩB‍∂2l∂zi2f(z;p)dz1⋯zn=∫-BB‍∂2l∂zi2f(zi;p)dzi=-2pf(0;p)+2f(0;p)1g(0;q)×(dg(|zi|;q)d|zi|)zi=0+∫-BB‍1g(|zi|;q)d2g(|zi|;q)d|zi|2f(zi;p)dzi-∫-BB‍1(g(|zi|;q))2(dg(|zi|;q)d|zi|)2×f(zi;p)dzi,
since (dg(|zi|;q)/d|zi|)zi=0=(dg(u;q)/du)u=0, where u=|zi|. The quantity ζi(p,g,q;B) does not depend on i, so we omit the index. Hence,
(94)E(ℋβ,n,m)=ζ(p,g,q;B)Xn,mTXn,m.
If g(z;q)=g1(z)=1, then
(95)ζ(p,g1;B)=-2p2Q(p,g1;B)=-p2(1-e-pB),
and so although the generalised Hessian ℋβ,n,m has infinite elements, its expected value is negative definite. By continuity, if g(z;q) is close enough to the constant map g1, the expected value E(ℋβ,n,m) will still be negative definite. Note that if g has negative slope at the origin, the peak of lz,n at the origin becomes sharper, compared to that for the case g=g1. If g has positive slope at the origin, we see the opposite effect.


## 4.7. The Generalized Variance-Covariance Matrix for the Model Coefficients

We use a generalized Taylor series expansion (in the coefficientsβj) to approximate lβ,n,m by a negative definite quadratic function about a local maximum. Although we know that, for finite n, this approximation is not exact, we show in Section 5 that we would expect it to become more accurate as n→∞. Assuming Xn,m and y∈ℝn are given, conditional on p and q, we could write a Taylor series expansion for lβ,n,m about a ML estimator β^=β^n∈ℬ as follows:
(96)lβ,n,m(β)≈lβ,n,m(β^n)+(β-β^n)T(∇lβ)β^n+(12)(β-β^n)Tℋβ^n(β-β^n)+⋯,
if f and hence Lβ,n,m and lβ,n,m were sufficiently smooth. Now our probability density function f and hence lβ,n,m are not sufficiently smooth, but we can replace the second derivative of lβ,n,m by its expected value using generalized functions. Let Δβ=β-β^n. This yields the approximation
(97)lβ,n,m(β)≈lβ,n,m(β^n)+(12)(Δβ)TE(ℋβ,n,m)Δβ.
Assuming E(ℋβ,n,m) is nonsingular, which is true when g=g1, we ignore higher order terms. Equation (97) provides an indication of the behaviour of lβ,n,m about a maximum since, for example, if g=g1, then E(ℋβ,n,m) is negative definite.Next we consider the score function∇lβ,n,m and use a Taylor series approximation incorporating generalized functions (about a local maximum β^n) to derive a relationship between the expected value of the generalized Hessian E(ℋβ,n,m), the information matrix 𝒥β,n,m, and the generalized variance-covariance matrix 𝒱β,n,m. This approximation is
(98)(∇lβ,n,m)β≈(∇lβ,n,,m)β^n+E(ℋβ,n,m)Δβ=E(ℋβ,n,m)Δβ.
Since E(∇lβ,n,m)=0 and E(ℋβ,n,m) has full rank, E(β^n)=β. We multiply each side by its own transpose and take expected values to obtain
(99)E[(∇lβ,n,m)β(∇lβ,n,m)βT]=(E(ℋβ,n,m))E[Δβ(Δβ)T](E(ℋβ,n,m))T.
Hence,
(100)𝒥β,n,m=(E(ℋβ,n,m))𝒱β,n,m(E(ℋβ,n,m))T=(E(ℋβ,n,m))𝒱β,n,m(E(ℋβ,n,m)),
so
(101)𝒱β,n,m=𝒱β,n,m(p,g,q,Xn,m;B)=(E(ℋβ,n,m))-1𝒥β,n,m(E(ℋβ,n,m))-1=ν(p,g,q;B)(ζ(p,g,q;B))2(Xn,mTXn,m)-1.
Here, we require Xn,m to have full rank and that ζ(p,g,q;B)≠0, which is certainly true when g=g1.


## 4.8. Generalized Statistical Expressions and Relations

We have shown that the expected value of the generalized HessianE(ℋβ,n,m) and the information matrix 𝒥β,n,m are multiples of Xn,mTXn,m. Specifically,
(102)E(ℋβ,n,m)=Xn,mTE(ℋz,n)Xn,m=E(∂2l(z1;p,g,q)∂z12)Xn,mTXn,m=ζ(p,g,q;B)Xn,mTXn,m,𝒥β,n,m=E((∂l(z1;p,g,q)∂z1)2)Xn,mTXn,m=ν(p,g,q;B)Xn,mTXn,m
(see (94) and (90)). Hence, if ν(p,g,q;B)≠0 (true if g=g1), then
(103)E(ℋβ,n,m)=E(∂2l(z1;p,g,q)/∂z12)E((∂l(z1;p,g,q)/∂z1)2)𝒥β,n,m=ζ(p,g,q;B)ν(p,g,q;B)𝒥β,n,m.
In addition, assuming that the scalar ζ(p,g,q;B) is also non-zero (true if g=g1) and that Xn,m has full rank m, we have proved the following relations:
(104)(𝒱β,n,m(p,g,q,Xn,m;B))-1=(ζ(p,g,q;B))2ν(p,g,q;B)Xn,mTXn,m(105)=ζ(p,g,q;B)ν(p,g,q;B)E(ℋβ,n,m)(106)=(ζ(p,g,q;B))2(ν(p,g,q;B))2𝒥β,n,m,
to be used in the derivation of the generalized log-likelihood ratio statistic (Section 4.11).


## 4.9. Statistical Relations for the Laplace Distribution

Ifg(z,q)=g1(z)=1, then using (21), (84), and (95):
(107)Q(p,g1;B)=2(1-e-pB),ν(p,g1;B)=p2,ζ(p,g1;B)=-2p2Q(p,g1;B)=-p2(1-e-pB).
Hence,
(108)𝒥β,n,m=p2Xn,mTXn,m,E(ℋβ,n,m)=(-p2)(1-e-pB)Xn,mTXn,m,
so
(109)𝒱β,n,m=(1-e-pB)2p2(Xn,mTXn,m)-1.
Here,
(110)𝒥β,n,m=-(1-e-pB)E(ℋβ,n,m),(111)𝒱β,n,m=(1-e-pB)2𝒥β,n,m-1,
where p>0, B>0, and Xn,m has full rank m. Hence, (75), derived for 𝒞2 distributions, does not hold for the Laplace distribution with bounded support. However, (75) describes the limiting behaviour, as B→∞.


## 4.10. Statistical Relations for a Laplace Distribution with Added Kurtosis

Now consider our motivating example, a Laplace distribution with bounded support[-1,1] amended by adding kurtosis. In this case g(z,q)=g2(z;q)=1+qH3(z) (see (14)), B=1,
(112)Q(p,q)=2[(p3-3qp2+6q)-e-p(p3(1-2q)+6pq+6q)]p3,f(z;p,q)=p4e(-p|z|)[1+qH3(|z|)]2[(p3-3qp2+6q)-e-p(p3(1-2q)+6pq+6q)].
Also, using (82) and (81), we obtain
(113)F(z;p,q)=F(z;p,g2,q)=-p+3q(z2-1)1+q(z3-3z),ν(p,g2,q;1)=2∫01‍(F(z;p,q))2f(z;p,q)dz,
and using (93), we obtain
(114)ζ(p,g2,q;1)=-2pf(0;p,q)+2f(0;p,q)×1g(0;q)(dg(|z|;q)d|z|)z=0+∫-11‍1g(|z|;q)d2g(|z|;q)d|z|2f(z;p,q)dz-∫-11‍1(g(|z|;q))2(dg(|z|;q)d|z|)2×f(z;p,q)dz=-2p2Q(p,q)+-6pqQ(p,q)+2∫01‍1g(z;q)d2g(z;q)dz2f(z;p,q)dz-2∫01‍1(g(z;q))2(dg(z;q)dz)2×f(z;p,q)dz=-2p2-6pqQ(p,q)+2∫01‍6qz(1+q(z3-3z))f(z;p,q)dz-2∫01‍1(1+q(z3-3z)2)×(3q(z2-1))2f(z;p,q)dz.
For typical values p=5.254 and q=0.025, numerical integration gives ν=28.3561 and ζ=-28.4957, to four decimal places. In this case,
(115)𝒥β,n,m=28.3561Xn,mTXn,m,E(ℋβ,n,m)=-28.4957Xn,mTXn,m,
and so
(116)𝒱β,n,m=28.3561(-28.4957)2(Xn,mTXn,m)-1.
Here
(117)𝒥β,n,m=28.3561-28.4957E(ℋβ,n,m),𝒱β,n,m=(28.3561)2(-28.4957)2𝒥β,n,m-1.
In this example, (75) could be used as an approximation but does not exactly describe the relationship between E(ℋβ,n,m), 𝒥β,n,m, and 𝒱β,n,m.


## 4.11. The Generalized Log-Likelihood Ratio Statistic

The log-likelihood ratio statistic enables us to assess the adequacy of a model. It enables us to compare a model withM coefficients (parameters) with a model of interest which differs only in that it has fewer coefficients, say P, with M>P≥1, see [1], for example. We wish to compare a linear model with m coefficients, the βj, for j=1,2,…,m with a lesser linear model with fewer coefficients. The aim is to decide whether or not the excluded coefficients are useful. Our comparison is conditional on the parameters p and q.LetLρ,n,M(ρ1,…,ρM;p,q;Xn,M,y) denote the likelihood function for the model with M parameters, and let Lψ,n,P(ψ1,…,ψP;p,q;Xn,P,y) denote the likelihood function for the lesser model. Let ρ^ (or ρ^n) be the maximum likelihood estimator of ρ=(ρ1,…,ρM)T, and let ψ^ (or ψ^n) be the maximum likelihood estimator of ψ=(ψ1,…,ψP)T. Then the likelihood ratio
(118)λ=Lρ,n,M(ρ^;p,q;Xn,M,y)Lψ,n,P(ψ^;p,q;Xn,P,y)
is a ratio of two probabilities and will be greater than one since the model with M parameters provides the more complete description of the model. In our application, M=m, ρ=β, and usually (but not necessarily) P=m-1. We show that a multiple (conditional on p and q) of
(119)log(λ)=log(Lρ,n,M(ρ^;p,q;Xn,M,y))-log(Lψ,n,P(ψ^;p,q;Xn,P,y))
has a chi-squared distribution as follows.The derivation of the log-likelihood ratio statistic (for smooth functions) may be found in the textbooks. For example, for generalized linear models (see [1]) where the log-likelihood function is smooth in a neighbourhood of its maximum, D=2log(λ) is distributed approximately as χ2(M-P,δD). Here δD is a noncentrality parameter, a positive constant which will be near zero if the lesser model fits the data almost as well as the model with more coefficients. Consider our situation in which the log-likelihood function is continuous at a maximum, but where this maximum occurs at a vertex or possibly on a ridge in (β,l)-space, and generalized calculus is required to consider Taylor series expansions. Then, as in (97), around ρ^, replacing the Hessian of log(LM)=lM by its expected value, we obtain the following approximation:
(120)log(Lρ,n,M(E(ρ^);p,q;Xn,M,y))-log(Lρ,n,M(ρ^;p,q;Xn,M,y))=log(Lρ,n,M(E(ρ^)))-log(Lρ,n,M(ρ^))≈12(E(ρ^)-ρ^)TE(ℋρ,n,M)(E(ρ^)-ρ^).So, by (105) assuming that the scalar factors ν(p,g,q;B) and ζ(p,g,q;B) are both non-zero and that Xn,M has full rank
(121)2E((∂2l(z1;p,g,q))/∂z12)E(((∂l(z1;p,g,q))/∂z1)2)(lρ,n,M(E(ρ^))-lρ,n,M(ρ^))≈(E(ρ^)-ρ^)Tζ(p,g,q;B)ν(p,g,q;B)E(ℋρ,n,M)(E(ρ^)-ρ^)=(E(ρ^)-ρ^)T𝒱ρ,n,M-1(E(ρ^)-ρ^),
which has the distribution χ2(M) if the MLE has a normal distribution. We show in Section 5 that the distribution of the MLE is asymptotically normal. Hence, if n is large enough, this will be a good approximation. Similarly
(122)2E((∂2l(z1;p,g,q))/∂z12)E(((∂l(z1;p,g,q))/∂z1)2)(lψ,n,P(E(ψ^))-lψ,n,P(ψ^))≈(E(ψ^)-ψ^)Tζ(p,g,q;B)ν(p,g,q;B)E(ℋψ,n,P)(E(ψ^)-ψ^)=(E(ψ^)-ψ^)T𝒱ψ,n,P-1(E(ψ^)-ψ^),
which has the distribution χ2(P), approximately. Noting that
(123)ζ(p,g,q;B)ν(p,g,q;B)=-1(1-e-pB)<0,
when g=g1, let
(124)Dgen=-2ζ(p,g,q;B)ν(p,g,q;B)log(λ).
Then Dgen, the log-likelihood ratio statistic calculated with generalized functions, a positive number, may be expressed as
(125)Dgen=-2ζ(p,g,q;B)ν(p,g,q;B)×(lρ,n,M(ρ^;p,q;Xn,M,y)-lψ,n,P(ψ^;p,q;Xn,P,y))=-2ζ(p,g,q;B)ν(p,g,q;B)×(lρ,n,M(ρ^;p,q;Xn,M,y)-lρ,n,M(E(ρ^);p,q;Xn,M,y))+2ζ(p,g,q;B)ν(p,g,q;B)×(lψ,n,P(ψ^;p,q;Xn,P,y)-lψ,n,P(E(ψ^);p,q;Xn,P,y))-2ζ(p,g,q;B)ν(p,g,q;B)×(lρ,n,M(E(ρ^);p,q;Xn,M,y)-lψ,n,P(E(ψ^);p,q;Xn,P,y)).
Hence, Dgen is the sum of three terms; the first (positive) has the distribution χ2(M) (approximately). The second (negative) has the distribution χ2(P) (approximately). The third is a positive constant (say δgen) that depends on p and q. Hence, Dgen is distributed approximately as χ2(M-P,δgen). If the lesser model gives a good description of the data, then δgen will be small. The generalized log-likelihood ratio statistic Dgen is easily calculated and is hence a potentially useful statistic for assessing our linear model for which the log-likelihood function is not 𝒞2. Note if g=g1, Dgen=(2log(λ))/(1-e-pB), and so Dgen→2log(λ), as B→∞.


## 5. The Maximum Likelihood Estimator Is Consistent and Asymptotically Normal

The generalized expressions derived in Section4 will be used to prove the asymptotic convergence of our MLE to a random variable with a normal distribution. Recall the following assumptions.  (i)
Our model is linear (see (4)).(ii)
The response functionf is a Laplace probability density function, generally perturbed and/or truncated, as given in (10).We make the following further assumptions.  (i)
There exists a unique true vector of coefficientsβ0∈ℝm whose value we are trying to estimate.(ii)
The matrixXn,m has full rank m, so that Xn,mTXn,m has full rank m.(iii)
Thelimn→∞(1/n)Xn,mTXn,m=Wm, a positive definite matrix.(iv)
Assuming fixedm, for n≥m, denote by β^n a (not necessarily unique) MLE of the true value β0, corresponding to the explanatory variables in Xn,m.Lemma 11.
The ML estimatesβ^n exist, for n≥m.Proof of Lemma11.
Since a continuous function on a compact set attains it maximum, the existence of a maximum of the log-likelihood functionlβ,n,m is guaranteed for finite bound B. Even if the domain is not bounded, we can work with finite bound B, and then let B→∞.Theorem 12.
The random variablen(β^n-β0) converges in distribution to an m-dimensional normally distributed random vector with mean 0 and covariance matrix (ν/ζ2)Wm-1, that is, as n→∞,
(126)n(β^n-β0)⟶DN(0,(νζ2)Wm-1).Lemma 13.
The random variable(1/n)∇lβ,n,m converges in distribution to an m-dimensional normally distributed random vector with mean 0 and covariance matrix νWm, that is, as n→∞,
(127)(1n)∇lβ,n,m⟶DN(0,νWm).
Hence, (1/n)∇lβ,n,m converges in probability to 0.Proof of Lemma13.
Consider the random variable∇lβ,n,m=-Xn,mT∇lz,n. If we repeat the sampling of n≥m data points, using Xn,m a total of t times, we may write
(128)∇lβ,nt,m=-Xnt,mT∇lz,nt,=-tXn,mT∇lz,n¯,
where ∇lz,n¯ denotes the average of t samples of the random vector ∇lz,n. Now ∇lz,n has mean 0 and covariance matrix νIn, and so by the multivariate Central Limit Theorem, as t→∞,
(129)t∇lz,n¯⟶DN(0,νIn),Xn,mTt∇lz,n¯⟶DN(0,Xn,mTνXn,m),(1t)∇lβ,nt,m⟶DN(0,Xn,mTνXn,m),(1nt)∇lβ,nt,m⟶DN(0,νnXn,mTXn,m).
Hence, as n→∞ and t→∞,
(130)(1nt)∇lβ,nt,m⟶DN(0,νWm).
Hence, as n→∞,
(131)(1n)∇lβ,n,m⟶DN(0,νWm).Lemma 14.
Although for finiten, the log-likelihood function lβ,n,m is not differentiable at a maximum; (1/n)lβ,n,m converges in distribution to a negative definite quadratic function centred at β0. Hence, the MLEs β^n are consistent, that is, they converge in probability to the true value β0.Proof of Lemma14.
Using a generalized Taylor series expansion aboutβ0 and noting that by Lemma 13, (1/n)∇lβ,n,m converges in probability to 0, we can write
(132)limn→∞lβ,n,m(β)-lβ,n,m(β0)n=limn→∞(β-β0)TE(ℋβ,n,m)n(β-β0)=(β-β0)T(limn→∞ζXn,mTXn,mn)(β-β0)=ζ(β-β0)TWm(β-β0).
This shows that, in the limit as n→∞, the log-likelihood function lβ,n,m has an isolated maximum at β0, and so a sequence of MLEs β^n must converge to β0 in probability (consistency).In the Proof of Lemma14 we ignored the third partial derivatives of lβ,n,m in the generalized Taylor series expansion. The justification is as follows. Since g is assumed to be 𝒞3 on an open interval which contains [0,B], |d3g/dz3| must be bounded above by some positive real number G3 on [0,B]. Hence, the absolute values of the third partial derivatives of l must be bounded above by some positive real number except at points z, where some zi=0 (data points). Note that E(∂3l/∂zi3)=0, i=1,2,…,n. This follows from (66), using the symmetry introduced by the modulus function and (62). Hence, the third partial derivative terms will be small compared to the second partial derivative terms near a critical point, and so we may ignore them in our generalized Taylor series expansion for l, when the expected value of the Hessian has full rank.Proof of Theorem12.
Consider the first degree approximation(133)∇lβ,n,m(β)~E(ℋβ,n,m)(β-β0).
Since E(ℋβ,n,m) has full rank,
(134)(β-β0)~(E(ℋβ,n,m))-1∇lβ,n,m(β),
so
(135)nn(β-β0)~(E(ℋβ,n,m)n)-1∇lβ,n,m(β)n.
By Lemma 13, as n→∞,
(136)(1n)∇lβ,n,m⟶DN(0,νWm).
Hence, as n→∞,
(137)n(β-β0)⟶DN(0,(ζWm)-1νWm(ζWm)-T).
Hence, as n→∞,
(138)n(β-β0)⟶DN(0,(νζ2)Wm-1).An alternative Proof of Theorem12 for the simplest case, that is, in the absence of truncation or perturbation, (i.e., g=g1 and B=∞) is in [2]. This case is also discussed in [11, page 451].


## 6. Real and Simulated Data Illustrations

### 6.1. Empirical Distribution of Methylation Proportion Deviations

Quantitative analysis of DNA methylation at specific genomic sites (known as CpG sites) was carried out with the Sequenom MassARRAY Compact System (http://www.sequenom.com/). Briefly, this involves accurate determination and comparison of the mass of transcription cleavage products following chemical modification of the DNA which is dependent upon the a priori methylation status, using MALDI-TOF mass spectrometry (Bruker-Sequenom) (see [12]). Quantitative CpG methylation was assessed using proprietary EpiTyper software v1.0.5 (http://www.sequenom.com/). Sequenom measurements of 1440 CpG sites in each of 41 human umbilical cord tissue samples were performed in duplicate, and the difference between the measurements was recorded. This difference represents the deviation in the measurement of CpG methylation due to sample nanodispensing and MALDI-TOF mass spectrometry detection. Figure 4 is a histogram of the deviations of these 1440 repeated measurements. Although not obvious from the histogram, about 1.04% of values are greater than 0.2 in absolute value. The histogram has heavy tails. The data was shown not to conform to a normal distribution using the test proposed by [13] (P value <0.00001).Figure 4
Deviations in methylation proportion.


### 6.2. Simulated Data Example Using Methylation Proportion Deviations

In order to illustrate the application of the theory developed, a sample of 40 deviations was chosen at random from the total pool of 1440 available CpG methylation proportion deviations. A constant value of 0.48 was added to 20 of these samples and designated treatment H, while a constant value of 0.45 was added to the other 20 samples and designated treatment L. A uniform random variable sampled between −0.01 and 0.01 was added to each value to simulate the additional differences expected to occur between individuals.We analysed the data using our amended Laplace distribution (13), with g=g2, p=53.41, q=0.0314 (machine characteristics), and B=1. The parameter MLEs p=53.41±3.45 and q=0.0314±0.0097 (estimates plus or minus standard errors) were given by Buckland's algorithm [14, 15], with the q being highly significantly different from zero, (with a log-likelihood ratio statistic of 135.6 distributed as χ2(1)). Assuming q=0.0 (LAE regression with truncation) yields p=36.22±1.167. Buckland's algorithm does not take into account the truncation, but since p is large and (1-e-pB) is so close to 1 in this example, the effect of truncation on the standard errors and Dgen (for P values) is negligible. If the machine characteristic p was smaller, say 1≤p≤3, assuming B=1, this effect would become more significant.We coded a low value treatment (L) by settingxi2=-1, i=1,2,…,20 and a high value treatment (H) by setting xi2=1, i=21,22,…,40. Estimates for β1 and β2 were calculated by maximum likelihood estimation, using the simplex method. The standard errors of the coefficient estimates were calculated using our generalized 𝒱. Hence, the estimated treatment means (β1±β2), and their standard errors were calculated. A P value was obtained using Dgen, which is assumed distributed χ2(1). This simulation was performed twice, and the results are given in Table 1. Note for each simulation, the P value is less than 0.01, indicating a significant difference between the means.Table 1
Two simulation results, adding high (H) or low (L) treatments (T) to DNA methylation proportions.
T
Amended Laplace distribution
Laplace distribution
Normal distribution
p
=
53.41, q=0.0314
p
=
36.22, q=0.0
MATLAB glmfit
Mean
Std. err.
P value
Mean
Std. err.
P value
Mean
Std. err.
P value
H
0.4817
0.0042
2
e
-
05
0.4812
0.0062
4
e
-
04
0.4658
0.0094
4
e
-
01
L
0.4532
0.0042
0.4519
0.0062
0.4540
0.0094
H
0.4803
0.0042
2
e
-
04
0.4809
0.0062
2
e
-
03
0.4887
0.0119
2
e
-
01
L
0.4592
0.0042
0.4586
0.0062
0.4673
0.0119For both simulationsν(p,g2,q;1)=2862.70 and ζ(p,g2,q;1)=-2862.71 (see Section 4.10). Also,
(139)XTX=(400040),
rounding to five significant figures,
(140)E(ℋ)=(-11451000-114510),𝒥=(11451000114510),𝒱=(8.7330e-06008.7330e-06).
Here,
(141)E(ℋ)=-1.0000𝒥,𝒱-1=1.0000𝒥,Dgen=2.0000log(λ).
These expressions correspond to the classical expressions for smooth functions, up to our numerical tolerance.For comparison, the results of a standard analysis of variance, assuming the deviations have a normal distribution, are also given in Table1, for each of the two simulations. Using the P values obtained, 4e-01 and 2e-01, we would not reject the null hypothesis (that the means are equal) at any usual level of significance, using a least squares approach. However, our algorithm based on the amended Laplace distribution correctly identified the structure of the simulated data set, separating two means which were fairly close. The standard least squares algorithm failed to do this. The two data sets analysed are given in Tables 3 and 4 in Appendix B.For comparison we also included LAE regression, estimating the model coefficients assuming the response function is the Laplace probability density function truncated to[-1,1], using p=36.22, given by Buckland's algorithm. We see that including the perturbing Hermite polynomial improves (decreases) the P values, meaning we can be even more confident, than when using LAE regression, that the means are different.


### 6.3. Primiparous versus Multiparous Effects on DNA Methylation Proportion at the Promoter of the H19 Gene

The CpG methylation at two CpG sites in the promoter of the H19 gene was measured by Sequenom in umbilical cord tissues collected as part of an ongoing prospective birth cohort study. Phenotype variables in this population include birth order or parity, defined as first born child (primiparous) or later born (multiparous). We analysed the relationship between H19 gene methylation status and birth order in this study, using our amended Laplace distribution (13), with g=g2 and B=1, and for comparison, the usual least squares algorithm. Here an amended Laplace distribution (13) has been used to model the measurement errors as opposed to the differences of errors in repeated measurements (the deviations in Section 6.2). In practice this works well, despite the observation that if the measurement errors have a distribution of the form (13), then the deviations would not be expected to have exactly this form of distribution. Modelling the deviations as the difference between two random variables both distributed Laplace also works well. We used the estimate of pd=36.22 for the deviations to estimate pe=51.22~36.22*2 for the errors. The factor 2 is chosen to halve the variance in the underlying Laplace distribution. For the error distribution, we estimated qe to be as close to 0.5 as we allow, say qe=0.4999. This was done by assuming the errors were distributed as an amended Laplace distribution (13) with parameters pe and qe, deriving a corresponding distribution for the deviations and using MLE with the deviation data to obtain parameter estimates. We also tried pe=75.53~53.41*2 with qe=0.4999. The latter (Sequenom) parameter estimates gave smaller standard errors and P values and were consistent with Infinium parameter estimates.The problem was coded by substitutingxi2=1 for primiparous and xi2=-1 for multiparous. Estimates for β1 and β2 and their standard errors were calculated. The data used is given in Table 5 in Appendix B. The estimated means (β1±β2) and their standard errors are given in Table 2. The small P values associated with the amended Laplace distribution (13), calculated using Dgen, identify a difference between the mean methylation proportions (primiparous versus multiparous) at a given CpG site. This demonstrates the power of accounting properly for the distributional properties of the methylation errors and hence enables clearer inference of the epigenetic mechanisms underlying these biological phenomena.Table 2
Primiparous (p) versus multiparous (m) effects on DNA methylation proportion at the promoter of the H19 gene.
Site
Parity
Amended Laplace distributionp=75.53, q=0.4999
Normal distribution
Mann-Whitney
Mean
Std. err.
P value
Mean
Std. err.
P value
P value
CpG9
p
0.180
0.0029
<1.0e-09
0.300
0.059
1
e
-
01
5
e
-
02
CpG9
m
0.450
0.0029
0.441
0.059
CpG13
p
0.230
0.0029
<1.0e-09
0.326
0.061
4
e
-
02
4
e
-
02
CpG13
m
0.560
0.0029
0.523
0.061Table 3
First simulation data, treatments either H (xi=1) or L (xi=-1) plus randomly sampled methylation deviances and randomly sampled uniformly distributed individual variation.
T
y
i
T
y
i
T
y
i
T
y
i
L
0.4579
L
0.4467
H
0.4841
H
0.4873
L
0.4243
L
0.4610
H
0.4735
H
0.4761
L
0.4993
L
0.4609
H
0.4878
H
0.2391
L
0.4131
L
0.4851
H
0.4823
H
0.4805
L
0.4463
L
0.4340
H
0.4462
H
0.4877
L
0.4317
L
0.4573
H
0.4664
H
0.4779
L
0.4473
L
0.4360
H
0.4845
H
0.4929
L
0.4347
L
0.4584
H
0.4817
H
0.4863
L
0.4760
L
0.4420
H
0.4861
H
0.4751
L
0.4776
L
0.4914
H
0.4668
H
0.4543Table 4
Second simulation data, treatments either H (xi=1) or L (xi=-1) plus randomly sampled methylation deviances and randomly sampled uniformly distributed individual variation.
T
y
i
T
y
i
T
y
i
T
y
i
L
0.5287
L
0.4416
H
0.4881
H
0.4829
L
0.5224
L
0.4547
H
0.4803
H
0.4246
L
0.5162
L
0.4496
H
0.4789
H
0.4609
L
0.4564
L
0.4568
H
0.4790
H
0.4233
L
0.4628
L
0.4574
H
0.4739
H
0.5412
L
0.5230
L
0.4599
H
0.4974
H
0.5193
L
0.3731
L
0.5124
H
0.7010
H
0.4921
L
0.4389
L
0.4592
H
0.4725
H
0.4702
L
0.4519
L
0.4458
H
0.4871
H
0.5520
L
0.4675
L
0.4685
H
0.4885
H
0.3600Table 5
CpG methylation measurements at sites 9 and 13 on the promoter of the H19 gene versus primiparous (p) or multiparous (m).
CpG9
p/m
CpG9
p/m
CpG13
p/m
CpG13
p/m
1.00
p
0.16
p
0.30
p
0.16
p
0.08
p
0.19
p
0.00
p
0.36
p
0.04
p
0.15
p
0.03
p
0.02
p
0.17
p
0.35
p
0.25
p
0.60
p
0.46
p
0.04
p
0.80
p
0.01
p
1.00
p
0.27
p
0.71
p
0.70
p
0.18
p
0.32
p
0.17
p
0.56
p
0.33
m
0.37
p
0.56
m
0.70
p
0.28
m
0.05
p
0.40
m
0.00
p
0.82
m
0.39
p
0.57
m
0.61
p
0.20
p
0.07
p
0.18
p
0.02
p
0.08
p
0.17
p
0.03
p
0.23
p
0.15
p
0.14
m
0.09
p
0.99
m
1.00
p
0.61
m
0.96
p
0.60
m
0.10
m
0.53
m
0.79
m
0.35
m
0.89
m
0.45
m
0.83
m
0.63
m
0.07
m
0.09
m
0.02
m
0.07
m
0.62
m
0.57
m
0.38
m
0.53
m
0.48
m
0.73
m
0.68
m
0.72
m
0.31
m
0.30
m
0.27
m
0.22
m
0.62
m
0.80
mIn this example, settingq=0.0 (LAE regression with truncation) yields similar means and variances to our MLE algorithm (the same to two decimal places) and also gives small P values. For the CpG13 example with p=75.53 and q=0.4999, Dgen=201.479. For this example with p=75.53 and q=0.0, Dgen=197.889. These values are very close as pB is large and so e-pB is small. For smaller values of pB, the value of taking into account the truncation and perturbation increases. A Mann-Whitney 2-sample U-test was also done for comparison.


## 6.1. Empirical Distribution of Methylation Proportion Deviations

Quantitative analysis of DNA methylation at specific genomic sites (known as CpG sites) was carried out with the Sequenom MassARRAY Compact System (http://www.sequenom.com/). Briefly, this involves accurate determination and comparison of the mass of transcription cleavage products following chemical modification of the DNA which is dependent upon the a priori methylation status, using MALDI-TOF mass spectrometry (Bruker-Sequenom) (see [12]). Quantitative CpG methylation was assessed using proprietary EpiTyper software v1.0.5 (http://www.sequenom.com/). Sequenom measurements of 1440 CpG sites in each of 41 human umbilical cord tissue samples were performed in duplicate, and the difference between the measurements was recorded. This difference represents the deviation in the measurement of CpG methylation due to sample nanodispensing and MALDI-TOF mass spectrometry detection. Figure 4 is a histogram of the deviations of these 1440 repeated measurements. Although not obvious from the histogram, about 1.04% of values are greater than 0.2 in absolute value. The histogram has heavy tails. The data was shown not to conform to a normal distribution using the test proposed by [13] (P value <0.00001).Figure 4
Deviations in methylation proportion.


## 6.2. Simulated Data Example Using Methylation Proportion Deviations

In order to illustrate the application of the theory developed, a sample of 40 deviations was chosen at random from the total pool of 1440 available CpG methylation proportion deviations. A constant value of 0.48 was added to 20 of these samples and designated treatment H, while a constant value of 0.45 was added to the other 20 samples and designated treatment L. A uniform random variable sampled between −0.01 and 0.01 was added to each value to simulate the additional differences expected to occur between individuals.We analysed the data using our amended Laplace distribution (13), with g=g2, p=53.41, q=0.0314 (machine characteristics), and B=1. The parameter MLEs p=53.41±3.45 and q=0.0314±0.0097 (estimates plus or minus standard errors) were given by Buckland's algorithm [14, 15], with the q being highly significantly different from zero, (with a log-likelihood ratio statistic of 135.6 distributed as χ2(1)). Assuming q=0.0 (LAE regression with truncation) yields p=36.22±1.167. Buckland's algorithm does not take into account the truncation, but since p is large and (1-e-pB) is so close to 1 in this example, the effect of truncation on the standard errors and Dgen (for P values) is negligible. If the machine characteristic p was smaller, say 1≤p≤3, assuming B=1, this effect would become more significant.We coded a low value treatment (L) by settingxi2=-1, i=1,2,…,20 and a high value treatment (H) by setting xi2=1, i=21,22,…,40. Estimates for β1 and β2 were calculated by maximum likelihood estimation, using the simplex method. The standard errors of the coefficient estimates were calculated using our generalized 𝒱. Hence, the estimated treatment means (β1±β2), and their standard errors were calculated. A P value was obtained using Dgen, which is assumed distributed χ2(1). This simulation was performed twice, and the results are given in Table 1. Note for each simulation, the P value is less than 0.01, indicating a significant difference between the means.Table 1
Two simulation results, adding high (H) or low (L) treatments (T) to DNA methylation proportions.
T
Amended Laplace distribution
Laplace distribution
Normal distribution
p
=
53.41, q=0.0314
p
=
36.22, q=0.0
MATLAB glmfit
Mean
Std. err.
P value
Mean
Std. err.
P value
Mean
Std. err.
P value
H
0.4817
0.0042
2
e
-
05
0.4812
0.0062
4
e
-
04
0.4658
0.0094
4
e
-
01
L
0.4532
0.0042
0.4519
0.0062
0.4540
0.0094
H
0.4803
0.0042
2
e
-
04
0.4809
0.0062
2
e
-
03
0.4887
0.0119
2
e
-
01
L
0.4592
0.0042
0.4586
0.0062
0.4673
0.0119For both simulationsν(p,g2,q;1)=2862.70 and ζ(p,g2,q;1)=-2862.71 (see Section 4.10). Also,
(139)XTX=(400040),
rounding to five significant figures,
(140)E(ℋ)=(-11451000-114510),𝒥=(11451000114510),𝒱=(8.7330e-06008.7330e-06).
Here,
(141)E(ℋ)=-1.0000𝒥,𝒱-1=1.0000𝒥,Dgen=2.0000log(λ).
These expressions correspond to the classical expressions for smooth functions, up to our numerical tolerance.For comparison, the results of a standard analysis of variance, assuming the deviations have a normal distribution, are also given in Table1, for each of the two simulations. Using the P values obtained, 4e-01 and 2e-01, we would not reject the null hypothesis (that the means are equal) at any usual level of significance, using a least squares approach. However, our algorithm based on the amended Laplace distribution correctly identified the structure of the simulated data set, separating two means which were fairly close. The standard least squares algorithm failed to do this. The two data sets analysed are given in Tables 3 and 4 in Appendix B.For comparison we also included LAE regression, estimating the model coefficients assuming the response function is the Laplace probability density function truncated to[-1,1], using p=36.22, given by Buckland's algorithm. We see that including the perturbing Hermite polynomial improves (decreases) the P values, meaning we can be even more confident, than when using LAE regression, that the means are different.


## 6.3. Primiparous versus Multiparous Effects on DNA Methylation Proportion at the Promoter of the H19 Gene

The CpG methylation at two CpG sites in the promoter of the H19 gene was measured by Sequenom in umbilical cord tissues collected as part of an ongoing prospective birth cohort study. Phenotype variables in this population include birth order or parity, defined as first born child (primiparous) or later born (multiparous). We analysed the relationship between H19 gene methylation status and birth order in this study, using our amended Laplace distribution (13), with g=g2 and B=1, and for comparison, the usual least squares algorithm. Here an amended Laplace distribution (13) has been used to model the measurement errors as opposed to the differences of errors in repeated measurements (the deviations in Section 6.2). In practice this works well, despite the observation that if the measurement errors have a distribution of the form (13), then the deviations would not be expected to have exactly this form of distribution. Modelling the deviations as the difference between two random variables both distributed Laplace also works well. We used the estimate of pd=36.22 for the deviations to estimate pe=51.22~36.22*2 for the errors. The factor 2 is chosen to halve the variance in the underlying Laplace distribution. For the error distribution, we estimated qe to be as close to 0.5 as we allow, say qe=0.4999. This was done by assuming the errors were distributed as an amended Laplace distribution (13) with parameters pe and qe, deriving a corresponding distribution for the deviations and using MLE with the deviation data to obtain parameter estimates. We also tried pe=75.53~53.41*2 with qe=0.4999. The latter (Sequenom) parameter estimates gave smaller standard errors and P values and were consistent with Infinium parameter estimates.The problem was coded by substitutingxi2=1 for primiparous and xi2=-1 for multiparous. Estimates for β1 and β2 and their standard errors were calculated. The data used is given in Table 5 in Appendix B. The estimated means (β1±β2) and their standard errors are given in Table 2. The small P values associated with the amended Laplace distribution (13), calculated using Dgen, identify a difference between the mean methylation proportions (primiparous versus multiparous) at a given CpG site. This demonstrates the power of accounting properly for the distributional properties of the methylation errors and hence enables clearer inference of the epigenetic mechanisms underlying these biological phenomena.Table 2
Primiparous (p) versus multiparous (m) effects on DNA methylation proportion at the promoter of the H19 gene.
Site
Parity
Amended Laplace distributionp=75.53, q=0.4999
Normal distribution
Mann-Whitney
Mean
Std. err.
P value
Mean
Std. err.
P value
P value
CpG9
p
0.180
0.0029
<1.0e-09
0.300
0.059
1
e
-
01
5
e
-
02
CpG9
m
0.450
0.0029
0.441
0.059
CpG13
p
0.230
0.0029
<1.0e-09
0.326
0.061
4
e
-
02
4
e
-
02
CpG13
m
0.560
0.0029
0.523
0.061Table 3
First simulation data, treatments either H (xi=1) or L (xi=-1) plus randomly sampled methylation deviances and randomly sampled uniformly distributed individual variation.
T
y
i
T
y
i
T
y
i
T
y
i
L
0.4579
L
0.4467
H
0.4841
H
0.4873
L
0.4243
L
0.4610
H
0.4735
H
0.4761
L
0.4993
L
0.4609
H
0.4878
H
0.2391
L
0.4131
L
0.4851
H
0.4823
H
0.4805
L
0.4463
L
0.4340
H
0.4462
H
0.4877
L
0.4317
L
0.4573
H
0.4664
H
0.4779
L
0.4473
L
0.4360
H
0.4845
H
0.4929
L
0.4347
L
0.4584
H
0.4817
H
0.4863
L
0.4760
L
0.4420
H
0.4861
H
0.4751
L
0.4776
L
0.4914
H
0.4668
H
0.4543Table 4
Second simulation data, treatments either H (xi=1) or L (xi=-1) plus randomly sampled methylation deviances and randomly sampled uniformly distributed individual variation.
T
y
i
T
y
i
T
y
i
T
y
i
L
0.5287
L
0.4416
H
0.4881
H
0.4829
L
0.5224
L
0.4547
H
0.4803
H
0.4246
L
0.5162
L
0.4496
H
0.4789
H
0.4609
L
0.4564
L
0.4568
H
0.4790
H
0.4233
L
0.4628
L
0.4574
H
0.4739
H
0.5412
L
0.5230
L
0.4599
H
0.4974
H
0.5193
L
0.3731
L
0.5124
H
0.7010
H
0.4921
L
0.4389
L
0.4592
H
0.4725
H
0.4702
L
0.4519
L
0.4458
H
0.4871
H
0.5520
L
0.4675
L
0.4685
H
0.4885
H
0.3600Table 5
CpG methylation measurements at sites 9 and 13 on the promoter of the H19 gene versus primiparous (p) or multiparous (m).
CpG9
p/m
CpG9
p/m
CpG13
p/m
CpG13
p/m
1.00
p
0.16
p
0.30
p
0.16
p
0.08
p
0.19
p
0.00
p
0.36
p
0.04
p
0.15
p
0.03
p
0.02
p
0.17
p
0.35
p
0.25
p
0.60
p
0.46
p
0.04
p
0.80
p
0.01
p
1.00
p
0.27
p
0.71
p
0.70
p
0.18
p
0.32
p
0.17
p
0.56
p
0.33
m
0.37
p
0.56
m
0.70
p
0.28
m
0.05
p
0.40
m
0.00
p
0.82
m
0.39
p
0.57
m
0.61
p
0.20
p
0.07
p
0.18
p
0.02
p
0.08
p
0.17
p
0.03
p
0.23
p
0.15
p
0.14
m
0.09
p
0.99
m
1.00
p
0.61
m
0.96
p
0.60
m
0.10
m
0.53
m
0.79
m
0.35
m
0.89
m
0.45
m
0.83
m
0.63
m
0.07
m
0.09
m
0.02
m
0.07
m
0.62
m
0.57
m
0.38
m
0.53
m
0.48
m
0.73
m
0.68
m
0.72
m
0.31
m
0.30
m
0.27
m
0.22
m
0.62
m
0.80
mIn this example, settingq=0.0 (LAE regression with truncation) yields similar means and variances to our MLE algorithm (the same to two decimal places) and also gives small P values. For the CpG13 example with p=75.53 and q=0.4999, Dgen=201.479. For this example with p=75.53 and q=0.0, Dgen=197.889. These values are very close as pB is large and so e-pB is small. For smaller values of pB, the value of taking into account the truncation and perturbation increases. A Mann-Whitney 2-sample U-test was also done for comparison.


## 7. Discussion

### 7.1. A Comparison with LAE Regression

The original MLE theory and methods in this paper were developed assuming the response function is a modified version of the Laplace probability density function, that is, assuming nontrivial perturbation and/or truncation to compact support[-B,B]. Such response functions have been observed in measurement data generated by nearly all the analytical platforms currently used to assess DNA methylation, including the Sequenom EpiTyper, Infinium Mass Array and Restricted Representation Bisulphite Sequencing platforms [16].In the absence of perturbation or truncation of the response function, the results in this paper correspond to the theory of LAE (or median) regression as found in [2–4]. That is, if g(z)=1 and B=∞, then our MLE method corresponds with LAE regression. In this case, Q=2, ν=p2, ζ=-p2, f(z)=(p/2)e-p|z| (in one dimension), and so
(142)νζ2=p2p4=1p2=1(2f(0))2,
which is the asymptotic variance of the ordinary sample median for f [2].We present an original and practical method of obtaining the covariance matrix for the model coefficients. This involves evaluating(Xn,mTXn,m)-1, where although n might be large, m generally is not, and using generalized functions to numerically evaluate two one-dimensional integrals (to find ν and ζ). The calculation of ν and ζ takes into account the characteristics of the response function (truncation and perturbation) which would be ignored if we used median regression. This is possible when we know the response function parameters, or have fairly accurate estimates, as in our epigenetic application, modelling DNA methylation proportion deviations.For LAE regression, other methods of determining approximations to this covariance matrix may be found in the literature. In particular, in the method of quantile regression [5] implemented in the statistical package R, the covariance matrix for the MLE is calculated by resampling techniques, by bootstrapping or by using hierarchical spline models [17].We prove that, even for truncated and perturbed Laplace response functions, subject to certain restrictions, the maximum of the log-likelihood function occurs at a data point. This result is well-known in the case of LAE regression. A proof that the LAE estimator passes through at leastrX data points may be found in references [3, 4].Three asymptotically equivalent test statistics for LAE regression may be found in [18], namely, a likelihood ratio test statistic, a Wald test statistic, and a Lagrange multiplier test statistic. Our likelihood ratio test statistic is an original modification of the former, applicable to our general case (not restricted to LAE regression), calculated using generalized functions. An F test statistic for LAE regression is found in [4].When working with a model for which the response function is assumed to be a truncated Laplace probability density function, we could ignore the truncation to[-B,B]. However, taking the truncation into account reduces the variance in the model coefficient estimates by a factor (1-e-pB)2 and increases the log-likelihood ratio statistic by a factor (1-e-pB)-1. This effect is small if e-pB is small but becomes more significant as pB decreases, that is, as more of the density function is truncated. Refer to Section 4.9, (111) and Section 4.11, (123), and (124) which show that, for example, when g=g1,
(143)Dgen=2(1-e-pB)log(λ).
Hence, by taking into account the truncation, we can be more confident of our coefficient estimates and the value of appropriate beta coefficients than the standard formulae for LAE regression indicate. This effect will also be seen for small perturbations of the density function.


### 7.2. Original Formulae Enable Model Comparison

The original formulae derived in Section4 for the covariance matrix of the model coefficients and the generalized log-likelihood ratio statistic enable us to do a generalized analysis of variance, in our case comparing the means of two different groups of data. We see that in Section 6.3, the primiparous and multiparous means could be judged distinct with very high probability using either the Laplace or an amended Laplace response function. The P-values are slightly lower when using the Hermite modification. In the simulation in Section 6.2, using the Hermite modification decreases the P-value by an order of magnitude. In a situation in which p is low or two means are less distinct, the Hermite modification may prove very useful.Preliminary results indicate the use of an amended Laplace distribution enables a clearer separation of means than that given by other more standard procedures, for example, beta regression, in cases where independent evidence suggests that the means are different [16].


### 7.3. Summary

The Laplace distribution is the basis of many mathematical models (see [19]). Our focus has been modelling the distributions of errors in the proportions of DNA methylation measured at genomic CpG sites.Molecular biology deals with complex interactions both in terms of the physiology of the processes of interest and in the instrumentation required to measure these effects. The non-linearity of these processes can result in frequency distributions that are far from normal, so that application of  “standard” methods of statistical inference based on least squares may be inadequate. Methods which deal with the form of the frequency distribution directly such as maximum likelihood are necessary for adequate inference to be made.The Laplace or double exponential distribution considered here has been observed in molecular biology studies, where a significant proportion of high deviations appears to occur regularly [16, 20]. The extension by Hermite polynomials considered here provides flexibility for describing the tails of the distribution. However, as noted, the use of the Laplace distribution as the “key” function introduces problems in finding maximum likelihood estimators and particularly their standard errors. This paper presents both a practical method for dealing with these problems and the underlying asymptotic theory.


## 7.1. A Comparison with LAE Regression

The original MLE theory and methods in this paper were developed assuming the response function is a modified version of the Laplace probability density function, that is, assuming nontrivial perturbation and/or truncation to compact support[-B,B]. Such response functions have been observed in measurement data generated by nearly all the analytical platforms currently used to assess DNA methylation, including the Sequenom EpiTyper, Infinium Mass Array and Restricted Representation Bisulphite Sequencing platforms [16].In the absence of perturbation or truncation of the response function, the results in this paper correspond to the theory of LAE (or median) regression as found in [2–4]. That is, if g(z)=1 and B=∞, then our MLE method corresponds with LAE regression. In this case, Q=2, ν=p2, ζ=-p2, f(z)=(p/2)e-p|z| (in one dimension), and so
(142)νζ2=p2p4=1p2=1(2f(0))2,
which is the asymptotic variance of the ordinary sample median for f [2].We present an original and practical method of obtaining the covariance matrix for the model coefficients. This involves evaluating(Xn,mTXn,m)-1, where although n might be large, m generally is not, and using generalized functions to numerically evaluate two one-dimensional integrals (to find ν and ζ). The calculation of ν and ζ takes into account the characteristics of the response function (truncation and perturbation) which would be ignored if we used median regression. This is possible when we know the response function parameters, or have fairly accurate estimates, as in our epigenetic application, modelling DNA methylation proportion deviations.For LAE regression, other methods of determining approximations to this covariance matrix may be found in the literature. In particular, in the method of quantile regression [5] implemented in the statistical package R, the covariance matrix for the MLE is calculated by resampling techniques, by bootstrapping or by using hierarchical spline models [17].We prove that, even for truncated and perturbed Laplace response functions, subject to certain restrictions, the maximum of the log-likelihood function occurs at a data point. This result is well-known in the case of LAE regression. A proof that the LAE estimator passes through at leastrX data points may be found in references [3, 4].Three asymptotically equivalent test statistics for LAE regression may be found in [18], namely, a likelihood ratio test statistic, a Wald test statistic, and a Lagrange multiplier test statistic. Our likelihood ratio test statistic is an original modification of the former, applicable to our general case (not restricted to LAE regression), calculated using generalized functions. An F test statistic for LAE regression is found in [4].When working with a model for which the response function is assumed to be a truncated Laplace probability density function, we could ignore the truncation to[-B,B]. However, taking the truncation into account reduces the variance in the model coefficient estimates by a factor (1-e-pB)2 and increases the log-likelihood ratio statistic by a factor (1-e-pB)-1. This effect is small if e-pB is small but becomes more significant as pB decreases, that is, as more of the density function is truncated. Refer to Section 4.9, (111) and Section 4.11, (123), and (124) which show that, for example, when g=g1,
(143)Dgen=2(1-e-pB)log(λ).
Hence, by taking into account the truncation, we can be more confident of our coefficient estimates and the value of appropriate beta coefficients than the standard formulae for LAE regression indicate. This effect will also be seen for small perturbations of the density function.


## 7.2. Original Formulae Enable Model Comparison

The original formulae derived in Section4 for the covariance matrix of the model coefficients and the generalized log-likelihood ratio statistic enable us to do a generalized analysis of variance, in our case comparing the means of two different groups of data. We see that in Section 6.3, the primiparous and multiparous means could be judged distinct with very high probability using either the Laplace or an amended Laplace response function. The P-values are slightly lower when using the Hermite modification. In the simulation in Section 6.2, using the Hermite modification decreases the P-value by an order of magnitude. In a situation in which p is low or two means are less distinct, the Hermite modification may prove very useful.Preliminary results indicate the use of an amended Laplace distribution enables a clearer separation of means than that given by other more standard procedures, for example, beta regression, in cases where independent evidence suggests that the means are different [16].


## 7.3. Summary

The Laplace distribution is the basis of many mathematical models (see [19]). Our focus has been modelling the distributions of errors in the proportions of DNA methylation measured at genomic CpG sites.Molecular biology deals with complex interactions both in terms of the physiology of the processes of interest and in the instrumentation required to measure these effects. The non-linearity of these processes can result in frequency distributions that are far from normal, so that application of  “standard” methods of statistical inference based on least squares may be inadequate. Methods which deal with the form of the frequency distribution directly such as maximum likelihood are necessary for adequate inference to be made.The Laplace or double exponential distribution considered here has been observed in molecular biology studies, where a significant proportion of high deviations appears to occur regularly [16, 20]. The extension by Hermite polynomials considered here provides flexibility for describing the tails of the distribution. However, as noted, the use of the Laplace distribution as the “key” function introduces problems in finding maximum likelihood estimators and particularly their standard errors. This paper presents both a practical method for dealing with these problems and the underlying asymptotic theory.



---

*Source: 496180-2013-11-02.xml*