Distributed sequential method for analyzing massive data

TL;DR

This paper introduces a distributed sequential estimation method for analyzing massive datasets, combining parallel processing, adaptive sampling, and variable selection to improve efficiency and accuracy.

Contribution

It proposes a novel parallel divide-and-conquer approach with adaptive sampling and shrinkage estimation, maintaining statistical properties for large-scale data analysis.

Findings

Method achieves efficient estimation on large datasets.

Adaptive sampling accelerates the estimation process.

Successfully applied to real-world energy and pollution data.

Abstract

To analyse a very large data set containing lengthy variables, we adopt a sequential estimation idea and propose a parallel divide-and-conquer method. We conduct several conventional sequential estimation procedures separately, and properly integrate their results while maintaining the desired statistical properties. Additionally, using a criterion from the statistical experiment design, we adopt an adaptive sample selection, together with an adaptive shrinkage estimation method, to simultaneously accelerate the estimation procedure and identify the effective variables. We confirm the cogency of our methods through theoretical justifications and numerical results derived from synthesized data sets. We then apply the proposed method to three real data sets, including those pertaining to appliance energy use and particulate matter concentration.

Tables

Table 1. Table 1: Simulation results with β 0 = ( − 1 , 1 ) subscript 𝛽 0 1 1 \beta_{0}=(-1,1) and (-1.0,1.0,0.7,0.5,0.2), denoted by S 1 subscript 𝑆 1 S_{1} and S 2 subscript 𝑆 2 S_{2} .

$\beta$	$d$	$M$	method	stopping time	coverage probability	computation time
$S_1$	0.5	1	SM	63.494(16.188)∗	0.95	0.086(0.026)
		2	PSM	65.108(16.325)	(0.956,0.95)+	0.051(0.016)
		5	PSM	71.174(13.888)	(0.952,0.922)	0.027(0.01)
	0.4	1	SM	99.23(19.957)	0.934	0.14(0.033)
		2	PSM	99.466(22.475)	(0.932,0.924)	0.08(0.021)
		5	PSM	104.104(18.914)	(0.948,0.94)	0.04(0.012)
	0.3	1	SM	173.472(25.825)	0.942	0.255(0.042)
		2	PSM	176.484(28.582)	(0.948,0.948)	0.141(0.028)
		5	PSM	179.474(27.596)	(0.93,0.924)	0.068(0.015)
	0.2	1	SM	392.276(39.207)	0.946	0.613(0.071)
		2	PSM	393.452(41.269)	(0.944,0.942)	0.318(0.043)
		5	PSM	400.028(39.049)	(0.956,0.954)	0.145(0.023)
$S_2$	0.5	1	SM	76.964(15.167)	0.914	0.123(0.03)
		2	PSM	88.196(15.5)	(0.934,0.91)	0.073(0.019)
		5	PSM	113.932(14.514)	(0.938,0.882)	0.037(0.01)
	0.4	1	SM	114.196(18.78)	0.934	0.192(0.039)
		2	PSM	125.012(18.163)	(0.944,0.94)	0.11(0.022)
		5	PSM	154.202(18.089)	(0.93,0.9)	0.055(0.013)
	0.3	1	SM	192.826(25.864)	0.95	0.344(0.056)
		2	PSM	205.438(22.773)	(0.942,0.93)	0.19(0.03)
		5	PSM	236.398(23.574)	(0.96,0.926)	0.092(0.017)
	0.2	1	SM	423.92(36.282)	0.944	0.822(0.091)
		2	PSM	430.794(36.571)	(0.95,0.948)	0.42(0.052)
		5	PSM	467.084(37.767)	(0.936,0.924)	0.19(0.03)
∗ Standard deviations are in parentheses.
+ Coverage probabilities in parentheses are for the regions (4) and (8).

Table 2. Table 2: Simulation results of the procedures with the D-optimality criterion, to locate observations with two sets of parameters: β 0 = ( − 1 , 1 ) subscript 𝛽 0 1 1 \beta_{0}=(-1,1) and ( − 1.0 , 1.0 , 0.7 , 0.5 , 0.2 ) 1.0 1.0 0.7 0.5 0.2 (-1.0,1.0,0.7,0.5,0.2) , denoted as S 1 subscript 𝑆 1 S_{1} and S 2 subscript 𝑆 2 S_{2} . The PSM results here are for M = 2 𝑀 2 M=2 .

parameter	$d$	method	stopping time	cover probability	computation time
$S_1$	0.5	SM	26.108(8.119)∗	0.912	4.194(1.65)
		PSM	27.988(6.71)	(0.932,0.928)+	1.016(0.373)
	0.4	SM	42.702(9.816)	0.928	7.522(1.985)
		PSM	42.084(9.597)	(0.936,0.936)	1.686(0.467)
	0.3	SM	75.554(14.305)	0.944	14.194(2.974)
		PSM	75.15(14.548)	(0.938,0.938)	3.236(0.685)
	0.2	SM	176.196(20.177)	0.944	35.207(4.456)
		PSM	176.446(19.616)	(0.958,0.958)	7.799(0.994)
$S_2$	0.5	SM	47.246(11.993)	0.902	7.441(2.575)
		PSM	48.86(9.748)	(0.912,0.896)	1.58(0.539)
	0.4	SM	72.83(12.541)	0.932	12.781(2.867)
		PSM	73.794(13.367)	(0.922,0.918)	2.784(0.677)
	0.3	SM	132.39(16.786)	0.938	26.128(4.532)
		PSM	131.668(17.604)	(0.924,0.924)	5.633(1.02)
	0.2	SM	299.536(25.686)	0.936	69.088(21.81)
		PSM	301.656(26.456)	(0.948,0.948)	14.36(1.959)
∗ Standard deviations are in parentheses.
+ Coverage probabilities in parentheses are for the regions (4) and (8).

Table 3. Table 3: Simulation results of the proposed procedure (PSM) and the common divide-and-conquer method (DC) with M = 2 𝑀 2 M=2 , to locate observations with two sets of parameters: β 0 = ( − 1 , 1 ) subscript 𝛽 0 1 1 \beta_{0}=(-1,1) and ( − 1.0 , 1.0 , 0.7 , 0.5 , 0.2 ) 1.0 1.0 0.7 0.5 0.2 (-1.0,1.0,0.7,0.5,0.2) , denoted as S 1 subscript 𝑆 1 S_{1} and S 2 subscript 𝑆 2 S_{2} .

			PSM
parameter	$\rho$		d=0.5	0.4	0.3	0.2	DC
S1	0.01	N	27.42(4.345)	38.176(8.119)	65.694(12.779)	150.636(17.864)	6000(0)
		SE	0.049(0.059)	0.031(0.036)	0.017(0.022)	0.007(0.009)	0.002(0.001)
		AD	0.237(0.139)	0.189(0.107)	0.136(0.081)	0.089(0.051)	0.055(0.019)
	0.05	N	27.636(4.459)	38.324(7.752)	66.352(11.996)	149.584(18.151)	6000(0)
		SE	0.053(0.06)	0.035(0.046)	0.018(0.026)	0.008(0.01)	0.043(0.007)
		AD	0.251(0.143)	0.194(0.121)	0.136(0.081)	0.093(0.052)	0.264(0.024)
	0.10	N	27.874(4.495)	38.378(7.892)	65.98(12.529)	150.598(19.862)	6000(0)
		SE	0.053(0.057)	0.037(0.046)	0.018(0.024)	0.008(0.01)	0.171(0.016)
		AD	0.253(0.139)	0.201(0.12)	0.137(0.083)	0.092(0.055)	0.53(0.029)
	0.15	N	27.872(4.512)	38.554(7.921)	66.194(12.023)	151.664(21.38)	6000(0)
		SE	0.052(0.059)	0.033(0.041)	0.018(0.025)	0.008(0.01)	0.387(0.026)
		AD	0.243(0.139)	0.188(0.114)	0.138(0.087)	0.09(0.055)	0.801(0.032)
S2	0.01	N	49.374(9.544)	74.348(13.115)	133.454(17.453)	301.856(26.503)	6000(0)
		SE	0.061(0.049)	0.036(0.028)	0.018(0.014)	0.009(0.007)	0.003(0.002)
		AD	0.424(0.156)	0.324(0.121)	0.233(0.084)	0.162(0.059)	0.093(0.032)
	0.05	N	48.326(9.108)	73.042(13.383)	130.872(18.231)	300.338(25.382)	6000(0)
		SE	0.054(0.049)	0.033(0.028)	0.018(0.014)	0.008(0.006)	0.046(0.01)
		AD	0.397(0.158)	0.312(0.122)	0.231(0.087)	0.158(0.058)	0.433(0.048)
	0.10	N	49.53(9.74)	74.774(13.441)	132.256(17.763)	302.33(25.535)	6000(0)
		SE	0.058(0.051)	0.033(0.027)	0.018(0.014)	0.008(0.005)	0.189(0.025)
		AD	0.411(0.169)	0.309(0.12)	0.229(0.08)	0.153(0.051)	0.888(0.059)
	0.15	N	49.112(9.981)	73.792(13.892)	132.136(18.651)	302.398(25.299)	6000(0)
		SE	0.058(0.046)	0.035(0.03)	0.019(0.014)	0.008(0.006)	0.434(0.047)
		AD	0.416(0.163)	0.321(0.127)	0.233(0.084)	0.159(0.057)	1.358(0.072)
∗ Standard deviations are in parentheses.

Table 4. Table 4: Simulation results of the common divide-and-conquer method (DC) with M = 2 , 5 , 10 , 15 , 20 𝑀 2 5 10 15 20 M=2,5,10,15,20 and total sample size N = 500 , 1000 , 2000 𝑁 500 1000 2000 N=500,1000,2000 .

N		M=2	5	10	15	20
500	SE	0.043(0.013)	0.05(0.014)	0.068(0.014)	0.11(0.016)	0.248(0.025)
	AD	0.739(0.125)	0.797(0.114)	0.93(0.1)	1.182(0.088)	1.767(0.089)
1000	SE	0.02(0.006)	0.022(0.006)	0.025(0.006)	0.029(0.007)	0.034(0.007)
	AD	0.511(0.083)	0.529(0.08)	0.563(0.076)	0.605(0.07)	0.657(0.066)
2000	SE	0.01(0.003)	0.01(0.003)	0.011(0.003)	0.012(0.003)	0.012(0.003)
	AD	0.358(0.059)	0.364(0.058)	0.374(0.056)	0.386(0.054)	0.399(0.053)
∗ Standard deviations are in parentheses.

Table 5. Table 5: Adaptive shrinkage estimation results with β 0 = ( − 2 , 1 , 1.5 , 2 , 0 , 0 , 0 , 0 , 0 , 0 ) subscript 𝛽 0 2 1 1.5 2 0 0 0 0 0 0 \beta_{0}=(-2,1,1.5,2,0,0,0,0,0,0) and ( − 2 , 2 , 2 , 2 , 0 , … , 0 ) 2 2 2 2 0 … 0 (-2,2,2,2,0,...,0) , denoted by S 1 subscript 𝑆 1 S_{1} and S 2 subscript 𝑆 2 S_{2} , where there are 46 zero cases in S 2 subscript 𝑆 2 S_{2} .

