s609078902 发表于 2021-12-7 10:14 
l 和 g 是什么类型的数据
###################################################fy的估计
R=function(fy,N,x){
  corr=matrix(0,N,N)
  for(m in 1:N){
    for(n in 1:N){
      corr[m,n]=exp(-(fy[1]*(x[m,1]-x[n,1])^2+fy[2]*(x[m,2]-x[n,2])^2+fy[3]*(x[m,3]-x[n,3])^2+fy[4]*(x[m,4]-x[n,4])^2))
    }
  }
  return(corr)
}
Rh=R(fy,n1,xh)
Rl=R(fy,n2,xl)
beta_h_hat=function(fy){
  solve(t(Fh)%*%solve(R(fy,n1,xh))%*%Fh)%*%(t(Fh)%*%solve(R(fy,n1,xh))%*%yh)
}
beta_l_hat=function(fy){
  solve(t(Fl)%*%solve(R(fy,n2,xl))%*%Fl)%*%(t(Fl)%*%solve(R(fy,n2,xl))%*%yl1)
}
sigma_h_hat=function(fy){
  (1/n1)*t(yh-Fh%*%beta_h_hat(fy))%*%solve(R(fy,n1,xh))%*%(yh-Fh%*%beta_h_hat(fy))
}
sigma_l_hat=function(fy){
  (1/n2)*t(yl1-Fl%*%beta_l_hat(fy))%*%solve(R(fy,n2,xl))%*%(yl1-Fl%*%beta_l_hat(fy))
}
#fy的似然函数
fun_fy=function(fy){
  n1*log(sigma_h_hat(fy))+n2*log(sigma_l_hat(fy))+log(det(R(fy,n1,xh)))+log(det(R(fy,n2,xl)))
}
#梯度函数
fun_fy_dao=function(fy){
  deriv(fun_fy(fy))
}
fy=optim(par=c(0.5,0.2,0.4,0.6),fun_fy,method="L-BFGS-B",hessian = T)
fy
optim的第一个参数是  1,是设置的需要估计的参数fy的初始值
应该将初始值设置为一个向量的
g是需要优化的函数,最后需要返回一个向量