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

[Rob C]

Nell,

I still might not be interpreting your question correctly, but I will try.

When you say that the sum of the probabilities of all distance classes
must equal one, I am going to treat distance as a class, instead of as a
continuous variable.

distance_class_probabilities <- c(9, 11, 10, 8.9, 8, 7, 6, 5.8, 5.1,
4, 3.9, 3.7, 3.4, 3.1, 2, 1.9, 1.6, 1.4, 1, 0.9, 0.8, 0.7, 0.4, 0.3,
0.1)/100
sum(distance_class_probabilities)
distance_classes <- factor(paste0("d", 1:25), ordered = TRUE,
levels=paste0("d", 1:25))
input_parameters <- c("dispersal_distance", "temperature", "pressure")
N <- 1000

plot(1:25, distance_class_probabilities, type="h", lwd=5)

set.seed(1)
require(lhs)
X <- randomLHS(N, length(input_parameters))
dimnames(X) <- list(NULL, input_parameters)
Y <- X
Y[,"dispersal_distance"] <-
approx(x=cumsum(distance_class_probabilities), y=1:25,
xout=X[,"dispersal_distance"], method="constant", yleft=0)$y + 1

hist(Y[,"dispersal_distance"], breaks=c(seq(0.5, 25.5, by=1)))
plot(table(distance_classes[Y[,"dispersal_distance"]]))


Is it still a Latin hypercube?

This is a more difficult question.  In some ways, the sample is still a Latin
hypercube since it was drawn that way.  But once the sample has been discretized
into the distance classes, then it loses the latin property of having only one
sample per "row".  It might be close enough for your purposes though.

Rob



On Tue, May 30, 2017 at 10:59 AM, Nelly Reduan <nell.redu at hotmail.fr> wrote:
> Thanks a lot Rob for your answer. I need to add a condition for the
> parameter “dispersal distance”. The sum of the probabilities of all distance
> classes must be equal to 1:
>
> y <- c(9, 11, 10, 8.9, 8, 7, 6, 5.8, 5.1, 4, 3.9, 3.7, 3.4, 3.1, 2, 1.9,
> 1.6, 1.4, 1, 0.9, 0.8, 0.7, 0.4, 0.3, 0.1)
>
> x <- seq(1, 25, by = 1)
>
> barplot(y/100, names.arg=x, ylab="Probability", xlab="Distance (km)")
>
>
>
> With this condition, is it possible to perform a LHS?
>
> Thanks a lot for your time.
>
> Nell
>
>
> ________________________________
> De : R-help <r-help-bounces at r-project.org> de la part de Rob C
> <bertcarnell at gmail.com>
> Envoyé : samedi 27 mai 2017 13:32:23
> À : r-help at r-project.org
> Objet : Re: [R] Latin Hypercube Sampling when parameters are defined
> according to specific probability distributions
>
>>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
>
> ______________________________________________
> R-help at r-project.org mailing list -- To UNSUBSCRIBE and more, see
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.