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-boun...@r-project.org> de la part de Rob C
<bertcarn...@gmail.com>
Envoyé : samedi 27 mai 2017 13:32:23
À : r-help@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@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.
[[alternative HTML version deleted]]
______________________________________________
R-help@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.
