A data augmentation approach for a class of statistical inference problems.

We present an algorithm for a class of statistical inference problems. The main idea is to reformulate the inference problem as an optimization procedure, based on the generation of surrogate (auxiliary) functions. This approach is motivated by the MM algorithm, combined with the systematic and iterative structure of the Expectation-Maximization algorithm. The resulting algorithm can deal with hidden variables in Maximum Likelihood and Maximum a Posteriori estimation problems, Instrumental Variables, Regularized Optimization and Constrained Optimization problems. The advantage of the proposed algorithm is to provide a systematic procedure to build surrogate functions for a class of problems where hidden variables are usually involved. Numerical examples show the benefits of the proposed approach.

1 Introduction

Problems in statistics and system identification often involve variables for which measurements are not available. Among others, real-life examples can be found in communication systems [1, 2] and systems with quantized data [3, 4]. In Maximum Likelihood (ML) estimation problems, the likelihood function is in general difficult to optimize by using closed-form expressions, and numerical approximations are usually cumbersome. These difficulties are traditionally avoided by the utilization of the Expectation-Maximization (EM) algorithm [5], where a surrogate (auxiliary) function is optimized instead of the main objective function. This surrogate function includes the complete data, i.e. the measurements and the random variables for which there are no measurements. The incorporation of such hidden data or latent variables is usually termed as data augmentation, where the main goal is to obtain, in general, simple and fast algorithms [6].

On the other hand, the MM (MM stands for Maximization-Minorization or Minimization-Majorization, depending on the optimization problem that needs to be solved) algorithm [7] is generally employed for solving more general optimization problems, not only for ML and Maximum a Posteriori (MAP) estimation problems. In general, the main motivation for using the MM algorithm is the lack of closed-form expressions for the solution of the optimization problem or dealing with objective cost functions that are not convex. Applications where the MM algorithm has been utilized include communication systems problems [8] and image processing [9]. For constrained optimization problems, an elegant solution is presented by Marks and Wright [10], where the constraints are incorporated via the formulation of surrogate functions. Surprisingly, Marks’ approach has not received the same attention from the scientific community when it comes to compare it with the EM and the MM algorithms. In fact, these three approaches are contemporary, but the EM algorithm has attracted most of the attention (out of the three methods), and it has been used for solving linear and nonlinear statistical inference problems in biology and engineering, see e.g. [11-16], amongst others. On the other hand, as shown in [7], the MM algorithm has obtained much less attention, while Marks’ approach is mostly known to a limited audience in the the Communication Systems community. These three approaches have important similarities: i) a surrogate function is defined and optimized in place of the original optimization problem, and ii) the solution is obtained iteratively. In general, these algorithms are “principles and recipes” [17] or a “philosophy” [7] for constructing solutions to a broad variety of optimization problems.

In this paper we adopt the ideas behind [5, 7, 10] to develop an algorithm for a special class of functions. Our approach generalizes the ones of [5, 7, 10] by reinterpreting the E-step in the EM algorithm and expressing the cost function in terms of an infinite mixture or kernel. This kind of problems can be interpreted as inverse problems that, in turn, can be posed as integral equations, such as channel modelling in wireless communications [18], estimation of the distribution of stellar rotational velocities [19], and mass estimation in particle physics problems [20, 21]. In particular cases, the kernel corresponds to a variance-mean Gaussian mixture (VMGM), see e.g. [22]. VMGMs, also referred to as normal variance-mean mixtures [23] and normal scale mixtures [24], have been considered in the literature for formulating EM-based approaches to solve ML [25] and MAP problems, including regularized sparse estimation problems [22, 26, 27]. Sparse estimation problems have been widely studied in the last two decades and several techniques have been developed that include different types of penalties or contraint, see e.g [28, 29] and different strategies for the formulation of those penalties/constraints [30]. In this paper we show that our proposal can also be considered for sparsity problems, however the analysis of the solution is out of the scope of the paper. Our approach is applicable to a wide class of functions, which allows for defining the likelihood function, the prior density function, and constraints as kernels, extending also the work in [10]. Thus, our work encompasses the following contributions: i) a systematic approach to constructing surrogate functions for a class of cost functions and constraints, ii) a class of kernels where the unknown quantities of the algorithm can be easily computed, and iii) a generalization of [5, 7, 10] by including the cost function and the constraints in one general expression. Our proposal is based, among other things, on a particular way to apply Jensen’s inequality [31]. In addition, we provide the details on how to construct quadratic surrogate functions for cost functions and constraints.