$\beta$	$d$	$M$	method	stopping time	coverage probability	computation time	${\hat{p}}_0$
$S_1$	0.5	1	SM	111.772(18.947)	0.898	0.155(0.035)	4.108(0.323)
		5	PSM	165.812(16.988)	0.976	0.038(0.011)	3.958(0.254)
	0.4	1	SM	168.122(25.127)	0.92	0.267(0.053)	4.076(0.273)
		5	PSM	220.206(22.532)	0.952	0.06(0.013)	3.98(0.178)
	0.3	1	SM	289.008(34.338)	0.91	0.577(0.105)	4.046(0.21)
		5	PSM	338.978(31.237)	0.964	0.106(0.017)	3.996(0.089)
	0.2	1	SM	635.088(53.877)	0.934	2.414(0.484)	4.02(0.14)
		5	PSM	684.934(51.767)	0.958	0.255(0.03)	4(0)
$S_2$	0.5	1	SM	155.836(21.272)	0.844	0.837(0.204)	4.25(0.587)
		5	PSM	358.108(17.462)	0.926	0.178(0.056)	3.992(0.089)
	0.4	1	SM	215.852(28.16)	0.816	1.45(0.307)	4.294(0.583)
		5	PSM	413.196(27.578)	0.932	0.298(0.08)	4(0)
	0.3	1	SM	340.644(35.802)	0.848	2.989(0.507)	4.194(0.444)
		5	PSM	540.934(36.577)	0.96	0.573(0.11)	4(0)
	0.2	1	SM	696.362(65.537)	0.86	10.419(2)	4.144(0.379)
		5	PSM	912.79(52.673)	0.956	1.398(0.198)	4(0)
∗ Standard deviations are in parentheses.
+ Coverage probabilities in parentheses are for the regions (4) and (8).

Table 6. Table 6: Stopping times and computation times for data sets: Energy and PM2.5 in Beijing and Shanghai. The upper and lower panels involve randomly selected samples and feature a D-optimal design.

				Random selection
			Stopping time					Computation time
Dataset	$M$	Method	$d=0.5$	0.4	0.3	0.2		$d=0.5$	0.4	0.3	0.2
Energy	1	SM	2053	3200	5370	11311		6	13	37	150
	5	PSM	1876	2805	4967	10829		1	2	3	11
Beijing	1	SM	4023	6182	11329	25193		15	33	116	557
	5	PSM	4088	6233	10981	25005		2	3	7	30
Shanghai	1	SM	7559	11989	23182	31880∗		59	146	475	856
	5	PSM	8673	13888	23835	31880∗		6	13	37	69
				D-optimal design
			Stopping time					Computation time
Dataset	$M$	Method	$d=0.5$	0.4	0.3	0.2		$d=0.5$	0.4	0.3	0.2
Energy	1	SM	287	522	1200	5461		393	1134	4772	67439
	5	PSM	371	601	1153	5597		11	31	87	932
Beijing	1	SM	576	1072	2656	9608		3219	7810	32572	320896
	5	PSM	594	1048	2572	9589		51	115	396	3222
Shanghai	1	SM	2122	3935	12480	31880∗		13913	39275	283890	813693
	5	PSM	2151	4090	14192	31880∗		210	583	5585	13793
∗ The stopping criterion is not satisfied, even when all samples are used.

Table 7. Table 7: Parameter estimation for the Energy data set.

	$M=1$ (SM)					$M=5$ (PSM)
Para.	d=0.5	0.4	0.3	0.2		d=0.5	0.4	0.3	0.2
T1	0.381(0.05)∗	0.336(0.04)	0.337(0.03)	0.304(0.02)		0.299(0.051)	0.367(0.04)	0.36(0.03)	0.316(0.02)
RH1	0.625(0.052)	0.578(0.041)	0.568(0.031)	0.551(0.021)		0.502(0.053)	0.562(0.042)	0.57(0.032)	0.559(0.021)
T2	-0.416(0.068)	-0.324(0.055)	-0.343(0.041)	-0.295(0.027)		-0.291(0.068)	-0.385(0.054)	-0.371(0.041)	-0.314(0.027)
RH2	-0.491(0.059)	-0.412(0.047)	-0.398(0.036)	-0.364(0.024)		-0.336(0.06)	-0.408(0.048)	-0.419(0.036)	-0.381(0.024)
T3	0.226(0.049)	0.18(0.038)	0.205(0.03)	0.167(0.02)		0.221(0.051)	0.241(0.04)	0.188(0.03)	0.159(0.02)
RH3	0.123(0.037)	0.084(0.029)	0.086(0.022)	0.046(0.015)		0.035(0.038)	0.067(0.03)	0.05(0.022)	0.043(0.015)
T4	-0.063(0.033)	-0.082(0.027)	-0.064(0.02)	-0.066(0.014)		-0.078(0.035)	-0.08(0.027)	-0.077(0.02)	-0.071(0.014)
RH4	-0.283(0.03)	-0.296(0.024)	-0.286(0.018)	-0.254(0.012)		-0.223(0.031)	-0.231(0.025)	-0.234(0.019)	-0.258(0.013)
Prec	-0.039(0.015)	-0.043(0.012)	-0.036(0.009)	-0.034(0.006)		-0.033(0.015)	-0.023(0.012)	-0.028(0.009)	-0.035(0.006)
WS	0.014(0.015)	0.016(0.012)	0.024(0.01)	0.022(0.006)		0.014(0.016)	0.025(0.013)	0.032(0.01)	0.026(0.007)
V	0.01(0.014)	0.013(0.011)	0.004(0.009)	0.002(0.006)		0.016(0.015)	-0.002(0.012)	-0.007(0.009)	0.002(0.006)
DEWP	-0.008(0.04)	-0.024(0.032)	-0.056(0.025)	-0.069(0.017)		-0.101(0.043)	-0.101(0.034)	-0.064(0.025)	-0.047(0.017)
∗ Standard variance estimation in parentheses.

Table 8. Table 8: Parameter estimation for the Beijing PM2.5 and Shanghai PM2.5 data sets.

			Beijing PM2.5
$M$∗	Method	$d$	DEWP	HUMI	PRES	TEMP	Iws	Prec	Iprec
1	SM	0.5	0.399(0.091)+	0.176(0.053)	-0.278(0.025)	-0.589(0.075)	-0.248(0.015)	-0.002(0.012)	-0.104(0.018)
		0.4	0.393(0.073)	0.196(0.042)	-0.247(0.021)	-0.577(0.06)	-0.236(0.012)	-0.001(0.009)	-0.113(0.016)
		0.3	0.291(0.055)	0.27(0.032)	-0.226(0.015)	-0.485(0.045)	-0.232(0.009)	-0.012(0.008)	-0.093(0.01)
		0.2	0.291(0.036)	0.27(0.021)	-0.211(0.01)	-0.476(0.03)	-0.234(0.006)	-0.025(0.006)	-0.085(0.006)
5	PSM	0.5	0.282(0.091)	0.29(0.052)	-0.243(0.026)	-0.484(0.075)	-0.224(0.016)	-0.092(0.044)	-0.24(0.042)
		0.4	0.283(0.073)	0.286(0.042)	-0.246(0.021)	-0.482(0.06)	-0.24(0.012)	-0.054(0.023)	-0.21(0.028)
		0.3	0.3(0.055)	0.258(0.031)	-0.251(0.015)	-0.504(0.045)	-0.249(0.009)	-0.074(0.017)	-0.135(0.016)
		0.2	0.35(0.036)	0.231(0.021)	-0.223(0.01)	-0.524(0.03)	-0.235(0.006)	-0.026(0.009)	-0.12(0.008)
			Shanghai PM2.5
$M$	Method	$d$	DEWP	HUMI	PRES	TEMP	Iws	Prec	Iprec
1	SM	0.5	-0.198(0.088)	-0.051(0.042)	-0.163(0.017)	-0.255(0.083)	-0.207(0.008)	-0.013(0.01)	-0.101(0.009)
		0.4	-0.135(0.07)	-0.086(0.034)	-0.176(0.014)	-0.325(0.066)	-0.21(0.006)	-0.016(0.008)	-0.096(0.008)
		0.3	-0.148(0.053)	-0.08(0.025)	-0.166(0.01)	-0.312(0.05)	-0.216(0.005)	-0.012(0.005)	-0.088(0.005)
		0.2	-0.11(0.045)	-0.101(0.021)	-0.164(0.008)	-0.345(0.042)	-0.216(0.004)	-0.017(0.004)	-0.084(0.004)
5	PSM	0.5	-0.209(0.088)	-0.047(0.042)	-0.145(0.016)	-0.235(0.083)	-0.213(0.007)	-0.028(0.011)	-0.101(0.01)
		0.4	-0.174(0.071)	-0.072(0.033)	-0.17(0.013)	-0.294(0.066)	-0.213(0.006)	-0.003(0.006)	-0.111(0.008)
		0.3	-0.128(0.053)	-0.095(0.025)	-0.173(0.01)	-0.34(0.05)	-0.216(0.005)	-0.01(0.005)	-0.093(0.005)
		0.2	-0.116(0.045)	-0.098(0.021)	-0.164(0.008)	-0.34(0.042)	-0.216(0.004)	-0.016(0.004)	-0.087(0.005)
∗ Num. stands for number of machines in parallel sequential estimation.
+ Standard variance estimation in parentheses.

Table 9. Table 9: Number of effective variables, and stopping and computation times, with ASE for the data sets: Energy and PM2.5 in Beijing and Shanghai.

			Stopping time					Computation time					${\hat{p}}_0$
Dataset	$M$	Method	d=0.5	0.4	0.3	0.2		d=0.5	0.4	0.3	0.2		d=0.5	0.4	0.3	0.2
Energy	1	SM	1680	2563	4743	10627		47	131	765	8583		9	10	11	12
	5	PSM	1439	2196	4001	8977		1	3	13	91		6	6	8	8
Beijing	1	SM	3666	5735	10305	25193		351	1241	6722	92463		7	7	7	8
	5	PSM	2931	4147	10539	24013		9	25	105	907		6	6	7	7
Shanghai	1	SM	7040	10703	20900	31880		2249	7149	51720	187753		7	7	7	8
	5	PSM	4878	8171	15355	28660		20	102	430	1427		4	4	4	6

Table 10. Table 10: Parameter estimation for the Energy data set.

	$M=1$ (SM)					$M=5$ (PSM)
