PUBDEV-8947: aic glm

h2o-algos/src/main/java/hex/gam/MetricBuilderGAM.java

@@ -7,6 +7,8 @@
 import water.util.ArrayUtils;
 import water.util.MathUtils;
 
+import java.util.Arrays;
+
 import static hex.glm.GLMModel.GLMParameters.Family.*;
 
 public class MetricBuilderGAM extends ModelMetricsSupervised.MetricBuilderSupervised<MetricBuilderGAM> {
@@ -51,7 +53,7 @@ public double[] perRow(double[] ds, float[] yact, double weight, double offset,
     if (weight == 0) return ds;
     _metricBuilder.perRow(ds, yact, weight, offset, m); // grab the generic terms
     if (_glmf._family.equals(GLMModel.GLMParameters.Family.negativebinomial))
-      _log_likelihood += m.likelihood(weight, yact[0], ds[0]);
+      _log_likelihood += m.likelihood(weight, yact[0], ds);
     if (!ArrayUtils.hasNaNsOrInfs(ds) && !ArrayUtils.hasNaNsOrInfs(yact)) {
       if (_glmf._family.equals(GLMModel.GLMParameters.Family.multinomial) || _glmf._family.equals(GLMModel.GLMParameters.Family.ordinal))
         add2(yact[0], ds[0], weight, offset);
@@ -86,8 +88,10 @@ public void reduce(MetricBuilderGAM other) {
       _metricBuilder.reduce(other._metricBuilder);
     _residual_deviance += other._residual_deviance;
     _null_deviance += other._null_deviance;
-    if (_glmf._family.equals(negativebinomial))
+    if (Arrays.asList(gaussian, binomial, quasibinomial, fractionalbinomial, 
+            poisson, negativebinomial, gamma, tweedie).contains(_glmf._family)) {
       _log_likelihood += other._log_likelihood;
+    }
     _nobs += other._nobs;
     _aic2 += other._aic2;
     _wcount += other._wcount;
@@ -157,24 +161,24 @@ public ModelMetrics makeModelMetrics(Model m, Frame f, Frame adaptedFrame, Frame
       }
       mm = new ModelMetricsBinomialGLM(m, f, mm._nobs, mm._MSE, _domain, metricsBinomial._sigma, 
               metricsBinomial._auc, metricsBinomial._logloss, residualDeviance(), _null_deviance, _aic, nullDOF(), 
-              resDOF(), gl, _customMetric);
+              resDOF(), gl, _customMetric, _log_likelihood);
     } else if (_glmf._family.equals(multinomial)) {
       ModelMetricsMultinomial metricsMultinomial = (ModelMetricsMultinomial) mm;
       mm = new ModelMetricsBinomialGLM.ModelMetricsMultinomialGLM(m, f, metricsMultinomial._nobs, 
               metricsMultinomial._MSE, metricsMultinomial._domain, metricsMultinomial._sigma, metricsMultinomial._cm, 
               metricsMultinomial._hit_ratios, metricsMultinomial._logloss, residualDeviance(),_null_deviance, _aic, 
-              nullDOF(), resDOF(), metricsMultinomial._auc,  _customMetric);
+              nullDOF(), resDOF(), metricsMultinomial._auc,  _customMetric, _log_likelihood);
     } else if (_glmf._family == GLMModel.GLMParameters.Family.ordinal) { // ordinal should have a different resDOF()
       ModelMetricsOrdinal metricsOrdinal = (ModelMetricsOrdinal) mm;
       mm = new ModelMetricsBinomialGLM.ModelMetricsOrdinalGLM(m, f, metricsOrdinal._nobs, metricsOrdinal._MSE, 
               metricsOrdinal._domain, metricsOrdinal._sigma, metricsOrdinal._cm, metricsOrdinal._hit_ratios, 
-              metricsOrdinal._logloss, residualDeviance(), _null_deviance, _aic, nullDOF(), resDOF(), _customMetric);
+              metricsOrdinal._logloss, residualDeviance(), _null_deviance, _aic, nullDOF(), resDOF(), _customMetric, _log_likelihood);
     } else {
       ModelMetricsRegression metricsRegression = (ModelMetricsRegression) mm;
       mm = new ModelMetricsRegressionGLM(m, f, metricsRegression._nobs, metricsRegression._MSE, 
               metricsRegression._sigma, metricsRegression._mean_absolute_error, 
               metricsRegression._root_mean_squared_log_error, residualDeviance(), 
-              residualDeviance() / _wcount, _null_deviance, _aic, nullDOF(), resDOF(), _customMetric);
+              residualDeviance() / _wcount, _null_deviance, _aic, nullDOF(), resDOF(), _customMetric, _log_likelihood);
     }
     return gamM.addModelMetrics(mm);
   }

