X-Git-Url: https://pintos-os.org/cgi-bin/gitweb.cgi?a=blobdiff_plain;f=src%2Flanguage%2Fstats%2Fglm.c;h=1eaf1a988a335264a5db70d9cd3bc137363c7602;hb=d725a9fcccd8c826cda4278294419024d04efc90;hp=d2fdde5621d0baf34fec0bc2190bb704f2dccd43;hpb=01c3f002d3a9ad239bb145c4a20dccdef75b0775;p=pspp diff --git a/src/language/stats/glm.c b/src/language/stats/glm.c index d2fdde5621..1eaf1a988a 100644 --- a/src/language/stats/glm.c +++ b/src/language/stats/glm.c @@ -39,6 +39,7 @@ #include "linreg/sweep.h" #include "math/categoricals.h" #include "math/covariance.h" +#include "math/interaction.h" #include "math/moments.h" #include "output/tab.h" @@ -53,12 +54,19 @@ struct glm_spec size_t n_factor_vars; const struct variable **factor_vars; + size_t n_interactions; + struct interaction **interactions; + enum mv_class exclude; /* The weight variable */ const struct variable *wv; + const struct dictionary *dict; + bool intercept; + + double alpha; }; struct glm_workspace @@ -81,29 +89,37 @@ static void output_glm (const struct glm_spec *, static void run_glm (struct glm_spec *cmd, struct casereader *input, const struct dataset *ds); + +static bool parse_design_spec (struct lexer *lexer, struct glm_spec *glm); + + int cmd_glm (struct lexer *lexer, struct dataset *ds) { + int i; struct const_var_set *factors = NULL; - const struct dictionary *dict = dataset_dict (ds); struct glm_spec glm; + bool design = false; + glm.dict = dataset_dict (ds); glm.n_dep_vars = 0; glm.n_factor_vars = 0; + glm.n_interactions = 0; + glm.interactions = NULL; glm.dep_vars = NULL; glm.factor_vars = NULL; glm.exclude = MV_ANY; glm.intercept = true; - glm.wv = dict_get_weight (dict); - + glm.wv = dict_get_weight (glm.dict); + glm.alpha = 0.05; - if (!parse_variables_const (lexer, dict, + if (!parse_variables_const (lexer, glm.dict, &glm.dep_vars, &glm.n_dep_vars, PV_NO_DUPLICATE | PV_NUMERIC)) goto error; lex_force_match (lexer, T_BY); - if (!parse_variables_const (lexer, dict, + if (!parse_variables_const (lexer, glm.dict, &glm.factor_vars, &glm.n_factor_vars, PV_NO_DUPLICATE | PV_NUMERIC)) goto error; @@ -163,16 +179,84 @@ cmd_glm (struct lexer *lexer, struct dataset *ds) } } } -#if 0 + else if (lex_match_id (lexer, "CRITERIA")) + { + lex_match (lexer, T_EQUALS); + if (lex_match_id (lexer, "ALPHA")) + { + if (lex_force_match (lexer, T_LPAREN)) + { + if (! lex_force_num (lexer)) + { + lex_error (lexer, NULL); + goto error; + } + + glm.alpha = lex_number (lexer); + lex_get (lexer); + if ( ! lex_force_match (lexer, T_RPAREN)) + { + lex_error (lexer, NULL); + goto error; + } + } + } + else + { + lex_error (lexer, NULL); + goto error; + } + } + else if (lex_match_id (lexer, "METHOD")) + { + lex_match (lexer, T_EQUALS); + if ( !lex_force_match_id (lexer, "SSTYPE")) + { + lex_error (lexer, NULL); + goto error; + } + + if ( ! lex_force_match (lexer, T_LPAREN)) + { + lex_error (lexer, NULL); + goto error; + } + + if ( ! lex_force_int (lexer)) + { + lex_error (lexer, NULL); + goto error; + } + + if (3 != lex_integer (lexer)) + { + msg (ME, _("Only type 3 sum of squares are currently implemented")); + goto error; + } + + lex_get (lexer); + + if ( ! lex_force_match (lexer, T_RPAREN)) + { + lex_error (lexer, NULL); + goto error; + } + } else if (lex_match_id (lexer, "DESIGN")) { - size_t n_des; - const struct variable **des; lex_match (lexer, T_EQUALS); - parse_const_var_set_vars (lexer, factors, &des, &n_des, 0); + if (! parse_design_spec (lexer, &glm)) + goto error; + + if ( glm.n_interactions == 0) + { + msg (ME, _("One or more design variables must be given")); + goto error; + } + + design = true; } -#endif else { lex_error (lexer, NULL); @@ -180,13 +264,18 @@ cmd_glm (struct lexer *lexer, struct dataset *ds) } } + if ( ! design ) + { + lex_error (lexer, _("/DESIGN is mandatory in GLM")); + goto error; + } { struct casegrouper *grouper; struct casereader *group; bool ok; - grouper = casegrouper_create_splits (proc_open (ds), dict); + grouper = casegrouper_create_splits (proc_open (ds), glm.dict); while (casegrouper_get_next_group (grouper, &group)) run_glm (&glm, group, ds); ok = casegrouper_destroy (grouper); @@ -195,14 +284,22 @@ cmd_glm (struct lexer *lexer, struct dataset *ds) const_var_set_destroy (factors); free (glm.factor_vars); + for (i = 0 ; i < glm.n_interactions; ++i) + interaction_destroy (glm.interactions[i]); + free (glm.interactions); free (glm.dep_vars); + return CMD_SUCCESS; error: const_var_set_destroy (factors); free (glm.factor_vars); + for (i = 0 ; i < glm.n_interactions; ++i) + interaction_destroy (glm.interactions[i]); + + free (glm.interactions); free (glm.dep_vars); return CMD_FAILURE; @@ -225,38 +322,33 @@ not_dropped (size_t j, size_t * dropped, size_t n_dropped) } static void -get_ssq (struct covariance *cov, gsl_vector * ssq, const struct glm_spec *cmd) +get_ssq (struct covariance *cov, gsl_vector *ssq, const struct glm_spec *cmd) { - const struct variable **vars; - gsl_matrix *small_cov = NULL; gsl_matrix *cm = covariance_calculate_unnormalized (cov); size_t i; size_t j; size_t k; - size_t n; - size_t m; - size_t *dropped; - size_t n_dropped; - - dropped = xcalloc (covariance_dim (cov), sizeof (*dropped)); - vars = xcalloc (covariance_dim (cov), sizeof (*vars)); - covariance_get_var_indices (cov, vars); + size_t *dropped = xcalloc (covariance_dim (cov), sizeof (*dropped)); + const struct categoricals *cats = covariance_get_categoricals (cov); - for (k = 0; k < cmd->n_factor_vars; k++) + for (k = 0; k < cmd->n_interactions; k++) { - n_dropped = 0; - for (i = 1; i < covariance_dim (cov); i++) + size_t n = 0; + size_t m = 0; + gsl_matrix *small_cov = NULL; + size_t n_dropped = 0; + for (i = cmd->n_dep_vars; i < covariance_dim (cov); i++) { - if (vars[i] == cmd->factor_vars[k]) + if (categoricals_get_interaction_by_subscript (cats, i - cmd->n_dep_vars) + == cmd->interactions[k]) { + assert (n_dropped < covariance_dim (cov)); dropped[n_dropped++] = i; } } small_cov = gsl_matrix_alloc (cm->size1 - n_dropped, cm->size2 - n_dropped); gsl_matrix_set (small_cov, 0, 0, gsl_matrix_get (cm, 0, 0)); - n = 0; - m = 0; for (i = 0; i < cm->size1; i++) { if (not_dropped (i, dropped, n_dropped)) @@ -282,7 +374,6 @@ get_ssq (struct covariance *cov, gsl_vector * ssq, const struct glm_spec *cmd) } free (dropped); - free (vars); gsl_matrix_free (cm); } @@ -301,7 +392,8 @@ run_glm (struct glm_spec *cmd, struct casereader *input, struct glm_workspace ws; struct covariance *cov; - ws.cats = categoricals_create (cmd->factor_vars, cmd->n_factor_vars, + + ws.cats = categoricals_create (cmd->interactions, cmd->n_interactions, cmd->wv, cmd->exclude, NULL, NULL, NULL, NULL); @@ -335,8 +427,6 @@ run_glm (struct glm_spec *cmd, struct casereader *input, } casereader_destroy (reader); - categoricals_done (ws.cats); - for (reader = input; (c = casereader_read (reader)) != NULL; case_unref (c)) { @@ -360,8 +450,8 @@ run_glm (struct glm_spec *cmd, struct casereader *input, reg_sweep (cm, 0); /* - Store the overall SSE. - */ + Store the overall SSE. + */ ws.ssq = gsl_vector_alloc (cm->size1); gsl_vector_set (ws.ssq, 0, gsl_matrix_get (cm, 0, 0)); get_ssq (cov, ws.ssq, cmd); @@ -398,7 +488,7 @@ output_glm (const struct glm_spec *cmd, const struct glm_workspace *ws) struct tab_table *t; const int nc = 6; - int nr = heading_rows + 4 + cmd->n_factor_vars; + int nr = heading_rows + 4 + cmd->n_interactions; if (cmd->intercept) nr++; @@ -427,8 +517,7 @@ output_glm (const struct glm_spec *cmd, const struct glm_workspace *ws) if (cmd->intercept) df_corr += 1.0; - for (f = 0; f < cmd->n_factor_vars; ++f) - df_corr += categoricals_n_count (ws->cats, f) - 1.0; + df_corr += categoricals_df_total (ws->cats); mse = gsl_vector_get (ws->ssq, 0) / (n_total - df_corr); @@ -452,13 +541,15 @@ output_glm (const struct glm_spec *cmd, const struct glm_workspace *ws) r++; } - for (f = 0; f < cmd->n_factor_vars; ++f) + for (f = 0; f < cmd->n_interactions; ++f) { - const double df = categoricals_n_count (ws->cats, f) - 1.0; + struct string str = DS_EMPTY_INITIALIZER; + const double df = categoricals_df (ws->cats, f); const double ssq = gsl_vector_get (ws->ssq, f + 1); const double F = ssq / df / mse; - tab_text (t, 0, r, TAB_LEFT | TAT_TITLE, - var_to_string (cmd->factor_vars[f])); + interaction_to_string (cmd->interactions[f], &str); + tab_text (t, 0, r, TAB_LEFT | TAT_TITLE, ds_cstr (&str)); + ds_destroy (&str); tab_double (t, 1, r, 0, ssq, NULL); tab_double (t, 2, r, 0, df, wfmt); @@ -467,8 +558,6 @@ output_glm (const struct glm_spec *cmd, const struct glm_workspace *ws) tab_double (t, 5, r, 0, gsl_cdf_fdist_Q (F, df, n_total - df_corr), NULL); - - r++; } @@ -534,3 +623,126 @@ dump_matrix (const gsl_matrix * m) printf ("\n"); } #endif + + + + +/* Match a variable. + If the match succeeds, the variable will be placed in VAR. + Returns true if successful */ +static bool +lex_match_variable (struct lexer *lexer, const struct glm_spec *glm, const struct variable **var) +{ + if (lex_token (lexer) != T_ID) + return false; + + *var = parse_variable_const (lexer, glm->dict); + + if ( *var == NULL) + return false; + return true; +} + +/* An interaction is a variable followed by {*, BY} followed by an interaction */ +static bool +parse_design_interaction (struct lexer *lexer, struct glm_spec *glm, struct interaction **iact) +{ + const struct variable *v = NULL; + assert (iact); + + switch (lex_next_token (lexer, 1)) + { + case T_ENDCMD: + case T_SLASH: + case T_COMMA: + case T_ID: + case T_BY: + case T_ASTERISK: + break; + default: + return false; + break; + } + + if (! lex_match_variable (lexer, glm, &v)) + { + interaction_destroy (*iact); + *iact = NULL; + return false; + } + + assert (v); + + if ( *iact == NULL) + *iact = interaction_create (v); + else + interaction_add_variable (*iact, v); + + if ( lex_match (lexer, T_ASTERISK) || lex_match (lexer, T_BY)) + { + // lex_error (lexer, "Interactions are not yet implemented"); return false; + return parse_design_interaction (lexer, glm, iact); + } + + return true; +} + +static bool +parse_nested_variable (struct lexer *lexer, struct glm_spec *glm) +{ + const struct variable *v = NULL; + if ( ! lex_match_variable (lexer, glm, &v)) + return false; + + if (lex_match (lexer, T_LPAREN)) + { + if ( ! parse_nested_variable (lexer, glm)) + return false; + + if ( ! lex_force_match (lexer, T_RPAREN)) + return false; + } + + lex_error (lexer, "Nested variables are not yet implemented"); return false; + return true; +} + +/* A design term is an interaction OR a nested variable */ +static bool +parse_design_term (struct lexer *lexer, struct glm_spec *glm) +{ + struct interaction *iact = NULL; + if (parse_design_interaction (lexer, glm, &iact)) + { + /* Interaction parsing successful. Add to list of interactions */ + glm->interactions = xrealloc (glm->interactions, sizeof *glm->interactions * ++glm->n_interactions); + glm->interactions[glm->n_interactions - 1] = iact; + return true; + } + + if ( parse_nested_variable (lexer, glm)) + return true; + + return false; +} + + + +/* Parse a complete DESIGN specification. + A design spec is a design term, optionally followed by a comma, + and another design spec. +*/ +static bool +parse_design_spec (struct lexer *lexer, struct glm_spec *glm) +{ + if (lex_token (lexer) == T_ENDCMD || lex_token (lexer) == T_SLASH) + return true; + + if ( ! parse_design_term (lexer, glm)) + return false; + + lex_match (lexer, T_COMMA); + + return parse_design_spec (lexer, glm); +} +