text:
\begin{document}

\begin{frontmatter}

\title{Local polynomial regression and variable selection}

\runtitle{Local variable selection}

\begin{aug}
\author{\fnms{Hugh} \snm{Miller}\ead[label=e1]{h.miller@ms.unimelb.edu.au}}
\and
\author{\fnms{Peter} \snm{Hall}\ead[label=e2]{halpstat@ms.unimelb.edu.au}}

\runauthor{H. Miller and P. Hall}

\affiliation{The University of Melbourne}

\address[a]{Department of Mathematics and Statistics, The University of Melbourne, Melbourne, VIC 3052, Australia, \printead{e1,e2}}

\end{aug}

\begin{abstract}
We propose a method for incorporating variable selection into local polynomial regression. This can improve the accuracy of the regression by extending the bandwidth in directions corresponding to those variables judged to be are unimportant. It also increases our understanding of the dataset by highlighting areas where these variables are redundant. The approach has the potential to effect complete variable removal as well as perform partial removal when a variable redundancy applies only to particular regions of the data. We define a nonparametric oracle property and show that this is more than satisfied by our approach under asymptotic analysis. The usefulness of the method is demonstrated through simulated and real data numerical examples.
\end{abstract}

\begin{keyword}[class=AMS]
\kwd[Primary ]{62G08}
\kwd[; secondary ]{62G20}
\end{keyword}

\begin{keyword}
\kwd{variable selection}
\kwd{local regression}
\kwd{adaptive bandwidth}
\kwd{local variable significance}
\end{keyword}

\end{frontmatter}

\section{Introduction}

The classical regression problem is concerned with predicting a noisy continuous response using a $d$-dimensional predictor vector with support on some $d$-dimensional subspace. This functional relationship is often taken to be smooth and methods for estimating it range from parametric models, which specify the form of the relationship between predictors and response, through to nonparametric models, which have fewer prior assumptions about the shape of the fit. An important consideration for fitting such a regression model is whether all $d$ predictors are in fact necessary. If a particular predictor has no relationship to the response, the model be made both simpler and more accurate by removing it. In recent years there has been strong interest in techniques that automatically generate such ``sparse'' models. Most attention has been given to parametric forms, and in particular the linear model, where the response is assumed to vary linearly with the predictors. There has also been some investigation into variable selection for nonlinear models, notably through the use of smoothing splines and local regression.

One common feature of the existing sparse methods is that the variable selection is ``global'' in nature, attempting to universally include or exclude a predictor. Such an approach does not naturally reconcile well with some nonparametric techniques, such as local polynomial regression, which focus on a ``local'' subset of the data to estimate the response. In this local context it would be more helpful to understand local variable influence, since predictors that are irrelevant in some regions may in fact be important elsewhere in the subspace. Just as in the global setting, such information would allow us to improve the accuracy and parsimony of a model, but at a local level.

However this approach to variable selection can be problematic. Most notably, variable significance  affects the definition of ``local''. To illustrate concretely, suppose that two data points were close in every dimension except one. In typical local regression these points would not be considered close, and so the response at one point would not impact the other. If, however, we establish that the one predictor they differ by is not influential over a range that includes both these points, then they should actually be treated as neighbouring, and be treated as such in the model. Any methodology seeking to incorporate local variable influence needs to accommodate such potential situations.

Understanding local variable significance can also give additional insight into a dataset. If a variable is not important in certain regions of the support, knowledge of this allows us to discount it in certain circumstances, simplifying our understanding of the problem. For example, if none of the variables are relevant in a region, we may treat the response as locally constant and so know that we can ignore predictor effects when an observation lies in this region.

A final consideration is theoretical performance. In particular we shall present on approach that is ``oracle''; that is, its performance is comparable to that of a particularly well-informed statistician, who has been provided in advance with the correct variables. It is interesting to note that variable interactions often cause sparse parametric approaches to fail to be oracle, but in the local nonparametric setting this is not an issue, because such interactions vanish as the neighbourhood of consideration shrinks.

In this paper we propose a flexible and adaptive approach to local variable selection using local polynomial regression. The key technique is careful adjustment of the local regression bandwidths to allow for variable redundancy. The method has been named LABAVS, standing for ``locally adaptive bandwidth and variable selection''. Section~2 will introduce the LABAVS algorithm, including a motivating example and possible variations. Section~3 will deal with theoretical properties and in particular establishes a result showing that the performance of LABAVS is better than oracle when the dimension remains fixed. Section~4 presents numerical results for both real and simulated data, showing that the algorithm can improve prediction accuracy and is also a useful tool in arriving at an intuitive understanding of the data. Technical details have been relegated to an appendix which may be found in the long version of this paper (Miller and Hall, 2010).

LABAVS is perhaps best viewed as an improvement to local polynomial regression, and will retain some of the advantages and disadvantages associated with this approach. In particular, it still suffers the ``curse of dimensionality,'' in that it struggles to detect local patterns when the dimension of genuine variables increases beyond a few. It is not the first attempt at incorporating variable selection into local polynomial regression; the papers by Lafferty and Wasserman (2008) and Bertin and Lecu\'{e} (2008) also do this. We compare our approach to these in some detail in Section~\ref{se:compare}. LABAVS can also be compared to other nonparametric techniques in use for low to moderate dimensions. These include generalised additive models, MARS and tree based methods (see Hastie \etalc 2001).

The earliest work on local polynomial regression dates back to that of Nadaraya (1964) and Watson (1964). General references on the subject include Wand and Jones (1995), Simonoff (1996) and Loader (1999). An adaptive approach to bandwidth selection may be found in Fan and Gijbels (1995), although this was not in the context of variable selection. Tibshirani (1996) studies the LASSO, one of the most popular sparse solutions for the linear model; more recent related work on the linear model includes that of Candes and Tao (2007) and Bickel \etal (2009). Zou (2006) created the adaptive version of the LASSO and proved oracle performance for it. Lin and Zhang (2006) and Yuan and Lin (2006) have investigated sparse solutions to smoothing spline models.

The LABAVS algorithm also bears some similarity to the approach adopted by Hall \etal (2004). There the aim was to estimate the conditional density of a response using the predictors. Cross-validation was employed and the bandwidths in irrelevant dimensions diverged, thereby greatly downweighting those components. In the present paper the focus is more explicitly on variable selection, as well as attempting to capture local variable dependencies.

\section{Model and Methodology}

\subsection{Model and definitions}

Suppose that we have a continuous response $Y_i$ and a $d$-dimensional random predictor vector $X_i = (X^{(1)}_i, \ldots, X^{(d)}_i)$ which has support on some subspace $\cC \subset \bR^d$. Further, assume that the observation pairs $(Y_i,X_i)$ are independent and identically distributed for $i=1,\ldots,n$, and that $X_i$ has density function $f$. The response is related to the predictors through a function $g$,
\beqn Y_i = g(X_i) + \ep_i \,, \label{eq:simple} \neqn
with the error $\ep_i$ having zero mean and fixed variance. Smoothness conditions for $f$ and $g$ will be discussed in the theory section.

Local polynomial regression makes use of a kernel and bandwidth to assign increased weight to neighbouring observations compared to those further away, which will often have zero weight. We take $K(u) = \prod_{1 \leq j \leq d} K^*(u^{(j)})$ to be the $d$-dimensional rectangular kernel formed from a one dimensional kernel $K^*$ such as the tricubic kernel,
$$ K^*(u^{(j)})=(35/32)(1-x^2)^3 \bI(|x|<1) \,.$$
\label{kernel} Assume $K^*$ is symmetric with support on $[-1,1]$. For $d \times d$ bandwidth matrix $H$ the kernel with bandwidth $H$, denoted $K_H$, is
\beqn K_H(u) = \frac{1}{|H|^{1/2}}K(H^{-1/2}u) \,. \label{eq:ker} \neqn
We assume that the bandwidth matrices are diagonal, $H = \diag(h_1^2,\ldots,h_d^2)$, with each $h_j>0$, and write $H(x)$ when $H$ varies as a function of $x$. Asymmetric bandwidths can be defined as having both a lower and an upper (diagonal) bandwidth matrix,  $H^L$ and $H^U$ respectively, for a given estimation point $x$, rather than a single bandwidth $H$ for all $x$. The kernel weight of an observation $X_i$ at estimation point $x$ with asymmetrical bandwidth matrices $H^L(x)$ and $H^U(x)$, is
\begin{eqnarray}
K_{H^U(x),H^L(x)}(X_i - x) &=& \prod_{j \,:\, X_i^{(j)} < x^{(j)}} \frac{1}{h^L_j(x)} K^*\left( \frac{X^{(j)}_i - x^{(j)}}{h^L_j(x)} \right) \\
&& \qquad \cdot \prod_{j\,:\, X_i^{(j)} \ge x^{(j)}} \frac{1}{h^U_j(x)} K^*\left( \frac{X^{(j)}_i - x^{(j)}}{h^U_j(x)} \right) \,. \nonumber
\end{eqnarray}
This amounts to having (possibly) different window sizes above and below $x$ in each direction. Although such unbalanced bandwidths would often lead to undesirable bias properties in local regression, here they will be used principally to extend bandwidths in dimensions considered redundant, so this issue is not a concern.

We also allow the possibility of infinite bandwidths $h_j = \infty$. In calculating the kernel in (\ref{eq:ker}) when $h_j$ is infinite, proceed as if the $j$th dimension did not exist (or equivalently, as if the $j$th factor in rectangular kernel product is always equal to 1). If all bandwidths are infinite, consider the kernel weight to be 1 everywhere. Although the kernel and bandwidth conditions above have been defined fairly narrowly to promote simplicity in exposition, many of these assumptions are easily generalised.

Local polynomial regression estimates of the response at point $x$, $\hg(x)$, are found by fitting a polynomial $q$ to the observed data, using the kernel and bandwidth to weight observations. This is usually done by minimising the weighted sum of squares,
\beqn \sum_{i=1}^n \{Y_i - q(X_i-x)\}^2 K_H(X_i-x) \,. \label{eq:locpoly} \neqn
Once the minimisation has been performed, $q(0)$ becomes the point estimate for $g(x)$. The polynomial is of some fixed degree $p$, with larger values of $p$ generally decreasing bias at the cost of increased variance. Of particular interest in the theoretical section will be the local linear fit, which minimises
\beqn  \sum_{i=1}^n \left[ \left\{ Y_i - \ga_0 - \sum_{j=1}^n (X_i^{(j)}-x^{(j)})\ga_j   \right\}^2 K_{H}(X_i-x) \right] \,,\label{eq:loclin} \neqn
over $\ga_0$ and $\ga=(\ga_1,\ldots,\ga_d)$.

