DRAFT! Logistic Kernel Ridge Regression Stan model. Parameters: alpha_hat = fitted coefficients, yhat2 = estimated train response; Arguments: N = number of training samples (bags), P = dimensions of Gram matrix or kernel (usually same as N), K = Gram or kernel matrix, y = response of training data (0,1), lambda = regularization coefficient in KR…

Logistic_Kernel_Ridge_Regression.stan

data {
  int<lower=0> N;
  int<lower=0> P;
  matrix[P, P] K;
  vector[N]    y;
  //real         lambda;
}
transformed data{
  // this block results verified with KRR_logit() analytical solution
  vector[N] q;
  vector[N] e_;
  vector[N] diagW;
  vector[N] pi_;
  vector[N] spec;
  vector[N] Kalpha;
  vector[N] alpha;
  //matrix[N, N] m;
  matrix[N, N] ident_N;
  ident_N = diag_matrix(rep_vector(1, rows(K)));
  //m = K + lambda * ident_N;
  for (n in 1:N) {
    alpha[n] = 1.0 / N;
  }
  Kalpha = K * alpha;
  spec = 1 + exp(-Kalpha);
  for (n in 1:N){
    pi_[n] = 1.0 / spec[n];
  }
  for (n in 1:N){
    diagW[n] = pi_[n] * (1.0 - pi_[n]);
  }
  for (n in 1:N){
    e_[n] = (y[n] - pi_[n]) / diagW[n];
  }
  for (n in 1:N){
    q[n] = Kalpha[n] + e_[n];
  }
  //print(m)
}
parameters {
  vector[N] alpha_hat;
  real<lower=0> lambda;
  real<lower=0> sigma;
  //real<lower=0> sigma2;
}
transformed parameters{
  // this block results verified with KRR_logit() analytical solution
}
model {
  // this block results approximate KRR_logit() analytical solution
  matrix[N, N] m;
  m = K + lambda * ident_N;

  q ~ normal(m * alpha_hat, sigma);
  
  alpha_hat ~ normal(0,10); // sigma is a wide guess
  lambda ~ normal(0,10); // could put a param on sigma2 here
  //sigma2 ~ cauchy(0,5);
  sigma ~ cauchy(0,5);
}
generated quantities{
  vector[N] yhat1;
  vector[N] yhat2;
  yhat1 =  (1.0 + exp(-(K * alpha_hat)));
  for (n in 1:N){
    yhat2[n] = 1.0 / yhat1[n];
  }
}


// KRR_logit <- function(K,y,lambda){
//   #### Logistic KRR
//    N = nrow(K)
//   alpha = rep(1/N, N) # initial values of alpha, transformed Parameters block
//   Kalpha = as.vector(K %*% alpha) # as 1D matrix of vector? stan wants a vector it seems
//   spec = 1 + exp(-Kalpha) # transformed Parameters block
//   pi = 1 / spec # transformed Parameters block
//   diagW = pi * (1 - pi) # transformed Parameters block
//   e = (y - pi) / diagW # transformed Parameters block // errors started here
//   q = Kalpha + e  # transformed Parameters block // errors continued here
//   ident.N <- diag(rep(1,N)) # constructed in model block
//   theSol = solve(K + lambda * ident.N, q) # the objective
//   log_pred <- 1 / (1 + exp(-as.vector(K %*% theSol))) # generated quantities
//   return(list(pred = log_pred, alphas = theSol))
// }

Edit 1) Move most of the calculations to get 'q' into the Transformed Data block. Moved the estimation of matrix 'm' into the Model block. Made 'lambda' (the Ridge coefficient) a random variable. Added parameter 'sigma'.