X-Git-Url: https://pintos-os.org/cgi-bin/gitweb.cgi?a=blobdiff_plain;f=src%2Flanguage%2Fstats%2Flogistic.c;fp=src%2Flanguage%2Fstats%2Flogistic.c;h=e3a66f4de004c81d1c18bd5b80f358fa27465a0b;hb=5937fc72919993c44ae9956a44fb53de797483b5;hp=4645ee14d208356d816ba386c73e3f45a8081033;hpb=53f7aa8f9d84203b0f62c6964ebf2ab6e0e51584;p=pspp

diff --git a/src/language/stats/logistic.c b/src/language/stats/logistic.c
index 4645ee14d2..e3a66f4de0 100644
--- a/src/language/stats/logistic.c
+++ b/src/language/stats/logistic.c
@@ -344,13 +344,13 @@ xt_times_y_pi (const struct lr_spec *cmd,
 	       struct casereader *input,
 	       const struct variable **x, size_t n_x,
 	       const struct variable *y_var,
-	       double *likelihood)
+	       double *llikelihood)
 {
   struct casereader *reader;
   struct ccase *c;
   gsl_vector *output = gsl_vector_calloc (res->beta_hat->size);
 
-  *likelihood = 1.0;
+  *llikelihood = 0.0;
   res->tn = res->tp = res->fn = res->fp = 0;
   for (reader = casereader_clone (input);
        (c = casereader_read (reader)) != NULL; case_unref (c))
@@ -363,7 +363,7 @@ xt_times_y_pi (const struct lr_spec *cmd,
 
       double y = map_dependent_var (cmd, res, case_data (c, y_var));
 
-      *likelihood *= pow (pi, weight * y) * pow (1 - pi, weight * (1 - y));
+      *llikelihood += (weight * y) * log (pi) + log (1 - pi) * weight * (1 - y);
 
       for (v0 = 0; v0 < res->beta_hat->size; ++v0)
 	{
@@ -578,10 +578,10 @@ run_lr (const struct lr_spec *cmd, struct casereader *input,
 
   bool converged = false;
 
-  /* Set the likelihoods to a negative sentinel value */
-  double likelihood = -1;
-  double prev_likelihood = -1;
-  double initial_likelihood = -1;
+  /* Set the log likelihoods to a sentinel value */
+  double log_likelihood = SYSMIS;
+  double prev_log_likelihood = SYSMIS;
+  double initial_log_likelihood = SYSMIS;
 
   struct lr_result work;
   work.n_missing = 0;
@@ -643,7 +643,7 @@ run_lr (const struct lr_spec *cmd, struct casereader *input,
       v = xt_times_y_pi (cmd, &work, input,
 			 cmd->predictor_vars, cmd->n_predictor_vars,
 			 cmd->dep_var,
-			 &likelihood);
+			 &log_likelihood);
 
       {
 	/* delta = M.v */
@@ -666,29 +666,29 @@ run_lr (const struct lr_spec *cmd, struct casereader *input,
 	gsl_vector_free (delta);
       }
 
-      if ( prev_likelihood >= 0)
+      if (i > 0)
 	{
-	  if (-log (likelihood) > -(1.0 - cmd->lcon) * log (prev_likelihood))
+	  if (-log_likelihood > -(1.0 - cmd->lcon) * prev_log_likelihood)
 	    {
 	      msg (MN, _("Estimation terminated at iteration number %d because Log Likelihood decreased by less than %g%%"), i + 1, 100 * cmd->lcon);
 	      converged = true;
 	    }
 	}
       if (i == 0)
-	initial_likelihood = likelihood;
-      prev_likelihood = likelihood;
+	initial_log_likelihood = log_likelihood;
+      prev_log_likelihood = log_likelihood;
 
       if (converged)
 	break;
     }
   casereader_destroy (input);
-  assert (initial_likelihood >= 0);
+
 
   if ( ! converged) 
     msg (MW, _("Estimation terminated at iteration number %d because maximum iterations has been reached"), i );
 
 
-  output_model_summary (&work, initial_likelihood, likelihood);
+  output_model_summary (&work, initial_log_likelihood, log_likelihood);
 
   if (work.cats)
     output_categories (cmd, &work);
@@ -1307,7 +1307,7 @@ output_variables (const struct lr_spec *cmd,
 /* Show the model summary box */
 static void
 output_model_summary (const struct lr_result *res,
-		      double initial_likelihood, double likelihood)
+		      double initial_log_likelihood, double log_likelihood)
 {
   const int heading_columns = 0;
   const int heading_rows = 1;
@@ -1329,15 +1329,15 @@ output_model_summary (const struct lr_result *res,
 
   tab_text (t,  0, 0, TAB_LEFT | TAT_TITLE, _("Step 1"));
   tab_text (t,  1, 0, TAB_CENTER | TAT_TITLE, _("-2 Log likelihood"));
-  tab_double (t,  1, 1, 0, -2 * log (likelihood), 0);
+  tab_double (t,  1, 1, 0, -2 * log_likelihood, 0);
 
 
   tab_text (t,  2, 0, TAB_CENTER | TAT_TITLE, _("Cox & Snell R Square"));
-  cox =  1.0 - pow (initial_likelihood /likelihood, 2 / res->cc);
+  cox =  1.0 - exp((initial_log_likelihood - log_likelihood) * (2 / res->cc));
   tab_double (t,  2, 1, 0, cox, 0);
 
   tab_text (t,  3, 0, TAB_CENTER | TAT_TITLE, _("Nagelkerke R Square"));
-  tab_double (t,  3, 1, 0, cox / ( 1.0 - pow (initial_likelihood, 2 / res->cc)), 0);
+  tab_double (t,  3, 1, 0, cox / ( 1.0 - exp(initial_log_likelihood * (2 / res->cc))), 0);
 
 
   tab_submit (t);