I want to add a constraint to my optimization.
Setting: Optimization of parameters b[]using a given likelihood function lf[i] <-. The values of vand are entered a1inv. As can be seen from the code below, the likelihood function uses omegavwhich is the product omegav <- s %*% t(s), and the matrix is sobtained from the parameters b[].
Problem: Using this code, I get a random b[]one that optimizes the likelihood function through omegav. However, I need those b[]that create a matrix omegavequal to the given matrix omegaoptim. After optimization, I want to have the optimal b[]one such that s %*% t(s) = omegaoptim.
I can also show my limitations as follows: I know, for example, that b[3]*b[2]+b[7]*b[6]=omegaoptim[3,2], which 0,1736304 = omegaoptim[3,2]in my case. This way I can create 25 constraints for each 5x5 matrix number omegaoptim. (see all restrictions below)
My question is how to add these additional restrictions to my optimization problem?
Function used to optimize the likelihood:
loglt <- function(b,v,a1inv) {
t <- nrow(v)
n <- ncol(v)
lf <- rep(0, t)
s55 <- -(a1inv[4,4]*b[15])/a1inv[4,5]
s <- matrix(c(b[1], 0, 0, 0, 0,
b[2], b[6], 0, 0, 0,
b[3], b[7], b[10], 0, 0,
b[4], b[8], b[11], b[13], b[15],
b[5], b[9], b[12], b[14], s55), nrow=n, byrow=T)
omegav <- s %*% t(s)
for (i in seq(t)) {
lf[i] <- -0.5*n*log(2*pi) - 0.5*log(det(omegav)) - 0.5*v[i,] %*% inv(omegav) %*% cbind(v[i,])
}
return(lf)
}
Ps: this omegaoptimcan also be defined asomegaoptim <- t(v) %*% v/nrow(v)
25 restrictions:
b[1]*b[1]==omegav[1,1]
b[1]*b[2]==omegav[1,2]
b[1]*b[3]==omegav[1,3]
b[1]*b[4]==omegav[1,4]
b[1]*b[5]==omegav[1,5]
b[2]*b[1]==omegav[2,1]
b[2]*b[2]+b[6]*b[6]==omegav[2,2]
b[2]*b[3]+b[6]*b[7]==omegav[2,3]
b[2]*b[4]+b[6]*b[8]==omegav[2,4]
b[2]*b[5]+b[6]*b[9]==omegav[2,5]
b[3]*b[1]==omegav[3,1]
b[3]*b[2]+b[7]*b[6]==omegav[3,2]
b[3]*b[3]+b[7]*b[7]+b[10]*b[10]==omegav[3,3]
b[3]*b[4]+b[7]*b[8]+b[10]*b[11]==omegav[3,4]
b[3]*b[5]+b[7]*b[9]+b[10]*b[12]==omegav[3,5]
b[4]*b[1]==omegav[4,1]
b[4]*b[2]+b[8]*b[6]==omegav[4,2]
b[4]*b[3]+b[8]*b[7]+b[11]*b[10]==omegav[4,3]
b[4]*b[4]+b[8]*b[8]+b[11]*b[11]+b[13]*b[13]+b[15]*b[15]==omegav[4,4]
b[4]*b[5]+b[8]*b[9]+b[11]*b[12]+b[13]*b[14]+b[15]*b55==omegav[4,5]
b[5]*b[1]==omegav[5,1]
b[5]*b[2]+b[9]*b[6]==omegav[5,2]
b[5]*b[3]+b[9]*b[7]+b[12]*b[10]==omegav[5,3]
b[5]*b[4]+b[9]*b[8]+b[12]*b[11]+b[14]*b[13]+s55*b[15]==omegav[5,4]
b[5]*b[5]+b[9]*b[9]+b[12]*b[12]+b[14]*b[14]+s55*s55==omegav[5,5]
Thank you so much!