Maximum Likelihood Estimation of a Semiparametric Two-component Mixture Model using Log-concave Approximation

TL;DR

This paper introduces a semiparametric maximum likelihood estimation method for a two-component mixture model with a known component, assuming the unknown component's density is log-concave, with applications in biology and astronomy.

Contribution

It develops a novel EM algorithm for joint estimation of mixing proportions and unknown distribution under log-concavity assumptions, with proven identifiability and consistency.

Findings

Estimator performs well in simulations

Method successfully applied to biological data

Approach outperforms existing estimators in certain scenarios

Abstract

Motivated by studies in biological sciences to detect differentially expressed genes, a semiparametric two-component mixture model with one known component is being studied in this paper. Assuming the density of the unknown component to be log-concave, which contains a very broad family of densities, we develop a semiparametric maximum likelihood estimator and propose an EM algorithm to compute it. Our new estimation method finds the mixing proportions and the distribution of the unknown component simultaneously. We establish the identifiability of the proposed semiparametric mixture model and prove the existence and consistency of the proposed estimators. We further compare our estimator with several existing estimators through simulation studies and apply our method to two real data sets from biological sciences and astronomy.

Figures

Tables

Table 1. Table 1 : Bias(MSE) of estimates of p / μ 𝑝 𝜇 p/\mu and mean of the classification error for model 1 when n = 1000 𝑛 1000 n=1000 .

	EM_logconcave	Patra	Bordes	Song EM	Song max $\pi$	Xiang
$p=0.2$
$p$	0.002(0.0004)	0.009(0.0007)	0.001(0.0009)	0.08(0.0066)	0.087(0.0122)	0.006(0.0006)
$\mu$	0.063(0.0180)	-0.152(0.0650)	-0.021(0.0426)	0.116(0.0446)	-0.675(0.5680)	0.116(0.0396)
Cla_error	0.0960	0.1056	0.1052	0.1102	0.1052	0.0973
$p=0.5$
$p$	-0.002(0.0004)	-0.025(0.0011)	0.001(0.0006)	-0.132(0.0177)	0.106(0.0149)	0.007(0.0006)
$\mu$	0.018(0.0042)	0.051(0.0073)	0.000(0.0046)	0.185(0.0375)	-0.322(0.1392)	0.013(0.0056)
Cla_error	0.1094	0.1219	0.1198	0.1352	0.1206	0.1104
$p=0.8$
$p$	0.001(0.0002)	-0.252(0.0020)	0.001(0.0003)	-0.107(0.0118)	0.063(0.0047)	0.009(0.0003)
$\mu$	0.005(0.0013)	0.066(0.0057)	0.000(0.0016)	0.118(0.0153)	-0.128(0.0220)	-0.002(0.0021)
Cla_error	0.0645	0.0739	0.0694	0.0834	0.0721	0.0664

Table 2. Table 2 : Bias(MSE) of estimates of p / μ 𝑝 𝜇 p/\mu and mean of the classification error for model 2 when n = 1000 𝑛 1000 n=1000 .

	EM_logconcave	Patra	Bordes	Song EM	Song max $\pi$	Xiang
$p=0.2$
$p$	-0.008(0.0014)	-0.023(0.0015)	-0.015(0.0012)	-0.15(0.0228)	0.382(0.1496)	0.017(0.0019)
$\mu$	-0.018(0.0015)	0.027(0.0016)	-0.029(0.0017)	-0.014(0.0007)	0.199(0.0401)	-0.007(0.0010)
Cla_error	0.1270	0.1520	0.1511	0.1676	0.1847	0.1339
$p=0.5$
$p$	0.001(0.0007)	-0.046(0.0030)	-0.040(0.0024)	-0.248(0.0811)	0.228(0.0548)	-0.047(0.0035)
$\mu$	-0.003(0.0001)	-0.011(0.0002)	-0.032(0.0011)	-0.038(0.0015)	0.077(0.0064)	-0.022(0.0012)
Cla_error	0.1609	0.1990	0.1974	0.2638	0.1887	0.1753
$p=0.8$
$p$	-0.001(0.0004)	-0.074(0.0059)	-0.070(0.0055)	-0.311(0.0974)	0.099(0.0105)	-0.060(0.0043)
$\mu$	-0.001(0.00004)	-0.019(0.0004)	-0.033(0.0011)	-0.040(0.0016)	0.025(0.0007)	-0.030(0.0014)
Cla_error	0.1000	0.1264	0.1261	0.2103	0.1129	0.1142

Table 3. Table 3 : Bias(MSE) of estimates of p / μ 𝑝 𝜇 p/\mu and mean of the classification error for model 3 when n = 1000 𝑛 1000 n=1000 .

	EM_logconcave	Patra	Bordes	Song EM	Song max $\pi$	Xiang
$p=0.2$
$p$	0.001(0.0002)	-0.001(0.0006)	NA	-0.060(0.0038)	0.410(0.1698)	0.024(0.0011)
$\mu$	0.006(0.0082)	-0.039(0.0152)	NA	0.048(0.0149)	-1.140(1.3094)	-0.105(0.0184)
Cla_error	0.0709	0.0851	NA	0.0879	0.1568	0.0790
$p=0.5$
$p$	0.000(0.0003)	-0.013(0.0006)	NA	-0.073(0.0057)	0.259(0.0681)	0.042(0.0028)
$\mu$	0.003(0.0021)	-0.011(0.0030)	NA	0.018(0.0030)	-0.502(0.2578)	-0.091(0.0157)
Cla_error	0.0595	0.0767	NA	0.0790	0.1166	0.0732
$p=0.8$
$p$	0.001(0.0002)	-0.228(0.0010)	NA	-0.231(0.0012)	0.104(0.0112)	0.071(0.0060)
$\mu$	-0.001(0.0013)	-0.002(0.0014)	NA	-0.002(0.0014)	-0.159(0.0283)	-0.104(0.0224)
Cla_error	0.0260	0.0325	NA	0.0322	0.0526	0.0617

Table 4. Table 4 : Bias(MSE) of estimates of p / μ 𝑝 𝜇 p/\mu and mean of the classification error for model 4 when n = 1000 𝑛 1000 n=1000 .

	EM_logconcave	Patra	Bordes	Song EM	Song max $\pi$	Xiang
$p=0.2$
$p$	-0.000(0.0002)	0.005(0.0005)	NA	0.006(0.0003)	0.106(0.0160)	0.056(0.0041)
$\mu$	-0.023(0.0321)	-0.286(0.1299)	NA	-0.304(0.1353)	-1.066(1.4174)	-0.738(1.0279)
Cla_error	0.0112	0.0139	NA	0.0137	0.0205	0.0215
$p=0.5$
$p$	0.000(0.0002)	-0.009(0.0004)	NA	0.014(0.0005)	0.067(0.0057)	0.049(0.0030)
$\mu$	-0.005(0.0074)	-0.148(0.0332)	NA	-0.185(0.0459)	-0.333(0.1439)	-0.676(0.6616)
Cla_error	0.0110	0.0157	NA	0.0163	0.0160	0.0207
$p=0.8$
$p$	-0.001(0.0001)	-0.023(0.0007)	NA	0.006(0.0002)	0.038(0.0019)	0.066(0.0048)
$\mu$	-0.002(0.0052)	-0.025(0.0077)	NA	-0.054(0.0098)	-0.147(0.0352)	-0.718(0.5396)
Cla_error	0.0045	0.0078	NA	0.0089	0.0108	0.0319

Table 5. Table 5 : Bias(MSE) of estimates of p / μ 𝑝 𝜇 p/\mu and mean of the classification error for model 5 when n = 1000 𝑛 1000 n=1000 .

	EM_logconcave	Patra	Bordes	Song EM	Song max $\pi$	Xiang
$p=0.2$
$p$	0.001(0.0002)	0.006(0.0006)	NA	0.018(0.0005)	0.110(0.0154)	0.037(0.0018)
$\mu$	0.004(0.0193)	-0.399(0.2021)	NA	-0.427(0.2126)	-1.129(1.4965)	-0.923(0.9073)
Cla_error	0.0012	0.0044	NA	0.0045	0.0105	0.0092
$p=0.5$
$p$	0.000(0.0003)	-0.009(0.0005)	NA	0.023(0.0008)	0.131(0.0208)	0.061(0.0044)
$\mu$	-0.001(0.0079)	-0.173(0.0384)	NA	-0.213(0.0538)	-0.681(0.5620)	-0.650(0.4645)
Cla_error	0.0007	0.0079	NA	0.0090	0.0202	0.0189
$p=0.8$
$p$	0.001(0.0001)	-0.223(0.0007)	NA	0.009(0.0002)	0.079(0.0068)	0.081(0.0071)
$\mu$	0.000(0.0047)	-0.027(0.0061)	NA	-0.051(0.0077)	-0.313(0.1123)	-0.473(0.2482)
Cla_error	0.0003	0.0022	NA	0.0030	0.0190	0.0426

Table 6. Table 6 : Bias(MSE) of estimates of p / μ 𝑝 𝜇 p/\mu and mean of the classification error for model 6 when n = 1000 𝑛 1000 n=1000 .

	EM_logconcave	Patra	Bordes	Song EM	Song max $\pi$	Xiang
