[R] Latin Hypercube Sampling when parameters are defined according to specific probability distributions

[Rob C]

>May 26, 2017; 11:41am  Nelly Reduan Latin Hypercube Sampling when parameters are >defined according to specific probability distributions
>Hello,
> I would like to perform a sensitivity analysis using a Latin Hypercube Sampling (LHS).
>Among the input parameters in the model, I have a parameter dispersal distance which is defined according to an exponential probability distribution.

>In the model, the user thus sets a default probability value for each distance class.

>For example, for distances ([0  2]; ]2  4]; ]4  6]; ]6  8]; ]8  10];; ]48  50],

>respective probabilities are 0.055; 0.090; 0.065; 0.035; 0.045;; 0.005.

 >Here is the code to represent an exponential probability
distribution for the parameter dispersal distance:

>set.seed(0)
>foo <- rexp(100, rate = 1/10)
>hist(foo, prob=TRUE, breaks=20, ylim=c(0,0.1), xlab ="Distance (km)")
>lines(dexp(seq(1, 100, by = 1), rate = 1/mean(foo)),col="red")
>1/mean(foo)

>When a parameter is defined according to a specific probability distribution, how can I perform a LHS ?
>For example, should I sample N values from a uniform distribution for each distance class (i.e., [0 � 2]; ]2 � 4]; ]4 � 6]; ]6 � 8]; ]8 � 10];��; ]48 � 50])
>or sample N values from exponential distributions with different rates ?

>Here is the code used to perform a LHS when the parameter �dispersal distance� is defined by one default value in the model:

>library(pse)
>factors <- c("distance")
>q <- c("qexp")
>q.arg <- list( list(rate=1/30) )
>uncoupledLHS <- LHS(model=NULL, factors, 50, q, q.arg)
>head(uncoupledLHS)

>Thanks a lot for your time.
>Have a nice day
>Nell

Nell,

I would like to suggest a slightly different method for generating the
sample using the lhs library,  then I will try using the pse library.
Generally when you have a package specific
question, you should try to contact the package maintainer first.

set.seed(1)
# I don't think your model has only one parameter, so I will include multiple
input_parameters <- c("dispersal_distance", "temperature", "pressure")
N <- 50
exponential_rate <- 1/30

library(lhs)
X <- randomLHS(N, length(input_parameters))
dimnames(X) <- list(NULL, input_parameters)
# X is now a uniformly distributed Latin hypercube
head(X)
hist(X[,1], breaks=5)
hist(X[,2], breaks=5)
hist(X[,3], breaks=5)
# now, transform the dispersal_distance paramter to an exponential sample
Y <- X
Y[,"dispersal_distance"] <- qexp(X[,"dispersal_distance"],
rate=exponential_rate)
hist(Y[,1], breaks=10)
# you can transform the other marginals as required and then assess
function sensitivity
model_function <- function(z) z[1]*z[2] + z[3]
apply(Y, 1, model_function)

# now, trying to use pse
library(pse)
q <- list("qexp", "qunif", "qunif")
q.arg <- list(list(rate=exponential_rate), list(min=0, max=1),
list(min=0, max=1))
uncoupledLHS <- LHS(model=model_function, input_parameters, N, q, q.arg)
hist(uncoupledLHS$data$dispersal_distance, breaks=10)

Rob