Regularization Methods Based on the L_q-Likelihood for Linear Models with Heavy-Tailed Errors

Abstract

We propose regularization methods for linear models based on the $L_q$-likelihood, which is a generalization of the log-likelihood using a power function. Regularization methods are popular for the estimation in the normal linear model. However, heavy-tailed errors are also important in statistics and machine learning. We assume q-normal distributions as the errors in linear models. A q-normal distribution is heavy-tailed, which is defined using a power function, not the exponential function. We find that the proposed methods for linear models with q-normal errors coincide with the ordinary regularization methods that are applied to the normal linear model. The proposed methods can be computed using existing packages because they are penalized least squares methods. We examine the proposed methods using numerical experiments, showing that the methods perform well, even when the error is heavy-tailed. The numerical experiments also illustrate that our methods work well in model selection and generalization, especially when the error is slightly heavy-tailed.

1. Introduction

We propose regularization methods based on the $L_q$-likelihood for linear models with heavy-tailed errors. These methods turn out to coincide with the ordinary regularization methods that are used for the normal linear model. The proposed methods work efficiently, and can be computed using existing packages.

Linear models are widely applied, and many methods have been proposed for estimation, prediction, and other purposes. For example, for estimation and variable selection in the normal linear model, the literature on sparse estimation includes the least absolute shrinkage and selection operator (LASSO) [1], smoothly clipped absolute deviation (SCAD) [2], Dantzig selector [3], and minimax concave penalty (MCP) [4]. The LASSO has been studied extensively and generalized to many models, including the generalized linear models [5]. As is well known, the regularization methods have many good properties. Many regularization methods are the penalized maximum likelihood estimators, that is, minimizing the sum of the negative log-likelihood and a penalty. The literature proposed various penalties. As described later, our regularization methods use another likelihood with existing penalties.

Because the regularization methods for the normal linear model are useful, they are sometimes used in linear models with non-normal errors. Here, popular errors include the Cauchy error and the t-distribution error, both of which are heavy-tailed errors. For example, Ref. [6] partly consider the Cauchy and t-distribution errors in their extensive experiments. These heavy-tailed distributions are known to be q-normal distributions, which are studied in the literature on statistical mechanics [7,8,9]. The q-normal model is also studied in the literature on the generalized Cauchy distribution. For example, see [10,11,12,13].

In this study, we consider the problem of a linear regression with a q-normal error. We propose sparse estimation methods based on the $L_q$-likelihood, which is a generalization of the log-likelihood using a power function. The maximizer of the $L_q$-likelihood, the maximum $L_q$-likelihood estimator (MLqE), is investigated by [14] as an extension of the ordinary maximum likelihood estimator (MLE). Ref. [14] studies the asymptotic properties of the MLqE. However, we are interested in the regularization, not in the MLqE, because regularization estimators can be better than the MLqE. We examine the proposed methods using numerical experiments. The experiments show that our methods perform well in model selection and generalization, even when the error is heavy-tailed. Moreover, we consider the effects of the sample size, dimension and sparseness of the parameter, and value of the nonzero elements in the numerical experiments.

We also find that the proposed methods for linear models with q-normal errors coincide with the ordinary regularization methods that are applied to the normal linear model. This finding partly justifies the use of the ordinary regularization methods for linear regressions with heavy-tailed errors. Moreover, the proposed methods are penalized least squares methods, and can be efficiently computed by existing packages.

The rest of the paper is organized as follows. In Section 2, we introduce several tools, including the normal linear model, regularization methods, $L_q$-likelihood, and q-normal models. In Section 3, we describe the problem under consideration, that is, estimations in linear models with q-normal errors. Moreover, we propose several regularization methods based on the $L_q$-likelihood. In Section 4, we evaluate the proposed methods using numerical experiments. Section 5 concludes the paper.

2. Preliminaries

2.1. Normal Linear Model and Sparse Estimation

First, we introduce the normal linear model, the estimation of which is a basic problem in statistics and machine learning [15]. Furthermore, we briefly describe some well-known regularization methods.