$p=0.2$
$p$	-0.010(0.0003)	-0.007(0.0006)	0.083(0.0012)	-0.022(0.0006)	0.077(0.0080)	0.001(0.0003)
$\mu$	0.171(0.0455)	0.029(0.0262)	-0.075(0.0843)	0.083(0.0257)	-0.414(0.2369)	0.001(0.0177)
Cla_error	0.0440	0.0450	0.0455	0.0457	0.0468	0.0435
$p=0.5$
$p$	0.001(0.0003)	-0.031(0.0015)	-0.002(0.0006)	-0.053(0.0031)	0.038(0.0028)	-0.002(0.0631)
$\mu$	-0.010(0.0115)	0.148(0.0264)	-0.003(0.0066)	0.182(0.0375)	-0.020(0.0185)	0.018(0.0066)
Cla_error	0.0094	0.0672	0.0658	0.0680	0.0656	0.0631
$p=0.8$
$p$	-0.001(0.0001)	-0.059(0.0037)	-0.002(0.0004)	-0.063(0.0043)	0.008(0.0006)	0.001(0.0034)
$\mu$	-0.004(0.0072)	0.169(0.0307)	-0.001(0.0024)	0.174(0.0321)	0.055(0.0073)	-0.003(0.0034)
Cla_error	0.0046	0.0637	0.0570	0.0643	0.0567	0.0545

Table 7. Table 7 : Estimates of p 𝑝 p for the prostate cancer data.

EM_logconcave	Patra	Bordes	Song EM	Song max  $\pi$	Xiang
0.0173	0.0817	0.1975	0.0076	0.6132	0.1915

Table 8. Table 8 : Estimates of p 𝑝 p for the Carina data.

EM_logconcave	Patra	Bordes	Song EM	Song max  $\pi$	Xiang
0.354	0.364	0.363	0.370	0.687	0.385

Table 9. Table 9 : Bias(MSE) of estimates of p / μ 𝑝 𝜇 p/\mu and mean of the classification error for model 1 when n = 250 𝑛 250 n=250 .

	EM_logconcave	Patra	Bordes	Song EM	Song max $\pi$	Xiang
$p=0.2$
$p$	0.008(0.0019)	-0.001(0.0021)	0.018(0.0046)	-0.071(0.0057)	0.106(0.0160)	0.021(0.0056)
$\mu$	0.057(0.0738)	-0.396(0.3286)	-0.166(0.2109)	-0.108(0.2243)	-0.846(0.9437)	0.243(0.1864)
Cla_error	0.1029	0.1094	0.1097	0.1138	0.1104	0.1058
$p=0.5$
$p$	0.000(0.0023)	-0.041(0.0036)	0.005(0.0025)	-0.130(0.0185)	0.100(0.0153)	0.014(0.0026)
$\mu$	0.017(0.0225)	0.023(0.0167)	-0.021(0.0198)	0.143(0.0393)	-0.344(0.1635)	0.032(0.0208)
Cla_error	0.1151	0.1259	0.1232	0.1379	0.1248	0.1138
$p=0.8$
$p$	-0.001(0.0011)	-0.070(0.0057)	-0.001(0.0014)	-0.104(0.0123)	0.056(0.0040)	0.016(0.0014)
$\mu$	0.003(0.0072)	0.059(0.0102)	-0.001(0.0085)	0.097(0.0158)	-0.147(0.0323)	-0.008(0.0079)
Cla_error	0.0670	0.0781	0.0722	0.0835	0.0752	0.0703

Table 10. Table 10 : Bias(MSE) of estimates of p / μ 𝑝 𝜇 p/\mu and mean of the classification error for model 1 when n = 500 𝑛 500 n=500 .

	EM_logconcave	Patra	Bordes	Song EM	Song max $\pi$	Xiang
$p=0.2$
$p$	0.000(0.0008)	-0.008(0.0014)	0.003(0.0021)	-0.077(0.0063)	0.086(0.0102)	0.011(0.0013)
$\mu$	0.059(0.0348)	-0.258(0.1430)	-0.054(0.1013)	0.001(0.0734)	-0.738(0.6525)	0.175(0.0843)
Cla_error	0.0972	0.1070	0.1060	0.1106	0.1051	0.0990
$p=0.5$
$p$	-0.003(0.0009)	-0.031(0.0021)	0.000(0.0012)	-0.132(0.0181)	0.107(0.0158)	0.011(0.0011)
$\mu$	0.019(0.0097)	0.035(0.0109)	-0.003(0.0098)	0.169(0.0363)	-0.346(0.1605)	0.020(0.0097)
Cla_error	0.1111	0.1239	0.1209	0.1359	0.1226	0.1120
$p=0.8$
$p$	0.003(0.0006)	-0.053(0.0033)	0.001(0.0007)	-0.104(0.0117)	0.056(0.0040)	0.014(0.0006)
$\mu$	-0.001(0.0032)	0.065(0.0073)	0.000(0.0041)	0.110(0.0155)	-0.121(0.0220)	-0.007(0.0040)
Cla_error	0.0644	0.0758	0.0693	0.0822	0.0711	0.0685

Table 11. Table 11 : Bias(MSE) of estimates of p / μ 𝑝 𝜇 p/\mu and mean of the classification aerror for model 2 when n = 250 𝑛 250 n=250 .

	EM_logconcave	Patra	Bordes	Song EM	Song max $\pi$	Xiang
$p=0.2$
$p$	0.004(0.0051)	-0.021(0.0033)	-0.006(0.0036)	-0.156(0.0248)	0.371(0.1443)	0.056(0.0087)
$\mu$	-0.022(0.0038)	0.061(0.0073)	-0.019(0.0029)	0.013(0.0032)	0.197(0.0401)	-0.011(0.0026)
Cla_error	0.1368	0.1554	0.1568	0.1746	0.1858	0.1437
$p=0.5$
$p$	0.004(0.0043)	-0.064(0.0071)	-0.037(0.0041)	-0.300(0.0916)	0.230(0.0576)	-0.013(0.0041)
$\mu$	-0.005(0.0007)	-0.004(0.0004)	-0.032(0.0013)	-0.034(0.0013)	0.080(0.0070)	-0.014(0.0011)
Cla_error	0.1678	0.2110	0.2031	0.2778	0.1964	0.1764
$p=0.8$
$p$	0.009(0.0019)	-0.105(0.0124)	-0.065(0.0057)	-0.312(0.1001)	0.081(0.0078)	-0.049(0.0056)
$\mu$	0.000(0.0002)	-0.020(0.0005)	-0.033(0.0012)	-0.039(0.0016)	0.010(0.0006)	-0.018(0.0010)
Cla_error	0.1030	0.1384	0.1266	0.2137	0.1124	0.1140

Table 12. Table 12 : Bias(MSE) of estimates of p / μ 𝑝 𝜇 p/\mu and mean of the classification error for model 2 when n = 500 𝑛 500 n=500 .

	EM_logconcave	Patra	Bordes	Song EM	Song max $\pi$	Xiang
$p=0.2$
$p$	-0.001(0.0028)	-0.021(0.0020)	-0.007(0.0022)	-0.154(0.0238)	0.379(0.1496)	0.035(0.0047)
$\mu$	-0.019(0.0024)	0.041(0.0030)	-0.025(0.0022)	-0.006(0.0008)	0.197(0.0393)	-0.009(0.0018)
Cla_error	0.1294	0.1544	0.1520	0.1697	0.1862	0.1369
$p=0.5$
$p$	0.002(0.0022)	-0.053(0.0045)	-0.039(0.0032)	-0.292(0.0860)	0.234(0.0583)	-0.034(0.0035)
$\mu$	-0.004(0.0003)	-0.008(0.0002)	-0.033(0.0013)	-0.037(0.0014)	0.080(0.0069)	-0.014(0.0009)
Cla_error	0.1638	0.2031	0.2011	0.2723	0.1940	0.1753
$p=0.8$
$p$	0.003(0.0010)	-0.086(0.0081)	-0.070(0.0058)	-0.312(0.0990)	0.097(0.0102)	-0.048(0.0038)
$\mu$	-0.001(0.0001)	-0.020(0.0005)	-0.034(0.0012)	-0.040(0.0016)	0.024(0.0007)	-0.022(0.0010)
Cla_error	0.1001	0.1307	0.1263	0.2119	0.1129	0.1139

Table 13. Table 13 : Bias(MSE) of estimates of p / μ 𝑝 𝜇 p/\mu and mean of the classification error for model 3 when n = 250 𝑛 250 n=250 .

	EM_logconcave	Patra	Bordes	Song EM	Song max $\pi$	Xiang
$p=0.2$
$p$	0.005(0.0011)	0.009(0.0021)	NA	-0.050(0.0034)	0.431(0.1889)	0.048(0.0042)
$\mu$	0.026(0.0400)	-0.041(0.0581)	NA	0.048(0.0591)	-1.115(1.2677)	-0.139(0.0493)
Cla_error	0.0737	0.0869	NA	0.0895	0.1718	0.0842
$p=0.5$
$p$	0.002(0.0013)	-0.019(0.0019)	NA	-0.069(0.0065)	0.269(0.0742)	0.081(0.0103)
$\mu$	-0.001(0.0096)	-0.001(0.0126)	NA	0.013(0.0120)	-0.495(0.2618)	-0.174(0.0590)
Cla_error	0.0623	0.0806	NA	0.0839	0.1271	0.0860
$p=0.8$
$p$	0.003(0.0007)	-0.046(0.0029)	NA	-0.225(0.0017)	0.107(0.0122)	0.087(0.0096)
$\mu$	0.001(0.0052)	0.003(0.0061)	NA	-0.004(0.0057)	-0.157(0.0316)	-0.150(0.0446)
Cla_error	0.0274	0.0351	NA	0.0354	0.0589	0.0734

Table 14. Table 14 : Bias(MSE) of estimates of p / μ 𝑝 𝜇 p/\mu and mean of the classification error for model 3 when n = 500 𝑛 500 n=500 .

	EM_logconcave	Patra	Bordes	Song EM	Song max $\pi$	Xiang
