multistagevariance: Expected variance after selection after k stages selection

Description

This function uses the algorithm described by Tallis (1961) to calculate the variance after multi-stage selection. The variance among candidates of y in the selected area $\textbf{S}_{Q}$ is defined as the second central moment, $\psi_n(y)=E(Y^2|\textbf{S}_{Q}) - \left[E(Y|\textbf{S}_{Q})\right]^2$, where $E(Y^2|\textbf{S}_{Q}) = \alpha^{-1} \int_{-\infty} ^\infty \int_{q_{1}}^\infty...\int_{q_{n}}^\infty y^2\, \phi_{n+1}(\textbf{x}^{*}; \bm{\Sigma}^{*}) \, d \textbf{x}^*$

Usage

multistagevariance(Q, corr, alg)

Arguments

are the coordinates of the truncation points, which are the output of the function multistagetp that we are going to introduce.

corr

is the correlation matrix of y and X, which is introduced in the function multistagecorr. The correlation matrix must be symmetric and positive-definite. If the estimated correlation matrix is negative-definite, it must be adjusted before using this funct

alg

is used to switch between two algorithms. If alg = GenzBretz(), which is by default, the quasi-Monte Carlo algorithm from Genz et al. (2009, 2013), will be used. If alg = Miwa(), the program will use the Miwa algorithm (Mi et al.

Value

The output is the value of $\psi_n(y|\textbf{S}_{Q})$.

References

A. Genz and F. Bretz. Computation of Multivariate Normal and t Probabilities. Lecture Notes in Statistics, Vol. 195, Springer-Verlag, Heidelberg, 2009. A. Genz, F. Bretz, T. Miwa, X. Mi, F. Leisch, F. Scheipl and T. Hothorn. mvtnorm: Multivariate normal and t distributions. R package version 0.9-9995, 2013. G.M. Tallis. Moment generating function of truncated multi-normal distribution. J. Royal Stat. Soc., Ser. B, 23(1):223-229, 1961. X. Mi, T. Miwa and T. Hothorn. Implement of Miwa's analytical algorithm of multi-normal distribution. R Journal, 1:37-39, 2009.

Examples

Run this code

# first example

Q =c(0.4308,0.9804,1.8603)

corr=matrix( c(1,       0.3508,0.3508,0.4979,
               0.3508,  1,     0.3016,0.5630,
               0.3508,  0.3016,1,     0.5630,
               0.4979,  0.5630,0.5630,1), 
              nrow=4  
)


multistagevariance(Q=Q,corr=corr,alg=Miwa)

# time comparsion

var.time.miwa=system.time (var.miwa<-multistagevariance(Q=Q,corr=corr,alg=Miwa))

var.time.bretz=system.time (var.bretz<-multistagevariance(Q=Q,corr=corr))



# second examples 


Q= c(0.9674216, 1.6185430)
corr=matrix( c(1,         0.7071068, 0.9354143,
               0.7071068, 1,         0.7559289,
               0.9354143, 0.7559289, 1), 
              nrow=3  
)


multistagevariance(Q=Q,corr=corr,alg=Miwa)

var.time.miwa=system.time (var.miwa<-multistagevariance(Q=Q, corr=corr, alg=Miwa))

var.time.bretz=system.time (var.bretz<-multistagevariance(Q=Q, corr=corr))




# third examples 

 alpha1<- 1/(24)^0.5
 alpha2<- 1/(24)^0.5
 Q=multistagetp(alpha=c(alpha1,alpha2),corr=corr)


corr=matrix( c(1,         0.7071068,0.9354143,
               0.7071068, 1,        0.7559289,
               0.9354143, 0.7559289,1), 
              nrow=3  
)

multistagevariance(Q=Q, corr=corr, alg=Miwa)