\subsection{The LABAVS Algorithm}

Below is the LABAVS algorithm that will perform local variable selection and vary the bandwidths accordingly. The choice of $H$ in the first step can be local or global and should be selected as for a traditional polynomial regression, using cross-validation, a plug-in estimator or some other standard technique. Methods for assessing variable significance in Step~2, and the degree of shrinkage needed in Step~4, are discussed below.

$$\begin{array}{p{13cm}}
\ni {\sl LABAVS Algorithm}
\begin{enumerate}
\item Find a starting $d\times d$ bandwidth $H = \diag(h^2,\ldots,h^2)$.
\item For each point $x$ of a representative grid in the data support, perform local variable selection to determine disjoint index sets $\chA^+(x),\chA^-(x)$, with $\chA^+(x) \cup \chA^-(x) = \{1,\ldots,d\}$,  for variables that are considered relevant and redundant respectively.
\item For any given $x$, derive new local bandwidth matrices $H^L(x)$ and $H^U(x)$ by extending the bandwidth in each dimension indexed in $\chA^-(x)$. The resulting space given nonzero weight by the kernel $K_{H^L(x),H^U(x)}(u-x)$ is the rectangle of maximal area with all grid points $x_0$ inside the rectangle satisfying $\chA^+(x_0) \subset \chA^+(x)$. Here $\chA^+(x)$ is calculated explicitly as in Step~2, or taken as the set corresponding the closest grid point to $x$.
\item Shrink the bandwidth slightly for those variables in $\chA^+(x)$ according to the amount that bandwidths have increased in the other variables.
\item Compute the local polynomial estimator at $x$, excluding variables in $\chA^-(x)$ and using adjusted assymetrical bandwidths $H^L(x)$ and $H^U(x)$. The expression to be minimised is
$$ \sum_{i=1}^n \left\{ Y_i - q(X_i-x)   \right\}^2 K_{H^L(x), H^U(x)}(X_i-x)  \,, $$
where the minimisation runs over all polynomials $q$ of appropriate degree. The value of $q(0)$ in the minimisation is the final local linear estimator.
\end{enumerate}
\end{array}$$

The key feature of the algorithm is that variable selection directly affects the bandwidth, increasing it in the direction of variables that have no influence on the point estimator. If a variable has no influence anywhere, it has the potential to be completely removed from the local regression, reducing the dimension of the problem. For variables that have no influence in certain areas, the algorithm achieves a partial dimension reduction. The increased bandwidths reduce the variance of the estimate and Step~4 swaps some of this reduction for a decrease in the bias to further improve the overall estimator.

\label{se:eg}

As a concrete example of the approach, define the following one-dimensional ``huberised'' linear function:
\beqn g(x) = x^2\bI(0 < x \le 0.4) + (0.8x-0.16)\bI(x>0.4) \,,\neqn
and let $g(X) = g\{([X^{(1)}]_+^2+[X^{(2)}]_+^2)^{1/2}\}$ for 2-dimensional random variable $X = (X^{(1)}, X^{(2)})$. Assume that $X$ is uniformly distributed on the space $[-2,2]\times[-2,2]$. Notice that when $X^{(1)}$, $X^{(2)} <0$ the response variable $Y$ in (\ref{eq:simple}) is independent of $X^{(1)}$ and $X^{(2)}$; when $X^{(1)} < 0$ and $X^{(2)} >0$ the response depends on  $X^{(2)}$ only; when $X^{(1)} > 0$ and $X^{(2)} <0$ the response depends on  $X^{(1)}$ only; when $X^{(1)}$, $X^{(2)} >0$ the response depends on both $X^{(1)}$ and $X^{(2)}$. Thus in each of these quadrants a different subset of the predictors is significant. A local approach to variable significance can capture these different dependencies, while a global variable redundancy test would not eliminate any variables.

\begin{figure}[htbp] \begin{center}
\includegraphics[width = 3.5in, height=3.5in]{bandwidths2.jpg}
\vspace{-8pt}
{\it \caption{ Bandwidth adjustments under ideal circumstances in illustrative example. \label{fig2}}}
\end{center} \end{figure}

Now consider how the algorithm applies to this example, starting with a uniform initial bandwidth of $h=0.5$ in both dimensions. Assuming that variable significance perfectly on a dense grid, figure~\ref{fig2} illustrates the adjusted bandwidths for each of the quadrants. The dots are four sample estimation points, the surrounding unit squares indicate the initial bandwidths and the dashed lines indicate how the bandwidths are modified. In the bottom left quadrant both variables are considered redundant, and so the bandwidth expands to cover the entire quadrant. This is optimal behaviour, since the true function is constant over this region, implying that the best estimator will be produced by including the whole area. In the bottom right quadrant the first dimension is significant while the second is not. Thus the bandwidth for the second dimension is ``stretched'', while the first is shrunken somewhat. Again, this is desirable for improving the estimator. The stretching in the second dimension improves the estimator by reducing the variance as more points are considered. Then the shrunken first dimension swaps some of this reduction in variance for decreased bias. Finally, in the top right quadrant, there is no change in the bandwidth since both variables are considered to be significant.

\subsection{Variable selection step} \label{se:approaches}
Below are three possible ways to effect variable selection at $x_0$ in Step~2 of the algorithm, presented in the context of local linear regression. They all make use of a tuning parameter $\la$ which controls how aggressive the model is in declaring variables as irrelevant. Cross validation can be used to select an appropriate level for $\la$. So that the tuning parameters are comparable at different points in the data domain, it is useful to consider a local standardisation of the data at $x_0$. Define $\bX_{x_0} = (\bX_{x_0}^{(1)}, \ldots, \bX_{x_0}^{(d)})$ and $\bY_{x_0}$ by
\begin{eqnarray*}
 \bX_{x_0}^{(j)} &=& \frac{\sum_{i=1}^n X_i^{(j)} K_H(X_i-x_0)}{ \sum_{i=1}^n K_H(X_i-x_0)} \,, \quad \bY_{x_0} = \frac{\sum_{i=1}^n Y_i K_H(X_i-x_0)}{ \sum_{i=1}^n K_H(X_i-x_0)} \,,
\end{eqnarray*}
 and define $\tX_i = (\tX_i^{(1)},\ldots,\tX_i^{(d)})$ and $\tY_i$ by
\begin{eqnarray}
 \tX_i^{(j)} &=& \frac{(X^{(j)}_i-\bX^{(j)}_{x_0})K_H(X_i-x_0)^{1/2}} {\big\{\sum_{i=1}^n \big(X_i^{(j)} - \bX^{(j)}_{x_0} \big)^2 K_H(X_i-x_0) \big\}^{1/2}} \,, \label{eq:scale} \\
 \tY_i &=& (Y_i-\bY_{x_0})K_H(X_i-x_0)^{1/2} \,. \nonumber
\end{eqnarray}
Notice that $\tX$ and $\tY$ incorporate the weight $K_H(X_i-x_0)$ into the expression.

\begin{enumerate}
\item {\bf Hard thresholding:} Choose parameters to minimise the weighted least squares expression,
\beqn \sum_{i=1}^n \left\{ \tY_i - \be_0 - \sum_{j =1}^n \tX_i^{(j)}\be_j   \right\}^2  \,, \label{eq:appr1} \neqn
and classify as redundant those variables for which $|\hbe_j|<\la$. This can be extended to higher degree polynomials, although performance tends to be more unstable.
\item {\bf Backwards stepwise approach:} For each individual $j$, calculate the percentage increase in the sum of squares if the $j$th variable is excluded from the local fit. Explicitly, if $\hq$ is the optimal local fit using all variables and $\hq_j$ is the fit using all except the $j$th, we classify the $j$th variable as redundant if
\beqn \frac{\sum_{i=1}^n \left\{ Y_i - \hq_j(\hat{X}_i)   \right\}^2 K_H(\hat{X}_i) - \sum_{i=1}^n \left\{ Y_i - \hq(\hat{X}_i)   \right\}^2 K_H(\hat{X}_i)}{\sum_{i=1}^n \left\{ Y_i - \hq(\hat{X}_i)   \right\}^2 K_H(\hat{X}_i)} < \la \,, \label{eq:appr2} \neqn
where $\hat{X}_i =X_i-x_0$.
\item {\bf Local lasso:} Minimise the expression
\beqn \sum_{i=1}^n  \left\{ \tY_i - \ga_0 - \sum_{j =1}^n \tX_i^{(j)}\ga_j   \right\}^2  + \la \sum_{j=1}^d|\ga_j| \,. \label{eq:appr3} \neqn
Those variables for which $\ga_j$ are set to zero in this minimisation are then classed as redundant. While the normal lasso can have consistency problems (Zou, 2006), this local version does not since variables are asymptotically independent as $h\to 0$. The approach also scales naturally to higher order polynomials, provided all polynomial terms are locally standardised; a variable is considered redundant if all terms that include it have corresponding parameters set to zero by the lasso.
\end{enumerate}

We have found that the first and second of the above approaches have produced the most compelling numerical results. The numerical work in Section~\ref{se:num} uses the first approach for linear polynomials, while the theoretical work in Section~\ref{theory} establishes uniform consistency for both of the first two methods, guaranteeing oracle performance.

\subsection{Variable shrinkage step}