$p=0.2$
$p$	0.004(0.0005)	0.004(0.0011)	NA	-0.055(0.0034)	0.415(0.1746)	0.038(0.0024)
$\mu$	0.014(0.0204)	-0.039(0.0316)	NA	0.047(0.0315)	-1.119(1.2657)	-0.131(0.0357)
Cla_error	0.0722	0.0862	NA	0.0855	0.1610	0.0819
$p=0.5$
$p$	0.001(0.0005)	-0.016(0.0010)	NA	-0.070(0.0057)	0.260(0.0692)	0.060(0.0059)
$\mu$	0.004(0.0047)	-0.007(0.0061)	NA	0.017(0.0064)	-0.489(0.2475)	-0.126(0.0338)
Cla_error	0.0604	0.0787	NA	0.0811	0.1189	0.0790
$p=0.8$
$p$	0.002(0.0003)	-0.036(0.0017)	NA	-0.029(0.0013)	0.106(0.0115)	0.080(0.0078)
$\mu$	0.001(0.0027)	0.001(0.0026)	NA	-0.002(0.0030)	-0.159(0.0294)	-0.117(0.0284)
Cla_error	0.0270	0.0334	NA	0.0341	0.0557	0.0677

Table 15. Table 15 : Bias(MSE) of estimates of p / μ 𝑝 𝜇 p/\mu and mean of the classification error for model 4 when n = 250 𝑛 250 n=250 .

	EM_logconcave	Patra	Bordes	Song EM	Song max $\pi$	Xiang
$p=0.2$
$p$	-0.001(0.0009)	0.019(0.0021)	NA	0.019(0.0013)	0.129(0.0198)	0.113(0.0179)
$\mu$	0.016(0.1454)	-0.574(0.5652)	NA	-0.569(0.5268)	-1.296(2.0126)	-0.728(1.0325)
Cla_error	0.0128	0.0168	NA	0.0168	0.0247	0.0438
$p=0.5$
$p$	0.003(0.0009)	-0.014(0.0014)	NA	0.022(0.0014)	0.085(0.0089)	0.086(0.0101)
$\mu$	0.025(0.0493)	-0.192(0.0898)	NA	-0.200(0.0878)	-0.379(0.2136)	-0.605(0.6314)
Cla_error	0.0096	0.0184	NA	0.0188	0.0204	0.0378
$p=0.8$
$p$	0.000(0.0006)	-0.039(0.0022)	NA	0.013(0.0007)	0.044(0.0026)	0.077(0.0076)
$\mu$	0.009(0.0279)	-0.024(0.0247)	NA	-0.073(0.0332)	-0.157(0.0572)	-0.621(0.4861)
Cla_error	0.0044	0.0093	NA	0.0115	0.0149	0.0468

Table 16. Table 16 : Bias(MSE) of estimates of p / μ 𝑝 𝜇 p/\mu and mean of the classification error for model 4 when n = 500 𝑛 500 n=500 .

	EM_logconcave	Patra	Bordes	Song EM	Song max $\pi$	Xiang
$p=0.2$
$p$	-0.002(0.0003)	0.009(0.0009)	NA	0.010(0.0005)	0.108(0.0141)	0.074(0.0073)
$\mu$	-0.008(0.0721)	-0.410(0.2935)	NA	-0.404(0.2668)	-1.109(1.4363)	-0.803(1.1190)
Cla_error	0.0116	0.0151	NA	0.0147	0.0201	0.0271
$p=0.5$
$p$	0.000(0.0005)	-0.011(0.0008)	NA	0.017(0.0009)	0.074(0.0066)	0.069(0.0062)
$\mu$	-0.011(0.0244)	-0.169(0.0525)	NA	-0.220(0.0720)	-0.374(0.1801)	-0.607(0.6089)
Cla_error	0.0095	0.0167	NA	0.0178	0.0178	0.0288
$p=0.8$
$p$	0.002(0.0004)	-0.031(0.0014)	NA	0.011(0.0005)	0.042(0.0023)	0.068(0.0054)
$\mu$	-0.007(0.0137)	-0.034(0.0151)	NA	-0.072(0.0190)	-0.159(0.0447)	-0.699(0.5384)
Cla_error	0.0048	0.0082	NA	0.0103	0.0126	0.0360

Table 17. Table 17 : Bias(MSE) of estimates of p / μ 𝑝 𝜇 p/\mu and mean of the classification error for model 5 when n = 250 𝑛 250 n=250 .

	EM_logconcave	Patra	Bordes	Song EM	Song max $\pi$	Xiang
$p=0.2$
$p$	-0.001(0.0007)	0.020(0.0022)	NA	0.030(0.0018)	0.142(0.0240)	0.047(0.0036)
$\mu$	-0.026(0.0943)	-0.679(0.6177)	NA	-0.716(0.6592)	-1.449(2.3985)	-0.895(0.9545)
Cla_error	0.0017	0.0090	NA	0.0087	0.0185	0.0127
$p=0.5$
$p$	0.001(0.0010)	-0.015(0.0014)	NA	0.031(0.0019)	0.116(0.0162)	0.074(0.0075)
$\mu$	0.009(0.0288)	-0.208(0.0734)	NA	-0.267(0.1017)	-0.613(0.4915)	-0.680(0.5574)
Cla_error	0.0011	0.0114	NA	0.0136	0.0202	0.0248
$p=0.8$
$p$	0.000(0.0006)	-0.042(0.0025)	NA	0.015(0.0008)	0.069(0.0056)	0.089(0.0096)
$\mu$	0.006(0.0216)	-0.036(0.0232)	NA	-0.089(0.0305)	-0.296(0.1188)	-0.488(0.3336)
Cla_error	0.0004	0.0037	NA	0.0055	0.0188	0.0544

Table 18. Table 18 : Bias(MSE) of estimates of p / μ 𝑝 𝜇 p/\mu and mean of the classification error for model 5 when n = 500 𝑛 500 n=500 .

	EM_logconcave	Patra	Bordes	Song EM	Song max $\pi$	Xiang
$p=0.2$
$p$	-0.001(0.0003)	0.012(0.0011)	NA	0.023(0.0009)	0.119(0.0176)	0.042(0.0025)
$\mu$	-0.007(0.0438)	-0.500(0.3260)	NA	-0.553(0.3655)	1.237(1.7273)	-0.928(0.9531)
Cla_error	0.0014	0.0058	NA	0.0061	0.0132	0.0107
$p=0.5$
$p$	-0.001(0.0006)	-0.010(0.0007)	NA	0.025(0.0012)	0.120(0.0179)	0.077(0.0073)
$\mu$	0.014(0.0170)	-0.194(0.0565)	NA	-0.233(0.0722)	-0.630(0.4822)	-0.726(0.5840)
Cla_error	0.0008	0.0094	NA	0.0111	0.0188	0.0240
$p=0.8$
$p$	0.001(0.0003)	-0.031(0.0013)	NA	0.013(0.0004)	0.078(0.0067)	0.088(0.0087)
$\mu$	0.006(0.0094)	-0.020(0.0105)	NA	-0.065(0.0146)	-0.324(0.1188)	-0.506(0.3185)
Cla_error	0.0003	0.0026	NA	0.0039	0.0199	0.0501

Table 19. Table 19 : Bias(MSE) of estimates of p / μ 𝑝 𝜇 p/\mu and mean of the classification error for model 6 when n = 250 𝑛 250 n=250 .

	EM_logconcave	Patra	Bordes	Song EM	Song max $\pi$	Xiang
$p=0.2$
$p$	-0.006(0.0015)	0.003(0.0021)	0.020(0.0036)	-0.012(0.0009)	0.109(0.0147)	0.027(0.0023)
$\mu$	0.127(0.1230)	-0.200(0.1172)	-0.685(0.8830)	-0.118(0.0777)	-0.654(0.5415)	-0.067(0.0877)
Cla_error	0.0464	0.0468	0.1738	0.0470	0.0509	0.0473
$p=0.5$
$p$	-0.015(0.0028)	-0.044(0.0037)	0.000(0.0024)	-0.045(0.0031)	0.057(0.0054)	0.031(0.0061)
$\mu$	0.073(0.0525)	0.123(0.0325)	-0.009(0.0298)	0.127(0.0329)	-0.099(0.0405)	-0.068(0.0676)
Cla_error	0.0688	0.0682	0.0676	0.0688	0.0683	0.0738
$p=0.8$
$p$	0.006(0.0016)	-0.077(0.0067)	-0.003(0.0016)	-0.059(0.0044)	0.011(0.0012)	0.006(0.0013)
$\mu$	-0.024(0.0205)	0.173(0.0369)	0.004(0.0104)	0.155(0.0320)	0.032(0.0125)	-0.002(0.0105)
Cla_error	0.0591	0.0660	0.0595	0.0647	0.0588	0.0572

Table 20. Table 20 : Bias(MSE) of estimates of p / μ 𝑝 𝜇 p/\mu and mean of the classification error for model 6 when n = 500 𝑛 500 n=500 .

	EM_logconcave	Patra	Bordes	Song EM	Song max $\pi$	Xiang
$p=0.2$
$p$	-0.008(0.0006)	-0.001(0.0008)	0.016(0.0021)	-0.017(0.0007)	0.095(0.0119)	0.011(0.0007)
$\mu$	0.132(0.0565)	-0.066(0.0452)	-0.161(0.1613)	0.000(0.0312)	-0.538(0.3834)	-0.025(0.0370)
Cla_error	0.0435	0.0457	0.0451	0.0450	0.0479	0.0439
$p=0.5$
$p$	-0.011(0.0014)	-0.033(0.0019)	-0.001(0.0010)	-0.049(0.0030)	0.043(0.0032)	0.011(0.0017)
$\mu$	0.069(0.0304)	0.137(0.0272)	-0.004(0.0109)	0.161(0.0339)	-0.038(0.0183)	-0.023(0.0203)
Cla_error	0.0655	0.0676	0.0666	0.0684	0.0665	0.0664
$p=0.8$
$p$	0.004(0.0005)	-0.065(0.0047)	-0.002(0.0008)	-0.062(0.0043)	0.010(0.0011)	0.002(0.0007)
$\mu$	-0.021(0.0069)	0.169(0.0325)	-0.002(0.0054)	0.165(0.0311)	0.043(0.0102)	0.002(0.0060)
Cla_error	0.0541	0.0646	0.0576	0.0642	0.0574	0.0546

