0

我一直试图用 CVXR 弄湿我的脚，试图查看我是否可以获得已知为高斯分布且均值为零的数据的方差或精度参数的最大似然估计。

# setup the data
set.seed(23)
n = 100
M = diag(n)
x = matrix(rnorm(n))

# show an R objective function which looks convex to me when plotted
r_log_det = function(x) {
  ldet = determinant(x, logarithm = TRUE)
  as.numeric(ldet$modulus)
}

r_quad_form = function(x, y) {
  as.numeric(t(x) %*% y %*% x)
}

# Here is proof that this function optimizes to nearly the same value as the maximum likelihood estimate
test = function(p) {
  P = p*M
  S = (1/p)*M
  ldet = r_log_det(S)
  qform = r_quad_form(x, P)
  ldet + qform
}

# Plotting this for positive precision shows this function appears convex 
xi = seq(from = 0.1, to = 30, by = 0.1)
plot(xi, sapply(xi, test), type = "l")
# I also tried plotting with another "data set" where there were only two points: in that case it was still convex as well.

# Calling this mostly agrees with the result of lm for the same problem
1/sqrt(optimize(test, c(0.1, 10))$minimum) # => 0.9375864
sigma(lm(x ~ 0)) # => 0.9375892
# which is just:
sqrt(t(x) %*% x/n) # => 0.9375892

# Now, I try to implement that with CVXR
p = Variable(1, pos = TRUE)
P = diag(p, nrow = n)
S = diag(inv_pos(p), nrow = n)
loss = sum(log_det(S), quad_form(x, P))
obj  = Minimize(loss)
cons = list(p >= 0.1, p <= 10)
prob = Problem(obj, cons)
sol  = solve(prob)

但这错误表明我的函数不遵循 DCP 规则。我错过了什么？

4

0 回答 0