The variable shrinkage step depends on whether the initial bandwidth, and thus the shrunken bandwidth $h'$, is chosen locally or globally. Define
\beqn V(x,H) = \frac{\tsum K_{H}(X_i-x)^2 }{ \{\tsum K_{H}(X_i-x) \}^2 } \,, \neqn
where the bandwidth term in the function $V$ is allowed to be asymmetrical. Then in the local case, letting $d'(x)$ denote the cardinality of $\chA(x)$, let
\beqn M(x) = V[x,\{H^L(x),H^L(x)\}] / V(x,H)  \label{eq:fac} \,. \neqn
The expression is asymptotically proportional to $\{h'(x)\}^{-d'(x)}$ and estimates the degree of variance stabilisation resulting from the bandwidth adjustment. Using this, the correct amount of bandwidth needed in step~4 is $h'(x) = h \{M(x) d'(x)/d\}^{1/4} $. Since both sides of this expression depend on $h'(x)$, shrinkage can be approximated in the following way. Let
$$ M^*(x) = V[x,\{\tilde{H}^L(x),\tilde{H}^L(x)\}] / V(x,H) \,,$$
where $\tilde{H}^L(x)$ and $\tilde{H}^U(x)$ are the bandwidth matrices immediately after step 3. Then the shrunken bandwidths are $h'(x) = h\{ M^*(x) d'(x)/d )^{1/(d'(x)+4)}$.

In the global bandwidth case, we define
\beqn M[\{H^L(X),H^U(X)\},H] = \frac{E\left(V[X,\{H^L(X),H^L(X)\}]\right)}{E\{V(X,H)\}} \,. \label{newfac1} \neqn
This expression measures the average variance stabilisation across the domain. In this case, the shrinkage factor should satisfy
\beqn h' = h \left( M[\{H^L(X),H^U(X)\},H] E\{d'(X)\}/d \right)^{1/4} \,. \label{newfac2} \neqn
The theoretical properties in Section~\ref{theory} deal with the global bandwidth scenario. The treatment for the local case is similar, except that care must be taken in regions of the domain where the function $g$ behaves in a way that is exactly estimable by a local polynomial and thus has potentially no bias.

\subsection{Further remarks}

\begin{enumerate}
\item The choice of distance between grid points in Step~2 is somewhat arbitrary, but should be taken as less than $h$ so that all data points are considered in calculations. In the asymptotic theory we let this length decrease faster that the rate of the bandwidth, and in numerical experimentation the choice impacts only slightly on the results.
\item Step 5 of the algorithm forces the estimate at point $x$ to exclude variables indexed in $\chA^-(x)$. An alternative is to still use all variables in the final fit. This may be advantageous in situations with significant noise, where variable admission and omission is more likely to have errors. Despite including these extra variables, the adjusted bandwidths still ensure that estimation accuracy is increased.
\item Finding the maximal rectangle for each representative point, as suggested in step~3 of the algorithm, can be a fairly intensive computational task. In our numerical work we simplified this by expanding the rectangle equally until the boundary met a ``bad'' grid point (i.e.~a point $x'$ such that $\chA^+(x')\nsubseteq \chA^+(x)$). The corresponding direction was then held constant while the others continue to increase uniformly. We continued until each dimension stopped expanding or grew to be infinite. This approach does not invalidate the asymptotic results in Section~\ref{theory}, but there may be some deterioration in numerical performance associated with this simplification.
\item If a variable is redundant everywhere, results in Section~\ref{theory} demonstrate that the algorithm is consistent; the probability that the variable is classified as redundant everywhere tends to 1 as $n$ grows. However, the exact probability is not easy to calculate and for fixed $n$ we may want greater control over the ability to exclude a variable completely. In such circumstances a global variable selection approach may be appropriate.
\item As noted at the start of Section~2.2, the initial bandwidth in Step~1 does not necessarily have to be fixed over the domain.  For instance, a nearest neighbour bandwidth, where $h$ at $x$ is roughly proportional to $f(x)^{-1}$, could be used. Employing this approach offers many practical advantages and the theoretical basis is similar to that for the constant bandwidth. The numerical work makes use of nearest neighbour bandwidths throughout. In addition, we could use an initial bandwidth that was allowed to vary for each variable, $H = \diag(h_1^2,\ldots,h_p^2)$. So long as, asymptotically, each $h_j$ was equal to $C_j h$ for some controlling bandwidth $h$ and constant $C_j$, the theory would hold, although details are not pursued here.
\end{enumerate}

\subsection{Comparison to other local variable selection approaches} \label{se:compare}

As mentioned in the introduction, two recent papers take a similar approach to this problem. Firstly Lafferty and Wasserman (2008) introduce the rodeo procedure. This attempts to assign adaptive bandwidths based on the derivative with respect to the bandwidth for each dimension, $\partial \hg(x)/\partial h_j$. This has the attractive feature of bypassing the actual local shape and instead focussing on whether an estimate is improved by shrinking the bandwidths. It is also a greedy approach, starting with large bandwidths in each direction and shrinking only those that cause a change in the estimator at a point. The second paper is by Bertin and Lecu\'{e} (2008), who implement a two step procedure to reduce the dimensionality of a local estimate. The first step fits a local linear estimate with an $L_1$ or lasso type penalty, which identifies the relevant variables. This is followed by a second local linear fit using this reduced dimensionality. The lasso penalty they use is precisely the same as the third approach suggested in Section~\ref{se:approaches}.

We comment on the similarities and differences of these two approaches compared to the current presentation, which are summarised in Table~\ref{tabcomp}. Firstly the theoretical framework of the two other papers focus exclusively on the performance at a single point, while the LABAVS approach ensures uniformly oracle performance on the whole domain. The framework for the other two also assumes that variables are either active on the whole domain or redundant everywhere, while we have already discussed the usefulness of an approach that can adapt to variables that are redundant on various parts of the data. We believe this is particularly important, since local tests of variable significance will give the same results everywhere. Related to this, our method does not require an assumption of nonzero gradients (whether with respect to the bandwidth or variables) to obtain adequate theoretical performance, in contrast to the other methods. On the other hand, ensuring both uniform performance while allowing $d$ to be increasing is a quite challenging, so our presentation assumes $d$ is fixed,  while the others do not. It is also worth noting that the greedy approach of Lafferty and Wasserman potentially gives it an advantage in higher dimensional situations.

While all approaches work in a similar framework, the above discussion demonstrates that there are significant differences. Our methodology may be viewed as a generalisation of the work of Bertin and Lecu\'{e}, save for imposing fixed dimensionality. It can also be viewed as a competitor to the rodeo, and some numerical examples comparing the two are provided.

\begin{table}[htbp]\begin{center}
\caption{ Summary of locally adaptive bandwidth approaches  \label{tabcomp}}
\begin{tabular}{|l|c|c|c|} \hline
  &  LABAVS \T& Rodeo &  Bertin and \\
  &         &       &  Lecu\'{e} (2008) \B\\ \hline
Oracle performance on entire domain \T \B& \tick & \cross & \cross \\ \hline
Allows for locally redundant variables \T \B& \tick & \cross & \cross \\ \hline
Relevant variables allowed to have zero gradient \T \B& \tick & \cross & \cross \\ \hline
Theory allows dimension $d$ to increase with $n$ \T \B& \cross & \tick & \tick \\ \hline
Greedy algorithm applicable for higher dimensions \T \B& \cross & \tick & \cross \\ \hline
\end{tabular} \end{center} \end{table}

With regards to computation time, for estimation at a single point the rodeo is substantially faster, since calculating variable significance on a large grid of points is not required. If however we need to make predictions at a reasonable number of points, then Labavs is likely to be more efficient, since the grid calculations need only be done once, while rodeo requires a new set of bandwidth calculations for each point.

\section{Theoretical properties} \label{theory}

As mentioned in the introduction, a useful means of establishing the power of a model that includes variable selection is to compare it with an oracle model, where the redundant variables are removed before the modelling is undertaken. In the linear (and the parametric) context, we interpret the oracle property as satisfying two conditions as $n \to \infty$:
\begin{enumerate}
    \item the probability that the correct variables are selected converges to 1, and
    \item the nonzero parameters are estimated at the same asymptotic rate as they would be if the correct variables were known in advance.
\end{enumerate}
We wish to extend this notion of an oracle property to the nonparametric setting, where some predictors may be redundant. Here there are no parameters to estimate, so attention should instead be given to the error associated with estimating $g$. Below we define weak and strong forms of these oracle properties:
\begin{defn}
The {\it weak} oracle property in nonparametric regression is:
\begin{enumerate}
    \item the probability that the correct variables are selected converges to 1, and
    \item at each point $x$ the error of the estimator $\hg(x)$ decreases at the same asymptotic rate as it would if the correct variables were known in advance.
\end{enumerate}
\end{defn}
\begin{defn}
The {\it strong} oracle property in nonparametric regression is:
\begin{enumerate}
    \item the probability that the correct variables are selected converges to 1, and
    \item at each point $x$ the error of the estimator $\hg(x)$ has the same first-order asymptotic properties as it would if the correct variables were known in advance.
\end{enumerate}
\end{defn}
Observe that the weak oracle property achieves the correct rate of estimation while the strong version achieves both the correct rate and the same asymptotic distribution. The first definition is most analogous to its parametric counterpart, while the second is more ambitious in scope.

Here we establish the strong version of the nonparametric oracle property for the LABAVS algorithm, with technical details found in the appendix (Miller and Hall, 2010). We shall restrict attention to the case of fixed dimension. To apply to a situation of increasing dimension, we could add an asymptotically consistent screening method to reduce it back to fixed $d$. The treatment here focuses on local linear polynomials, partly for convenience but also recognising that the linear factors dominate higher order terms in the asymptotic local fit. Thus our initial fit is found by minimising the expression (\ref{eq:loclin}). We impose further conditions on the kernel $K$:
\beqn \begin{array}{p{12.7cm}}
$\int K(z)dz = 1$, $\int z^{(j)}K(z)dz = 0$ for each $j$,  $\int z^{(j)}z^{(k)} K(z)dz = 0$ when $j \neq k$ and $\int (z^{(j)})^2 K(z)dz = \mu_2(K)>0$, with  $\mu_2(K)$ independent of $j$. \label{eq:c4}
\end{array} \neqn
The useful quantity $R(K)$, depending on the choice of kernel, is defined as
$$ R(K) = \int K(z)^2 dz = \left\{ \int K^*(z^{(j)})^2 dz^{(j)}\right\}^d \,, $$
where $K^*$ is the univariate kernel introduced on page \pageref{kernel}. Let $a_n \asymp b_n$ denote the property that $a_n = O(b_n)$ and $b_n = O(a_n)$. We also require the following conditions (\ref{eq:cond}), needed to ensure uniform consistency of our estimators.
\beqn  \begin{array}{p{12cm}} \label{eq:cond}
\begin{enumerate}
\item  The support $\cC = \{x:f(x)>0\}$ of the random variable $X$ is compact. Further, $f$ and its first order partial derivatives are bounded and uniformly continuous on $\cC$, and $\inf_{x\in \cC} f(x) > 0$.
\item The kernel function $K$ is bounded with compact support and satisfies $|p(u) K(u) - p(v) K(v)| \le C_1||u-v||$ for some $C_1 > 0$ and all points $u,v$ in $\cC$. Here $p(u)$ denotes a single polynomial term of the form $\prod (u^{(j)})^{a_j}$ with the nonnegative integers $a_j$ satisfying $\sum a_j \le 4$. The bound $C_1$ should hold for all such choices of $p$.
\item The function $g$ has bounded and uniformly continuous partial derivatives up to order 2. If $(D^k g)(x)$ denotes the partial derivative
$$ \frac{\partial^{|k|}g(x)}{\partial (x^{(1)})^{k_1} \cdots \partial (x^{(d)})^{k_d}} \,,$$
with $|k| = \sum k_j$, then we assume that these derivatives satisfy, for some constant $C_2$,
$$ |h(u) - h(v)| \le C_2 || u - v || \,.$$
\end{enumerate}
\end{array} \neqn
$$ \begin{array}{p{12cm}}
\begin{enumerate}
\setcounter{enumi}{3}
\item $E(|Y|^\si) < \infty$ for some $\si > 2$.
\item The conditional density $f_{X|Y}(x|y)$ of $X_i$, conditional on $Y$, exists and is bounded.
\item $$ \mbox{For some $0<\xi<1$, } \frac{n^{1-2/\si} h^d}{\log n\{ \log n (\log\log n)^{1+\de}\}^{2/\si}} \to \infty \,.$$
\item The Hessian of $g$, $\cH_g$, is nonzero on a set of nonzero measure in $\cC$.
\end{enumerate}
\end{array} $$
The conditions in (\ref{eq:cond}), except perhaps the first, are fairly natural and not overly constrictive. For example, the sixth will occur naturally for any reasonable choice of $h$, while the second follows easily if $K$ has a bounded derivative. The last condition is purely for convenience in the asymptotics; if $\cH_g$ was zero almost everywhere then $g$ would be linear and there would be no bias in the estimate, improving accuracy. The first condition will not apply if the densities trail off to zero, rather than experiencing a sharp cutoff at the boundaries of $\cC$. However, in such circumstances our results apply to a subset of the entire domain, chosen so that the density did not fall below a specified minimum. Performance inside this region would then conform to the optimal accuracies presented, while estimation outside this region would be poorer. This distinction is unavoidable, since estimation in the tails is usually problematic and it would be unusual to guarantee uniformly good performance there.