Our algorithm is tested by two examples. In the first one we considered the problem of estimating the rotational velocities of stars. The system model corresponds to the convolution of two probability density functions (pdf’s) and thus it is an infinite mixture. We show that our reinterpretation of the EM algorithm allows for the direct application of our proposal for the correct estimation of the parameter of interest. In the second example, we considered the estimation of the channel in wireless communications, where the true distribution can be either Rayleigh or Rice, depending on environment where the electromagnetic waves propagate. The problem is solved considering a sum of a Rayleigh and a Rice term, allowing for a more complex channel distribution. To select the more representative distribution, Akaike’s Information Criterion [32] was also considered in order to obtain the least complex model that exhibits the best possible fitting.

2 Rudiments of the proposed approach

2.1 The EM algorithm

Let us consider an estimation problem and its corresponding log-likelihood function defined as ℓ() = log p(|), where p(|) is the likelihood function, , and . Denoting the complete data by ∈ Ω(), and using Bayes’ theorem, we can obtain: Let us assume that at the ith iteration we have the estimate . By integrating at both sides of (1) with respect to we obtain , where

Using Jensen’s inequality [31], it is possible to show that for any value of , the function is decreasing. Hence, the optimization is only carried out on the auxiliary function because, by maximizing , the new parameter is such that the likelihood function increases (see e.g. [5, 33]).

In general, the EM method can be summarised as follows:

E-step: Compute the expected value of the joint likelihood function for the complete data (measurements and hidden variables) based on a given parameter estimate, . Thus, we have (see e.g. [5]):

M-step: Maximize the function (4), with respect to :

This succession of estimates converges to a stationary point of the log-likelihood function [34].

2.2 The MM algorithm

The idea behind the MM algorithm [7] is to construct a surrogate function , that majorizes (for minimization problems) or minorizes (for maximization problems) a given cost functions f() [7] at such that, where is an estimate of . Then, the surrogate function is iteratively optimized until convergence. Hence, for maximizing f() we have [35] For the construction of the surrogate function, popular techniques include the second order Taylor approximation, the quadratic upper bound principle and Jensen’s inequality for convex functions, see, e.g., [35].

Remark 1. The iterative strategy utilized in the MM algorithm converges to a local optimum since

2.3 Data augmentation in inference problems

Data augmentation algorithms are based on the construction of the augmented data and its many-to-one mapping Ω(y). This augmented data is assumed to describe a model from which the observed data y is obtained via marginalization [36]. That is, a system with a likelihood function p(y|) can be understood to arise from where the augmented data corresponds to (y, x) and x is the latent data [6, 36]. This idea has been utilized for supervised learning [37] and the development of the Bayesian Lasso [38], to mention a few examples. In those problems, the Laplace distribution is expressed as a two-level hierarchical-Bayes model. This equivalence is obtained from the representation of the Laplace distribution as a VMGM: In fact, there are several pdf’s than can be expressed as VMGMs, as shown in Table 1 [22], where g(θ) is the penalty term that can be expressed as a pdf. In addition, in [26] it was developed an early version of the methodology presented in this paper, exploring the estimation of a sparse parameter vector utilizing the ℓ-(pseudo)norm, with 0 < q < 1.

Table 1

Selection of mean-variance mixture representations for penalty functions.

.

Penalty function	g(θ)	u	μ	p(λ)
Ridge	(θ/τ)2	0	0	λ = 1
Lasso	|θ/τ|	0	0	Exponential
Bridge	|θ/τ|α	0	0	Stable
Generalized Double-Pareto	[(1+α)τ]log(1+|θ|(ατ))	0	0	Exp-Gamma

.

3 A systematic approach to construct surrogate functions for a class of inference problems

