Smoothing Spline Semiparametric Density Models

TL;DR

This paper introduces a unified framework for semiparametric density estimation using reproducing kernel Hilbert spaces, providing new models, estimation methods, and theoretical guarantees validated through simulations and real data.

Contribution

It develops a systematic, kernel-based approach to semiparametric density modeling, estimation, and theory, unifying existing methods and extending their applicability.

Findings

Proposed general semiparametric density models encompassing existing ones.

Established joint consistency and convergence rates for estimators.

Validated methods through simulations and real data application.

Abstract

Density estimation plays a fundamental role in many areas of statistics and machine learning. Parametric, nonparametric and semiparametric density estimation methods have been proposed in the literature. Semiparametric density models are flexible in incorporating domain knowledge and uncertainty regarding the shape of the density function. Existing literature on semiparametric density models is scattered and lacks a systematic framework. In this paper, we consider a unified framework based on the reproducing kernel Hilbert space for modeling, estimation, computation and theory. We propose general semiparametric density models for both a single sample and multiple samples which include many existing semiparametric density models as special cases. We develop penalized likelihood based estimation methods and computational methods under different situations. We establish joint consistency and derive convergence rates of the proposed estimators for both the finite dimensional Euclidean parameters and an infinite-dimensional functional parameter. We validate our estimation methods empirically through simulations and an application.

Tables

Table 1. Table 1: Kullback-Leibler distances of different methods. Numbers inside parentheses of KLs are bias and variances. All numbers are in 10 − 2 superscript 10 2 10^{-2} .

$a$	Method	KL	KL	KL
		$n=100$	$n=200$	$n=500$
0.25	Kernel	2.56(0.04, 2.52)	1.65(0.03, 1.62)	0.96(0.03, 0.93)
	HG	2.48(0.47, 2.01)	1.21(0.18, 1.03)	0.62(0.12, 0.50)
	Cubic	1.53(0.31, 1.22)	0.81(0.21, 0.60)	0.41(0.11, 0.30)
	Quintic	1.08(0.04, 1.04)	0.43(0.01, 0.42 )	0.20(0.00, 0.20)
	SEMI	1.08(0.04, 1.04)	0.43(0.01, 0.42)	0.20(0.00, 0.20)
1	Kernel	2.82(0.04, 2.78)	1.84(0.07, 1.77)	1.25(0.04, 1.21)
	HG	2.30(0.14, 2.16)	1.38(0.43, 0.95)	0.67(0.11, 0.56)
	Cubic	1.51(0.28, 1.23)	0.80(0.15, 0.65)	0.45(0.12, 0.33)
	Quintic	0.93(0.03, 0.90)	0.46(0.02, 0.44)	0.22(0.01, 0.21)
	SEMI	0.93(0.03, 0.90)	0.47(0.02, 0.45)	0.24(0.01, 0.23)
4	Kernel	7.79(0.83, 6.96)	6.24(0.74, 5.50)	4.78(0.67, 4.11)
	HG	4.32(2.19, 2.13)	2.77(1.73, 1.04)	1.78(1.29, 0.49)
	Cubic	1.54(0.27, 1.27)	0.79(0.18, 0.61)	0.34(0.07, 0.27)
	Quintic	1.27(0.19, 1.08)	0.63(0.14 ,0.49)	0.33(0.15, 0.18)
	SEMI	1.27(0.19, 1.08)	0.64(0.15, 0.49)	0.35(0.16, 0.19)

Table 2. Table 2: Kullback-Leibler distances of different methods. Numbers inside parentheses of KLs are bias and variances. All numbers are in 10 − 2 superscript 10 2 10^{-2} .

$a$	Method	KL	KL	KL
		$n=100$	$n=200$	$n=500$
0.25	Kernel	3.50(0.41, 3.09)	2.59(0.36, 2.23)	1.72(0.21, 1.51)
	HG	3.20(1.12, 2.08)	1.89(0.68, 1.21)	1.10(0.49, 0.61)
	Cubic	2.12(0.59, 1.53)	1.32(0.42, 0.90)	0.63(0.19, 0.44)
	SEMI	1.58(0.16, 1.42)	0.64(0.01, 0.63)	0.45(0.04, 0.41)
1	Kernel	4.32(0.47, 3.85)	3.22(0.49, 2.73)	2.35(0.31, 2.04)
	HG	3.07(1.00, 2.07)	2.04(1.02, 1.02)	1.19(0.61, 0.58)
	Cubic	2.17(0.69, 1.48)	1.21(0.52, 0.69)	0.64(0.21, 0.43)
	SEMI	1.28(0.11, 1.17)	0.64(0.06, 0.58)	0.45(0.05, 0.40)
4	Kernel	11.47(1.86, 9.61)	9.68(1.88, 7.80)	7.68(1.53, 6.15)
	HG	7.83(5.85, 1.98)	5.27(4.08, 1.19)	3.50(3.04, 0.46)
	Cubic	2.00(0.65, 1.35)	1.08(0.37, 0.71)	0.52(0.19, 0.33)
	SEMI	2.33(0.58, 1.75)	1.30(0.32, 0.98)	0.66(0.13, 0.53)

Table 3. Table 3: Kullback-Leibler distances and mean squared errors of different methods. Numbers inside parentheses of KLs are bias and variances. Numbers inside parentheses of MSEs are squared biases and variances. All numbers are in 10 − 2 superscript 10 2 10^{-2} .

$\gamma$	Method	KL	MSE ($\hat{\gamma }$)	KL	MSE ($\hat{\gamma }$)
		$n=100$		$n=200$
1	Cubic	0.65(0.01, 0.64)		0.43(0.01, 0.42)
	TransKernel	1.87(0.69, 1.18)	13.40(11.44, 1.96)	1.57(0.74, 0.83)	15.95(15.05, 0.90)
	SEMI	1.09(0.02, 1.07)	3.07(0.01, 3.06)	0.61(0.00, 0.61)	1.33(0.01, 1.32)
2	Cubic	1.90(0.67, 1.23)		1.26(0.43, 0.83)
	TransKernel	2.05(0.52, 1.53)	58.35(50.44, 7.91)	1.64(1.04, 0.60)	61.80(58.02, 3.78)
	SEMI	0.94(0.02, 0.92)	8.63(0.17, 8.46)	0.62(0.01, 0.61)	6.02(0.18, 5.84)
3	Cubic	1.72(0.55, 1.17)		0.96(0.35, 0.61)
	TransKernel	1.92(0.74, 1.18)	123.05(106.61, 16.44)	1.37(0.71, 0.66)	135.24(128.98, 6.26)
	SEMI	0.95(0.03, 0.92)	21.07(0.48, 20.59)	0.48(0.01, 0.47)	10.75(0.34, 10.41)

Table 4. Table 4: Kullback-Leibler distances and mean squared errors of different methods. Numbers inside parentheses of KLs are bias and variances. Numbers inside parentheses of MSEs are squared biases and variances. All numbers are in 10 − 2 superscript 10 2 10^{-2} .

$n_1$	$n_2$	Method	Overall KL	${\text{KL}}_s$	MSE($\hat{\mu }$)	MSE($\hat{\sigma }$)
100	100	SEMI	5.36(0.92, 4.44)	2.56(0.49, 2.07)	3.58(0.00, 3.58)	2.34(0.00, 2.34)
		CHAR	13.39(10.06, 3.33)	6.86(3.53, 3.33)	2.24(0.00, 2.24)	1.95(0.00, 1.95)
		TP	5.81(1.95, 3.86)
	200	SEMI	4.27(0.79, 3.48)	2.56(0.41, 2.15)	2.31(0.01, 2.30)	1.41(0.00, 1.41)
		CHAR	11.07(8.82, 2.25)	6.26(4.02, 2.24)	1.58(0.00, 1.58)	1.34(0.00, 1.34)
		TP	4.95(1.61, 3.34)
200	100	SEMI	3.84(0.87, 2.97)	1.60(0.41, 1.19)	2.63(0.11, 2.52)	1.29(0.07, 1.22)
		CHAR	10.67(8.83, 1.84)	4.69(2.86, 1.83)	1.52(0.04, 1.48)	0.97(0.03, 0.94)
		TP	4.58(1.57, 3.01)
	200	SEMI	3.21(0.85, 2.36)	1.46(0.32, 1.14)	1.78(0.00, 1.78)	0.86(0.00, 0.86)
		CHAR	8.34(6.82, 1.52)	4.07(2.55, 1.52)	1.65(0.00, 1.65)	0.74(0.01, 0.73)
		TP	3.60(1.16, 2.44)

Table 5. Table 5: Kullback-Leibler distances and mean squared errors of different methods. Numbers inside parentheses of KLs are bias and variances. Numbers inside parentheses of MSEs are squared biases and variances. All numbers are in 10 − 2 superscript 10 2 10^{-2} .

$n_1$	$n_2$	Method	Overall KL	${\text{KL}}_s$	MSE($\hat{\mu }$)	MSE($\hat{\sigma }$)
100	100	SEMI	3.19(0.29, 2.90)	1.49(0.13, 1.36)	7.25(0.10, 7.15)	1.58(0.00, 1.58)
		CHAR	10.52(7.48, 3.04)	5.24(2.20, 3.04)	6.03(0.09, 5.94)	1.41(0.00, 1.41)
		TP	3.71(0.59, 3.12)
	200	SEMI	2.41(0.15, 2.26)	1.53(0.04, 1.49)	4.67(0.15, 4.52)	1.44(0.00, 1.44)
		CHAR	8.20(5.92, 2.28)	4.74(2.49, 2.25)	3.94(0.15, 3.79)	1.36(0.00, 1.36)
		TP	2.69(0.36, 2.33)
200	100	SEMI	2.17(0.12, 2.05)	0.82(0.09, 0.73)	5.53(0.00, 5.53)	1.16(0.04, 1.12)
		CHAR	8.20(6.41, 1.79)	3.70(1.91, 1.79)	5.34(0.01, 5.33)	1.04(0.02, 1.02)
		TP	2.67(0.34, 2.33)
	200	SEMI	1.54(0.17, 1.37)	0.78(0.05, 0.73)	3.70(0.01, 3.69)	0.58(0.01, 0.57)
		CHAR	6.10(4.87, 1.23)	3.15(1.92, 1.23)	3.53(0.01, 3.52)	0.54(0.00, 0.54)
		TP	1.86(0.36, 1.50)

Equations

$$f(x)=\frac{\exp \left\{\eta (x;\theta ,h)\right\}}{\int \exp \left\{\eta (x;\theta ,h)\right\}dx},$$

$$f_l(x,\theta ,h)=\frac{\exp \left\{{\eta }_l(x;\theta ,h)\right\}}{\int \exp \left\{{\eta }_l(x;\theta ,h)\right\}dx},$$

$$\sum _{l=1}^m\left\{-\frac{1}{n_l}\sum _{i=1}^{n_l}{\eta }_l(X_{li};\theta ,h)+\log {\int }_{{\mathcal{X}}_l}{e}^{{\eta }_l(x;\theta ,h)}dx\right\}+\frac{\lambda }{2}J(h),$$

$${\eta }_l(x;\theta ,h)={\alpha }_l(x;\theta )+h(x),$$

$$\sum _{l=1}^m\left\{-\frac{1}{n_l}\sum _{i=1}^{n_l}{\alpha }_l(X_{li};\theta )-\frac{1}{n_l}\sum _{i=1}^{n_l}h(X_{li})+\log {\int }_{{\mathcal{X}}_l}\exp \{{\alpha }_l(x;\theta )+h(x)\}dx\right\}+\frac{\lambda }{2}J(h).$$

$$\sum _{l=1}^m\left\{-\frac{1}{n_l}\sum _{i=1}^{n_l}h(X_{li})+\log {\int }_{{\mathcal{X}}_l}{w}_l(x)e^{h(x)}dx\right\}+\frac{\lambda }{2}J(h),$$

$${\hat{h}}_{\theta ,\lambda }(x)=\sum _{\nu =1}^k{d}^{\nu }{\psi }^{\nu }(x)+\sum _{j=1}^q{c}^j{R}_J(Z_j,x)={\psi }^T(x)d+{\xi }^T(x)c,$$

$${\bar{\mu }}_{\tilde{h}}(\psi )=\sum _{l=1}^m{\mu }_{l,\tilde{h}}(\psi )={\left[\sum _{l=1}^m{\mu }_{l,\tilde{h}}({\psi }^i)\right]}_{i=1}^k,V_{\psi ,\tilde{h}}=\sum _{l=1}^m{V}_{l,\tilde{h}}(\psi ,\tilde{h})={\left[\sum _{l=1}^m{V}_{l,\tilde{h}}({\psi }^i,\tilde{h})\right]}_{i=1}^k,$$

$$V_{\psi ,\psi }=\sum _{l=1}^m{V}_{l,\tilde{h}}(\psi ,{\psi }^T)={\left[\sum _{l=1}^m{V}_{l,\tilde{h}}({\psi }^i,{\psi }^j)\right]}_{i,j=1}^k,V_{\psi ,\xi }=\sum _{l=1}^m{V}_{l,\tilde{h}}(\psi ,{\xi }^T)={\left[\sum _{l=1}^m{V}_{l,\tilde{h}}({\psi }^i,{\xi }^j)\right]}_{i=1,j=1}^{k,q},$$

$$\left(\begin{array}[]{cc}V_{\psi,\psi}&V_{\psi,\xi}\\ V_{\xi,\psi}&V_{\xi,\xi}+\lambda Q\end{array}\right)\left(\begin{array}[]{c}d\\ c\end{array}\right)=\left(\begin{array}[]{c}S^{T}W\bm{1}_{N}-\bar{\mu}_{\tilde{h}}({\psi})+V_{\psi,\tilde{h}}\\ R^{T}W\bm{1}_{N}-\bar{\mu}_{\tilde{h}}({\xi})+V_{\xi,\tilde{h}}\end{array}\right).$$

$$\text{RKL}({\eta }_l(x;\theta ,h),{\eta }_l(x;\theta ,{\hat{h}}_{\theta ,\lambda }),l=1,\dots ,m)=\sum _{l=1}^m\left\{\log {\int }_{{\mathcal{X}}_l}{e}^{{\eta }_l(x;\theta ,{\hat{h}}_{\theta ,\lambda })}dx-{\mu }_{l,h}({\eta }_l(x;\theta ,{\hat{h}}_{\theta ,\lambda }))\right\}.$$

$$V(\lambda ,\theta )=\sum _{l=1}^m\left\{-\frac{1}{n_l}\sum _{i=1}^{n_l}{\eta }_l(X_i;\theta ,{\hat{h}}_{\theta ,\lambda })+\log {\int }_{{\mathcal{X}}_l}{e}^{{\eta }_l(x;\theta ,{\hat{h}}_{\theta ,\lambda })}dx+\alpha \frac{\text{tr}(P_{l,1}^{\perp }\breve{R_l}{H}^{-1}{\breve{R_l}}^T{P}_{l,1}^{\perp })}{n_l(n_l-1)}\right\},$$

$$\sum _{l=1}^m\left\{-\frac{1}{n_l}\sum _{i=1}^{n_l}({\alpha }_l(X_{li};\theta )+{\hat{h}}_{\theta }(X_{li}))+\log {\int }_{{\mathcal{X}}_l}\exp \{{\alpha }_l(x;\theta )+{\hat{h}}_{\theta }(x)\}dx\right\}+\frac{{\lambda }_{\theta }}{2}J({\hat{h}}_{\theta }).$$

$${\eta }_l(x;\theta ,h)\approx {\eta }_l(x;\theta ,\tilde{h})+{\mathcal{L}}_{l,x}(h-\tilde{h})\equiv r_l(x;\theta ,\tilde{h})+{\mathcal{L}}_{l,x}h,$$

$${\eta }_l(x;\theta ,h)=h({\alpha }_l(x;\theta )),$$

$$\sum _{l=1}^m\left\{-\frac{1}{n_l}\sum _{i=1}^{n_l}h(Y_{li})+\log \int \exp \left\{h(y)\right\}{w}_l(y)dy\right\}+\frac{\lambda }{2}J(h)$$

