我正在通过教科书"时间序列的隐马尔可夫模型:使用 R 的简介"来自学如何在 R 中运行一些马尔可夫模型。我有点纠结于如何实施文本中提到的内容。
所以,我有以下功能:
f <- function(samples,lambda,delta) -sum(log(outer(samples,lambda,dpois)%*%delta))
我可以使用以下方法针对lambda进行优化:
optim(par, fn=f, samples=x, delta=d)
where "par" is the initial guess for lambda, for some x and d.
这工作得很好。但是,在与我尝试重现的示例相对应的文本部分中,他们指出:"对于 i=1,参数 delta 和 lambda 受 sum(delta_i(=1 的约束,...m、delta_i>0 和 lambda_i>0。因此,如果希望使用无约束的优化器(如 nlm(,则有必要重新参数化"。一种可能性是最大化相对于2m-1无约束参数的可能性"。
不受约束的参数由下式给出
eta<-log(lambda)
tau<-log(delta/(1-sum(delta))
我不完全了解如何实现这一点。如何编写函数来优化此转换后的参数空间?
当使用没有参数的optim()
时,如下所示:
simpleFun <- function(x)
(x-3)^2
out = optim(par=5,
fn=simpleFun)
参数估计集将通过out$par
获得,3
正如您所料,这种情况。 或者,您可以包装函数 在转换中f
参数如下:
out = optim(par=5,
fn=function(x)
# apply the transformation x -> x^3
simpleFun(x^3))
现在的诀窍是将正确的最佳参数集分配给您的您需要对参数应用相同转换的函数估计如下:
(out$par)^3
#> 2.99741
(是的,参数估计值略有不同。为此,人为的 例如,您可以设置method="BFGS"
以获得稍微好一点的估计值。 无论如何,这表明转型的选择确实很重要有些情况,但那是另一个讨论...
要完成答案,听起来您想使用这样的包装器
# the function to be optimized
f <- function(samples,lambda,delta)
-sum(log(outer(samples,lambda,dpois)%*%delta))
out <- optim(# par it now a 2m vector
par = c(eta1 = 1,
eta2 = 1,
eta3 = 1,
tau1 = 1,
tau2 = 1,
tau3 = 1),
# a wrapper that applies the constraints
fn=function(x,samples){
# exp() guarantees that the values of lambda are > 0
lambda = exp(x[1:3])
# delta is also > 0
delta = exp(x[4:6])
# and now it sums to 1
delta = delta / sum(delta)
f(samples,lambda,delta)
},
samples=samples)
以上保证传递给f()
的参数具有正确的约束,只要对out$par
应用相同的变换,optim()
就会估计f()
的最优参数集。