h2o-algos/src/main/java/hex/generic/GenericModelOutput.java

@@ -102,7 +102,7 @@ auc, binomial._logloss, convertTable(binomial._gains_lift_table),
                             convertTable(binomial._thresholds_and_metric_scores), convertTable(binomial._max_criteria_and_metric_scores),
                             convertTable(binomial._confusion_matrix), glmBinomial._nullDegreesOfFreedom, glmBinomial._residualDegreesOfFreedom,
                             glmBinomial._resDev, glmBinomial._nullDev, glmBinomial._AIC, convertTable(modelAttributesGLM._coefficients_table),
-                            glmBinomial._r2, glmBinomial._description);
+                            glmBinomial._r2, glmBinomial._description, glmBinomial._loglikelihood);
                 } else {
                     return new ModelMetricsBinomialGeneric(null, null, mojoMetrics._nobs, mojoMetrics._MSE,
                             _domains[_domains.length - 1], binomial._sigma,
@@ -126,7 +126,7 @@ glmMultinomial._logloss, customMetric(mojoMetrics),
                             glmMultinomial._mean_per_class_error, glmMultinomial._nullDegreesOfFreedom, glmMultinomial._residualDegreesOfFreedom,
                             glmMultinomial._resDev, glmMultinomial._nullDev, glmMultinomial._AIC, convertTable(modelAttributesGLM._coefficients_table),
                             glmMultinomial._r2, convertTable(glmMultinomial._multinomial_auc), convertTable(glmMultinomial._multinomial_aucpr),
-                            MultinomialAucType.valueOf((String)modelAttributes.getParameterValueByName("auc_type")), glmMultinomial._description);
+                            MultinomialAucType.valueOf((String)modelAttributes.getParameterValueByName("auc_type")), glmMultinomial._description, glmMultinomial._loglikelihood);
                 } else {
                     final MojoModelMetricsMultinomial multinomial = (MojoModelMetricsMultinomial) mojoMetrics;
                     return new ModelMetricsMultinomialGeneric(null, null, mojoMetrics._nobs, mojoMetrics._MSE,
@@ -147,7 +147,7 @@ multinomial._mean_per_class_error, multinomial._r2, convertTable(multinomial._mu
                             regressionGLM._sigma, regressionGLM._mae, regressionGLM._root_mean_squared_log_error, regressionGLM._mean_residual_deviance,
                             customMetric(regressionGLM), regressionGLM._r2,
                             regressionGLM._nullDegreesOfFreedom, regressionGLM._residualDegreesOfFreedom, regressionGLM._resDev,
-                            regressionGLM._nullDev, regressionGLM._AIC, convertTable(modelAttributesGLM._coefficients_table));
+                            regressionGLM._nullDev, regressionGLM._AIC, regressionGLM._loglikelihood, convertTable(modelAttributesGLM._coefficients_table));
                 } else {
                     MojoModelMetricsRegression metricsRegression = (MojoModelMetricsRegression) mojoMetrics;
 
@@ -174,7 +174,7 @@ multinomial._mean_per_class_error, multinomial._r2, convertTable(multinomial._mu
                             ordinalMetrics._domain, ordinalMetrics._sigma, convertTable(ordinalMetrics._cm), ordinalMetrics._hit_ratios,
                             ordinalMetrics._logloss, customMetric(ordinalMetrics),
                             ordinalMetrics._r2, ordinalMetrics._nullDegreesOfFreedom, ordinalMetrics._residualDegreesOfFreedom, ordinalMetrics._resDev,
-                            ordinalMetrics._nullDev, ordinalMetrics._AIC, convertTable(modelAttributesGLM._coefficients_table),
+                            ordinalMetrics._nullDev, ordinalMetrics._AIC, ordinalMetrics._loglikelihood, convertTable(modelAttributesGLM._coefficients_table),
                             convertTable(ordinalMetrics._hit_ratio_table), ordinalMetrics._mean_per_class_error, ordinalMetrics._description);
                 } else {
                     MojoModelMetricsOrdinal ordinalMetrics = (MojoModelMetricsOrdinal) mojoMetrics;

h2o-algos/src/main/java/hex/glm/DispersionUtils.java

@@ -264,15 +264,13 @@ public void reduce(CalculateInitialTheta mrt) {
         }
     };
 
