Во-первых, давайте создадим некоторые смоделированные данные (все специальные функции в конце ответа):
set.seed(12345) # updated to T=20 and L=40 for comparative purposes.
T = 20 # number of years
L = 40 # number of sites
N0 = 100 # average initial pop (to simulate data)
sd_env = 0.8 # to simulate the env (assumed mean 0)
env = matrix(rnorm(T*L, mean=0, sd=sd_env), nrow=T, ncol=L)
# 'real' parameters
alpha = 0.1
beta = 0.05
sd = 0.4
gamma = rnorm(T-1, mean=0, sd=sd)
mu_ini = log(rpois(n=L, lambda=N0)) #initial means
par_real = list(alpha=alpha, beta=beta, gamma=gamma,
sd=sd, mu_ini=mu_ini)
mu = dynamics(par=par_real, x=env, T=T, L=L)
# observed abundances
n = matrix(rpois(length(mu), lambda=mu), nrow=T, ncol=L)
Теперь для заданного набора параметров мы можем смоделировать ожидаемое количество особей на каждом участке и в каждом году. А поскольку у нас есть наблюдаемое количество особей, мы можем написать функцию правдоподобия для наблюдений (с распределением Пуассона) и добавить штраф за годовые отклонения скорости роста (чтобы сделать ее нормально распределенной). Для этого функция dynamics будет моделировать модель, а функция .getLogLike вычислять целевую функцию. Теперь нам нужно оптимизировать целевую функцию. Параметры для оценки: alpha, beta, годовые отклонения (gamma) и начальное ожидаемое количество особей (mu_ini) и, возможно, sigma.
Для первой попытки мы можем предоставить 0 для всех параметров в качестве начальных предположений, но для начальных ожидаемых чисел, для которых мы можем использовать начальные наблюдаемые содержания (в любом случае это MLE).
fit0 = fitModel0(obs=n, env=env, T=T, L=L)
Optimal parameters:
alpha beta gamma1 gamma2 gamma3
0.28018842 0.05464360 -0.12904373 -0.15795001 -0.04502903
gamma4 gamma5 gamma6 gamma7 gamma8
0.05045117 0.08435066 0.28864816 0.24111786 -0.80569709
gamma9 gamma10 gamma11 gamma12 gamma13
0.22786951 0.10326086 -0.50096088 -0.08880594 -0.33392310
gamma14 gamma15 gamma16 gamma17 gamma18
0.22664634 -0.47028311 0.11782381 -0.16328820 0.04208037
gamma19 mu_ini1 mu_ini2 mu_ini3 mu_ini4
0.17648808 4.14267523 4.19187205 4.05573114 3.90406443
mu_ini5 mu_ini6 mu_ini7 mu_ini8 mu_ini9
4.08975038 4.17185883 4.03679049 4.23091760 4.04940333
mu_ini10 mu_ini11 mu_ini12 mu_ini13 mu_ini14
4.19355333 4.05543081 4.15598515 4.18266682 4.09095730
mu_ini15 mu_ini16 mu_ini17 mu_ini18 mu_ini19
4.17922360 3.87211968 4.04509178 4.19385641 3.98403521
mu_ini20 mu_ini21 mu_ini22 mu_ini23 mu_ini24
4.08531659 4.19294203 4.29891769 4.21025211 4.16297457
mu_ini25 mu_ini26 mu_ini27 mu_ini28 mu_ini29
4.19265543 4.28925869 4.10752810 4.10957212 4.14953247
mu_ini30 mu_ini31 mu_ini32 mu_ini33 mu_ini34
4.09690570 4.34234547 4.18169575 4.01663339 4.32713905
mu_ini35 mu_ini36 mu_ini37 mu_ini38 mu_ini39
4.08121891 3.98256819 4.08658375 4.05942834 4.06988174
mu_ini40
4.05655031
Это работает, но обычно некоторые параметры могут быть коррелированы, и их сложнее идентифицировать из данных, поэтому можно использовать последовательный подход (можно прочитать Bolker et al. 2013 для получения дополнительной информации). В этом случае мы постепенно увеличиваем количество параметров, улучшая исходное предположение для оптимизации на каждом этапе калибровки. В этом примере на первом этапе оцениваются только alpha и beta и используются предположения для линейной модели скорости роста и окружающей среды. Затем на втором этапе мы используем оценки из первой оптимизации и добавляем годовые отклонения в качестве параметров (gamma). Наконец, мы используем оценки второй оптимизации и добавляем начальные ожидаемые значения в качестве параметров. Мы добавляем начальные ожидаемые значения последними, предполагая, что начальные наблюдаемые значения уже очень близки и являются хорошим предположением для начала, но, с другой стороны, у нас нет четкого представления о значениях остальных параметров.
fit = fitModel(obs=n, env=env, T=T, L=L)
Phase 1: alpha and beta only
Optimal parameters:
alpha beta
0.18172961 0.06323379
neg-LogLikelihood: -5023687
Phase 2: alpha, beta and gamma
Optimal parameters:
alpha beta gamma1 gamma2 gamma3
0.20519928 0.06238850 -0.35908716 -0.21453015 -0.05661066
gamma4 gamma5 gamma6 gamma7 gamma8
0.18963170 0.17800563 0.34303170 0.28960181 -0.72374927
gamma9 gamma10 gamma11 gamma12 gamma13
0.28464203 0.16900331 -0.40719047 -0.01292168 -0.25535610
gamma14 gamma15 gamma16 gamma17 gamma18
0.28806711 -0.38924648 0.19224527 -0.07875934 0.10880154
gamma19
0.24518786
neg-LogLikelihood: -5041345
Phase 3: alpha, beta, gamma and mu_ini
Optimal parameters:
alpha beta gamma1 gamma2
0.1962334008 0.0545361273 -0.4298024242 -0.1984379386
gamma3 gamma4 gamma5 gamma6
0.0240318556 0.1909639571 0.1116636126 0.3465693397
gamma7 gamma8 gamma9 gamma10
0.3478695629 -0.7500599493 0.3600551021 0.1361405398
gamma11 gamma12 gamma13 gamma14
-0.3874453347 -0.0005839263 -0.2305008546 0.2819608670
gamma15 gamma16 gamma17 gamma18
-0.3615273177 0.1792020332 -0.0685681922 0.1203006457
gamma19 mu_ini1 mu_ini2 mu_ini3
0.2506129351 4.6639314468 4.7205977429 4.5802529032
mu_ini4 mu_ini5 mu_ini6 mu_ini7
4.4293994068 4.6182382472 4.7039110982 4.5668031666
mu_ini8 mu_ini9 mu_ini10 mu_ini11
4.7610910879 4.5844180026 4.7226353021 4.5823048717
mu_ini12 mu_ini13 mu_ini14 mu_ini15
4.6814189824 4.7130039559 4.6135420745 4.7100006841
mu_ini16 mu_ini17 mu_ini18 mu_ini19
4.4080115751 4.5758092977 4.7209394881 4.5150790425
mu_ini20 mu_ini21 mu_ini22 mu_ini23
4.6171948847 4.7141188899 4.8303375556 4.7392110431
mu_ini24 mu_ini25 mu_ini26 mu_ini27
4.6893526309 4.7237687961 4.8234804183 4.6333012324
mu_ini28 mu_ini29 mu_ini30 mu_ini31
4.6392335265 4.6817044754 4.6260620666 4.8713345071
mu_ini32 mu_ini33 mu_ini34 mu_ini35
4.7107116580 4.5471434622 4.8540773708 4.6129553933
mu_ini36 mu_ini37 mu_ini38 mu_ini39
4.5134108799 4.6231016082 4.5823454113 4.5969785420
mu_ini40
4.5835763300
neg-LogLikelihood: -5047251
Сравнивая обе калибровки модели, мы видим, что вторая обеспечивает более низкое значение целевой функции. Кроме того, сравнивая корреляцию между «реальными» годовыми отклонениями и расчетными, мы имеем более высокую корреляцию для второй калибровки:
cor(gamma, fit0$par$gamma)
[1] 0.8708379
cor(gamma, fit$par$gamma)
[1] 0.9871758
И, глядя на результаты, мы видим, что у нас есть некоторые проблемы с оценкой начальных ожидаемых значений (заниженных для всех сайтов) в первой калибровке (с реальными данными обычно многофазная калибровка работает намного лучше):
par(mfrow=c(3,2), mar=c(3,5,1,1), oma=c(1,1,1,1))
for(i in 1:4) {
ylim=c(0, 1.1*log(max(fit$fitted, n)))
plot(log(fit$fitted[,i]), type="l", col="blue", ylim=ylim,
ylab="mu (log)")
lines(log(fit0$fitted[,i]), col="green")
points(log(mu[,i]), col="red")
mtext(paste("Site ", i), 3, adj=0.05, line=-2)
if(i==3) {
mtext(c("observed", "fitModel0", "fitModel1"), 1, adj=0.95,
line=-1.5:-3.5, col=c("red", "green", "blue"), cex=0.8)
}
}
mus = rbind(mu_ini, fit$par$mu_ini, fit0$par$mu_ini)
barplot(mus, beside=TRUE, col=c("red", "blue", "green"),
ylab="Initial expected population",
xlab="Sites", border=NA)
gam = rbind(gamma, fit$par$gamma, fit0$par$gamma)
barplot(gam, beside=TRUE, col=c("red", "blue", "green"),
ylab="Annual deviates", border=NA)
Окончательно,
system.time(fitModel(obs=n, env=env, T=T, L=L))
user system elapsed
9.85 0.00 9.85
Что примерно в четыре раза медленнее, чем решение, предложенное @Thierry с использованием INLA (из summary(model)):
Time used:
Pre-processing Running inla Post-processing Total
0.1070 2.3131 0.0460 2.4661
Однако после байтовой компиляции моих функций мы получаем:
user system elapsed
7.53 0.00 7.53
Это на 24% быстрее, и теперь только в 3 раза медленнее, чем метод INLA. Тем не менее, я думаю, что это разумно даже для тысяч экспериментов (моей собственной модели требуется 5 дней только для одной оптимизации, так что, возможно, у меня есть предвзятость), и поскольку мы используем смоделированные данные, мы можем сравнить надежность оценок параметров в плюс время за компьютером.
# The functions -----------------------------------------------------------
require(compiler)
dynamics = function(par, obs, x, T, L) {
alpha = par$alpha
beta = par$beta
gamma = if(!is.null((par$gamma))) par$gamma else rep(0, T-1)
mu_ini = if(!is.null(par$mu_ini)) exp(par$mu_ini) else obs[1,]
mu = matrix(nrow=T, ncol=L)
mu[1,] = mu_ini
for(t in seq_len(T-1)) {
log_mu_new = log(mu[t,]) + alpha + beta*x[t,] + gamma[t]
mu[t+1, ] = exp(log_mu_new)
}
return(mu)
}
dynamics = cmpfun(dynamics)
reListPars = function(par) {
out = list()
out$alpha = as.numeric(par["alpha"])
out$beta = as.numeric(par["beta"])
if(!is.na(par["sd"])) out$sd = as.numeric(par["sd"])
gammas = as.numeric(par[grepl("gamma", names(par))])
if(length(gammas)>0) out$gamma = gammas
mu_inis = as.numeric(par[grepl("mu_ini", names(par))])
if(length(mu_inis)>0) out$mu_ini = mu_inis
return(out)
}
reListPars = cmpfun(reListPars)
.getLogLike = function(par, obs, env, T, L) {
par = reListPars(par)
if(is.null(par$sd)) {
par$sd = if(!is.null(par$gamma)) sd(par$gamma)+0.01 else 1
}
mu = dynamics(par=par, obs=obs, x=env, T=T, L=L)
logLike = sum(obs*log(mu) - mu) - sum(par$gamma^2/(2*par$sd^2))
return(-logLike)
}
.getLogLike = cmpfun(.getLogLike)
.getUpper = function(par) {
par$alpha = 10*par$alpha + 1
par$beta = 10*abs(par$beta) + 1
if(!is.null(par$gamma)) {
if(!is.null(par$sd)) sd = par$sd else sd=sd(par$gamma)
if(sd==0) sd=1
par$gamma = rep(qnorm(0.999, sd=sd), length(par$gamma))
}
if(!is.null(par$mu_ini)) par$mu_ini = 5*par$mu_ini
if(!is.null(par$sd)) par$sd = 10*par$sd
par = unlist(par)
return(par)
}
.getUpper = cmpfun(.getUpper)
.getLower = function(par) {
par$alpha = 0 # alpha>0?
par$beta = -10*abs(par$beta) - 1
if(!is.null(par$gamma)) {
if(!is.null(par$sd)) sd = par$sd else sd=sd(par$gamma)
if(sd==0) sd=1
par$gamma = rep(qnorm(1-0.999, sd=sd), length(par$gamma))
}
if(!is.null(par$mu_ini)) par$mu_ini = 0.2*par$mu_ini
if(!is.null(par$sd)) par$sd = 0.0001*par$sd
par = unlist(par)
return(par)
}
.getLower = cmpfun(.getLower)
fitModel = function(obs, env, T, L) {
r = log(obs[-1,]/obs[-T,])
guess = data.frame(rate=as.numeric(r), env=as.numeric(env[-T,]))
mod1 = lm(rate ~ env, data=guess)
output = list()
output$par = NULL
# Phase 1: alpha an beta only
cat("Phase 1: alpha and beta only\n")
par = list()
par$alpha = as.numeric(coef(mod1)[1])
par$beta = as.numeric(coef(mod1)[2])
opt = optim(par=unlist(par), fn=.getLogLike, gr=NULL,
obs=obs, env=env, T=T, L=L, method="L-BFGS-B",
upper=.getUpper(par), lower=.getLower(par))
opt$bound = data.frame(par=unlist(par), low=.getLower(par),
upp=.getUpper(par))
output$phase1 = opt
cat("Optimal parameters: \n")
print(opt$par)
cat("\nneg-LogLikelihood: ", opt$value, "\n")
# phase 2: alpha, beta and gamma
cat("Phase 2: alpha, beta and gamma\n")
optpar = reListPars(opt$par)
par$alpha = optpar$alpha
par$beta = optpar$beta
par$gamma = rep(0, T-1)
opt = optim(par=unlist(par), fn=.getLogLike, gr=NULL,
obs=obs, env=env, T=T, L=L, method="L-BFGS-B",
upper=.getUpper(par), lower=.getLower(par))
opt$bound = data.frame(par=unlist(par), low=.getLower(par),
upp=.getUpper(par))
output$phase2 = opt
cat("Optimal parameters: \n")
print(opt$par)
cat("\nneg-LogLikelihood: ", opt$value, "\n")
# phase 3: alpha, beta, gamma and mu_ini
cat("Phase 3: alpha, beta, gamma and mu_ini\n")
optpar = reListPars(opt$par)
par$alpha = optpar$alpha
par$beta = optpar$beta
par$gamma = optpar$gamma
par$mu_ini = log(obs[1,])
opt = optim(par=unlist(par), fn=.getLogLike, gr=NULL,
obs=obs, env=env, T=T, L=L, method="L-BFGS-B",
upper=.getUpper(par), lower=.getLower(par),
control=list(maxit=1000))
opt$bound = data.frame(par=unlist(par), low=.getLower(par),
upp=.getUpper(par))
output$phase3 = opt
cat("Optimal parameters: \n")
print(opt$par)
cat("\nneg-LogLikelihood: ", opt$value, "\n")
output$par = reListPars(opt$par)
output$fitted = dynamics(par=output$par, obs=obs, x=env, T=T, L=L)
output$observed = obs
output$env = env
return(output)
}
fitModel = cmpfun(fitModel)
fitModel0 = function(obs, env, T, L) {
output = list()
output$par = NULL
par = list()
par$alpha = 0
par$beta = 0
par$gamma = rep(0, T-1)
par$mu_ini = log(obs[1,])
opt = optim(par=unlist(par), fn=.getLogLike, gr=NULL,
obs=obs, env=env, T=T, L=L, method="L-BFGS-B",
upper=.getUpper(par), lower=.getLower(par))
opt$bound = data.frame(par=unlist(par), low=.getLower(par),
upp=.getUpper(par))
output$phase1 = opt
cat("Optimal parameters: \n")
print(opt$par)
cat("\nneg-LogLikelihood: ", opt$value, "\n")
output$par = reListPars(opt$par)
output$fitted = dynamics(par=output$par, obs=obs, x=env, T=T, L=L)
output$observed = obs
output$env = env
return(output)
}
fitModel0 = cmpfun(fitModel0)
person
Ricardo Oliveros-Ramos
schedule
24.03.2014