Equations

$$g(x)=(1-p)f_0(x)+pf(x),$$

$$P_i=\frac{pf(t_i)}{(1-p)f_0(t_i)+pf(t_i)}.$$

$$g(x)=(1-p){\varphi }_{\sigma }(x)+pf(x)$$

$$g(x)=(1-p)f_0(x;\xi )+pf(x-\mu ),$$

$$\underset{x\to a^{+}}{\lim }\frac{f(x)}{f_0(x)}=0\text{or}\underset{x\to a^{-}}{\lim }\frac{f(x)}{f_0(x)}=0.$$

$$\text{log}(f_0(x))=\text{log}(\Gamma (\frac{\nu +1}{2}))-\frac{1}{2}\text{log}(\nu \pi )-\text{log}(\Gamma (\frac{\nu }{2}))-\frac{\nu +1}{2}\text{log}(1+\frac{x^2}{\nu }).$$

$$\varphi (x)-\text{log}{f}_0(x)$$

$$L(p,\varphi ,Q_n)=n\int \text{log}(g)d{Q}_n=\sum _{i=1}^n\text{log}((1-p)f_0(X_i)+p{e}^{\varphi (X_i)}),$$

$${\omega }_i^{(k+1)}=\frac{(1-p^{(k)})f_0(x_i)}{(1-p^{(k)})f_0(x_i)+p^{(k)}{f}^{(k)}(x_i)},$$

$$p^{(k+1)}=\frac{1}{n}\sum _{i=1}^n(1-{\omega }_i^{(k+1)}),$$

$${\varphi }^{(k+1)}=\stackrel{[}{\varphi }\in {\Phi }^1,\int e^{\varphi (x)}dx=1]\text{argmax}\sum _{i=1}^n(1-{\omega }_i^{(k+1)})\varphi (x_i),$$

$$f^{(k+1)}=e^{{\varphi }^{(k+1)}}.$$

$$l^{(k+1)}\ge l^{(k)},$$

$$L(p,\varphi ,Q)$$

$$\text{csupp}(Q)=\bigcap \{C:C\subseteq R^d\text{closed and convex},Q(C)=1\}.$$

$$\int ||x|{|}^{k}Q(dx)<\infty \text{and}\text{interior}(\text{csupp}(Q))\ne \varnothing .$$

$$L(Q)=\stackrel{[}{p\in [0,1],\varphi \in {\stackrel{~}{\Phi }}^d}]\text{sup}L(p,\varphi ,Q)$$

$$(p_0,{\varphi }_0)\in \stackrel{[}{p\in [0,1],\varphi \in {\stackrel{~}{\Phi }}^d}]\text{argmax}L(p,\varphi ,Q).$$

$$\text{interior}(\text{csupp}(Q))\subseteq \text{dom}({\varphi }_0)=\{x\in R^d:{\varphi }_0(x)>-\infty \}\subseteq \text{csupp}(Q).$$

$$\int \text{log}(g_0(x)/g_1(x))g_0(x)dx=0.$$

$${\mathcal{Q}}^k$$

$${\mathcal{Q}}_0$$

$$D_k(Q,Q^{\prime })=\stackrel{[}{\stackrel{[}{X}\sim Q,X^{\prime }\sim Q^{\prime }]X,X^{\prime }}]\text{inf}\{E||X-X^{\prime }|{|}^k{\}}^{1/k}.$$

$$\stackrel{[}{n}\to \infty ]\text{lim}L(Q_n)=L(Q).$$

$$\stackrel{[}{n}\to \infty ]\text{lim}{p}_n$$

$$\stackrel{[}{n}\to \infty ,x\to y]\text{lim}{f}_n(x)$$

$$\stackrel{[}{n}\to \infty ,x\to y]\text{limsup}{f}_n(x)$$

$$\stackrel{[}{n}\to \infty ]\text{lim}\int |f_n(x)-f^{\ast }(x)|dx$$

$$\widehat{w_i}=\frac{(1-\hat{p})f_0(x_i)}{(1-\hat{p})f_0(x_i)+\hat{p}\hat{f}(x_i)}.$$

$$\widehat{w_i}=\frac{2(1-\hat{p})f_0(x_i)}{(1-\hat{p})f_0(x_i)+\hat{h}(x_i)},$$

$$\hat{\mu }=\frac{{\sum }_{i=1}^n(1-\widehat{w_i})X_i}{{\sum }_{i=1}^n(1-\widehat{w_i})}.$$

$$\text{Cla\_error}=\frac{1}{n}\sum _{i=1}^n(\widehat{w_i}-w_i{)}^2,$$

Full text

1 Introduction

In this paper, we consider the following two-component mixture model,

[TABLE]

where the probability density function (pdf) $f_{0}(x)$ is known, whereas the mixing proportion $p\in[0,1]$ and the pdf $f$ are unknown. Model (1.1) is motivated by studies in the biological sciences to cluster differentially expressed genes in microarray data Bordes et al., (2006). Typically we build a test statistic, say $T_{i}$, for each gene $i$. Under the null hypothesis, which presumes no difference in expression levels under two or more conditions, $T_{i}$ is assumed to have a known distribution (in general Student’s or Fisher). Under the alternative hypothesis, the distribution is unknown. Thus, the distribution of the test statistic is modelled by (1.1) where $p$ is the proportion of non-null statistics. The estimation of $p$ and the pdf $f$ can tell us the probability $P_{i}$ that gene $i$ is differentially expressed given $T_{i}=t_{i}$:

[TABLE]

Bordes et al., (2006) considered model (1.1) where $f$ is assumed to be symmetric. They obtained some identifiability results under moment and symmetry conditions, and proposed to estimate this model under the symmetry of $f$. In addition, they proved the consistency of their estimator under mild conditions.

Song et al., (2010) considered another special case,

[TABLE]

where $f_{0}=\phi_{\sigma}$ is a normal density with mean [math] and unknown standard deviation $\sigma$. This model was inspired by sequential clustering (Song & Nicolae,, 2009), which finds candidates for centers of clusters first, then carries out a local search to find the objects that belong to those clusters, and finally selects the best cluster. Song et al., (2010) proposed an EM-type estimator and a maximizing $\pi$-type estimator for their model which can be easily extended to models where $f_{0}$ is not normal.

A slightly different model is considered by Xiang et al., (2014),

[TABLE]

where $\xi$ is a possibly unknown parameter, and $\mu$ is a non-null location parameter for $f$. They proposed a new effective estimator based on the minimum profile Hellinger distance (MPHD). They established the existence and uniqueness of their estimator and also proved its consistency under some regularity conditions. Their method actually does not require $f$ to be symmetric and thus can be applied to a more general model. For some other alternative estimators, see, for example Patra & Sen, (2016); Ma et al., (2015).

In this paper, we propose to estimate (1.1) using a new approach by imposing a log-concave assumption on $f$, i.e. $\text{log}(f)\in\Phi^{1}$; here $\Phi^{d}$ denotes the family of concave functions $\phi$ on $R^{d}$ which are upper semicontinuous and coercive in the sense that $\phi(x)\rightarrow-\infty,\ \text{as}\ ||x||\rightarrow\infty$. Note that $\text{log}(f)$ needs to be coercive in order for $f$ to be a density function. Many common parametric families of distributions belong to the family of log-concave densities, for example, normal distribution, exponential distribution, logistic distribution, etc. We propose to estimate the new model by maximizing a semiparametric mixture likelihood. Compared to the kernel density estimation of $f$ used by many existing methods (Bordes et al.,, 2006; Xiang et al.,, 2014; Ma et al.,, 2015), the new method does not require the choice of one or more bandwidths (Samworth, (2017)). We establish the identifiability of the proposed semiparametric mixture model and prove the existence and consistency of the proposed estimators. We further compare our estimator with several existing estimators through simulation studies and apply our method to two real data sets from biological sciences and astronomy.

The rest of the paper is organized as follows. In Section 2 we discuss some identifiability issues for model (1.1). Section 3 introduces our maximum likelihood estimator and a detailed EM type algorithm is provided. Existence and consistency properties of our estimator are established with detailed proofs given in the Appendix (Section 7). Section 4 demonstrates the finite sample performance of our proposed estimator by comparing with many other existing algorithms. Two real data applications are given in Section 5. Section 6 gives a brief discussion.

2 Identifiability

Note that the model (1.1) is non-identifiable without any constraint on the density $f$, see e.g., Bordes et al., (2006), and Patra & Sen, (2016). However a parametric model for $f$ might create biased or even misleading statistical inference when the model assumption is incorrect. In this paper, we assume $f(x)$ to be log-concave, i.e. $f(x)=e^{\phi(x)}$, where $\phi(x)$ is a concave function. Log-concave densities attracted lots of attention in the recent years since it is very flexible and can be estimated by nonparametric maximum likelihood estimator without requiring the choice of any tuning parameter. For more details, see Cule et al., (2010a), Dümbgen et al., (2009), Walther et al., (2009), Dümbgen et al., (2011) and the review of the recent progress in log-concave density estimation by Samworth, (2017).

Proposition 2.1**.**

Assume $f_{0}>0$ around a neighborhood of $a$, then model (1.1) is identifiable if,

[TABLE]

Remark 2.1**.**

