# R/expected.R In GUILDS: Implementation of Sampling Formulas for the Unified Neutral Model of Biodiversity and Biogeography, with or without Guild Structure

```sort_aux <- function(A) {
output <- rep(0,13)

for (k in seq_along(output)) {
start <- 2^(k - 1)
end <- -1 + 2^(k)
if(end > length(A)) {
end <- length(A)
X <- sum(A[start:end])
output[k] <- X
break
}
X <- sum(A[start:end])
output[k] <- X
}
return(output)
}

expected.SAD <- function(theta, m, J) {
if (theta < 1) {
"theta can not be below one")
}
if (m < 0) {
"m can not be below zero")
}
if (m > 1) {
"m can not be above 1")
}
if (J < 0) {
"J can not be below zero")
}

I <- (J - 1) * m / (1 - m)
}

J, n_replicates = 100) {
if (theta < 1) {
"theta can not be below one")
}
if (alpha_x < 0) {
"alpha_x can not be below zero")
}
if (alpha_x > 1) {
"alpha_x can not be above 1")
}
if (alpha_y < 0) {
"alpha_y can not be below zero")
}
if (alpha_y > 1) {
"alpha_y can not be above 1")
}
if (J < 1) {
"J can not be below one")
}

meanx <- rep(0, J)
meany <- rep(0, J)

for (r in seq_len(n_replicates)) {
M <- draw_local(theta, alpha_x, alpha_y, J)
for (m in seq_along(M\$guildX)) {
meanx[m] <- meanx[m] + M\$guildX[m]
}
for (m in seq_along(M\$guildY)) {
meany[m] <- meany[m] + M\$guildY[m]
}
}
meanx <- meanx / n_replicates
meany <- meany / n_replicates

gx <- sort_aux(meanx)
gy <- sort_aux(meany)

output <- list( guildX = gx, guildY = gy)
return(output)
}

alpha_x,
alpha_y,
Jx,
Jy,
n_replicates = 100) {
if (theta < 1) {
"theta can not be below one")
}
if (alpha_x < 0) {
"alpha_x can not be below zero")
}
if (alpha_x > 1) {
"alpha_x can not be above 1")
}
if (alpha_y < 0) {
"alpha_y can not be below zero")
}
if (alpha_y > 1) {
"alpha_y can not be above 1")
}
if (Jx < 1) {
"Jx can not be below one")
}
if (Jy < 1) {
"Jy can not be below one")
}

meanx <- rep(0,Jx)
meany <- rep(0,Jy)

for (r in seq_len(n_replicates)) {
M <- draw_local_cond(theta, alpha_x, alpha_y, Jx, Jy)
for (m in seq_along(M\$guildX)) {
meanx[m] <- meanx[m] + M\$guildX[m]
}
for (m in seq_along(M\$guildY)) {
meany[m] <- meany[m] + M\$guildY[m]
}
}
meanx <- meanx / n_replicates
meany <- meany / n_replicates

gx <- sort_aux(meanx)
gy <- sort_aux(meany)

output <- list( guildX = gx, guildY = gy)
return(output)
}
```

## Try the GUILDS package in your browser

Any scripts or data that you put into this service are public.

GUILDS documentation built on May 2, 2019, 8:25 a.m.