cov.R

#mysource("pv.R")

N=7;
K=10;
X = matrix(runif(K*N,0,1),nrow = N);
mu = colMeans(X);
X_ = X-matrix(rep(colMeans(X),N),nrow=N,byrow=TRUE);

reg = function(lambda,sig) {
  (1-lambda)*sig + lambda*diag(K)
}

sig = t(X_) %*% X_ / N;

sig_k = function(k) { t(X_[-k,]) %*% X_[-k,] / (N-1) };

vectorize = function(f) { function(X) { sapply(X,f) } }

F = function(lambda) {
  s = sum(sapply(1:N,
    function(i) {
      sig_i = reg(lambda,sig_k(i));
      s1 = X_[i,] %*% solve(sig_i) %*% X_[i,];
      s2 = log(det(sig_i));
      ## pv0(s1);
      ## pv0(s2);
      s1+s2
    }));
  -s/2
}

Fv = vectorize(F);