Para.	d=0.5	0.4	0.3	0.2		d=0.5	0.4	0.3	0.2
T1	0.376(0.055)	0.378(0.045)	0.348(0.032)	0.305(0.021)		0(0)	0.326(0.045)	0.37(0.033)	0.322(0.023)
RH1	0.615(0.058)	0.608(0.046)	0.573(0.033)	0.546(0.022)		0.518(0.062)	0.541(0.048)	0.575(0.035)	0.565(0.023)
T2	-0.411(0.076)	-0.376(0.061)	-0.353(0.044)	-0.291(0.028)		0(0)	-0.346(0.06)	-0.386(0.045)	-0.321(0.03)
RH2	-0.48(0.066)	-0.444(0.052)	-0.398(0.038)	-0.359(0.025)		-0.353(0.069)	0(0)	-0.414(0.039)	-0.394(0.027)
T3	0.227(0.055)	0.206(0.043)	0.214(0.031)	0.166(0.021)		0.248(0.056)	0.225(0.045)	0.213(0.032)	0.151(0.022)
RH3	0.12(0.041)	0.096(0.033)	0.09(0.024)	0.046(0.016)		0(0)	0(0)	0(0)	0(0)
T4	-0.062(0.036)	-0.085(0.029)	-0.064(0.021)	-0.073(0.014)		0(0)	0(0)	-0.083(0.022)	-0.074(0.015)
RH4	-0.261(0.033)	-0.288(0.027)	-0.284(0.02)	-0.259(0.013)		-0.225(0.034)	-0.234(0.028)	-0.23(0.02)	-0.253(0.014)
Prec	0(0)	-0.043(0.013)	-0.033(0.01)	-0.035(0.006)		0(0)	0(0)	0(0)	0(0)
WS	0(0)	0(0)	0(0)	0.023(0.007)		0(0)	0(0)	0(0)	0(0)
V	0(0)	0(0)	0(0)	0(0)		0(0)	0(0)	0(0)	0(0)
DEWP	0(0)	0(0)	-0.059(0.026)	-0.065(0.017)		-0.098(0.049)	0(0)	0(0)	0(0)
∗ Standard variance estimation in parentheses.

Table 11. Table 11: Parameter estimation for the Beijing PM2.5 and Shanghai PM2.5 data sets.

			Beijing PM2.5
$M$	Method	$d$	DEWP	HUMI	PRES	TEMP	Iws	Prec	Iprec
	SM	0.5	0.397(0.096)	0.174(0.055)	-0.275(0.027)	-0.593(0.079)	-0.248(0.016)	0(0)	-0.099(0.018)
		0.4	0.401(0.077)	0.193(0.044)	-0.246(0.021)	-0.577(0.063)	-0.24(0.013)	0(0)	-0.112(0.016)
		0.3	0.317(0.057)	0.251(0.033)	-0.224(0.016)	-0.507(0.047)	-0.232(0.009)	0(0)	-0.091(0.01)
		0.2	0.291(0.036)	0.27(0.021)	-0.211(0.01)	-0.476(0.03)	-0.234(0.006)	-0.025(0.006)	-0.085(0.006)
5	PSM	0.5	0(0)	0.354(0.062)	-0.253(0.03)	-0.41(0.09)	-0.221(0.018)	0(0)	-0.271(0.036)
		0.4	0(0)	0.29(0.051)	-0.255(0.025)	-0.475(0.074)	-0.235(0.015)	0(0)	-0.237(0.032)
		0.3	0.288(0.056)	0.263(0.032)	-0.25(0.016)	-0.492(0.046)	-0.25(0.009)	0(0)	-0.129(0.011)
		0.2	0.346(0.037)	0.233(0.021)	-0.224(0.01)	-0.52(0.031)	-0.237(0.006)	0(0)	-0.124(0.008)
			Shanghai PM2.5
$M$	Method	$d$	DEWP	HUMI	PRES	TEMP	Iws	Prec	Iprec
1	SM	0.5	-0.222(0.092)	-0.042(0.044)	-0.165(0.017)	-0.237(0.086)	-0.208(0.008)	0(0)	-0.1(0.01)
		0.4	-0.164(0.074)	-0.073(0.035)	-0.178(0.014)	-0.299(0.069)	-0.21(0.006)	0(0)	-0.107(0.009)
		0.3	-0.156(0.055)	-0.079(0.026)	-0.169(0.01)	-0.306(0.052)	-0.212(0.005)	0(0)	-0.093(0.006)
		0.2	-0.11(0.045)	-0.101(0.021)	-0.164(0.008)	-0.345(0.042)	-0.216(0.004)	-0.017(0.004)	-0.084(0.004)
5	PSM	0.5	0(0)	0(0)	-0.133(0.021)	0(0)	-0.221(0.01)	0(0)	-0.12(0.011)
		0.4	0(0)	0(0)	-0.151(0.016)	0(0)	-0.212(0.008)	0(0)	-0.104(0.009)
		0.3	0(0)	0(0)	-0.166(0.012)	0(0)	-0.214(0.006)	0(0)	-0.098(0.006)
		0.2	0(0)	-0.092(0.023)	-0.17(0.009)	-0.333(0.045)	-0.215(0.004)	0(0)	-0.089(0.005)
∗ Standard variance estimation in parentheses.

Equations

$$Y={\mathbf{X}}^{\top }\mathbf{\beta }+ϵ,$$

$${\hat{\sigma }}_{jn}^2=\frac{1}{n-p}\sum _{i=1}^n(y_{ji}-{\mathbf{x}}_{ji}^{\top }{\hat{\mathbf{\beta }}}_{jn}{)}^2$$

$$\sum _{j=1}^M{\tilde{a}}_j^2=a^2,$$

$$N_j=N_{jd}\equiv \inf \left\{n\ge n_0:({\hat{\sigma }}_{jn}^2+\frac{1}{n})\le \frac{d^{2}n}{{\tilde{a}}_j^2{\mu }_{jn}}\right\}.$$

$$R_{N_j}=\left\{\mathbf{z}\in R^p:\frac{S_{N_j}}{N_j}\le \frac{d^2}{{\mu }_{j{N}_j}}\right\}.$$

$$\hat{\mathbf{\beta }}\equiv {\hat{\mathbf{\beta }}}_{N^{\ast }}=\sum _{j=1}^M{\rho }_j{\hat{\mathbf{\beta }}}_{j{N}_j}$$

$$R_{N^{\ast }}=\left\{\mathbf{z}\in R^p:(\mathbf{z}-\hat{\mathbf{\beta }}{)}^{\top }{\left\{\sum _{j=1}^M{\rho }_j^2({\mathbf{X}}_{j{N}_j}{\mathbf{X}}_{j{N}_j}^{\top }{)}^{-1}\right\}}^{-1}(\mathbf{z}-\hat{\mathbf{\beta }})\le \frac{N^{\ast }{d}^2}{{\mu }_{N^{\ast }}}\right\},$$

$$N\equiv \inf \left\{n\ge n_0:\text{mbox}and({\hat{\sigma }}_n^2+\frac{1}{n})\le \frac{d^{2}n}{a^2{\mu }_n}\right\}.$$

$$\underset{n}{\sup }E(|{ϵ}_n{|}^{\alpha }|{\mathcal{F}}_{n-1})<\infty \text{mbox}almostsurelyforsome\alpha >2,$$

$${\lambda }_{min}(\sum _{i=1}^n{x}_i{x}_i^T)\to \infty \text{ and }\log ({\lambda }_{max}(\sum _{i=1}^n{x}_i{x}_i^T))=o({\lambda }_{min}(\sum _{i=1}^n{x}_i{x}_i^T)).$$

$$(\mathbf{A}+\mathbf{U}{\mathbf{U}}^T{)}^{-1}={\mathbf{A}}^{-1}-\frac{({\mathbf{A}}^{-1}\mathbf{U})({\mathbf{U}}^T{\mathbf{A}}^{-1})}{1+{\mathbf{U}}^T{\mathbf{A}}^{-1}\mathbf{U}}.$$

$$det(\mathbf{A}+\mathbf{U}{\mathbf{U}}^T)=det(\mathbf{A})(1+{\mathbf{U}}^T{\mathbf{A}}^{-1}\mathbf{U}).$$

$${\left(\sum _{j=1}^M\frac{{\rho }_j^2}{N_j}\right)}^{-1}\left\{\sum _{j=1}^M{\rho }_j^2({\mathbf{X}}_{j{N}_j}{\mathbf{X}}_{j{N}_j}^{\top }{)}^{-1}\right\}-N^{\ast }{\left(\sum _{j=1}^M{\mathbf{X}}_{j{N}_j}{\mathbf{X}}_{j{N}_j}^{\top }\right)}^{-1}⟶0,$$

$$R_{N^{\ast }}^{{}^{\prime }}=\left\{Z\in R^p:(\hat{\mathbf{\beta }}-\mathbf{z}{)}^{\top }\left(\sum _{j=1}^M{\mathbf{X}}_{j{N}_j}{\mathbf{X}}_{j{N}_j}^{\top }\right)(\hat{\mathbf{\beta }}-\mathbf{z})\le \frac{N^{\ast }{d}^2}{{\mu }_{N^{\ast }}}\right\}.$$

$$\left(\sum _{j=1}^M{X}_{j{N}_j}{X}_{j{N}_j}^T\right)\left\{\sum _{j=1}^M{\rho }_j^2(X_{j{N}_j}{X}_{j{N}_j}^T{)}^{-1}\right\}$$

$$\ge \left(\sum _{j=1}^M{\rho }_j^2\right)I+\sum _{j\ne l}{\rho }_l{\rho }_{j}I=I,$$

$${\lambda }_{max}\left\{{\left(\sum _{j=1}^M{X}_{j{N}_j}{X}_{j{N}_j}^T\right)}^{-1}\right\}\le {\lambda }_{max}\left[\left\{\sum _{j=1}^M{\rho }_j^2(X_{j{N}_j}{X}_{j{N}_j}^T{)}^{-1}\right\}\right].$$

$$n^{1/2}\lambda \to 0\text{ and }{n}^{1/2+\gamma \delta }\lambda ⟶\infty ,\text{mbox}asn\to \infty .$$

$${\tilde{N}}_j=N_d\equiv \inf \left\{k:k\ge n_0\text{mbox}and({\hat{\sigma }}_{jk}^2+\frac{1}{k})\le \frac{d^{2}k}{{\tilde{a}}_{jk}^2{\tilde{\mu }}_{jk}}\right\},$$

$${\mathbf{I}}^{\ast }=\prod _{j=1}^M{\mathbf{I}}_j(ϵ)$$

$$\tilde{\mathbf{\beta }}={\mathbf{I}}^{\ast }(\sum _{j=1}^M{\rho }_j{\mathbf{\beta }}_j^{\ast }),$$

$$\sqrt{{\tilde{N}}_j}({\mathbf{\beta }}_j^{\ast }-{\mathbf{\beta }}_0)⟶N(0,{\sigma }^2{\mathbf{I}}_0{\Sigma }^{-1}{\mathbf{I}}_0)\text{in distribution as }d\to 0,$$

$$\sqrt{\tilde{N}}(\tilde{\mathbf{\beta }}-{\mathbf{\beta }}_0)⟶N(0,{\sigma }^2{\mathbf{I}}_0{\Sigma }^{-1}{\mathbf{I}}_0)\text{in distribution as }d\to 0.$$

