Nothing
lambdamu = function(n,pars,ddep)
{
lnn = length(n)
zeros = rep(0,lnn)
ones = rep(1,lnn)
la = pars[1]
mu = pars[2]
K = pars[3]
r = pars[4]
lavec = la * ones
muvec = mu * ones
n0 = (ddep == 2 | ddep == 4)
if(ddep == 1)
{
lavec = pmax(zeros,la - (la - mu) * n / K)
} else if(ddep == 1.3)
{
lavec = pmax(zeros,la * (1 - n/K))
} else if(ddep == 2 | ddep == 2.1 | ddep == 2.2)
{
y = -(log(la/mu)/log(K+n0))^(ddep != 2.2)
lavec = pmax(zeros,la * (n + n0)^y)
} else if(ddep == 2.3)
{
y = -K
lavec = pmax(zeros,la * (n + n0)^y)
} else if(ddep == 3)
{
lavec = la * ones
muvec = mu + (la - mu) * n/K
} else if(ddep == 4 | ddep == 4.1 | ddep == 4.2)
{
y = (log(la/mu)/log(K+n0))^(ddep != 4.2)
muvec = mu * (n + n0)^y
} else if(ddep == 5)
{
lavec = pmax(zeros,la - 1/(r + 1)*(la - mu)/K * n)
muvec = muvec = mu + r/(r + 1)*(la - mu)/K * n
}
return(list(lavec,muvec))
}
dd_loglik_M_aux = function(pars,lx,k,ddep)
{
nvec = 0:(lx - 1);
lambdamu_nk = lambdamu(nvec + k,pars,ddep)
lambda_nk = lambdamu_nk[[1]]
mu_nk = lambdamu_nk[[2]];
i1 = 2:lx; j1 = 1:(lx - 1); x1 = lambda_nk[1:(lx - 1)] * (nvec[1:(lx - 1)] + 2 * k)
i2 = 1:(lx - 1); j2 = 2:lx; x2 = mu_nk[2:lx] * nvec[2:lx];
i3 = 1:lx; j3 = i3; x3 = -(lambda_nk + mu_nk) * (nvec + k)
MM = Matrix::sparseMatrix(i = c(i1,i2,i3),j = c(j1,j2,j3),x = c(x1,x2,x3), dims = c(lx,lx));
MM[lx,lx] = -mu_nk[lx] * (nvec[lx] + k)
return(MM)
}
changepars = function(pars)
{
if(length(pars) <= 3)
{
sel = 1:2
} else {
sel = c(1:2,4:min(length(pars),5))
}
if(sum(pars[sel] == Inf) > 0)
{
pars[which(pars[sel] == Inf)] = 1E+10
}
if(sum(pars[sel] == 0) > 0)
{
pars[which(pars[sel] == 0)] = 1E-14
}
if(sum(pars[sel] == 1) > 0)
{
pars[which(pars[sel] == 1)] = pars[which(pars[sel] == 1)] - 1E-14
}
return(pars)
}
dd_loglik_M = function(pars,lx,k,ddep,tt,p)
{
pars = changepars(pars)
MM = dd_loglik_M_aux(pars,lx,k,ddep)
#p = expoRkit::expv(x = MM,v = p,t = tt,m = 50L)
p = expm::expAtv(A = MM,v = p,t = tt,m.max = 50)[[1]]
return(p)
}
dd_loglik_M_bw_aux = function(pars,lx,k,ddep)
{
nvec = 0:(lx - 1);
lambdamu_nk = lambdamu(nvec + k,pars,ddep)
lambda_nk = lambdamu_nk[[1]]
mu_nk = lambdamu_nk[[2]];
i1 = 2:lx; j1 = 1:(lx - 1); x1 = mu_nk[1:(lx - 1) + 1] * nvec[1:(lx - 1) + 1];
i2 = 1:(lx - 1); j2 = 2:lx; x2 = lambda_nk[2:lx - 1] * (nvec[2:lx] + 2 * k - 1);
i3 = 1:lx; j3 = i3; x3 = -(lambda_nk + mu_nk) * (nvec + k)
MM = Matrix::sparseMatrix(i = c(i1,i2,i3),j = c(j1,j2,j3),x = c(x1,x2,x3), dims = c(lx,lx));
MM[lx,lx] = -mu_nk[lx] * (nvec[lx] + k)
return(MM)
}
dd_loglik_M_bw = function(pars,lx,k,ddep,tt,p)
{
pars = changepars(pars)
MM = dd_loglik_M_bw_aux(pars,lx,k,ddep)
#p = expoRkit::expv(x = MM,v = p,t = tt,m = 50L)
p = expm::expAtv(A = MM,v = p,t = tt,m.max = 50)[[1]]
return(p)
}
lambdamu2 = function(nxt,pars,ddep)
{
lnn <- length(nxt[1,])
laM = pars[1]
muM = pars[2]
KM = pars[3]
n0 = (ddep == 2 | ddep == 4)
if(ddep == 1)
{
lavec = pmax(matrix(0,lnn,lnn),laM - (laM-muM)/KM * nxt)
muvec = muM * matrix(1,lnn,lnn)
} else if(ddep == 1.3)
{
lavec = pmax(matrix(0,lnn,lnn),laM * (1 - nxt/KM))
muvec = muM * matrix(1,lnn,lnn)
} else if(ddep == 2 | ddep == 2.1 | ddep == 2.2)
{
x = -(log(laM/muM)/log(KM+n0))^(ddep != 2.2)
lavec = pmax(matrix(0,lnn,lnn),laM * (nxt + n0)^x)
muvec = muM * matrix(1,lnn,lnn)
} else if(ddep == 2.3)
{
x = -KM
lavec = pmax(matrix(0,lnn,lnn),laM * (nxt + n0)^x)
muvec = muM * matrix(1,lnn,lnn)
} else if(ddep == 3)
{
lavec = laM * matrix(1,lnn,lnn)
muvec = muM + (laM - muM)/KM * nxt
} else if(ddep == 4 | ddep == 4.1 | ddep == 4.2)
{
lavec = laM * matrix(1,lnn,lnn)
x = (log(laM/muM)/log(KM+n0))^(ddep != 4.2)
muvec = (nxt + n0)^x
}
return(list(lavec,muvec))
}
dd_loglik_M2_aux = function(pars,lx,ddep)
{
n = 0:(lx - 1);
nx1 = rep(n,lx)
dim(nx1) = c(lx,lx) # row index = number of species in first group
nx2 = t(nx1) # column index = number of species in second group
nxt = nx1 + nx2
lambdamu_n = lambdamu2(nxt,pars,ddep)
lambda_n = lambdamu_n[[1]]
mu_n = lambdamu_n[[2]];
ly = lx^2
# set elements
auxM1 = rep(0:(lx - 1),times = lx) + rep(0:(lx - 1),each = lx)
auxM2 = auxM1; auxM2[seq(lx,ly,by = lx)] = 0;
auxM3 = rep(1,ly); auxM3[seq(lx,ly,by = lx)] = 0;
# subdiagonal
i1 = 2:ly; j1 = 1:(ly - 1); x1 = lambda_n[1 + auxM1[1:(ly - 1)]] * (auxM2[1:(ly - 1)] > 0) * rep(0:(lx - 1),lx)[-ly];
# sublxdiagonal
i2 = (lx + 1):ly; j2 = 1:(ly - lx); x2 = lambda_n[1 + auxM1[1:(ly - lx)]] * (auxM2[1:(ly - lx)] > 0) * rep(0:(lx - 2),each = lx);
# superdiagonal
i3 = 1:(ly - 1); j3 = 2:ly; x3 = mu_n[1 + auxM1[2:ly]] * rep(0:(lx - 1), times = lx)[-1];
# superlxdiagonal
i4 = 1:(ly - lx); j4 = (lx + 1):ly; x4 = (mu_n[1 + auxM1[1:ly]] * rep(0:(lx - 1),each = lx))[-c(1:lx)];
i5 = 1:ly; j5 = i5; x5 = -(lambda_n[1 + auxM1[1:ly]] + mu_n[1 + auxM1[1:ly]]) * auxM1[1:ly];
MM = Matrix::sparseMatrix(i = c(i1,i2,i3,i4,i5),j = c(j1,j2,j3,j4,j5),x = c(x1,x2,x3,x4,x5));
return(MM)
}
dd_loglik_M2 = function(pars,lx,ddep,tt,p)
{
pars = changepars(pars)
MM = dd_loglik_M2_aux(pars,lx,ddep)
#p = expoRkit::expv(x = MM,v = p,t = tt,m = 50L)
p = expm::expAtv(A = MM,v = p,t = tt,m.max = 50)[[1]]
return(p)
}
lambdamu3 = function(n,pars,ddep,kM,kS)
{
lnn = length(n) + kM + kS
laM = pars[1]
muM = pars[2]
K = pars[3]
laS = pars[4]
muS = pars[5]
n0 = (ddep == 2 | ddep == 4)
nx1 = rep(0:(lnn - 1),lnn)
dim(nx1) = c(lnn,lnn) # row index = number of species in first group
nx2 = t(nx1) # column index = number of species in second group
nxt = nx1 + nx2 + kM + kS
muvecM = muM * matrix(1,lnn,lnn)
muvecS = muS * matrix(1,lnn,lnn)
if(ddep == 1.3)
{
lavecM = pmax(matrix(0,lnn,lnn),laM * (1 - nxt/K))
lavecS = pmax(matrix(0,lnn,lnn),laS * (1 - nxt/K))
} else if(ddep == 2.3)
{
x = -K
lavecM = pmax(matrix(0,lnn,lnn),laM * (nxt + n0)^x)
lavecS = pmax(matrix(0,lnn,lnn),laS * (nxt + n0)^x)
}
return(list(lavecM,muvecM,lavecS,muvecS))
}
dd_loglik_M3_aux = function(pars,lx,ddep,kM,kS)
{
nvec = 0:(lx - 1);
lambdamu_n = lambdamu3(nvec,pars,ddep,kM,kS)
lambdaM_n = lambdamu_n[[1]]
muM_n = lambdamu_n[[2]];
lambdaS_n = lambdamu_n[[3]]
muS_n = lambdamu_n[[4]];
ly = lx^2
# set elements
auxM1 = rep(0:(lx - 1),times = lx) + rep(0:(lx - 1),each = lx)
auxM2 = auxM1; auxM2[seq(lx,ly,by = lx)] = 0;
auxM3 = rep(1,ly); auxM3[seq(lx,ly,by = lx)] = 0;
# subdiagonal
i1 = 2:ly; j1 = 1:(ly - 1); x1 = lambdaS_n[1 + auxM1[1:(ly - 1)]] * (kS + auxM2[1:(ly - 1)] > 0) * (rep(0:(lx - 1),lx)[-ly] + 2 * kS);
# sublxdiagonal
i2 = (lx + 1):ly; j2 = 1:(ly - lx); x2 = lambdaM_n[1 + auxM1[1:(ly - lx)]] * (kM + auxM2[1:(ly - lx)] > 0) * (rep(0:(lx - 2),each = lx) + 2 * kM);
# superdiagonal
i3 = 1:(ly - 1); j3 = 2:ly; x3 = muS_n[1 + auxM1[2:ly]] * rep(0:(lx - 1), times = lx)[-1];
# superlxdiagonal
i4 = 1:(ly - lx); j4 = (lx + 1):ly; x4 = (muM_n[1 + auxM1[1:ly]] * rep(0:(lx - 1),each = lx))[-c(1:lx)];
# diagonal
i5 = 1:ly; j5 = i5; x5 = -(lambdaM_n[1 + auxM1[1:ly]] + muM_n[1 + auxM1[1:ly]]) * (rep(0:(lx - 1),each = lx) + kM) - (lambdaS_n[1 + auxM1[1:ly]] + muS_n[1 + auxM1[1:ly]]) * (rep(0:(lx - 1),times = lx) + kS)
MM = Matrix::sparseMatrix(i = c(i1,i2,i3,i4,i5),j = c(j1,j2,j3,j4,j5),x = c(x1,x2,x3,x4,x5));
return(MM)
}
dd_loglik_M3 = function(pars,lx,ddep,tt,p,kM,kS)
{
pars = changepars(pars)
MM = dd_loglik_M3_aux(pars,lx,ddep,kM,kS)
#p = expoRkit::expv(x = MM,v = p,t = tt,m = 50L)
p = expm::expAtv(A = MM,v = p,t = tt,m.max = 50)[[1]]
return(p)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.