1 |
n_species |
|
n_years |
|
n_stations |
|
phi |
|
n_factors |
|
SpatialScale |
|
SD_O |
|
SD_E |
|
SD_extra |
|
rho |
|
logMeanDens |
|
Lmat |
|
Loc |
|
RandomSeed |
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 | ##---- Should be DIRECTLY executable !! ----
##-- ==> Define data, use random,
##-- or do help(data=index) for the standard data sets.
## The function is currently defined as
function (n_species, n_years, n_stations = 20, phi = NULL, n_factors = 2,
SpatialScale = 0.1, SD_O = 0.5, SD_E = 0.2, SD_extra = 0.1,
rho = 0.8, logMeanDens = 1, Lmat = NULL, Loc = NULL, RandomSeed = NA)
{
if (!is.na(RandomSeed))
set.seed(RandomSeed)
if (is.null(Lmat)) {
Lmat = matrix(rnorm(n_factors * n_species), nrow = n_species,
ncol = n_factors)
for (i in 1:ncol(Lmat)) {
Lmat[seq(from = 1, to = i - 1, length = i - 1), i] = 0
if (Lmat[, i][which.max(abs(Lmat[, i]))] < 0) {
Lmat[, i] = -1 * Lmat[, i]
}
}
}
if (is.null(phi))
phi = rnorm(n_factors, mean = 0, sd = 1)
Beta = rep(logMeanDens, n_species)
if (is.null(Loc))
Loc = cbind(x = runif(n_stations, min = 0, max = 1),
y = runif(n_stations, min = 0, max = 1))
model_O <- RMgauss(var = SD_O^2, scale = SpatialScale)
model_E <- RMgauss(var = SD_E^2, scale = SpatialScale)
Omega = matrix(NA, nrow = n_stations, ncol = n_factors)
for (i in 1:n_factors) {
Omega[, i] = RFsimulate(model = model_O, x = Loc[, "x"],
y = Loc[, "y"])@data[, 1]
}
Epsilon = array(NA, dim = c(n_stations, n_factors, n_years))
for (i in 1:n_factors) {
Epsilon[, i, 1] = RFsimulate(model = model_E, x = Loc[,
"x"], y = Loc[, "y"])@data[, 1]
for (t in 2:n_years) {
Epsilon[, i, t] = rho * Epsilon[, i, t - 1] + RFsimulate(model = model_E,
x = Loc[, "x"], y = Loc[, "y"])@data[, 1]
}
}
Psi = array(NA, dim = c(n_stations, n_factors, n_years))
for (i in 1:n_factors) {
for (t in 1:n_years) {
Psi[, i, t] = phi[i] * rho^t + Epsilon[, i, t] +
Omega[, i]/(1 - rho)
}
}
Theta = array(NA, dim = c(n_stations, n_species, n_years))
for (s in 1:n_stations) {
for (t in 1:n_years) {
Theta[s, , t] = Lmat %*% Psi[s, , t]
}
}
DF = NULL
for (s in 1:n_stations) {
for (p in 1:n_species) {
for (t in 1:n_years) {
Tmp = c(sitenum = s, spp = p, year = t, catch = rpois(1,
lambda = exp(Theta[s, p, t] + logMeanDens +
SD_extra * rnorm(1))), waterTmpC = 0)
DF = rbind(DF, Tmp)
}
}
}
DF = data.frame(DF, row.names = NULL)
DF[, "spp"] = factor(letters[DF[, "spp"]])
if (n_species > 26)
stop("problem with using letters")
Sim_List = list(DF = DF, Psi = Psi, Lmat = Lmat, phi = phi,
Loc = Loc, Omega = Omega, Epsilon = Epsilon, Theta = Theta,
Psi = Psi)
return(Sim_List)
}
|
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.