Step 1 of the LABAVS Algorithm allows the initial bandwidth to be chosen globally or locally. Here we shall focus on the global case, where an initial bandwidth $H = \diag (h^2,\ldots,h^2)$ is used. Further, we assume that this $h$ is chosen to minimise the mean integrated squared error (MISE):
$$ E \left[ \int \{\hg(x) - g(x) \}^2 f(x) dx \right] \,,$$
where the outer expectation runs over the estimator $\hg$. It is possible to show that under our assumptions that
\beqn  h = \left[ \frac{d\si^2 R(K) A_\cC}{n \mu_2(K)^2 A_{\cH_g}} \right] ^{-1/(d+4)} \,. \label{eq:swan} \neqn
Notice in particular that $h \asymp n^{1/(d+4)}$. Details are given in Lemma~\ref{le:ace} in the appendix (Miller and Hall, 2010).

A key result in establishing good performance, in Theorem~\ref{th:1} below, is  uniform consistency of the local polynomial parameter estimates. It is a simplified version of a result by Masry (1996), and no proof is included.

\begin{thm} Suppose the conditions in (\ref{eq:cond}) hold and we use parameter estimates from a degree $p$ polynomial regression to estimate the partial derivatives of $g$. Then for each $k$ with $0 \le |k| \le p$ we have
$$ \sup_{x\in \cC} | (\widehat{D^k g})(x) - (D^k g)(x)| = O \left[ \left( \frac{\log n}{n h^{d+2|k|}} \right)^{1/2} \right] + O(h^{p-|k|+1}) \mbox{ almost surely.}$$ \label{th:1} \end{thm}

Since the partial derivative estimate at $x$ is proportional to the corresponding local polynomial coefficient, Theorem~1 ensures that the local polynomial coefficients are consistently estimated uniformly for suitable $h$. The scaling applied in (\ref{eq:scale}) does not impact on this, as the proof of Theorem~\ref{th:appr1} demonstrates.

Let $\cC^-$ denote the points $x\in\cC$ satisfying $\pa g(x)/\pa x^{(j)} = 0$ and $\pa^2 g(x)/\pa x^{(j)2} \ne 0$. That is, $\cC^-$ denotes the points where the true set of relevant variables changes. Notice that in the illustrative example in Section~\ref{se:eg} we had $\cC^- = \{ x \,|\, x^{(1)}=0 \,, x^{(2)}=0 \}$. The smoothness assumed of $g$ implies that $\cC^-$ has Lebesgue measure 0. Let $\de>0$ and let $\cO_\de$ be an open set in $\cC^-$ such that
\beqn \inf_{x \in \cC^+ \setminus \cO_\de,\, j \in \cA^+(x)} \ga_j = \de \,. \label{eq:duck} \neqn
Intuitively this means that on the set $\cC \setminus \cO_\de$ the relevant variables have the absolute value of their corresponding parameters $|\ga_j|$ bounded below by $\de>0$, while irrelevant variables have $\ga_j=0$. Thus we have a ``gap'' between the true and irrelevant variables in this region that we may exploit. The volume of $\cO_\de$ may be made arbitrarily small by choosing $\de$ small. Call the set $\chA^+(x)$ in the algorithm correct if the variables in it are the same as the set of variables $j$ with $\partial g(x)/\partial x^{(j)} \neq 0$. Denote the latter correct set by $\cA^+(x)$.

\begin{thm} Suppose $\de$ is given, $h$ is chosen to minimise squared error as in (\ref{eq:swan}), $\chA^+(x)$ is formed using the first approach in Section~\ref{se:approaches}, and $\la$ has a growth rate between arbitrary constant multiples of $h^2 (n\log n)^{1/2}$ and $hn^{1/2}$. If $f$ has bounded and uniformly continuous derivatives of degree 2, then the probability that $\chA^+(x)$ is correct on $\cC \setminus \cO_\de$ tends to 1 as $n \to \infty$. Furthermore, variables that are genuinely redundant everywhere will be correctly classified as such with probability tending to 1.
\label{th:appr1}  \end{thm}

The property (\ref{eq:duck}) ensures that the coefficients in the local linear fit are consistently estimated with error of order $O\{h(\log n)^{1/2}\}$. The adjustment in (\ref{eq:scale}) means that the actual coefficients estimated are of order $hn^{1/2}$ times this, so the range of $\la$ given is correct for separating true and redundant variables. The definition of $\cO_\de$ ensures that the classification is correct on $\cC \setminus \cO_\de$, while variables that are redundant everywhere will be recognised as such.

The next result ensures consistency for the second approach in Section~\ref{se:approaches}. We make one further assumption, concerning the error $\ep_i$.  Observe that this holds trivially if $\ep_i$ is bounded. Assume that:
\beqn \begin{array}{p{14cm}}
there exists $C_3$ such that $|E(\ep_i^\al)| \le C_3^{\al-2}\si^2$ for $\al=3,4,\ldots$.
\end{array} \label{eq:extra} \neqn

\begin{thm} Suppose $\de$ is given, $h$ is chosen to minimise squared error as in (\ref{eq:swan}), and $\chA^+(x)$ is formed using the second approach in Section~\ref{se:approaches}. Provided that $\la = o(h^2)$ and $h^4 \log n = o(\la)$, the probability that $\chA^+(x)$ is correct on $\cC \setminus \cO_\de$ tends to 1 as $n \to \infty$. Furthermore, variables that are genuinely redundant everywhere will be correctly classified as such with probability tending to 1.
\label{th:appr2} \end{thm}

The previous two results ensure that we have consistent variable selection for the first two approaches in Section~\ref{se:approaches}. Finally we can state and prove the strong oracle property for $\cC \setminus \cO_\de$. Although the result does not cover the whole space $\cC$, recall that we may make the area $\cO_\de$ arbitrarily small by decreasing $\de$. Furthermore, the proof implies that if we restricted attention to removing only those variables that are redundant everywhere, we would actually have the oracle property on the whole of $\cC$; however we sacrifice this performance on $\cO_\de$ to improve the fit elsewhere by adjusting for locally redundant variables. In the following theorem the matrix $\ddot{H}$ is the diagonal bandwidth matrix with bandwidth $\infty$ for globally redundant variables and $\ddot{h}$ for the other variables, where
$$ \ddot{h} = h \left( M(\ddot{H},H) \ddot{d}/d \right)^{1/4} \,.$$
Here $\ddot{d}$ denotes the number of variables that are not globally redundant.

\begin{thm} The estimates produced by the algorithm, where variable selection is performed using the first or second approach in Section~\ref{se:approaches}, satisfy the strong definition of the nonparametric oracle property on $\cC$. Further,  when there are locally redundant variables, squared estimation error is actually less that the oracle performance by a factor of $M[\{H^L(X),H^U(X)\},\ddot{H}] <1$. \label{th:big} \end{thm}

\section{Numerical properties} \label{se:num}

The examples presented in this section compare the performance of two versions of the LABAVS algorithm with ordinary least squares, a traditional local linear fit,  generalised additive models, tree-based gradient boosting and MARS. Table~\ref{tab:types} describes the approaches used. The implementations of the latter four methods were from the R packages {\sl locfit}, {\sl gam}, {\sl gbm} and {\sl polspline} respectively. Tuning parameters such as bandwidths for local methods, $\la$ in LABAVS, number of trees in boosting, and MARS model complexity, were chosen to give best performance for each method. The LABAVS models used the first variable selection approach of Section~\ref{se:approaches}. All the local methods used nearest neighbour bandwidths. The OLS linear model was included as a standard benchmark, but obviously will fail to adequately detect nonlinear features of a dataset.

\begin{table}[htbp]\begin{center}
\caption{ Approaches included in computational comparisons  \label{tab:types}}
\begin{tabular}{llll}
\hline
 Name &  Description \\
\hline
LABAVS-A & LABAVS with linear fit, all vars in final fit \\
LABAVS-B & LABAVS with linear fit, relevant vars only in final fit \\
LOC1      & Local linear regression     \\
OLS       & Ordinary least squares linear regression \\
GBM       & Boosting with trees, depth equal to three \\
GAM       & Generalised additive models with splines  \\
MARS      & Multivariate adaptive regression splines  \\
\hline
\end{tabular} \end{center} \end{table}

\begin{figure}[htb] \begin{center}
\includegraphics[width = 3.4in, height=3.4in]{sim3heat.pdf}
\vspace{-8pt}
{\it \caption{ Plot of detected variable significance across subspace in Example 1. \label{fig3}}}
 \end{center} \end{figure}

