Nothing
#' Orthogonal Latin Hypercube Design
#'
#' \code{OLHD.Y1998} returns a \code{2^m+1} by \code{2m-2} orthogonal Latin hypercube design generated by the construction method of Ye (1998)
#'
#' @param m A positive integer, and it must be greater than or equal to 2.
#'
#' @return If all inputs are logical, then the output will be an orthogonal LHD with the following run size: \code{n=2^m+1} and the following factor size: \code{k=2m-2}.
#'
#' @references Ye, K.Q. (1998) Orthogonal column Latin hypercubes and their application in computer experiments. \emph{Journal of the American Statistical Association}, \strong{93}(444), 1430-1439.
#'
#' @examples
#' #create an orthogonal LHD with m=3. So n=2^m+1=9 and k=2*m-2=4
#' OLHD.Y1998(m=3)
#'
#' #create an orthogonal LHD with m=4. So n=2^m+1=17 and k=2*m-2=6
#' OLHD.Y1998(m=4)
#'
#' @export
OLHD.Y1998=function(m){
if(m<2){
stop("m must be greater than or equal to 2")
}
#m >= 2
q=2^(m-1)
#construction of M starts
e=sample(seq(1,q),q) #first column
e=matrix(e,ncol=1)
I=matrix(c(1,0,0,1),ncol=2,nrow=2,byrow=T)
R=matrix(c(0,1,1,0),ncol=2,nrow=2,byrow=T)
AL=array(0,dim=c(q,q,m-1)) #there are m-1 of AL's
if(m==2){
AL[,,m-1]=R
M=cbind(e,AL[,,m-1]%*%e)
}
if(m>=3){
for (i in 1:(m-1-1)) { #This is the loop for L. However, When L=m-1, m-1-L=0, and AL cannot be
#calcultated within the loop. So L=1,...,m-2, and AL with L=m-1 will be calculated separately
a=1
for (j in 1:(m-1-i)) { #This is the loop for m-1-L, the first half of AL
a=kronecker(a,I)
}
b=1
for (k in 1:i) {
b=kronecker(b,R)
}
AL[,,i]=kronecker(a,b)
}
c=1
for (l in 1:(m-1)) {
c=kronecker(c,R)
}
AL[,,m-1]=c
M=e
for (i in 1:(m-1)) {
M=cbind(M,AL[,,i]%*%e)
}
for (i in 1:(m-2)) {
M=cbind(M,AL[,,i]%*%AL[,,m-1]%*%e)
}
}
#construction of M ends
#construction of S starts
j=rep(1,q)
j=matrix(j,ncol=1)
ak=array(0,dim=c(q,1,m-1)) #there are m-1 of ak's
B=array(1,dim=c(2,1,m-1))
if(m==2){
B[1,,m-1]=-1 #B_{m-k}
ak[,,m-1]=B[,,1]
S=cbind(j,ak[,,m-1])
}
if(m>=3){
for (i in 1:(m-1)) { #This is the loop for k.
temp=B
temp[1,,m-i]=-1 #B_{m-k}
d=1
for (k in 1:(m-1)) {
d=kronecker(d,temp[,,k])
}
ak[,,i]=d
}
S=j
for (i in 1:(m-1)) {
S=cbind(S,ak[,,i])
}
for (i in 2:(m-1)) {
S=cbind(S,ak[,,1]*ak[,,i])
}
}
#construction of S ends
#construction of T starts
T0=matrix(0,nrow=q,ncol=(2*m-2))
for (i in 1:q) {
for (k in 1:(2*m-2)) {
T0[i,k]=M[i,k]*S[i,k]
}
}
#construction of T ends
CP=matrix(0,ncol=(2*m-2),nrow=1) #center point
X=rbind(T0,CP,-T0)
X
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.