stotastic

// generates a pair of iid standard normals using Box-Muller
function boxMuller() {
  var x = new Array(2);
  var u1 = Math.random();
  var u2 = Math.random();
 
  var R = Math.sqrt(-2*Math.log(u1));
  var theta = 2*Math.PI*u2;
 
  x[0] = R*Math.cos(theta);
  x[1] = R*Math.sin(theta);
  return(x);
}
 
 
// generates a vector of iid standard normals
// n is the number of normals to generate
function rnorm(n) {
  var x = new Array(n);
  var y = new Array(2);
  var k = n % 2;
 
  // if n is even
  if(k==0){
    for(var i=0; i<(n/2); i++){
      y = boxMuller();
      x[i*2]     = y[0];
      x[i*2 + 1] = y[1];
    }
  }
 
  // if n is odd
  if(k==1){
    for(var i=0; i<Math.round((n-1)/2); i++){
      y = boxMuller();
      x[i*2]     = y[0];
      x[i*2 + 1] = y[1];
    }
    y = boxMuller();
    x[n-1] = y[0];
  }  
  return(x);
}
 
// returns a cholesky decomposed matrix
// A is the square matrix to decompose
function cholesky(A) {
  var n = A.length;
  // create square matrix L
  var L = new Array(n)
  for(var i=0; i<n; i++){
    L[i] = new Array(n)
    for(var j=0; j<n; j++){
      L[i][j] = 0;
    }
  }
 
  for(var i=0; i<n; i++){
    for(var j=0; j<=i; j++){
      var ss = 0;
      for(var k=0; k<j; k++){
        ss = ss + L[i][k] * L[j][k];
      }
      if(i==j){
        L[i][i] = Math.sqrt(A[i][i] - ss);
      } else{
        L[i][j] = (A[i][j] - ss) / L[j][j];
      }
    }
  }
  return(L);  
}
 
// generates multivariate normals
// n is the number of multivariate normals to generate
// mu is an array of means
// sig is an array of standard deviations
// rho is a 2 dim correlation array
function rmultinorm(n, mu, sig, rho){
  var N = mu.length;
  var A = cholesky(rho);
 
  // create x which will be the return values
  var x = new Array(n)
  for(var i=0; i<n; i++){
    x[i] = new Array(N);
    for(var j=0; j<N; j++){
      x[i][j] = 0;
    }
  }
 
  for(var k=0; k<n; k++){
    // generate standard normals
    var z = rnorm(N);
 
   // tranform to correlated standard normals
    for(var i=0; i<N; i++){
      var ss = 0;
      for(var j=0; j<N; j++){
        ss = ss + z[j]*A[i][j];
      }
      x[k][i] = ss*sig[i] + mu[i];
    }
  }
  return(x);
}