Proposition 2.1 also holds if $a=\pm\infty$, and this result is much more general (require much weaker condition) than the result of Proposition 3(i) of Bordes et al., (2006).

Remark 2.2**.**

Proposition 2.1 guarantees that model (1.1) is identifiable if the support of $f$ is strictly contained in the support of $f_{0}$ and the two supports have different Legesgue measure.

Now, for a log-concave density, it is easy to derive the following more specific result,

Proposition 2.2**.**

Assume $f_{0}>0$ and $\text{log}(f)\in\Phi^{1}$. Model (1.1) is identifiable if $\textup{log}f_{0}(x)=O(x^{k})$, for some $0<k<1$, as $x\rightarrow+\infty$ or $x\rightarrow-\infty$.

Example 2.1**.**

If $f_{0}(x)$ is the density of a $t$ distribution with $\nu$ degrees of freedom, then model (1.1) is identifiable.

Proof.

Since $f_{0}(x)=\frac{\Gamma(\frac{\nu+1}{2})}{\sqrt{\nu\pi}\Gamma(\frac{\nu}{2})}(1+\frac{x^{2}}{\nu})^{-\frac{\nu+1}{2}}$, we have,

[TABLE]

Thus for any $0<k<1$, $\text{log}(f_{0}(x))/x^{k}\rightarrow 0$, as $x\rightarrow+\infty$. And by Proposition 2.2, we can conclude that model (1.1) is identifiable when $\text{log}(f)\in\Phi^{1}$. ∎

Remark 2.3**.**

Similarly, one can check that when $f_{0}$ is the pdf of an $F$ distribution, log-normal distribution, or Pareto distribution, then model (1.1) is identifiable under the condition that $\text{log}(f)\in\Phi^{1}$.

In general, we need $\phi(x)-\text{log}f_{0}(x)\rightarrow-\infty$ as $x\rightarrow+\infty$ or $x\rightarrow-\infty$ to make model (1.1) identifiable. Since $f_{0}(x)$ is known, it is not difficult to give some sufficient conditions with respect to different $f_{0}$’s. Here we give two examples of easy to apply conditions.

Example 2.2**.**

Suppose $f_{0}(x)$ is the density of a normal distribution with mean $\mu$ and variance $\sigma^{2}$, then model (1.1) is identifiable if $\lim_{x\rightarrow+\infty}\frac{\phi(x)}{x^{2}}<-\frac{1}{2\sigma^{2}}$, or $\lim_{x\rightarrow-\infty}\frac{\phi(x)}{x^{2}}<-\frac{1}{2\sigma^{2}}$, or the condition of Remark 2.2 holds.

Proof.

Suppose $\lim_{x\rightarrow+\infty}\frac{\phi(x)}{x^{2}}<-\frac{1}{2\sigma^{2}}$, or $\lim_{x\rightarrow-\infty}\frac{\phi(x)}{x^{2}}<-\frac{1}{2\sigma^{2}}$. Since

[TABLE]

Hence $\displaystyle\frac{f(x)}{f_{0}(x)}\rightarrow 0$ as $x\rightarrow+\infty$ or $x\rightarrow-\infty$, and Proposition 2.1 asserts the identifiability of model (1). ∎

Remark 2.4**.**

Under the constraints set by Example 2.2, Model 1, 4, and 5 from Section 4 are identifiable.

Example 2.3**.**

Suppose $f_{0}(x)$ is the density of an exponential distritution with rate $\lambda$, then model (1.1) is identifiable if $\lim_{x\rightarrow+\infty}\frac{\phi(x)}{x}<-\lambda$, or the condition of Remark 2.2 holds.

Proof.

Suppose $\lim_{x\rightarrow+\infty}\frac{\phi(x)}{x}<-\lambda$. Since,

[TABLE]

Hence $\displaystyle\frac{f(x)}{f_{0}(x)}\rightarrow 0$ as $x\rightarrow+\infty$, and again Proposition 2.1 ensures the identifiability of model (1.1). ∎

Remark 2.5**.**

Under the constraints set by Example 2.3, Model 3 from Section 4 is identifiable.

3 Maximum Likelihood Estimation

Suppose we have a random sample of $n$ i.i.d. observations $X_{1},X_{2},\cdots,X_{n}$ from the density $g(x)=(1-p)f_{0}(x)+pf(x)$, $p\in[0,1]$, and $f=e^{\phi}$ is a log-concave density, i.e., $\phi\in\Phi^{1}$. Then with the empirical distribution $Q_{n}=\displaystyle\frac{1}{n}\sum_{i=1}^{n}\delta_{X_{i}}$ , where $\delta_{X_{i}}$ is the degenerate distribution function at $\{X_{i}\}$, the log likelihood of our random sample can be written as,

[TABLE]

subject to the condition that $\int e^{\phi(x)}dx=1$. A natural approach to estimate $p$ and $\phi$ is to find the maximizer of $L(p,\phi,Q_{n})$.

3.1 Algorithm

We propose to estimate $p$ and $f$ by maximizing $L(p,\phi,Q_{n})$, using the EM algorithm that consists of iterating an E step and M step until convergence:

“E step”: Given $p^{(k)}$ and $f^{(k)}$,

[TABLE]

“M step”:

[TABLE]

We find $\phi^{(k+1)}$ using an active set algorithm, which is described in Dümbgen et al., (2007) and implemened in the R package logcondens by Rufibach & Duembgen, (2010). Through out this paper, we use “EM_logconcave” to represent our method. The following result establishes the monotone properties of EM_logconcave algorithm.

Proposition 3.1**.**

Let $l^{(k)}=\stackrel{{\scriptstyle[}}{{i}}=1]{n}{\sum}\textup{log}((1-p^{(k)})f_{0}(x_{i})+p^{(k)}e^{\phi^{(k)}(x_{i})})$, then

[TABLE]

for any $k\geq 0$.

3.2 Theoretical Properties

In general, for any distribution $Q$ on $R^{d}$, we define,

[TABLE]

For the existence of a maximizer of $L(p,\phi,Q)$, we follow the approach of Dümbgen et al., (2011). We define the convex support of $Q$ as,

[TABLE]

Theorem 3.1**.**

For fixed $f_{0}$, assume $\textup{supp}\{f_{0}\}\subseteq\textup{csupp}(Q)$, and there exist some integer $k\geq 1$, such that,

[TABLE]

Let $\tilde{\Phi}^{d}=\{\phi\in\Phi^{d}:\int e^{\phi(x)}dx=1\ \text{and}\ f_{0}(x)\leq m(x)e^{\phi(x)}\ \text{for \ some \ }m(x)=c_{0}e^{c_{1}||x||^{k}},c_{0}\geq 0,c_{1}\geq 0\}$, then

[TABLE]

is real. In that case, there exists,

[TABLE]

Moreover,

[TABLE]

The proof of Theorem 3.1 is given in the Appendix (Section 7).

Example 3.1**.**

Assume $Q$ represents the distribution of Model 1: $g(x)=(1-p)*N(\mu=0,\sigma=2)+p*N(\mu=3,\sigma=1)$, hence $f_{0}$ represents the pdf of $N(\mu=0,\sigma=2)$ distribution. For any integer $k\geq 2$, $\tilde{\Phi}^{1}$ contains all the pdfs of normal distribution, logistic distribution, and Laplace distribution, etc. And Theorem 3.1 implies that the maximum of $L(p,\phi,Q)$ exists over $p\in[0,1]$ and $\phi\in\tilde{\Phi}^{1}$.

In general, the maximizer of $L(p,\phi,Q)$ is not unique. But if $Q$ has density $g_{0}(x)=(1-p_{0})f_{0}(x)+p_{0}e^{\phi_{0}(x)}$, where $g_{0}(x)$ is identifiable, then $L(p_{0},\phi_{0},Q)=\int\text{log}(g_{0}(x))g_{0}(x)dx$, and this $(p_{0},\phi_{0})$ is the unique maximizer. This is because as noted by Dümbgen et al., (2011), if we have $(p_{1},\phi_{1})$, such that $L(Q)=L(p_{0},\phi_{0},Q)=L(p_{1},\phi_{1},Q)$, let $g_{1}(x)=(1-p_{1})f_{0}(x)+p_{1}e^{\phi_{1}(x)}$, then

[TABLE]

Note the above integral is exactly the Kullback-Leibler divergence which is positive and equals [math] iff $g_{0}=g_{1}$ almost everywhere. Thus $(p_{0},\phi_{0})=(p_{1},\phi_{1})$ except that $\phi_{0}$ and $\phi_{1}$ may differ on a set of Lebesgue measure zero.

Next we establish the consistency of our maximum likelihood estimator. First, we introduce some notations,

[TABLE]

In what follows, we consider the convergence of distributions under Mallows’ distance $D_{1}$(Mallows, (1972)). Specifically, for two distributions $Q$, $Q^{\prime}\in\mathcal{Q}^{k}$,

[TABLE]

It is known that $\stackrel{{\scriptstyle[}}{{n}}\rightarrow\infty]{}{\text{lim}}D_{k}(Q_{n},Q)\rightarrow 0$ is equivalent to $Q_{n}\rightarrow_{w}Q$ and $\int||x||^{k}Q_{n}(dx)\rightarrow\int||x||^{k}Q(dx)$(Bickel & Freedman, (1981); Mallows, (1972)). Here $Q_{n}\rightarrow_{w}Q$ means weak convergence, or convergence in distribution.

Now we are ready to state our consistency theorem.

Theorem 3.2**.**

