NOTE: Need sn library. See
par.skewpar.meths.CI <-function(seal.mat, prey.mat, cal.mat, dist.meas, noise, nprey, R.p,
R.ps, R,p.mat, alpha, FC, ext.fa = ext.common.fa.list) {
# COMPUTES THE MIX AND SKEW-MIX INTERVALS GIVEN A SAMPLE OF PREDATOR FA
# SIGNATURES AND A PREY BASE.
I <- length(unique(prey.mat[, 1]))
ns <- nrow(seal.mat)
par.list <- vector("list", 1)
par.list.sk <- vector("list", 1)
data.star.seals <- seal.mat
for(r.p in 1:R.p) {
cat("r.p")
print(r.p)
# RESAMPLE CALIBRATION FACTORS
if(is.vector(unlist(cal.mat)) || (nrow(cal.mat) == 1) || (ncol(cal.mat) == 1)) {
data.star.cal <- matrix(cal.mat, ncol = 1)
}
else {
cal.seq <- seq(1, ncol(cal.mat), 1)
cal.sample <- sample(cal.seq, size = ncol(cal.mat),replace = T)
data.star.cal <- cal.mat[, cal.sample]
}
# RESAMPLE PREY
prey.seq <- seq(1, nrow(prey.mat), 1)
prey.sample <- tapply(prey.seq, prey.mat[, 1], sample, replace = T)
data.star.FC <- FC[unlist(prey.sample)]
data.star.prey <- prey.mat[unlist(prey.sample), ]
data.star.prey.ext <- as.data.frame(data.star.prey) [c("Species",ext.fa)]
data.star.prey.ext[, -1] <- data.star.prey.ext[, -1]/apply(data.star.prey.ext[, -1],1,sum)
p.mat.star <- as.matrix(p.prey(data.star.seals, MEANmeth(data.star.prey.ext),
as.vector(apply(data.star.cal,1, mean)), dist.meas, tapply(data.star.FC,
prey.mat[ , 1], mean, na.rm = T)))
#data.star.FC <- data.star.FC[!is.na(data.star.FC)]
# ESTIMATE NUISANCE PARAMETERS BY GENERATING R.ps PSEUDOSEALS
R.ps.mod <- floor(R.ps/ns) * ns
seq.vec <- rep(seq(1, ns, 1), floor(R.ps/ns))
p.mat.initboot <- p.mat.star[seq.vec, ]
data.initboot <- nonparsimseals.sim(data.star.prey, p.mat.initboot, data.star.cal, data.star.FC, 0,
nprey)
seals.initboot <- data.initboot$seals.pseudo
cal.initboot <- data.initboot$cal.mod
FC.initboot <- data.initboot$fat.mod
p.boot.mat <- p.prey(seals.initboot, MEANmeth(data.star.prey.ext),cal.initboot,
dist.meas,FC.initboot)
prob.zero.vec <- apply(p.boot.mat, 2, prob.zero.par.skewpar)
# FOR MIX METHOD
data.trans <- apply(p.boot.mat, 2, logistrans.nozero2.par.skewpar)
var.est <- apply(data.trans, 2, var.nozero)
# SKEW MIX METHOD
var.skest <- rep(0, I)
shape.est <- rep(0, I)
for(k in 1:I) {
if(round(prob.zero.vec[k], 6) != 1) {
par.sn <- sn.mle(y = logistrans.nozero(
p.boot.mat[, k]), plotit = F)$cp
par.sn <- cp.to.dp(param = par.sn)
var.skest[k] <- par.sn[2]
shape.est[k] <- par.sn[3]
}
} # END k
if(r.p == 1) {
par.list.prey <- list(rbind(prob.zero.vec, var.est))
par.list.sk.prey <- list(rbind(prob.zero.vec, var.skest,shape.est))
}
else {
par.list.prey <- append(par.list.prey, list(rbind(prob.zero.vec, var.est)))
par.list.sk.prey <- append(par.list.sk.prey, list(
rbind(prob.zero.vec, var.skest, shape.est)))
}
} # END R.p
par.list[[1]] <- par.list.prey
par.list.sk[[1]] <- par.list.sk.prey
cat("PAR Intervals")
CI.L.1 <- rep(NA, I)
CI.U.1 <- rep(NA, I)
CI.L.2 <- rep(NA, I)
CI.U.2 <- rep(NA, I)
alpha1 <- alpha/I
alpha2 <- alpha
for(k in 1:I) {
CI <- bisect.parmeth.2.lim(alpha1, alpha2, par.list, R, p.mat,k)
CI.L.1[k] <- CI$CI.L.1
CI.U.1[k] <- CI$CI.U.1
CI.L.2[k] <- CI$CI.L.2
CI.U.2[k] <- CI$CI.U.2
}
cat("Skew-par Intervals")
CI.L.sk.1 <- rep(NA, I)
CI.U.sk.1 <- rep(NA, I)
CI.L.sk.2 <- rep(NA, I)
CI.U.sk.2 <- rep(NA, I)
for(k in 1:I) {
CI <- bisect.parmeth.2.skew.lim.onealpha(alpha1, par.list.sk,R, p.mat, k)
CI.L.sk.1[k] <- CI$CI.L.1
CI.U.sk.1[k] <- CI$CI.U.1
CI.L.sk.2[k] <- CI$CI.L.2
CI.U.sk.2[k] <- CI$CI.U.2
}
return(CI.L.1, CI.U.1, CI.L.2, CI.U.2, CI.L.sk.1, CI.U.sk.1, CI.L.sk.2,CI.U.sk.2)
}
p.prey <- function(seal.mat, prey.mat, cal.mat, dist.meas, FC = rep(1., nrow(prey.mat)),
start.val = rep(1., nrow(prey.mat)), ext.fa = ext.common.fa.list){
# COMPUTES THE DIET ESTIMATE FOR EACH SEAL IN seal.mat USING KL (1) OR
# AIT DISTANCE MEASURE (2).
# NOTE: seal.mat MUST BE MATRIX AND prey.mat CONTAINS A
# REPRESENTATIVE FA SIGNATURE (SUCH AS THE MEAN) FROM
# EACH SPECIES.
# INPUT:
# seal.mat --> MATRIX CONTAINING FA SIGNATURES OF SEALS
# prey.mat --> MATRIX CONTAINING A REPRESENTATIVE FA SIGNATURE
FROM EACHPREY TYPE.
# *** ASSUMES THAT FIRST COLUMN CONTAINS NAME OF SPECIES ***
# cal.mat --> MATRIX OF CALIBRATION FACTORS WHERE THE iTH COLUMN IS TO
# BE USED WITH THE itH SEAL.
# dist.meas --> DISTANCE MEASURE TO USE: 1 = KL OR 2 = AITCHISON
# OUTPUT:
# p.mat --> VECTOR THAT MINIMIZED THE FUNCTION (NORMALIZED)
if(is.vector(unlist(cal.mat)) || (nrow(cal.mat) == 1.) || (ncol(cal.mat ) == 1.)) {
# ONLY ONE SEAL
# REMOVE EXTENDED DIETARY
seal.mat <- t(t(seal.mat)/as.vector(unlist(cal.mat)))
seal.mat <- as.data.frame(seal.mat)[ext.fa]
seal.mat <- seal.mat/apply(seal.mat, 1., sum)
seal.mat <- as.matrix(seal.mat)
}
else {
# REMOVE EXTENDED DIETARY
seal.mat <- seal.mat/t(cal.mat)
seal.mat <- as.data.frame(seal.mat)[ext.fa]
seal.mat <- seal.mat/apply(seal.mat, 1., sum)
}
# eps USED TO ADJUST DIET ESTIMATES THAT ARE 1.
eps <- 0.0005
I <- nrow(prey.mat)
ns <- nrow(seal.mat)
p.mat <- matrix(rep(0., nrow(prey.mat) * nrow(seal.mat)), byrow = T, nrow(seal.mat),
nrow(prey.mat))
if(!(is.matrix(start.val))) {
start.val <- matrix(rep(start.val, ns), byrow = T, ns, I)
}
if(dist.meas == 1.) {
for(i in 1.:nrow(seal.mat)) {
p.all <- nlminb(start = start.val[i, ], objective =optquantile.obj2, gradient = optquantile.grad2,
seal = seal.mat[i, ], prey.quantiles =prey.mat, lower =
rep(0.,nrow(prey.mat)),upper = rep(1., nrow(prey.mat)), control =
nlminb.control(iter.max = 1000., eval.max =1000.))
p.mat[i, ] <- p.all$parameters
p.mat[i, ] <- p.mat[i, ]/sum(p.mat[i, ])
eq.one.bool <- p.mat[i, ] == 1.
if(any(eq.one.bool)) {
p.mat[i, ][eq.one.bool] <- 1. - eps
p.mat[i, ][!eq.one.bool] <- eps/I
}
} # END for
}
else {
for(i in 1.:nrow(seal.mat)) {
p.all <- nlminb(start = start.val[i, ], objective = optcompdiff.obj,
gradient = optcompdiff.grad,seal = seal.mat[i, ],
prey.quantiles =prey.mat, lower = rep(0., nrow(prey.mat)),
upper = rep(1., nrow(prey.mat)), control =
nlminb.control(iter.max = 1000., eval.max =
1000.))
p.mat[i, ] <- p.all$parameters
p.mat[i, ] <- p.mat[i, ]/sum(p.mat[i, ])
eq.one.bool <- p.mat[i, ] == 1.
if(any(eq.one.bool)) {
p.mat[i, ][eq.one.bool] <- 1. - eps
p.mat[i, ][!eq.one.bool] <- eps/I
}
} # END for
}
FC.mat <- matrix(rep(FC, nrow(seal.mat)), byrow = T, nrow(seal.mat), I)
p.mat <- p.mat/FC.mat
p.mat <- p.mat/apply(p.mat, 1., sum)
return(p.mat)
}
optquantile.obj2 <- function(alpha, seal, prey.quantiles){
# USED IN optquantile (IN CALL TO nlminb) AS THE OBJECTIVE FUNCTION TO BE
# MINIMIZED
# INPUT:
# alpha --> VECTOR OVER WHICH MINIMIZATION TAKES PLACE.
# seal --> VECTOR OF FATTY ACID COMPOSITIONS OF SEAL.
# prey.quantiles --> MATRIX OF FATTY ACID COMPOSITION OF PREY.
# EACH ROW CONTAINS AN INDIVIDUAL PREY.
# FROM A DIFFERENT SPECIES.
no.zero <- sum(seal == 0.)
seal[seal == 0.] <- 1e-05
seal[seal > 0.] <- (1. - no.zero * 1e-05) * seal[seal > 0.]
sealhat <- t(as.matrix(alpha)) %*% prey.quantiles
no.zero <- sum(sealhat == 0.)
sealhat[sealhat == 0.] <- 1e-05
sealhat[sealhat > 0.] <- (1. - no.zero * 1e-05) * sealhat[sealhat >0.]
sealhat <- sealhat/sum(alpha)
return(KLdist(seal, sealhat))
}
optquantile.grad2 <- function(alpha, seal, prey.quantiles){
# USED IN optquantile (IN CALL TO nlminb) AS THE GRADIENT OF THE OBJECTIVE
# FUNCTION TO BE MINIMIZED
# INPUT:
# alpha --> VECTOR OVER WHICH MINIMIZATION TAKES PLACE.
# seal --> VECTOR OF FATTY ACID COMPOSITIONS OF SEAL.
# prey.quantiles --> MATRIX OF FATTY ACID COMPOSITION OF PREY. EACH ROW
# CONTAINS AN INDIVIDUAL PREYFROM A DIFFERENT SPECIES.
no.zero <- sum(seal == 0.)
seal[seal == 0.] <- 1e-05
seal[seal > 0.] <- (1. - no.zero * 1e-05) * seal[seal > 0.]
sealhat <- t(as.matrix(alpha)) %*% prey.quantiles
no.zero <- sum(sealhat == 0.)
sealhat[sealhat == 0.] <- 1e-05
sealhat[sealhat > 0.] <- (1. - no.zero * 1e-05) * sealhat[sealhat > 0.]
sumalpha <- sum(alpha)
sealhat <- sealhat/sumalpha
alphadum <- sealhat/(sumalpha)
alphadum.mat <- matrix(rep(alphadum, nrow(prey.quantiles)), byrow = T,
nrow(prey.quantiles), length(alphadum))
prey.quantilesdum <- prey.quantiles/sumalpha - alphadum.mat
Term1 <- (prey.quantilesdum) %*% t(matrix(log(seal/sealhat), ncol = length(seal)))
Term2 <- (prey.quantilesdum) %*% t(matrix((seal - sealhat)/sealhat,ncol = length(seal)))
return( - Term1 - Term2)
}
optcompdiff.obj <- function(alpha, seal, prey.quantiles){
# USED IN optquantile (IN CALL TO nlminb) AS THE OBJECTIVE FUNCTION TO BE
# MINIMIZED
# INPUT:
# alpha --> VECTOR OVER WHICH MINIMIZATION TAKES PLACE
# seal --> VECTOR OF FATTY ACID COMPOSITIONS OF SEAL
# prey.quantiles --> MATRIX OF FATTY ACID COMPOSITION OF PREY.
# EACH ROW CONTAINS AN INDIVIDUAL PREY
# FROM A DIFFERENT SPECIES.
no.zero <- sum(seal == 0.)
seal[seal == 0.] <- 1e-05
seal[seal > 0.] <- (1. - no.zero * 1e-05) * seal[seal > 0.]
sealhat <- t(as.matrix(alpha)) %*% prey.quantiles
no.zero <- sum(sealhat == 0.)
sealhat[sealhat == 0.] <- 1e-05
sealhat[sealhat > 0.] <- (1. - no.zero * 1e-05) * sealhat[sealhat > 0.]
sealhat <- sealhat/sum(alpha)
return(compdiff(seal, sealhat))
}
optcompdiff.grad <- function(alpha, seal, prey.quantiles){
# USED IN optquantile (IN CALL TO nlminb) AS THEGRADIENT FUNCTION TO BE MINIMIZED
# INPUT:
# alpha --> VECTOR OVER WHICH MINIMIZATION TAKES PLACE
# seal --> VECTOR OF FATTY ACID COMPOSITIONS OF SEAL
# prey.quantiles --> MATRIX OF FATTY ACID COMPOSITION OF PREY.
# EACH ROW CONTAINS AN INDIVIDUAL PREY
# FROM A DIFFERENT SPECIES.
nfa <- ncol(prey.quantiles)
no.zero <- sum(seal == 0.)
seal[seal == 0.] <- 1e-05
seal[seal > 0.] <- (1. - no.zero * 1e-05) * seal[seal > 0.]
sealhat <- t(as.matrix(alpha)) %*% prey.quantiles
no.zero <- sum(sealhat == 0.)
sealhat[sealhat == 0.] <- 1e-05
sealhat[sealhat > 0.] <- (1. - no.zero * 1e-05) * sealhat[sealhat > 0.]
sealhat <- sealhat/sum(alpha)
sumalpha <- sum(alpha)
alphadum <- sealhat/(sumalpha)
alphadum.mat <- matrix(rep(alphadum, nrow(prey.quantiles)), byrow = T,
nrow(prey.quantiles), length(alphadum))
prey.quantilesdum <- t(t((prey.quantiles/sumalpha - alphadum.mat))/as.vector(sealhat))
dist.dum <- log(seal/mean.geometric(seal)) - log(sealhat/mean.geometric(sealhat))
dist.dum <- matrix(dist.dum, byrow = T, 1., length(seal))
alphadum.mat2 <- matrix(rep(apply(prey.quantilesdum, 1., sum), nfa),byrow = F,
nrow(prey.quantiles), nfa)
Term1 <- prey.quantilesdum %*% t(dist.dum)
Term2 <- -1./nfa * alphadum.mat2 %*% t(dist.dum)
return(-2. * (Term1 + Term2))
}
mean.geometric <- function(x) {
# RETURNS GEOMETRIC MEAN
# INPUT:
# x --> VECTOR
D <- length(x)
return(prod(x)^(1./D))
}
MEANmeth <- function(prey.mat){
# RETURNS THE MULTIVARIATE MEAN FA SIGNATURE FROM EACH PRETY TYPE
# INPUT:
# prey.mat --> MATRIX CONTAINING FA SIGNATURES OF THE PREY. NOTE THAT
# FIRST COLUMN INDEXES PREY TYPE.
prey.means <- apply(prey.mat[, -1.], 2., tapply, prey.mat[, 1.], mean)
return(prey.means)
}
nonparsimseals.sim <- function(prey.mat, diet.null.mat, cal.mat, fat.cont, noise, nprey){
# GENERATES PSEUDO SEALS WITH iTH PSEUDO SEAL HAVING TRUE DIET
# GIVEN BY iTH ROW OF diet.null.mat
ns <- nrow(diet.null.mat)
I <- length(unique(prey.mat[, 1.]))
nFA <- ncol(prey.mat) - 1.
fat.sim <- rep(0., I)
fat.mod <- matrix(rep(0., I * ns), byrow = T, ns, I)
if(is.vector(unlist(cal.mat)) || (nrow(cal.mat) == 1.) || (ncol(cal.mat) == 1.)) {
seals.pseudo <- matrix(rep(0., ns * nFA), byrow = T, ns, nFA)
for(i in 1.:ns) {
prey.out <- split.prey(prey.mat)
prey.sim <- prey.out$prey.sim
prey.mod <- prey.out$prey.mod
for(k in 1.:I) {
fat.split <- split.fatcont(fat.cont[prey.mat[ , 1.] == unique(prey.mat[, 1.])[k]])
fat.sim[k] <- fat.split$fat.1.mean
fat.mod[i, k] <- fat.split$fat.2.mean
}
seals.pseudo[i, ] <- pseudo.seal(prey.sim, diet.null.mat[i, ], noise, nprey, cal.mat,
fat.sim, rep(F, I))
}
cal.mod <- cal.mat
seals.pseudo <- as.data.frame(seals.pseudo)
dimnames(seals.pseudo)[[2.]] <- dimnames(prey.mat[, -1.])[[2.]]
seals.pseudo <- as.matrix(seals.pseudo)
}
else {
ind.sim <- sample(seq(1., ncol(cal.mat), 1.), ns, replace = T)
cal.sim <- as.matrix(cal.mat[, ind.sim])
seals.pseudo <- matrix(rep(0., ns * nFA), byrow = T, ns, nFA)
cal.mod <- matrix(rep(0., ns * nFA), byrow = T, nFA, ns)
for(i in 1.:ns) {
prey.out <- split.prey(prey.mat)
prey.sim <- prey.out$prey.sim
prey.mod <- prey.out$prey.mod
for(k in 1.:I) {
fat.split <- split.fatcont(fat.cont[prey.mat[ , 1.] == unique(prey.mat[, 1.])[k]])
fat.sim[k] <- fat.split$fat.1.mean
fat.mod[i, k] <- fat.split$fat.2.mean
}
seals.pseudo[i, ] <- pseudo.seal(prey.sim, diet.null.mat[i, ], noise, nprey, cal.sim[
, i], fat.sim, rep(F, I))
cal.mod[, i] <- apply(cal.mat[, - ind.sim[i]], 1.,
mean)
}
seals.pseudo <- as.data.frame(seals.pseudo)
dimnames(seals.pseudo)[[2.]] <- dimnames(prey.mat[, -1.])[[
2.]]
cal.mod <- as.data.frame(cal.mod)
dimnames(cal.mod)[[1.]] <- dimnames(prey.mat[, -1.])[[2.]]
seals.pseudo <- as.matrix(seals.pseudo)
cal.mod <- as.matrix(cal.mod)
}
return(seals.pseudo, cal.mod, fat.mod, prey.mat)
}
split.prey <- function(prey.mat){
# RANDOMLY SPLITS EACH PREY TYPE INTO TWO SETS .
# IF nk<=5, prey.mod AND prey.sim WILL BE THE SAME. (ie. NOT SPLIT.)
# (RESULTS TO BE USED IN pseudo.seal.)
# INPUT:
# prey.mat --> MATRIX OF FATTY ACID SIGNATURES WHERE THE FIRST COLUMN
# DENOTES THE PREY TYPE
# OUPUT:
# prey.sim --> PREY USED TO SIMULATE SEALS (USING pseudo.seal)
# prey.mod --> MODELING PREY
numprey <- tapply(prey.mat[, 1.], prey.mat[, 1.], length)
#n.sim <- floor(numprey/2)
#n.mod <- ceiling(numprey/2)
I <- length(numprey)
n.sim <- floor(numprey/3.)
n.mod <- numprey - n.sim
split.list <- tapply(seq(1., nrow(prey.mat), 1.), prey.mat[, 1.],sample)
prey.mat.rand <- prey.mat[unlist(split.list), ]
if(numprey[1.] > 5.) {
split.bool.sim <- c(rep(T, n.sim[1.]), rep(F, n.mod[1.]))
split.bool.mod <- c(rep(F, n.sim[1.]), rep(T, n.mod[1.]))
}
else {
split.bool.sim <- c(rep(T, numprey[1.]))
split.bool.mod <- c(rep(T, numprey[1.]))
}
for(i in 2.:I) {
if(numprey[i] > 5.) {
split.bool.sim <- c(split.bool.sim, c(rep(T, n.sim[i]), rep(F, n.mod[i])))
split.bool.mod <- c(split.bool.mod, c(rep(F, n.sim[i]), rep(T, n.mod[i])))
}
else {
split.bool.sim <- c(split.bool.sim, rep(T, numprey[i]))
split.bool.mod <- c(split.bool.mod, rep(T, numprey[i]))
}
}
prey.sim <- prey.mat.rand[split.bool.sim, ]
prey.mod <- prey.mat.rand[split.bool.mod, ]
return(prey.sim, prey.mod)
}
split.fatcont <- function(fat.cont.k){
# USED TO RANDOMLY SPLIT FAT CONTENT FOR SPECIES K IN TWO
fat.cont.k <- fat.cont.k[!is.na(fat.cont.k)]
fat.sample <- sample(fat.cont.k)
fat.1 <- fat.sample[1.:round(length(fat.cont.k)/2.)]
fat.1.mean <- mean(fat.1)
fat.2 <- fat.sample[(round(length(fat.cont.k)/2.) + 1.):length(fat.cont.k)]
fat.2.mean <- mean(fat.2)
return(fat.1.mean, fat.2.mean)
}
pseudo.seal <- function(prey.sim, diet, noise, nprey, cal, fat.cont, specify.noise){
# GENERATES PSEUDO-SEALS
# THIS IS THE NEW pseudo.seal FUNCTION THAT ALLOWS 1) FAT CONTENT
# TO BE INCLUDED IN THE GENERATED SEALS AND 2) SOME SPECIES TO BE
# TRULY ZERO (THAT IS,”ZERO SPECIES” DO NOT HAVE TO BE INCLUDED IN THE
# "NOISE" )
# NOTE: IT IS ASSUMED THAT SUM(DIET) IS 1-NOISE
# INPUT:
# prey.sim --> OUTPUT OF split.prey
# diet --> DIET COMPOSITION VECTOR (NOTE: THIS VECTOR SHOULD SUM TO
# 1-NOISE. THE NOISE WILL BE ADDED TO THE diet VECTOR.)
# noise --> AMOUNT OF NOISE
# nprey --> TOTAL NUMBER OF PREY TO BE SAMPLED
# cal --> CALIBRATION FACTORS
# fat.cont --> VECTOR OF FAT CONTENT OF LENGTH=I (# OF SPECIES)
# specify.noise --> A BOOLEAN VECTOR WITH TRUES DENOTING SPECIES
# TO USE IN NOISE.
# OUTPUT:
# seal.star --> SIMULATED SEAL FA SIGNATURE.
# MODIFYING DIET TO INCLUDE FAT CONTENT
diet <- fat.cont * diet
# WANT sum(diet) + noise = 1
# MODIFYING NOISE TO INCLUDE FAT CONTENT
if(noise != 0.) {
fat.cont.mean.noise <- mean(fat.cont[specify.noise])
noise <- fat.cont.mean.noise * noise
}
diet.mod <- diet/(sum(diet) + noise)
noise <- noise/(sum(diet) + noise)
diet <- diet.mod
numprey.diet <- round(nprey * diet)
numprey.noise <- round(nprey * noise)
# FIRST SELECT FROM prey.sim, THE PREY INCLUDED IN THE DIET
I <- length(unique(prey.sim[, 1.]))
numprey.sim <- tapply(prey.sim[, 1.], prey.sim[, 1.], length)
# SAMPLE WITH REPLACEMENT numprey.diet FROM prey.sim.in
prey.sim.in <- prey.sim[rep(numprey.diet != 0., numprey.sim), ]
numprey.diet.bool <- numprey.diet != 0.
nonzero.ind <- seq(1., I, 1.)[numprey.diet.bool]
ind.list <- list(sample(seq(1., numprey.sim[nonzero.ind[1.]], 1.),numprey.diet[nonzero.ind[1.]],
replace = T))
if(sum(numprey.diet != 0.) > 1.) {
for(i in 2.:sum(numprey.diet != 0.)) {
ind.list <- append(ind.list, list(sample(seq(1.,numprey.sim[nonzero.ind[i]], i),
numprey.diet[nonzero.ind[i]], replace = T)))
}
}
numprey.sim.in <- numprey.sim[numprey.diet.bool]
ind <- unlist(ind.list)
ind <- ind + rep(cumsum(numprey.sim.in), numprey.diet[numprey.diet.bool]) –
rep(numprey.sim.in, numprey.diet[numprey.diet.bool])
# SELECT FROM prey.sim, THE PREY INCLUDED IN NOISE
prey.star <- prey.sim.in[ind, ]
if(noise != 0.) {
prey.sim.out <- prey.sim[rep(specify.noise, numprey.sim), ]
}
if((noise != 0.) & (numprey.noise != 0.)) {
sample.repl <- sample(seq(1., nrow(prey.sim.out), 1.), numprey.noise, replace = T)
if(length(sample.repl) == 1.) {
noise.star <- prey.sim.out[sample.repl, ]
seal.star <- (apply(prey.star[, -1.], 2., sum) + noise.star[1., -1.])/(nprey) * cal
}
else {
noise.star <- prey.sim.out[sample.repl, ]
seal.star <- (apply(prey.star[, -1.], 2., sum) +
apply(noise.star[, -1.], 2., sum))/(nprey) * cal
}
}
else {
seal.star <- apply(prey.star[, -1.], 2., sum)/(nprey) * cal
}
seal.star <- seal.star/sum(seal.star)
seal.star <- as.vector(seal.star)
names(seal.star) <- dimnames(prey.sim[, -1.])[[2.]]
return(seal.star)
}
prob.zero.par.skewpar <- function(x){
# RETURNS THE PROPORTION OF ELEMENTS OF VECTOR x THAT ARE ZERO.
# NOTE: IF ONLY ONE ELEMENT IS NON ZERO, CAN'T COMPUTE A VARIANCE SO JUST
# REMOVE.
x.bool <- x == 0.
if(sum(!x.bool) == 1.) {
return(1.)
}
else {
return(sum(x.bool)/length(x))
}
}
logistrans.nozero2.par.skewpar <- function(x){
# SIMILAR TO logistrans.nozero2 BUT RETURNS A VECTOR OF ZEROS
# TRANSFORMS EVERY NONZERO ELEMENT OF A VECTOR x BY THE
# LOGISTIC TRANSFORMATION
if(all(x == 0.)) {
return(x)
}
x.bool <- x == 0.
if(sum(!x.bool) == 1.) {