#include <math/design-matrix.h>
#include <math/coefficient.h>
#include <math/linreg/linreg.h>
+#include <math/moments.h>
#include <output/table.h>
#include "gettext.h"
/* (functions) */
static struct cmd_regression cmd;
+/*
+ Moments for each of the variables used.
+ */
+struct moments_var
+{
+ struct moments1 *m;
+ struct variable *v;
+};
+
/* Linear regression models. */
static pspp_linreg_cache **models = NULL;
/*
P values for the test statistic above.
*/
- pval = 2 * gsl_cdf_tdist_Q (fabs (t_stat), 1.0);
+ pval = 2 * gsl_cdf_tdist_Q (fabs (t_stat), (double) (c->n_obs - c->n_coeffs));
tab_float (t, 6, j + 1, 0, pval, 10, 2);
}
tab_title (t, _("Coefficients"));
for (i = 0; i < n_vars; i++)
{
- size_t n_categories = cat_get_n_categories (varlist[i]);
- size_t j;
+ int n_categories = cat_get_n_categories (varlist[i]);
+ int j;
fprintf (fp, "%s.name = \"%s\";\n\t",
var_get_name (varlist[i]),
/*
Mark missing cases. Return the number of non-missing cases.
+ Compute the first two moments.
*/
static size_t
mark_missing_cases (const struct casefile *cf, struct variable *v,
- int *is_missing_case, double n_data)
+ int *is_missing_case, double n_data, struct moments_var *mom)
{
struct casereader *r;
struct ccase c;
size_t row;
const union value *val;
+ double w = 1.0;
for (r = casefile_get_reader (cf, NULL);
casereader_read (r, &c); case_destroy (&c))
row = casereader_cnum (r) - 1;
val = case_data (&c, v);
+ if (mom != NULL)
+ {
+ moments1_add (mom->m, val->f, w);
+ }
cat_value_update (v, val);
- if (var_is_value_missing (v, val))
+ if (var_is_value_missing (v, val, MV_ANY))
{
if (!is_missing_case[row])
{
static int
prepare_data (int n_data, int is_missing_case[],
struct variable **indep_vars,
- struct variable *depvar, const struct casefile *cf)
+ struct variable *depvar, const struct casefile *cf, struct moments_var *mom)
{
int i;
int j;
cat_stored_values_create (v_variables[i]);
}
n_data =
- mark_missing_cases (cf, v_variables[i], is_missing_case, n_data);
+ mark_missing_cases (cf, v_variables[i], is_missing_case, n_data, mom + i);
}
}
/*
Mark missing cases for the dependent variable.
*/
- n_data = mark_missing_cases (cf, depvar, is_missing_case, n_data);
+ n_data = mark_missing_cases (cf, depvar, is_missing_case, n_data, NULL);
return n_data;
}
pspp_coeff_init (c->coeff + 1, dm);
}
+/*
+ Put the moments in the linreg cache.
+ */
+static void
+compute_moments (pspp_linreg_cache *c, struct moments_var *mom, struct design_matrix *dm, size_t n)
+{
+ size_t i;
+ size_t j;
+ double weight;
+ double mean;
+ double variance;
+ double skewness;
+ double kurtosis;
+ /*
+ Scan the variable names in the columns of the design matrix.
+ When we find the variable we need, insert its mean in the cache.
+ */
+ for (i = 0; i < dm->m->size2; i++)
+ {
+ for (j = 0; j < n; j++)
+ {
+ if (design_matrix_col_to_var (dm, i) == (mom + j)->v)
+ {
+ moments1_calculate ((mom + j)->m, &weight, &mean, &variance,
+ &skewness, &kurtosis);
+ gsl_vector_set (c->indep_means, i, mean);
+ gsl_vector_set (c->indep_std, i, sqrt (variance));
+ }
+ }
+ }
+}
static bool
run_regression (const struct ccase *first,
const struct casefile *cf, void *cmd_ UNUSED, const struct dataset *ds)
struct ccase c;
struct variable **indep_vars;
struct design_matrix *X;
+ struct moments_var *mom;
gsl_vector *Y;
pspp_linreg_opts lopts;
}
is_missing_case = xnmalloc (n_cases, sizeof (*is_missing_case));
-
+ mom = xnmalloc (n_variables, sizeof (*mom));
+ for (i = 0; i < n_variables; i++)
+ {
+ (mom + i)->m = moments1_create (MOMENT_VARIANCE);
+ (mom + i)->v = v_variables[i];
+ }
lopts.get_depvar_mean_std = 1;
for (k = 0; k < cmd.n_dependent; k++)
}
n_data = prepare_data (n_cases, is_missing_case, indep_vars,
cmd.v_dependent[k],
- (const struct casefile *) cf);
+ (const struct casefile *) cf, mom);
Y = gsl_vector_alloc (n_data);
X =
Find the least-squares estimates and other statistics.
*/
pspp_linreg ((const gsl_vector *) Y, X->m, &lopts, models[k]);
+ compute_moments (models[k], mom, X, n_variables);
subcommand_statistics (cmd.a_statistics, models[k]);
subcommand_export (cmd.sbc_export, models[k]);
free (lopts.get_indep_mean_std);
casereader_destroy (r);
}
-
+ for (i = 0; i < n_variables; i++)
+ {
+ moments1_destroy ((mom + i)->m);
+ }
+ free (mom);
free (is_missing_case);
return true;