The normal linear model is defined as follows. A response is represented by a linear combination of explanatory variables $x_1,x_2,\cdots ,x_d$ as

$$y^a={\theta }^0+\sum _{i=1}^d{x}_i^a{\theta }^i+{\epsilon }^a(a=1,2,\cdots ,n),$$

(1)

where $y^a$ is the response of the a-th sample, n is the sample size, d is the number of explanatory variables, $x_i^a$ is the i-th explanatory variable of the a-th sample, ${\epsilon }^a$ is a normal error with mean zero and known variance, and the regression coefficient $\mathit{\theta }={({\theta }^0,{\theta }^1,\cdots ,{\theta }^d)}^{\top }$ is the parameter to be estimated. The normal linear model is equivalently given by

$$\mathit{\mu }=X\mathit{\theta },$$

where ${\mu }^a=\mathrm{E}\left[y^a\right]$ is the expectation of the response $y^a$, $\mathit{\mu }=\left({\mu }^a\right)$, and $X=\left(x_i^a\right)$ is a design matrix of size $n\times (d+1)$, with $x_0^a=1(a=1,2,\cdots ,n)$. Moreover, we define a row vector ${\mathit{x}}^a(a=1,2,\cdots ,n)$ as ${\mathit{x}}^a=(1,x_1^a,x_2^a,\cdots ,x_d^a)$, and a column vector ${\mathit{x}}_i(i=0,1,2,\cdots ,d)$ as ${\mathit{x}}_i={(x_i^1,x_i^2,\cdots ,x_i^n)}^{\top }$, which results in $X={({\mathit{x}}^{1\top },{\mathit{x}}^{2\top },\cdots ,{\mathit{x}}^{n\top })}^{\top }=({\mathit{x}}_0,{\mathit{x}}_1,{\mathit{x}}_2,\cdots ,{\mathit{x}}_d)$. Let $\mathit{y}=\left(y^a\right)$ be the response vector of length n. We assume that each column vector ${\mathit{x}}_i(i=1,2,\cdots ,d)$ is standardized, as follows: ${\sum }_{a=1}^n{x}_i^a=0$ and $\parallel {\mathit{x}}_i\parallel =1$, for $i=1,2,\cdots ,d$.

As is well known, some regularization methods for the normal linear model are formulated as an optimization problem in the form of

$$\underset{\mathit{\theta }\in {\mathbb{R}}^{d+1}}{min}\left\{\frac{1}{2n}{\parallel \mathit{y}-X\mathit{\theta }\parallel }^2+{\rho }_{\lambda }\left(\mathit{\theta }\right)\right\},$$

(2)

where ${\rho }_{\lambda }\left(\mathit{\theta }\right)$ is a penalty term, and $\lambda \ge 0$ is a regularization parameter. The LASSO [1] uses ${\rho }_{\lambda }\left(\mathit{\theta }\right)={\lambda \parallel \mathit{\theta }\parallel }_1=\lambda {\sum }_{i=1}^d\left|{\theta }^i\right|$. The path of the LASSO estimator when $\lambda$ varies can be made by the least angle regression (LARS) algorithm [16]. The SCAD [2] uses