\ni {\bf Example 1:} The example introduced in Section~\ref{se:eg} was simulated with $n=500$. The error for $Y_i$ was normal with standard deviation $0.3$. We first compare LABAVS to the rodeo and the methodology of Bertin and Lecu\'{e} (2008) at the four representative points in Figure~\ref{fig2}. Table~\ref{tabcomp1} shows the mean squared error of the prediction compared to the true value over 100 simulations. In all cases parameters were chosen to minimise this average error. At all points the LABAVS approach performed strongest. The method of Bertin and Lecu\'{e} (2008) performed poorly in situations where at least one variable is redundant; this is to be expected, since it excludes the variable completely and so will incorporate regions where it is actually important, causing significant bias. The rodeo also did not perform as well; we found it tended to overestimate the optimal bandwidths in redundant directions.

\begin{table}[htb]\begin{center}
\caption{ Mean squared prediction error on sample points in Example 1 \label{tabcomp1}}
\begin{tabular}{lcccc}
\hline
 Test Point & LABAVS-A & LABAVS-B & rodeo & Bertin and\\
&&&& Lecu\'{e} \\ \hline
(1,1)  & 0.0022 & 0.0022 & 0.0065 & 0.0023 \\
(1,-1) & 0.0011 & 0.0013 & 0.0015 & 0.0018 \\
(-1,1) & 0.0009 & 0.0011 & 0.0015 & 0.0013 \\
(-1,-1)& 0.0006 & 0.0007 & 0.0008 & 0.0013 \\
\hline
\end{tabular} \end{center} \end{table}

We then compared LABAVS with the other model approaches which are designed to make multiple predictions, rather than a specific point. For each simulation all the models were fitted and the average squared error was estimated using a separate test set of 500 observations. The simulation was run 100 times and the average error and its associated standard deviation for each model are recorded in Table~\ref{tab1}.

\begin{table}[htbp]\begin{center}
\caption{ Mean squared error sum of test dataset in Example 1 \label{tab1}}
\begin{tabular}{lcc}
\hline
\quad Approach &  Error & Std Dev\\
\hline
LABAVS-A  &    2.18   &       (0.71) \\
LABAVS-B  &    1.87   &       (0.65) \\
LOC1   &       2.31   &       (0.73) \\
OLS    &       42.85  &       (2.64) \\
GBM    &       2.47   &       (0.67) \\
GAM    &       5.93   &       (0.57) \\
MARS   &       2.35   &       (0.90) \\
\hline
\end{tabular} \end{center} \end{table}

Inspection of the results shows that the LABAVS models performed best, able to allow for the different dependencies on the variables. In particular the algorithm improved on the performance of the local linear model on which it is based. The local linear regression, the boosted model and MARS also performed reasonably, while the GAM struggled with the nonadditive nature of the problem, and a strict linear model is clearly unsuitable here.

To show how effectively variable selection is for  LABAVS, Figure~\ref{fig3} graphically represents the sets $\chA^+$ at each grid point for one of the simulations, with the darkest representing $\{\}$, the next darkest $\{1\}$, the next darkest $\{2\}$ and finally the lightest $\{1,2\}$. Here the variable selection has performed well; there is some encroachment of irrelevant variables into the wrong quadrants but the selection pattern is broadly correct. The encroachment is more prevalent near the boundaries since the bandwidths are slightly larger there, to cover the same number of neighbouring points.

\ni {\bf Example 2:} We next show that LABAVS can effectively remove redundant variables completely. Retain the setup of Example~1, except that we add $d^*=d-2$ variables similarly distributed (uniform on [-2,2]), which have no influence on the response. Also, keep the parameters relating to the LABAVS fit the same as the previous example, except that the cutoff for hard threshold variable selection, $\la$ is permitted to vary. Table~\ref{tabex2} shows the proportion of times from 500 simulations that LABAVS effected complete removal of the redundant dimensions, for various $\la$ and $p^*$. Note that the cutoff level of 0.55 is that used in the previous example, and the two genuine variables were never completely removed in any of the simulations. The results suggest that to properly exclude redundant variables, a higher threshold is needed than would otherwise be the case. This causes the final model to be slightly underfit when compared to the oracle model, but this effect is not too severe; Figure~\ref{figex2} shows how the variable significance plots change for a particular simulation with different values of the cutoff. It is clear that the patterns are still broadly correct, and the results still represent a significant improvement to traditional linear regression.

\begin{table}[htbp]\begin{center}
\caption{ Proportion of simulations where redundant variables completely removed by LABAVS \label{tabex2}}
\begin{tabular}{lcccc}
\hline
\quad  & \multicolumn{4}{c}{Number of redundant dimensions}  \\
\hline
$\la$ & 1&2&3&4 \\ \hline
0.55& 0.394&0.086&0.034&0.038\\
0.65& 0.800&0.542&0.456&0.506\\
0.75& 0.952&0.892&0.874&0.864\\
0.85& 0.996&0.984&0.994&0.974\\
0.95& 0.998&1.000&1.000&0.992\\
\hline
\end{tabular} \end{center} \end{table}

\begin{figure}[tbhp] \begin{center}
\includegraphics[width = 5in, height=1.7in]{sim4.pdf}
\vspace{-8pt}
{\it \caption{ Plot of detected variable significance across subspace in Example 2, under various choices for $\la$. \label{figex2}}}
\end{center} \end{figure}

\ni {\bf Example 3:} The first real data example used is the ozone dataset from Hastie~\etal (2001), p.175. It is the same as the air dataset in S-PLUS, up to a cube root transformation in the response. The dataset contains meteorological measurements for New York collected from May to September 1973. There are 111 observations in the dataset, a fairly moderate size. Our aim here is to predict the ozone concentration using two of the other variables, temperature and wind and scaled to unit variance when fitting the models. The smoothed perspective plot of the data in Figure~\ref{fig4} shows strong dependence on each of the two variables in some parts of the domain, but some sections appear flat in one or both directions in other parts. For example, the area surrounding a temperature of 70 and wind speed of 15 appears to be flat, implying that for reasonably low winds and high temperatures the ozone concentration is fairly stable. This suggests that LABAVS, by expanding the bandwidths here, could be potentially useful in reducing error. We performed a similar comparative analysis to that in Example 1, except that error rates were calculated using leave-one-out cross validation, where an estimate for each individual observations was made after using all other observations to build the model. The resulting mean squared errors and corresponding standard deviations are presented in Table~\ref{tab:Oz}. The results suggest that the data is best modelled using local linear methods, and that LABAVS offers a noticeable improvement over a traditional local fit, due to its ability to improve the estimate in the presence of redundant variables. The perspective plot in left panel of Figure~\ref{fig4} suggests a highly non-additive model, which may explain why GAM performs poorly. There is also a large amount of local curvature, which hinder the OLS, GBM and MARS fits. The right panel of Figure~\ref{fig4} shows the variable selection results for the linear version of LABAVS across the data support, using the same shading as in Figure~\ref{fig2}. We see that variable dependence is fairly complex, with all combinations of variables being significant in different regions. In particular, notice that the procedure has labelled both variables redundant in the region around $(70,15)$, confirming our initial suspicions. This plot is also highly suggestive, revealing further interesting features. For instance, there is also little dependence on wind when temperatures are relatively high. Such observations are noteworthy and potentially useful.

\begin{table}[htb]\begin{center}
\caption{ Cross-validated mean squared error sum for the ozone dataset  \label{tab:Oz}}
\begin{tabular}{lcc}
\hline
\quad Approach &  Error & Std Dev\\
\hline
LABAVS-A  &   277  &       (53) \\
LABAVS-B  &   284  &       (55) \\
LOC1      &   290  &       (55) \\
OLS       &   491  &       (110) \\
GBM       &   403  &       (118) \\
GAM       &   391  &       (98) \\
MARS      &   457  &       (115) \\
\hline
\end{tabular} \end{center} \end{table}

\begin{figure}[tbh] \begin{center}
\includegraphics[width =5in, height=2.25in]{Ozone3.jpg}
\vspace{-8pt}
{\it \caption{ Ozone dataset smoothed perspective plot and variable selection plot. \label{fig4}}}
 \end{center} \end{figure}

\ni {\bf Example 4:} As a second low-dimensional real data example, we use the ethanol dataset which has been studied extensively, for example by Loader (1999). The response is the amount of a certain set of pollutants emitted by an engine, with two predictors: the compression ratio of the engine and the equivalence ratio of air to petrol. There are 88 observations, a fairly moderate size. Inspection of the data shows strong dependence on the equivalence ratio, but the case for the compression ratio is less clear. This suggests LABAVS could be potentially useful in reducing error. We performed a similar analysis to that in Example 1, with the results are presented in Table~\ref{tab:Eth}.

\begin{table}[htb]\begin{center}
\caption{ Cross-validated mean squared error sum for the ethanol dataset  \label{tab:Eth}}
\begin{tabular}{lcc}
\hline
\quad Approach &  Error & Std Dev\\
\hline
LABAVS-A  &   0.075  &       (0.011) \\
LABAVS-B  &   0.085  &       (0.014) \\
LOC1      &   0.090  &       (0.012) \\
OLS       &   1.348  &       (0.128) \\
GBM       &   0.104  &       (0.020) \\
GAM       &   0.098  &       (0.012) \\
MARS      &   0.045  &       (0.008) \\
\hline
\end{tabular} \end{center} \end{table}

The results in Table~\ref{tab:Eth} show that this problem is particularly suited to MARS, which performed the best. After MARS, LABAVS produced the next strongest result, again improving on the traditional local linear model. The GBM and GAM models were inferior to the local linear fit.

config: arxiv


TITLE: Continuity of a function such that $\lambda f(x) = f(\lambda x) $
QUESTION [1 upvotes]: Given $f:\mathbb{R} \rightarrow \mathbb{R} $ such that $\lambda f(x) = f(\lambda x) \mbox{   }\forall \lambda, x \in \mathbb{R}$. Show that $f$ is continuous.
I know one shall prove that $f$ is continuous in all its domain. I tried to use the $\epsilon - \delta$ definition and prove that $\lambda f(x)$ is continuous, so for the case $\lambda=1$, $f$ is continuous. But I got stuck.

