, Wenping Wang [ctb],
Bill Denney [ctb] | Title: | Random Number Generation Functions for 'rxode2' |
|---|---|
| Description: | Provides the random number generation (in parallel) needed for 'rxode2' (Wang, Hallow and James (2016) <doi:10.1002/psp4.12052>) and 'nlmixr2' (Fidler et al (2019) <doi:10.1002/psp4.12445>). This split will reduce computational burden of recompiling 'rxode2'. |
| Authors: | Matthew L. Fidler [aut, cre]
, Wenping Wang [ctb],
Bill Denney [ctb] |
| Maintainer: | Matthew L. Fidler <[email protected]> |
| License: | GPL (>= 3) |
| Version: | 2.1.1 |
| Built: | 2024-09-29 06:05:16 UTC |
| Source: | https://github.com/cranhaven/cranhaven.r-universe.dev |
cbind Ome
.cbindOme(et, mat, n).cbindOme(et, mat, n)
et |
The theta data frame |
mat |
The full matrix simulation from omegas |
n |
number of subject simulated |
data frame with et combined with simulated omega matrix values
Matthew Fidler
Convert numeric vector to repeated data.frame
.vecDf(vec, n).vecDf(vec, n)
vec |
Named input vector |
n |
Number of columns |
Data frame with repeated vec
Matthew Fidler
Note this Inverse wishart rescaled to match the original scale of the covariance matrix.
cvPost( nu, omega, n = 1L, omegaIsChol = FALSE, returnChol = FALSE, type = c("invWishart", "lkj", "separation"), diagXformType = c("log", "identity", "variance", "nlmixrSqrt", "nlmixrLog", "nlmixrIdentity") )cvPost( nu, omega, n = 1L, omegaIsChol = FALSE, returnChol = FALSE, type = c("invWishart", "lkj", "separation"), diagXformType = c("log", "identity", "variance", "nlmixrSqrt", "nlmixrLog", "nlmixrIdentity") )
nu |
Degrees of Freedom (Number of Observations) for covariance matrix simulation. |
omega |
Either the estimate of covariance matrix or the estimated standard deviations in matrix form each row forming the standard deviation simulated values |
n |
Number of Matrices to sample. By default this is 1.
This is only useful when |
omegaIsChol |
is an indicator of if the omega matrix is in
the Cholesky decomposition. This is only used when |
returnChol |
Return the Cholesky decomposition of the
covariance matrix sample. This is only used when |
type |
The type of covariance posterior that is being simulated. This can be:
|
diagXformType |
Diagonal transformation type. These could be:
The nlmixr transformations only make sense when there is no off-diagonal correlations modeled. |
If your covariance matrix is a 1x1 matrix, this uses an scaled inverse chi-squared which is equivalent to the Inverse Wishart distribution in the uni-directional case.
In general, the separation strategy is preferred for diagonal
matrices. If the dimension of the matrix is below 10, lkj
is numerically faster than separation method. However, the
lkj method has densities too close to zero (XXXX) when the
dimension is above 10. In that case, though computationally more
expensive separation method performs better.
For matrices with modeled covariances, the easiest method to use is the inverse Wishart which allows the simulation of correlation matrices (XXXX). This method is more well suited for well behaved matrices, that is the variance components are not too low or too high. When modeling nonlinear mixed effects modeling matrices with too high or low variances are considered sub-optimal in describing a system. With these rules in mind, it is reasonable to use the inverse Wishart.
a matrix (n=1) or a list of matrices (n > 1)
Matthew L.Fidler & Wenping Wang
Alvarez I, Niemi J and Simpson M. (2014) Bayesian Inference for a Covariance Matrix. Conference on Applied Statistics in Agriculture.
Wang1 Z, Wu Y, and Chu H. (2018) On Equivalence of the LKJ distribution and the restricted Wishart distribution. <doi:10.48550/arXiv.1809.047463
## Sample a single covariance. draw1 <- cvPost(3, matrix(c(1, .3, .3, 1), 2, 2)) ## Sample 3 covariances set.seed(42) draw3 <- cvPost(3, matrix(c(1, .3, .3, 1), 2, 2), n = 3) ## Sample 3 covariances, but return the cholesky decomposition set.seed(42) draw3c <- cvPost(3, matrix(c(1, .3, .3, 1), 2, 2), n = 3, returnChol = TRUE) ## Sample 3 covariances with lognormal standard deviations via LKJ ## correlation sample cvPost(3, sapply(1:3, function(...) { rnorm(10) }), type = "lkj") ## or return cholesky decomposition cvPost(3, sapply(1:3, function(...) { rnorm(10) }), type = "lkj", returnChol = TRUE ) ## Sample 3 covariances with lognormal standard deviations via separation ## strategy using inverse Wishart correlation sample cvPost(3, sapply(1:3, function(...) { rnorm(10) }), type = "separation") ## or returning the cholesky decomposition cvPost(3, sapply(1:3, function(...) { rnorm(10) }), type = "separation", returnChol = TRUE )## Sample a single covariance. draw1 <- cvPost(3, matrix(c(1, .3, .3, 1), 2, 2)) ## Sample 3 covariances set.seed(42) draw3 <- cvPost(3, matrix(c(1, .3, .3, 1), 2, 2), n = 3) ## Sample 3 covariances, but return the cholesky decomposition set.seed(42) draw3c <- cvPost(3, matrix(c(1, .3, .3, 1), 2, 2), n = 3, returnChol = TRUE) ## Sample 3 covariances with lognormal standard deviations via LKJ ## correlation sample cvPost(3, sapply(1:3, function(...) { rnorm(10) }), type = "lkj") ## or return cholesky decomposition cvPost(3, sapply(1:3, function(...) { rnorm(10) }), type = "lkj", returnChol = TRUE ) ## Sample 3 covariances with lognormal standard deviations via separation ## strategy using inverse Wishart correlation sample cvPost(3, sapply(1:3, function(...) { rnorm(10) }), type = "separation") ## or returning the cholesky decomposition cvPost(3, sapply(1:3, function(...) { rnorm(10) }), type = "separation", returnChol = TRUE )
This uses simulations to match the rse
dfWishart(omega, n, rse, upper, totN = 1000, diag = TRUE, seed = 1234)dfWishart(omega, n, rse, upper, totN = 1000, diag = TRUE, seed = 1234)
omega |
represents the matrix for simulation |
n |
This represents the number of subjects/samples this comes from (used to calculate rse). When present it assumes the rse= sqrt(2)/sqrt(n) |
rse |
This is the rse that we try to match, if not specified,
it is derived from |
upper |
The upper boundary for root finding in terms of degrees of freedom. If not specified, it is n*200 |
totN |
This represents the total number of simulated inverse wishart deviates |
diag |
When |
seed |
to make the simulation reproducible, this represents the seed that is used for simulating the inverse Wishart distribution |
output from uniroot() to find the right estimate
Matthew L. Fidler
dfWishart(lotri::lotri(a+b~c(1, 0.5, 1)), 100)dfWishart(lotri::lotri(a+b~c(1, 0.5, 1)), 100)
Cumulative distribution of standard normal
phi(q)phi(q)
q |
vector of quantiles |
cumulative distribution of standard normal distribution
Matthew Fidler
# phi is equivalent to pnorm(x) phi(3) # See pnorm(3) # This is provided for NONMEM-like compatibility in rxode2 models# phi is equivalent to pnorm(x) phi(3) # See pnorm(3) # This is provided for NONMEM-like compatibility in rxode2 models
Scaled Inverse Chi Squared distribution
rinvchisq(n = 1L, nu = 1, scale = 1)rinvchisq(n = 1L, nu = 1, scale = 1)
n |
Number of random samples |
nu |
degrees of freedom of inverse chi square |
scale |
Scale of inverse chi squared distribution (default is 1). |
a vector of inverse chi squared deviates.
rinvchisq(3, 4, 1) ## Scale = 1, degrees of freedom = 4 rinvchisq(2, 4, 2) ## Scale = 2, degrees of freedom = 4rinvchisq(3, 4, 1) ## Scale = 1, degrees of freedom = 4 rinvchisq(2, 4, 2) ## Scale = 2, degrees of freedom = 4
Care should be taken with this method not to encounter the
birthday problem, described
https://www.johndcook.com/blog/2016/01/29/random-number-generator-seed-mistakes/.
Since the sitmo threefry, this currently generates
one random deviate from the uniform distribution to seed the
engine threefry and then run the code.
rxbeta(shape1, shape2, n = 1L, ncores = 1L)rxbeta(shape1, shape2, n = 1L, ncores = 1L)
shape1, shape2
|
non-negative parameters of the Beta distribution. |
n |
number of observations. If |
ncores |
Number of cores for the simulation
|
Therefore, a simple call to the random number generated followed by a second call to random number generated may have identical seeds. As the number of random number generator calls are increased the probability that the birthday problem will increase.
The key to avoid this problem is to either run all simulations in the
rxode2 environment once (therefore one seed or series of seeds
for the whole simulation), pre-generate all random variables
used for the simulation, or seed the rxode2 engine with rxSetSeed()
Internally each ID is seeded with a unique number so that the results do not depend on the number of cores used.
beta random deviates
## Use threefry engine rxbeta(0.5, 0.5, n = 10) # with rxbeta you have to explicitly state n rxbeta(5, 1, n = 10, ncores = 2) # You can parallelize the simulation using openMP rxbeta(1, 3)## Use threefry engine rxbeta(0.5, 0.5, n = 10) # with rxbeta you have to explicitly state n rxbeta(5, 1, n = 10, ncores = 2) # You can parallelize the simulation using openMP rxbeta(1, 3)
Care should be taken with this method not to encounter the
birthday problem, described
https://www.johndcook.com/blog/2016/01/29/random-number-generator-seed-mistakes/.
Since the sitmo threefry, this currently generates
one random deviate from the uniform distribution to seed the
engine threefry and then run the code.
rxbinom(size, prob, n = 1L, ncores = 1L)rxbinom(size, prob, n = 1L, ncores = 1L)
size |
number of trials (zero or more). |
prob |
probability of success on each trial. |
n |
number of observations. If |
ncores |
Number of cores for the simulation
|
Therefore, a simple call to the random number generated followed by a second call to random number generated may have identical seeds. As the number of random number generator calls are increased the probability that the birthday problem will increase.
The key to avoid this problem is to either run all simulations in the
rxode2 environment once (therefore one seed or series of seeds
for the whole simulation), pre-generate all random variables
used for the simulation, or seed the rxode2 engine with rxSetSeed()
Internally each ID is seeded with a unique number so that the results do not depend on the number of cores used.
binomial random deviates
## Use threefry engine rxbinom(10, 0.9, n = 10) # with rxbinom you have to explicitly state n rxbinom(3, 0.5, n = 10, ncores = 2) # You can parallelize the simulation using openMP rxbinom(4, 0.7)## Use threefry engine rxbinom(10, 0.9, n = 10) # with rxbinom you have to explicitly state n rxbinom(3, 0.5, n = 10, ncores = 2) # You can parallelize the simulation using openMP rxbinom(4, 0.7)
Care should be taken with this method not to encounter the
birthday problem, described
https://www.johndcook.com/blog/2016/01/29/random-number-generator-seed-mistakes/.
Since the sitmo threefry, this currently generates
one random deviate from the uniform distribution to seed the
engine threefry and then run the code.
rxcauchy(location = 0, scale = 1, n = 1L, ncores = 1L)rxcauchy(location = 0, scale = 1, n = 1L, ncores = 1L)
location, scale
|
location and scale parameters. |
n |
number of observations. If |
ncores |
Number of cores for the simulation
|
Therefore, a simple call to the random number generated followed by a second call to random number generated may have identical seeds. As the number of random number generator calls are increased the probability that the birthday problem will increase.
The key to avoid this problem is to either run all simulations in the
rxode2 environment once (therefore one seed or series of seeds
for the whole simulation), pre-generate all random variables
used for the simulation, or seed the rxode2 engine with rxSetSeed()
Internally each ID is seeded with a unique number so that the results do not depend on the number of cores used.
Cauchy random deviates
## Use threefry engine rxcauchy(0, 1, n = 10) # with rxcauchy you have to explicitly state n rxcauchy(0.5, n = 10, ncores = 2) # You can parallelize the simulation using openMP rxcauchy(3)## Use threefry engine rxcauchy(0, 1, n = 10) # with rxcauchy you have to explicitly state n rxcauchy(0.5, n = 10, ncores = 2) # You can parallelize the simulation using openMP rxcauchy(3)
Care should be taken with this method not to encounter the
birthday problem, described
https://www.johndcook.com/blog/2016/01/29/random-number-generator-seed-mistakes/.
Since the sitmo threefry, this currently generates
one random deviate from the uniform distribution to seed the
engine threefry and then run the code.
rxchisq(df, n = 1L, ncores = 1L)rxchisq(df, n = 1L, ncores = 1L)
df |
degrees of freedom (non-negative, but can be non-integer). |
n |
number of observations. If |
ncores |
Number of cores for the simulation
|
Therefore, a simple call to the random number generated followed by a second call to random number generated may have identical seeds. As the number of random number generator calls are increased the probability that the birthday problem will increase.
The key to avoid this problem is to either run all simulations in the
rxode2 environment once (therefore one seed or series of seeds
for the whole simulation), pre-generate all random variables
used for the simulation, or seed the rxode2 engine with rxSetSeed()
Internally each ID is seeded with a unique number so that the results do not depend on the number of cores used.
chi squared random deviates
## Use threefry engine rxchisq(0.5, n = 10) # with rxchisq you have to explicitly state n rxchisq(5, n = 10, ncores = 2) # You can parallelize the simulation using openMP rxchisq(1)## Use threefry engine rxchisq(0.5, n = 10) # with rxchisq you have to explicitly state n rxchisq(5, n = 10, ncores = 2) # You can parallelize the simulation using openMP rxchisq(1)
Care should be taken with this method not to encounter the
birthday problem, described
https://www.johndcook.com/blog/2016/01/29/random-number-generator-seed-mistakes/.
Since the sitmo threefry, this currently generates
one random deviate from the uniform distribution to seed the
engine threefry and then run the code.
rxexp(rate, n = 1L, ncores = 1L)rxexp(rate, n = 1L, ncores = 1L)
rate |
vector of rates. |
n |
number of observations. If |
ncores |
Number of cores for the simulation
|
Therefore, a simple call to the random number generated followed by a second call to random number generated may have identical seeds. As the number of random number generator calls are increased the probability that the birthday problem will increase.
The key to avoid this problem is to either run all simulations in the
rxode2 environment once (therefore one seed or series of seeds
for the whole simulation), pre-generate all random variables
used for the simulation, or seed the rxode2 engine with rxSetSeed()
Internally each ID is seeded with a unique number so that the results do not depend on the number of cores used.
exponential random deviates
## Use threefry engine rxexp(0.5, n = 10) # with rxexp you have to explicitly state n rxexp(5, n = 10, ncores = 2) # You can parallelize the simulation using openMP rxexp(1)## Use threefry engine rxexp(0.5, n = 10) # with rxexp you have to explicitly state n rxexp(5, n = 10, ncores = 2) # You can parallelize the simulation using openMP rxexp(1)
Care should be taken with this method not to encounter the
birthday problem, described
https://www.johndcook.com/blog/2016/01/29/random-number-generator-seed-mistakes/.
Since the sitmo threefry, this currently generates
one random deviate from the uniform distribution to seed the
engine threefry and then run the code.
rxf(df1, df2, n = 1L, ncores = 1L)rxf(df1, df2, n = 1L, ncores = 1L)
df1, df2
|
degrees of freedom. |
n |
number of observations. If |
ncores |
Number of cores for the simulation
|
Therefore, a simple call to the random number generated followed by a second call to random number generated may have identical seeds. As the number of random number generator calls are increased the probability that the birthday problem will increase.
The key to avoid this problem is to either run all simulations in the
rxode2 environment once (therefore one seed or series of seeds
for the whole simulation), pre-generate all random variables
used for the simulation, or seed the rxode2 engine with rxSetSeed()
Internally each ID is seeded with a unique number so that the results do not depend on the number of cores used.
f random deviates
## Use threefry engine rxf(0.5, 0.5, n = 10) # with rxf you have to explicitly state n rxf(5, 1, n = 10, ncores = 2) # You can parallelize the simulation using openMP rxf(1, 3)## Use threefry engine rxf(0.5, 0.5, n = 10) # with rxf you have to explicitly state n rxf(5, 1, n = 10, ncores = 2) # You can parallelize the simulation using openMP rxf(1, 3)
Care should be taken with this method not to encounter the
birthday problem, described
https://www.johndcook.com/blog/2016/01/29/random-number-generator-seed-mistakes/.
Since the sitmo threefry, this currently generates
one random deviate from the uniform distribution to seed the
engine threefry and then run the code.
rxgamma(shape, rate = 1, n = 1L, ncores = 1L)rxgamma(shape, rate = 1, n = 1L, ncores = 1L)
shape |
The shape of the gamma random variable |
rate |
an alternative way to specify the scale. |
n |
number of observations. If |
ncores |
Number of cores for the simulation
|
Therefore, a simple call to the random number generated followed by a second call to random number generated may have identical seeds. As the number of random number generator calls are increased the probability that the birthday problem will increase.
The key to avoid this problem is to either run all simulations in the
rxode2 environment once (therefore one seed or series of seeds
for the whole simulation), pre-generate all random variables
used for the simulation, or seed the rxode2 engine with rxSetSeed()
Internally each ID is seeded with a unique number so that the results do not depend on the number of cores used.
gamma random deviates
## Use threefry engine rxgamma(0.5, n = 10) # with rxgamma you have to explicitly state n rxgamma(5, n = 10, ncores = 2) # You can parallelize the simulation using openMP rxgamma(1)## Use threefry engine rxgamma(0.5, n = 10) # with rxgamma you have to explicitly state n rxgamma(5, n = 10, ncores = 2) # You can parallelize the simulation using openMP rxgamma(1)
Care should be taken with this method not to encounter the
birthday problem, described
https://www.johndcook.com/blog/2016/01/29/random-number-generator-seed-mistakes/.
Since the sitmo threefry, this currently generates
one random deviate from the uniform distribution to seed the
engine threefry and then run the code.
rxgeom(prob, n = 1L, ncores = 1L)rxgeom(prob, n = 1L, ncores = 1L)
prob |
probability of success in each trial. |
n |
number of observations. If |
ncores |
Number of cores for the simulation
|
Therefore, a simple call to the random number generated followed by a second call to random number generated may have identical seeds. As the number of random number generator calls are increased the probability that the birthday problem will increase.
The key to avoid this problem is to either run all simulations in the
rxode2 environment once (therefore one seed or series of seeds
for the whole simulation), pre-generate all random variables
used for the simulation, or seed the rxode2 engine with rxSetSeed()
Internally each ID is seeded with a unique number so that the results do not depend on the number of cores used.
geometric random deviates
## Use threefry engine rxgeom(0.5, n = 10) # with rxgeom you have to explicitly state n rxgeom(0.25, n = 10, ncores = 2) # You can parallelize the simulation using openMP rxgeom(0.75)## Use threefry engine rxgeom(0.5, n = 10) # with rxgeom you have to explicitly state n rxgeom(0.25, n = 10, ncores = 2) # You can parallelize the simulation using openMP rxgeom(0.75)
Get the rxode2 seed
rxGetSeed()rxGetSeed()
rxode2 seed state or -1 when the seed isn't set
rxSetSeed, rxWithSeed, rxWithPreserveSeed
# without setting seed rxGetSeed() # Now set the seed rxSetSeed(42) rxGetSeed() rxnorm() rxGetSeed() # don't use the rxode2 seed again rxSetSeed(-1) rxGetSeed() rxnorm() rxGetSeed()# without setting seed rxGetSeed() # Now set the seed rxSetSeed(42) rxGetSeed() rxnorm() rxGetSeed() # don't use the rxode2 seed again rxSetSeed(-1) rxGetSeed() rxnorm() rxGetSeed()
Care should be taken with this method not to encounter the
birthday problem, described
https://www.johndcook.com/blog/2016/01/29/random-number-generator-seed-mistakes/.
Since the sitmo threefry, this currently generates
one random deviate from the uniform distribution to seed the
engine threefry and then run the code.
rxnbinom(size, prob, n = 1L, ncores = 1L) rxnbinomMu(size, mu, n = 1L, ncores = 1L)rxnbinom(size, prob, n = 1L, ncores = 1L) rxnbinomMu(size, mu, n = 1L, ncores = 1L)
size |
target for number of successful trials, or dispersion parameter (the shape parameter of the gamma mixing distribution). Must be strictly positive, need not be integer. |
prob |
probability of success in each trial. |
n |
number of observations. If |
ncores |
Number of cores for the simulation
|
mu |
alternative parametrization via mean: see ‘Details’. |
Therefore, a simple call to the random number generated followed by a second call to random number generated may have identical seeds. As the number of random number generator calls are increased the probability that the birthday problem will increase.
The key to avoid this problem is to either run all simulations in the
rxode2 environment once (therefore one seed or series of seeds
for the whole simulation), pre-generate all random variables
used for the simulation, or seed the rxode2 engine with rxSetSeed()
Internally each ID is seeded with a unique number so that the results do not depend on the number of cores used.
negative binomial random deviates. Note that rxbinom2
uses the mu parameterization an the rxbinom uses the prob
parameterization (mu=size/(prob+size))
## Use threefry engine rxnbinom(10, 0.9, n = 10) # with rxbinom you have to explicitly state n rxnbinom(3, 0.5, n = 10, ncores = 2) # You can parallelize the simulation using openMP rxnbinom(4, 0.7) # use mu parameter rxnbinomMu(40, 40, n=10)## Use threefry engine rxnbinom(10, 0.9, n = 10) # with rxbinom you have to explicitly state n rxnbinom(3, 0.5, n = 10, ncores = 2) # You can parallelize the simulation using openMP rxnbinom(4, 0.7) # use mu parameter rxnbinomMu(40, 40, n=10)
Simulate random normal variable from threefry/vandercorput generator
rxnorm(mean = 0, sd = 1, n = 1L, ncores = 1L)rxnorm(mean = 0, sd = 1, n = 1L, ncores = 1L)
mean |
vector of means. |
sd |
vector of standard deviations. |
n |
number of observations |
ncores |
Number of cores for the simulation
|
normal random number deviates
## Use threefry engine rxnorm(n = 10) # with rxnorm you have to explicitly state n rxnorm(n = 10, ncores = 2) # You can parallelize the simulation using openMP rxnorm(2, 3) ## The first 2 arguments are the mean and standard deviation## Use threefry engine rxnorm(n = 10) # with rxnorm you have to explicitly state n rxnorm(n = 10, ncores = 2) # You can parallelize the simulation using openMP rxnorm(2, 3) ## The first 2 arguments are the mean and standard deviation
Get the MD5 hash of the current rxode2random revision
rxode2randomMd5()rxode2randomMd5()
md5 hash of rxode2random revision
Matthew L. Fidler
rxode2randomMd5()rxode2randomMd5()
Simulate ordinal value
rxord(...)rxord(...)
... |
the probabilities to be simulated. These should sum up to a number below one. |
The values entered into the 'rxord' simulation will simulate the probability of falling each group. If it falls outside of the specified probabilities, it will simulate the group (number of probabilities specified + 1)
A number from 1 to the (number of probabilities specified + 1)
Matthew L. Fidler
# This will give values 1, and 2 rxord(0.5) rxord(0.5) rxord(0.5) rxord(0.5) # This will give values 1, 2 and 3 rxord(0.3, 0.3) rxord(0.3, 0.3) rxord(0.3, 0.3)# This will give values 1, and 2 rxord(0.5) rxord(0.5) rxord(0.5) rxord(0.5) # This will give values 1, 2 and 3 rxord(0.3, 0.3) rxord(0.3, 0.3) rxord(0.3, 0.3)
Care should be taken with this method not to encounter the
birthday problem, described
https://www.johndcook.com/blog/2016/01/29/random-number-generator-seed-mistakes/.
Since the sitmo threefry, this currently generates
one random deviate from the uniform distribution to seed the
engine threefry and then run the code.
rxpois(lambda, n = 1L, ncores = 1L)rxpois(lambda, n = 1L, ncores = 1L)
lambda |
vector of (non-negative) means. |
n |
number of random values to return. |
ncores |
Number of cores for the simulation
|
Therefore, a simple call to the random number generated followed by a second call to random number generated may have identical seeds. As the number of random number generator calls are increased the probability that the birthday problem will increase.
The key to avoid this problem is to either run all simulations in the
rxode2 environment once (therefore one seed or series of seeds
for the whole simulation), pre-generate all random variables
used for the simulation, or seed the rxode2 engine with rxSetSeed()
Internally each ID is seeded with a unique number so that the results do not depend on the number of cores used.
poission random number deviates
## Use threefry engine rxpois(lambda = 3, n = 10) # with rxpois you have to explicitly state n rxpois(lambda = 3, n = 10, ncores = 2) # You can parallelize the simulation using openMP rxpois(4) ## The first arguments are the lambda parameter## Use threefry engine rxpois(lambda = 3, n = 10) # with rxpois you have to explicitly state n rxpois(lambda = 3, n = 10, ncores = 2) # You can parallelize the simulation using openMP rxpois(4) ## The first arguments are the lambda parameter
Simulate a from a Poisson process
rxPp( n, lambda, gamma = 1, prob = NULL, t0 = 0, tmax = Inf, randomOrder = FALSE )rxPp( n, lambda, gamma = 1, prob = NULL, t0 = 0, tmax = Inf, randomOrder = FALSE )
n |
Number of time points to simulate in the Poisson process |
lambda |
Rate of Poisson process |
gamma |
Asymmetry rate of Poisson process. When gamma=1.0, this simulates a homogenous Poisson process. When gamma<1.0, the Poisson process has more events early, when gamma > 1.0, the Poisson process has more events late in the process. When gamma is non-zero, the tmax should not be infinite but indicate the end of the Poisson process to be simulated. In most pharamcometric cases, this will be the end of the study. Internally this uses a rate of: l(t) = lambdagamma(t/tmax)^(gamma-1) |
prob |
When specified, this is a probability function with one argument, time, that gives the probability that a Poisson time t is accepted as a rejection time. |
t0 |
the starting time of the Poisson process |
tmax |
the maximum time of the Poisson process |
randomOrder |
when |
This returns a vector of the Poisson process times; If the dropout is >= tmax, then all the rest of the times are = tmax to indicate the dropout is equal to or after tmax.
Matthew Fidler
## Sample homogenous Poisson process of rate 1/10 rxPp(10, 1 / 10) ## Sample inhomogenous Poisson rate of 1/10 rxPp(10, 1 / 10, gamma = 2, tmax = 100) ## Typically the Poisson process times are in a sequential order, ## using randomOrder gives the Poisson process in random order rxPp(10, 1 / 10, gamma = 2, tmax = 10, randomOrder = TRUE) ## This uses an arbitrary function to sample a non-homogenous Poisson process rxPp(10, 1 / 10, prob = function(x) { 1 / x })## Sample homogenous Poisson process of rate 1/10 rxPp(10, 1 / 10) ## Sample inhomogenous Poisson rate of 1/10 rxPp(10, 1 / 10, gamma = 2, tmax = 100) ## Typically the Poisson process times are in a sequential order, ## using randomOrder gives the Poisson process in random order rxPp(10, 1 / 10, gamma = 2, tmax = 10, randomOrder = TRUE) ## This uses an arbitrary function to sample a non-homogenous Poisson process rxPp(10, 1 / 10, prob = function(x) { 1 / x })
This is simulated with the fast, thread-safe threefry simulator and can use multiple cores to generate the random deviates.
rxRmvn( n, mu = NULL, sigma, lower = -Inf, upper = Inf, ncores = 1, isChol = FALSE, keepNames = TRUE, a = 0.4, tol = 2.05, nlTol = 1e-10, nlMaxiter = 100L )rxRmvn( n, mu = NULL, sigma, lower = -Inf, upper = Inf, ncores = 1, isChol = FALSE, keepNames = TRUE, a = 0.4, tol = 2.05, nlTol = 1e-10, nlMaxiter = 100L )
n |
Number of random row vectors to be simulated OR the matrix to use for simulation (faster). |
mu |
mean vector |
sigma |
Covariance matrix for multivariate normal or a list
of covariance matrices. If a list of covariance matrix, each
matrix will simulate |
lower |
is a vector of the lower bound for the truncated multivariate norm |
upper |
is a vector of the upper bound for the truncated multivariate norm |
ncores |
Number of cores used in the simulation |
isChol |
A boolean indicating if |
keepNames |
Keep the names from either the mean or covariance matrix. |
a |
threshold for switching between methods; They can be tuned for maximum speed; There are three cases that are considered: case 1: a < l < u case 2: l < u < -a case 3: otherwise where l=lower and u = upper |
tol |
When case 3 is used from the above possibilities, the tol value controls the acceptance rejection and inverse-transformation; When abs(u-l)>tol, uses accept-reject from randn |
nlTol |
Tolerance for newton line-search |
nlMaxiter |
Maximum iterations for newton line-search |
If n==integer (default) the output is an (n x d) matrix
where the i-th row is the i-th simulated vector.
If is.matrix(n) then the random vector are store in n,
which is provided by the user, and the function returns
NULL invisibly.
Matthew Fidler, Zdravko Botev and some from Matteo Fasiolo
John K. Salmon, Mark A. Moraes, Ron O. Dror, and David E. Shaw (2011). Parallel Random Numbers: As Easy as 1, 2, 3. D. E. Shaw Research, New York, NY 10036, USA.
The thread safe multivariate normal was inspired from the mvnfast package by Matteo Fasiolo https://CRAN.R-project.org/package=mvnfast
The concept of the truncated multivariate normal was taken from Zdravko Botev Botev (2017) doi:10.1111/rssb.12162 and Botev and L'Ecuyer (2015) doi:10.1109/WSC.2015.7408180 and converted to thread safe simulation;
## From mvnfast ## Unlike mvnfast, uses threefry simulation d <- 5 mu <- 1:d # Creating covariance matrix tmp <- matrix(rnorm(d^2), d, d) mcov <- tcrossprod(tmp, tmp) set.seed(414) rxRmvn(4, 1:d, mcov) set.seed(414) rxRmvn(4, 1:d, mcov) set.seed(414) rxRmvn(4, 1:d, mcov, ncores = 2) # r.v. generated on the second core are different ###### Here we create the matrix that will hold the simulated # random variables upfront. A <- matrix(NA, 4, d) class(A) <- "numeric" # This is important. We need the elements of A to be of class "numeric". set.seed(414) rxRmvn(A, 1:d, mcov, ncores = 2) # This returns NULL ... A # ... but the result is here ## You can also simulate from a truncated normal: rxRmvn(10, 1:d, mcov, lower = 1:d - 1, upper = 1:d + 1) # You can also simulate from different matrices (if they match # dimensions) by using a list of matrices. matL <- lapply(1:4, function(...) { tmp <- matrix(rnorm(d^2), d, d) tcrossprod(tmp, tmp) }) rxRmvn(4, setNames(1:d, paste0("a", 1:d)), matL)## From mvnfast ## Unlike mvnfast, uses threefry simulation d <- 5 mu <- 1:d # Creating covariance matrix tmp <- matrix(rnorm(d^2), d, d) mcov <- tcrossprod(tmp, tmp) set.seed(414) rxRmvn(4, 1:d, mcov) set.seed(414) rxRmvn(4, 1:d, mcov) set.seed(414) rxRmvn(4, 1:d, mcov, ncores = 2) # r.v. generated on the second core are different ###### Here we create the matrix that will hold the simulated # random variables upfront. A <- matrix(NA, 4, d) class(A) <- "numeric" # This is important. We need the elements of A to be of class "numeric". set.seed(414) rxRmvn(A, 1:d, mcov, ncores = 2) # This returns NULL ... A # ... but the result is here ## You can also simulate from a truncated normal: rxRmvn(10, 1:d, mcov, lower = 1:d - 1, upper = 1:d + 1) # You can also simulate from different matrices (if they match # dimensions) by using a list of matrices. matL <- lapply(1:4, function(...) { tmp <- matrix(rnorm(d^2), d, d) tcrossprod(tmp, tmp) }) rxRmvn(4, setNames(1:d, paste0("a", 1:d)), matL)
This sets the seed for the rxode2 parallel random number generation. If set, then whenever a seed is set for the threefry or vandercorput simulation engine, it will use this seed, increment for the number of seeds and continue with the sequence the next time the random number generator is called.
rxSetSeed(seed)rxSetSeed(seed)
seed |
An integer that represents the rxode2 parallel and internal random number generator seed. When positive, use this seed for random number generation and increment and reseed any parallel or new engines that are being called. When negative, turn off the rxode2 seed and generate a seed from the R's uniform random number generator. Best practice is to set this seed. |
In contrast, when this is not called, the time that the vandercorput or threefry simulation engines are seeded it comes from a uniform random number generated from the standard R random seed. This may cause a duplicate seed based on the R seed state. This means that there could be correlations between simulations that do not exist This will avoid the birthday problem picking exactly the same seed using the seed state of the R random number generator. The more times the seed is called, the more likely this becomes.
Nothing, called for its side effects
Matthew Fidler
JD Cook. (2016). Random number generator seed mistakes. https://www.johndcook.com/blog/2016/01/29/random-number-generator-seed-mistakes/
rxGetSeed, rxWithSeed, rxWithPreserveSeed
rxSetSeed(42) # seed with generator 42 rxnorm() # Use R's random number generator rnorm(1) rxSetSeed(42) # reproduces the same number rxnorm() # But R's random number is not the same rnorm(1) # If we reset this to use the R's seed # (internally rxode2 uses a uniform random number to span seeds) # This can lead to duplicate sequences and seeds rxSetSeed(-1) # Now set seed works for both. # This is not recommended, but illustrates the different types of # seeds that can be generated. set.seed(42) rxnorm() rnorm(1) set.seed(42) rxnorm() rnorm(1)rxSetSeed(42) # seed with generator 42 rxnorm() # Use R's random number generator rnorm(1) rxSetSeed(42) # reproduces the same number rxnorm() # But R's random number is not the same rnorm(1) # If we reset this to use the R's seed # (internally rxode2 uses a uniform random number to span seeds) # This can lead to duplicate sequences and seeds rxSetSeed(-1) # Now set seed works for both. # This is not recommended, but illustrates the different types of # seeds that can be generated. set.seed(42) rxnorm() rnorm(1) set.seed(42) rxnorm() rnorm(1)
Care should be taken with this method not to encounter the
birthday problem, described
https://www.johndcook.com/blog/2016/01/29/random-number-generator-seed-mistakes/.
Since the sitmo threefry, this currently generates
one random deviate from the uniform distribution to seed the
engine threefry and then run the code.
rxt(df, n = 1L, ncores = 1L)rxt(df, n = 1L, ncores = 1L)
df |
degrees of freedom ( |
n |
number of observations. If |
ncores |
Number of cores for the simulation
|
Therefore, a simple call to the random number generated followed by a second call to random number generated may have identical seeds. As the number of random number generator calls are increased the probability that the birthday problem will increase.
The key to avoid this problem is to either run all simulations in the
rxode2 environment once (therefore one seed or series of seeds
for the whole simulation), pre-generate all random variables
used for the simulation, or seed the rxode2 engine with rxSetSeed()
Internally each ID is seeded with a unique number so that the results do not depend on the number of cores used.
t-distribution random numbers
## Use threefry engine rxt(df = 3, n = 10) # with rxt you have to explicitly state n rxt(df = 3, n = 10, ncores = 2) # You can parallelize the simulation using openMP rxt(4) ## The first argument is the df parameter## Use threefry engine rxt(df = 3, n = 10) # with rxt you have to explicitly state n rxt(df = 3, n = 10, ncores = 2) # You can parallelize the simulation using openMP rxt(4) ## The first argument is the df parameter
Care should be taken with this method not to encounter the
birthday problem, described
https://www.johndcook.com/blog/2016/01/29/random-number-generator-seed-mistakes/.
Since the sitmo threefry, this currently generates
one random deviate from the uniform distribution to seed the
engine threefry and then run the code.
rxunif(min = 0, max = 1, n = 1L, ncores = 1L)rxunif(min = 0, max = 1, n = 1L, ncores = 1L)
min, max
|
lower and upper limits of the distribution. Must be finite. |
n |
number of observations. If |
ncores |
Number of cores for the simulation
|
Therefore, a simple call to the random number generated followed by a second call to random number generated may have identical seeds. As the number of random number generator calls are increased the probability that the birthday problem will increase.
The key to avoid this problem is to either run all simulations in the
rxode2 environment once (therefore one seed or series of seeds
for the whole simulation), pre-generate all random variables
used for the simulation, or seed the rxode2 engine with rxSetSeed()
Internally each ID is seeded with a unique number so that the results do not depend on the number of cores used.
uniform random numbers
## Use threefry engine rxunif(min = 0, max = 4, n = 10) # with rxunif you have to explicitly state n rxunif(min = 0, max = 4, n = 10, ncores = 2) # You can parallelize the simulation using openMP rxunif()## Use threefry engine rxunif(min = 0, max = 4, n = 10) # with rxunif you have to explicitly state n rxunif(min = 0, max = 4, n = 10, ncores = 2) # You can parallelize the simulation using openMP rxunif()
Care should be taken with this method not to encounter the
birthday problem, described
https://www.johndcook.com/blog/2016/01/29/random-number-generator-seed-mistakes/.
Since the sitmo threefry, this currently generates
one random deviate from the uniform distribution to seed the
engine threefry and then run the code.
rxweibull(shape, scale = 1, n = 1L, ncores = 1L)rxweibull(shape, scale = 1, n = 1L, ncores = 1L)
shape, scale
|
shape and scale parameters, the latter defaulting to 1. |
n |
number of observations. If |
ncores |
Number of cores for the simulation
|
Therefore, a simple call to the random number generated followed by a second call to random number generated may have identical seeds. As the number of random number generator calls are increased the probability that the birthday problem will increase.
The key to avoid this problem is to either run all simulations in the
rxode2 environment once (therefore one seed or series of seeds
for the whole simulation), pre-generate all random variables
used for the simulation, or seed the rxode2 engine with rxSetSeed()
Internally each ID is seeded with a unique number so that the results do not depend on the number of cores used.
Weibull random deviates
## Use threefry engine # with rxweibull you have to explicitly state n rxweibull(shape = 1, scale = 4, n = 10) # You can parallelize the simulation using openMP rxweibull(shape = 1, scale = 4, n = 10, ncores = 2) rxweibull(3)## Use threefry engine # with rxweibull you have to explicitly state n rxweibull(shape = 1, scale = 4, n = 10) # You can parallelize the simulation using openMP rxweibull(shape = 1, scale = 4, n = 10, ncores = 2) rxweibull(3)
Preserved seed and possibly set the seed
rxWithSeed( seed, code, rxseed = rxGetSeed(), kind = "default", normal.kind = "default", sample.kind = "default" ) rxWithPreserveSeed(code)rxWithSeed( seed, code, rxseed = rxGetSeed(), kind = "default", normal.kind = "default", sample.kind = "default" ) rxWithPreserveSeed(code)
seed |
R seed to use for the session |
code |
Is the code to evaluate |
rxseed |
is the rxode2 seed that is being preserved |
kind |
character or |
normal.kind |
character string or |
sample.kind |
character string or |
returns whatever the code is returning
rxGetSeed, rxSetSeed
rxGetSeed() rxWithSeed(1, { print(rxGetSeed()) rxnorm() print(rxGetSeed()) rxnorm() }, rxseed=3)rxGetSeed() rxWithSeed(1, { print(rxGetSeed()) rxnorm() print(rxGetSeed()) rxnorm() }, rxseed=3)
Swaps the matrix list with a cube
swapMatListWithCube(matrixListOrCube)swapMatListWithCube(matrixListOrCube)
matrixListOrCube |
Either a list of 2-dimensional matrices or a cube of matrices |
A list or a cube (opposite format as input)
Matthew L. Fidler
# Create matrix list matLst <- cvPost(10, lotri::lotri(a+b~c(1, 0.25, 1)), 3) print(matLst) # Convert to cube matCube <- swapMatListWithCube(matLst) print(matCube) # Convert back to list matLst2 <- swapMatListWithCube(matCube) print(matLst2)# Create matrix list matLst <- cvPost(10, lotri::lotri(a+b~c(1, 0.25, 1)), 3) print(matLst) # Convert to cube matCube <- swapMatListWithCube(matLst) print(matCube) # Convert back to list matLst2 <- swapMatListWithCube(matCube) print(matLst2)