Assume, (a). $\textup{supp}\{f_{0}\}\subseteq\textup{csupp}(Q)$; (b). for some fixed integer $k\geq 1$, the unknown density $f$ satisfies the following condition: $\exists m(x)=c_{0}e^{c_{1}||x||^{k}}$, where $c_{i}\geq 0$, $i=0,1$, such that, $f_{0}(x)\leq m(x)f(x)=m(x)e^{\phi(x)}$. Let $\{Q_{n}\}$ be a sequence of distributions in $\mathcal{Q}_{0}\bigcap\mathcal{Q}^{k}$ such that $\stackrel{{\scriptstyle[}}{{n}}\rightarrow\infty]{}{\textup{lim}}D_{k}(Q_{n},Q)=0$ for some $Q\in\mathcal{Q}_{0}\bigcap\mathcal{Q}^{k}$. Suppose $f_{0}$ is upper semi-continuous and $\displaystyle\frac{\text{log}(f_{0})}{1+||x||}$ is bounded. Then

[TABLE]

Assume there exist maximizers $(p_{n},\phi_{n})$ of $L(p,\phi,Q_{n})$, and a unique maximizer $(p^{*},\phi^{*})$ of $L(p,\phi,Q)$, where $p_{n},p^{*}\in[0,1],\phi_{n},\phi^{*}\in\tilde{\Phi}^{d}$. Let $f_{n}=exp(\phi_{n})$, $f^{*}=exp(\phi^{*})$, then

[TABLE]

Practically, $Q_{n}$ will be the empirical distribution function which automatically satisfies the above assumption.

4 Simulation

In this section, we investigate the finite sample performance of our algorithm and compare it to the estimator proposed by Patra & Sen, (2016) ($\hat{\alpha}_{0}^{0.1k_{n}}$ from their paper), the Symmetrization estimator by Bordes et al., (2006), the EM-type estimator, Maximizing-$\pi$ type estimator by Song et al., (2010), and the Minimum profile Hellinger distance estimator by Xiang et al., (2014).

In order to test our method under different settings, we simulate $K=200$ samples of $n$ i.i.d. random variables with the common distribution given by the following six models:

• •

Model 1: $g(x)=(1-p)*N(\mu=0,\sigma=2)+p*N(\mu=3,\sigma=1)$,

• •

Model 2: $g(x)=(1-p)*\text{unif}(0,1)+p*\text{beta}(\alpha=1,\beta=5)$,

• •

Model 3: $g(x)=(1-p)*\text{exp}(\lambda=1)+p*(\text{exp}(\lambda=1)+2)$,

• •

Model 4: $g(x)=(1-p)*N(0,1)+p*(\chi^{2}(3)+2)$,

• •

Model 5: $g(x)=(1-p)*N(0,1)+p*(\text{exp}(\lambda=0.5)+3)$,

• •

Model 6: $g(x)=(1-p)*N(0,1)+p*(t(d.f.=5)+3)$.

For each sample we estimate $p$, the mean ($\mu$) of the unknown component $f$ and the classification error. The detailed calculation is explained bellow.

For our algorithm and the algorithm by Xiang et al., (2014), final estimators $\hat{p}$ and $\hat{f}$ are always produced, thus the estimated probability $\hat{w}_{i}$ that the $i$-th observation is from the known component $f_{0}(x)$, given $X_{i}=x_{i}$, can be calculated by

[TABLE]

For other methods, $\hat{f}$ may not always be given directly. Suggested by Song et al., (2010), we estimate $\hat{w_{i}}$ by the following,

[TABLE]

where $\hat{h}$ is the kernel density estimator of $g$ with Gaussian kernel and Silverman’s “rule of thumb” bandwidth (Silverman,, 1986). Note that the algorithm proposed by Patra & Sen, (2016) actually can estimate $\hat{f}$ when $f$ is non-increasing. But we find the algorithm works best when $f_{0}$ and $f$ has the same support and often produces unreliable estimates when two supports differ from each other. Thus, we do not use $\hat{f}$ to estimate $\hat{w_{i}}$ for Patra & Sen, (2016)’s algorithm even when the true $f$ does decrease on its support, instead, we follow Song et al., (2010)’s recommendation to get $\hat{w_{i}}$.

The algorithms by Xiang et al., (2014) and Bordes et al., (2006) give a final mean estimator $\hat{\mu}$ directly. For other methods, after we get $\hat{w_{i}}$, we estimate $\mu$ by the following weighted sum,

[TABLE]

Last, we report the classification error (Cla_error) based on $\hat{w_{i}}$ as the mean squared error between $\hat{w_{i}}$ and the true $w_{i}$, i.e.,

[TABLE]

where $w_{i}=1$ if $x_{i}$ is from the known component $f_{0}(x)$ and [math] if $x_{i}$ is from the unknown component $f(x)$.

For model 1, Table 1 reports the bias and MSE of the estimates of $p$, the bias and MSE of the estimates of $\mu$, and the mean of the classification error for different methods over $K=200$ repetitions when $p=0.2$, $p=0.5$, and $p=0.8$, with sample size $n=1000$. Similar reports of other models can be found in Tables 1 — 6. Simulation results for sample sizes $n=250$ and $n=500$ are reported in the Appendix(Section 7). We report the results of Bordes et al., (2006)’s algorithm only for model 1, model 2 and model 6, since this method fails to estimate $p$ for other models, in which the real $f(x)$’s are not symmetric on their supports.

All the simulation results strongly suggest that our method is very competitive and often outperforms all other methods. Moreover, our method is even more favorable when the sample size $n$ gets larger.

To better display our simulation results, we also plot the MSE of point estimates of $p$ and $\mu$ vs. different models for all the methods we mentioned above when $p=0.2$ and $n=1000$, except for the method by Bordes et al., (2006) as their method fails to estimate $p$ and $\mu$ for half of the models we discussed here. Figure 1 shows that the curve representing our method always lies at the bottom which demonstrates the effectiveness of our algorithm, while the Maximizing-$\pi$ type estimator by Song et al., (2010) gives the worst results in terms of MSE.

5 Real Data Application

5.1 Prostate Data

In this section we consider the prostate data consisting of genetic expression levels related to prostate cancer patients of Efron, (2012). The data set is a $6033\times 102$ matrix, with entries $x_{ij}=$ expression level for gene $i$ on patient $j$, $i=1,\cdots,n$, $j=1,\cdots,m$, here $n=6033$, $m=102$. Among the $m=102$ patients, $m_{1}=50$ of them are normal control subjects (corresponding to $j=1,\cdots,m_{1}$) and $m_{2}=52$ of them are prostate cancer patients (corresponding to $j=m_{1}+1,\cdots,m_{2}$). The goal of the study is to discover the potential genes that are differentially expressed between normal control and prostate cancer patients.

Two-sample $t$-test is performed to test the significance of each gene $i$ by,

[TABLE]

where $\bar{x}_{i}(1)=(\stackrel{{\scriptstyle[}}{{j}}=1]{m_{1}}{\sum}x_{ij})/m_{1}$, $\bar{x}_{i}(2)=(\stackrel{{\scriptstyle[}}{{j}}=m_{1}+1]{m_{2}}{\sum}x_{ij})/m_{2}$, $s_{i}^{2}=(1/m_{1}+1/m_{2})\left\{\stackrel{{\scriptstyle[}}{{j}}=1]{m_{1}}{\sum}(x_{ij}-\bar{x}_{i}(1))^{2}+\right.$ $\left.\stackrel{{\scriptstyle[}}{{j}}=m_{1}+1]{m_{2}}{\sum}(x_{ij}-\bar{x}_{i}(2))^{2}\right\}/(m-2)$. These two-sided $t$-tests produce $n=6033$ $p$-values, and the distribution of these $p$-values under the null hypothesis (i.e., gene $i$ is not differentially expressed) has a uniform density, while under the alternative hypothesis (i.e., gene $i$ is differentially expressed) has a non-increasing density.

The estimation of $p$ is reported in Table 7. We can see that the estimate by Bordes et al., (2006) and the Maximizing-$\pi$ type estimate by Song et al., (2010) give a relatively big estimate. The estimate procedure by Bordes et al., (2006) assumes the density function under the alternative hypothesis to be symmetric, while in our example this density is non-increasing, which implies the violation of their symmetric assumption. It is known that the Maximizing-$\pi$ type estimator by Song et al., (2010) tends to overestimate the $p$ value, which can also been seen in Table 7. We also want to point out that several approaches have been proposed by Efron, (2012) to estimate $p$ as well, the estimator based on central matching method gives $\hat{p}=0.020$ (please see Efron, (2012) and Efron et al., (2007) for detailed description of those estimators), and Table 7 shows that our estimator gives a closest value to Efron’s result.

Figure 2 shows that our estimate of the density $\hat{f}$ under the alternative tends to have a much smaller support comparing to the one given by Patra & Sen, (2016). Again, as we noted before, the method by Patra & Sen, (2016) assumes that $f$ is decreasing on the whole support of $f_{0}$. While in reality, smaller $p$-values tends to indicate the alternative hypothesis, hence it actually makes sense that the support of $f$ for this prostate data may be much smaller than $(0,1)$. The estimate produced by Bordes et al., (2006) is not very reliable, since the density $f$ is not symmetric. For this example, if we apply Bordes et al., (2006)’s method to the original $t$ statistics directly, the estimate is $\hat{p}=0.0072$.

5.2 Carina Data

Carina is one of the seven dwarf spheroidal (dSph) satellite galaxies of Milkey Way. Here we consider the data consisting of radial velocities (RV) of $n=1266$ stars from Carina galaxy. The data is obtained by Magellan and MMT telescopes (Walker et al., (2007)). The stars of Milkey Way contribute contamination to this data set. We assume the distribution $f_{0}$ of RV from stars of Milkey Way is known from the Besancon Milky Way model (Robin et al., (2003)). Now we would like to analyze this data set to better understand the mixture distribution of the RV of stars in Carina galaxy.

