[R] latin hypercube sampling

[Rob Carnell]

Aimee Kopolow <alj27 <at> georgetown.edu> writes:
> 
> Hi all,
> 
> I am attempting to use latin hypercube sampling to sample different
> variable functions in a series of simultaneous differential equations.
> There is very little code online about lhs or clhs, so from different
> other help threads I have seen, it seems I need to create a
> probability density function for each variable function, and then use
> latin hypercube sampling on this pdf.
> 
> So far, I have created a data frame consisting of the "y" output of
> density(functionX) for each of the different functions I wish to
> sample. [examples of functions include T1(t), WL1(T1,t),
> BE1(WL1,T1,t)] The dataframe consists of 512 rows/vectors for each
> function.
> 
> I tried running
> res <- clhs(df, size = 500, iter = 2000, progress = FALSE, simple = TRUE)
> 
> and it returned a single series of 500 samples, rather than a series
> of 500 samples per function.
> 
> I ultimately need a sample of each variable function that I can run
> through my model, putting each individual variable function as a
> constant instead, and then performing PRCC. Is there anyone who can
> advise on how to do this, or failing that, where I should look for
> sample code?
> 
> Thank you for any help you are able to give,
> Aimee.
> 
> 

Aimee,

I'm the package maintainer for the lhs package.  Unfortunately, I'm not familiar
with the functions you mentioned (reproducible code would help us answer your
post).  I will try to show something parallel to what you described.

require(lhs)

# functions you described
T1 <- function(t) t*t
WL1 <- function(T1, t) T1*t
BE1 <- function(WL1, T1, t) WL1*T1*t

# t is distributed according to some pdf (e.g. normal)
require(lhs)
# draw a lhs with 512 rows and 3 columns (one for each function)
y <- randomLHS(512, 3)
# transform the three columns to a normal distribution (these could be any
# distribution)
t <- apply(y, 2, function(columny) qnorm(columny, 2, 1))
# transform t using the functions provided
result <- cbind(
  T1(t[,1]),
  WL1(T1(t[,2]), t[,2]),
  BE1(WL1(T1(t[,3]), t[,3]), T1(t[,3]), t[,3])
)
# check the results
# these should be approximately uniform
windows()
par(mfrow=c(2,2))
apply(y, 2, hist, breaks=50)
# these should be approximately normal
windows()
par(mfrow=c(2,2))
apply(t, 2, hist, breaks=50)
# these should be the results of the functions
windows()
par(mfrow=c(2,2))
apply(result, 2, hist, breaks=50)

Please feel free to contact me as the package maintainer if you need additional
help with lhs.

Rob