$${\hat{\sigma }}^{-2}(\tilde{\mathbf{\beta }}-{\mathbf{\beta }}_0{)}^{\top }{\left[{\mathbf{I}}^{\ast }{\left\{\sum _{j=1}^M{\mathbf{X}}_{j{N}_j}{\mathbf{X}}_{j{N}_j}^{\top }\right\}}^{-1}{\mathbf{I}}^{\ast }\right]}^{-}(\tilde{\mathbf{\beta }}-{\mathbf{\beta }}_0)⟶{\chi }_{p_0}^2,$$

$${\hat{\sigma }}^{-2}(\tilde{\mathbf{\beta }}-{\mathbf{\beta }}_0{)}^{\top }{\left[{\mathbf{I}}^{\ast }{\left\{\sum _{j=1}^M{\mathbf{X}}_{j{N}_j}{\mathbf{X}}_{j{N}_j}^{\top }\right\}}^{-1}{\mathbf{I}}^{\ast }\right]}^{-}(\tilde{\mathbf{\beta }}-{\mathbf{\beta }}_0)$$

$$={\hat{\sigma }}^{-2}({\tilde{\mathbf{\beta }}}_1-{\mathbf{\beta }}_{01}{)}^{\top }{\tilde{\mathbf{\Sigma }}}_{11}({\tilde{\mathbf{\beta }}}_1-{\mathbf{\beta }}_{01}),$$

$$R_{\tilde{N}}=\left\{\mathbf{Z}\in R^p:\frac{S_{\tilde{N}}}{\tilde{N}}\le \frac{d^2}{{\nu }_{\tilde{N}}}\text{mbox}and{z}_j=0\text{mbox}for{I}_{jj}^{\ast }=0,1\le j\le p\right\},$$

$$\displaystyle\tilde{\bm{{\Sigma}}}=\left(\begin{array}[]{cc}\tilde{\bm{{\Sigma}}}_{11}&0\\ 0&0\\ \end{array}\right).$$

$$Y={\mathbf{X}}^T{\mathbf{\beta }}_0+ϵ,$$

$$\frac{(\mathbf{z}-{\hat{\mathbf{\beta }}}_{j{N}_j}{)}^{\top }({\mathbf{X}}_{j{N}_j}{\mathbf{X}}_{j{N}_j}^{\top })(\mathbf{z}-{\hat{\mathbf{\beta }}}_{j{N}_j})}{N_j}\le \frac{d^2}{{\mu }_{j{N}_j}}.$$

$$D=2{\left(\frac{N_j{d}^2}{{\mu }_{j{N}_j}}\right)}^{1/2}{\lambda }_{max}^{1/2}\left[\left({\mathbf{X}}_{j{N}_j}{\mathbf{X}}_{j{N}_j}^{\top }{)}^{-1}\right\}\right]=2{\left(\frac{d^2}{{\mu }_{j{N}_j}}\right)}^{1/2}{\lambda }_{min}^{1/2}\left[\left({\mathbf{X}}_{j{N}_j}{\mathbf{X}}_{j{N}_j}^{\top })/N_j\right\}\right],$$

Full text

Distributed sequential method for analyzing massive data

Zhanfeng Wang  Yuan-chin Ivan Chang111Corresponding author: E-mail: ycchang@gate.sinica.edu.tw; ivan.chang.1@gmail.com

Abstract: To analyse a very large data set containing lengthy variables, we adopt a sequential estimation idea and propose a parallel divide-and-conquer method. We conduct several conventional sequential estimation procedures separately, and properly integrate their results while maintaining the desired statistical properties. Additionally, using a criterion from the statistical experiment design, we adopt an adaptive sample selection, together with an adaptive shrinkage estimation method, to simultaneously accelerate the estimation procedure and identify the effective variables. We confirm the cogency of our methods through theoretical justifications and numerical results derived from synthesized data sets. We then apply the proposed method to three real data sets, including those pertaining to appliance energy use and particulate matter concentration.

Keywords: Sequential sampling; Stopping rule; Confidence set; Distributed/Parallel computation

AMS Subject Classification (2000): Primary 62F12; Secondary 62E20

1 Introduction

While the development of modern measurement and communication technologies has frequently made data collection procedures more efficient, we as researchers have been hard-pressed to analyse and extract information from large data sets and to keep up with our data collection capacity. Although we can leverage concepts from the divide-and-conquer algorithm to analyse very large data sets—such that we can apply longstanding statistical methods without modifying existing procedures and computation facilities—the number of partitions, the size of each partition, and some tuning parameters will soon become follow-up issues. From a statistical perspective, the more important issue will be how to legitimately and effectively merge the individual result from each partition into an informative one. This is an essential issue when we apply such divide-and-conquer thinking, especially when there is a lack of consistency among the partition results.

From a computation perspective, the distributed/parallel computation method is a powerful way of accelerating the computation procedure, when statistical procedures can be ‘parallelized’ without overly modifying their current algorithms. Leveraging the divide-and-conquer concept, many researchers have applied distributed computation methods to statistical hypothesis testing and estimation; one can reference Chen and Xie (2012); Zhang et al. (2015); Battey et al. (2015); Lu et al. (2016), and the references therein. These studies are proposed under fixed sample size scenarios, and there is currently a dearth of research about how to integrate sequential procedures into a distributed computation method. The current study looks to fill this research gap, at least in part.

When analysing a very large data set featuring lengthy variables, the computation issue becomes an essential one. Computer scientists may want to resolve these issues from algorithm and hardware perspectives, but these usually require a complicated software setup and/or modern computation facilities. In the current study, we adopt sequential estimation methods for regression models and allow the model for each ‘partition’ to sequentially choose from the data pool its own new subjects, until its stopping criterion is fulfilled. With a suitable selected stopping criterion, we are able to combine and integrate the models and maintain good statistical properties in such an integrated estimate. Conventionally, how to partition a large data set into several small ones may affect the overall analytical performance. Here, the data sizes of the individual partitions differ and depend on the corresponding performance of the model of each partition. In particular, we adopt the method of the fixed-size confidence set estimation (Siegmund, 1985), such that when data recruiting is stopped, the estimates will have a prescribed accuracy. Because the accuracy of the coefficient estimates from each partition is under control, we are able to merge these estimates into one and retain the required statistical properties of the original sequential estimation.

The method of sequential analysis was established by Wald (1945); Ward (1947), and it has been applied to many areas since then, including clinical trials, finance, engineering control, and educational and psychological tests, inter alia. The major feature of the sequential method is that it allows the sample size to be random and depend on the observed information (Chow and Robbins, 1965; Woodroofe, 1982; Whitehead, 1997; Bartroff et al., 2013). In addition to the random sample sizes, some sequential methods allow users to recruit new observations, based on information gathered while analysing the current-stage data. This type of sequential method is common in the literature on stochastic control, educational test, active learning in machine learning, and the like (Lai and Wei, 1982; Lord, 1971; Wainer, 2000; Deng et al., 2009). The identification of important variables is a critical feature in applications, especially where the ability to interpret model results is essential and the data set has a lengthy list of variables. To this end, we adopt the adaptive shrinkage estimate (ASE) of Wang and Chang (2013), such that we can effectively detect high-impact variables during the sequential modelling procedure.

To demonstrate our method, we conduct several sequential fixed-sized confidence set estimation procedures at once on different machines (or central processing units (CPUs)); then, once all the ‘sampling’ procedures are stopped, we combine their estimates. This kind of computation scheme is similar to ‘distributed’ computing in the computer science literature; it can also be viewed as a naive parallel procedure. Please note that when we already have very large data sets in hand—as in many modern data analysis scenarios—we simply acquire new observations from an existing data set and ‘re-estimate’ the regression coefficients of the corresponding model; this situation is different from traditional sequential analysis applications, where we need to actually ‘collect’ new samples.

The remainder of this paper is organized as follows. In Section 2, we review the general sequential fixed-sized confidence region estimate, introduce our distributed sequential estimation procedure, discuss how to combine individual estimates into a final estimate when all procedures have finished, and study the method’s asymptotic properties. Furthermore, in Section 3, we present a sequential procedure with an ASE. Section 4 summarizes the numerical results by using the synthesized data and some real examples, including a data set pertaining to appliance energy use, and two fine particle/particulate matter (PM2.5) concentration data sets. Technical details are presented in the Appendix.

2 Distributed sequential estimation in linear models

Consider a linear regression model,

[TABLE]

where $Y\in R$ is a response variable, $\bm{{X}}\in\bm{{R}}^{p}$ is a covariate vector with length $p$, $\bm{{\beta}}$ is an unknown vector of parameters to be estimated, and $\epsilon$ denotes the random error with mean $\hbox{E}(\epsilon)=0$ and variance $\hbox{Var}(\epsilon)=\sigma^{2}>0$. Building a fixed-sized confidence set for $\bm{{\beta}}$ with a prespecified coverage probability is a classical problem, and some authors have extended the thinking of Chow and Robbins (1965) to linear regression models Albert (1966); Gleser (1965); Srivastava (1971). Since then, many papers have appeared in the literature that feature various setups and perspectives Grambsch (1989); Chang and Martinsek (1992); Muthu Poruthotage (2013); Vorobeichikov and Konev (2017).

Assume that there is already a very large data set available for analysis, similar to that seen in most ‘big data’ scenarios. When its size is too large—such that using the data set all at once is impractical—the divide-and-conquer method is an economical, hands-on approach that requires fewer changes to existing tools (e.g. hardware and software). If we simultaneously run $M$ independent estimation procedures (e.g. use $M$ independent machines), then we will definitely reduce the computation time; at this point, how to partition the data set appropriately into $M$ partitions becomes an issue. Most practitioners follow the conventional sampling method by treating the available data set as a pseudo-population. Besides those studies that examine how to divide the original data set without introducing biases, as well as partition sizes, there has been a lack of discussion in the literature regarding how to legitimately integrate all the results from each partition into a statistically meaningful estimate. This will be an important issue when there is inconsistency among the results derived from partitioning. One of the advantages inherent in using a sequential methodology is that we do not need to specify the sample size first. For this reason, we can conduct $M$ sequential fixed-width confidence set estimation procedures on these machines. With the properties of the fixed-sized confidence set estimation procedure, the estimates of various sequential estimation procedures will have a similar estimation accuracy and coverage probabilities; this will allow us to combine them into an estimate while retaining the desired statistical properties.

2.1 Sequential fixed-sized confidence set estimate for regression model

Our present goal is to estimate regression parameters and build with them a confidence set of a prescribed size. We independently run $M$ estimation procedures; each of them will independently and sequentially recruit observations from the same data set, without replacement. When all sequential procedures are stopped, we combine their results into one estimate, such that it can retain the desired statistical properties of the conventional sequential confidence set estimation procedure with less computation time and only some of the data. We run each procedure independently, and so we first describe below an individual estimation procedure.

