R自定义约束优化函数

问题描述

目标:使用最佳"功能估算sigma1和sigma2，而sigma2必须大于sigma1

Goal: Estimate sigma1 and sigma2 with the "optim" function, while sigma2 must be greater than sigma1

模拟数据(y)

我有以下类型的数据y:

I have the following kind of data y:

N<-50

delta<-matrix(rep(0, N*N), nrow=N, ncol=N)
 for(i in 1:(N )){
  for (j in 1:N)
    if (i == j+1 | i == j-1){
    delta[i,j] <- 1;
    }
 }

sigma1<-5
sigma2<-10
diagonal=2*sigma1^2+sigma2^2
nondiag<--sigma1^2*delta
Lambda_i<-(diag(diagonal,N)+-nondiag)/diagonal
sig<-as.matrix(diagonal*Lambda_i)
sig

mu<-rep(0, N)
y<-as.vector(mvnfast::rmvn(1,mu, sig))

创建最大似然函数

mle<-function(par){
  sigma1<-par[1]
  sigma2<-par[2]
  diagonal=2*sigma1^2+sigma2^2
  nondiag<--sigma1^2*delta
  Lambda_i<-(diag(diagonal,N)+-nondiag)/diagonal
  sig<-as.matrix(diagonal*Lambda_i)


  #lokli
  loglik<--as.numeric(mvnfast::dmvn(matrix(y, byrow=T, ncol=N),mu, sig, log=T))
  loglik
}

优化

par <- c(5,5)
fit<-optim(par,mle,hessian=T,
       method="L-BFGS-B",lower=c(0.01,0.01),
       upper=c(30,30))
fit$par

问题:如何在优化过程中设置约束:"sigma2始终大于sigma1"?

Question: How I can set the constraint: "sigma2 always greater sigma1" in the optimization procedure?

推荐答案

请继续关注我的评论.我们可以使用两个技巧:

Just to follow-up on my comment. We can use two tricks:

1. 用(sigma1 + sigma2)替换所有出现的sigma2，以确保sigma2现在代表添加到sigma1的数量.
2. 在sigma2上使用exp以确保它不是负数.

1. Replace all occurrences of sigma2 in your likelihood with (sigma1 + sigma2) to ensure that sigma2 now represents the amount that is added to sigma1
2. Use exp on sigma2 to ensure that it is non-negative.

可能性变为

mlenew<-function(par){
  sigma1<-par[1]
  sigma2<-par[2]
  diagonal=2*sigma1^2+(sigma1 + exp(sigma2))^2
  nondiag<--sigma1^2*delta
  Lambda_i<-(diag(diagonal,N)+-nondiag)/diagonal
  sig<-as.matrix(diagonal*Lambda_i)
  #lokli
  loglik<--as.numeric(mvnfast::dmvn(matrix(y, byrow=T, ncol=N),mu, sig, log=T))
  loglik
}

如果我运行您的代码，我会得到

If I run you code I get

> fit<-optim(par,mle,hessian=T,
+        method="L-BFGS-B",lower=c(0.01,0.01),
+        upper=c(30,30))
> fit$par
[1]  1.738656 12.672040

有了新代码，我得到了

> fit<-optim(par,mlenew,hessian=T,
+        method="L-BFGS-B",lower=c(0.01,0.01),
+        upper=c(30,30))
> fit$par
[1] 1.737843 2.391921

然后需要反向转换":使用新代码的sigma2旧版本的实际值为

and then you need to "back-transform": The actual value of the old version of sigma2 using the new code is

> exp(2.391921) + 1.737843
[1] 12.67232

希望这会有所帮助.

这篇关于R自定义约束优化函数的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持！