2

目标：我有一个已知的矩阵 A 和一个矩阵 C。C 是某个未知矩阵 BC 的函数，A 应该几乎相等。所以我尝试用 B 作为我想要找到的未知数来最小化 AC。我为此选择了 fminsearch。

1）在我的成本函数中，我收到一个向量 B 并用它来计算成本函数。计算出的成本函数是一个矩阵，我将其转换为向量并返回该向量。

2）就在调用/使用 fminsearch 之前，我将 B 定义为向量并初始化 B（定义矩阵并转换为向量）。我调用 fminsearch 并将 B 作为算法的起点。问题是不平滑和不连续的。

5 次迭代后出现错误：

Iteration   Func-count     min f(x)         Procedure

 received B 0.1 0.1 0.1 0.1
 returning  0.163476621025508 0.293733517856535 0.127416423140963 -0.00271296077015587

 received B 0.105 0.1 0.1 0.1
 returning  0.16316485573574 0.293716885638579 0.127398072233325 -0.0027135923399549

 received B 0.1 0.105 0.1 0.1
 returning  0.163460015346026 0.293389431663331 0.127415792328966 -0.00273253871794082

 received B 0.1 0.1 0.105 0.1
 returning  0.163460015346026 0.293732946124453 0.127036780812566 -0.00273253871794082

 received B 0.1 0.1 0.1 0.105
 returning  0.1634760485781 0.293715744158374 0.127396812798301 -0.00313260869222808
     0            1        0.1634766                             
     1            5        0.1631649         initial simplex     

 received B 0.1025 0.1025 0.1025 0.1025
 returning  0.163303791321532 0.293543554409175 0.127206829827164 -0.00294271094775549
Error in this$fv[ive, 1] <- fv : 
  number of items to replace is not a multiple of replacement length
In addition: Warning messages:

代码片段是：

require("neldermead")
require("fBasics")
  objective_f<-function(B){
  cat("\n received B",B)
  P_k_l_cts=matrix(0,nrow=regimes_optimum,ncol=regimes_optimum)
  B2=matrix(0,nrow=regimes_optimum,ncol=regimes_optimum)
  for(k in 1:(regimes_optimum*regimes_optimum))
  {
    m=ceiling(k/regimes_optimum)
    p=k%%regimes_optimum
    if(p==0)
    {p=regimes_optimum}
    B2[m,p]=B[k]
  }
  for(k in 1:regimes_optimum)
  {
    for(l in 1:regimes_optimum)
    {
      P_k_l_cts[k,l]=((omega%*%expm((B2-omega))*(1/lamda_cts[k]))[k,l])*(1/lamda_cts[k])
    }
  }
  P_k_l_cts_vec=vec(t(P_k_l_cts))
  P_k_l_d=matrix(0,nrow=regimes_optimum,ncol=regimes_optimum)
  for(k in 1:regimes_optimum)
  {
  for(l in 1:regimes_optimum)
  {
   if(k==l)
  {
    P_k_l_d[k,l]=(A_j_k_optimum[k,l])*(lamda_optimum[k])*(exp(-lamda_optimum[k]*    (1/lamda_optimum[k])))*(1/lamda_optimum[k])
   }
   else
   {
        P_k_l_d[k,l]=(lamda_optimum[k]/(lamda_optimum[k]-lamda_optimum[l]))*(A_j_k_optimum[k,l])*((exp(-lamda_optimum[l]*(1/lamda_optimum[k])))-(exp(-lamda_optimum[k]*(1/lamda_optimum[k]))))
      }
    }
  }

  P_k_l_d_vec=vec(t(P_k_l_d))
  cat("\n returning ",P_k_l_d_vec-P_k_l_cts_vec)
  cat("\n")
  return(P_k_l_d_vec-P_k_l_cts_vec)
}



    regime_cts_param<-function(A_j_k_cts,lamda_cts,p_cts,pi_cts,sigma_2_log_cts,regimes_optimum){      sigma_cts<<-matrix(0,nrow=regimes_optimum,ncol=1)
        for(j in 1:regimes_optimum)
     {
        sigma_cts[j]=sqrt(((1-p_cts[j])/lamda_cts[j])*exp(sigma_2_log_cts[j]))      
     }

  opt <- optimset(Display="iter",OutputFcn=outfun)
  output_fsolve=fminsearch(f=objective_f,x0=rep(0.1,regimes_optimum*regimes_optimum),opt)  
  return(list(sigma_cts=sigma_cts,B=output_fsolve$x,val=output_fsolve$fval))
}

里面的变量是全局的，我提供的是实际值：

regimes_optimum =2

sigma_2_log_optimum<-structure(c(-12.0947707925948, -12.0017466436778), .Dim = c(2L, 
1L))

lamda_optimum<-structure(c(2.1264772727783, 1.92731793847921), .Dim = c(2L, 
1L))

pi_optimum<-structure(c(0.792902540088065, 0.207097459911935), .Dim = c(2L, 
1L))

p_optimum<-structure(c(0.964696592041726, 0.393390784503164), .Dim = c(2L, 
1L))

A_j_k_optimum<-structure(c(0.613756531711541, 0.386243468288459, 0.779456839468863, 
0.220543160531137), .Dim = c(2L, 2L))

omega<-structure(c(2.1264772727783, 0, 0, 1.92731793847921), .Dim = c(2L, 
2L))

lamda_cts<-lamda_optimum
4

0 回答 0