docs: private double simulateParameters(double[] curparams) {
      // Initialize alpha, beta
      for(int i = 0; i < numfit; i++) {
        Arrays.fill(alpha[i], 0.);
      }
      Arrays.fill(beta, 0.);
  
      double newchisq = 0.0;
  
      // Simulation loop over all data
      for(int di = 0; di < x.length; di++) {
        FittingFunctionResult res = func.eval(x[di], curparams);
        // compute inverse squared standard deviation of the point (confidence?)
        double sigma2inv = 1.0 / (s[di] * s[di]);
        double deltay = y[di] - res.y;
        // i2 and j2 are the indices that only count the params with dofit true!
        for(int i = 0, i2 = 0; i < numfit; i++) {
          if(dofit[i]) {
            double wt = res.gradients[i] * sigma2inv;
            // fill only half of the matrix, use symmetry below to complete the
            // remainder.
            for(int j = 0, j2 = 0; j <= i; j++) {
              if(dofit[j]) {
                alpha[i2][j2++] += wt * res.gradients[j];
              }
            }
            beta[i2++] += deltay * wt;
          }
        }
        newchisq += deltay * deltay * sigma2inv;
      }
      // fill symmetric side of matrix
      for(int i = 1; i < numfit; i++) {
        for(int j = i + 1; j < numfit; j++) {
          alpha[i][j] = alpha[j][i];
        }
      }
  
      return newchisq;
    }