At stage $n$ of procedure $j$, let notation $C_{jn}=\{(y_{ji},\bm{{x}}_{ji}),i=1,...,n\}$ be the data set with the observations recruited from the original data set—say, $\mathcal{D}$. Suppose that $\hat{\bm{{\beta}}}_{jn}$ is a conventional least square estimate (LSE) of $\bm{{\beta}}_{0}$, based on $C_{jn}$, of model (1); additionally, let

[TABLE]

be an estimate of $\sigma^{2}$. With the notations defined above, we now describe the $j$th sequential estimation procedure as follows.

Let $\hat{\bm{{\beta}}}_{jn_{0}}$ be an LSE based on the initial data set $C_{jn_{0}}=\{(y_{ji},\bm{{x}}_{ji}),i=1,\ldots,n_{0}\}$ of size $n_{0}$. Let $\mu_{jn}=\lambda_{min}[(\bm{{X}}_{jn}\bm{{X}}_{jn}^{\top})/n]$; $\bm{{X}}_{jn}$ is a design matrix based on a set of covariate vectors, $\{\bm{{x}}_{j1},...,\bm{{x}}_{jn}\}$, of the observations in $C_{jn}$. (Notation $\lambda_{min}(A)$ denotes the minimum eigenvalue of matrix $A$.) Let $a^{2}$ be the $1-\alpha$ quantile of the chi-square distribution with degree of freedom $p$ and

[TABLE]

where $\tilde{a}_{j}>0$. For a given $d>0$, define a stopping time

[TABLE]

Hence, $N_{j}$ is the smallest number of observations, such that the inequality in (2) is fulfilled. Let $S_{N_{j}}=(\bm{{z}}-\hat{\bm{{\beta}}}_{jN_{j}})^{\top}(\bm{{X}}_{jN_{j}}\bm{{X}}_{jN_{j}}^{\top})(\bm{{z}}-\hat{\bm{{\beta}}}_{jN_{j}})$, and $\bm{{z}}=(z_{1},\cdots,z_{p})^{\top}$ and

[TABLE]

When the $j$th procedure is stopped according to (2), it is shown that $R_{N_{j}}$ is a confidence ellipsoid of $\bm{{\beta}}_{0}$ with the length of a maximum axis no greater than $2d$ and the coverage probability equals approximately $1-\tilde{\alpha}_{j}$, where $\tilde{\alpha}_{j}$ is the probability corresponding to $\tilde{a}_{j}$, which is the $1-\tilde{\alpha}_{j}$ quantile of $\chi^{2}(p)$. The coverage probability of $R_{N_{j}}$ depends on both $\tilde{a}_{j}$ and $d$. The constant $d>0$ serves as a precision index, which restricts the size of the confidence ellipsoid; this can be determined by users, based on their practical needs. When $d$ becomes smaller, we will need more observations, such that the estimate can achieve the desired precision. There have been many discussions in the literature regarding this conventional sequential estimation procedure; hence we only briefly touch upon technical matters, in Appendix A.

In a conventional sequential confidence set estimation procedure (i.e. $M=1$), we first specify the coverage, and then use its corresponding quantile and the size of the confidence set to define the stopping time. In the current method, we do not specify the coverage probability for each procedure $j$, $j=1,\ldots,M$. Instead, we only specify the coverage probability for the integrated one, which is a combination of all $M(>1)$ procedures. Thus, $\tilde{a}_{j}$ depends on the total number of procedures $M$, and for each $j$, the coverage of $R_{N_{j}}$ is less than $1-\alpha$. We show that the integrated one will have the desired coverage probability. This is why we demand few observations in each procedure; it also requires less computation time.

2.2 Combining several sequential estimation procedures

Suppose we have enough computing capacity available to simultaneously support $M$ sequential estimation procedures. When conducting a sequential procedure, we usually start with a small initial data set; technically, each procedure can have a different size of initial data set. However, to simplify our discussion, we assume that all initial data sets are of the same size ($n_{0}$), and we let $C_{jn_{0}}$, $j=1,\ldots,M$ be the initial data set for each of them.

Assume that there are already $k-1$ observations used in the estimation procedure $j$ ($k>n_{0}$). If the estimates obtained from procedure $j$ based on $C_{j{k-1}}$ cannot fulfil the inequality in (2), then we will select a new observation from $\mathcal{D}$; from there, we will calculate the new estimates, $\hat{\bm{{\beta}}}_{jk}$ and $\hat{\sigma}^{2}_{jk}$, based on $C_{j{k}}$. We repeat this sequential recruiting process until the inequality in (2) is satisfied.

Let $N^{*}=\sum_{j=1}^{M}N_{j}$ and $\rho_{j}=N_{j}/N^{*}$, with $j=1,\ldots,M$. When all $M$ procedures are stopped, let

[TABLE]

be a weighted average of $\{\hat{\bm{{\beta}}}_{jN_{j}}:j=1\ldots,M\}$, where $\sum_{j=1}^{M}\rho_{j}=1$. By definition, we know that $N_{j}$ is a function of $d$ and $\tilde{a}$, which are the same for all procedures in our current setup. Hence, $N_{j}$ values are random variables with the same mean (see Appendix B, (A.5) for the proof of Theorem 1). Thus, the actual number of observations used in each of the various procedures is different. Hence, we take into account the number of real observations used in each procedure, and we define $\hat{\bm{{\beta}}}$ as the weighted average of estimates obtained from each procedure with weights equal to $\rho_{j}=N_{j}/N^{*}$.

Then, based on the total number of observations $N^{*}$, we can define a confidence set for $\bm{{\beta}}_{0}$ as follows:

[TABLE]

where $\mu_{N^{*}}=\sum_{j=1}^{M}\rho_{j}\mu_{jN_{j}}$, $\mu_{jN_{j}}=\lambda_{min}[(\bm{{X}}_{jN_{j}}\bm{{X}}_{jN_{j}}^{\top})/N_{j}]$, $j=1,...,M$. We show that $R_{N^{*}}$ has the desired coverage probability and estimation precision. (Detailed arguments are provided in Appendix B.) We also prove that this distributed sequential estimation procedure retains some of the conventional statistical properties of sequential fixed-sized confidence set estimation, and summarize these properties as per the following theorem. (The proof of Theorem 1 is in Appendix B.)

Theorem 1

*Assume that for each $j$, $\{(x_{ji},y_{ji}),i\geq 1\}$ follows the linear regression model (1), and

(A1) $\lim_{n\rightarrow\infty}\sum_{i=1}^{n}\bm{{x}}_{ji}\bm{{x}}_{ji}^{\top}/n=\Sigma$, where the matrix $\Sigma$ is positive-definite.

Let $N_{j}$ be defined as in (2). Then,*

• (1)

$\lim_{d\rightarrow 0}\frac{d^{2}N^{*}}{a^{2}\sigma^{2}\mu}=1~{}~{}~{}~{}\text{\rm almost surely}$,

• (2)

$\lim_{d\rightarrow 0}P(\bm{{\beta}}_{0}\in R_{N}^{*})=1-\alpha$,

• (3)

$\lim_{d\rightarrow 0}\frac{d^{2}E(N^{*})}{a^{2}\sigma^{2}\mu}=1$,

where $\mu$ is the minimum eigenvalue of matrix $\Sigma$.

Theorems 1 (2) and (3) state that the proposed method has the prespecified coverage probability $1-\alpha$, and that the ratio of total data used in all procedures to the ‘theoretical’ one—the best but unknown one—is asymptotically equal to 1. Chow and Robbins (1965) first named these two properties asymptotic consistency and asymptotic efficiency, respectively, to describe the asymptotic properties of a sequential estimation procedure.

Under the current setup, the conventional sequential estimation is a case with $M=1$. The stopping time in this case is

[TABLE]

Obviously, $a^{2}>\tilde{a}_{j}^{2}$ implies that $N$ is stochastically larger than $N_{j}$ for each $j$, and because we run several procedures concurrently, the computation time will be smaller (as expected). However, it is worthwhile to note that when the number of variables increases or the procedure involves variable selection, the differences in computation times will increase. In addition, because the proposed procedure combines the estimation results of each procedure (each of which features a different data size), it is more stable than that of a single procedure in terms of coverage probabilities and the like. In our numerical studies, we will compare performance levels by using different $M$ values.

2.3 Active learning with D-optimal design criterion

In the machine learning literature, the term ‘active learning’ usually refers to some learning algorithms that can interactively query users. From statistical modelling perspectives, this concept approximates subject selection. The work of Deng et al. (2009)—which extends that of Wu (1985) and closely relates to that of Robbins and Monro (1951); Lai and Robbins (1979)—gives an example of how the sequential experimental design can play a role in an active learning procedure. In particular, Lai and Wei (1982) shows that the LSE for linear regression models with adaptive designs has some nice asymptotic properties under rather general design conditions, as follows.

• (A2) The random error $\{\epsilon_{n}\}$ is a martingale difference sequence with respect to an increasing sequence of $\sigma$-fields $\{{\cal{F}}_{n}=\sigma\{(x_{j},y_{j}):j=1,\ldots,n\}\}$ with

[TABLE]

• (A3) The maximum and minimum eigenvalues of design matrix $\sum_{i=1}^{n}x_{i}x_{i}^{T}$ satisfy, with a probability of 1, that

[TABLE]

These conditions are very general, and in their paper, they do not refer to any particular design. Following their results, we will show later that if we adopt the D-optimality of statistical experimental design methods to locate new observations from the data pool, these conditions will be satisfied. Hence, with the previously defined stopping criterion, the parameter estimates will have asymptotic consistency and efficiency, as stated below.

Theorem 2

Assume that conditions in Theorem 1, (A2), and (A3) hold. Let $N_{j}$ be defined as in (2), where new samples are selected by the aforementioned D-optimality. Then, we have $\lim_{d\rightarrow 0}(d^{2}N^{*})/(a^{2}\sigma^{2}\mu)=1$ almost surely, $\lim_{d\rightarrow 0}P(\beta_{0}\in R_{N}^{*})=1-\alpha$ and $\lim_{d\rightarrow 0}[d^{2}E(N^{*})]/(a^{2}\sigma^{2}\mu)=1$, where $\mu$ is the minimum eigenvalue of matrix $\Sigma$.

Theorem 2 simply says that with adaptive designs that use the $D$-optimality criterion, we still have asymptotic properties similar to those in Theorem 1. (The proof of this theorem is in Appendix C.)

We know that the properties of the estimates will depend on the design method. The D-optimal design criterion selects observations which maximize the determinant of the Fisher information matrix (or, equivalently, minimize the volume of the confidence ellipsoid for the parameters) (Neyer, 1994; Silvey, 1980). Because the selection is based on the current estimates of the unknown parameter for each procedure, the candidate observations for each procedure are usually different. In addition, we adopt a ‘without replacement’ principal, where each observation will be recruited by one procedure only.

2.4 Computation matter

