向量的约束优化

Question

我有一个（非对称）概率矩阵和一个观察到的整数结果向量。给定转移矩阵，我想找到一个使结果概率最大化的向量。简而言之，我试图根据它们在陆地上的最终分布来估计海上粒子的分布，以及从海洋中给定点释放的粒子最终到达陆地上给定点的概率矩阵。

我要找到的向量受制于所有分量必须在0-1之间，并且分量之和必须等于1的约束。我试图找出解决问题的最佳优化方法。

我的转换矩阵和数据集都很大，但我在这里创建了一个较小的：

我使用了一个模拟的已知海上分布 msim<-c(.3,.2,.1,.3,.1,0)和一个模拟的概率矩阵（t）来得出一个估计的海岸矩阵（Datasim2），如下：

t<-matrix (c(0,.1,.1,.1,.1,.2,0,.1,0,0,.3,0,0,0,0,.4,.1,.3,0,.1,0,.1,.4,0,0,0,.1,0,.1,.1),
nrow=5,ncol=6, byrow=T)
rownames(t)<-c("C1","C2","C3","C4","C5") ### locations on land
colnames(t)<-c("S1","S2","S3","S4","S5","S6") ### locations at sea

Datasim<-as.numeric (round((t %*% msim)*500))

Datasim2<-c(rep("C1",95), rep("C2",35), rep("C3",90),rep("C4",15),rep("C5",30))
M <-c(0.1,0.1,0.1,0.1,0.1,0.1) ## starting M

我从一个简单的函数开始，如下所示：

EstimateSource3<-function(M,Data,T){

EstEndProbsall<-M%*%T
  TotalLkhd<-rep(NA, times=dim(Data)[1])

  for (j in 1:dim(Data)[1]){

ObsEstEndLkhd<-0
    ObsEstEndLkhd<-1-EstEndProbsall[1,]  ## likelihood of particle NOT ending up at locations other than the location of interest

      IndexC<-which(colnames(EstEndProbsall)==Data$LocationCode[j], arr.ind=T) ## likelihood of ending up at location of interest

      ObsEstEndLkhd[IndexC]<-EstEndProbsall[IndexC]

      #Total likelihood
      TotalLkhd[j]<-sum(log(ObsEstEndLkhd)) 
  }

  SumTotalLkhd<-sum(TotalLkhd)


  return(SumTotalLkhd)
}

DistributionEstimate <- optim(par = M, fn = EstimateSource3, Data = Datasim2, T=t, 
control = list(fnscale = -1, trace=5, maxit=500), lower = 0, upper = 1)

为了将总和限制为 1，我尝试使用此处发布的一些建议：How to set parameters' sum to 1 in constrained optimization

例如将M<-M/sum(M) 或SumTotalLkhd<-SumTotalLkhd-(10*pwr) 添加到函数体中，但都没有产生类似msim 的任何东西，事实上，第二个解决方案出现了错误“L-BFGS-B 需要'fn'的有限值”

我认为 quadprog 包可能会有所帮助，但我认为我没有对称正定矩阵......

提前感谢您的帮助！

Answer 1

那又如何：令 D = 陆地分布，M = 海上分布，T 为转移矩阵。你知道D，T，你要计算M。你有

D' = M' T

因此 D' T' = M' (T T')

因此 D'T'(T T')^(-1) = M'

基本上你在做线性回归时解决它（似乎SO不支持数学符号：'是转置，^（-1）是普通矩阵逆。）

或者，D 可能是粒子数，现在您可以提出以下问题：海洋中最可能的粒子分布是什么。不过，这需要一种不同的方法。

Answer 2

嗯，我从来没有做过这样的模型，但按照以下思路思考。设 M 的长度为 3，D 的长度为 2，因此 T 为 3x2。我们知道 T，我们在位置 1 观察到 D_1 粒子，在位置 2 观察到 D_2 粒子。

您在位置 D_1 观察到一个粒子的可能性有多大？ Pr(D = 1) = M_1 T_11 + M_2 T_21 + M_3 T_32。类似地，Pr(D = 2) = M_1 T_12 + M_2 T_22 + M_3 T_32。现在您可以轻松编写在位置 1 和 2 观察 D_1 和 D_2 粒子的对数似然。代码可能如下所示：

loglik <- function(M) {
   if(M[1] < 0 | M[1] > 1)
      return(NA)
   if(M[2] < 0 | M[2] > 1)
      return(NA)
   M3 <- 1 - M[1] - M[2]
   if(M3 < 0 | M3 > 1)
      return(NA)
   D[1]*log(T[1,1]*M[1] + T[2,1]*M[2] + T[3,1]*M3) +
       D[2]*log(T[1,2]*M[1] + T[2,2]*M[2] + T[3,2]*M3)
}
T <- matrix(c(0.1,0.2,0.3,0.9,0.8,0.7), 3, 2)
D <- c(100,200)
library(maxLik)
m <- maxLik(loglik, start=c(0.4,0.4), method="BFGS")
summary(m)

我在估计时得到答案(0, 0.2, 0.8)，但标准误差非常大。

正如我所说，我从未这样做过，所以我不知道这是否有意义。