projects
/
pspp
/ blobdiff
commit
grep
author
committer
pickaxe
?
search:
re
summary
|
shortlog
|
log
|
commit
|
commitdiff
|
tree
raw
|
inline
| side by side
Yet more sealage
[pspp]
/
src
/
language
/
stats
/
logistic.c
diff --git
a/src/language/stats/logistic.c
b/src/language/stats/logistic.c
index 4645ee14d208356d816ba386c73e3f45a8081033..b3cf843400f6460b077de5b396e2e7c10cada8d1 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)
{
@@
-478,6
+478,11
@@
initial_pass (const struct lr_spec *cmd, struct lr_result *res, struct casereade
double weight = dict_get_case_weight (cmd->dict, c, &res->warn_bad_weight);
const union value *depval = case_data (c, cmd->dep_var);
double weight = dict_get_case_weight (cmd->dict, c, &res->warn_bad_weight);
const union value *depval = case_data (c, cmd->dep_var);
+ if (var_is_value_missing (cmd->dep_var, depval, cmd->exclude))
+ {
+ missing = true;
+ }
+ else
for (v = 0; v < cmd->n_indep_vars; ++v)
{
const union value *val = case_data (c, cmd->indep_vars[v]);
for (v = 0; v < cmd->n_indep_vars; ++v)
{
const union value *val = case_data (c, cmd->indep_vars[v]);
@@
-578,10
+583,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;
@@
-589,11
+594,12
@@
run_lr (const struct lr_spec *cmd, struct casereader *input,
work.warn_bad_weight = true;
work.cats = NULL;
work.beta_hat = NULL;
work.warn_bad_weight = true;
work.cats = NULL;
work.beta_hat = NULL;
+ work.hessian = NULL;
/* Get the initial estimates of \beta and their standard errors.
And perform other auxilliary initialisation. */
if (! initial_pass (cmd, &work, input))
/* Get the initial estimates of \beta and their standard errors.
And perform other auxilliary initialisation. */
if (! initial_pass (cmd, &work, input))
-
return false
;
+
goto error
;
for (i = 0; i < cmd->n_cat_predictors; ++i)
{
for (i = 0; i < cmd->n_cat_predictors; ++i)
{
@@
-607,7
+613,7
@@
run_lr (const struct lr_spec *cmd, struct casereader *input,
msg (ME, _("Category %s does not have at least two distinct values. Logistic regression will not be run."),
ds_cstr(&str));
ds_destroy (&str);
msg (ME, _("Category %s does not have at least two distinct values. Logistic regression will not be run."),
ds_cstr(&str));
ds_destroy (&str);
-
return false
;
+
goto error
;
}
}
}
}
@@
-623,6
+629,12
@@
run_lr (const struct lr_spec *cmd, struct casereader *input,
NULL,
NULL);
NULL,
NULL);
+ input = casereader_create_filter_missing (input,
+ &cmd->dep_var,
+ 1,
+ cmd->exclude,
+ NULL,
+ NULL);
work.hessian = gsl_matrix_calloc (work.beta_hat->size, work.beta_hat->size);
work.hessian = gsl_matrix_calloc (work.beta_hat->size, work.beta_hat->size);
@@
-643,7
+655,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
+678,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;
}
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);
@@
-696,12
+708,20
@@
run_lr (const struct lr_spec *cmd, struct casereader *input,
output_classification_table (cmd, &work);
output_variables (cmd, &work);
output_classification_table (cmd, &work);
output_variables (cmd, &work);
+ casereader_destroy (input);
gsl_matrix_free (work.hessian);
gsl_vector_free (work.beta_hat);
gsl_matrix_free (work.hessian);
gsl_vector_free (work.beta_hat);
-
categoricals_destroy (work.cats);
return true;
categoricals_destroy (work.cats);
return true;
+
+ error:
+ casereader_destroy (input);
+ gsl_matrix_free (work.hessian);
+ gsl_vector_free (work.beta_hat);
+ categoricals_destroy (work.cats);
+
+ return false;
}
struct variable_node
}
struct variable_node
@@
-730,6
+750,7
@@
lookup_variable (const struct hmap *map, const struct variable *var, unsigned in
int
cmd_logistic (struct lexer *lexer, struct dataset *ds)
{
int
cmd_logistic (struct lexer *lexer, struct dataset *ds)
{
+ int i;
/* Temporary location for the predictor variables.
These may or may not include the categorical predictors */
const struct variable **pred_vars;
/* Temporary location for the predictor variables.
These may or may not include the categorical predictors */
const struct variable **pred_vars;
@@
-1081,6
+1102,10
@@
cmd_logistic (struct lexer *lexer, struct dataset *ds)
ok = proc_commit (ds) && ok;
}
ok = proc_commit (ds) && ok;
}
+ for (i = 0 ; i < lr.n_cat_predictors; ++i)
+ {
+ interaction_destroy (lr.cat_predictors[i]);
+ }
free (lr.predictor_vars);
free (lr.cat_predictors);
free (lr.indep_vars);
free (lr.predictor_vars);
free (lr.cat_predictors);
free (lr.indep_vars);
@@
-1089,6
+1114,10
@@
cmd_logistic (struct lexer *lexer, struct dataset *ds)
error:
error:
+ for (i = 0 ; i < lr.n_cat_predictors; ++i)
+ {
+ interaction_destroy (lr.cat_predictors[i]);
+ }
free (lr.predictor_vars);
free (lr.cat_predictors);
free (lr.indep_vars);
free (lr.predictor_vars);
free (lr.cat_predictors);
free (lr.indep_vars);
@@
-1289,10
+1318,14
@@
output_variables (const struct lr_spec *cmd,
if (cmd->print & PRINT_CI)
{
if (cmd->print & PRINT_CI)
{
+ int last_ci = nr;
double wc = gsl_cdf_ugaussian_Pinv (0.5 + cmd->confidence / 200.0);
wc *= sqrt (sigma2);
double wc = gsl_cdf_ugaussian_Pinv (0.5 + cmd->confidence / 200.0);
wc *= sqrt (sigma2);
- if (idx < cmd->n_predictor_vars)
+ if (cmd->constant)
+ last_ci--;
+
+ if (row < last_ci)
{
tab_double (t, 8, row, 0, exp (b - wc), 0);
tab_double (t, 9, row, 0, exp (b + wc), 0);
{
tab_double (t, 8, row, 0, exp (b - wc), 0);
tab_double (t, 9, row, 0, exp (b + wc), 0);
@@
-1307,7
+1340,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
+1362,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);