-    public static double estimateNegBinomialDispersionMomentMethod(GLMModel model, double[] beta, DataInfo dinfo, Vec weights, Vec response) {
-        DispersionTask.GenPrediction gPred = new DispersionTask.GenPrediction(beta, model, dinfo).doAll(
-                1, Vec.T_NUM, dinfo._adaptedFrame);
-        Vec mu = gPred.outputFrame(Key.make(), new String[]{"prediction"}, null).vec(0);
+    public static double estimateNegBinomialDispersionMomentMethod(GLMModel model, double[] beta, DataInfo dinfo, Vec weights, Vec response, Vec mu) {
         class MomentMethodThetaEstimation extends MRTask<MomentMethodThetaEstimation> {
             double _muSqSum;
             double _sSqSum;
             double _muSum;
             double _wSum;
+
             @Override
             public void map(Chunk[] cs) {
                 // mu, y, w
@@ -292,10 +290,11 @@ public void reduce(MomentMethodThetaEstimation mrt) {
                 _muSum += mrt._muSum;
                 _wSum += mrt._wSum;
             }
-        };
+        }
+        ;
         MomentMethodThetaEstimation mm = new MomentMethodThetaEstimation().doAll(mu, response, weights);
 
-        return mm._muSqSum/(mm._sSqSum - mm._muSum/mm._wSum);
+        return mm._muSqSum / (mm._sSqSum - mm._muSum / mm._wSum);
     }
 
 

h2o-algos/src/main/java/hex/glm/GLM.java

@@ -1281,6 +1281,29 @@ public void init(boolean expensive) {
         if (_parms._fix_tweedie_variance_power && !_parms._fix_dispersion_parameter)
           _tweedieDispersionOnly = true;
 
+        // likelihood calculation for gaussian, gamma, negativebinomial and tweedie families requires dispersion parameter estimation
+        // _dispersion_parameter_method: gaussian - pearson (default); gamma, negativebinomial, tweedie - ml.
+        if(!_parms._HGLM && _parms._calc_like) {
+          switch (_parms._family) {
+            case gaussian:
+              _parms._compute_p_values = true;
+              _parms._remove_collinear_columns = true;
+              break;
+            case gamma:
+            case negativebinomial:
+              _parms._compute_p_values = true;
+              _parms._remove_collinear_columns = true;
+            case tweedie:
+              // dispersion value estimation for tweedie family does not require 
+              // parameters compute_p_values and remove_collinear_columns
+              _parms._dispersion_parameter_method = ml;
+              // disable regularization as ML is supported only without regularization
+              _parms._lambda = new double[] {0.0};
+            default:
+              // other families does not require dispersion parameter estimation
+          }
+        }
+
       if (_parms.hasCheckpoint()) {
         if (!Family.gaussian.equals(_parms._family))  // Gaussian it not iterative and therefore don't care
           _checkPointFirstIter = true;  // mark the first iteration during iteration process of training
@@ -2446,62 +2469,67 @@ private boolean updateNegativeBinomialDispersion(int iterCnt, double[] betaCnd,
       double delta;
       double theta;
       boolean converged = false;
-      if (iterCnt == 1) {
-        theta = estimateNegBinomialDispersionMomentMethod(_model, betaCnd, _dinfo, weights, response);
-      } else {
-        theta = _parms._theta;
+      try {
+        Scope.enter();
         DispersionTask.GenPrediction gPred = new DispersionTask.GenPrediction(betaCnd, _model, _dinfo).doAll(
                 1, Vec.T_NUM, _dinfo._adaptedFrame);
-        Vec mu = gPred.outputFrame(Key.make(), new String[]{"prediction"}, null).vec(0);
-
-        NegativeBinomialGradientAndHessian nbGrad = new NegativeBinomialGradientAndHessian(theta).doAll(mu, response, weights);
-        delta = _parms._dispersion_learning_rate * nbGrad._grad / nbGrad._hess;
-        double bestLLH = Math.max(-previousNLLH, nbGrad._llh);
-        double bestTheta = theta;
+        Vec mu = Scope.track(gPred.outputFrame(Key.make(), new String[]{"prediction"}, null)).vec(0);
 
-        delta = Double.isFinite(delta) ? delta : 1; // NaN can occur in extreme datasets so try to get out of this neighborhood just by linesearch
-
-        // Golden section search for the optimal size of delta
-        // Set lowerbound to -10 or lowest value that will keep theta > 0 which ever is bigger
-        // Negative value here helps with datasets where we use to diverge, I'm not sure yet if it's caused by some
-        // numerical issues or if the likelihood can get multimodal for some cases.
-        double lowerBound = (theta + 10 * delta < 0) ? (1 - 1e-15) * theta / delta : -10;
-        double upperBound = (theta - 1e3 * delta < 0) ? (1 - 1e-15) * theta / delta : 1e3;
-        double d = upperBound - lowerBound;
-
-        for (int i = 0; i < _parms._max_iterations_dispersion; i++) {
-          d *= 0.618;  // division by golden ratio
-          final double lowerBoundProposal = upperBound - d;
-          final double upperBoundProposal = lowerBound + d;
-          NegativeBinomialGradientAndHessian nbLower = new NegativeBinomialGradientAndHessian(theta - lowerBoundProposal * delta).doAll(mu, response, weights);
-          NegativeBinomialGradientAndHessian nbUpper = new NegativeBinomialGradientAndHessian(theta - upperBoundProposal * delta).doAll(mu, response, weights);
-
-          if (nbLower._llh >= nbUpper._llh) {
-            upperBound = upperBoundProposal;
-            if (nbLower._llh > bestLLH) {
-              bestLLH = nbLower._llh;
-              bestTheta = nbLower._theta;
+        if (iterCnt == 1) {
+          theta = estimateNegBinomialDispersionMomentMethod(_model, betaCnd, _dinfo, weights, response, mu);
+        } else {
+          theta = _parms._theta;
+          NegativeBinomialGradientAndHessian nbGrad = new NegativeBinomialGradientAndHessian(theta).doAll(mu, response, weights);
+          delta = _parms._dispersion_learning_rate * nbGrad._grad / nbGrad._hess;
+          double bestLLH = Math.max(-previousNLLH, nbGrad._llh);
+          double bestTheta = theta;
+
+          delta = Double.isFinite(delta) ? delta : 1; // NaN can occur in extreme datasets so try to get out of this neighborhood just by linesearch
+
+          // Golden section search for the optimal size of delta
+          // Set lowerbound to -10 or lowest value that will keep theta > 0 which ever is bigger
+          // Negative value here helps with datasets where we use to diverge, I'm not sure yet if it's caused by some
+          // numerical issues or if the likelihood can get multimodal for some cases.
+          double lowerBound = (theta + 10 * delta < 0) ? (1 - 1e-15) * theta / delta : -10;
+          double upperBound = (theta - 1e3 * delta < 0) ? (1 - 1e-15) * theta / delta : 1e3;
+          double d = upperBound - lowerBound;
+
+          for (int i = 0; i < _parms._max_iterations_dispersion; i++) {
+            d *= 0.618;  // division by golden ratio
+            final double lowerBoundProposal = upperBound - d;
+            final double upperBoundProposal = lowerBound + d;
+            NegativeBinomialGradientAndHessian nbLower = new NegativeBinomialGradientAndHessian(theta - lowerBoundProposal * delta).doAll(mu, response, weights);
+            NegativeBinomialGradientAndHessian nbUpper = new NegativeBinomialGradientAndHessian(theta - upperBoundProposal * delta).doAll(mu, response, weights);
+
+            if (nbLower._llh >= nbUpper._llh) {
+              upperBound = upperBoundProposal;
+              if (nbLower._llh > bestLLH) {
+                bestLLH = nbLower._llh;
+                bestTheta = nbLower._theta;
+              }
+            } else {
+              lowerBound = lowerBoundProposal;
+              if (nbUpper._llh > bestLLH) {
+                bestLLH = nbUpper._llh;
+                bestTheta = nbUpper._theta;
+              }
             }
-          } else {
-            lowerBound = lowerBoundProposal;
-            if (nbUpper._llh > bestLLH) {
-              bestLLH = nbUpper._llh;
-              bestTheta = nbUpper._theta;
+            if (Math.abs((upperBoundProposal - lowerBoundProposal) * Math.max(1, delta / Math.max(_parms._theta, bestTheta))) < _parms._dispersion_epsilon || _job.stop_requested()) {
+              break;
             }
           }
-          if (Math.abs((upperBoundProposal - lowerBoundProposal) * Math.max(1, delta / Math.max(_parms._theta, bestTheta))) < _parms._dispersion_epsilon || _job.stop_requested()) {
-            break;
-          }
+
+          theta = bestTheta;
+          converged = (nbGrad._llh + previousNLLH) <= _parms._objective_epsilon || !Double.isFinite(theta);
         }
+        delta = _parms._theta - theta;
+        converged = converged && (Math.abs(delta) / Math.max(_parms._theta, theta) < _parms._dispersion_epsilon);
 
-        theta = bestTheta;
-        converged = (nbGrad._llh + previousNLLH) <= _parms._objective_epsilon || !Double.isFinite(theta);
+        updateTheta(theta);
+        return converged;
+      } finally {
+        Scope.exit();
       }
-      delta = _parms._theta - theta;
-      converged = converged &&  (Math.abs(delta) / Math.max(_parms._theta, theta) < _parms._dispersion_epsilon);
-
-      updateTheta(theta);
-      return converged;
     }
 
     private void fitLBFGS() {
@@ -2773,7 +2801,7 @@ else if (gaussian.equals(_parms._family) && Link.identity.equals(_parms._link))
         }
       }
 
-      if (_parms._compute_p_values) { // compute p-values, standard error, estimate dispersion parameters...
+      if (_parms._compute_p_values) { // compute p-values, standard error, estimate dispersion parameters... 
         double se = _parms._init_dispersion_parameter;
         boolean seEst = false;
         double[] beta = _state.beta();  // standardized if _parms._standardize=true, original otherwise
@@ -2802,6 +2830,8 @@ else if (gaussian.equals(_parms._family) && Link.identity.equals(_parms._link))
               se = estimateTweedieDispersionOnly(_parms, _model, _job, beta, _state.activeData());
             }
           }
+          // save estimation to the _params, so it is available for params.likelihood computation
+          _parms._dispersion_estimated = se;
         }
         double[] zvalues = MemoryManager.malloc8d(_state.activeData().fullN() + 1);
        // double[][] inv = cholInv(); // from non-standardized predictors
@@ -3457,8 +3487,11 @@ private void doCompute() {
         _model.setVcov(_vcov);
         _model.update(_job._key);
       }
-      if (!_parms._HGLM)  // no need to do for HGLM
+      if (!_parms._HGLM) {  // no need to do for HGLM
+        _model._finalScoring = true; // enables likelihood calculation while scoring
         scoreAndUpdateModel();
+        _model._finalScoring = false; // avoid calculating likelihood in case of further updates
+      }
 
       if (dfbetas.equals(_parms._influence))
         genRID();