$$(\hat{\theta },\hat{h})\equiv \underset{(\theta ,h)\in \mathcal{Q}}{\mathrm{argmin}}\sum _{l=1}^m{\mathfrak{L}}_{n_l,\lambda }(\theta ,h)$$

$$\{\begin{array}{l}{D}_{\theta }{\sum }_{l=1}^m{\mathfrak{L}}_{n_l,\lambda }(\theta ,h)=0\\ D_h{\sum }_{l=1}^m{\mathfrak{L}}_{n_l,\lambda }(\theta ,h)=0\end{array}$$

$$\sum _{l=1}^m{\tilde{\mathfrak{L}}}_{n_l,\lambda }(\theta ,h)$$

$$\equiv \sum _{l=1}^m\{-\frac{1}{n_l}\sum _{i=1}^{n_l}\left[D_{\theta }{\eta }_l(X_{li};{\tau }_0)\theta +D_h{\eta }_l(X_{li};{\tau }_0)h\right]+{\mu }_{l,{\tau }_0}\left[D_{\theta }{\eta }_l({\tau }_0)\theta +D_h{\eta }_l({\tau }_0)h\right]$$

$$+\frac{1}{2}{V}_{l,{\tau }_0}\left[D_{\theta }{\eta }_l({\tau }_0)(\theta -{\theta }_0)+D_h{\eta }_l({\tau }_0)(h-h_0)\right]\}+\frac{\lambda }{2}J(h),$$

$$\{\begin{array}{l}{D}_{\theta }{\sum }_{l=1}^m{\tilde{\mathfrak{L}}}_{n_l,\lambda }(\theta ,h)=0\\ D_h{\sum }_{l=1}^m{\tilde{\mathfrak{L}}}_{n_l,\lambda }(\theta ,h)=0.\end{array}$$

$$\mathrm{SKL}({\tau }_0,\hat{\tau })$$

$$={\mu }_{{\tau }_0}[\eta ({\theta }_0,h_0)-\eta (\hat{\theta },\hat{h})]+{\mu }_{\hat{\tau }}[\eta (\hat{\theta },\hat{h})-\eta ({\theta }_0,h_0)].$$

$$M_1{\Vert g\Vert }^2\le \sum _{l=1}^m{V}_{l,{\tau }_0}(L_{l,h}g)+\lambda J(g)\le M_2{\Vert g\Vert }^2.$$

$$\sum _{l=1}^m{V}_{l,{\tau }_0}(L_{l,\theta }\zeta -L_{l,h}g)=\sum _{l=1}^m{V}_{l,{\tau }_0}(L_{l,\theta }\zeta -L_{l,h}g,L_{l,\theta }\zeta -L_{l,h}g)>c_{\delta }.$$

$$⟨({\zeta }_1,g_1),({\zeta }_2,g_2){⟩}_{\mathcal{Q}}\equiv \sum _{l=1}^m{V}_{l,{\tau }_0}(L_{l,\theta }{\zeta }_1+L_{l,h}{g}_1,L_{l,\theta }{\zeta }_2+L_{l,h}{g}_2)+\lambda J(g_1,g_2),$$

$$C_1\sum _{l=1}^m{V}_{l,{\tau }_0}(L_{l,\theta }\zeta +L_{l,h}g)$$

$$\sum _{l=1}^m{V}_{l,{\tau }_0}\left[(D_{\theta }{\eta }_l({\theta }_1,h_1)-D_{\theta }{\eta }_l({\theta }_2,h_2))\zeta +(D_h{\eta }_l({\theta }_1,h_1)-D_h{\eta }_l({\theta }_2,h_2))g\right]\le C_d\sum _{l=1}^m{V}_{l,{\tau }_0}(L_{l,\theta }\zeta +L_{l,h}g).$$

$$C_3\sum _{l=1}^m{V}_{l,{\tau }_0}(f)\le \sum _{l=1}^m{V}_{l,\tau }(f)\le C_4\sum _{l=1}^m{V}_{l,{\tau }_0}(f)$$

$$\sum _{l=1}^m{V}_{l,{\tau }_0}(L_{l,h}{\varphi }_{\nu },L_{l,h}{\varphi }_{\mu })={\delta }_{\nu \mu }\text{ and }J({\varphi }_{\mu },{\varphi }_{\nu })={\rho }_{\nu }{\delta }_{\nu \mu }$$

$$\sum _{l=1}^m{V}_{l,{\tau }_0}(D_h{\eta }_l({\theta }_{\ast },h_{\ast }){\varphi }_{\ast ,\nu },D_h{\eta }_l({\theta }_{\ast },h_{\ast }){\varphi }_{\ast ,\mu })={\delta }_{\nu \mu }\text{ and }J({\varphi }_{\ast ,\mu },{\varphi }_{\ast ,\nu })={\rho }_{\ast ,\nu }{\delta }_{\nu \mu }$$

Full text

$\hbox{University of California, Santa Barbara}^{1}$

$\hbox{University of Massachusetts, Amherst}^{2}$

1 Introduction

Density estimation plays a fundamental role in many areas of statistics and machine learning. The estimated density functions are useful for model building and diagnostics, inference, prediction and clustering. Traditionally there are two distinct approaches for density estimation: maximum likelihood method within a parametric family and nonparametric methods. Parametric models are in general parsimonious with meaningful and interpretable parameters [pearson1902systematic, pearson1902systematic2, kendall1946advanced, fisher1997absolute]. Nonparametric methods based on minimal assumptions are in general more flexible [Silverman, izenman1991review, Gu2013bk]. Often in practice it is desirable to model some components of the density function parametrically while leaving other components unspecified. Many semiparametric density models have been proposed for different purposes, including the partially linear semiparametric density models [yang2009penalized], the mixture models [olkin1987semiparametric, bordes2006semiparametric, ma2015flexible], mixture models with exponential tilts for multiple samples [Anderson1972, Qin1999, ZouFineYandell2002, Tan2009], combinations of parametric and nonparametric functions [hjort1995nonparametric, efron1996using], semiparametric models for Bayesian testing [lenk2003bayesian], and transformation models including the location-scale family distributions [wand1991transformations, potgieter2012nonparametric]. However, a systematic framework with a broad formulation of the semiparametric density models is lacking. In this paper, we propose a general framework that includes the above mentioned semiparametric density models as special cases. We develop a smoothing spline based estimation method for the general model and prove the asymptotic results based on our estimation method. The general framework provides a unified platform for the developments of estimation, computation, and theory.

For a single random sample, $X_{i}$, $i=1,\ldots,n$, from a common probability density $f(x)$ on a general domain ${\cal X}$, we consider the following general semiparametric density model

[TABLE]

where $\eta$ is a known function of $x$ given $\theta$ and $h$, which will be referred to as the logistic transformation of $f$. The parameter $\theta\in\mathbb{R}^{p}$ and the nonparametric function $h$ are unknown and need to be estimated. Often certain conditions, which depend on the form of $\eta(x;\theta,h)$, are necessary to make model (1) identifiable. We assume that model (1) is identifiable, and discuss identifiability conditions for specific models in the following sections.

Many existing semiparametric models are special cases of model (1). \citeasnounolkin1987semiparametric proposed a mixture of a parametric and a nonparametric density function, namely, $f(x)=\theta_{0}f_{1}(x,\theta_{1})+(1-\theta_{0})f_{2}(x)$, where $f_{1}(x,\theta_{1})$ is a known density function up to parameters $\theta_{1}\in\mathbb{R}^{p}$, $\theta_{0}\in[0,1]$ is an unknown weight parameter, and $f_{2}(x)$ is a nonparametric density function. They showed that the semiparametric density estimate provides a compromise between the parametric and nonparametric estimates. It is a special case of model (1) with $\eta(x)=\log\{\theta_{0}f_{1}(x,\theta_{1})+(1-\theta_{0})f_{2}(x)\}$, where $\theta_{0}$ and $\theta_{1}$ are the parameters and $f_{2}$ is the nonparametric function. \citeasnounhjort1995nonparametric proposed a density estimation procedure by starting out with a parametric density estimate $f_{1}(x,\hat{\theta})$, and then multiplying with a nonparametric kernel type estimate of a correction function $r(x)=f(x)/f_{1}(x,\hat{\theta})$. It was shown that their semiparametric density model can perform better than a nonparametric fit when the true density is in the neighborhood of the initial parametric density. Their model is a special case of model (1) with $\eta(x)=\log\{f_{1}(x,\theta)\}+h(x)$, where $\theta$ is the parameter and $h(x)=\log\{r(x)\}$ is the nonparametric function. \citeasnounefron1996using proposed a specially designed exponential family for density estimation, $f(x)=f_{0}(x)\exp\{\theta_{0}+t^{T}(x)\theta_{1}\}$, where $f_{0}(x)$ is a carrier density and estimated by kernel density estimation, $t(x)$ is a $p$-dimensional vector of sufficient statistics, $\theta_{1}$ is a $p$-dimensional vector of parameters, and $\theta_{0}$ is a normalizing parameter making $f(x)$ integrate to 1 over $\mathcal{X}$. The proposed method matches the estimated expectation of $t(x)$ with its sample expectation. The model is a special case of model (1) with $\eta(x)=\log\{f_{0}(x)\}+\theta_{0}+t^{T}(x)\theta_{1}$, where $f_{0}(x)=\exp\{h(x)\}/\int{\exp\{h(x)\}dx}$ is the carrier density given by nonparametric function $h(x)$, and $\theta_{0}$ and $\theta_{1}$ are the parameters. \citeasnounlenk2003bayesian proposed a flexible semiparametric model for Bayesian testing of $f(x|\theta,Z)={\exp\{\alpha^{T}(x)\theta+Z(x)\}}/{\int_{\mathcal{X}}\exp\{\alpha^{T}(x)\theta+Z(x)\}dG(x)}$, where $\alpha(x)=[\alpha^{1}(x),\ldots,\alpha^{m}(x)]^{T}$ is a vector of $m$ nonconstant functions, $Z$ is a zero mean, second-order Gaussian process with bounded, continuous covariance function, and $G$ is a known dominating measure on the support $\mathcal{X}$. The semiparametric model allows the predictive distribution to deviate from the parametric family. If the parametric family is inadequate, the semiparametric predictive density coherently adapts to the data. \citeasnounyang2009penalized also used the logistic transformation of density function as in \citeasnounlenk2003bayesian with $Z(x)$ modeled as an unknown smooth function. The model in \citeasnounyang2009penalized is a special case of model (1) with $\eta(x)=\alpha^{T}(x)\theta+h(x)$. \citeasnounwand1991transformations considered density estimation of data with local features, and proposed a data transformation technique so that global smoothing is appropriate. This transformation model is a special case of model (1) with $\eta(x)=h(t(x;\theta))$ where $t(x;\theta)$ is the transformation, $\theta$ is the parameter and $h$ is the nonparametric function.

In the case of multiple samples, assume there are $m$ independent groups, and in each group $l=1,\cdots,m$, there are $n_{l}$ iid observations such that $X_{l1},\ldots,X_{ln_{l}}\overset{iid}{\sim}f_{l}(x;\theta,h)$ on domain ${\cal X}_{l}$. We consider the following general semiparametric density model

[TABLE]

where $\eta_{l}$, the logistic transformation of $f_{l}$, is a known function of the parameter $\theta$ and the nonparametric function $h$. We are interested in the estimation of $\theta\in\mathbb{R}^{p}$, $h$, and ultimately the overall density function $f_{l}(x,\theta,h)$ from the observations.

\citeasnoun

potgieter2012nonparametric considered a two sample transformation model. Suppose that $X_{11},\ldots,X_{1n_{1}}\overset{iid}{\sim}f_{1}$ and $X_{21},\ldots,X_{2n_{2}}\overset{iid}{\sim}f_{2}$, and $X_{2}$’s have the same distribution as $X_{1}$’s after a certain invertible transformation parametrized by $\theta$, denoted as $t(x;\theta)$. They considered the nonparametric estimation of the density function based on asymptotic likelihood, and showed that the estimators are often near optimal when compared to fully parametric methods. The model is a special case of our model (2) with $\eta_{1}(x;\theta,h)=h(x)$, and $\eta_{2}(x;\theta,h)=h(t(x;\theta))\left\lvert t^{\prime}(x;\theta)\right\rvert$, where $h$ is the logistic transformation of $f_{1}$. \citeasnounAnderson1972, \citeasnounQin1999, \citeasnounZouFineYandell2002, and \citeasnounTan2009 considered exponential tilt mixture models for biased sampling and case control studies, where multiple samples are collected and each sample follows an exponential tilt mixture model, that is, for sample $l$, $X_{l1},\ldots,X_{ln_{1}}\overset{iid}{\sim}\lambda_{l}f(x)+(1-\lambda_{l})g(x)$, and $\log(f(x)/g(x))=\theta^{1}+\theta^{2}x$. We can see that this model is also a special case of model (2) with $\eta_{l}(x;\theta,h)=\log\{\lambda_{l}+(1-\lambda_{l})\exp(\theta^{1}+\theta^{2}x)\}+h(x)$, where $\theta=[\theta^{1},\theta^{2}]$ is the parameter and $h$ is the logistic transformation of $f$.

We consider the smoothing spline based estimation method for models (1) and (2), that is, we assume that $h\in{\cal H}$, where ${\cal H}$ is a reproducing kernel Hilbert space (RKHS). Our work builds upon the existing literature and extends it to include semiparametric nonlinear density models. We develop computation methods based on profile penalized likelihood and backfitting, and joint asymptotic theory of the parametric and nonparametric components. For the single sample case, smoothing spline based nonparametric density estimation has been considered by many authors, including \citeasnounsilverman_1982, \citeasnounCoxOsullivan, and \citeasnounGu2013bk. \citeasnounyang2009penalized extended such methods for nonparametric models to a semiparametric density estimation model, where $\eta$ is linear in both the parametric and nonparametric components $\theta$ and $h$, as discussed above. Great progress has been made toward joint asymptotic theory in semiparametric models recently, with the seminal work by \citeasnouncheng2015joint and follow-up papers [ZhaoChengLiu16, ChaoVolgushevCheng17]. Prior to this, semiparametric asymptotic theory focused on the convergence of the parametric component, whereas the nonparametric component was usually considered as a nuisance parameter [Bickelbk, Tsiatisbk, Kosorokbk]. We note that one approach is to extend the framework for joint asymptotics of the semiparametric linear regression model in \citeasnouncheng2015joint to the semiparametric linear density model, and the same joint rate of convergence as derived by \citeasnounyang2009penalized can be obtained. This approach relies heavily on the extension of the reproducing kenel of ${\cal H}$, where the nonparametric component $h$ lives, to the product parameter space $\mathbb{R}^{p}\times{\cal H}$ using the linear form of $\eta$ in $\theta$ and $h$. However, such an extension for a general nonlinear function $\eta(\theta,h)$ is nontrivial to obtain. In order to develop joint asymptotics for model (1) with general $\eta$ and for model (2) with general $\eta_{l}$, we extend an approach from \citeasnoungu_qiu_1993 for nonparametric density estimation to the general nonlinear semiparametric model through a first order linearization technique motivated by the study of the nonlinear nonparametric regression model in \citeasnouno’sullivan_1990. In particular, we specify the additional assumptions needed for our general framework and we introduce an appropriate metric in which the rate of convergence is derived for the joint estimator. Additional assumptions include, for example, smoothness criteria for $\eta$ and the existence of certain bounded linear operators related to the first Fréchet partial derivatives of $\eta$, all of which are redundant when $\eta$ is linear in $(\theta,h)$. In addition to the rate of convergence of the joint estimator, we also obtain the convergence rate of the parametric component in the standard $l^{2}$ norm as well as the convergence rate of the overall density function in the symmetrized Kullback-Leibler distance.

In Section 3, we consider models (1) and (2) in three scenarios and develop different computation procedures. Asymptotic properties of the proposed estimators are considered in Section 4. Simulations are conducted to evaluate the proposed estimation procedures in Section 5. Section 6 shows applications to real dataset.

2 Some notations

We begin by introducing some notations that will be used throughout the rest of the paper. Unless specified otherwise, let ${a=\left[a^{k}\right]_{k=1}^{p}}$ denote any $p\times 1$ vector, whose the $k$th element is $a^{k}$. The standard $l^{2}$ norm of $a\in\mathbb{R}^{p}$ is denoted $\left\lVert a\right\rVert_{l^{2}}=[\sum_{k=1}^{p}(a^{k})^{2}]^{1/2}$. We also denote any $p\times p$ matrix ${M=\left[M^{i,j}\right]_{i,j=1}^{p}}$, where $M^{i,j}$ represent the $(i,j)$th entry of the matrix. If there exist positive constants $c_{1},c_{2}$, such that $c_{1}A\leq B\leq c_{2}A$, we write $A\sim B$. We use the notation $\operatorname{E}(\cdot)$ to represent the expectation taken over the joint sample distribution, whose density function is $f(x_{1},\ldots,x_{m})=\prod_{l=1}^{m}f_{l}(x_{l})$. For simplicity, we will also sometimes use $\tau$ to represent the pair $(\theta,h)$, i.e., $\tau=(\theta,h)$ and $\eta_{l}(\tau)=\eta_{l}(\theta,h)$.

Denote the product parameter space as $\mathcal{Q}\equiv\mathbb{R}^{p}\times{\cal H}$. Let $D_{h}$ be the Fréchet partial differential operator with respect to $h$, and let $D_{\theta}$ be the Fréchet partial differential operator with respect to $\theta$. If $X$ and $Y$ are any (real) Banach spaces, $\mathscr{L}(X,Y)$ represents the space of bounded linear operators from $X$ to $Y$. Note that for any function $f:\mathcal{Q}\rightarrow Y$, the Fréchet partial derivatives of $f$ are maps $D_{\theta}f:\mathcal{Q}\rightarrow\mathscr{L}(\mathbb{R}^{p},Y)$ and $D_{h}f:\mathcal{Q}\rightarrow\mathscr{L}({\cal H},Y)$ by definition.

3 Penalized likelihood estimation

In this section, we will focus on model (2). Model (1) is a special case with $m=1$. We propose to estimate $\theta=[\theta^{i}]_{i=1}^{p}\in\mathbb{R}^{p}$ and $h\in{\cal H}$ in (2) as the minimizer of the following penalized likelihood with respect to $(\theta,h)\in\mathcal{Q}$:

[TABLE]

where the first part of (3) is the joint negative log likelihood, $J(h)$ is a square seminorm penalty, and $\lambda$ is the smoothing parameter. Assume that $\mathcal{H}_{0}=\{h:J(h)=0\}$ is a $k$-dimensional space with basis functions $\psi(x)=[\psi^{i}(x)]_{i=1}^{k}$, then ${\cal H}={\cal H}_{0}\oplus{\cal H}_{1}$, where ${\cal H}_{1}$ is an RKHS with reproducing kernel (RK) denoted by $R_{J}$.

The minimizer of the penalized likelihood (3) does not fall in a finite dimensional space. In the following, we consider estimation of the model (2) under the following three scenarios: additive, general, and transformation models.

3.1 The additive model

The model (2) is additive when $\eta_{l}$ is the sum of a term involving the parameters and the nonparametric component, that is,

[TABLE]

where $\alpha_{l}$ is a known function of $x$ with unknown parameters $\theta$. The model proposed by \citeasnounefron1996using, \citeasnounlenk2003bayesian, \citeasnounyang2009penalized, and \citeasnounhjort1995nonparametric are special cases of $\eqref{additive}$ with $m=1$. In particular, the first three of them considered the linear additive model in which $\alpha_{l}(x;\theta)=\theta^{T}a(x)$, where $a(x)=[a^{i}(x)]_{i=1}^{p}$. The model proposed in \citeasnounZouFineYandell2002 is an example of the nonlinear additive model with $\alpha_{l}(x)=\log\{\lambda_{l}+(1-\lambda_{l})\exp(\theta^{1}+\theta^{2}x)\}$, where $0\leq\lambda_{l}\leq 1$ and $\lambda_{1}\neq\cdots\neq\lambda_{m}$.

For the additive semiparametric density model (4), the penalized likelihood objective (3) becomes the following:

[TABLE]

We propose a profiled penalized likelihood based approach to optimize (5). First, with a fixed ${\theta}$, optimizing (5) with respect to $h\in{\cal H}$ is equivalent to optimizing the following weighted penalized likelihood:

[TABLE]

where $w_{l}(x)=\exp\{\alpha_{l}(x,\theta)\}$. As in Gu (2013), we approximate the solution by

[TABLE]

where $\{Z_{j},j=1,\ldots,q\}$ is a random sample of $\{X_{li},i=1,\ldots,n_{l},l=1,\ldots,m\}$, $c=[c^{j}]_{j=1}^{q}$, ${d=[d^{\nu}]_{\nu=1}^{k}}$, $\psi(x)=[\psi^{\nu}(x)]_{\nu=1}^{k}$, and $\xi(x)=[R_{J}(Z_{1},x),\ldots,R_{J}(Z_{q},x)]^{T}$.

For $l=1,\ldots,m$ and $i=1,\ldots,n_{l}$, let $N=\sum_{l=1}^{m}n_{l}$, $W$ be the $N\times N$ diagonal matrix whose $r$th diagonal entry is $n_{l}$ where $r=(l-1)n_{l}+i$, $S$ be the $N\times k$ matrix whose $(r,\nu$)th entry is $\psi^{\nu}(X_{li})$ for $\nu=1,\ldots,k$, and $R$ be the $N\times q$ matrix whose $(r,j)$th entry is $\xi^{j}(X_{li})=R_{J}(Z_{j},X_{li})$ for $j=1,\ldots,q$. Let $Q$ be the $q\times q$ matrix whose $(i,j)$th entry is $R_{J}(Z_{i},Z_{j})$. Define for $l=1,\ldots,m$, $\mu_{l,f}(g)={\int gw_{l}e^{f}dx}/{\int w_{l}e^{f}dx}$, $V_{l,f}(g,h)=\mu_{l,f}(gh)-\mu_{l,f}(g)\mu_{l,f}(h)$, and $V_{l,f}(g)=V_{l,f}(g,g)$. These notations extend in the obvious way to vectors and matrices of functions.

We plug (7) back into (6) and compute minimizers $c$ and $d$ using the Newton iterative method. Let $\tilde{h}_{\theta,\lambda}(x)=\psi^{T}(x)\tilde{d}+\xi^{T}(x)\tilde{c}$, where $\tilde{c}=[\tilde{c}^{j}]_{j=1}^{q}$ and $\tilde{d}=[\tilde{d}^{\nu}]_{\nu=1}^{k}$ are $c$ and $d$ vectors calculated in the previous step of the iteration. Denote

[TABLE]

[TABLE]

and note $V_{\xi,\psi}$ is the transpose of $V_{\psi,\xi}$. We can also define $\bar{\mu}_{\tilde{h}}({\xi})$, $V_{\xi,\tilde{h}}$, and $V_{\xi,\xi}$ similarly to the first three equations above by replacing $\psi$ by $\xi$ and $k$ by $q$. Similar to \citeasnounGu2013bk, the Newton updating equation is

[TABLE]

At convergence of the Newton method, we denote the estimate of $h$ with fixed $\theta$ and $\lambda$ as $\hat{h}_{\theta,\lambda}$. Similar to \citeasnounGu2013bk, the smoothing parameter $\lambda$ is selected by optimizing the relative Kullback-Leibler (KL) distance between the density associated with the estimate, $\eta_{l}(x;\theta,\hat{h}_{\theta,\lambda})$, and the true density corresponding to $\eta_{l}(x;\theta,h)$ for $l=1,\ldots,m$,

[TABLE]

Following similar derivations in \citeasnounGu2013bk, an approximated cross-validation estimate of the relative Kullback-Leibler distance is

[TABLE]

where $P_{l,\bf 1}^{\perp}=I_{n_{l}}-\mathbb{1}_{n_{l}}\mathbb{1}_{n_{l}}^{T}/n_{l}$, $\breve{\xi}=(\psi^{T},\xi^{T})^{T}$, ${\breve{R}_{l}^{T}=(\breve{\xi}(X_{l1}),\ldots,\breve{\xi}(X_{ln_{l}}))}$, $H=\sum_{l=1}^{m}V_{l,\tilde{h}}(\breve{\xi},\breve{\xi}^{T})+\text{diag}(0,\lambda Q)$, and $\alpha=1$. We may set a larger $\alpha$ to prevent under-smoothing, such as $\alpha=1.4$ as in \citeasnounGu2013bk.

The optimal smoothing parameter for a given $\theta$, $\lambda_{\theta}$, is defined as the $\lambda$ that minimizes the cross validation criterion (10). The corresponding nonparametric function estimate $\hat{h}_{\theta,\lambda_{\theta}}$ is denoted as $\hat{h}_{\theta}$ for simplicity. Plugging $\hat{h}_{\theta}$ into (5), we have the profiled penalized likelihood

[TABLE]

Then the estimate of $\theta$, $\hat{\theta}$, is the minimizer of (11). The minimization is achieved by the line search algorithm in \citeasnounnelder1965simplex. The final estimate of $h$ is $\hat{h}_{\hat{\theta}}$.

In fact, for a given $\theta$, minimizing the cross validation criterion (10) to obtain $\lambda_{\theta}$ and minimizing (5) to obtain $\hat{h}_{\theta}$ can be easily implemented by utilizing the ssden function in the R package gss with the option bias to specify the weights [gss].

3.2 The general nonlinear model

In this section, we consider the general form of $\eta_{l}$ in model (2) that is not additive. For example, when $m=1$, a mixture model of two densities assumes $\eta(x;\theta,h)=\log\left\{\theta_{0}f_{1}(x,\theta_{1})+(1-\theta_{0})h(x)\right\}$, where $0\leq\theta_{0}\leq 1$ and $\theta=[\theta_{0},\theta_{1}^{1},\ldots,\theta_{1}^{p}]^{T}$. Originally considered by \citeasnounolkin1987semiparametric, this model is often referred to as the two-component mixture model. Different estimation methods and applications can be found in \citeasnounbordes2006semiparametric and \citeasnounma2015flexible.

We will extend the Gauss-Newton procedure for the estimation of $\theta$ and $h$ to the general nonlinear semiparametric density model. For fixed $\theta$, let $\tilde{h}$ be the current estimate of $h$. We assume that the Fréchet derivative of $\eta_{l}$ with respect to $h$ at $(\theta,\tilde{h})$ exists, and we denote it by $D_{h}\eta_{l}(x;\theta,\tilde{h})\equiv{\cal L}_{l,x}$. Approximating $\eta_{l}$ by its first order Taylor expansion at $\theta$ and $\tilde{h}$, we have

[TABLE]

where $r_{l}(x;\theta,\tilde{h})=\eta_{l}(x;\theta,\tilde{h})-{\cal L}_{l,x}\tilde{h}$.

Through this linearization, the nonlinear semiparametric density is approximated by an additive model as in Section 3.1 but with the evaluational functional replaced by the bounded linear functional ${\cal L}_{l,x}$. Note that the method in Section 3.1 can be extended to bounded linear functionals with $\psi^{\nu}(X_{li})$ replaced by ${\cal L}_{l,X_{li}}\psi^{\nu}$ and $\xi^{j}(X_{li})$ replaced by ${\cal L}_{l,X_{li}}R_{J}(Z_{j},\cdot)$, $j=1,\cdots,q$.

In summary, for the general nonlinear semiparametric density model, there are two loops of iterations. In the inner loop, for a fixed $\theta$, we estimate $h$ by the extended Gaussian-Newton iteration until convergence, where smoothing parameters are estimated at each iteration. In the outer loop, the profile penalized likelihood (11) is optimized with respect to $\theta$.

3.3 The transformation model

Model (2) is called a transformation model if

[TABLE]

where $\alpha_{l}$’s are known differentiable and invertible functions with unknown parameters $\theta$. When $\theta$ is fixed, the logistic density of $Y_{l}=\alpha_{l}(X_{l};\theta)$ is $\eta_{l}(y)=h(y)-\log\left(|\alpha_{l}^{\prime}(\alpha_{l}^{-1}(y;\theta))|\right)$. Then the problem reduces to nonparametric density estimation with different weights. The two sample location-scale family density model [potgieter2012nonparametric] belongs to this case with $\alpha_{1}(x;\theta)=x$ and $\alpha_{2}(x;\theta)=(x-\mu)/\sigma$, where $\mu$ and $\sigma$ are the location and scale parameters and $\theta=(\mu,\sigma)$.

We note that the profiled penalized likelihood method developed for the general nonlinear model in Section 3.2 applies to the transformation model. Nevertheless, the following backfitting procedure works better for the transformation model. First, provide initial value $\theta_{0}$ for the parameter $\theta$. Next we iterate between the following two steps until convergence:

Based on current estimates $\tilde{\theta}$, transform the data using $Y_{li}=\alpha_{l}(X_{li};\tilde{\theta})$ and we have $\eta_{Y}(y)=h(y)-\log(|\alpha^{\prime}(\alpha^{-1}(y;\tilde{\theta}))|)$. With the transformed data, update $h$ by minimizing the penalized likelihood

[TABLE]

with the weight function $w_{l}(y)=1/|\alpha_{l}^{\prime}(\alpha_{l}^{-1}(y;\tilde{\theta}))|$. 2. 2.

Based on current estimate $\tilde{h}$ of $h$, update $\theta$ as the MLE based on the likelihood of $X_{11},\ldots,X_{mn_{m}}$.

The minimization of the weighted penalized likelihood (14) is solved by the Newton method proposed in \citeasnounGu2013bk.

4 Joint asymptotic theory

In this section, we develop a joint asymptotic theory for the penalized likelihood estimators for models (1) and (2).

For $1\leq l\leq m$, $f_{l}(x)$ is a density function over the sample space $\mathcal{X}_{l}$ of the form given by (2). We consider $(\theta,h)\in\mathcal{Q}=\mathbb{R}^{p}\times{\cal H}$ and a known function $\eta_{l}:\mathcal{Q}\to L^{2}_{0}(\mathcal{X}_{l})$, where $L_{0}^{2}(\mathcal{X}_{l})=L_{\tau_{0}}^{2}(\mathcal{X}_{l})\ominus\{1\}$, $L^{2}_{\tau_{0}}(\mathcal{X}_{l})$ is the space of functions on $\mathcal{X}_{l}$ with finite second moment with respect to the measure given by the true density $f_{l}(x;\theta_{0},h_{0})$, and $\tau_{0}=(\theta_{0},h_{0})$ is the true parameter. The constant functions have been removed for identifiability. In general, $\eta_{l}(\theta,h)$ need not be one-to-one on $\mathcal{Q}$. However, for identifiability, we require that for all $l=1,\ldots,m$, $\eta_{l}(\theta,h)$ is one-to-one in a neighborhood $\mathcal{N}_{\theta_{0}}\times\mathcal{N}_{h_{0}}$ of the true parameter. We specify the precise assumptions on this neighborhood in Section 4.2.2.

Recall the penalized likelihood for estimating the parameter $(\theta,h)$ given by (3), and denote it as $\sum_{l=1}^{m}\mathfrak{L}_{n_{l},\lambda}(\theta,h)$. Let $\mathfrak{L}_{n_{l}}(\theta,h)\equiv{-n_{l}^{-1}\sum_{i=1}^{n_{l}}\eta_{l}(X_{li};\theta,h)+\log\int_{\mathcal{X}_{l}}\exp\{\eta_{l}(x;\theta,h)\}dx}$. If the joint negative log likelihood $\sum_{l=1}^{m}\mathfrak{L}_{n_{l}}(\theta,h)$ is convex, under certain conditions (see Theorem $2.9$ in \citeasnounGu2013bk), the existence of the penalized semiparametric estimator given by

[TABLE]

is guaranteed. In general, the functions $\eta_{l}(\theta,h)$ need not satisfy these conditions, so for our analysis on joint consistency, we assume that $\sum_{l=1}^{m}\mathfrak{L}_{n_{l}}(\theta,h)$ is convex and achieves a unique local minimum at $(\hat{\theta},\hat{h})$ in $\mathcal{N}_{\theta_{0}}\times\mathcal{N}_{h_{0}}$. We study the consistency of this estimator and obtain the rate of convergence of $(\hat{\theta},\hat{h})$ to the true parameters $(\theta_{0},h_{0})$ as $n=\min_{1\leq l\leq m}\{n_{l}\}\to\infty$ and $\lambda\to 0$. In our analysis, we may drop the subscript $l=1$ when $m=1$ for convenience.

Note that by adapting the methods in \citeasnouncox_1988, \citeasnounCoxOsullivan, and \citeasnouno’sullivan_1990, we have established local existence and uniqueness of $(\hat{\theta},\hat{h})$. A detailed proof can be found in the Appendix.

4.1 Outline of the proof of consistency

By our assumption for $\sum_{l=1}^{m}\mathfrak{L}_{n_{l},\lambda}(\theta,h)$, the estimator $\hat{\tau}=(\hat{\theta},\hat{h})$ is taken to be the unique solution of

[TABLE]

in $\mathcal{N}_{\theta_{0}}\times\mathcal{N}_{h_{0}}$. To study the asymptotic behavior of the estimator $\hat{\tau}=(\hat{\theta},\hat{h})$, we follow a similar outline as in the nonparametric case [gu_qiu_1993]. We first introduce an approximation of $\hat{\tau}$, denoted $\tilde{\tau}=(\tilde{\theta},\tilde{h})$, which is the minimizer of

[TABLE]

where $\mu_{l,\tau}(\varphi)={\int_{\mathcal{X}_{l}}\varphi(t)\exp(\eta_{l}(t;\tau))dt}/{\int_{\mathcal{X}_{l}}\exp(\eta_{l}(t;\tau))dt}$, $V_{l,\tau}(\varphi_{1},\varphi_{2})=\mu_{l,\tau}(\varphi_{1}\varphi_{2})-\mu_{l,\tau}(\varphi_{1})\mu_{l,\tau}(\varphi_{2})$, and $V_{l,\tau}(\varphi)=V_{l,\tau}(\varphi,\varphi)$. We see that for each $l=1,\ldots,m$, $\tilde{\mathfrak{L}}_{n_{l}}(\theta,h)$ is almost like a quadratic approximation of $\mathfrak{L}_{n_{l}}(\theta,h)$ at $\tau_{0}=(\theta_{0},h_{0})$, ignoring terms independent of $(\theta,h)$ and terms involving second derivatives of $\eta_{l}(\theta,h)$. Since $V_{l,\tau_{0}}(\cdot)$ and $J(\cdot)$ are quadratic functionals, one can check that $\sum_{l=1}^{m}\tilde{\mathfrak{L}}_{n_{l},\lambda}(\theta,h)$ is convex with respect to $(\theta,h)$, and attains its minimum at $\tilde{\tau}=(\tilde{\theta},\tilde{h})$ if and only if $\tilde{\tau}$ is the solution for the system

[TABLE]

Using a quadratic approximation of the penalized likelihood to attain an approximation of the estimator is a common intermediate step in the study of convergence of smoothing spline estimators (see \citeasnounsilverman_1982, \citeasnounCoxOsullivan, \citeasnouno’sullivan_1990, and \citeasnoungu_qiu_1993). In particular, when $m=1$, our definition of $\tilde{\mathfrak{L}}_{n_{1},\lambda}(\theta,h)$ can be seen as a generalized version of the quadratic approximation given by equation (5.2) in \citeasnoungu_qiu_1993.

In Section 4.2, we define a norm related to $\sum_{l=1}^{m}V_{l,\tau_{0}}(\cdot)$ and $J(\cdot)$, in which the rate of convergence for $\hat{\tau}-\tau_{0}$ can be derived. For readability, we restrict our analysis to the special case when $m=1$ and establish the rate of convergence for $\tilde{\tau}-\tau_{0}$ in Section 4.3. Together with a bound for the approximation error $\hat{\tau}-\tilde{\tau}$, the rate of convergence for $\hat{\tau}-\tau_{0}$ then follows by the triangle inequality (Section 4.4). In addition, we also establish the rate of convergence of the estimate $\hat{\tau}$ in terms of the symmetrized Kullback-Leibler distance, which is defined as

[TABLE]

Section 4.5 states results for $m\geq 1$ without proofs since the proofs are straightforward extensions of those in Sections 4.3 and 4.4.

4.2 Norm and assumptions

In this section, we introduce an appropriate norm in which the convergence rate of our estimator is measured and we provide some assumptions needed for our analysis.

4.2.1 Parameter space

Assumption 1**.**

The penalty $J(h)$ is a square seminorm in ${\cal H}$ with a finite-dimensional null space ${\cal H}_{0}\subset{\cal H}$. Therefore, $J((\theta,h))\equiv J(h)$ extends $J$ to a square seminorm on $\mathcal{Q}$, and its null space $\mathbb{R}^{p}\times{\cal H}_{0}$ is again finite-dimensional. Denote by $J(g,h)$ the semi-inner product associated with the seminorm $J(h)$. We also assume that $J(h_{0})<\infty$.

Assumption 2**.**

For $l=1,\ldots,m$, there are bounded linear operators $L_{l,\theta}:\mathbb{R}^{p}\to L^{2}_{0}(\mathcal{X}_{l})$ and $L_{l,h}:{\cal H}\to L^{2}_{0}(\mathcal{X}_{l})$, with zero nullspaces, which satisfy the following conditions:

1. (i)

Suppose ${\cal H}$ is a real Hilbert space of functions, equipped with norm $\left\lVert\cdot\right\rVert$. For any $g\in{\cal H}$, there exist positive constants $M_{1},M_{2}$, such that

[TABLE] 2. (ii)

For any $\zeta\in\mathbb{R}^{p}$ satisfying $\left\lVert\zeta\right\rVert_{l^{2}}=1$ and for any $g\in{\cal H}$, there exists a positive constant $c_{\delta}$ such that

[TABLE]

By Assumptions 1 and 2$(i)$, we see that $\langle g_{1},g_{2}\rangle_{{\cal H}}\equiv\sum_{l=1}^{m}V_{l,\tau_{0}}(L_{l,h}g_{1},L_{l,h}g_{2})+\lambda J(g_{1},g_{2})$ is an inner product on ${\cal H}$, and its induced norm, denoted by $\left\lVert\cdot\right\rVert_{{\cal H}}$, is complete on ${\cal H}$. One can also see that $L_{l,\theta}$ can be represented by the $p\times 1$ vector of $L^{2}_{0}(\mathcal{X}_{l})$ functions $[L_{l,\theta}^{k}(x)]_{k=1}^{p}$. We may use $L_{l,\theta}$ to denote the linear operator or its vector form, i.e., $L_{l,\theta}\zeta=\zeta^{T}L_{l,\theta}$. Denote by $V_{l,\tau_{0}}(L_{l,\theta},L_{l,\theta})$ the $p\times p$ matrix in which the $(i,j)$th entry is $V_{l,\tau_{0}}(L_{l,\theta}^{i}(x),L_{l,\theta}^{j}(x))$, and note that one can write $\sum_{l=1}^{m}V_{l,\tau_{0}}(L_{l,\theta}\zeta,L_{l,\theta}\zeta)=\zeta^{T}\sum_{l=1}^{m}V_{l,\tau_{0}}(L_{l,\theta},L_{l,\theta})\zeta$. When $g=0$, (17) implies that $\sum_{l=1}^{m}V_{l,\tau_{0}}(L_{l,\theta},L_{l,\theta})$ is positive definite. Therefore, $\langle\zeta_{1},\zeta_{2}\rangle_{\mathbb{R}^{p}}\equiv\sum_{l=1}^{m}V_{l,\tau_{0}}(L_{l,\theta}\zeta_{1},L_{l,\theta}\zeta_{2})$ is an inner product on $\mathbb{R}^{p}$ and we use $\left\lVert\cdot\right\rVert_{\mathbb{R}^{p}}$ to denote its induced norm.

For any $(\zeta_{1},g_{1}),(\zeta_{2},g_{2})\in\mathcal{Q}$, we define an inner product as follows,

[TABLE]

and we denote the norm induced by this inner product by $\left\lVert\cdot\right\rVert_{\mathcal{Q}}$. In the Appendix, we prove that $\langle\cdot,\cdot\rangle_{\mathcal{Q}}$ as defined above is, in fact, an inner product, and that its induced norm is complete.

Remark 1**.**

1. (i)

We need to use the auxiliary operators $L_{l,\theta}$ and $L_{l,h}$ to define the norm in which the rate of convergence can be conveniently measured. We point out here that $L_{l,\theta}$ and $L_{l,h}$ are related to the Fréchet partial derivatives of $\eta_{l}$ by Assumption 4 in Section 4.2.2. This is analogous to the approach in \citeasnouno’sullivan_1990. 2. (ii)

It is well known that all norms on $\mathbb{R}^{p}$, including $\left\lVert\cdot\right\rVert_{\mathbb{R}^{p}}$ as defined above, are equivalent to the Euclidean norm, namely $\left\lVert\zeta\right\rVert_{l^{2}}=(\sum_{i=1}^{p}\left\lvert\zeta^{i}\right\rvert^{2})^{1/2}$ (Theorem 3.1 in \citeasnounconway1990). Therefore $\left\lVert\cdot\right\rVert_{\mathbb{R}^{p}}$ is complete. Moreover, the convergence of $\hat{\theta}$ to $\theta_{0}$ in $\left\lVert\cdot\right\rVert_{\mathbb{R}^{p}}$ implies the convergence in the Euclidean norm (see Corollary 1). 3. (iii)

Assumption 2(ii) is a regularity condition similar to Assumption A3. in \citeasnouncheng2015joint. This assumption guarantees that $\left\langle\cdot,\cdot\right\rangle_{\mathcal{Q}}$ is indeed an inner product, and its induced norm is complete (see Theorem A.1 in the Appendix). 4. (iv)

We note that under Assumption 2, for any $(\zeta,g)\in\mathcal{Q}$, $\left\lVert(\zeta,g)\right\rVert_{\mathcal{Q},1}\equiv\left\lVert\zeta\right\rVert_{\mathbb{R}^{p}}+\left\lVert g\right\rVert_{{\cal H}}$ defines a norm on the product space $\mathcal{Q}$. By a similar argument as in the proof of Theorem A.1 in the Appendix, a sequence $\{(\zeta_{i},g_{i})\}_{i=1}^{\infty}\in\mathcal{Q}$ converges in norm $\left\lVert\cdot\right\rVert_{\mathcal{Q}}$ if and only if it converges in norm $\left\lVert\cdot\right\rVert_{\mathcal{Q},1}$. This implies that $\left\lVert\cdot\right\rVert_{\mathcal{Q}}$ and $\left\lVert\cdot\right\rVert_{\mathcal{Q},1}$ are equivalent norms on $\mathcal{Q}$.

4.2.2 Properties of $\eta(\theta,h)\text{ and }(\theta_{0},h_{0})$

Recall $\tau_{0}=(\theta_{0},h_{0})$ is the true parameter in $\mathcal{Q}$. We assume there are neighborhoods $\mathcal{N}_{\theta_{0}}\subset\mathbb{R}^{p}$ of $\theta_{0}$ and $\mathcal{N}_{h_{0}}\subset{\cal H}$ of $h_{0}$ such that the following Assumptions 3 to 6 hold.

Assumption 3**.**

For all $l=1,\ldots,m$, $\eta_{l}(\theta,h)$ is one-to-one and three times continuously Fréchet differentiable with respect to $(\theta,h)$ in $\mathcal{N}_{\theta_{0}}\times\mathcal{N}_{h_{0}}$. Moreover, $\tau_{0}=(\theta_{0},h_{0})$ is the unique root of ${\operatorname{E}\left[D_{\theta}\sum_{l=1}^{m}\mathfrak{L}_{n_{l}}(\theta,h)\right]=0}$ and ${\operatorname{E}\left[D_{h}\sum_{l=1}^{m}\mathfrak{L}_{n_{l}}(\theta,h)\right]=0}$ in $\mathcal{N}_{\theta_{0}}\times\mathcal{N}_{h_{0}}$.

Assumption 4**.**

For any $\theta_{*}\in\mathcal{N}_{\theta_{0}},h_{*}\in\mathcal{N}_{h_{0}},l\in\{1,\ldots,m\}$, $D_{\theta}\eta_{l}(\theta_{*},h_{*})$ is a bounded linear operator from $\mathbb{R}^{p}$ to $L^{2}_{0}(\mathcal{X}_{l})$, $D_{h}\eta_{l}(\theta_{*},h_{*})$ is a bounded linear operator from ${\cal H}$ to $L^{2}_{0}(\mathcal{X}_{l})$, and there exist positive constants $C_{1},C_{2}$, such that for all $(\zeta,g)\in\mathcal{Q}$,

[TABLE]

Assumption 5**.**

For any $(\theta_{1},h_{1}),(\theta_{2},h_{2})\in\mathcal{N}_{\theta_{0}}\times\mathcal{N}_{h_{0}}$, there exists a positive constant $C_{d}<2C_{1}$, where $C_{1}$ is as defined in Assumption 4, such that for any $(\zeta,g)\in\mathcal{N}_{\theta_{0}}\times\mathcal{N}_{h_{0}}$,

[TABLE]

Assumption 6**.**

$\mathcal{N}_{\theta_{0}}\times\mathcal{N}_{h_{0}}$* is a convex set that contains $\hat{\tau}$, and there exist $C_{3},C_{4}>0$ such that*

[TABLE]

holds uniformly for any $\tau=(\theta,h)\in\mathcal{N}_{\theta_{0}}\times\mathcal{N}_{h_{0}}$ and $f\in L^{2}_{0}(\mathcal{X}_{l})$.

Remark 2**.**

We compare the assumptions above to the existing literature on the single sample ($m=1$) nonparametric model and the semiparametric additive model.

1. (i)

For the nonparametric model where $\eta(x;\theta,h)=h(x)$, and for the semiparametric additive model where $\eta(x;\theta,h)=\alpha(x;\theta)+h(x)$, we see that ${\cal H}\subset L_{0}^{2}(\mathcal{X})$ and $L_{h}$ can be chosen to be the inclusion operator $\iota:{\cal H}\to L_{0}^{2}(\mathcal{X})$. It is easy to see that $V_{\tau_{0}}(g_{1},g_{2})$ defines an inner product on $L_{0}^{2}(\mathcal{X})$. Moreover, the norm $\left\lVert g\right\rVert_{{\cal H}}=[V_{\tau_{0}}(g,g)+\lambda J(g,g)]^{1/2}$ and its analog for regression models have been widely used in the smoothing spline literature for nonparametric models, including \citeasnounsilverman_1982, \citeasnouncox_1988, \citeasnoungu_qiu_1993, \citeasnounshang2010convergence, etc. 2. (ii)

Consider the linear additive model $\eta(x;\theta,h)=\theta^{T}a(x)+h(x)$, where $a(x)=[a^{k}(x)]_{k=1}^{p}$ is a vector of bounded $L^{2}_{0}(\mathcal{X})$ functions. One may choose $L_{\theta}=D_{\theta}\eta(\theta,h)=a(x)$ and $L_{h}=D_{h}\eta(\theta,h)=\iota$, where $\iota$ is the inclusion operator from ${\cal H}$ to $L^{2}_{0}(\mathcal{X})$. Since $V_{\tau_{0}}(\theta^{T}a-h)$ measures the variance of the difference between the parametric and nonparametric components, identifiability of this model follows from Assumption 2(ii). 3. (iii)

As discussed above, by choosing the appropriate operators $L_{\theta}$ and $L_{h}$ for either the nonparametric model or the semiparametric linear additive model, our analysis is reduced to that studied by \citeasnoungu_qiu_1993 or \citeasnounyang2009penalized, respectively. In both cases, due to the linearity of $\eta$, the neighborhood $\mathcal{N}_{\theta}\times\mathcal{N}_{h}$ can be taken to be the whole parameter space. Furthermore, Assumptions 3 to 5 are satisfied automatically. One may see that such assumptions are simply redundant when $\eta$ is linear in $\theta$ and $h$. 4. (iv)

Note that Assumption 6 is similar to Assumption A.3 in \citeasnoungu_qiu_1993, Assumption A.3 in \citeasnoungu1995smoothing, and Condition 3 in \citeasnounyang2009penalized. As mentioned by \citeasnoungu1995smoothing, $V_{\tau}(f-g)$ can be viewed as a kind of weighted mean square error between functions $f$ and $g$ with weight function $\exp\{\eta(x;\tau)\}$. Since $\eta$ is continuous in $\mathcal{N}_{\theta_{0}}\times\mathcal{N}_{h_{0}}$, a small change in $\tau$ yields a small change in the weight function, and Assumption 6 simply guarantees that this also yields a relatively small change in the weighted mean square error.

4.2.3 Spectral decomposition

We now construct an eigensystem for the functionals $g\mapsto\sum_{l=1}^{m}V_{l,\tau_{0}}(L_{l,h}g)$ and $J$.

Assumption 7**.**

The quadratic functional $g\mapsto\sum_{l=1}^{m}V_{l,\tau_{0}}(L_{l,h}g)$ is completely continuous with respect to the quadratic functional $J$.

Under Assumption 7, Theorem 3.1 of \citeasnounweinberger_1974 yields a sequence $\{\phi_{\nu}:\nu=1,2,\ldots\}$ of eigenfunctions and a sequence $\{\rho_{\nu}:\nu=1,2,\ldots\}$ of eigenvalues such that

[TABLE]

for all pairs $\nu,\mu$ of positive integers, where $\delta_{\nu\mu}$ is the Kronecker delta and $0\leq\rho_{\nu}\rightarrow\infty$.

Assumption 8**.**

$\rho_{\nu}=\kappa_{\nu}\nu^{r}$, where $r>1$ and $\kappa_{\nu}\in(\beta_{1},\beta_{2})\subset(0,\infty)$.

By Assumptions 4 and 7, for any $(\theta_{*},h_{*})\in\mathcal{N}_{\theta_{0}}\times\mathcal{N}_{h_{0}}$, there exist sequences of eigenfunctions $\{\phi_{*,\nu}:\nu=1,2,\ldots\}$ and eigenvalues $\{\rho_{*,\nu}:\nu=1,2,\ldots\}$ such that

[TABLE]

for all pairs $\nu,\mu$ of positive integers, where $0\leq\rho_{*,\nu}\rightarrow\infty$. Assumption 4 implies that there exist positive constants $c_{1},c_{2}$ such that for $\nu$ large enough, $c_{1}\rho_{\nu}\leq\rho_{*,\nu}\leq c_{2}\rho_{\nu}$. By Assumption 8, $\rho_{*,\nu}\sim\nu^{r}$ for large enough $\nu$, where $r>1$.

Note that for every $h\in{\cal H}$ and any $(\theta_{*},h_{*})\in\mathcal{N}_{\theta_{0}}\times\mathcal{N}_{h_{0}}$, we can have a Fourier expansion

[TABLE]

4.3 Linear approximation

For this section and the next, we present our analysis and provide the rate of convergence for the single sample case when $m=1$.

Consider the eigensystem discussed in Section 4.2.3 with $(\theta_{*},h_{*})=(\theta_{0},h_{0})$. We have the Fourier series expansions $h=\sum_{\nu}h_{\nu}\phi_{0,\nu}$ and $h_{0}=\sum_{\nu}h_{0,\nu}\phi_{0,\nu}$, where

[TABLE]

are the Fourier coefficients of $h$ and $h_{0}$ with respect to the base $\phi_{0,\nu}$. Plugging these into equation (15), we get

[TABLE]

In equation (19), using the fact that $D_{\theta}\eta(x;\tau_{0})$ can be represented by a $p\times 1$ vector whose $k$th entry is $D_{\theta^{k}}\eta(x;\tau_{0})$, we denote by $D_{\theta}\eta(x;\tau_{0})$ both the linear operator and its vector form, i.e., $D_{\theta}\eta(x;\tau_{0})\theta=\theta^{T}D_{\theta}\eta(x;\tau_{0})$. We also have $\mu_{\tau_{0}}[D_{\theta}\eta(\tau_{0})]=[\mu_{\tau_{0}}[D_{\theta^{1}}\eta(\tau_{0})],\ldots,\mu_{\tau_{0}}[D_{\theta^{p}}\eta(\tau_{0})]]^{T}$, $V_{\tau_{0}}[D_{\theta}\eta(\tau_{0})]$ is the $p\times p$ covariance matrix whose $(i,j)$th entry is $V_{\tau_{0}}[D_{\theta^{i}}\eta(\tau_{0}),D_{\theta^{j}}\eta(\tau_{0})]$, and $V_{\tau_{0}}[D_{\theta}\eta(\tau_{0}),D_{h}\eta(\tau_{0})\phi_{0,\nu}]$ is the $p\times 1$ vector whose $i$th entry is $V_{\tau_{0}}[D_{\theta^{i}}\eta(\tau_{0}),D_{h}\eta(\tau_{0})\phi_{0,\nu}]$. Let

[TABLE]

The Fourier coefficients $\tilde{h}_{\nu}$ and $\tilde{\theta}$ that minimize equation (19) are therefore given by setting the partial derivatives of equation (19) to 0 and solving the resulting system. We get

[TABLE]

where

[TABLE]

and $a^{\otimes 2}=aa^{T}$ ($a$ is a vector or matrix). Therefore, $(\tilde{\theta},\tilde{h})$, where $\tilde{h}=\sum_{\nu}\tilde{h}_{\nu}\phi_{0,\nu}$, is the minimizer of equation (15). Using the terminology from the nonparametric setting [Gu2013bk], $\tilde{h}$ can be called a linear approximation of $\hat{h}$ since $\tilde{h}_{\nu}$ is linear in $\phi_{0,\nu}(X_{i})$ given $\tilde{\theta}$. We will slightly abuse the terminology for our nonlinear case here and also call $\tilde{\tau}=(\tilde{\theta},\tilde{h})$ the linear approximation of $\hat{\tau}=(\hat{\theta},\hat{h})$.

Using the fact that $\operatorname{E}(\bar{\alpha}_{n})=\operatorname{E}(\bar{\beta}_{\nu,n})=0$, $\operatorname{E}(\|\bar{\alpha}_{n}\|_{l^{2}}^{2})=O(n^{-1})$, $\operatorname{E}(\bar{\beta}_{\nu,n}^{2})=n^{-1}$, and after tedious calculation, we get the following lemma and theorem, whose proofs are given in the Appendix.

Lemma 1**.**

Under Assumptions 1 to 3, 7, and 8, as $n\rightarrow\infty$ and $\lambda\rightarrow 0$,

[TABLE]

Theorem 1**.**

Under Assumptions 1 to 4, 7, and 8, as $n\rightarrow\infty$ and $\lambda\rightarrow 0$,

[TABLE]

Note that in addition to $n\rightarrow\infty$ and $\lambda\rightarrow 0$, if also $n^{-1}\lambda^{-\frac{1}{r}}\rightarrow 0$, then the probability that $\tilde{\tau}\in\mathcal{N}_{\theta_{0}}\times\mathcal{N}_{h_{0}}$ tends to $1$ (since $\tilde{\tau}\rightarrow\tau_{0}$ under these conditions). We can restrict our attention to this event for the rest of our analysis, or for simplicity, we assume that $\tilde{\tau}\in\mathcal{N}_{\theta_{0}}\times\mathcal{N}_{h_{0}}$.

4.4 Approximation error and main results for case $m=1$

In this section, we first find a bound for the approximation error $\hat{\tau}-\tilde{\tau}=(\hat{\theta}-\tilde{\theta},\hat{h}-\tilde{h})$, which will then imply the convergence of $\hat{\tau}-\tau_{0}=(\hat{\theta}-\theta_{0},\hat{h}-h_{0})$.

Define

[TABLE]

By straightforward calculations, we have

[TABLE]

Equation (20) (or equation (21)) can be understood as the Fréchet partial derivative of $\mathfrak{L}_{n,\lambda}(\theta,h)$ with respect to $\theta$ (or $h$) at $\tau_{1}=(\theta_{1},h_{1})$ in the direction of $\theta_{2}$ (or $h_{2}$). Likewise, equation (22) (or equation (23)) is the Fréchet partial derivative of $\tilde{\mathfrak{L}}_{n,\lambda}(\theta,h)$ with respect to $\theta$ (or $h$) at $\tau_{1}=(\theta_{1},h_{1})$ in the direction of $\theta_{2}$ (or $h_{2}$). Therefore, for any $\tau_{2}$, we note that $\dot{A}^{\theta}_{\hat{\tau},\tau_{2}}(0)=\dot{A}^{h}_{\hat{\tau},\tau_{2}}(0)=0$, and $\dot{B}^{\theta}_{\tilde{\tau},\tau_{2}}(0)=\dot{B}^{h}_{\tilde{\tau},\tau_{2}}(0)=0$.

Setting $\tau_{1}=\hat{\tau}$, $\tau_{2}=\hat{\tau}-\tilde{\tau}$ in (20) + (21), we get

[TABLE]

and setting $\tau_{1}=\tilde{\tau}$, $\tau_{2}=\hat{\tau}-\tilde{\tau}$ in (22) + (23) implies

[TABLE]

Equating (24) and (25), rearranging terms, and subtracting ${\mu_{\tilde{\tau}}[D_{\theta}\eta(\hat{\tau})(\hat{\theta}-\tilde{\theta})+D_{h}\eta(\hat{\tau})(\hat{h}-\tilde{h})]}$ on both sides, we have

[TABLE]

Note that for any function $f\in L^{2}_{0}(\mathcal{X})$ and $\tau_{1},\tau_{2}\in\mathcal{N}_{\theta_{0}}\times\mathcal{N}_{h_{0}}$, it is easy to show using the mean value theorem that

[TABLE]

for some $0\leq\alpha\leq 1$. We are now ready to prove the following theorem.

Theorem 2**.**

Under Assumptions 1 to 8, as $n\rightarrow\infty$, $\lambda\rightarrow 0$, and $n\lambda^{\frac{1}{r}}\rightarrow\infty$,

[TABLE]

Therefore,

[TABLE]

Proof.

We first obtain a lower bound for the first two terms in the LHS of (26) and an upper bound for the second term in the RHS of (26). For some $0\leq\alpha\leq 1$, let $\tau_{*}=\tilde{\tau}+\alpha(\hat{\tau}-\tilde{\tau})$. By Assumptions 4, 5, and 6, and equation (27), the first two terms in the LHS of (26) yield

[TABLE]

for some $0\leq c_{3}\leq C_{3}(2C_{1}-C_{d})$. For the second term in the RHS of (26), let $\tau_{u}=\tilde{\tau}+u(\tau_{0}-\tilde{\tau})$ for some $0\leq u\leq 1$. By equation (27), the Cauchy-Schwarz inequality, and Assumptions 4 and 6, we have

[TABLE]

Next, for any random sample $X_{1},\ldots,X_{n}$ and any function $f\in L^{2}_{0}(\mathcal{X})$, it is well known that $\operatorname{E}\left\{[n^{-1}\sum_{i=1}^{n}f(X_{i})-\operatorname{E}(f(X))]^{2}\right\}=n^{-1}\operatorname{Var}[f(X)]$, and this implies $n^{-1}\sum_{i=1}^{n}f(X_{i})-\operatorname{E}[f(X)]=O_{p}(n^{-\frac{1}{2}})\operatorname{Var}^{\frac{1}{2}}[f(X)]$. Thus, the last two terms in the RHS of (26) are

[TABLE]

respectively.

Putting everything together, as $\lambda\rightarrow 0$, $n\rightarrow\infty$, and $n\lambda^{\frac{1}{r}}\rightarrow\infty$, by all assumptions and Theorem 1, we get

[TABLE]

The result then follows after trivial manipulation. ∎

The following corollaries as direct results from Theorem 2.

Corollary 1**.**

Under Assumptions 1 to 8, as $n\rightarrow\infty$, $\lambda\rightarrow 0$, and $n\lambda^{\frac{1}{r}}\rightarrow\infty$, we have

[TABLE]

Corollary 2**.**

If $\eta(\theta,h)=\alpha(x;\theta)+h(x)$, $L_{h}$ can be chosen to be the inclusion operator from ${\cal H}$ to $L_{0}^{2}(\mathcal{X})$. Under the same conditions as in Theorem 2, we have

[TABLE]

Corollaries 1 and 2 hold since $\left\lVert\cdot\right\rVert_{\mathbb{R}^{p}}$ is equivalent to the Euclidean norm $\left\lVert\cdot\right\rVert_{l^{2}}$ on $\mathbb{R}^{p}$, and $\left\lVert\cdot\right\rVert_{\mathcal{Q}}$ is equivalent to the product norm $\left\lVert\cdot\right\rVert_{\mathcal{Q},1}$ on the joint parameter space $\mathcal{Q}$ (see Remark 1$(ii)$ and $(iv)$). We now derive a convergence rate of the overall density function, measured by the symmetrized Kullback-Leibler distance, defined in (16).

Theorem 3**.**

Under Assumptions 1 to 8, as $n\rightarrow\infty$, $\lambda\rightarrow 0$, and $n\lambda^{\frac{1}{r}}\rightarrow\infty$,

[TABLE]

Proof.

By the definition of the Fréchet derivative of $\eta(\theta,h)$ at $(\theta_{0},h_{0})$, we have

[TABLE]

By Theorem 2, $\|(\hat{\theta}-\theta_{0},\hat{h}-h_{0})\|_{\mathcal{Q}}\rightarrow 0$ with probability tending to 1 as $\lambda\rightarrow 0$, and $n\lambda^{\frac{1}{r}}\rightarrow\infty$. Therefore,

[TABLE]

By equation (27), for some $0\leq u\leq 1$, we have

[TABLE]

By the Cauchy-Schwarz inequality, Assumptions 4 and 6, and Theorem 2, (28) is bounded above by

[TABLE]

∎

Remark 3**.**

Note that results in Corollary 1, Corollary 2, and Theorem 3 are independent of the linear operators $L_{\theta}$ and $L_{h}$.

4.5 Joint asymptotics for case $m\geq 1$

In this section, we will discuss the joint consistency for $m\geq 1$. The proof for general $m$ is a simple extension of the proof for $m=1$ discussed in Sections 4.3 and 4.4.

The linear approximation estimate $\tilde{\tau}=(\tilde{\theta},\tilde{h})$ of $\hat{\tau}=(\hat{\theta},\hat{h})$, where $\tilde{h}=\sum_{\nu}\tilde{h}_{\nu}\phi_{0,\nu}$, can be obtained by minimizing equation (15). We get

[TABLE]

where

[TABLE]

Recall that $n=\min_{1\leq l\leq m}\{n_{l}\}$. Since

[TABLE]

one can get the following results, which are similar to Lemma 1 and Theorems 1 – 3. Note that one also has corollaries to Theorem 5 which are analogous to Corollaries 1 and 2, but for $m>1$. They are omitted to save space.

Lemma 2**.**

Under Assumptions 1 to 3, 7, and 8, as ${n=\min_{1\leq l\leq m}\{n_{l}\}\rightarrow\infty}$ and $\lambda\rightarrow 0$,

[TABLE]

Theorem 4**.**

Under Assumptions 1 to 4, 7, and 8, as ${n=\min_{1\leq l\leq m}\{n_{l}\}\rightarrow\infty}$ and $\lambda\rightarrow 0$,

[TABLE]

Theorem 5**.**

Under Assumptions 1 to 8, as ${n=\min_{1\leq l\leq m}\{n_{l}\}\rightarrow\infty}$, $\lambda\rightarrow 0$, and $n\lambda^{\frac{1}{r}}\rightarrow\infty$,

[TABLE]

Therefore,

[TABLE]

Theorem 6**.**

Under Assumptions 1 to 8, as ${n=\min_{1\leq l\leq m}\{n_{l}\}\rightarrow\infty}$, $\lambda\rightarrow 0$, and $n\lambda^{\frac{1}{r}}\rightarrow\infty$,

[TABLE]

5 Simulations

To evaluate the proposed estimation procedures and algorithms, we conduct several simulations for both additive and nonadditive cases. We use the KL divergence between the true density and the estimated density to evaluate the performance of density estimation, $\text{KL}(f,\hat{f})=\int_{\cal X}f(x)\log\{f(x)/\hat{f}(x)\}dx$. Furthermore, we will use the generalized decomposition

[TABLE]

proposed by \citeasnounheskes1998bias to evaluate the bias-variance trade-off where $\bar{f}=\exp\{\hbox{E}(\log\hat{f})\}/Z$ and $Z$ is a normalization constant.

We emphasize that in all the following tables, all numbers are in $10^{-2}$. The numbers in the parentheses for the KL measures are the biases and the variances in the decomposition of the KL distances. The numbers in the parentheses for the mean squared errors (MSEs) are the squared biases and the variances of the parameter estimators.

5.1 Additive cases

5.1.1 Near Normal distribution

Consider the density function

[TABLE]

where $x\in[0,1]$, $\alpha(x;\mu,\sigma)=-\frac{x^{2}}{2\sigma^{2}}+\frac{\mu x}{\sigma^{2}}$, $h(x)=ax^{3}$ and $a$ is a constant. The function $\alpha(x;\mu,\sigma)$ is the logistic transformation of the truncated Normal density function with mean $\mu$ and standard deviation $\sigma$. The constant $a$ controls the departure from the truncated Normal distribution. The choice of this density function is motivated by \citeasnounhjort1995nonparametric, with a truncated Normal as the starting parametric density. We will consider the additive model (4) with $\alpha(x;\mu,\sigma)=-\frac{x^{2}}{2\sigma^{2}}+\frac{\mu x}{\sigma^{2}}$ and $h\in W_{2}^{3}[0,1]\ominus\{1,x,x^{2}\}$, where $W_{2}^{s}[0,1]$ is the $s$th order Sobolev space on ${\cal X}=[0,1]$ defined as

[TABLE]

We consider the third order Sobolev space because the null space with quadratic functions corresponds to the probability density functions of the truncated Normal distributions. This null space has been removed from the model space of $h$ to make the model identifiable. We use the profile likelihood approach in Section 3.1 to compute estimates of $\theta=(\mu,\sigma)$ and $h$.

For comparison, we also implemented the nonparametric estimation methods based on kernels and polynomial splines. Given a random sample $X_{1},\ldots,X_{n}$ that has density $f(x)$ on $[0,1]$, the polynomial spline based nonparametric estimation method estimates the logistic transformation of $f$, $\eta$, by minimizing the penalized likelihood

[TABLE]

over $\eta\in W_{2}^{s}[0,1]$. We implemented the cubic spline ($s=2$) and the quintic spline ($s=3$). In the implementation of kernel density estimation and the method by \citeasnounhjort1995nonparametric, we use the Gaussian kernel and compare different methods for bandwidth selection, including the approach in \citeasnounscott1992curse, the unbiased and biased cross-validation procedures, and the method using pilot estimation of derivatives in \citeasnounsheather1991reliable. We find that the bandwidth obtained by

[TABLE]

provides the best overall performance. We report results with bandwidth determined by (30) only.

We set $\mu=0.5$ and $\sigma=0.2$. We consider three choices of $a$: $a=0.25,1,4$, and three sample sizes, $n=100,200$ and $500$. For each simulation setting, we generate 100 simulation data sets. For each simulated data set, we compute density estimates using the kernel method, the method by \citeasnounhjort1995nonparametric (HG), cubic spline, quintic spline and our proposed semiparametric (SEMI) method. Table 1 summarizes the performances of these methods in terms of the density estimation. As expected, the performance improves as sample size increases, and the semiparametric and quintic spline estimates have the smallest KL than other nonparametric estimates. Our semiparametric approach performs better than the semiparametric approach in \citeasnounhjort1995nonparametric. There is virtually no difference between the quintic spline approach. This is expected since they are fitting the same model using different methods.

5.1.2 Near Gumbel distribution

Consider the density function

[TABLE]

where $x\in{\cal X}=[0,1]$, $\alpha(x;\mu,\sigma)=-{(x-\mu)}/{\sigma}-\exp({(x-\mu)}/{\sigma})$, $h(x)=ax^{3}$, and $a$ is a constant. This is equivalent to start with a truncated Gumbel distribution whose logistic transformation is $\alpha(x;\mu,\sigma)$, and the constant $a$ controls the departure from the truncated Gumbel distribution. We will consider the additive model (4) with $\alpha(x;\mu,\sigma)=-{(x-\mu)}/{\sigma}-\exp({(x-\mu)}/{\sigma})$ and $h\in W_{2}^{2}[0,1]\ominus\{1\}$. Note that $\alpha$ is non-linear in $\theta$. Therefore, the method in \citeasnounyang2009penalized does not apply. We compare our proposed method with the kernel, cubic spline and HG’s method, where $f_{0}$ in the HG’s approach is a truncated Gumbel distribution. We use the profile likelihood approach as described in Section 3.1 to compute estimates of $\theta=(\mu,\sigma)$ and $h$.

We set $\mu=0.5$ and $\sigma=0.2$. We consider three choices of $a$: $a=0.25,1,4$, and three sample sizes $n=100,200,500$. For each simulation setting, we generate 100 simulated data sets. Table 2 shows that the semiparametric methods has smaller KL distances when the true density is close to truncated Gumbel (i.e. $a=0.25$ and $1$). When the true density if far away from the truncated Gumbel ($a=4$), cubic spline has smaller KL distances.

5.2 Nonadditive cases

5.2.1 Power transformation

The truncated Weibull distribution has density function

[TABLE]

where $x\in{\cal X}=[0,1]$, $s>0$ is the scale parameter and $\gamma>0$ is the shape parameter. Since $Y=X^{\gamma}\sim\text{truncated Exp}(s)$, we consider the transformation model in Section 3.3 with

[TABLE]

and $h\in W_{2}^{2}[0,1]\ominus\{1\}$. The algorithm in Section 3.3 is used to compute the estimates of $\theta$ and $h$. We also compute the cubic spline estimate of the density function (Cubic) and the kernel density estimation after transformation (TransKernel) [wand1991transformations], where the transformation form is known as in (32) up to the parameter $\theta$. In the simulations, we set $s=1$, and consider $\gamma=1,2,3$ where $\gamma=1$ corresponding to the truncated exponential distribution for which the cubic spline is well-suited since the null space corresponds to the truncated exponential density. We consider two sample sizes $n=100$ and $200$. For each simulation setting, we generate 100 data sets. The simulation results are shown in Table 3. As expected, cubic spline performs best when $\gamma=1$. Our semiparametric approach provides smaller KL distances in other cases. As sample size increases, both the performances of density and parameter estimation improved, except for the transformation kernel density estimation, the biases of the parameter estimates actually increased although the variances reduced, leading to slightly larger MSEs.

5.2.2 Two-sample density estimation

The Gumbel and logistic distributions are members of location scale family with the location parameter $\mu$ and scale parameter $\sigma$. They have been extensively used in many different areas.

We generate two independent samples, $X_{1},\ldots,X_{n_{1}}\overset{iid}{\sim}f(x;0,1)$ and $Y_{1},\ldots,Y_{n_{2}}\overset{iid}{\sim}f(x;\mu,\sigma)$ where $f$ is either the Gumbel distribution or the logistic distribution. We set $\mu=2$, $\sigma=1$, and consider four sample sizes $(n_{1},n_{2})=(100,100),(100,200),(200,100)$ and $(200,200)$. We fit model (13) with $\alpha_{1}(x;\theta)=x$ and $\alpha_{2}(x;\theta)=\frac{x-\mu}{\sigma}$ and $h$ belongs to the RKHS for the univariate thin-plate splines with ${\cal X}=\mathbb{R}$ and

[TABLE]

We estimate the density functions using the procedure described in Section 3.3, and report KL$(f(\cdot;0,1),\hat{f}(\cdot;0,1))+\text{KL}(f(\cdot;\mu,\sigma),\hat{f}(\cdot;\hat{\mu},\hat{\sigma}))$ for the overall performance of density estimation (overall KL). For comparison, we consider the characteristic function based two sample transformation density estimation method (CHAR) in \citeasnounpotgieter2012nonparametric, and another method which estimates the densities for $X_{i}$’s and $Y_{i}$’s using thin-plate spline models separately (TP), with logistic of densities belong to ${\cal H}$ in (33). Denote the separate thin-plate estimates for $X$ and $Y$ samples as $\hat{f}_{1}$ and $\hat{f}_{2}$ respectively. We report KL$(f(\cdot;0,1),\hat{f}_{1})+\text{KL}(f(x;\mu,\sigma),\hat{f}_{2})$ for the overall performance (overall KL) of the density estimation. For the method CHAR and our method SEMI, we report the MSEs, squared biases and variances as measures for the performance of parameter estimations. We report $\hbox{KL}_{s}=\hbox{KL}(f(\cdot;0,1),\exp(\hat{h})/\int\exp(\hat{h})dx)$ as a measure of the performance for estimating the function $h$.

Gumbel distribution The Gumbel distribution has the density $f(x)={\frac{1}{\sigma}}e^{{-({(x-\mu)}/{\sigma}}+e^{{-{(x-\mu)}/{\sigma}}})}$, $\quad x\in\mathbb{R}$. Table 4 shows although that the CHAR method estimates the parameters slightly better, our semiparametric method provides better density estimation than the method in \citeasnounpotgieter2012nonparametric and the separate thin-spline estimates. All methods improved as the sample size increases.

Logistic distribution The pdf of the logistic distribution is given by $f(x;\mu,\sigma)={\frac{e^{-{{(x-\mu)}/{\sigma}}}}{\sigma\left(1+e^{-{{(x-\mu)}/{\sigma}}}\right)^{2}}}$, $\quad x\in\mathbb{R}.$ The simulation results are in Table 5. We observe the same relative behavior as in the case of two-sample Gumbel distributions.

6 Application to the Old Faithful data

The data set consists of waiting time between eruptions and the duration of the eruption gathered from 272 consecutive eruptions of the Old Faithful geyser in Yellowstone National Park, Wyoming, USA [hardle2012smoothing]. Histograms of these two variables suggest that both of the two variables have two modes. Therefore, we consider a semiparametric model with a mixture Normal as the parametric component:

[TABLE]

where $h\in W_{2}^{2}[40,100]\ominus\{1\}$ for the waiting time and $h\in W_{2}^{2}[1.5,5.5]\ominus\{1\}$ for the duration time, and $f_{0}$ is a mixture Normal density

[TABLE]

Model (34) is a special case of the additive model (4). We compute density estimates for both waiting time and duration using the profiled likelihood estimation procedure in Section 3.1. We also compute the estimates of the mixture Normal density (35) and cubic spline for both the waiting time and the duration.

The estimated densities are shown in Figures 1 and 2 for the waiting time and the duration respectively. For the waiting time, the estimates from the semiparametric model and mixture Normal are quite close suggesting that the parametric mixture Normal fits data adequately. Both of them follow the data on the left slightly better than the cubic spline estimate. For the duration, the estimates from the semiparametric model and cubic spline are close. And both of them portray both modes more precisely.

Acknowledgements

We gratefully acknowledge support from the National Science Foundation (DMS-1507078 for Anna Liu and DMS-1507620 for Yuedong Wang). We acknowledge support from the Center for Scientific Computing from the CNSI, MRL: an NSF MRSEC (DMR-1720256).

Appendix

Verification of inner product $\left\langle\cdot,\cdot\right\rangle_{\mathcal{Q}}$

We now turn to the discussion of the validity of (18). For any $p_{1}\times 1$ and $p_{2}\times 1$ vectors of functions in ${\cal H}$, say $G_{1}=[G_{1}^{k}]_{k=1}^{p_{1}}$ and $G_{2}=[G_{2}^{k}]_{k=1}^{p_{2}}$, we use the vector form of the inner product $\langle G_{1},G_{2}\rangle_{{\cal H}}$ to denote a $p_{1}\times p_{2}$ matrix in which the $(i,j)$th entry is $\langle G_{1}^{i},G_{2}^{j}\rangle_{{\cal H}}$. For any $g\in{\cal H}$, let $\mathcal{F}^{k}g=\sum_{l=1}^{m}V_{l,\tau_{0}}[L_{\l,\theta}^{k},L_{l,h}g]$. Since

[TABLE]

for some positive constant $C$, $\mathcal{F}^{k}$ is a bounded linear functional on ${\cal H}$. By the Riesz representation theorem, there exists a $F^{k}\in{\cal H}$ such that for any $g\in{\cal H}$,

[TABLE]

Let $F=\left[F^{k}\right]_{k=1}^{p}$. We define $V_{l,\tau_{0}}(L_{l,\theta},L_{l,h}g)$, $V_{l,\tau_{0}}(L_{l,h}F,L_{l,h}g)$, and $J(F,g)$ to be $p\times 1$ vectors whose $k$th entries are $V_{l,\tau_{0}}(L_{l,\theta}^{k},L_{l,h}g)$, $V_{l,\tau_{0}}(L_{l,h}F^{k},L_{l,h}g)$, and $J(F^{k},g)$, respectively. Therefore,

[TABLE]

We also define the $p\times p$ matrix $\Omega_{F}=\sum_{l=1}^{m}V_{l,\tau_{0}}(L_{l,\theta}-L_{l,h}F,L_{l,\theta}-L_{l,h}F)$, whose $(i,j)$th entry is $\sum_{l=1}^{m}V_{l,\tau_{0}}(L^{i}_{l,\theta}-L_{l,h}F^{i},L^{j}_{l,\theta}-L_{l,h}F^{j})$.

Lemma A.1.

Under Assumption 2, $\Omega_{F}$ is positive definite and the eigenvalues of $\Omega_{F}$ are greater than $c_{\delta}$, which is as defined in Assumption 2$(ii)$.

Proof.

Fix a non-zero vector $\zeta\in\mathbb{R}^{p}$ and write $\zeta_{*}=\zeta/\left\lVert\zeta\right\rVert_{l^{2}}$. We have

[TABLE]

where the last inequality holds by Assumption 2$(ii)$ because $\left\lVert\zeta_{*}\right\rVert_{l^{2}}=1$ and $\zeta_{*}^{T}F\in{\cal H}$. Therefore, $\Omega_{F}$ is positive definite.

Let $\delta$ be any eigenvalue of $\Omega_{F}$, and let $\zeta_{\delta}\in\mathbb{R}^{p}$ be a unit eigenvector associated with $\delta$. By definition, we have $\Omega_{F}\zeta_{\delta}=\delta\zeta_{\delta}$. We have

[TABLE]

∎

Theorem A.1.

Suppose Assumption 2 holds. Then $\langle\cdot,\cdot\rangle_{\mathcal{Q}}$ given by (18) is a well-defined inner product on $\mathcal{Q}$, and $\mathcal{Q}$ is complete with respect to the norm $\left\lVert\cdot\right\rVert_{\mathcal{Q}}$ induced by this inner product. Hence, $\mathcal{Q}$ is a Hilbert space.

Proof.

It is easy to check that (18) satisfies symmetry, linearity and positive semi-definiteness for an inner product. If $(\zeta,g)=0$, $\langle(\zeta,g),(\zeta,g)\rangle_{\mathcal{Q}}=0$ is obvious. We will now show that $\langle(\zeta,g),(\zeta,g)\rangle_{\mathcal{Q}}=0$ implies $(\zeta,g)=0$. We see that

[TABLE]

and every term in (36) is non-negative. If $\langle(\zeta,g),(\zeta,g)\rangle=0$, the first term in (36) implies $\zeta=0$ by Lemma 1. This further implies that

[TABLE]

Therefore, $g=0$ because $\langle\cdot,\cdot\rangle_{{\cal H}}$ is an inner product on ${\cal H}$. Hence, $\langle\cdot,\cdot\rangle_{\mathcal{Q}}$ is a well-defined inner product on $\mathcal{Q}$.

Next, we want to show that $\mathcal{Q}$ is complete with respect to the norm $\left\lVert\cdot\right\rVert_{\mathcal{Q}}$. Let $\{(\zeta_{i},g_{i})\}_{i=1}^{\infty}\subset\mathcal{Q}$ be a Cauchy sequence. For any $\epsilon>0$, there exist a positive integer $M$ such that for all $i,j>M$, we have

[TABLE]

This implies that

[TABLE]

where $(\zeta_{i}-\zeta_{j})^{*}=(\zeta_{i}-\zeta_{j})/\left\lVert\zeta_{i}-\zeta_{j}\right\rVert_{l^{2}}$, and $(g_{i}-g_{j})^{*}=(g_{i}-g_{j})/\left\lVert\zeta_{i}-\zeta_{j}\right\rVert_{l^{2}}$. By Assumption 2$(ii)$, $\sum_{l=1}^{m}V_{l,\tau_{0}}[L_{l,\theta}(\zeta_{i}-\zeta_{j})^{*}+L_{l,h}(g_{i}-g_{j})^{*}]>c_{\delta}$ for some positive constant $c_{\delta}$ as defined in (17). Therefore,

[TABLE]

and $\{\zeta_{i}\}_{i=1}^{\infty}$ is a Cauchy sequence in $\mathbb{R}^{p}$ under the Euclidean norm, which therefore converges to some limit $\zeta_{\infty}\in\mathbb{R}^{p}$.

To find a limit for the sequence $\{g_{i}\}_{i=1}^{\infty}$ in ${\cal H}$, we consider

[TABLE]

where the first inequality follows from the Cauchy-Schwarz inequality, and the second inequality follows from the triangle inequality. For $a>0,b>0$, we have $(1/4)a+b-a^{1/2}b^{1/2}=[(1/2)a^{1/2}-b^{1/2}]^{2}\geq 0$, and it follows that $2a^{1/2}b^{1/2}\leq(1/2)a+2b$. For $a=\left\lVert g_{i}-g_{j}\right\rVert_{{\cal H}}^{2}$ and $b=\left\lVert\zeta_{i}-\zeta_{j}\right\rVert_{\mathbb{R}^{p}}^{2}$, (38) becomes

[TABLE]

Since $\left\lVert\cdot\right\rVert_{\mathbb{R}^{p}}$ is equivalent to $\left\lVert\cdot\right\rVert_{l^{2}}$ on $\mathbb{R}^{p}$, $\left\lVert\zeta_{i}-\zeta_{j}\right\rVert_{\mathbb{R}^{p}}^{2}\leq C\epsilon$ for some positive constant $C$ by (37). Therefore, after rearranging (39), we get $\left\lVert g_{i}-g_{j}\right\rVert_{{\cal H}}^{2}\leq(2+C)\epsilon$. Hence, $\{g_{i}\}_{i=1}^{\infty}$ is a Cauchy sequence in ${\cal H}$ under the norm $\left\lVert\cdot\right\rVert_{{\cal H}}$. By Assumption 2$(i)$, this sequence converges to some limit $g_{\infty}\in{\cal H}$.

Lastly, we show that $(\zeta_{i},g_{i})$ converges to $(\zeta_{\infty},g_{\infty})$ in $\left\lVert\cdot\right\rVert_{\mathcal{Q}}$. By the Cauchy-Schwarz inequality and triangle inequality, as $i\to\infty$, we have

[TABLE]

Therefore, we conclude that $\mathcal{Q}$ is a Hilbert space with respect to the inner product $\langle\cdot,\cdot\rangle_{\mathcal{Q}}$. ∎

Proof of Lemma 1

Recall that $\theta=\left[\theta^{k}\right]_{k=1}^{p}$. Let $\mathcal{G}^{k}g=V_{\tau_{0}}[D_{\theta^{k}}\eta(\tau_{0}),D_{h}\eta(\tau_{0})g]$, for any $g\in{\cal H}$. Since

[TABLE]

for some positive constant $C$, $\mathcal{G}^{k}$ is a bounded linear functional on ${\cal H}$. By the Riesz representation theorem, there exists a $G^{k}\in{\cal H}$ such that for any $g\in{\cal H}$,

[TABLE]

Let $G=\left[G^{k}\right]_{k=1}^{p}$. We define $V_{\tau_{0}}[D_{\theta}\eta(\tau_{0}),D_{h}\eta(\tau_{0})g]$, $V_{\tau_{0}}[D_{h}\eta(\tau_{0})G,D_{h}\eta(\tau_{0})g]$, and $J(G,g)$ to be $p\times 1$ vectors whose $k$th entries are $V_{\tau_{0}}[D_{\theta^{k}}\eta(\tau_{0}),D_{h}\eta(\tau_{0})g]$, $V_{\tau_{0}}[D_{h}\eta(\tau_{0})G^{k},D_{h}\eta(\tau_{0})g]$, and $J(G^{k},g)$, respectively. Therefore,

[TABLE]

The Fourier expansion of $G^{k}$ given the eigensystem discussed in Section 4.2.3 is

[TABLE]

Simple calculation shows that

[TABLE]

and hence,

[TABLE]

Recall from Section 4.3 that

[TABLE]

where

[TABLE]

The $(i,j)$th entry of this $p\times p$ matrix $\Omega_{\lambda}$ can be written as

[TABLE]

Let $\Omega=V_{\tau_{0}}[D_{\theta}\eta(\tau_{0})-D_{h}\eta(\tau_{0})G]$ and $\Sigma_{\lambda}=\lambda J(G)$ be the matrices such that

[TABLE]

and $\Sigma_{\lambda}^{i,j}=\lambda J(G^{i},G^{j})$. Thus, $\Omega_{\lambda}=\Omega+\Sigma_{\lambda}$.

We now prove some properties of $\Omega$ and $\Sigma_{\lambda}$, which will be used to establish the bound for $\operatorname{E}[(\tilde{\theta}-\theta_{0})^{T}(\tilde{\theta}-\theta_{0})]$.

Lemma A.2.

$\Sigma_{\lambda}\rightarrow 0$* as $\lambda\rightarrow 0$.*

Proof.

By (40), we get that the $(i,j)$th entry of $\Sigma_{\lambda}$ is

[TABLE]

By square summability of $\{V_{\tau_{0}}[D_{\theta^{i}}\eta(\tau_{0}),D_{h}\eta(\tau_{0})\phi_{0,\nu}]\}_{\nu\in\mathbb{N}}$ and the dominated convergence theorem, the above sum converges to 0 as $\lambda\rightarrow 0$. ∎

Lemma A.3.

Under Assumption 2, $\Omega$ is positive definite. Let $c_{\delta}$ and $C_{1}$ be constants defined in Assumptions 2 and 4 respectively, and let $\delta$ be any eigenvalue of $\Omega$. Then $\delta>C_{1}c_{\delta}=\tilde{c}_{\delta}$.

Proof.

Using Assumption 4, the proof is similar to the one for Lemma 1. ∎

Note that by Lemma 3, the eigenvalues of $\Omega$ have a uniform lower bound independent of $\lambda$. Then as $\lambda\rightarrow 0$, we have

[TABLE]

where

[TABLE]

and $a_{n}^{i}$ is the $i$th entry of $a_{n}$.

Before we proceed to derive the bound for $\operatorname{E}\left[(a_{n}^{i})^{2}\right]$, we also need the following lemma, which is the same as Lemma 5.2 in \citeasnoungu_qiu_1993.

Lemma A.4.

Under Assumption 8, as $\lambda\rightarrow 0$,

[TABLE]

We now return to our analysis for $\operatorname{E}\left[(a_{n}^{i})^{2}\right]$. We have

[TABLE]

Note that $\operatorname{E}\left[(\bar{\alpha}_{n}^{i})^{2}\right]=O(n^{-1})$, and by square summability of $\{V_{\tau_{0}}[D_{\theta^{i}}\eta(\tau_{0}),D_{h}\eta(\tau_{0})\phi_{0,\nu}]\}_{\nu\in\mathbb{N}}$ and $\{h_{0,\nu}\}_{\nu\in\mathbb{N}}$, the dominated convergence theorem, the Cauchy-Schwarz inequality, $\operatorname{E}(\bar{\beta}^{2}_{\nu,n})=n^{-1}$, and Lemma 4, we have

[TABLE]

Therefore, we conclude that as $\lambda\rightarrow 0$, $\operatorname{E}\left[(a_{n}^{i})^{2}\right]=O(n^{-1}\lambda^{-\frac{1}{r}})$, which implies that

[TABLE]

This concludes the proof for the first bound in Lemma 1.

Now for the second bound in Lemma 1, we see that

[TABLE]

Plugging in the formula of $\tilde{h}_{\nu}$ given in Section 4.3, we get

[TABLE]

Since $\operatorname{E}(\bar{\beta}_{\nu,n})=0$, $\operatorname{E}(\bar{\beta}_{\nu,n}^{2})=\frac{1}{n}$, and $\sum_{\nu}\rho_{0,\nu}h^{2}_{0,\nu}=J(h_{0})<\infty$, by Lemma 4, we have

[TABLE]

If $a,b$ are $p\times 1$ vectors, then $a^{T}b$ is $1\times 1$, which implies that $a^{T}b=b^{T}a$ and $(a^{T}b)^{2}=b^{T}aa^{T}b$. Using this fact and the bound for $\operatorname{E}[(\tilde{\theta}-\theta_{0})^{T}(\tilde{\theta}-\theta_{0})]$, we have

[TABLE]

Therefore,

[TABLE]

Similar analysis shows that

[TABLE]

Hence, the second bound in Lemma 1 is established.

Proof of Theorem 1

By Assumption 4, Lemma 1 implies that

[TABLE]

By the Cauchy-Schwarz inequality and completing the square, we have

[TABLE]

Together with $\lambda J(\tilde{h}-h_{0})=O_{p}(n^{-1}\lambda^{-\frac{1}{r}}+\lambda)$, we have the desired result.

Proof of existence and uniqueness of $(\hat{\theta},\hat{h})$

In this section, adapting the frameworks in \citeasnouncox_1988, \citeasnounCoxOsullivan, \citeasnouno’sullivan_1990, and \citeasnounke2004existence, we establish the existence and uniqueness of the semiparametric estimator of the penalized likelihood given by (3) for the case when $m=1$ in the main text of the paper. The proof for the case when $m>1$ can be carried out in a similar manner. Note that we have the same assumptions as given in Section 4.2 for the proof of consistency in the main text, with the exception of assuming the existence of $(\hat{\theta},\hat{h})$ in $\mathcal{N}_{\theta_{0}}\times\mathcal{N}_{h_{0}}$. For convenience, we drop the subscript $l=1$. Denote the $k$th order partial Fréchet derivative operators by $D^{k}_{a_{1}\ldots a_{k}}=D_{a_{1}}\ldots D_{a_{k}}$, where $a_{i}\in\{{\theta},h\}$ for $i=1,\ldots,k$.

Linearization

We extend the linearization technique used to approximate the systematic and stochastic components of the estimation error in \citeasnounCoxOsullivan and \citeasnouno’sullivan_1990 to our semiparametric setting by using a bivariate Taylor series expansion for a nonlinear operator. Note that \citeasnounke2004existence used a similar analysis to show existence and uniqueness of a smoothing spline nonlinear nonparametric regression model. We first state the following proposition, whose proof is provided in \citeasnounke2004existence.

Proposition A.1.

Let $f:D(f)\subset X\times Y\to Z$, where $X$, $Y$ and $Z$ are Banach spaces. If $f^{\prime\prime}$ exists at $(x,y)$, then the partial Fréchet derivatives $f_{xx}$, $f_{xy}$, $f_{yx}$ and $f_{yy}$ exist at $(x,y)$. For any $h,a\in X$, $k,b\in Y$,

[TABLE]

By above theorem and the Taylor formula given in Chapter 1 Section 4 in \citeasnounambrosetti1995primer, we can write the first order Taylor series expansion of $f(x,y)$ as

[TABLE]

where $R$ is the remainder, given by

[TABLE]

Linear expansions

Since $D_{\theta}\mathfrak{L}_{\lambda}(\theta,h)$ is a bounded linear functional on $\mathbb{R}^{p}$, by the Riesz representation theorem, there exists $Z_{\theta}(\theta,h)\in\mathbb{R}^{p}$ such that for any $a\in\mathbb{R}^{p}$,

[TABLE]

Similarly, we can denote the Riesz representer of $D_{h}\mathfrak{L}_{\lambda}(\theta,h)$ in ${\cal H}$ by $Z_{h}(\theta,h)\in{\cal H}$. For convenience, we use

[TABLE]

to represent either the functionals or their Riesz representers in $\mathbb{R}^{p}$ and ${\cal H}$, respectively. For any $\theta_{0}+a\in\mathcal{N}_{\theta_{0}}$ and $h_{0}+g\in\mathcal{N}_{h_{0}}$, the first order Taylor series expansions of $Z_{h},Z_{\theta}$ at the true parameter $(\theta_{0},h_{0})$ are

[TABLE]

where

[TABLE]

For $Z_{h}(\theta,h)$, $u,v\in{\cal H}$, $a\in\mathbb{R}^{p}$, we have

[TABLE]

For $Z_{\theta}(\theta,h)$, $a,b\in\mathbb{R}^{p}$, $u\in{\cal H}$, we have

[TABLE]

For any $u,v\in{\cal H}$, ${a,b\in\mathbb{R}^{p}}$, and $\tau=(\theta,h)\in\mathcal{N}_{\theta_{0}}\times\mathcal{N}_{h_{0}}$, define the operators $U_{\theta}(\theta,h)$ and $U_{h}(\theta,h)$ on $\mathbb{R}^{p}$ and ${\cal H}$, respectively, such that

[TABLE]

Note that these operators are well-defined by the Riesz representation theorem applied to the linear functionals

[TABLE]

which are bounded in the corresponding norms on ${\cal H}$ and $\mathbb{R}^{p}$, respectively. Similarly, we also define $U_{h\theta}(\theta,h):{\cal H}\to\mathbb{R}^{p}$ and $U_{\theta h}(\theta,h):\mathbb{R}^{p}\to{\cal H}$ by

[TABLE]

By Lemma S.2 in the supplement of \citeasnouncheng2015joint, there exists a bounded linear operator $W_{\lambda}$ on ${\cal H}$ such that

[TABLE]

Therefore,

[TABLE]

where $G_{h}(\theta,h)=U_{h}(\theta,h)+W_{\lambda}$. We provide the presentations of the remainder terms $R_{h}(\theta,h)ag$ and $R_{\theta}(\theta,h)ag$ in Section Remainder terms.

Suppose $(\theta_{\lambda},h_{\lambda})$ is a solution for $Z_{\theta}(\theta,h)=Z_{h}(\theta,h)=0$. We define the systematic error as $(\theta_{\lambda}-\theta_{0},h_{\lambda}-h_{0})$. Ignoring the remainder terms, we can get an approximation to the systematic error by setting the system of equations (41) to 0 and solving for $\bar{\theta}_{\lambda}-\theta_{0}$, and $\bar{h}_{\lambda}-h_{0}$, i.e.,

[TABLE]

By the Lax-Milgram theorem (Section 3.6 of \citeasnounaubin_1979) and Assumptions 2 and 4 in the main text, for any $(\theta,h)\in\mathcal{N}_{\theta_{0}}\times\mathcal{N}_{h_{0}}$, the operators $G_{h}(\theta,h)$, $U_{\theta}(\theta,h)$ have bounded inverses on ${\cal H}$ and $\mathbb{R}^{p}$, respectively. Let

[TABLE]

Assuming both operators above have bounded inverses for any ${(\theta,h)\in\mathcal{N}_{\theta_{0}}\times\mathcal{N}_{h_{0}}}$, we get

[TABLE]

Next, we define the stochastic error as $(\hat{\theta}-\theta_{\lambda},\hat{h}-h_{\lambda})$. Similar to the definition of $Z_{\theta}$ and $Z_{h}$, we let

[TABLE]

The approximation of the stochastic errors can be obtained by the linearizations of $Z_{n\theta}$ and $Z_{nh}$. Since for $u,v\in{\cal H}$ and $a,b\in\mathbb{R}^{p}$, we have

[TABLE]

for any $\theta_{\lambda}+a\in\mathcal{N}_{\theta_{0}},h_{\lambda}+g\in\mathcal{N}_{h_{0}}$, the first order Taylor series expansions of $Z_{n\theta}$ and $Z_{nh}$ at $(\theta_{\lambda},h_{\lambda})\in\mathcal{N}_{\theta_{0}}\times\mathcal{N}_{h_{0}}$ can be written as

[TABLE]

where the error terms are given by

[TABLE]

and $R_{nh},R_{n\theta}$ are defined similarly to $R_{h},R_{\theta}$ by replacing $Z_{h},Z_{\theta}$ with $Z_{nh},Z_{n\theta}$, respectively. Recall that $(\hat{\theta},\hat{h})$ is the solution for $Z_{n\theta}(\theta,h)=Z_{nh}(\theta,h)=0$. Dropping the error terms and the remainder terms, we can get an approximation to the stochastic error ${(\hat{\theta}-\theta_{\lambda},\hat{h}-h_{\lambda})}$ by setting the linearizations (42) to 0 and solving for $\bar{\theta}_{n\lambda}-\theta_{\lambda}$, and $\bar{h}_{n\lambda}-h_{\lambda}$. We get

[TABLE]

Remainder terms

By definition of $R_{h}(\theta,h)$ and $R_{\theta}(\theta,h)$ given in Section Linear expansions, we need to find the second partial Fréchet derivatives of $Z_{\theta}(\theta,h)$ and $Z_{h}(\theta,h)$.

For $u,v,w\in{\cal H}$, $a,b,c\in\mathbb{R}^{p}$, we have

[TABLE]

[TABLE]

[TABLE]

[TABLE]

By replacing the terms $\mu_{\tau_{0}}[\cdot(x)]$ with $\frac{1}{n}\sum_{i=1}^{n}\cdot(x_{i})$ in each term above, we can have the second partial Fréchet derivatives of $Z_{n\theta}(\theta,h)$ and $Z_{nh}(\theta,h)$ for the remainder terms $R_{nh}(\theta,h)$ and $R_{n\theta}(\theta,h)$.

Bounds for the remainders

We see that the magnitude of the remainder terms $R_{\theta}$, $R_{h}$, $R_{n\theta}$, $R_{nh}$, $e_{\theta}$, and $e_{h}$ determine how accurate $(\bar{\theta}_{\lambda}-\theta_{0},\bar{h}_{\lambda}-h_{0})$ and $(\bar{\theta}_{n\lambda}-\theta_{\lambda},\bar{h}_{n\lambda}-h_{\lambda})$ are as approximations of the systematic error and the stochastic error, respectively. To obtain measures of these terms, we first define that for $\lambda>0$, $\theta_{1},\theta_{2}\in\mathcal{N}_{\theta_{0}}$, $h_{1},h_{2}\in\mathcal{N}_{h_{0}}$, and unit elements $u_{1},u_{2}\in\mathbb{R}^{p}$ and $v_{1},v_{2}\in{\cal H}$,

[TABLE]

For ${i=1,2,3,4}$, we also define $K^{i}_{nh},K^{i}_{n\theta}$ by replacing $Z_{\theta},Z_{h}$ with $Z_{n\theta},Z_{nh}$ in $K^{i}_{h}$ and $K^{i}_{\theta}$, respectively. In addition, for $\lambda>0$, $\theta_{1},\theta_{2}\in\mathcal{N}_{\theta_{0}}$, $h_{1},h_{2}\in\mathcal{N}_{h_{0}}$, and unit elements $u_{1},u_{2}\in\mathbb{R}^{p}$ and $v_{1},v_{2}\in{\cal H}$, we define

[TABLE]

Therefore, for any $a\in\mathbb{R}^{p}$ and $g\in{\cal H}$, standard analysis allows us to have the following bounds for the remainder terms for the systematic error and the stochastic error,

[TABLE]

[TABLE]

[TABLE]

where $E_{nh}^{1}=E_{nh}^{11}+E_{nh}^{21}$ and $E_{nh}^{2}=E_{nh}^{12}+E_{nh}^{22}$, and

[TABLE]

where $E_{n\theta}^{1}=E_{n\theta}^{21}+E_{n\theta}^{11}$ and $E_{n\theta}^{2}=E_{n\theta}^{22}+E_{n\theta}^{12}$.

Proof of Existence and Uniqueness

We are now ready to show the existence and uniqueness of $(\theta_{\lambda},h_{\lambda})$ and $(\hat{\theta},\hat{h})$ in the neighborhood $\mathcal{N}_{\theta_{0}}\times\mathcal{N}_{h_{0}}$. Let

[TABLE]

One can get the following theorem for the existence and uniqueness of $(\theta_{\lambda},h_{\lambda})$ via a contraction mapping argument.

Theorem A.2.

(Existence of $(\theta_{\lambda},h_{\lambda})$ and the bias approximation) If $d_{\theta}(\lambda)\rightarrow 0,d_{h}(\lambda)\rightarrow 0,r_{\theta}(\lambda)\rightarrow 0,r_{h}(\lambda)\rightarrow 0$ as $\lambda\rightarrow 0$, there exist $\lambda_{0}>0$ such that for $\lambda\in[0,\lambda_{0}]$, there are unique $\theta_{\lambda}\in S_{\theta,\theta_{0}}(2d_{\theta}(\lambda))$ and $h_{\lambda}\in S_{h,h_{0}}(2d_{h}(\lambda))$ satisfying $Z_{\theta}(\theta_{\lambda},h_{\lambda})=0,Z_{h}(\theta_{\lambda},h_{\lambda})=0$, and $(\theta_{\lambda},h_{\lambda})\in\mathcal{N}_{\theta_{0}}\times\mathcal{N}_{h_{0}}$. In addition, as $\lambda\to 0$,

[TABLE]

Proof.

Let $t_{\theta\lambda}=2d_{\theta}(\lambda)$, $t_{h\lambda}=2d_{h}(\lambda)$. Define

[TABLE]

Let $\vec{F}(\zeta,g)=(F_{\theta}(\zeta,g),F_{h}(\zeta,g))$ be a function on $\mathcal{Q}=\mathbb{R}^{p}\times{\cal H}$, and for any subset $\mathcal{Q}_{1}\subset\mathcal{Q}$, denote by $\vec{F}(\mathcal{Q}_{1})$ the image of $\mathcal{Q}_{1}$ under $\vec{F}$. The proof has three steps:

$\vec{F}(S_{\theta}(t_{\theta\lambda})\times S_{h}(t_{h\lambda}))\subset S_{\theta}(t_{\theta\lambda})\times S_{h}(t_{h\lambda})$. 2. 2.

$\vec{F}$ is a contraction map on $S_{\theta}(t_{\theta\lambda})\times S_{h}(t_{h\lambda})$. 3. 3.

Obtaining the estimate for the bias approximation, $\left\|\bar{\theta}_{\lambda}-\theta_{\lambda}\right\|_{\mathbb{R}^{p}}+\left\|\bar{h}_{\lambda}-h_{\lambda}\right\|_{{\cal H}}$.

For (a), by our assumption, we can choose $\lambda_{0}$ small enough that $S_{\theta,\theta_{0}}(t_{\theta\lambda})\subset\mathcal{N}_{\theta_{0}}$, $S_{h,h_{0}}(t_{h\lambda})\subset\mathcal{N}_{h_{0}}$, and $r_{\theta}(\lambda)<1/2$ for all $\lambda\in(0,\lambda_{0}]$. For every $(\theta,h)\in\mathcal{Q}$, we denote

[TABLE]

For $(\zeta,g)\in S_{\theta}(t_{\theta\lambda})\times S_{h}(t_{h\lambda})$, we have

[TABLE]

For $\left\lVert F_{\theta}(\zeta,g)\right\rVert_{\mathbb{R}^{p}}$, by triangle inequality, we have

[TABLE]

By the definition of $\bar{\theta}_{\lambda}-\theta_{0}$ and $G_{\theta\theta}(\theta,h)$, Taylor series expansion of $Z_{\theta}(\theta_{0}+\zeta,h_{0}+g)$ and $Z_{h}(\theta_{0}+\zeta,h_{0}+g)$, and the remainder bound (44), we can get

[TABLE]

Similarly, by the definition of $\bar{h}_{\lambda}-h_{0}$ and $G_{hh}(\theta,h)$, Taylor series expansion of $Z_{\theta}(\theta_{0}+\zeta,h_{0}+g)$ and $Z_{h}(\theta_{0}+\zeta,h_{0}+g)$, and the remainder bound (43), we also have

[TABLE]

and

[TABLE]

Since $t_{\theta\lambda}=2\left\lVert\bar{\theta}_{\lambda}-\theta_{0}\right\rVert_{\mathbb{R}^{p}}$, $t_{h\lambda}=2\left\lVert\bar{h}_{\lambda}-h_{0}\right\rVert_{{\cal H}}$, $r_{h}(\lambda)<1/2$, and $r_{\theta}(\lambda)<1/2$, for $(\zeta,g)\in{S_{\theta}(t_{\theta\lambda})\times S_{h}(t_{h\lambda})}$, we have

[TABLE]

Now for step (b): For $\zeta_{1},\zeta_{2}\in S_{\theta}(t_{\theta\lambda})$, $g_{1},g_{2}\in S_{h}(t_{h\lambda})$, by Taylor expansion, we get

[TABLE]

Applying Taylor expansion again to the terms inside the integral and letting $\zeta^{*}=\zeta_{1}+t(\zeta_{2}-\zeta_{1})$, ${g^{*}=g_{1}+t(g_{2}-g_{1})}$, we have

[TABLE]

Since

[TABLE]

and for $0\leq t\leq 1$, $\zeta^{*}=\zeta_{1}+t(\zeta_{2}-\zeta_{1})\in S_{\theta}(t_{\theta\lambda})$, $g^{*}=g_{1}+t(g_{2}-g_{1})\in S_{h}(t_{h\lambda})$ by convexity of $S_{\theta}(t_{\theta\lambda})$ and $S_{h}(t_{h\lambda})$, similar algebraic manipulations as in the proof of (a) show that

[TABLE]

Similarly for $F_{h}$, we get

[TABLE]

Therefore,

[TABLE]

where $0<C_{1}<1$, $0<C_{2}<1$, so $\vec{F}=(F_{\theta},F_{h})$ is a contraction map on $S_{\theta}(t_{\theta\lambda})\times S_{h}(t_{h\lambda})$. By the contraction mapping theorem (see Theorem 9.23 in \citeasnounrudin1976principles), there exists a unique ${(\zeta_{\lambda},g_{\lambda})\in S_{\theta}(t_{\theta\lambda})\times S_{h}(t_{h\lambda})}$ such that

[TABLE]

Let $\theta_{\lambda}=\theta_{0}+\zeta_{\lambda}$, $h_{\lambda}=h_{0}+g_{\lambda}$. Then $\theta_{\lambda}\in S_{\theta,\theta_{0}}(t_{\theta\lambda})$, $h_{\lambda}\in S_{h,h_{0}}(t_{\theta\lambda})$, and $(\theta_{\lambda},h_{\lambda})$ are the unique solutions to $Z_{\theta}(\theta_{\lambda},h_{\lambda})=0$, $Z_{h}(\theta_{\lambda},h_{\lambda})=0$.

For part (c), note that

[TABLE]

Thus,

[TABLE]

This completes the proof of Theorem A.2. ∎

Next, we consider the existence of $(\hat{\theta},\hat{h})\in\mathcal{N}_{\theta_{0}}\times\mathcal{N}_{h_{0}}$. Define

[TABLE]

We get the following existence theorem for $(\hat{\theta},\hat{h})\in\mathcal{N}_{\theta_{0}}\times\mathcal{N}_{h_{0}}$.

Theorem A.3.

(Existence of $(\hat{\theta},\hat{h})$ and the variability approximation) Suppose $\lambda_{n}$ is a sequence such that for all $n$ sufficiently large, $\theta_{\lambda_{n}}\in\mathcal{N}_{\theta_{0}}$, $h_{\lambda_{n}}\in\mathcal{N}_{h_{0}}$, and

[TABLE]

Then, with probability tending to unity as $n\to\infty$, there is a unique root $(\hat{\theta},\hat{h})$ of $Z_{n\theta}(\hat{\theta},\hat{h})=0$, $Z_{nh}(\hat{\theta},\hat{h})=0$ in ${S_{\theta,\theta_{\lambda_{n}}}(2d_{n\theta}(\lambda_{n}))\times S_{h,h_{\lambda_{n}}}(2d_{nh}(\lambda_{n}))}\subset\mathcal{N}_{\theta_{0}}\times\mathcal{N}_{h_{0}}$. In addition, as $n\to\infty$, $\lambda_{n}\to 0$,

[TABLE]

Proof.

For convenience, we drop the subscript on $\lambda_{n}$ and let $t_{n\theta\lambda}=2d_{n\theta}(\lambda)$, $t_{nh\lambda}=2d_{nh}(\lambda)$. Let

[TABLE]

The proof proceeds in three steps, similar to the proof of Theorem 2, with additional terms introduced in approximating $D_{\theta}Z_{n\theta}$ and $D_{h}Z_{nh}$ by $D_{\theta}Z_{\theta}$ and $D_{h}Z_{h}$, respectively. Take $n$ large enough so that $S_{\theta,\theta_{\lambda}}(t_{n\theta\lambda})\subset\mathcal{N}_{\theta_{0}}$, $S_{h,h_{\lambda}}(t_{nh\lambda})\subset\mathcal{N}_{h_{0}}$ and $r_{n\theta}(\lambda)<\frac{1}{2}$, $r_{nh}(\lambda)<\frac{1}{2}$.

First, we show that $\vec{F}_{n}(\zeta,g)=(F_{n\theta}(\zeta,g),F_{nh}(\zeta,g))$ maps $S_{\theta}(t_{n\theta\lambda})\times S_{h}(t_{nh\lambda})$ to itself, i.e.,

[TABLE]

By definition, for $(\zeta,g)\in S_{\theta}(t_{n\theta\lambda})\times S_{h}(t_{nh\lambda})$, we have

[TABLE]

For $F_{n\theta}$, by the triangle inequality, we get

[TABLE]

Using the definition of $\bar{\theta}_{n\lambda}-\theta_{\lambda}$, $G_{\theta\theta}(\theta,h)$, Taylor expansion of $Z_{n\theta}(\theta_{\lambda}+\zeta,h_{\lambda}+g)$ and $Z_{nh}(\theta_{\lambda}+\zeta,h_{\lambda}+g)$, and the remainder bound (46), we get that

[TABLE]

Similarly, we have

[TABLE]

and

[TABLE]

Thus, for $(\zeta,g)\in S_{\theta}(t_{n\theta\lambda})\times S_{h}(t_{nh\lambda})$,

[TABLE]

Therefore, we have shown that $\vec{F}_{n}(S_{\theta}(t_{n\theta\lambda})\times S_{h}(t_{nh\lambda}))\subset S_{\theta}(t_{n\theta\lambda})\times S_{h}(t_{nh\lambda})$.

Next, we show that $\vec{F}_{n}$ is a contraction map. By similar calculations as in the proof for Theorem 1, after applying Taylor expansion twice, for $\zeta_{1},\zeta_{2}\in S_{\theta}(t_{n\theta\lambda})$, $g_{1},g_{2}\in S_{h}(t_{nh\lambda})$, we get

[TABLE]

[TABLE]

Thus,

[TABLE]

Since $r_{n\lambda}(\lambda)<\frac{1}{2}$, $r_{nh}(\lambda)<\frac{1}{2}$, we have shown that $\vec{F}(\zeta,g)$ is a contraction on $S_{\theta}(t_{n\theta\lambda})\times S_{h}(t_{nh\lambda})$. By the contraction mapping theorem, there exists a unique $(\zeta_{n\lambda},g_{n\lambda})\in S_{\theta}(t_{n\theta\lambda},\alpha)\times S_{h}(t_{nh\lambda},\alpha)$ such that

[TABLE]

Let $\hat{\theta}=\theta_{\lambda}+\zeta_{n\lambda}\in S_{\theta,\theta_{\lambda}}(2d_{n\theta}(\lambda))$ and $\hat{h}=h_{\lambda}+g_{n\lambda}\in S_{h,h_{\lambda}}(2d_{nh}(\lambda))$. Then $(\hat{\theta},\hat{h})$ is the unique root of $Z_{n\theta}(\hat{\theta},\hat{h})=0$ and $Z_{nh}(\hat{\theta},\hat{h})=0$.

To get the upper bound, we observe that

[TABLE]

Therefore,

[TABLE]

This completes the proof. ∎