REPLY [0 votes]: I think I got it:
for a given $\epsilon > 0, \exists \delta = \frac{\epsilon}{|f(1)|+1}$ such that $|x-a|<\delta $ implies $|f(x) -f(a)|=|f(1\cdot x) -f(1 \cdot a)| = |f(1)\cdot x -f(1) \cdot a| = |f(1)| \cdot|x -a| \leq |f(1)|\cdot \frac{\epsilon}{|f(1)|+1} < \epsilon$.
Hence, $f$ is continuous.


TITLE: Inequality for recursive-defined values
QUESTION [1 upvotes]: $f_{n+2} = \frac{6}{5}f_{n+1}-f_{n}, f_0 = 0,  f_1 = 1$
I need to prove that $f_n < 5/4$
I found that $f_{n} = \frac{1}{8} i 5^{1-n} \left((3-4 i)^n-(3+4 i)^n\right)$ and spend much time for efforts to solve this equality, but end in failure..

REPLY [2 votes]: If we take $z=\frac{3-4i}{5}$ we have that:
$$ f_n = \frac{5i}{8}(z^n-\bar{z}^n) $$
but since $z$ belongs to the unit circle so does $z^n$ and $|z^n-\bar{z}^n|$ cannot exceed $2$.


TITLE: Prove $f$ is constant in $\Omega$
QUESTION [0 upvotes]: They ask me:

$\Omega$ is a domain and $f$ is holomorphic in $\Omega$, prove that if
  $\lvert f(z) \rvert$ is constant in $\Omega$ then $f$ is also constant
  in $\Omega$

My idea to prove this is use the maximum modulus principle like:

I say as $\lvert f(z) \rvert$ is constant then $\lvert f(z) \rvert = M$, $M$ arbitrary value
As $\lvert f(z) \rvert = M$ I noticed that $\lvert f(z) \rvert$ reach the maximum
So using the the maximum modulus principle I can affirm that $f$ is constant.

I don't know if it is correct or not.

REPLY [0 votes]: Your solution is correct. Note that this property is usually proved before proving the maximum modulus principle, using the Cauchy-Riemann equations:
$$|f(z)|=const\Rightarrow u^2(x,y)+v^2(x,y)=const$$
Taking the partial derivatives and using Cauchy-Riemann we get:
$$2uu_x+2vv_x=0\Rightarrow2uu_x-2vu_y=0$$
$$2uu_y+2vv_y=0\Rightarrow2uu_y+2vu_x=0$$
We get a homogeneous system of 2 equations with variables $u_x,u_y$. For a non-trivial solution (which naturally gives a constant function since $u_x=u_y=0$ and similarly $v_x=v_y=0$), we required the determinant to be $0$ and so:
$$4u^2+4v^2=0\Rightarrow u=v=0$$
Which once again gives a constant solution. Either way - $f$ is constant.


TITLE: How to solve for the flow of two layers of viscous sandwiched between two moving plates
QUESTION [1 upvotes]: I've come across a problem which modifies the traditional problem of Couette flow by adding two layers of viscous fluid with different dynamic viscosities $\mu_{1}$ and $\mu_{2}$. Suppose we have two infinite plates, distance $h$ apart, with the plate at $h = 0$ stationary and the plate at $h$ moving at velocity $U$. Suppose we have a layer of fluid $A$ which occupies the region $0 \leq y \leq h_{1}$, and fluid $B$ is between $h_{1} \leq y \leq h$, each with dynamic viscosities $\mu_{1}$ and $\mu_{2}$ respectively.
I want to solve for the velocity profile. Now it is clear to me that we can assume the velocities $\mathbf{u_{1}}$ and $\mathbf{u_{2}}$ are directed only in the $x$ direction. So we assume that:
$$\mathbf{u_{i}} = (u_{i}(y), 0, 0). $$
Now, plugging into Navier Stokes we get that $u_{i}''(y) = 0$. The obvious boundary conditions are that
\begin{align}
u_{1}(0) &= 0 \\
u_{2}(h_{1}) &= u_{1}(h_{1}) \\
u_{2}(h) &= U
\end{align}
And, clearly we have that $u_{i}(y) = A_{i} + B_{i}y$ where the constants are determined by these boundary conditions. However, upon substitution I found that these conditions aren't enough to determine the flow completely. My question is, do we have another boundary condition at the interface of the fluid regarding the shear stresses in between the two fluids being equal? If so, I suspect it looks like
$$ \mu_{1} \frac{\partial u_{1}}{\partial y} = - \mu_{2}\frac{\partial u_{2}}{\partial y}$$
at $y = h_{1}$.
Having done this, I've determined the system in full, for example finding that $$ u_{1} = y\frac{U\mu _{2}}{h_{1}(\mu_{1} + \mu_{2}) - (h-h_{1})\mu_{1}}.$$ Is this the right track in proceeding the problem?

REPLY [1 votes]: You have the right idea, but

you have a sign error at the boundary, as mentioned by Chester in a comment
you left $y$ out of your equation for $u_1$

You should find
$$u_1 = Uy\frac{ \mu_2 }{h_1 \mu_2 + \mu_1 (h-h_1)}$$
It should be possible to detect these sorts of errors with a little care. For example, suppose $h = h_1$. Then the second fluid doesn't exist and $\mu_2$ has to disappear from your equation. Does it with your original answer? Does it with mine?
Or if $h_1 =0$, only fluid 2 exists; does your solution for $\mu_2$ do the right thing?
Also, you can detect things like forgetting the $y$ by doing a quick dimensional check.


TITLE: How to negate the satement "If $3y$ is a perfect cube then there exists an even number $n$ so that $y = 3n^3$"
QUESTION [1 upvotes]: Let $x \in \mathbb{Z}$. If $3y$ is a perfect cube then there exists an even number $n$ so that $y = 3n^3$.

Question: I am a little bit confused here. I know this statement is false and I would like to find 2 things.
(1) What is the negation of this statement? Important: If I want to disprove a statement, do I need to take the negation of the WHOLE implication like $\neg (a \implies b)$? Or just take the negation by parts like $a \implies b \equiv \neg a \implies \neg b$?
I need to find a counter example, but according to a professor I need to find it using the negation.
(2) In what cases would you use $\neg(a \implies b)$ versus $\neg a \implies \neg b$?

REPLY [1 votes]: (1) What is the negation of this statement? Important: If I want to disprove a statement, do I need to take the negation of the WHOLE implication like $\neg (a \implies b)$? Or just take the negation by parts like $a \implies b \equiv \neg a \implies \neg b$?

This is important, yes! You ALWAYS take the negation of the whole statement, not the parts. The negation of an implication is NOT the negation of its parts. What we do is
$$
\lnot (a \implies b) = a \land (\lnot b).
$$
In the end you should get the statement: "$3y$ is a perfect cube and for all even numbers $n$, $y \ne 3n^3$."

(2) In what cases would you use $\neg(a \implies b)$ versus $\neg a \implies \neg b$?

Always the first one. Never the second, when negating a statement.
$\lnot a \implies \lnot b$ may come up if you are doing a proof and you want to prove the converse of a statement $a \implies b$. But don't worry about that.


TITLE: De Rham cohomology and antiderivatives
QUESTION [3 upvotes]: A couple of recent questions about antiderivatives reminded me of the following, which I can't recall seeing tackled explicitly anywhere to my satisfaction and that I sketched to an ambitious calculus student long ago. I hope it is not completely trivial...
Consider $f: \mathbb{R} \rightarrow \mathbb{R}$ nice. An antiderivative $F$ of $f$ is defined by $f \ dx = dF$. Let $\Omega = [a,b]$. The Stokes theorem is basically that $d$ and $\partial$ are adjoint w/r/t the standard pairing (a/k/a the integral), viz. $\langle dF,\Omega \rangle = \langle F, \partial \Omega \rangle$.
This is obviously just an overly fancy formulation of the fundamental theorem of calculus (à la Spivak, for example). But its generality is a virtue in motivating the following question (which also betrays my long-decayed knowledge of the most trivial bits of algebraic topology): is there a general formulation of the notion of antiderivative that incorporates basic information about de Rham (co)homology? As a followup, is it necessary to bring bundles into the picture?

REPLY [2 votes]: If we assume that the higher dimensional generalization of the fundamental theorem of calculus,
$$ \int_a^b f'(x)\,dx = f(b) - f(a)$$
is Stokes' theorem,
$$ \int_{M} d\omega = \int_{\partial M} \omega, $$
then it appears to me that the higher dimensional generalization of an "antiderivative" is the "anti-exterior-derivative" of an exact $n$-form and the higher dimensional analogue of an "arbitrary constant" is a closed $(n-1)$-form.
So it appears that the $(n-1)$-th deRham cohomology group of $M$ should say something about the space of "arbitrary constants", but I'm not sure what.


TITLE: Square root of an integer has only even digits
QUESTION [12 upvotes]: Is there a non-square positive integer $n$, that $\sqrt{n}$ has only even digits in its decimal representation ?

REPLY [3 votes]: No such $n$ is known.
If one were found, it would be the biggest shock in Mathematics since, well, maybe since ever; certainly, since Godel's incompleteness results.
No proof is known that such an $n$ does not exist.


TITLE: Finding solutions to a non homogeneous differential equation knowing the solutions of its homogenous counterpart
QUESTION [2 upvotes]: Let $f(x)$ and $xf(x)$ be the particular solutions of a differential
  equation $$y''+R(x)y'+S(x)y=0$$ Find the solution of the differential
  equation $$y''+R(x)y'+S(x)y=f(x)$$ in terms of $f(x)$.

I was trying to find a solution if the form $(ax^2+bx+c)f(x)$. Now as $(bx+c)f(x)$ is a solution to the homogeneous equation, I just need to assume that $ax^2f(x)$ is a particular solution of the non homogeneous equation, and find the suitable value of $a$. But I am getting stuck. Please help.

