

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
[[alternative HTML version deleted]]
______________________________________________
[hidden email] mailing list  To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/rhelpPLEASE do read the posting guide http://www.Rproject.org/postingguide.htmland provide commented, minimal, selfcontained, reproducible code.


>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
______________________________________________
[hidden email] mailing list  To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/rhelpPLEASE do read the posting guide http://www.Rproject.org/postingguide.htmland provide commented, minimal, selfcontained, reproducible code.


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 : Rhelp < [hidden email]> de la part de Rob C < [hidden email]>
Envoyé : samedi 27 mai 2017 13:32:23
À : [hidden email]
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
______________________________________________
[hidden email] mailing list  To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/rhelpPLEASE do read the posting guide http://www.Rproject.org/postingguide.htmland provide commented, minimal, selfcontained, reproducible code.
[[alternative HTML version deleted]]
______________________________________________
[hidden email] mailing list  To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/rhelpPLEASE do read the posting guide http://www.Rproject.org/postingguide.htmland provide commented, minimal, selfcontained, reproducible code.


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 < [hidden email]> 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 : Rhelp < [hidden email]> de la part de Rob C
> < [hidden email]>
> Envoyé : samedi 27 mai 2017 13:32:23
> À : [hidden email]
> 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
>
> ______________________________________________
> [hidden email] mailing list  To UNSUBSCRIBE and more, see
> https://stat.ethz.ch/mailman/listinfo/rhelp> PLEASE do read the posting guide http://www.Rproject.org/postingguide.html> and provide commented, minimal, selfcontained, reproducible code.
______________________________________________
[hidden email] mailing list  To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/rhelpPLEASE do read the posting guide http://www.Rproject.org/postingguide.htmland provide commented, minimal, selfcontained, reproducible code.


Thank you very much Rob for your answer. I have some difficulties to understand how to apply my agentbased model to each parameter combination generated by the LHS, in particular when parameters are defined by probability distributions. Indeed, I have multiple parameters in my model: parameters which are defined by a single value (like “temperature", "pressure”) and parameters which are defined by probability distributions (like “dispersal distance”). It’s correct for me to treat distance as a class. When all parameters are defined by a single value, I need first to create a data frame in which each column represents a different parameter, and each row represents a different combination of parameter values. Then, I apply my model to each row of the data frame. But, it’s not clear for me how to do this when parameters are defined from probability distributions? In particular, how can I use your code to apply my model to each of the 50 rows of the data frame “tabLHS”? Given that one row corresponds to one model simulation, I should have a value generated by the LHS for all distance classes at the first line of the data frame.
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)
tabLHS < get.data(uncoupledLHS)
Sorry, it’s the first time that I perform a sensitivity analysis using the LHS.
Thank you very much for your time.
Have a nice day
Nell
________________________________
De : Rob C < [hidden email]>
Envoyé : mardi 30 mai 2017 16:26:08
À : Nelly Reduan
Cc : [hidden email]
Objet : Re: [R] Latin Hypercube Sampling when parameters are defined according to specific probability distributions
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 < [hidden email]> 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 : Rhelp < [hidden email]> de la part de Rob C
> < [hidden email]>
> Envoyé : samedi 27 mai 2017 13:32:23
> À : [hidden email]
> 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
>
> ______________________________________________
> [hidden email] mailing list  To UNSUBSCRIBE and more, see
> https://stat.ethz.ch/mailman/listinfo/rhelp> PLEASE do read the posting guide http://www.Rproject.org/postingguide.html> and provide commented, minimal, selfcontained, reproducible code.
[[alternative HTML version deleted]]
______________________________________________
[hidden email] mailing list  To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/rhelpPLEASE do read the posting guide http://www.Rproject.org/postingguide.htmland provide commented, minimal, selfcontained, reproducible code.


I think you should take this conversation private or seek local
statistical expertise. This is about strategies for analysis, not
about programming in R.
Cheers,
Bert
Bert Gunter
"The trouble with having an open mind is that people keep coming along
and sticking things into it."
 Opus (aka Berkeley Breathed in his "Bloom County" comic strip )
On Wed, May 31, 2017 at 9:49 PM, Nelly Reduan < [hidden email]> wrote:
> Thank you very much Rob for your answer. I have some difficulties to understand how to apply my agentbased model to each parameter combination generated by the LHS, in particular when parameters are defined by probability distributions. Indeed, I have multiple parameters in my model: parameters which are defined by a single value (like “temperature", "pressure”) and parameters which are defined by probability distributions (like “dispersal distance”). It’s correct for me to treat distance as a class. When all parameters are defined by a single value, I need first to create a data frame in which each column represents a different parameter, and each row represents a different combination of parameter values. Then, I apply my model to each row of the data frame. But, it’s not clear for me how to do this when parameters are defined from probability distributions? In particular, how can I use your code to apply my model to each of the 50 rows of the data frame “tabLHS”? Given that one row corresponds to one model simulation, I should have a value generated by the LHS for all distance classes at the first line of the data frame.
>
>
>
> 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)
>
> tabLHS < get.data(uncoupledLHS)
>
>
>
> Sorry, it’s the first time that I perform a sensitivity analysis using the LHS.
>
>
> Thank you very much for your time.
>
> Have a nice day
>
> Nell
>
>
> ________________________________
> De : Rob C < [hidden email]>
> Envoyé : mardi 30 mai 2017 16:26:08
> À : Nelly Reduan
> Cc : [hidden email]
> Objet : Re: [R] Latin Hypercube Sampling when parameters are defined according to specific probability distributions
>
> 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 < [hidden email]> 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 : Rhelp < [hidden email]> de la part de Rob C
>> < [hidden email]>
>> Envoyé : samedi 27 mai 2017 13:32:23
>> À : [hidden email]
>> 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
>>
>> ______________________________________________
>> [hidden email] mailing list  To UNSUBSCRIBE and more, see
>> https://stat.ethz.ch/mailman/listinfo/rhelp>> PLEASE do read the posting guide http://www.Rproject.org/postingguide.html>> and provide commented, minimal, selfcontained, reproducible code.
>
> [[alternative HTML version deleted]]
>
> ______________________________________________
> [hidden email] mailing list  To UNSUBSCRIBE and more, see
> https://stat.ethz.ch/mailman/listinfo/rhelp> PLEASE do read the posting guide http://www.Rproject.org/postingguide.html> and provide commented, minimal, selfcontained, reproducible code.
______________________________________________
[hidden email] mailing list  To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/rhelpPLEASE do read the posting guide http://www.Rproject.org/postingguide.htmland provide commented, minimal, selfcontained, reproducible code.