Here, we consider a general optimization cost defined as: where is a parameter vector, is a given data (i.e. measurements), is the complete data (comprised of the observed data and the hidden variables (unobserved data), Ω() is a mapping from the sample space of to the sample space of , K(⋅, ⋅) is a (positive) kernel function, and μ(⋅) is a measure, see e.g [31]. The definition in (9) is based on the definition of the auxiliary function in the EM algorithm [5], where it is assumed throughout the paper that there is a mapping that relates the not observed data to the observed data, and that the complete data lies in Ω() [5]. Notice that in (9) the kernel function may not be a pdf. However, several functions can be expressed in terms of a pdf. The most common cases are Gaussian kernels (yielding VMGMs) [23] and Laplace kernels (yielding Laplace mixtures) [39].

Remark 2. Notice that, as explained in Section 1, once the hidden data has been selected, the data augmentation procedure comes with the definition of in (9). From here, we follow the systematic nature of the EM and MM algorithms in terms of the iterative nature of the technique.

3.1 Constructing the surrogate function

Since we are considering the optimization of the function , we can also consider the optimization of the function Without modifying the cost function in (10), we can multiply and divide by the logarithm of the kernel function, obtaining: Let us assume that at the ith iteration we have the estimate . Then, we can multiply by and integrate on both sides of (11) with respect to dμ(), obtaining: where: are auxiliary functions. As in the EM algorithm, for any , and using Jensen’s inequality [31], we have: Hence, for any value of , the function in (14) is a decreasing function.

Remark 3. The kernel function K(, ) satisfies the standing assumption K(, ) > 0 since the proposed scheme is built, among others, on the logarithm of the kernel function K(, ). The definition of the kernel function in (9) allows for kernels that are not pdf’s. On the other hand, some kernels may correspond to a scaled version of a pdf. In that sense, for the cost function in (9) we can define a new kernel and a new measure as

Remark 4. In the proposed methodology, it is possible to optimize the surrogate function defined by since in (13) does not depend on the parameter . Thus, the proposed method corresponds to a variation of the EM algorithm that is not limited to probability density functions (e.g. the likelihood function) for solving ML and MAP estimation problems. Instead, our version considers general measures (μ()), where the mapping over the measurement data Ω() is a given set.

The idea behind using a surrogate function is to obtain a simpler algorithm for the optimization of the objective function when compared to the original optimization problem. This can be achieved iteratively if the Fisher Identity for the surrogate function and the objective function is satisfied. That is,

Lemma 1. For the class of objective functions in (9), the surrogate function in (13) satisfies the Fisher identity defined in (17).

Proof. From (12) we have: Next, let us consider the gradient of the auxiliary function : Hence, (17) holds.

Remark 5. Note that the Fisher identity in Lemma 1 is well known in the EM-framework. However, we have specialized this result for the problem in this paper (i.e. when K(, ) is not necessarily a probability density function)

Lemma 2. The surrogate function in (13) can be utilized to obtain an adequate surrogate function that satisfies the properties in Mark’s approach in (38)–(40).

Proof. Notice that in Mark’s approach the optimization problem corresponds to the minimization of the objective function. Hence, to maximize, we have . From (12) we can construct the surrogate functions and since The function satisfies ≥ 0, which implies that , satisfying (38). From we can obtain , satisfying (39). Finally, given that the auxiliary function satisfies (17), satisfies (40).

Remark 6. Since , it is simpler to consider the function instead of in penalized (regularized) and MAP estimation problems, as shown in [26] and [27].

We summarize our proposed algorithm in Table 2.

Table 2

Proposed algorithm.

Step 1: Find a kernel that satisfies (9).

Step 2: i = 0.

Step 3: Obtain an initial guess θ^^(i).

Step 4: Compute Q(θ,θ^^(i)) as in (13).

Step 5: Compute Q˜(θ,θ^^(i)).

Step 6: Incorporate Q˜(θ,θ^^(i)) in the optimization problem and solve.

Step 7: i = i + 1 and back to Step 4 until convergence.

3.2 Surrogate functions for inverse problems

The objective function in (9) can be understood as an integral equation [40], where the unknown function of the integral equation corresponds to the kernel K(, ). In this kind of problems, samples from are available, whilst dμ() is assumed known. Several problems can be posed as this kind of problems. In the following, we explain how to use the approach presented in this paper to solve different inverse problems that arise in the integral equation form.

3.2.1 Stellar rotational velocity estimation

One of the many problems in Astronomy deals with is the estimation of rotational velocities of stars. This particular problem is of great importance, since it allows astronomers to describe and model the stars formation, their internal structure and evolution, as well as how they interact with other stars, see e.g. [19, 41, 42].

Modern telescopes allow for the measurement of the rotational velocities from the telescope point of view, that is, a projection of the true rotational velocity. This is modelled (spatially) as the convolution of the true rotational velocity pdf and a uniform distribution over the sphere (for more details see e.g. [19]): where p(y|σ) is the uniform projected rotational velocity pdf and p(x|σ) is the true rotational velocity pdf to be estimated, and σ a hyperparameter. Thus, we can define the kernel function as K(x, σ) = p(x|σ) and . This definition allows for the direct utilization of the expressions in (13) in order to estimate the parameter σ that defines the unknown rotational velocity pdf. We illustrate with one example in Section 6.1.

3.2.2 Channel estimation in wireless communications

It is well known in the Communications community that the wireless channel corresponds to the superposition of different copies of the transmitted signal that have been reflected, refracted and scattered. Thus, those copies arrive at the receiver with different phase. Those components are referred to as multipath components [43]. On the other hand, it has been shown empirically that a good model for the multipath channel corresponds to either Rayleigh or Rice distributions [44]. However, there are cases when those distributions do not provide a good model for the channel. One of those cases corresponds to the presense of different channel models in the vicinity of the transmitter and the vicinity of the receiver. This situation ocurrs particularly in the so called urban scenarios, where the channel exhibits different behaviouirs in different places. For this scenario, an adequate model that takes into consideration different models and a transition from a local distribution and a global distribution is a continuous mixture of the form [18] where K(, ) is a pdf in a local area, μ() is a distribution of which depends on the path from transmitter to the local cluster, and is a vector in the parameter space. Our approach can be directly used in order to estimate the true nature of the channel when expressed as a mixture. This can be done since the four most common chanel distributions are Rayleigh, Rice, log-normal and Nakagami-m [45]. Hence, assuming that the local and global distribution families are known, the attainment of the auxiliary function is straightforward and thus the ML estimate of .

3.2.3 Neutrino mass search in particle physics

In the Particle Physics community there is a plethora of works dealing with the estimation of masses of neutrinos, see e.g. [20] and the references therein. Among the methods that are generally used to determine the absolute masses of neutrinos we find the β-decay and the direct determination of neutrino mass, see e.g [21]. In β-decay methods, the neutrinos are analized based on their energy spectrums, where the measurements, corresponding to the observed β-spectra, are associated with an integral equation of the form [20, 21] where b is a constant that represents the measurement noise, R(E) is the emitted β-spectra, and T′(E, U) is the impulse response of the equipment. Again, by regarding T′(E, U) as the kernel function and R(E)dE as a measure, the attainment of the auxiliary function is straightforward.

3.2.4 Estimation of mixture distributions

Mixture distributions have been widely studied in the literature, particularly finite mixtures, see e.g. [46, 47]. Their representation can be expressed in a general fashion by the notation [48] where K(, ) is a suitable function that may be either a pdf (for continuous random variables) or a probability function (for discrete random variables). The expression in (22) represents both the sum for a finite mixture, or the integral for an infinite mixture.

Our approach can be tailored for the estimation of the parameters of a finite mixture of the form where we have assumed that , , and λ = [λ1,…, λ] are the mixing weights, ϕ(, θ) are the components densities parametrised by θ. The kernel function defined in our approach can be utilized to represent jth component in the discrete mixture in (25) as where the = [1,…, ] is the vector of parameters to be estimated, with = [λ, ]. Notice that the dependence with respect to the variable is implicit. Utilizing the expression derived in (13) we obtain the following E-step

Notice that, as shown in [49], the auxiliary function in (27) is the same that we obtain when we consider that the data are fully categorized, i.e. each y, k = 1,…, N, is assumed to be drawn from only one distribution of the mixture. This assumption yields a data augmentation problem that is solved using the EM algorithm [46].

Remark 7. If we consider a combination of infinite mixtures and finite mixtures of the form with we can utlize the same approach described here to solve the problem of estimating the parameters in (28), . In this case, the jth kernel is defined as [50] and measure Then, the log-likelihood function can be expressed as with This choice of functions leads to the direct implementation of our proposal, from which it is obtained and the ML estimator can be locally obtained from

4 Marks’ approach for constrained optimization

4.1 Constrained problems in statistical inference

Statistical Inference and System Identification techniques include a variety of methods that can be used in order to obtain a model of a system from data. Classical methods, such as Least Squares, ML, MAP [51], Prediction Error Method, Instrumental Variables [52], and Stochastic Embedding [53] have been considered in the literature for such task. However, the increasing complexity of modern system models has motivated researchers to revisit and reconsider those techniques for some problems. This has resulted in the incorporation of constraints and penalties, yielding an often more complicated optimization problem. For instance, it has been shown that the incorporation of linear equality constraints may improve the accuracy of the parameter estimates, see e.g [54]. On the other hand, the incorporation of regularization terms (or penalties) also improves the accuracy of the estimates, reducing the effect of noise and eliminating spurious local minima [55]. Regularization can be mainly incorporated in two ways: by adding regularizing constraints (a penalty function) or by including a probability density function (pdf) as a prior distribution for the parameters, see e.g. [27]. Another way to improve the estimation is by incorporating inequality constraints, where certain functions of the parameters may be required, for physical reasons amongst others, to lie between certain bounds [56]. From this point of view, it is possible to consider the classical methods with constraints or penalties, as in [53, 55–57].

Perhaps one of the most utilized approaches for penalized estimation (with complicated non-linear expressions) is the MM algorithm—for details on the MM algorithm see Section 2.2. This technique allows for the utilization of a surrogate function that is simple to handle, in terms of derivatives and optimization techniques, and that is, in turn, iteratively solved. However, its inequality constraint counterpart, here referred to as Marks’ approach [10], is not as well known as the MM algorithm. Moreover, there is no straightforward manner to obtain such surrogate function. In this paper we focus on a systematic way to obtain the corresponding surrogate function using Marks’ approach for a class of constraints.

4.2 Mark’s approach

The approach in [10] deals with inequality constraints by using a similar approach to the EM and MM algorithms. The basic idea is, again, to generate a surrogate function that allows for an iterative procedure whose optimum value is the optimum value of the original optimization problem.

Let us consider the following constrained optimization problem: where f() is the objective function and g() encodes the constraint of the optimization problem. In particular, let us focus on the case where g() is not a convex function. This implies that the optimization problem cannot be solved directly using standard techniques, such as quadratic programming or fractional programming. This difficulty can be overcome by utilizing a surrogate function at a given estimate , such that

Provided the above properties are satisfied, then the following approximation of (37): iteratively converges to the solution of the optimization problem (37). As shown in [10], the optimization problem in (41) is equivalent to the original problem (37), since the solution of (41) converges to a point that satisfies the Karush-Kuhn-Tucker conditions of the original optimization problem.

Remark 8. Mark’s approach can be considered as a generalization of the MM algorithm, since the latter can be derived (for a broad class of problems) from the former. Let us consider the following problem: Using the epigraph representation of (42) [58], we obtain the equivalent problem Using Mark’s approach (41), we can iteratively find a local optimum of (42) via where in (44) is a surrogate function for f() in (42). From the epigraph representation we then obtain which is the definition of the MM algorithm (see 2.2) for more details.

5 A quadratic surrogate function for a class of kernels

In this section we focus on a special class of the kernel functions K(, ). For this particular class, the following is satisfied: where A() is a matrix and b is a vector, both of adequate dimensions. Then, we have that

Remark 9. Notice that the previous expression is linear with respect to . This implies that the function is quadratic with respect to the parameter vector .

From the Fisher Identity in (17) we have that from which we can solve for R in some cases. In other cases, the matrix R can also be computed using Monte Carlo algorithms. In particular, if A() is a diagonal matrix, then R is also a diagonal matrix defined by R = diag[r1, r2,…]. Thus, we have where θ is the ith component of the parameter vector , is the ith component of the estimate , r is the ith element of the diagonal of R, and b is the ith element of the vector b. Hence, when optimizing the auxiliary function we obtain Equivalently, This implies that in our approach, it is not necessary to obtain the auxiliary function and optimize it. Instead, by computing R and b at every iteration, the new estimate can be obtained.

Remark 10. The computation of the matrix can be cumbersome when the matrix () is not diagonal. In those cases, the integral that defines can be computed utilizing Markov Chain Monte Carlo, quasi-Monte Carlo [59] or quadrature methods [60].

The class defined in (46) arises naturally when dealing with VMGM, because the corresponding kernel function is normal and, thus, its logarithm is a quadratic function.

In particular, the utilization of VMGM encompasses different expressions commonly used for parameter estimation. Indeed, we have:

(i) Lasso: The Lasso [61], expressed as a Laplace pdf, is represented by a VMGM [38], with p(|) a zero-mean Gaussian distribution (of iid terms) as the kernel and p() an exponential distribution with parameter γ2/2. That is,

(ii) Elastic-Net: The Elastic-Net penalty [62] is interpreted as a pdf if it corresponds to the product of two pdf’s, a Laplacian (as in the Lasso case) and a Gaussian pdf. In this sense, by grouping those pdf’s, we obtain [63] with .

(iii) Group-Lasso: The Group-Lasso penalty is obtained via VMGM-representation as where is the chi-squared distribution with l degrees of freedom.

For the class of kernels here described, the proposed method for constructing surrogate functions can also be understood as part of sequential quadratic programming (SQP) methods [64] when, for instance, the above penalties are utilized as constraints in a constrained ML estimation problem. Indeed, the general case of equality and inequality-constrained minimization problems is defined as [65]: which is solved by iteratively defining quadratic functions that approximate the objective function and the inequality constraint around a current iterate . In the same way, our proposal generates an algorithm with quadratic surrogate functions, where an auxiliary function in (18) must be constructed for f() and/or g() in (53).

6 Numerical examples

In this section, we illustrate our proposed algorithm with two numerical examples.

6.1 Example 1: Estimation of the distribution of stellar rotational velocities

A common model for p(x|σ) found in the Astronomy literature is the Maxwellian distribution (see e.g. [19, 66]) In practice, the measurements correspond to realizations of p(y|σ) [19], from which the likelihood function can be defined as: where y = [y1,…, y], x is Maxwellian distributed, and N is the number of measurement points. Hence, the log-likelihood function can be expressed as: If we define the complete data = (x, y), the kernel function K(⋅, ⋅) and the measure μ(⋅) in (9) can be defined as and Then, the log-likelihood function in (56) can be written as: with Thus, the ML estimator is obtained from: Since the parameter that is needed to be estimated is part of the convolution in (19), the optimization problem in (61) cannot be solved in a straightforward manner. Instead, we utilize the re-interpretation of the EM algorithm that we propose for solving (61).

First, notice that from the surrogate function can be expressed as: with For convenience, we can differentiate the auxiliary function in (62) with respect to 1/σ obtaining: Then, equating to zero and solving for σ we finally obtain where

In Table 3 we summarized the specialisation of our proposed algorithm for this example.

Table 3

Proposed algorithm for Maxwellian distribution estimation in Example 1.

Step 1: i = 0.

Step 2: Obtain an initial guess σ^^(i).

Step 3: Compute the integral given by (66).

Step 4: Compute σ^^(i+1) (65)

Step 5: i = i + 1 and back to Step 3 until convergence.

For the numerical simulation, we have considered the problem solved in [19], with the true dispersion parameter σ0 = 8. The measurement data y = [y1,…, y] was generated using the Slice Sampler (see e.g. [67]) applied to (19). The simulation setup is as follows:

The data length is given by N = 10000.

The number of Monte Carlo (MC) simulations is 50.

The stopping criterion is given by: or the maximum number of iterations of 100 has been reached.

The results are shown in Fig 1, were the estimated p(x|σ) for each MC simulation is shown. It is clear that the estimated Maxwellian distributions are very similar to the true density distribution. The mean value of the estimated parameter was . The estimation from each MC simulation is shown in Fig 2. It can be clearly seen that the estimated parameter is close to the true value.

Fig 1

Estimated distribution for the stellar rotational velocity.

Fig 2

Convergence of the proposed approach to the global optimum.

6.2 Example 2: Channel estimation in wireless communications

When modelling the wireless channel, a popular technique that is commonly used corresponds to the transmission of a sine tone at a given frequency, see e.g. [68]. The received power is then modelled as a random variable. The corresponding distribution has been widely studied in the literature from measurements, and the empirical data have shown that the two most common distributions are Rayleigh an Rice [44, 45]. Hence, using similar ideas as in [50], in this example we formulate the channel distribution as a discrete sum based on a Rayleigh and a Rice component to determine the nature of the wireless channel.

First, the discrete mixture that we want to adjust from data is given by where and I0(⋅) is the modified Bessel function of zeroth order. In addition, we must also include the constraint λ1+ λ2 = 1 so p(x) is a pdf. Thus, we can directly apply our proposed approach by assuming that each measurement point can be associated with a hidden variable that describes if such data point comes from the Rayleigh component or the Rice component, as it is traditionally formulated when dealing with discrete mixtures [46]. Hence, the auxiliary function is given by where is the current estimate, ζ is the unobserved (hidden) data and ζ is an indicator parameter such that ζ = 1 if the t-th observation comes from component j and is zero otherwise. It is given by: The estimate at the next iteration is then given by: from which we obtain the following expressions: with We also consider the utilization of Akaike’s Information Criterion (AIC) in order to obtain an accurate yet simple model and, thus, discriminating from a Rayleigh channel, a Rice channel, and a mixture of both.

With the above formulation, we consider two cases: a Rayleigh distributed channel and a Rice distributed channel.

6.2.1 Rayleigh distributed channel

In this example, the random variable x is drawn from the Rayleigh distribution with σ2 = 1, using the Slice Sampler [69]. The best estimated model corresponds to the single Rayleigh component in the mixture. The corresponding estimation of σ2 yields . Fig 3 shows the true Rayleigh distribution and the mean estimated pdf from the 50 MC simulations. We observe an important agreement between the true pdf and the estimated model.

Fig 3

Rayleigh distribution estimation using a Rayleigh-Rice mixture.

6.2.2 Rice distributed channel

In this case, the data is drawn from the Rice distribution with v = 4 and σ2 = 1, using the Slice Sampler. The best model is selected as a single Rician component. The corresponding estimated parameters are and . In Fig 4 we show the true Rice distribution and the mean estimated pdf from 50 MC simulations. We can observe that the estimator exhibits a good performance for the estimation of a Rice distribution.

Fig 4

Rice distribution estimation using a Rayleigh-Rice mixture.

Conclusions and future work

In this paper we have presented a systematic approach for constructing surrogate functions in a wide range of inference problems. Our approach can be utilized for constructing surrogate functions for both the cost function and the constraints, generalizing the popular EM and MM algorithms. Our approach is based on the utilization of data augmentation and kernel functions, yielding simple optimization algorithms when the kernel can be expressed as VMGM. We have shown how our proposal can be utilized to solve inverse problems that are expressed as integral equations and mixture distributions.

In addition, we have shown that our approach can be utilised for constrained/penalized ML and MAP estimations problems. In particular, common problems in statistical inference can directly be solved using our proposal since they can be posed as Variance Mean Gaussian Mixtures (VMGM), yielding quadratic surrogate functions.

In the last two decades the problem of sparse estimation has attracted a lot of attention. Since our approach can be utilized in those problems, and since it is based on the principles of the MM algorithm, a detailed analysis can be done in terms of accuracy and convergence of our technique, and compared against other techniques, such as the ones in [28, 30], and [29], where different Lasso-type problems are compared, the MM algorithm is utilized in constrained problems for ML estimation in generalized linear model regression, and the MM algorithm is used for (unconstrained) sparse estimation under non-convex penalties, respectively.

Monte Carlo simulations for Examples 1 and 2.

(ZIP)

Click here for additional data file.