########## x: one-dimensional data vector

REg1 <- function(x,g0=0.5,nstep=500,ediff=10^(-3)){
	mu0 = median(x); sd0 = mad(x);
	for(istep in 1:nstep){
		weight = dnorm(x,mu0,sd0)^g0;
		mu1 = weighted.mean(x,weight);
		var1 = ( weighted.mean(x^2,weight) - mu1^2 ) * (1+g0);
		sd1 = sqrt(var1);
		diff = abs(mu1-mu0) + abs(sd1-sd0);
		if( diff < ediff ) break;
		mu0 = mu1; sd0 = sd1;
	}
	if( istep == nstep ){
		print("NOT CONVERGENCE");
		print(diff);
	}
	c(mu1,sd1)
}



########## estimation of asymptotic variance of hlambda

AsyVarForLambdaS <- function(x,mu,sigma,g0=0.5){
	n = length(x);
	wc = x - mu; 
	v1 = mean(wc^2*dnorm(x,mu,sigma)^g0);
	v2 = mean(wc^2*dnorm(x,mu,sigma)^(2*g0));
	avar = (v2/v1^2)*(sigma^2)^2;
	return(avar)
}




########## estimation of variance sigma

AdjustMean <- function(x){
	wx = x;
	wa = apply(x,1,mean); wx = wx - wa;
	wa = apply(x,2,mean); wx = t( t(wx) - wa );
	wa = mean(x); wx = wx + wa;
	return(wx)
}


EstSigmaDirect <- function(x){
	wx = x; nj = nrow(x); nk = ncol(x)/2;
	wa = AdjustMean(x[,1:nk]);
	wb = AdjustMean(x[,nk+1:nk]);
	wc = sum(wa^2+wb^2)/(2*(nj-1)*(nk-1));
	sigma = sqrt(wc);
	return(sigma)
}



########## multiple test
##########   for a probe set

MulTest <- function(y,hlambda,hsigy,hsigma,nk,g0=0.5){
	##### necessary variables
	nj = length(y);
	####### t-value
	avar = AsyVarForLambdaS(y,hlambda,hsigy,g0);
	tv.gene = hlambda/sqrt(avar/nj);
	##### pvalue
	pvalue = NULL; 
	##### test for each probe
	for(ij in 1:nj){
		########## hlambdaj
		hlambdaj = y[ij];
		##### asymptotic variance of hlambdaj
		varj = 2*hsigma^2/nk;
		##### test statistic
		tj = (hlambdaj-hlambda)/sqrt(avar/nj+varj);
		##### p-value (upper: two-side test, lower: one-side test)
		per = 1 - pnorm(tj);
		pvalue = append( pvalue, per );
	}
	names(pvalue) <- NULL
	names(tv.gene) <- NULL
	c(pvalue,tv.gene)
}



########## multiple test
##########   for all the probes

MulTestForGenes <- function(xa,npset,g0=0.5,nstep=500,ediff=10^(-3)){
	pvalue = NULL;
	tv.gene = NULL;
	for(i in 1:npset){
		##### probe set data 
		ip0 = ProbeFirstNum[i]; ip1 = ProbeFirstNum[i+1];
		nj = ip1 - ip0;
		ij = ip0+0:(nj-1);
		x <- xa[ij,];
		##### necesarry variables
		np = ncol(x); nk = np/2;
		##### simple variable
		wa = apply(x[,1:nk],1,mean);
		wb = apply(x[,nk+1:nk],1,mean);
		y = wa - wb;
		##### estimation of parameter
		hsigma = EstSigmaDirect(x);
		est.diff = REg1(y,g0,nstep,ediff); 
		hlambda = est.diff[1];
		hsigy = max( est.diff[2], hsigma*sqrt(2/nk) );
		##### multiple test
		w = MulTest(y,hlambda,hsigy,hsigma,nk,g0);
		wa = length(w);
		pvalue = append(pvalue,w[-wa]);
		tv.gene = append(tv.gene,w[wa]);
	}
	c(pvalue,tv.gene)
}