Active learning is a special case of machine learning, and it is usually operated in a sequential way. From a statistical viewpoint, this type of method approximates a general sequential method and a statistical experimental design. When applying active learning methods under ‘big data’ scenarios, we usually assume that the data set used in the analysis already exists. Hence, how to efficiently find informative observations sequentially, such that we can efficiently and effectively conduct the analysis, is an important research problem. In this situation, design criteria are usually used to select the promising candidates from the existing data, rather than randomly select them. The work of Deng et al. (2009) is an example.

Here we do not actually conduct a sampling process when recruiting a new observation; instead, we only draw observations from a data storage device. The sampling cost here is reduced; however, the computation cost is still noticeable, since we need to iteratively conduct the estimation process for each procedure time until its stopping criterion is fulfilled. That is, we need to compute the determinant of the matrix, repeatedly. To this end, we consider the following two computation strategies.

Iterative formulas

It is clear that in such a sequential estimation procedure, the calculations of the determinant and inverse of the design matrix demand considerable computation power. However, with help of simple algebra, both can be calculated based on their predecessors. The following properties allow us to update the determinant of the new design matrix, when we sequentially add a new observation, with simple linear algebra. These two formulas are common in textbooks, such as Rao (1973); for convenience, we describe these below.

Let $\mathbf{A}$ be a nonsingular $p\times p$ matrix and $\mathbf{U}\in R^{p}$ be a column vector. Then

[TABLE]

Moreover, we also know that

[TABLE]

Equation (6) suggests that we can calculate the inverse of $(\mathbf{A}+\mathbf{U}\mathbf{U}^{T})$ when we know the inverse of $\mathbf{A}$; Equation (7) suggests a way of calculating the determinant of $\mathbf{A}+\mathbf{U}\mathbf{U}^{T}$ when we know both the determinant and inverse of $\mathbf{A}$. For the sequential estimation procedures discussed here, if we are at the $k$-th stage, then $\mathbf{A}=\left(\sum_{i=1}^{k}\bm{{x}}_{ji}\bm{{x}}_{ji}^{\top}\right)$ and $\mathbf{U}=\bm{{x}}_{j(k+1)}$. Because we need to repeatedly calculate the inverse and determinant of the design matrix when we recruit a new observation into the procedure at each stage, these two formulas are useful, as they save us a considerable amount of computation time.

Approximation method

Let $\rho_{j}=N_{j}/N^{*}$, as before; then, with a probability of 1

[TABLE]

as $d$ goes to 0. Based on this fact, we define an alternative confidence set for $\beta_{0}$ as follows.

[TABLE]

Obviously, it is easy to calculate (8) when the number of dimensions $p$ is large, compared to (4). In addition, we show that the length of the maximum axis of the ellipsoid $R^{{}^{\prime}}_{N^{*}}$ is no greater than $2d$, and that its coverage probability is approximately equal to the nominated $1-\alpha$ when the sample size is sufficiently large (or, equivalently, when $d$ is small) (see Appendix B). Hence, besides the aforementioned iterative formulas, these asymptotic results suggest that we can use $R^{{}^{\prime}}_{N^{*}}$ as an alternative confidence set, and that they provide us another way of reducing the computation cost of calculating the inverse of the high-dimensional matrix.

On the other hand, we know that

[TABLE]

where $I$ is an identical matrix with rank $p$. This implies that

[TABLE]

To have the same coverage probability, the confidence set based on (8) requires more observations. Therefore, we suggest that users use a confidence set based on (8) only for models with a lengthy parameter list, and/or when the cost of computation is high. Otherwise, we recommend using the confidence region (4), especially when $p$ is not overly large.

Remark 1

In most ‘big data’ analysis scenarios, data have already been collected and are stored in some digital device, such that we can easily manage and access the data: we simply ‘pick out’ those desired observations from the available data pool, based on a predetermined selection criterion (in our method, $D$-optimality). The ‘sampling’ cost here is lower than that in conventional situations, since there is no need to conduct new experiments or data collection processes.

3 Adaptive shrinkage estimation

Model interpretation is essential in most data analysis scenarios. Users always prefer to have solutions that are interpretable and understandable. Additionally, researchers can always benefit from the interpretability of a model to validate and/or further improve their hypotheses. By adding the ASE feature to the distributed sequential estimation procedure, we can better detect variables that significantly affect the model, and thereby increase the interpretability of the model. This is especially the case when the data set at hand contains a long list of variables. In this section, we provide a definition for ASE, such that the proposed procedure can detect the effective variables for the model during its estimation process.

Suppose that only $p_{0}$ of $p$ components of $\bm{{\mathbf{\beta}}}_{0}$ are effective to the model, and that we want to simultaneously identify these variables and construct a fixed-size confidence set with a prescribed coverage probability for them. Let $\hat{\bm{{\beta}}}_{j}$ be the LSE of $\bm{{\beta}}_{0}$ in procedure $j$, as before, and let $\lambda_{k}=\lambda|\hat{\beta}_{j}^{k}|^{-\gamma}$, where $\hat{\beta}_{j}^{k}$ is the $k$th element of $\hat{\bm{{\beta}}}_{j}$. In addition, let $\lambda\equiv\lambda(n)$ be a nonrandom function of $n$, such that for some $0<\delta<1/2$ and $\gamma>0$,

[TABLE]

Let $\epsilon>0$ be constant; then, define $I_{j}^{k}(\epsilon)=I(\sqrt{n}\lambda_{k}<\epsilon)$, for $k=1,\ldots,p$ as an indicator function, and let $\bm{{I}}_{j}(\epsilon)=\text{diag}\{I_{j}^{1}(\epsilon),\cdots,I_{j}^{p}(\epsilon)\}$ be a $p\times p$ diagonal matrix. Then, $\bm{{\beta}}_{j}^{*}\equiv\bm{{I}}_{j}(\epsilon)\hat{\bm{{\beta}}}_{j}$ is an ASE of $\mathbf{\beta}_{0}$, where ${\beta}_{j}^{*k}=0$ if $I_{j}^{k}(\epsilon)=0$; otherwise, it remains the same as $\hat{\beta}_{j}^{k}$.

Suppose that $\hat{\sigma}_{jk}^{2}$ and $\hat{p}_{0}(k)=\sum_{l=1}^{p}I_{jl}(\epsilon)$ are estimates of $\sigma_{jk}^{2}$ and $p_{0}$, respectively, based on the observations in $C_{jk}$ at the $k$th stage. Let $\chi_{\hat{p}_{0}(k)}(\alpha)$ be an $1-\alpha$ quantile of the chi-square distribution with $\hat{p}_{0}(k)$ degrees of freedom, and let $\tilde{a}_{jk}^{2}=\chi_{\hat{p}_{0}(k)}(\alpha)/M$. Then, define a stopping rule for procedure $j$ as follows:

[TABLE]

where $\tilde{\mu}_{jk}=\lambda_{max}[k\bm{{I}}_{j}(\epsilon)(\bm{{X}}_{jk}\bm{{X}}_{jk}^{\top})^{-1}\bm{{I}}_{j}(\epsilon)]$, and $\bm{{X}}_{jk}$ is the design matrix under the current stage.

We conduct procedure $j$ as before, until the stopping criterion $\tilde{N}_{j}$ is satisfied. Please note that $\tilde{a}_{jk}^{2}$ in (12) is now random and depends on the estimated number of effective variables at the current stage. Because ${a}_{k}^{2}$ in (2) is a constant calculated with a total number ($p$) of variables, this implies that $\tilde{a}_{jk}^{2}\leq{a}_{k}^{2}$ almost surely. Therefore, ${N}_{j}$ is stochastically larger than $\tilde{N}_{j}$; this is especially the case when $p_{0}$ is much smaller than $p$. These phenomena will also appear in our numerical studies.

3.1 Adaptive shrinkage estimation in distributed sequential estimation

When we adopt the ASE in all procedures, the major concern will be whether or not the selected variables from each procedure will be consistent.

For each $j=1,\ldots,M$, let $\tilde{N}_{j}$, $\bm{{\beta}}_{j}^{*}$, and $\bm{{I}}_{j}(\epsilon)$ be defined as before. Let $\tilde{N}=\sum_{j=1}^{M}\tilde{N}_{j}$ be the sum of the stopping time: when all procedures are stopped, it is equal to the sum of observations used in each procedure. Let

[TABLE]

be an element-wise product of the indicator vectors and define the combined ASE

[TABLE]

as a weighted average of the ASE $\bm{{\beta}}_{j}^{*}$ of each procedure with random weights $\rho_{j}=\tilde{N}_{j}/\tilde{N}$, $j=1,\ldots,M$. We know that for each $j$,

[TABLE]

where $\bm{{I}}_{0}=\text{diag}\{I(\beta_{01}\neq 0),...,I(\beta_{0p}\neq 0)\}$. Because each observation will be recruited into one procedure only, $\bm{{\beta}}_{j}^{*},j=1,...,M$, are independent. This implies that

[TABLE]

Hence, if $\beta_{0l}=0$, for some $1\leq l\leq p$, $\sqrt{\tilde{N}}(\tilde{\beta}_{l}-\beta_{0l})$ will degenerate to the constant 0 when the sample size is large ($d$ tends to 0); otherwise, it retains the same distribution as that of its corresponding component in the LSE. Thus, Equation (13) implies that as $d\rightarrow 0$,

[TABLE]

where notation $\bm{{A}}^{-}$ denotes a general inverse of matrix $\bm{{A}}$.

Let $\bm{{O}}$ be an orthonormal matrix, depending on the current ASE estimate, such that $(\tilde{\bm{{\beta}}}_{1},\tilde{\bm{{\beta}}}_{2})^{\top}=\bm{{O}}\tilde{\bm{{\beta}}}$ and $\bm{{O}}\bm{{O}}^{\top}=I_{p}$, where $\tilde{\bm{{\beta}}}_{2}=0$ and $\tilde{\bm{{\beta}}}_{1}$ includes all nonzero elements of $\tilde{\bm{{\beta}}}$. We show in Appendix D that

[TABLE]

where $\tilde{\bm{{\Sigma}}}_{11}$ is defined as (A.13) in Appendix D. Based on this, we define a confidence set of $\tilde{\bm{{\beta}}}$ as follows.

[TABLE]

where $S_{\tilde{N}}=(\bm{{O}}\bm{{Z}}-\bm{{O}}\tilde{\bm{{\beta}}})^{\top}\tilde{\bm{{\Sigma}}}(\bm{{O}}\bm{{Z}}-\bm{{O}}\tilde{\bm{{\beta}}})$, $\nu_{\tilde{N}}=\lambda_{max}\left(\tilde{N}\bm{{I}}^{*}\left\{\sum_{j=1}^{M}\bm{{X}}_{jN_{j}}\bm{{X}}_{jN_{j}}^{\top}\right\}^{-1}\bm{{I}}^{*}\right)$ and

[TABLE]

Using similar arguments, it is easy to see that the maximum axis of the confidence set $R_{\tilde{N}}$ is less than $2d$, and that the distributed ASE-based sequential procedure has the following properties. (The proofs for these properties are given in Appendix D.)

