projects
/
pspp
/ blobdiff
commit
grep
author
committer
pickaxe
?
search:
re
summary
|
shortlog
|
log
|
commit
|
commitdiff
|
tree
raw
|
inline
| side by side
Logistic Regression: Calculate the logs of the likelihood instead of the raw value.
[pspp]
/
src
/
language
/
stats
/
logistic.c
diff --git
a/src/language/stats/logistic.c
b/src/language/stats/logistic.c
index 4645ee14d208356d816ba386c73e3f45a8081033..e3a66f4de004c81d1c18bd5b80f358fa27465a0b 100644
(file)
--- 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,
struct casereader *input,
const struct variable **x, size_t n_x,
const struct variable *y_var,
- double *likelihood)
+ double *l
l
ikelihood)
{
struct casereader *reader;
struct ccase *c;
gsl_vector *output = gsl_vector_calloc (res->beta_hat->size);
{
struct casereader *reader;
struct ccase *c;
gsl_vector *output = gsl_vector_calloc (res->beta_hat->size);
- *l
ikelihood = 1
.0;
+ *l
likelihood = 0
.0;
res->tn = res->tp = res->fn = res->fp = 0;
for (reader = casereader_clone (input);
(c = casereader_read (reader)) != NULL; case_unref (c))
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));
double y = map_dependent_var (cmd, res, case_data (c, y_var));
- *l
ikelihood *= pow (pi, weight * y) * pow (1 - pi, weight * (1 - y)
);
+ *l
likelihood += (weight * y) * log (pi) + log (1 - pi) * weight * (1 - y
);
for (v0 = 0; v0 < res->beta_hat->size; ++v0)
{
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;
bool converged = false;
- /* Set the l
ikelihoods to a negative
sentinel value */
- double l
ikelihood = -1
;
- double prev_l
ikelihood = -1
;
- double initial_l
ikelihood = -1
;
+ /* Set the l
og likelihoods to a
sentinel value */
+ double l
og_likelihood = SYSMIS
;
+ double prev_l
og_likelihood = SYSMIS
;
+ double initial_l
og_likelihood = SYSMIS
;
struct lr_result work;
work.n_missing = 0;
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,
v = xt_times_y_pi (cmd, &work, input,
cmd->predictor_vars, cmd->n_predictor_vars,
cmd->dep_var,
- &likelihood);
+ &l
og_l
ikelihood);
{
/* delta = M.v */
{
/* delta = M.v */
@@
-666,29
+666,29
@@
run_lr (const struct lr_spec *cmd, struct casereader *input,
gsl_vector_free (delta);
}
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)
{
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_l
ikelihood =
likelihood;
- prev_l
ikelihood =
likelihood;
+ initial_l
og_likelihood = log_
likelihood;
+ prev_l
og_likelihood = log_
likelihood;
if (converged)
break;
}
casereader_destroy (input);
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 );
if ( ! converged)
msg (MW, _("Estimation terminated at iteration number %d because maximum iterations has been reached"), i );
- output_model_summary (&work, initial_l
ikelihood,
likelihood);
+ output_model_summary (&work, initial_l
og_likelihood, log_
likelihood);
if (work.cats)
output_categories (cmd, &work);
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,
/* Show the model summary box */
static void
output_model_summary (const struct lr_result *res,
- double initial_l
ikelihood, double
likelihood)
+ double initial_l
og_likelihood, double log_
likelihood)
{
const int heading_columns = 0;
const int heading_rows = 1;
{
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_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"));
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, 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);
tab_submit (t);