functions{
  real sigmoidal_neural_network(vector x, matrix w_1, vector threshold_1, vector w_2, real threshold_2){
    int M; //dimension for X
    int H;
    vector[rows(w_1)] hidden_unit;
    vector[rows(w_1)] connection_between_input_hidden;
    real connection_between_hidden_output;
    real fw;
    H <- rows(w_1);
    
    connection_between_input_hidden <- w_1*x;
    for(l in 1:H){
      hidden_unit[l] <- inv_logit(connection_between_input_hidden[l]+threshold_1[l]);
    }
    
    connection_between_hidden_output <- w_2'*hidden_unit;
    fw <- inv_logit(connection_between_hidden_output + threshold_2);
    
    return fw;
  }
}

data{
  int<lower=0> n; //number of samples
  int<lower=0> M; //dimension of input
  vector[n] y; //samples for output
  matrix[n,M] x; //samples for input
  int<lower=0> H; //number of hidden unit
}
transformed data{
  real<lower=0> lambda;
  lambda <- 0.00002;
}
parameters{
  matrix[H,M] w_1; //weights between input and hidden
  vector[H] threshold_1; // threshold for hidden unit
  
  vector[H] w_2; //weights between hidden and output
  real threshold_2; // threshold for output
}
transformed parameters{
  real<lower=0> squared_error;
  squared_error <- 0;
  for(i in 1:n){
    squared_error <- squared_error + pow(y[i]-sigmoidal_neural_network(x[i]', w_1, threshold_1, w_2, threshold_2),2);
  }
}
model{
  //for lasso
//   for(i in 1:H){
//     for(j in 1:M){
//       increment_log_prob(-lambda*fabs(w_1[i,j]));
//     }
//     increment_log_prob(-lambda*fabs(threshold_1[i]));
//     increment_log_prob(-lambda*fabs(w_2[i]));
//   }
//   increment_log_prob(-lambda*fabs(threshold_2));

  //for ridge
  for(i in 1:H){
    for(j in 1:M){
      increment_log_prob(-lambda*pow(w_1[i,j],2));
    }
    increment_log_prob(-lambda*pow(threshold_1[i],2));
    increment_log_prob(-lambda*pow(w_2[i],2));
  }
  increment_log_prob(-lambda*pow(threshold_2,2));

  increment_log_prob(-squared_error);
}