Theorem 3

Let $\tilde{N}_{j}$ be the stopping time defined in (12). Assume that the conditions of Theorem 1 are satisfied; then, (1) $\lim_{d\rightarrow 0}P(\bm{{\beta}}_{0}\in R_{\tilde{N}})=1-\alpha$, (2) $\lim_{d\rightarrow 0}{(d^{2}\tilde{N})}/{(a^{2}\sigma^{2}\tilde{\mu})}=1$ almost surely, and (3) $\lim_{d\rightarrow 0}{\{(d^{2}E(\tilde{N})\}}/{(a^{2}\sigma^{2}\tilde{\mu})}=1$, where $a^{2}$ is the $1-\alpha$ quantile of $\chi_{p_{0}}^{2}$ and $\tilde{\mu}$ is the maximum eigenvalue of matrix ${\bm{{I}}_{0}}\Sigma^{-1}{\bm{{I}}_{0}}$. Moreover, $\lim_{d\rightarrow 0}\hat{p}_{0}=p_{0}$ almost surely and $\lim_{d\rightarrow 0}E(\hat{p}_{0})=p_{0}$.

Remark 2

The divide-and-conquer method is a general concept and it can usually be applied to most procedures. However, there is no fixed sample size procedure that can concurrently guarantee the accuracy (i.e. size of confidence set) and the coverage probability of the estimate. Hence, we cannot apply the same idea to any fixed sample size procedure.

4 Numerical studies

We apply our methods to an appliance energy consumption data set and two PM2.5 data sets collected in two major cities in China (i.e. Beijing and Shanghai). Before discussing the real-data analysis, we first use some synthesized data sets containing 500 runs for each case, to show the performance of the proposed methods and compare their results to that of the conventional sequential method (i.e. $M=1$), the latter of which we refer to as ‘SM’. Additionally, we offer the results of multiple parallel sequential estimation procedures ($M>1$), which we refer to as ‘PSM’.

4.1 Simulation studies

Simple distributed sequential estimation

We first consider $M=2,5$ with $d=0.5,0.4,0.3$ and $0.2$, without using the ASE feature. We use two different true parameter vectors: $\bm{{\beta}}_{0}=(-1,1)$ and $(-1.0,1.0,0.7,0.5,0.2)$, denoted by $S_{1}$ and $S_{2}$, and generate data using the following model.

[TABLE]

where $\mathbf{X}=(X_{1},...,X_{p})$ with $X_{1}=1$ and $X_{i},i=2,...,p$, generated from normal distributions with mean $1$ and 0.2 for $p=2$ and $5$, and variance $1$; additionally, each component of $\epsilon$ follows the standard normal distribution.

Table 1 summarizes the total number of observations used (i.e. stopping times), empirical coverage probabilities, and computation time (i.e. CPU time). For PSM, $N^{*}$ is equal to the sum of the numbers of observations used in all procedures; the computation time is the maximum computation time for all procedures. There are two empirical coverage probabilities for PSM. The first is from the original confidence region (defined in (4)); the second is based on the alternative formula (defined in (8)). In this simple case, we can see from this table that the empirical coverage probabilities of all cases are close to the nominated $95\%$ coverage probability; we also see that they get closer as $d$ becomes smaller. For a given $d$, the total numbers of observations used are similar among SM and PSMs with $M=2,5$. The case with a higher dimensional vector ($\bm{{\beta}}_{0}$) tends to use larger data sizes and the computation times are as expected. We basically see that in this simple case, the proposed method can have the same desired properties as the conventional method.

Using the D-optimality criterion in subject selection

To show the procedure with D-optimality, we set only $M=2$, in addition to the conventional method. We generate data under the same setup as in the earlier case. The purpose of this part of the simulation study is to see the performance of using the D-optimality criterion for subject selection. To mimic the aforementioned ‘big data’ scenario, we first generate a set with $3000\times 2$ data points, as the data pools used in the PSM procedure (i.e. 3000 data points for each of them); the SM procedure will select the joint one with $6000$ data points. We summarize the results of this study in Table 2.

When we compare the current results to the previous ones, we see that the procedures featuring D-optimality as its selection criterion during the recruiting process can save more computation time: in this case, the ratio of the computation time of SM to that of PSM exceeds 3 1. This is due to the change in subject selection. The coverage is lower than the nominated $95\%$ for larger $d$ values; however, in the cases of both $S_{1}$ and $S_{2}$, they approach the nominated one as $d$ becomes smaller. The observation sizes in all cases are smaller than before; this too is due to the $D$-optimality subject selection criterion. Please note that in conventional experiments, after we know the ‘theoretically best’ design points, we need to collect new observations under such a specific setup; in some applications, this means conducting new experiments under some particular setup. However, in the current situation, we simply find the ‘best one’ among the existing data points. This is often the case in ‘big data’ scenarios, and this motivates us to develop such a method.

Distributed Sequential method versus divide-and-conquer

It is known that for a linear model when the variance is unknown there is no fixed sample size solution for constructing a confidence set for regression parameters with prescribed coverage probability and precision(Siegmund, 1985); this also implies that we cannot apply a simple divide and conquer (DC) method for such a problem. On the other hand, because the existence of nonhomogeneous data in modern large data collections is common, thus in spite of the fact that there is no solution available to the problem above via using the thinking of divide and conquer, we still conduct a numerical study using a contaminated data set below to compare the performances of the proposed method under such a messy data situation in terms of estimation errors.

We generate the major part of data set using the same setup as that in the the D-optimality criterion subsection. In addition, we add a small portions of data generated from a “wrong” regression model such that the rate of the contaminated data size to the regular data size, say $\rho$, is equal to $0.01,0.05,0.1$ and $0.15$. For S1 and S2, we respectively generate two sets of noisy data using linear models with parameters $\bm{{\beta}}=(-5,5)$ and $(-5,5,5,5,5)$. For the estimate baed on the divide and conquer method, we always use the whole data to estimate the regression parameters. Thus, the data size used in DC procedure is always equal to $6000$. In this case, we require no data selection scheme.

We report both the square error $SE=\sum_{i=1}^{p}(\hat{\beta}_{i}-\beta_{0i})^{2}$ and absolute deviation $AD=\sum_{i=1}^{p}|\hat{\beta}_{i}-\beta_{0i}|$ to assess the performances of PSM and DC. Table 3 summarizes the results of PSM and DC, where $N$ is the total sample sizes used. Although, when $\rho=0.01$, DC has smaller SE and AD than PSM, the performance of DC decays very fast when $\rho$ increases. Table 3 shows that SEs and ADs of PSM are much smaller than those of DC cases. Please note that in Table 3, when there is no or only small portion of noisy data, the PSMs use only much smaller data sizes than that of DCs, and this is the reason why PM has larger SE and AD. However, even in such a contaminated data situation, the PSM can locate the “most informative” points from a data pool via analyzing the current observations in hands, which largely reduces the chance of adding “wrong” data into analysis and boost the analysis efficiency. Table 4 reports the performances of DC with different $M$. For a total sample size $N$ in data pool, each partition has $N/M$ samples. SE and AD (see Table 4). For fixed $N$, SE and AD from DC become larger when $M$ increases, especially with small $N=500$.

Adaptive shrinkage estimation for detecting variables

We use two vectors to show the variable-selection ability of ASE: \color[rgb]{0,0,0}\bm{{\beta}}_{0}=(-2,1,1.5,2,0,0,0,0,0,0) and $\bm{{\beta}}_{0}=(-2,2,2,2,0,\ldots,0)$. The second vector has a total of 50 elements. Only the first four elements of this vector are nonzero; the other 46 are zero elements. In a real-data situation, the so-called noneffective variables will not precisely equal [math]; hence, as mentioned in the previous section, we need to choose a cutting parameter, as described in Section 3 (see also Wang and Chang, 2013).

In this simulation study, we take $M=5$ and then compare the performance of PSM to that of SM (which also has the ASE feature). In addition to the previous tables, we will report $\hat{p}_{0}$ here. From Table 5, we first note that the estimates of $p_{0}$ from both SM and PSM are similar. However, the total number of PSM observations is larger than that with SM. This is because when there are only a few observations in the beginning of a sequential process, it is difficult to derive stable variable selection results. (PSM needs to use more observations to derive a stable estimate of the number of effective variables $p_{0}$.) Despite the increase in data size in the PSM, the computation time used with PSM is around one-fifth that used with SM. The empirical coverage probability of PSM is closer to the nominated $95\%$ level than that of $SM$. It is worth noting that in the second case—where the true $\bm{{\beta}}_{0}$ has 50 elements, of which only four are nonzero—the coverage frequencies of SM are often less than $90\%$ for all $d$ values. This confirms that use of the PSM is more advantageous in high-dimensional cases.

4.2 Real-data examples

We apply the proposed method to appliance energy use (Energy) data and PM2.5 data collected in Shanghai and Beijing, China. For simplicity, when we apply the proposed method to these data sets, we first randomly permute the whole data set and then divide it into $M$ partitions. In this analysis, we set $M=5$—that is, we will simultaneously conduct five sequential estimation procedures with $M$ data partitions, and we will apply the SM procedure to the data at once without partitioning. It is easy for us to undertake the partitioning process as we conduct our numerical study, and we will recruit data sequentially in the previously described manner.

Appliance energy prediction data

Data concerning appliance energy use are reported in Luis et al. (2017), who recorded energy use and other variables such as house temperature and humidity conditions at 10-min intervals for about 4.5 months. This data set merges the records from an automatic detecting device, some wireless sensors, and weather data downloaded from a public data set. This data set is also available from a UCI machine learning repository (Dua and Karra Taniskidou, 2017). There are 19735 records in this data set. We refer readers to their paper (Luis et al., 2017) for further details. We study the relationship between appliance energy use and 12 other variables—namely, temperature in the kitchen area (T1), humidity in the kitchen area (RH1), temperature in the living room area (T2), humidity in the living room area (RH2), temperature outside the building (T3), humidity outside the building (RH3), temperature in the ironing room (T4), humidity in the ironing room (RH4), pressure (Pres), wind speed (WS), visibility (V), and dew point (DEWP).

Particulate matter 2.5 data

Particulate matter (PM) is a general term that describes the mixture of solid particles and liquid droplets in the ambient air. The terminology ‘PM2.5’ refers to fine particulate matter—in other words, the mass per cubic metre of air of particles with a size (diameter) generally smaller than 2.5 micrometres ($\mu m$; 2.5 micrometres is equal to 1/400 of a millimetre). Recent studies show that long-term exposure to PM2.5 might increase age-specific mortality risk, particularly from cardiovascular causes.

The PM2.5 data sets used in the current study were collected in Beijing and Shanghai; they derive from hourly-based records from January 1, 2010 to December 31, 2015 (Liang et al., 2016). After deleting missing data values, there are 49,579 and 31,880 records remaining in the Beijing and Shanghai data sets, respectively. We apply the proposed methods to these two data sets while using PM2.5 concentration (ug/$m^{3}$) as a response variable with dew point (DEWP), humidity (HUMI), pressure (PRES), temperature (TEMP), cumulated wind speed (Iws), hourly precipitation (Prec), and cumulated precipitation (Iprec) as the model covariates.

Results

Tables 6–8 contain the results from using the aforementioned real data sets with procedures that use different ways of selecting new observations during sequential estimation processes—including random and $D$-optimal criterion-based selection—with $M=1$ or $5$, and $d=0.2,0.3,0.4,$ or $0.5$. Table 6 reports the numbers of observations and the computation time with different computation setups. Tables 7 and 8 state the parameter estimates for Energy and PM2.5, respectively.

When we applied our method to these real data sets, we found that both the PSM and SM methods use about the same number of observations. From Table 7, we see that the differences in the regression parameter estimates between SM and PSM are not statistically significant. As both SM and PSM use the same model with similar numbers of observations, these results are reasonable.

On the other hand, in these real-data cases, the PSM uses less computation time than does the SM. These results suggest that the sequential estimation procedure can readily benefit from the divide-and-conquer strategy: by putting them together, we can save a considerable amount of computation time and accelerate our data analysis (or model-fitting, in this case).

When we use the D-optimality criterion to select new observations in the estimation process, both methods use fewer observations than their random selection counterparts. However, the differences in the computation time increase. Both methods use similar total numbers of observations under this case; hence, each estimation process in the PSM procedure, with $M=5$, uses only around one-fifth of the observations that the SM procedure uses. Therefore, the growing differences are due to the time spent in searching for new observations close to the D-optimality criterion. It is for this reason that PSM with the D-optimality feature is more efficient than the SM method under the same selection scheme.

In Tables 9–11, we report the results when we add the ASE feature to the procedures. The total number of observations used in PSM is slightly smaller than that used in SM, and the computation time of PSM is dramatically lower than that of SM (see 9). Moreover, we can see that PSM with the ASE feature tends to select fewer variables, which is one of the reasons why PSM requires less computation time. From Table 9, we found that the total amount of data used by PSM is slightly smaller than that used by SM; however, the computation time of PSM is much smaller than that of SM. We can see that PSM with ASE functionality tends to choose fewer variables, which is one of the reasons why PSM requires less computation time. Tables 10 and 11 contain parameter estimates for all three data sets when we apply the ASE method. When applying the ASE method, we need to select the cutting parameters. The choice of cutting parameters should depend on actual and application needs. From a model interpretation perspective, having a parsimonious model is often advantageous, because it can provide clear information for use in further research or related applications.

5 Discussion and closing remarks

In this study, we applied the divide-and-conquer method to a sequential confidence estimation procedure for linear models. Together with the adaptive shrinkage estimation (ASE) methods, we were able to decide the effective variables for the models. In addition, we also used the D-optimal design criterion for selecting informative subjects for model building, both sequentially and adaptively, as in the active learning methods found in the machine learning literature. To use this kind of method, we chose the ‘most informative’ observations, based on analysis of the current data on hand without doing extra experiments or data collection. The proposed sequential procedure is a synergy of many useful ideas raised separately in other studies. Our numerical results show that the proposed method can effectively detect important variables for a model that requires estimation accuracy and less computation time; additionally, its adaptive sample selection feature makes it useful in scenarios where a large, precollected data set is available. Such situations are common in modern ‘big data’ analysis. We fully exploited the adaptive sequential sampling features, such that the proposed method performed better than conventional methods in terms of estimation accuracy and computation time, while making only a minor and flexible demand on computation facilities. We applied the proposed method to analyse one real-world data set concerning appliance energy consumption and two real-world data sets pertaining to particulate matter 2.5. Our method can detect within the data sets the important variables (from a lengthy variable list) and estimate those detected variables at a prescribed accuracy level. The features in the proposed methods are essential in addressing such problems and will provide useful information that can inform future research and/or policymaking. It is clear that we can apply this methodology to many other sequential procedures, such as sequential methods in generalized linear models and classification problems. We will report on those results elsewhere, at a later date.

Appendix

A. Properties of sequential estimation for the machine $j$

At the stopping time $N_{j}$ defined in (2), a confidence ellipsoid of $\bm{{\beta}}_{0}$ is

[TABLE]

The length of the maximum axis of this ellipsoid is

[TABLE]

where $\lambda_{max}(A)$ is the maximum eigenvalue of matrix $A$. By definition, $\mu_{jN_{j}}=\lambda_{min}[(\bm{{X}}_{jN_{j}}\bm{{X}}_{jN_{j}}^{\top})/{N_{j}}]$. Hence, $D=2d$. For the stopping time and confidence set for the machine $j$, we show some statistical properties as presented in the following Lemma.

Lemma 1

Assume that the conditions of Theorem 1 are satisfied, and ${N}_{j}$ is defined in (2). Then,

[TABLE]

where $\tilde{\alpha}_{j}$ satisfies $P(\chi_{p}^{2}>\tilde{a}_{j}^{2})=\tilde{\alpha}_{j}$, and $\mu$ is the minimum eigenvalue of matrix $\Sigma$.

Proof. Similar to Wang and Chang (2013), the proof of this lemma is straightforward and omitted here.

B. Properties of distributed sequential estimation

Since derivatives of distributed sequential estimation with $M>1$ machines are similar to those with $M=2$, hereafter in this section—and without loss of generality—let $M=2$.

Maximum axis of $R_{N^{*}}$: From the definition of the ellipsoid $R_{N^{*}}$ defined in (4), the length of the maximum axis of the set defined by

[TABLE]

is

[TABLE]

This easily shows that

[TABLE]

Hence, with $\rho_{j}=N_{j}/N^{*}$, we have

[TABLE]

Proof of Theorem 1: At the stopping time $N_{j}$, from Lemma 4 and for the confidence set $R_{N_{j}}$ defined in (3), we have

[TABLE]

where $\tilde{\alpha}_{j}$ satisfies $P(\chi_{p}^{2}>\tilde{a}_{j}^{2})=\tilde{\alpha}_{j}$.

From (A.3) and (A.5), we have for each $j$, as $d\rightarrow 0$,

[TABLE]

which, combined with $\tilde{a}^{2}_{1}+\tilde{a}^{2}_{2}=a^{2}$, show that

[TABLE]

It follows that

[TABLE]

As $d$ tends to 0, we know that $\sqrt{N_{j}}(\hat{\bm{{\beta}}}_{N_{j}}-\bm{{\beta}}_{0})$ has an asymptotic normal distribution and the variance of $\hat{\bm{{\beta}}}_{N_{j}}$ can be estimated with $(\bm{{X}}_{N_{j}}\bm{{X}}_{N_{j}}^{T})^{-1}\sigma^{2}$. Since the samples for the two machines are independent, the variance of $\hat{\bm{{\beta}}}$ can be estimated with

[TABLE]

It follows that as $d\rightarrow 0$,

[TABLE]

By definition, $\mu_{N^{*}}=\rho_{1}\mu_{1N_{1}}+\rho_{2}\mu_{2N_{2}}$, $\mu_{jN_{j}}=\lambda_{min}[(\bm{{X}}_{jN_{j}}\bm{{X}}_{jN_{j}}^{\top})/N_{j}]$, $j=1,2$, and $\mu=\lambda_{min}(\Sigma)$. Hence, as $d\rightarrow 0$, we have $\mu_{N^{*}}\rightarrow\mu$ and from ((1)), ${N^{*}d^{2}}/{(\sigma^{2}\mu_{N^{*}})}\rightarrow a^{2}$ almost surely. Therefore, (A.6) implies that

[TABLE]

Therefore, Theorem 1 is proved.

Maximum axis of $R^{{}^{\prime}}_{N^{*}}$: The length of the maximum axis of the ellipsoid $R^{{}^{\prime}}_{N^{*}}$ is

[TABLE]

Denoted by $E=\bm{{X}}_{1N_{1}}\bm{{X}}_{1N_{1}}^{\top}$ and $F=\bm{{X}}_{2N_{2}}\bm{{X}}_{2N_{2}}^{\top}$. It shows that

[TABLE]

which indicates that when $\rho_{j}=N_{j}/N^{*}$,

[TABLE]

Hence, $D^{{}^{\prime}}\leq 2d$.

C. Proof of Theorem 2

Under D-optimality, the new samples are selected with the maximized determinant of the information matrix. Then, from condition (A1) of Theorem 1, this implies that the minimum eigenvalue of $\bm{{X}}_{jn}\bm{{X}}_{jn}^{\top}$ still has an order of $n$, where $\bm{{X}}_{jn}$ is the design matrix under the current status (i.e. $(\bm{{x}}_{j1},...,\bm{{x}}_{jn})_{p\times n}$). Therefore, the proof of Theorem 2 follows arguments similar to those for Theorem 1 and Wang and Chang (2013, Theorem 8). Hence, the details are omitted here.

D. Properties of distributed sequential estimation with adaptive shrinkage estimation

Partition of design matrix: By simple computation, we have that

[TABLE]

According to $\tilde{\bm{{\beta}}}_{1}$ and $\tilde{\bm{{\beta}}}_{2}$, partition the matrix $\sum_{j=1}^{M}\bm{{O}}\bm{{X}}_{jN_{j}}(\bm{{O}}\bm{{X}}_{jN_{j}})^{\top}$ as follows,

[TABLE]

where $\bm{{\Sigma}}_{11}$, $\bm{{\Sigma}}_{12}$ and $\bm{{\Sigma}}_{22}$ are ${\hat{p}_{0}\times\hat{p}_{0}}$, ${\hat{p}_{0}\times(p-\hat{p}_{0})}$ and ${(p-\hat{p}_{0})\times(p-\hat{p}_{0})}$ matrices, respectively. Then, we have

[TABLE]

where

[TABLE]

Hence,

[TABLE]

is a general inverse matrix of the matrix defined on the left-hand side of (A.12). Consequently, we have

[TABLE]

Proof of Theorem 3: At the stopping time $\tilde{N}_{j}$, similar to Lemma 4, we have

[TABLE]

which indicates that for each $j$, as $d\rightarrow 0$,

[TABLE]

Hence, from $\tilde{a}^{2}_{1}+\tilde{a}^{2}_{2}=a^{2}$, it shows that

[TABLE]

From Wang and Chang (2013), we have as $d\rightarrow 0$,

[TABLE]

where $\bm{{I}}_{0}={\rm{diag}}\{I(\beta_{01}\neq 0),\cdots,\beta_{0p}\neq 0)\}$. Since the sequential sampling procedures for $M$ machines are independent, as $d\rightarrow 0$,

[TABLE]

which gives that

[TABLE]

Furthermore, it is easily shown that $\lim_{d\rightarrow 0}E(\hat{p}_{0})=p_{0}$. Similar to proof of Theorem 1, we have

[TABLE]

The Proof of Theorem 3 is finished.