$$\begin{array}{c}\hfill {\rho }_{\lambda }\left(\mathit{\theta }\right)=\left\{\begin{array}{c}{\displaystyle \sum _{i=1}^d\lambda |{\theta }^i\left|\right(|{\theta }^i|\le \lambda ),}\hfill \\ {\displaystyle -\sum _{i=1}^d\frac{|{\theta }^i{|}^2-2a\lambda \left|{\theta }^i\right|+{\lambda }^2}{2(a-1)}(\lambda <|{\theta }^i|\le a\lambda ),}\hfill \\ {\displaystyle \sum _{i=1}^d\frac{(a+1){\lambda }^2}{2}(a\lambda <|{\theta }^i\left|\right),}\hfill \end{array}\right.\end{array}$$

(3)

and the MCP [4] uses

$$\begin{array}{c}\hfill {\rho }_{\lambda }\left(\mathit{\theta }\right)=\lambda \sum _{i=1}^d{\int }_0^{|{\theta }^i|}{\left(1-\frac{u}{\gamma \lambda }\right)}_{+}\mathrm{d}u,\end{array}$$

(4)

where $a(>2)$ and $\gamma (>0)$ are tuning parameters.

The regularization problem given in (2) can be represented by

$$\underset{\mathit{\theta }\in {\mathbb{R}}^{d+1}}{min}\left\{-\frac{1}{n}logf\left(\mathit{y}\right|\mathit{\theta })+{\rho }_{\lambda }\left(\mathit{\theta }\right)\right\},$$

(5)

where $f\left(\mathit{y}\right|\mathit{\theta })$ is the probability density function of the statistical model. Note that $logf\left(\mathit{y}\right|\mathit{\theta })$ is the log-likelihood.

2.2. $L_q$-Likelihood

The $L_q$-likelihood is a generalization of the log-likelihood that uses a power function instead of the logarithmic function. Let $\mathit{y}={(y^1,y^2,\cdots ,y^n)}^{\top }$ be a vector of independent and identically distributed (i.i.d.) observations, and let $\mathit{\theta }$ be a parameter of a statistical model. For $q>0(q\ne 1)$, the $L_q$-likelihood function is defined as

$$\begin{array}{c}\hfill L_q\left(\mathit{\theta }\right|\mathit{y})=\sum _{a=1}^n{log}_{q}f\left(y^a\right|\mathit{\theta }),\end{array}$$

(6)

where $f(·|\mathit{\theta })$ is a probability density function of the statistical model, and

$$\begin{array}{c}\hfill {log}_q\left(u\right)=\frac{1}{1-q}(u^{1-q}-1)(u>0)\end{array}$$

is the q-logarithmic function [9]. For $q=1$, we define

$$\begin{array}{c}\hfill {log}_1\left(u\right)=logu(u>0),\end{array}$$

which is the ordinary logarithmic function. When $q=1$, the $L_q$-likelihood is the log-likelihood.

The MLqE is defined as the estimator that maximizes the $L_q$-likelihood. [14] studied the asymptotic performance of the MLqE, showing that it enjoys good asymptotic properties (e.g., asymptotic normality).

2.3. q-Normal Model

Before defining the q-normal distribution [7,8,9], we introduce the q-exponential function. For $q>0(q\ne 1)$, the q-exponential function is the inverse function of the q-logarithmic function, and is given by

$$\begin{array}{c}\hfill {exp}_q\left(u\right)={\{1+(1-q)u\}}^{\frac{1}{1-q}}(u<-1/(1-q)).\end{array}$$

For $q=1$, the 1-exponential function is the ordinary exponential function

$$\begin{array}{c}\hfill {exp}_1\left(u\right)=expu(u\in \mathbb{R}).\end{array}$$

Using the q-exponential function, the q-normal model is given by

$$\begin{array}{cc}\hfill {\mathcal{S}}_q& =\left\{f_q\left(y\right|\xi ,\sigma )\right|\xi \in \mathsf{\Xi },\sigma >0\},\hfill \\ \hfill f_q\left(y\right|\xi ,\sigma )& =\frac{1}{Z_q}{exp}_q\left\{-\frac{1}{3-q}{\left(\frac{y-\xi }{\sigma }\right)}^2\right\}\hfill \\ \hfill & =\frac{1}{Z_q}{\left\{1-\frac{1-q}{3-q}{\left(\frac{y-\xi }{\sigma }\right)}^2\right\}}^{\frac{1}{1-q}},\hfill \end{array}$$

where $\xi$ is a location parameter, $\mathsf{\Xi }\subset \mathbb{R}$ is the parameter space, and $\sigma$ is a dispersion parameter. The constant $Z_q$ is a normalizing constant.

We assume that $1\le q<3$, which ensures that the sample space is the real line itself, not just part of it. Moreover, the parameter space is $\mathsf{\Xi }=\mathbb{R}$ when $1\le q<3$.

For example, the 1-normal model is the ordinary normal model. Another example is the Cauchy distribution for $q=2$:

$$\begin{array}{c}\hfill f_2\left(y\right|\mu ,\sigma )=\frac{1}{\sigma B(\frac{1}{2},\frac{1}{2})}{\left(1+\frac{{(y-\xi )}^2}{{\sigma }^2}\right)}^{-1},\end{array}$$

where $B(·,·)$ is the beta function. Furthermore, the t-distribution of the degree of freedom $\nu$ is obtained for $q=1+2/(\nu +1)$:

$$\begin{array}{c}\hfill f_{1+2/(\nu +1)}\left(y\right|\mu ,\sigma )=\frac{1}{\sqrt{\nu }\sigma B(\frac{\nu }{2},\frac{1}{2})}{\left(1+\frac{{(y-\mu )}^2}{\nu {\sigma }^2}\right)}^{-\frac{\nu +1}{2}}.\end{array}$$

3. Problem and Estimation Method

3.1. Linear Model with q-Normal Error

In this subsection, we formulate our problem, that is, a linear regression with a heavy-tailed error. The errors of the Cauchy and t-distributions in linear models have been studied by researchers in the context of heavy-tailed errors [17,18,19,20]. However, they focused mainly on the least squares methods, whereas we are interested in sparse estimators. Moreover, our approach is based on the $L_q$-likelihood, not the ordinary log-likelihood.

We examine the problem of estimating the linear model given in (1) with i.i.d. errors from a q-normal distribution; henceforth, we refer to this as the q-normal linear model. In terms of probability distributions, we wish to estimate the parameter $\mathit{\theta }$ of the q-normal linear model ${\mathcal{M}}_q$:    

$$\begin{array}{cc}\hfill {\mathcal{M}}_q& =\left\{f(·|\mathit{\theta })\right|\mathit{\theta }\in {\mathbb{R}}^{d+1}\},\hfill \\ \hfill f\left(\mathit{y}\right|\mathit{\theta })& =\frac{1}{Z_q^n}\prod _{a=1}^n{exp}_q\left\{-\frac{{\left(y^a-{\mathit{x}}^a\mathit{\theta }\right)}^2}{3-q}\right\}\hfill \\ \hfill & =\frac{1}{Z_q^n}\prod _{a=1}^n{\left\{1-\frac{1-q}{3-q}{\left(y^a-{\mathit{x}}^a\mathit{\theta }\right)}^2\right\}}^{\frac{1}{1-q}},\hfill \end{array}$$

(7)

where the dispersion parameter is assumed to be known ($\sigma =1$). The 1-normal linear model is identical to the normal linear model, as described in Section 2.1.

3.2. $L_q$-Likelihood-Based Regularization Methods

We propose regularization methods based on the $L_q$-likelihood. For q-normal linear models, the proposed methods coincide with the original regularization methods for the normal linear model. In other words, we apply the ordinary regularization methods as if the error distribution were a normal distribution. The literature describes how to compute the proposed methods efficiently. Moreover, our method calculates the MLqE.

We define the $L_q$-likelihood for the q-normal linear model in (7) as (6), where $\mathit{\theta }$ is the regression coefficient. Note that the components of $\mathit{y}$ are not assumed to be identically distributed because their distributions are dependent on the explanatory variables.

The $L_q$-likelihood for the q-normal linear model is

$$\begin{array}{cc}\hfill L_q\left(\mathit{\theta }\right|\mathit{y})& =\sum _{a=1}^n{log}_{q}f\left(y^a\right|\mathit{\theta })\hfill \\ \hfill & =\sum _{a=1}^n{log}_q\left[\frac{1}{Z_q}{exp}_q\left\{-\frac{{\left(y^a-{\mathit{x}}^a\mathit{\theta }\right)}^2}{3-q}\right\}\right]\hfill \\ \hfill & =-\frac{Z_q^{q-1}}{3-q}{\parallel \mathit{y}-X\mathit{\theta }\parallel }^2-n{log}_q\left(Z_q\right),\hfill \end{array}$$

(8)

where the second term is a constant. The MLqE of the parameter $\mathit{\theta }$ is defined as the maximizer of the $L_q$-likelihood. In the q-normal linear model, the MLqE is equal to the ordinary least square, the MLE for the normal linear model.

We propose a LASSO, SCAD, and MCP based on the $L_q$-likelihood by replacing the log-likelihood with the $L_q$-likelihood in the optimization problem in (5). That is, the $L_q$-likelihood-based regularization methods are given in the form of

$$\underset{\mathit{\theta }\in {\mathbb{R}}^{d+1}}{min}\left\{-\frac{1}{n}{L}_q\left(\mathit{\theta }\right|\mathit{y})+{\rho }_{\lambda }\left(\mathit{\theta }\right)\right\}.$$

(9)

The penalty ${\rho }_{\lambda }$ is ${\lambda \parallel \mathit{\theta }\parallel }_1$ for the LASSO, (3) for the SCAD, and (4) for the MCP. Note that the estimator for $\lambda =0$ is the MLqE. As a special case, the proposed methods are the ordinary regularization methods when $q=1$.

Because of (8) and (9), for the q-normal linear models, the $L_q$-likelihood-based regularization methods are essentially the same as the penalized least square (2). In other words, we implicitly use the $L_q$-likelihood-based regularization methods when we apply the ordinary LASSO, SCAD, and MCP to data with heavy-tailed errors.

4. Numerical Experiments

In this section, we describe the results of our numerical experiments and compare the proposed methods. Here, we focus on model selection and generalization.

Our methods do not require additional implementations because the LASSO, SCAD, and MCP are already implemented in software packages. In the experiments, we use the ncvreg package of the software R.

4.1. Setting

The procedure for the experiments is as follows. We fix the value q of the q-normal linear model and the $L_q$-likelihood, the dimension d of the parameter $\mathit{\theta }$, the ratio of nonzero components $r_{\mathrm{nz}}$ of $\mathit{\theta }$, the true value ${\theta }_0$ of the nonzero components of $\mathit{\theta }$, and the sample size n. The value of q is selected from $1,13/11,3/2,5/3,2,2.01,2.1$, and $2.5$, where $q=13/11$ is the t-distribution with $\nu =10$ degrees of freedom, $q=3/2$ is the t-distribution with $\nu =3$ degrees of freedom, and $q=5/3$ is the t-distribution with $\nu =2$ degrees of freedom. The sample size is $n=100$ or $n=1000$. The true parameter consists of $d\times r_{\mathrm{nz}}$${\theta }_0$s and $d\times (1-r_{\mathrm{nz}})$ zeros. All cases are illustrated in

Table 1. All cases in the experiments. Each case is studied for the values of q and n.

${\mathit{\theta }}_0$	${\mathit{r}}_{\mathbf{nz}}=0.2$		${\mathit{r}}_{\mathbf{nz}}=0.4$		${\mathit{r}}_{\mathbf{nz}}=0.6$		${\mathit{r}}_{\mathbf{nz}}=0.8$
	$\mathit{d}=\mathbf{10}$	$\mathit{d}=\mathbf{100}$	$\mathit{d}=\mathbf{10}$	$\mathit{d}=\mathbf{100}$	$\mathit{d}=\mathbf{10}$	$\mathit{d}=\mathbf{100}$	$\mathit{d}=\mathbf{10}$	$\mathit{d}=\mathbf{100}$
${10}^0$	1	5	9	13	17	21	25	29
${10}^1$	2	6	10	14	18	22	26	30
${10}^2$	3	7	11	15	19	23	27	31
${10}^3$	4	8	12	16	20	24	28	32

.

For each of $m=1000$ trials, we create the design matrix X using the rnorm() function in R. The response $\mathit{y}$ is generated as q-normal random variables using the qGaussian package. For the estimation, we apply the ncvreg() function to $(\mathit{y},X)$ with the default options; for example, the values of the tuning parameters are $a=3.7$ and $\gamma =3$.

To select one model and one estimate from a sequence of parameter estimates generated by a method, we use the AIC and BIC:

$$\begin{array}{c}\hfill \mathrm{AIC}=-2logp\left(\mathit{y}\right|\widehat{\mathit{\theta }})+2d^{\prime },\end{array}$$

(10)

$$\begin{array}{c}\hfill \mathrm{BIC}=-2logp\left(\mathit{y}\right|\widehat{\mathit{\theta }})+d^{\prime }logn,\end{array}$$

(11)

where $d^{\prime }$ is the dimension of parameters of the model under consideration. Moreover, we use other criteria based on the $L_q$-likelihood:

$$\begin{array}{c}\hfill L_q\text{-AIC}=-2L_{q}p\left(\mathit{y}\right|\widehat{\mathit{\theta }})+2d^{\prime },\end{array}$$

(12)

$$\begin{array}{c}\hfill L_q\text{-BIC}=-2L_{q}p\left(\mathit{y}\right|\widehat{\mathit{\theta }})+d^{\prime }logn.\end{array}$$

(13)

For a sequence $\left({\widehat{\mathit{\theta }}}_{\left(k\right)}\right)$ made by each of the methods, let $I_{\left(k\right)}=\left\{i\right|{\widehat{\theta }}_{\left(k\right)}^i\ne 0\}$ and ${\widehat{\mathit{\theta }}}_{\mathrm{MLE}}^{\left(k\right)}$ the MLE of the model ${\mathcal{M}}_{\left(k\right)}=\left\{p(·|\mathit{\theta })\right|{\theta }^j=0(j\notin I_{\left(k\right)})\}$. We call (10) with $\widehat{\mathit{\theta }}={\widehat{\mathit{\theta }}}_{\mathrm{MLE}}^{\left(k\right)}$ AIC1, and (10) with $\widehat{\mathit{\theta }}={\widehat{\mathit{\theta }}}_{\left(k\right)}$ AIC2. Similarly, (11) with $\widehat{\mathit{\theta }}={\widehat{\mathit{\theta }}}_{\mathrm{MLE}}^{\left(k\right)}$ is BIC1, and (11) with $\widehat{\mathit{\theta }}={\widehat{\mathit{\theta }}}_{\left(k\right)}$ is BIC2. The $L_q$-AIC and $L_q$-BIC are referred to in the same manner; for example, (12) with $\widehat{\mathit{\theta }}={\widehat{\mathit{\theta }}}_{\mathrm{MLE}}^{\left(k\right)}$ is $L_q$-AIC1. Note that AIC1, BIC1, $L_q$-AIC1, and $L_q$-BIC1 are available only when the MLE exists; AIC2, BIC2, $L_q$-AIC2, and $L_q$-BIC2 are always applicable. Finally, we used cross-validation (CV) in addition to these information criteria.

4.2. Result

The results are presented in Figure 1, Figure 2, Figure 3, Figure 4, Figure 5, Figure 6, Figure 7, Figure 8, Figure 9, Figure 10, Figure 11, Figure 12, Figure 13, Figure 14, Figure 15, Figure 16, Figure 17, Figure 18, Figure 19, Figure 20, Figure 21 and Figure 22, which report the best result for each method based on the various information criteria. We present the tables of the results of the numerical experiments in the Supplementary Material. In the figures, white bars represent LASSO, gray bars represent SCAD, and black bars represent MCP.

The model selection results are reported in Figure 1, Figure 2, Figure 3, Figure 4, Figure 5, Figure 6, Figure 7, Figure 8, Figure 9, Figure 10, Figure 11, Figure 12, Figure 13 and Figure 14. The vertical axis indicates the number of trials (among $m=1000$ trials) where a method selects the true model. Here, a larger value is better. The horizontal axis shows the value of ${\theta }_0$.

The generalization results are reported in Figure 15, Figure 16, Figure 17, Figure 18, Figure 19, Figure 20, Figure 21 and Figure 22. To evaluate the generalization error of the proposed methods, we newly make $m=1000$ independent copies $\{({\mathit{y}}_1^{\prime },X_1^{\prime }),\cdots ,({\mathit{y}}_m^{\prime },X_m^{\prime })\}$ in each trial. We computed the difference between $({\mathit{y}}_1^{\prime },\cdots ,{\mathit{y}}_m^{\prime })$ and the m predictions using each of the methods. The vertical axis indicates the average prediction error over m trials. In this case, a smaller value is better. The scaling of Figure 21 and Figure 22 $(q=5/3)$ is different from that of the other figures.

Our first concern is whether the proposed methods work well. The results for $q=1$ can be regarded as a reference for the other values of q. The figures show that the proposed methods work well in both model selection and generalization, especially for $q<2$. The methods also perform well in terms of model selection for $q=2,2.01$, and $2.1$. However, they perform poorly for $q=2.5$ in terms of model selection and for $q\ge 2$ in terms of generalization. As anticipated, a large q makes the problem difficult.

Second, we evaluate the performance of the proposed methods, finding that the MCP performs best in most cases. In a few cases, the MCP performed similarly to or slightly worse than the other methods. For model selection, the cases with $q=1,n=1000$ and large ${\theta }_0$ are exceptions. Furthermore, the LASSO performed worse than the SCAD and MCP.

Third, we consider the effect of $r_{\mathrm{nz}}$, ${\theta }_0$, d, and n, in addition to q. The cases with large $r_{\mathrm{nz}}$ and/or small ${\theta }_0$ are difficult. Moreover, a large d makes the problems difficult. However, if we have a small q ($1\le q<2$), large ${\theta }_0$ (${\theta }_0={10}^2,{10}^3$) and small $r_{\mathrm{nz}}$, the problems with large d can be easier than those with small d. Furthermore, a small n makes the problems difficult in a similar manner to a large d. These observations imply that, for $1\le q<2$, small-sample problems can be easier than large-sample problems if $r_{\mathrm{nz}}$ is small and ${\theta }_0$ is large.

Fourth, the choice of information criterion changes the methods’ performance. In terms of model selection, BIC2 was mostly the best for many values of q. For $3/2\le q\le 2.1$, BIC1 was a little better than BIC2 if BIC1 was available. For $q=1$ and $13/11$, BIC2 was better than BIC1. AIC1 and AIC2 were as good as BICs for $2\le q\le 2.1$. Moreover, the $L_q$-BIC1 and -BIC2 were best only for $q=3/2$, when BIC1 and BIC2 performed just as well. Overall, the $L_q$-information criteria performed poorly.

Furthermore, in terms of generalization, BIC2 was mostly the best. AIC2 was as good as BIC2, whereas AIC2 was sometimes a little worse than BIC2. The information criteria using the $L_q$-likelihood were poor for $q=13/11$. For $q=1,3/2$, and $5/3$, the $L_q$-information criteria worked as well as the ordinary criteria and CV, except for some cases. The performance of CV was mostly good, but was occasionally very poor.

In summary, using an appropriate criterion, the proposed methods perform well for linear models with slightly heavy-tailed errors ($1\le q<2$). Moreover, the proposed methods work in terms of model selection, even if the error is heavy-tailed ($2\le q<2.5$). Overall, we recommend using the MCP and BIC2.

5. Conclusions

We proposed regularization methods for q-normal linear models based on the $L_q$-likelihood. The proposed methods coincide with the ordinary regularization methods. Our methods perform well for slightly heavy-tailed errors ($1\le q<2$) in terms of model selection and generalization. Moreover, they work well in terms of model selection for heavy-tailed errors ($2\le q<2.5$). A theoretical analysis of the proposed methods is left to future work.