REPLY [2 votes]: If $y(x) = u(x) f(x)$ is a solution of the non-homogeneous DE, substituting that in should give you
$$ u''(x) f(x) + u'(x) (2 f'(x) + R(x) f(x)) + u(x) (f''(x) + R(x) f'(x) + S(x) f(x)) = f(x) $$
From the fact that $y=f(x)$ and $y=x f(x)$ are solutions, you should get that the coefficients of $u(x)$ and $u'(x)$ are $0$.


TITLE: Velocity distribution for Molecular beam vs. gas
QUESTION [2 upvotes]: In the book: https://g.co/kgs/JRBLeA (Atomic Physics by C.J. Foot) on page 152 there's a table in the end that compares gas' velocity distribution and that for a molecular beam. I am unable to account for the extra velocity/speed factor in the latter.

I understand that the former has a Maxwell Boltzmann distribution.
The text mention, on the same page, that this extra factor is attributed to the observation that faster atoms a more likely to effuse through a hole of area A compared to slower atoms. I don't understand why this should be so. First of all, if you have an atom beam (say) traveling with a particular velocity, how do you associate the concept of temperature to this?

REPLY [2 votes]: You probably know how to derive an expression for the pressure in the context of the kinetic gas theory.  If we use this line of argument, but we were interested only in the collisions between the molecules and the small area (=hole) within the wall, we have to ask: If a particle exists the hole during the time interval $[0, dt]$, where does it come from?
Assuming that the particle has velocity $v$, and that no collision takes place in the time interval $[0, dt]$ the particles must come from the volume which is ${\bf ds} = {\bf v}\cdot dt$ away from the hole. Or, more precisely, they must come from a volume $V = {\bf A} \cdot {\bf ds} = A \, v dt \cdot \cos{\theta}$, where $\theta$ is the angle of the veclocity to the normal direction of the hole. The key point here is that particles with a "larger" velocity are allowed to be further way at time $t=0s$. Thus, "very fast" particles have a non-zero probability to exit the hole, no matter where they are at time $t=0s$, while "slow" particles must come from the direct surrounding of the hole. Proper derivations can be found in statistical physics textbooks.

[...] if you have an atom beam [...] traveling with a particular velocity, how do you associate the concept of temperature to this?

If we take about a gas having temperature $T$ we implicitly assume a random velocity distribution, where the atoms possess different velocities. Also, temperature uses the concept of an ensemble of particles -- we are unable to associate a temperature with a single particle, which happens to have the kinetic energy $k_B T$. So, if we use a velocity filter to generate a beam of particles with the velocity range $v + [0, dv]$, we should not use the term temperature to describe it.
However, it is conceptional fine to use the small hole in the wall of the container to generate a beam of particles, which is used probe the distribution of the ensemble inside the container.


TITLE: Uniform lamina (equilateral triangle) in a hoop.
QUESTION [0 upvotes]: A uniform lamina is in the shape of an equilateral triangle $ABC$, and rests in equilibrium inside a vertical circular hoop with $B$ and $C$ in contact with the hoop, and $A$ at the centre of the circle. Contact is smooth at $C$ and rough at $B$, where the coefficient of friction is $μ$. The lamina is in limiting equilibrium with $B$ lower than $C$, and the axis of symmetry through $A$ makes an angle with the vertical. By first taking moments about $A$, show that
$$tanθ=\frac{μ√3}{√3-2μ}$$
The lamina $ABC$ is placed so that $θ = 0$, and the hoop is slowly rotated about its centre in a vertical plane so that $θ$ increases. Find the value, $θ_o$, of $θ$ beyond which the lamina will topple about $B$, if it has not slid already. Also find the least value of $μ$ for which $θ$ will reach the value of $θ_o$ without the lamina sliding.
The axis of symmetry, which is from $A$ to the base of the triangle (which forms a chord of the circle (hoop)) has length $\frac{√3r}{2}$. The normal reaction is in line with $A$ and the frictional force is at right angles to this at $B$. The lamina will have weight $W$.

REPLY [0 votes]: I have modeled it slightly differently.
I have a normal reaction $R$ acting along $AB$ towards $A$. The frictional force $μR$ orthogonal to this pointing to the left and $W$ vertically downwards from the centre of $BC$.
So I have:
Moments about $A$:  $$\frac{r√3}{2}W=rμRsinθ$$
So $$W=\frac{2rμRsinθ}{r√3}$$
Resolving vertically (with R pointing towards the right and μR to the left (at an angle $θ$)): $$W=Rsinθ - μRcosθ$$
So $$\frac{2rμRsinθ}{r√3}= Rsinθ - μRcosθ$$
Which gives $$2rμsinθ - r√3sinθ = -r√3μcosθ$$
Which rearranges to give $$tanθ = \frac{μ√3}{√3-2μ}$$
For the value of $μ$, I have:
$$tan(60)=\frac{μ√3}{√3-2μ}$$
$$√3(√3-2μ)=μ√3$$
$$μ=\frac{3}{√3+2√3}=\frac{√3}{3}$$


TITLE: Compute this particular determinant
QUESTION [2 upvotes]: I want to know if there is a good way to compute the following determinant
$$D_n= \begin{vmatrix}a_1b_1& a_1b_2 & a_1b_3&\cdots & a_1b_n \\a_1b_2 & a_2b_2 & a_2b_3 &\cdots& a_2b_n \\ a_1b_3 &  a_2b_3 &  a_3b_3&\cdots& a_3b_n \\\vdots&\vdots&\vdots&\ddots&\vdots\\a_1b_n & a_2b_n & a_3b_n &\cdots&a_nb_n \end{vmatrix} $$
It's not hard to show that $D_n$ can be reduced to the following
$$D_n=a_1 \begin{vmatrix}b_1& a_1b_2-a_2b_1 & a_1b_3-a_3b_1&\cdots & a_1b_n-b_1a_n \\b_2 & 0 & a_2b_3-a_3b_2 &\cdots& a_2b_n-b_2a_n \\ b_3 &  0 &  0&\cdots& a_3b_n-b_3a_n \\\vdots&\vdots&\vdots&\ddots&\vdots\\b_n & 0 & 0&\cdots&0\end{vmatrix}. $$
Then is there a good way to compute this determinant efficiently?
Also, if $n=2$, the corresponding matrix can be written as multiplication by
two simpler matrix. But I don't know if one can do similar things in higher dimensions

REPLY [2 votes]: If you expand your determinant with respect to the last row, you get:
$$D=(-1)^{n+1}a_1b_n\prod_{i=1}^{n-1}(a_ib_{i+1}-a_{i+1}b_i)$$


TITLE: Expected Value Problem including permutations
QUESTION [0 upvotes]: Please help!
Let F denote some permutation that maps from A to A. Where A={1,2,...2n} is a set. For such a permutation, F, let p denote the number of indices i, belonging to A, such that $F(i) > 2i$
Suppose that the permutation F is chosen randomly, from the set of all permutations of the set F, meaning that each permutation is equally liked to be chosen.
How would I go about calculating the expected value of p. I have no idea even how to go about it.

REPLY [2 votes]: The probability that $F(i)>2i$ is $E(i)=\frac{2n-2i}{2n}$.
Then the expected value of $p$ is the sum of the $E(i)$ for $1\le i \le n-1$.
This is $\frac{1+2+ ... +n-1}{n}=\frac{n-1}{2}.$


TITLE: You are stopping by Timmy’s to buy $12$ donuts. There are $4$ varieties to choose from.
QUESTION [1 upvotes]: You are stopping by Timmy’s to buy $12$ donuts. There are $4$ varieties to choose from. You need to
have at least one donut from each variety. How many different ways can you select your $12$ donuts?

REPLY [2 votes]: Taking $1$ donut firstly from each variety to fulfill The condition for at least one donut from each variety, so we have $4$ donuts selected. $12-4 = 8$ donuts remain to be selected.

Now we can select $8$ donuts with any number of donuts from each variety(including zero selections), consider the $4$ varieties to be named $a,b,c,d$ then we have the situation as $a+b+c+d = 8$. Here $r = 4$ (four varieties) and $n = 8 $

So, to find the number of ways($w$) in this situation use $w = \binom{n+r-1}{r-1}$. applying this, we get the result as $\binom{11}{3}$ which is equal to $165$.


TITLE: Munkres, §22, Exercise 4 (a). Proof verification.
QUESTION [2 upvotes]: Define an equivalence relation on the plane $X = \mathbb{R}^2$ as follows:
\begin{align*}
x_0 \times y_0 \sim x_1 \times y_1 \quad \mbox{if} \quad x_0 + y_0^2 = x_1 + y_1^2.
\end{align*}
Let $X^*$ be the corresponding quotient space. It is homeomorphic to a familiar space; what is it? [Hint: Set $g(x \times y) = x + y^2$.]
I want to know if the following resolution is ok.
Let $g : \mathbb{R}^2 \to \mathbb{R}$ such that $g(x \times y) = x + y^2$. The map $g$ is continuous, because the addition and product are continuous functions. Furthermore $g$ is a surjective map, because for every $x \in \mathbb{R}$, there is $x \times 0 \in \mathbb{R}^2$ such that $g(x \times 0) = x$. Now, $X^*$ is the collection of the subsets of $X = \mathbb{R}^2$ consisting of all the points $x \times y$ such that $z = x + y^2$, for every $z \in Z = \mathbb{R}$. Then
\begin{align*}
X^* = \{ g^{-1}(\{ z \}) \; | \; z \in \mathbb{R} \}.
\end{align*}
We prove that $g$ is a quotient map. Consider all the points $x \times 0$ of the plane $\mathbb{R} \times \mathbb{R}$; these points belong to the line $\mathbb{R}$ (First coordinate of $\mathbb{R}^2$), taken as a subset of $X = \mathbb{R}^2$. Since $g(x \times 0) = x$, we can see that $g$ maps every point $x$ belonging to the line $\mathbb{R}$ to the same point $x$, and therefore $g$ is a retraction of $X = \mathbb{R}^2$ onto $\mathbb{R}$. And then, by exercise 2 (b), the function $g$ is a quotient map. Hence, by Corollary 22.3, pag. 142, of the Munkres book, it follows that there exists an homeomorphism $f : X^* \to \mathbb{R}$; and the familiar space is $\mathbb{R}$.

REPLY [1 votes]: You cannot have a retraction from $\mathbb{R}^2$ onto $\mathbb{R}$, because $\mathbb{R}$ is not a subset of $\mathbb{R}^2$.
But $\mathbb{R} \times \{0\}$ is, by $g'(x,y) = (g(x),0)$.
So $g'$ is a quotient map (using exercise 2(b) if you like), and as $(x,0) \rightarrow x$ (the projection) is a homeomorphism between $\mathbb{R} \times \{0\}$ and $\mathbb{R}$, the composition of these maps, which is $g$, is also quotient.
Now you apply the Corollary.


text:
\begin{document}

\title[Burnside's formula]{Sharp inequalities related with Burnside's formula}
\author{necdet batir}
\address{department of  mathematics, nev{\c{s}}ehir hbv university, nev{\c{s}}ehir, 50300 turkey}
\email{nbatir@hotmail.com}

\subjclass{33B15, 26D07}
\keywords{Burnside's formula, Stirling's formula, factorial function, inequalities}

\begin{abstract} We prove the following double inequality related with Burnside's formula for $n!$
\begin{equation*}
\sqrt{2\pi}\left(\frac{n+a_*}{e}\right)^{n+a_*}<n!<\sqrt{2\pi}\left(\frac{n+a^*}{e}\right)^{n+a^*}\,(n\in\mathbb{N}),
\end{equation*}
where the constants $a_*=0.428844044...$ and $a^*=0.5$ are the best possible. We believe that the method we used in the proof gives insight to undergraduate students to understand  how simple inequalities can be established.
\end{abstract}
\maketitle

\section{introduction}
In a paper written in 1733 the French mathematician Abraham de Moivre developped the formula
\begin{equation*}
n!\sim C\cdot\sqrt{n}n^ne^{-n},
\end{equation*}
where $C$ is a constant. He was unable, however, to evaluate numerical value of this constant; this task befell a Scot mathematician James Stirling (1692-1770), who found  $C=\sqrt{2\pi}$. Formula
\begin{equation}\label{e:1}
n!\sim n^ne^{-n}\sqrt{2\pi n}
\end{equation}
is  known as Stirling's formula today. It is known that this formula has many applications in statistical physics, probability theory and number theory. The most well known approximation formula for factorial function after Stirling formula is Burnside's formula \cite{2}, which is given by
\begin{equation}\label{e:2}
n!\sim\sqrt{2\pi}\left(\frac{n+\frac{1}{2}}{e}\right)^{n+\frac{1}{2}}.
\end{equation}
It is known that Burnside's formula is more accurate than Stirling's formula. The gamma function $\Gamma$ is defined by the improper integral, for $x>0$:
\begin{equation*}
\Gamma(x)=\int\limits_0^\infty t^{x-1}e^{-t}dt.
\end{equation*}
An important function related to $\Gamma$ is the digamma function $\psi$, which is defined by $\psi(x)=\Gamma'(x)/\Gamma(x)$ for $x>0$. The gamma function $\Gamma$ and factorial function $n!$ are related with  $\Gamma(n+1)=n!$, for all $n\in\mathbb{N}$. In the literature there are many inequalities related to Burnside's formula, for example
\begin{equation}\label{e:3}
\sqrt{2e}\left(\frac{n+\frac{1}{2}}{e}\right)^{n+\frac{1}{2}}<n!<e^{(1+\gamma)e^{-\gamma}}\left(\frac{n+e^{-\gamma}}{e}\right)^{n+e^{-\gamma}}.
\end{equation}
see \cite[Theorem 1.4]{1}, and
\begin{equation}\label{e:4}
\sqrt{2e}\left(\frac{n+\frac{1}{2}}{e}\right)^{n+\frac{1}{2}}<n!<\sqrt{2\pi}\left(\frac{n+\frac{1}{2}}{e}\right)^{n+\frac{1}{2}},
\end{equation}
see  \cite[Theorem 1.5]{1}. Lu \cite{4}  provided an asymptotic expansion for the gamma function starting from Burnside's formula. Chen \cite{3} found the following asymptotic expansion starting from Burnside's formula
\begin{equation*}
n!\sim\sqrt{2\pi}\left(\frac{n+\frac{1}{2}-\frac{1}{24}\frac{1}{(n+1/2)^2}+\frac{19}{5760}\frac{1}{(n+1/2)^3}+\cdots}{e}\right)^{n+\frac{1}{2}}.
\end{equation*}
In this short note, motivated by the inequalities (\ref{e:3}) and (\ref{e:4}), we determine the largest number $\alpha$ and the smallest number $\beta$ in such away that the following inequalities hold for all $n\in\mathbb{N}$:
 \begin{equation*}
\sqrt{2\pi}\left(\frac{n+\alpha}{e}\right)^{n+\alpha}<n!<\sqrt{2\pi}\left(\frac{n+\beta}{e}\right)^{n+\beta}.
\end{equation*}
We prove that the best possible constants $\alpha$ and $\beta$ satisfying these inequalities are $\alpha=0.428844044...$ and $\beta=1/2$. The right hand side inequality here is already known, but we want to emphasize that we have shown here that the scaler $a^*=1/2$ here can not be replaced by a smaller quantity. We believe that the methods we used in the proof help undergraduate students to gain insight in establishing simple inequalities.

\section{main result}
Our main result is the following theorem.
\begin{thm} For all $n\in\mathbb{N}$ we have
\begin{equation}\label{e:5}
\sqrt{2\pi}\left(\frac{n+a_*}{e}\right)^{n+a_*}<n!<\sqrt{2\pi}\left(\frac{n+a^*}{e}\right)^{n+a^*},
\end{equation}
where the constants $a_*=0.428844044...$ and $a^*=0.5$ are the best possible.
\end{thm}
\begin{proof} The right hand side of (\ref{e:5}) is already known (see \cite{1}). We assume that the right side of (\ref{e:5}) holds. Then we have for all $n\in\mathbb{N}$:
$$
\sqrt{2\pi}\left(\frac{a_*+n}{e}\right)^{a_*+n}\geq n!
$$
or taking logarithms of both sides, a simple calculation gives
$$
\log \frac{n!}{\sqrt{2\pi n}n^ne^{-n}}+(n+1/2)\log n+a^*-(n+a^*)\log(n+a^*)<0.
$$
Letting $n\to\infty$ and applying Stirling's formula yields
\[
\lim\limits_{n\to\infty}\log\frac{n^{\frac{1}{2}-a^*}e^{a^*}}{\left(1+\frac{a^*}{n}\right)^{n+a^*}}<0 \quad \mbox{or}\quad \lim\limits_{n\to\infty}\left(\frac{1}{2}-a^*\right)\log n<0,
\]
which implies that $a^*\geq \frac{1}{2}$. Now, we assume that  the left side inequality of (\ref{e:5}) is valid for all $n=1,2,,\cdots$. Then if we substitute  $n=1$ there we obtain
\begin{equation*}
\sqrt{2\pi}\left(\frac{a_*+1}{e}\right)^{a_*+1}\leq1
\end{equation*}
or
\begin{equation}\label{e:6}
\frac{1}{2}\log (2\pi)+(a_*+1)\log(a_*+1)-a_*-1\leq0.
\end{equation}
Let us define  for $t\geq0$
$$
g(t)=(t+1)\log(t+1)+\frac{1}{2}\log(2\pi)-t-1.
$$
Then clearly, (\ref{e:6}) is equivalent to
\begin{equation}\label{e:7}
g(a_*)\leq0.
\end{equation}
Since
$$
g(0)=\frac{1}{2}\log (2\pi)-1=-0.0810615...<0
$$
and
$$
g(1)=2\log2+\frac{1}{2}\log (2\pi)-2=0.305233...>0,
$$
and $g$ is strictly increasing on $[0,\infty)$, we conclude that $g$ has only one real root on $(0,1)$, which is $t_0=0.428844\cdots$. We therefore obtain from (\ref{e:7}) $g(a_*)\leq0=g(t_0)$, which implies $a_*\leq t_0$ from fact that $g$ is strictly increasing.

Now we shall show that the left hand side of (\ref{e:5}) holds for $a_*=0.428844044...$. For this reason we define
$$
G_b(t)=\frac{1}{2}\log(2\pi)+(t+b)\log(t+b)-t-b-\log\Gamma(t+1), \quad t\geq 1.
$$
Differentiation gives
$$
G_b'(t)=\log(t+b)-\psi(t+1),
$$
where $\psi$ is the digamma function. By \cite[Lemma 1.7]{1}  $G_{a_*}$ is strictly decreasing. Thus we get for $t\geq 1$
$$
G_{a_*}(t)\leq G_{a_*}(1)=-6.58087\times10^{-11}<0,
$$
which gives the left hand side of (\ref{e:5}).

We conclude that the constants $a_*=0.428844044...$ and $a^*=0.5$ are the best possible constants.
\end{proof}

\bibliographystyle{amsplain}


TITLE: Riesz transform does not preserve continuity
QUESTION [1 upvotes]: I've read somewhere that the Riesz operator $R_j$ defined by $$R_j f(t) := c(n) \, \text{pv} \int_{\mathbb{R}^n} \frac{x_j}{|x|^{n+1}} f(t-x) \, dx$$doesn't preserve the continuity, but I can't produce a counterexample (actually not even when $n=1$). My idea was to take a continuous function $f$ (maybe with compact support) such that its transformed explodes at $0$, but I couldn't find such $f$. Do you have any hint on how to construct such function?
Thanks in advance!

REPLY [1 votes]: Let $g:[0,\infty) \to (-\infty,0]$ be continuous and compactly supported, with $g(0) = 0$ and $g(x) = 1/\ln x, 0< x < 1/2.$ On $\mathbb R^n,$ define
$$f(x) = (x_1/|x|)g(|x|), x \ne 0, \,f(0) = 0.$$
Then $f$ is continuous with compact support on $\mathbb R^n.$ Let's get at $R_1f(0):$
$$\int_{\epsilon<|x|<\infty} \frac{x_1}{|x|^{n+1}}f(-x)\,dx =  \int_{\epsilon<|x|<\infty} \frac{-x_1^2}{|x|^{n+2}}g(x)\,dx \ge \int_{\epsilon<|x|<1/2} \frac{-x_1^2}{|x|^{n+2}\ln |x|}\,dx. $$
Let' reduce the domain of integration to $U_\epsilon = \{\epsilon <|x|<1/2\} \cap \{x_1>|x|/2\}.$ The last integral above is then greater than
$$\tag 1\int_{U_\epsilon}(1/4) \frac{-1}{|x|^{n}\ln |x|}\,dx.$$
Note that $U_\epsilon$ can be written as the set of points $r\zeta,$ where $0<r<1/2$ and $\zeta \in \kappa,$ where $\kappa$ is an open "spherical cap" on the unit sphere centered at $(1,0,\dots,0).$ Using the polar coordinates formula* in $\mathbb R^n,$ we see the integral in $(1)$ equals
$$\tag 2 \sigma(\kappa)\int_\epsilon^{1/2} \frac{dr}{r\ln r},$$
where $\sigma$ is $(n-1)$-dimensional area measure on the unit sphere. The limit of $(2)$ is $\infty$ as $\epsilon \to 0^+$ and thus we've shown $R_1f(0) = \infty.$
*This formula can be found in Rudin's Real and Complex Analysis.