From: John Darrington Date: Thu, 14 Jul 2011 21:43:11 +0000 (+0200) Subject: Extended the glm command to accept interactions. Unfortunately the ssqs are wrong X-Git-Tag: v0.7.9~193 X-Git-Url: https://pintos-os.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=d18270eb353b0f8ed23814be18f1de69c4001262;p=pspp-builds.git Extended the glm command to accept interactions. Unfortunately the ssqs are wrong --- diff --git a/src/language/stats/glm.c b/src/language/stats/glm.c index 3194fd81..8e3b07f4 100644 --- a/src/language/stats/glm.c +++ b/src/language/stats/glm.c @@ -322,7 +322,7 @@ 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) { gsl_matrix *cm = covariance_calculate_unnormalized (cov); size_t i; @@ -519,8 +519,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_interactions; ++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); @@ -547,7 +546,7 @@ output_glm (const struct glm_spec *cmd, const struct glm_workspace *ws) for (f = 0; f < cmd->n_interactions; ++f) { struct string str = DS_EMPTY_INITIALIZER; - const double df = categoricals_n_count (ws->cats, f) - 1.0; + const double df = categoricals_df (ws->cats, f); const double ssq = gsl_vector_get (ws->ssq, f + 1); const double F = ssq / df / mse; interaction_to_string (cmd->interactions[f], &str); diff --git a/src/language/stats/oneway.c b/src/language/stats/oneway.c index 96713280..e15bff74 100644 --- a/src/language/stats/oneway.c +++ b/src/language/stats/oneway.c @@ -770,7 +770,7 @@ run_oneway (const struct oneway_spec *cmd, pvw->ssa = pvw->sst - pvw->sse; - pvw->n_groups = categoricals_total (cats); + pvw->n_groups = categoricals_n_total (cats); pvw->mse = (pvw->sst - pvw->ssa) / (pvw->n - pvw->n_groups); @@ -783,8 +783,8 @@ run_oneway (const struct oneway_spec *cmd, categoricals_done (cats); - if (categoricals_total (cats) > ws.actual_number_of_groups) - ws.actual_number_of_groups = categoricals_total (cats); + if (categoricals_n_total (cats) > ws.actual_number_of_groups) + ws.actual_number_of_groups = categoricals_n_total (cats); } casereader_destroy (input); @@ -1015,7 +1015,7 @@ show_descriptives (const struct oneway_spec *cmd, const struct oneway_workspace if ( v > 0) tab_hline (t, TAL_1, 0, n_cols - 1, row); - for (count = 0; count < categoricals_total (cats); ++count) + for (count = 0; count < categoricals_n_total (cats); ++count) { double T; double n, mean, variance; @@ -1104,7 +1104,7 @@ show_descriptives (const struct oneway_spec *cmd, const struct oneway_workspace tab_double (t, 9, row + count, 0, ws->dd_total[v]->maximum, fmt); } - row += categoricals_total (cats) + 1; + row += categoricals_n_total (cats) + 1; } tab_submit (t); diff --git a/src/math/categoricals.c b/src/math/categoricals.c index 218ac6c1..c955ff95 100644 --- a/src/math/categoricals.c +++ b/src/math/categoricals.c @@ -28,42 +28,100 @@ #include "libpspp/hmap.h" #include "libpspp/pool.h" #include "libpspp/str.h" +#include "libpspp/hash-functions.h" #include "gl/xalloc.h" struct value_node { struct hmap_node node; /* Node in hash map. */ - struct ccase *ccase; - double cc; /* The total of the weights of cases with this value */ + + union value val; /* The value */ void *user_data; /* A pointer to data which the caller can store stuff */ +}; + +struct interaction_value +{ + struct hmap_node node; /* Node in hash map. */ + + struct ccase *ccase; /* A case (probably the first in the dataset) which matches this value */ - int subscript; /* A zero based integer, unique within the variable. - Can be used as an index into an array */ + /* Total of the weights of cases matching this interaction */ + double cc; }; +static struct value_node * +lookup_value (const struct hmap *map, const union value *val, unsigned int hash, int width) +{ + struct value_node *vn = NULL; + HMAP_FOR_EACH_WITH_HASH (vn, struct value_node, node, hash, map) + { + if (value_equal (&vn->val, val, width)) + break; + } + + return vn; +} + +struct variable_node +{ + struct hmap_node node; /* Node in hash map. */ + const struct variable *var; /* The variable */ + + struct hmap valmap; /* A map of value nodes */ +}; + +static void +dump_interaction (const struct interaction *iact) +{ + struct string str = DS_EMPTY_INITIALIZER; + interaction_to_string (iact, &str); + printf ("Interaction: %s\n", ds_cstr (&str)); + ds_destroy (&str); +} + + +static struct variable_node * +lookup_variable (const struct hmap *map, const struct variable *var, unsigned int hash) +{ + struct variable_node *vn = NULL; + HMAP_FOR_EACH_WITH_HASH (vn, struct variable_node, node, hash, map) + { + if (vn->var == var) + break; + + fprintf (stderr, "Warning: Hash table collision\n"); + } + + return vn; +} + + struct interact_params { - /* A map indexed by a union values */ - struct hmap map; + /* A map indexed by a interaction_value */ + struct hmap ivmap; const struct interaction *iact; int base_subscript_short; int base_subscript_long; - /* The number of distinct values of this variable */ + /* The number of distinct values of this interaction */ int n_cats; - /* A map of values indexed by subscript */ - const struct value_node **reverse_value_map; + /* The degrees of freedom for this interaction */ + int df; - /* Total of the weights of this variable */ - double cc; + /* A map of interaction_values indexed by subscript */ + struct interaction_value **reverse_interaction_value_map; + + double cc; }; +#if 0 /* Comparison function to sort the reverse_value_map in ascending order */ static int compare_value_node (const void *vn1_, const void *vn2_, const void *aux) @@ -74,7 +132,7 @@ compare_value_node (const void *vn1_, const void *vn2_, const void *aux) return interaction_case_cmp_3way (vp->iact, (*vn1)->ccase, (*vn2)->ccase); } - +#endif struct categoricals { @@ -84,6 +142,9 @@ struct categoricals /* An array of interact_params */ struct interact_params *iap; + /* Map whose members are the union of the variables which comprise IAP */ + struct hmap varmap; + /* The size of IAP. (ie, the number of interactions involved.) */ size_t n_iap; @@ -91,6 +152,8 @@ struct categoricals In the absence of missing values, this will be equal to N_IAP */ size_t n_vars; + size_t df_sum; + /* A map to enable the lookup of variables indexed by subscript. This map considers only the N - 1 of the N variables. */ @@ -118,113 +181,81 @@ struct categoricals }; -void -categoricals_destroy ( struct categoricals *cat) +static void +categoricals_dump (const struct categoricals *cat) { int i; - if (cat != NULL) +#if 1 + printf ("Reverse Variable Map (short):\n"); + for (i = 0; i < cat->df_sum; ++i) { - for (i = 0 ; i < cat->n_iap; ++i) - { - struct hmap *map = &cat->iap[i].map; - struct value_node *nn; - - HMAP_FOR_EACH (nn, struct value_node, node, map) - { - case_unref (nn->ccase); - } - - hmap_destroy (map); - } - - pool_destroy (cat->pool); - free (cat); + printf (" %d", cat->reverse_variable_map_short[i]); } -} - - -#if 0 -void -categoricals_dump (const struct categoricals *cat) -{ - int v; + printf ("\n"); - for (v = 0 ; v < cat->n_iap; ++v) + printf ("Number of interactions %d\n", cat->n_iap); + for (i = 0 ; i < cat->n_iap; ++i) { - const struct interact_params *vp = &cat->iap[v]; - const struct hmap *m = &vp->map; - struct hmap_node *node ; - int x; - struct string str ; + int v; + struct string str; + const struct interact_params *iap = &cat->iap[i]; + const struct interaction *iact = iap->iact; + ds_init_empty (&str); + interaction_to_string (iact, &str); - interaction_to_string (vp->iact, &str); - printf ("\n%s (%d) CC=%g n_cats=%d:\n", - ds_cstr (&str), - vp->base_subscript_long, vp->cc, vp->n_cats); + printf ("\nInteraction: %s (n: %d; df: %d ); ", ds_cstr (&str), iap->n_cats, iap->df); + ds_destroy (&str); + printf ("Base subscript: %d\n", iap->base_subscript_short); -#if 0 - printf ("Reverse map\n"); - for (x = 0 ; x < vp->n_cats; ++x) + printf ("\t("); + for (v = 0; v < hmap_count (&iap->ivmap); ++v) { - struct string s; - const struct value_node *vn = vp->reverse_value_map[x]; - ds_init_empty (&s); - var_append_value_name (vp->var, &vn->value, &s); - printf ("Value for %d is %s\n", x, ds_cstr(&s)); - ds_destroy (&s); - } - - printf ("\nForward map\n"); - for (node = hmap_first (m); node; node = hmap_next (m, node)) - { - struct string s; - const struct value_node *vn = HMAP_DATA (node, struct value_node, node); - ds_init_empty (&s); - var_append_value_name (vp->var, &vn->value, &s); - printf ("Value: %s; Index %d; CC %g\n", - ds_cstr (&s), - vn->subscript, vn->cc); - ds_destroy (&s); + int vv; + const struct interaction_value *iv = iap->reverse_interaction_value_map[v]; + + if (v > 0) printf (" "); + printf ("{"); + for (vv = 0; vv < iact->n_vars; ++vv) + { + const struct variable *var = iact->vars[vv]; + const union value *val = case_data (iv->ccase, var); + + printf ("%g", val->f); + if (vv < iact->n_vars - 1) + printf (", "); + } + printf ("}"); } -#endif + printf (")\n"); } +#endif +} - assert (cat->n_vars <= cat->n_iap); - - printf ("\n"); - printf ("Number of interactions: %d\n", cat->n_iap); - printf ("Number of non-empty categorical variables: %d\n", cat->n_vars); - printf ("Total number of categories: %d\n", cat->n_cats_total); - - printf ("\nReverse variable map (short):\n"); - for (v = 0 ; v < cat->n_cats_total - cat->n_vars; ++v) - printf ("%d ", cat->reverse_variable_map_short[v]); - - printf ("\nReverse variable map (long):\n"); - for (v = 0 ; v < cat->n_cats_total; ++v) - printf ("%d ", cat->reverse_variable_map_long[v]); - printf ("\n"); +void +categoricals_destroy (struct categoricals *cat) +{ + free (cat); } -#endif -static struct value_node * + +static struct interaction_value * lookup_case (const struct hmap *map, const struct interaction *iact, const struct ccase *c) { - struct value_node *nn; + struct interaction_value *iv = NULL; size_t hash = interaction_case_hash (iact, c); - HMAP_FOR_EACH_WITH_HASH (nn, struct value_node, node, hash, map) + HMAP_FOR_EACH_WITH_HASH (iv, struct interaction_value, node, hash, map) { - if (interaction_case_equal (iact, c, nn->ccase)) + if (interaction_case_equal (iact, c, iv->ccase)) break; fprintf (stderr, "Warning: Hash table collision\n"); } - return nn; + return iv; } @@ -252,13 +283,29 @@ categoricals_create (const struct interaction **inter, size_t n_inter, cat->aux1 = aux1; cat->aux2 = aux2; - cat->iap = pool_calloc (cat->pool, cat->n_iap, sizeof *cat->iap); + hmap_init (&cat->varmap); for (i = 0 ; i < cat->n_iap; ++i) { - hmap_init (&cat->iap[i].map); + int v; + hmap_init (&cat->iap[i].ivmap); cat->iap[i].iact = inter[i]; + cat->iap[i].cc = 0.0; + for (v = 0; v < inter[i]->n_vars; ++v) + { + const struct variable *var = inter[i]->vars[v]; + unsigned int hash = hash_pointer (var, 0); + struct variable_node *vn = lookup_variable (&cat->varmap, var, hash); + if (vn == NULL) + { + vn = pool_malloc (cat->pool, sizeof *vn); + vn->var = var; + hmap_init (&vn->valmap); + + hmap_insert (&cat->varmap, &vn->node, hash); + } + } } return cat; @@ -269,67 +316,93 @@ categoricals_create (const struct interaction **inter, size_t n_inter, void categoricals_update (struct categoricals *cat, const struct ccase *c) { - size_t i; - + int i; + struct variable_node *vn = NULL; const double weight = cat->wv ? case_data (c, cat->wv)->f : 1.0; assert (NULL == cat->reverse_variable_map_short); assert (NULL == cat->reverse_variable_map_long); + /* Interate over each variable, and add the value of that variable + to the appropriate map, if it's not already present. */ + HMAP_FOR_EACH (vn, struct variable_node, node, &cat->varmap) + { + const int width = var_get_width (vn->var); + const union value *val = case_data (c, vn->var); + unsigned int hash = value_hash (val, width, 0); + + struct value_node *valn = lookup_value (&vn->valmap, val, hash, width); + if (valn == NULL) + { + valn = pool_malloc (cat->pool, sizeof *valn); + value_init (&valn->val, width); + value_copy (&valn->val, val, width); + hmap_insert (&vn->valmap, &valn->node, hash); + } + } + + for (i = 0 ; i < cat->n_iap; ++i) { const struct interaction *iact = cat->iap[i].iact; - size_t hash; - struct value_node *node ; - if ( interaction_case_is_missing (iact, c, cat->exclude)) - continue; + // if ( interaction_case_is_missing (iact, c, cat->exclude)) + // continue; - hash = interaction_case_hash (iact, c); - node = lookup_case (&cat->iap[i].map, iact, c); + size_t hash = interaction_case_hash (iact, c); + struct interaction_value *node = lookup_case (&cat->iap[i].ivmap, iact, c); if ( NULL == node) { node = pool_malloc (cat->pool, sizeof *node); node->ccase = case_ref (c); - node->cc = 0.0; - - hmap_insert (&cat->iap[i].map, &node->node, hash); - cat->n_cats_total++; - - if ( 0 == cat->iap[i].n_cats) - cat->n_vars++; + node->cc = weight; - node->subscript = cat->iap[i].n_cats++ ; + hmap_insert (&cat->iap[i].ivmap, &node->node, hash); - if (cat->user_data_create) - node->user_data = cat->user_data_create (cat->aux1, cat->aux2); + // if (cat->user_data_create) + // node->user_data = cat->user_data_create (cat->aux1, cat->aux2); + } + else + { + node->cc += weight; } - - node->cc += weight; cat->iap[i].cc += weight; - if (cat->update) - cat->update (node->user_data, cat->exclude, cat->wv, NULL, c, cat->aux1, cat->aux2); + // if (cat->update) + // cat->update (node->user_data, cat->exclude, cat->wv, NULL, c, cat->aux1, cat->aux2); + } } -/* Return the number of categories (distinct values) for variable N */ +/* Return the number of categories (distinct values) for interction N */ size_t categoricals_n_count (const struct categoricals *cat, size_t n) { - return hmap_count (&cat->iap[n].map); + return hmap_count (&cat->iap[n].ivmap); +} + + +size_t +categoricals_df (const struct categoricals *cat, size_t n) +{ + return cat->iap[n].df; } /* Return the total number of categories */ size_t -categoricals_total (const struct categoricals *cat) +categoricals_n_total (const struct categoricals *cat) { return cat->n_cats_total; } +size_t +categoricals_df_total (const struct categoricals *cat) +{ + return cat->df_sum; +} /* This function must be called *before* any call to categoricals_get_*_by subscript and *after* all calls to categoricals_update */ @@ -344,46 +417,78 @@ categoricals_done (const struct categoricals *cat_) */ struct categoricals *cat = CONST_CAST (struct categoricals *, cat_); int v; + int i; int idx_short = 0; int idx_long = 0; + cat->df_sum = 0; + cat->n_cats_total = 0; + + /* Calculate the degrees of freedom, and the number of categories */ + for (i = 0 ; i < cat->n_iap; ++i) + { + const struct interaction *iact = cat->iap[i].iact; + + cat->iap[i].df = 1; + cat->iap[i].n_cats = 1; + + for (v = 0 ; v < iact->n_vars; ++v) + { + const struct variable *var = iact->vars[v]; + + struct variable_node *vn = lookup_variable (&cat->varmap, var, hash_pointer (var, 0)); + + cat->iap[i].df *= hmap_count (&vn->valmap) - 1; + cat->iap[i].n_cats *= hmap_count (&vn->valmap); + } + + cat->df_sum += cat->iap[i].df; + cat->n_cats_total += cat->iap[i].n_cats; + } + + cat->reverse_variable_map_short = pool_calloc (cat->pool, - cat->n_cats_total - cat->n_vars, + cat->df_sum, sizeof *cat->reverse_variable_map_short); cat->reverse_variable_map_long = pool_calloc (cat->pool, cat->n_cats_total, sizeof *cat->reverse_variable_map_long); - - for (v = 0 ; v < cat->n_iap; ++v) + + for (i = 0 ; i < cat->n_iap; ++i) { - int i; - struct interact_params *vp = &cat->iap[v]; - int n_cats_total = categoricals_n_count (cat, v); - struct hmap_node *node ; + struct interaction_value *ivn = NULL; + int x = 0; + int ii; + struct interact_params *iap = &cat->iap[i]; + + iap->base_subscript_short = idx_short; + iap->base_subscript_long = idx_long; - vp->reverse_value_map = pool_calloc (cat->pool, n_cats_total, sizeof *vp->reverse_value_map); + iap->reverse_interaction_value_map = pool_calloc (cat->pool, iap->n_cats, + sizeof *iap->reverse_interaction_value_map); - vp->base_subscript_short = idx_short; - vp->base_subscript_long = idx_long; - for (node = hmap_first (&vp->map); node; node = hmap_next (&vp->map, node)) + HMAP_FOR_EACH (ivn, struct interaction_value, node, &iap->ivmap) { - const struct value_node *vn = HMAP_DATA (node, struct value_node, node); - vp->reverse_value_map[vn->subscript] = vn; + iap->reverse_interaction_value_map[x++] = ivn; } +#if 0 /* For some purposes (eg CONTRASTS in ONEWAY) the values need to be sorted */ - sort (vp->reverse_value_map, vp->n_cats, sizeof (const struct value_node *), + sort (vp->reverse_interaction_value_map, vp->n_cats, sizeof (const struct interaction_value *), compare_value_node, vp); +#endif /* Populate the reverse variable maps. */ - for (i = 0; i < vp->n_cats - 1; ++i) - cat->reverse_variable_map_short[idx_short++] = v; + for (ii = 0; ii < iap->df; ++ii) + cat->reverse_variable_map_short[idx_short++] = i; - for (i = 0; i < vp->n_cats; ++i) - cat->reverse_variable_map_long[idx_long++] = v; + for (ii = 0; ii < iap->n_cats; ++ii) + cat->reverse_variable_map_long[idx_long++] = i; } + + assert (cat->n_vars <= cat->n_iap); } @@ -393,7 +498,7 @@ reverse_variable_lookup_short (const struct categoricals *cat, int subscript) { assert (cat->reverse_variable_map_short); assert (subscript >= 0); - assert (subscript < cat->n_cats_total - cat->n_vars); + assert (subscript < cat->df_sum); return cat->reverse_variable_map_short[subscript]; } @@ -401,6 +506,7 @@ reverse_variable_lookup_short (const struct categoricals *cat, int subscript) static int reverse_variable_lookup_long (const struct categoricals *cat, int subscript) { + printf ("%s\n", __FUNCTION__); assert (cat->reverse_variable_map_long); assert (subscript >= 0); assert (subscript < cat->n_cats_total); @@ -425,8 +531,7 @@ categoricals_get_case_by_subscript (const struct categoricals *cat, int subscrip { int vindex = reverse_variable_lookup_short (cat, subscript); const struct interact_params *vp = &cat->iap[vindex]; - const struct value_node *vn = vp->reverse_value_map [subscript - vp->base_subscript_short]; - + const struct interaction_value *vn = vp->reverse_interaction_value_map [subscript - vp->base_subscript_short]; return vn->ccase; } @@ -446,10 +551,10 @@ categoricals_get_sum_by_subscript (const struct categoricals *cat, int subscript int vindex = reverse_variable_lookup_short (cat, subscript); const struct interact_params *vp = &cat->iap[vindex]; - const struct value_node *vn = vp->reverse_value_map [subscript - vp->base_subscript_short]; - return vn->cc; -} + const struct interaction_value *iv = vp->reverse_interaction_value_map [subscript - vp->base_subscript_short]; + return iv->cc; +} /* Returns unity if the value in case C at SUBSCRIPT is equal to the category for that subscript */ @@ -468,28 +573,32 @@ categoricals_get_binary_by_subscript (const struct categoricals *cat, int subscr size_t categoricals_get_n_variables (const struct categoricals *cat) { + printf ("%s\n", __FUNCTION__); return cat->n_vars; } + /* Return a case containing the set of values corresponding to SUBSCRIPT */ const struct ccase * categoricals_get_case_by_category (const struct categoricals *cat, int subscript) { int vindex = reverse_variable_lookup_long (cat, subscript); const struct interact_params *vp = &cat->iap[vindex]; - const struct value_node *vn = vp->reverse_value_map [subscript - vp->base_subscript_long]; + // const struct interaction_value *vn = vp->reverse_interaction_value_map [subscript - vp->base_subscript_long]; - return vn->ccase; + // return vn->ccase; + return NULL; } - void * categoricals_get_user_data_by_category (const struct categoricals *cat, int subscript) { int vindex = reverse_variable_lookup_long (cat, subscript); const struct interact_params *vp = &cat->iap[vindex]; - const struct value_node *vn = vp->reverse_value_map [subscript - vp->base_subscript_long]; - return vn->user_data; + //const struct value_node *vn = vp->reverse_interaction_value_map [subscript - vp->base_subscript_long]; + //return vn->user_data; + assert (0); + return NULL; } diff --git a/src/math/categoricals.h b/src/math/categoricals.h index d9ac1f5e..ef9c37c2 100644 --- a/src/math/categoricals.h +++ b/src/math/categoricals.h @@ -50,9 +50,14 @@ void categoricals_update (struct categoricals *cat, const struct ccase *c); /* Return the number of categories (distinct values) for variable N */ size_t categoricals_n_count (const struct categoricals *cat, size_t n); +size_t categoricals_df (const struct categoricals *cat, size_t n); /* Return the total number of categories */ -size_t categoricals_total (const struct categoricals *cat); +size_t categoricals_n_total (const struct categoricals *cat); + +/* Return the total degrees of freedom */ +size_t categoricals_df_total (const struct categoricals *cat); + /* Return the total number of variables which participated in these categoricals. diff --git a/src/math/covariance.c b/src/math/covariance.c index 1a393351..b37da3d2 100644 --- a/src/math/covariance.c +++ b/src/math/covariance.c @@ -350,8 +350,7 @@ covariance_accumulate_pass2 (struct covariance *cov, const struct ccase *c) cov->dim = cov->n_vars; if (cov->categoricals) - cov->dim += categoricals_total (cov->categoricals) - - categoricals_get_n_variables (cov->categoricals); + cov->dim += categoricals_df_total (cov->categoricals); cov->n_cm = (cov->dim * (cov->dim - 1) ) / 2; cov->cm = xcalloc (sizeof *cov->cm, cov->n_cm); @@ -363,8 +362,8 @@ covariance_accumulate_pass2 (struct covariance *cov, const struct ccase *c) cov->moments[i] = resize_matrix (cov->moments[i], cov->dim); } - if (cov->categoricals) - categoricals_done (cov->categoricals); + // if (cov->categoricals) + // categoricals_done (cov->categoricals); /* Populate the moments matrices with the categorical value elements */ for (i = cov->n_vars; i < cov->dim; ++i)