The estimation of $p$ is reported in Table 8. Again we see that the estimation by Song et al., (2010)’s Maximizing-$\pi$ type estimator gives a relatively big estimate. Other estimates are relatively close.

Next in Figure 3 we plot the histogram of the RV data overlaid with our estimated two components of the mixture density, and we can see that our estimation approximates the real data fairly well. The component corresponding to the stars of Carina looks very symmetric, and in fact astronomers usually assume the distribution to be Gaussian, which causing the density estimation proposed by Patra & Sen, (2016) does not work here.

6 Discussion

In this paper we study the two-component mixture model with one component completely known. A nonparametric maximum likelihood estimator is developed via EM algorithm and log-concave approximation. Unlike most existing estimation procedures, our new method finds the mixing proportion and the distribution of the unknown component simultaneously without any selection of a tuning parameter and the proposed EM algorithm satisfies the non-decreasing property of the traditional EM algorithm. Simulation results show that our method is more favorable than many other competing estimation methods.

We are able to prove the existence and consistency of our maximum likelihood estimator for a general distribution $Q$. But we do require some extra conditions on $f_{0}$ and $f$. A possible future research direction is trying to ease these assumptions and make it more general. In addition, it would be also our interest to apply our method to a more general model where the component $f_{0}$ also contains some unknown parameter.

7 Appendix

7.1 Theoretical Proof

Proof of Proposition 2.1.

According to Patra & Sen, (2016), if we let $G,\ F_{0}$, and $F$ be the cumulative distribution functions of $g,\ f_{0}$ and $f$ respectively, define $p_{0}=\text{inf}\{\gamma\in(0,1]:[G-(1-\gamma)F_{0}]/\gamma\ \text{is a CDF}\}$, then

[TABLE]

where $\text{essinf}(h)=\text{sup}\{t\in R:\mathfrak{m}\{x:h(x)<t\}=0\}$, and here $\mathfrak{m}$ represents the Lebesgue measure. Now if $\text{essinf}\frac{f}{f_{0}}>0$, there must exist some $t>0$, such that, $\mathfrak{m}\{x:\frac{f(x)}{f_{0}(x)}<t\}=0$, i.e., $\frac{f(x)}{f_{0}(x)}\geq t$ almost everywhere, which contradicts to the fact that $\lim_{x\rightarrow a^{+}}\frac{f(x)}{f_{0}(x)}=0\ \text{or}\ \lim_{x\rightarrow a^{-}}\frac{f(x)}{f_{0}(x)}=0$. Hence we can conclude that $\text{essinf}\frac{f}{f_{0}}=0$, and consequently $p_{0}=p$, which means if we can write $g(x)=(1-p)f_{0}(x)+pf(x)$, this $p$ is fixed and equals $p_{0}$. Consequently $f(x)=(g(x)-(1-p)f_{0}(x))/p$ is fixed as well, and our model (1.1) is identifiable. ∎

Proof of Proposition 2.2.

Since $\displaystyle f(x)=e^{\phi(x)}$ is a log-concave density, there exist constants $a$ and $b>0$, such that $\phi(x)\leq a-b|x|$ (see Cule et al., (2010b)), which implies

[TABLE]

Now if $\text{log}f_{0}(x)=O(x^{k})$, for some $0<k<1$, apparently,

[TABLE]

Hence $\phi(x)-\text{log}f_{0}(x)\rightarrow-\infty$ as $x\rightarrow+\infty$ or $x\rightarrow-\infty$, which shows $\lim_{x\rightarrow+\infty}\frac{f(x)}{f_{0}(x)}=0$. Thus, model (1.1) is identifiable from Proposition 2.1. ∎

Proof of Theorem 3.1.

Suppose $\int||x||^{k}Q(dx)<\infty$, $\text{interior}(\text{csupp}(Q))\neq\emptyset$, $\int e^{\phi(x)}dx=1$ and $f_{0}(x)\leq m(x)e^{\phi(x)}$. For any concave function $\phi$ satisfying the above conditions, there exist $(a_{0},b_{0})$, such that $\phi(x)\leq a_{0}-b_{0}||x||$, thus for any $p\neq 0$, $L(p,-b_{0}||x||-\text{log}(\int e^{-b_{0}||x||}dx),Q)\geq\text{log}\frac{p}{\int e^{-b_{0}||x||}dx}-b_{0}\int||x||Q(dx)>-\infty$, thus we have $L(Q)>-\infty$. When maximizing $L(p,\phi,Q)$ over all $\phi\in\tilde{\Phi}^{d}$, we may restrict our attention to functions $\phi$ such that $\text{dom}(\phi)=\{x\in R^{d}:\phi(x)>-\infty\}\subseteq\text{csupp}(Q)$. For if $\text{dom}(\phi)\nsubseteq\text{csupp}(Q)$, replacing $\phi(x)$ with $-\infty$ for all $x\notin\text{csupp}(Q)$, then the value of $L(p,\phi-\text{log}(\int e^{\phi(x)}dx),Q)$ would be greater or equal to the original $L(p,\phi,Q)$. Note that since $\text{csupp}(f_{0})\subseteq\text{csupp}(Q)$, the new concave function $\phi^{\prime}=\phi-\text{log}(\int e^{\phi(x)}dx)$ still satisfies the conditions above, i.e., $\int e^{\phi^{\prime}(x)}dx=1$, $f_{0}(x)\leq m(x)e^{\phi^{\prime}(x)}$ and $\text{dom}(\phi^{\prime})=\{x\in R^{d}:\phi(x)>-\infty\}\subseteq\text{csupp}(Q)$. We denote $\Phi(Q)$ to be the family of all $\phi\in\Phi^{d}$ with these properties.

Now we show that $L(Q)<\infty$. Suppose that $\phi\in\Phi(Q)$ is such that $M=\text{max}_{x\in R^{d}}\phi(x)>0$. Let $D_{t}=\{\phi\geq t\}$, hence $D_{t}$ is closed and convex. For any $\alpha>0$, we have the following estimate,

[TABLE]

Note that $\int\text{log}(m(x)+1)Q(dx)$ exists since $\int||x||^{k}Q(dx)<\infty$. By Lemma 4.1 of Dümbgen et al., (2011), for any fixed $\alpha$,

[TABLE]

Lemma 2.1 of Dümbgen et al., (2011) says that for sufficiently large $\alpha$ and sufficiently small $\delta>0$, there exist some sufficiently small $\epsilon>0$, such that,

[TABLE]

which implies that $L(p,\phi,Q)\rightarrow-\infty$, as $M\rightarrow\infty$. Since for any $\phi\in\Phi(Q)$, we also have $L(p,\phi,Q)\leq\int\text{log}((1-p)m(x)+p)dQ+M$, $\Rightarrow L(Q)<\infty$ and there exist constants $M_{0}$ and $M_{*}$, such that

$L(Q)=\stackrel{{\scriptstyle[}}{{\stackrel{{\scriptstyle[}}{{M}}_{0}\leq\text{max}(\phi(x))\leq M_{*}]{}{p\in[0,1],\ \phi\in\Phi(Q)}}}]{}{\text{sup}}L(p,\phi,Q)$.

Now that we know $L(Q)$ is real, we are ready to prove the existence of a maximizer $(p_{0},\phi_{0})$ of $L(Q)$. Let $(p_{n},\phi_{n})$ be a sequence such that $p_{n}\in[0,1]$, $\phi_{n}\in\Phi(Q)$, $M_{n}=\text{max}(\phi_{n}(x))\in[M_{0},M_{*}]$, and $-\infty<L(p_{n},\phi_{n},Q)\uparrow L(Q)$ as $n\rightarrow\infty$. Here we assume $\{p_{n}\}$ is a convergent sequence, say $p_{n}\rightarrow p_{0}\in[0,1]$, as $n\rightarrow\infty$. If $\{p_{n}\}$ is not convergent, since it is bounded, it must have a convergent subsequence $\{p_{n_{k}}\}$, and the sequence $\{p_{n_{k}},\phi_{n_{k}}\}$ would satisfy all those properties above and we can just simply replace the original sequence with this subsequence.

Next, we show that,

[TABLE]

For any $x_{0}\in\text{interior}(\text{csupp}(Q))$, if $\phi_{n}(x_{0})<M_{n}$, then $x_{0}$ can not be an interior point of $\{\phi_{n}\geq\phi_{n}(x_{0})\}$, hence,

[TABLE]

where $h(Q,x)=\text{sup}\{Q(C):C\subseteq R^{d}\ \text{closed and convex},\ x\notin\text{interior}(C)\}<1$ by Lemma 2.13 of Dümbgen et al., (2011). And the above inequalities still hold even if $\phi_{n}(x_{0})=M_{n}$. Thus we have,

[TABLE]

which establishes (7.1). Since $\phi_{n}\leq M_{*}$, together with (7.1), Lemma 3.3 of Schuhmacher et al., (2011) implies that there exist constants $a$ and $b>0$ such that,

[TABLE]

Let $C=\{x\in R^{d}:\stackrel{{\scriptstyle[}}{{n}}\rightarrow\infty]{}{\text{liminf}}\phi_{n}(x)>-\infty\}\supseteq\text{interior}(\text{csupp}(Q))$ and $\bar{\phi}(x)=a-b||x||$, using Lemma 4.2 of Dümbgen et al., (2011), together with (7.1) and (7.2) we can conclude that there exist $\phi_{0}\in\Phi^{d}$ and a subsequence $\phi_{n_{k}}$ such that $C\subseteq\text{dom}(\phi_{0})\subseteq\text{csupp}(Q)$ and,

[TABLE]

Since $\text{dom}(\phi_{n_{k}})\subseteq\text{csupp}(Q)$, we have $\phi_{n_{k}}$ converges to $\phi_{0}$ almost everywhere as the Lebesgue measure of the boundary of $\text{csupp}(Q)$ is zero, then we can conclude $\int e^{\phi_{0}(x)}dx=1$ by dominated convergence. Thus, $\phi_{0}\in\Phi(Q)$. Next, we apply Fatou’s Lemma to the nonnegative functions $x\mapsto\int\text{log}(m(x)+1)Q(dx)+a-b||x||-\text{log}((1-p_{n_{k}})f_{0}+p_{n_{k}}e^{\phi_{n_{k}}})$, and we get,

[TABLE]

Hence,

[TABLE]

which shows $(p_{0},\phi_{0})$ is the maximizer that we are looking for.

∎

Proof of Theorem 3.2.

Since $\stackrel{{\scriptstyle[}}{{n}}\rightarrow\infty]{}{\text{lim}}D_{k}(Q_{n},Q)\rightarrow 0$, hence

[TABLE]

Suppose $\stackrel{{\scriptstyle[}}{{n}}\rightarrow\infty]{}{\text{limsup}}L(Q_{n})=\lambda\in[-\infty,\infty]$, thus there exist a subsequence $\{Q_{n_{k}}\}$, such that $L(Q_{n_{k}})\rightarrow\lambda$. If we let $h(x)=-b_{0}||x||-\text{log}(\int e^{-b_{0}||x||}dx)$ as we did in the proof of Theorem 3.1, then, $h\in\tilde{\Phi}^{d}$, and for any $p>0$,

[TABLE]

Note that in the above inequalities, we used the fact that $\stackrel{{\scriptstyle[}}{{n}}\rightarrow\infty]{}{\text{lim}}\int||x||Q_{n}(dx)=\int||x||Q(dx)$ by Lemma 4.6 of Schuhmacher et al., (2011).

Let $M_{n}=\text{max}_{x\in R^{d}}\phi_{n}(x)$. Since $\stackrel{{\scriptstyle[}}{{n}}\rightarrow\infty]{}{\text{lim}}\int\text{log}(m(x)+1)Q_{n}(dx)=\int\text{log}(m(x)+1)Q(dx)$ by Lemma 4.6 of Schuhmacher et al., (2011), similar to the proof of Theorem 3.1, one can show that for $n$ sufficiently large, we have $L(p_{n},\phi_{n},Q_{n})\rightarrow-\infty$, if $M_{n}\rightarrow\infty$ as $n\rightarrow\infty$, and $L(p_{n},\phi_{n},Q_{n})\leq\int\text{log}(m(x)+1)Q(dx)+M_{n}$, provided that

[TABLE]

Hence there exist some suitable constants $M_{0}$ and $M_{*}$, such that $M_{0}<M_{n_{k}}<M_{*}$ for $k$ sufficiently large and thus $\lambda<\infty$.

Here we explain how (7.3) is derived. As in the proof of Lemma 2.1 of Schuhmacher, Schuhmacher et al., (2011), there exist a simplex $\tilde{\Delta}=\text{conv}(\tilde{x}_{0},\cdots,\tilde{x}_{d})$ with positive Lebesgue measure and open sets $U_{0}$, $U_{1}$, $\cdots$, $U_{d}$ with $Q(U_{j})\geq\eta>0$, for $0\leq j\leq d$, here $\eta=\stackrel{{\scriptstyle[}}{{0}}\leq j\leq d]{}{\text{min}}Q(U_{j})>0$. For any convex and closed set $C$ with $C\cap U_{j}\neq\emptyset$ for all $j$, we have $\tilde{\Delta}\subseteq C$. By Theorem 4.4.4 of Chung, (2001), $\stackrel{{\scriptstyle[}}{{n}}\rightarrow\infty]{}{\text{liminf}}Q_{n}(U_{j})\geq Q(U_{j})\geq\eta$ for all $j$. Thus if $\stackrel{{\scriptstyle[}}{{n}}\rightarrow\infty]{}{\text{lim}}\text{Leb}(C_{n})=0$, then for any $n$ sufficiently large, $\text{Leb}(C_{n})<\text{Leb}(\tilde{\Delta}),\Rightarrow\tilde{\Delta}\nsubseteq C_{n},\Rightarrow$ there exist some $j$, such that $C_{n}\cap U_{j}=\emptyset,\Rightarrow Q_{n}(C_{n})\leq 1-Q_{n}(U_{j})\leq 1-\stackrel{{\scriptstyle[}}{{1}}\leq j\leq d]{}{\text{min}}Q_{n}(U_{j})$. Since,

[TABLE]

thus $\stackrel{{\scriptstyle[}}{{1}}\leq j\leq d]{}{\text{min}}Q_{n}(U_{j})\geq\eta+o(1)$, which shows that $Q_{n}(C_{n})\leq 1-\eta+o(1)$, and hence (7.3) is established.

Now that we know $M_{n_{k}}$ is bounded for $k$ sufficiently large, and $L(p_{n_{k}},\phi_{n_{k}},Q_{n_{k}})\rightarrow\lambda\in R$ as $k\rightarrow\infty$, we may assume $\{p_{n_{k}}\}$, is a convergent sequence, say $p_{n_{k}}\rightarrow p_{*}\in[0,1]$, as $k\rightarrow\infty$. For if $\{p_{n_{k}}\}$ is not convergent, since it is bounded, it must have a convergent subsequence $\{p_{n_{k_{l}}}\}$, and the sequence $\{p_{n_{k_{l}}},\phi_{n_{k_{l}}}\}$ would satisfy all those properties above and we can just simply replace the original sequence with this subsequence.

Again, as in the proof of Theorem 3.1, for any $x_{0}\in\text{interior}(\text{csupp}(Q))$, we have,

[TABLE]

As Lemma 2.13 of Schuhmacher et al., (2011) states that $\stackrel{{\scriptstyle[}}{{n}}\rightarrow\infty]{}{\text{limsup}}h(Q_{n_{k}},x)\leq h(Q,x)$ for any $x\in R^{d}$, we have,

[TABLE]

Hence, for k large enough,

[TABLE]

Again, we can deduce from (7.4) and the boundedness of $M_{n_{k}}$ that there exist constants $a$ and $b>0$ such that,

[TABLE]

Similar as before we conclude that there exist $\phi_{*}\in\Phi^{d}$ and a subsequence $\{\phi_{n_{k_{l}}}\}$ such that $\text{interior}(\text{csupp}(Q))\subseteq\text{dom}(\phi_{*})\subseteq\text{csupp}(Q)$ and,

[TABLE]

Then $\int e^{\phi_{*}(x)}dx=1$ by dominated convergence, which implies that $\phi_{*}\in\tilde{\Phi}^{d}$.

By Skorohod’s theorem, there exist a probability space $(\Omega,\mathcal{F},P)$ and random variables $X_{n_{k_{l}}}\sim Q_{n_{k_{l}}}$, $X\sim Q$, such that $\stackrel{{\scriptstyle[}}{{l}}\rightarrow\infty]{}{\text{lim}}X_{n}=X$ almost surely. Let $H_{n_{k_{l}}}=\int\text{log}(m(x)+1)Q(dx)+a-b||X_{n_{k_{l}}}||-\text{log}\{(1-p_{n_{k_{l}}})f_{0}(X_{n_{k_{l}}})+p_{n_{k_{l}}}\text{exp}(\phi_{n_{k_{l}}}(X_{n_{k_{l}}}))\}$. By Fatou’s Lemma, we have,

[TABLE]

In order to show that $\lambda\geq L(Q)$, we use the approximations $\phi^{*}\leq{\phi^{*}}^{(\epsilon)}\leq{\phi^{*}}^{(1)}$, $0<\epsilon\leq 1$ from Lemma 4.4 of Schuhmacher et al., (2011), since ${\phi^{*}}^{(\epsilon)}\in\Phi^{d}$ is Lipschitz continuous, one can show that $\frac{|{\phi^{*}}^{(\epsilon)}|}{1+||x||}$ is bounded, and hence by Lemma 4.6 of Schuhmacher et al., (2011), we have,

[TABLE]

The last step above is by applying dominated convergence on $e^{{\phi^{*}}^{(\epsilon)}}$ and monotone convergence on $(1-p^{*})f_{0}\int e^{{\phi^{*}}^{(1)}(x)}dx+p^{*}e^{{\phi^{*}}^{(1)}}-(1-p^{*})f_{0}\int e^{{\phi^{*}}^{(\epsilon)}(x)}dx-p^{*}e^{{\phi^{*}}^{(\epsilon)}}$. Thus we have shown that $\lambda=L(Q)$, and $(p^{*},\phi^{*})=(p_{*},\phi_{*})$ is the unique maximizer.

With exactly the same argument, we can show that $\stackrel{{\scriptstyle[}}{{n}}\rightarrow\infty]{}{\text{liminf}}L(Q_{n})=L(Q)$ as well, and hence $L(Q_{n})\rightarrow L(Q)$, as $n\rightarrow\infty$.

Also, if we let $f^{*}=\text{exp}\circ\phi^{*}$, $f_{n}=\text{exp}\circ\phi_{n}$, we have shown that,

[TABLE]

In particular, $\{f_{n_{k_{l}}}\}$ converges to $f^{*}$ almost everywhere w.r.t. Lebesgue measure and hence $\int|f_{n_{k_{l}}}(x)-f^{*}(x)|dx\rightarrow 0$, as $l\rightarrow\infty$ by dominated convergence. Our proof actually shows that for any subsequence of $\{Q_{n}\}$, we can further find a subsequence with the above convergence properties. That means the original sequence must satisfy those properties as well, otherwise we would arrive at contradictions and that completes the proof.

∎

7.2 More Simulation Result

Acknowledgment

Yao’s research is supported by NSF grant DMS-1461677 and Department of Energy with the award No: 10006272.