/* PSPP - a program for statistical analysis.
- Copyright (C) 2005 Free Software Foundation, Inc.
+ Copyright (C) 2005, 2009 Free Software Foundation, Inc.
This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
#include <libpspp/taint.h>
#include <math/design-matrix.h>
#include <math/coefficient.h>
-#include <math/linreg/linreg.h>
+#include <math/linreg.h>
#include <math/moments.h>
#include <output/table.h>
assert (c != NULL);
rsq = c->ssm / c->sst;
adjrsq = 1.0 - (1.0 - rsq) * (c->n_obs - 1.0) / (c->n_obs - c->n_indeps);
- std_error = sqrt ((c->n_indeps - 1.0) / (c->n_obs - 1.0));
+ std_error = sqrt (pspp_linreg_mse (c));
t = tab_create (n_cols, n_rows, 0);
tab_dim (t, tab_natural_dimensions);
tab_box (t, TAL_2, TAL_2, -1, TAL_1, 0, 0, n_cols - 1, n_rows - 1);
int this_row;
double t_stat;
double pval;
- double coeff;
double std_err;
double beta;
const char *label;
tab_float (t, 2, 1, 0, c->intercept, 10, 2);
std_err = sqrt (gsl_matrix_get (c->cov, 0, 0));
tab_float (t, 3, 1, 0, std_err, 10, 2);
- beta = c->intercept / c->depvar_std;
- tab_float (t, 4, 1, 0, beta, 10, 2);
+ tab_float (t, 4, 1, 0, 0.0, 10, 2);
t_stat = c->intercept / std_err;
tab_float (t, 5, 1, 0, t_stat, 10, 2);
pval = 2 * gsl_cdf_tdist_Q (fabs (t_stat), 1.0);
/*
Regression coefficients.
*/
- coeff = c->coeff[j]->estimate;
- tab_float (t, 2, this_row, 0, coeff, 10, 2);
+ tab_float (t, 2, this_row, 0, c->coeff[j]->estimate, 10, 2);
/*
Standard error of the coefficients.
*/
std_err = sqrt (gsl_matrix_get (c->cov, j + 1, j + 1));
tab_float (t, 3, this_row, 0, std_err, 10, 2);
/*
- 'Standardized' coefficient, i.e., regression coefficient
+ Standardized coefficient, i.e., regression coefficient
if all variables had unit variance.
*/
- beta = gsl_vector_get (c->indep_std, j + 1);
- beta *= coeff / c->depvar_std;
+ beta = pspp_coeff_get_sd (c->coeff[j]);
+ beta *= c->coeff[j]->estimate / c->depvar_std;
tab_float (t, 4, this_row, 0, beta, 10, 2);
/*
Test statistic for H0: coefficient is 0.
*/
- t_stat = coeff / std_err;
+ t_stat = c->coeff[j]->estimate / std_err;
tab_float (t, 5, this_row, 0, t_stat, 10, 2);
/*
P values for the test statistic above.
int n_cols = 7;
int n_rows = 4;
const double msm = c->ssm / c->dfm;
- const double mse = c->sse / c->dfe;
+ const double mse = pspp_linreg_mse (c);
const double F = msm / mse;
const double pval = gsl_cdf_fdist_Q (F, c->dfm, c->dfe);
/* Degrees of freedom */
- tab_float (t, 3, 1, 0, c->dfm, 4, 0);
- tab_float (t, 3, 2, 0, c->dfe, 4, 0);
- tab_float (t, 3, 3, 0, c->dft, 4, 0);
+ tab_text (t, 3, 1, TAB_RIGHT | TAT_PRINTF, "%g", c->dfm);
+ tab_text (t, 3, 2, TAB_RIGHT | TAT_PRINTF, "%g", c->dfe);
+ tab_text (t, 3, 3, TAB_RIGHT | TAT_PRINTF, "%g", c->dft);
/* Mean Squares */
-
tab_float (t, 4, 1, TAB_RIGHT, msm, 8, 3);
tab_float (t, 4, 2, TAB_RIGHT, mse, 8, 3);
tab_title (t, _("ANOVA"));
tab_submit (t);
}
+
static void
reg_stats_outs (pspp_linreg_cache * c)
{
assert (c != NULL);
}
+
static void
reg_stats_zpp (pspp_linreg_cache * c)
{
assert (c != NULL);
}
+
static void
reg_stats_label (pspp_linreg_cache * c)
{
assert (c != NULL);
}
+
static void
reg_stats_sha (pspp_linreg_cache * c)
{
Gets the predicted values.
*/
static int
-regression_trns_pred_proc (void *t_, struct ccase *c,
+regression_trns_pred_proc (void *t_, struct ccase **c,
casenumber case_idx UNUSED)
{
size_t i;
n_vals = (*model->get_vars) (model, vars);
vals = xnmalloc (n_vals, sizeof (*vals));
- output = case_data_rw (c, model->pred);
- assert (output != NULL);
+ *c = case_unshare (*c);
+ output = case_data_rw (*c, model->pred);
for (i = 0; i < n_vals; i++)
{
- vals[i] = case_data (c, vars[i]);
+ vals[i] = case_data (*c, vars[i]);
}
output->f = (*model->predict) ((const struct variable **) vars,
vals, model, n_vals);
Gets the residuals.
*/
static int
-regression_trns_resid_proc (void *t_, struct ccase *c,
+regression_trns_resid_proc (void *t_, struct ccase **c,
casenumber case_idx UNUSED)
{
size_t i;
n_vals = (*model->get_vars) (model, vars);
vals = xnmalloc (n_vals, sizeof (*vals));
- output = case_data_rw (c, model->resid);
+ *c = case_unshare (*c);
+ output = case_data_rw (*c, model->resid);
assert (output != NULL);
for (i = 0; i < n_vals; i++)
{
- vals[i] = case_data (c, vars[i]);
+ vals[i] = case_data (*c, vars[i]);
}
- obs = case_data (c, model->depvar);
+ obs = case_data (*c, model->depvar);
output->f = (*model->residual) ((const struct variable **) vars,
vals, obs, model, n_vals);
free (vals);
for (lc = models; lc < models + cmd.n_dependent; lc++)
{
- assert (*lc != NULL);
- assert ((*lc)->depvar != NULL);
- if (cmd.a_save[REGRESSION_SV_RESID])
- {
- reg_save_var (ds, "RES", regression_trns_resid_proc, *lc,
- &(*lc)->resid, n_trns);
- }
- if (cmd.a_save[REGRESSION_SV_PRED])
+ if (*lc != NULL)
{
- reg_save_var (ds, "PRED", regression_trns_pred_proc, *lc,
- &(*lc)->pred, n_trns);
+ if ((*lc)->depvar != NULL)
+ {
+ if (cmd.a_save[REGRESSION_SV_RESID])
+ {
+ reg_save_var (ds, "RES", regression_trns_resid_proc, *lc,
+ &(*lc)->resid, n_trns);
+ }
+ if (cmd.a_save[REGRESSION_SV_PRED])
+ {
+ reg_save_var (ds, "PRED", regression_trns_pred_proc, *lc,
+ &(*lc)->pred, n_trns);
+ }
+ }
}
}
}
for (i = 0; i < n_variables; i++)
if (!is_depvar (i, depvar))
indep_vars[n_indep_vars++] = v_variables[i];
- if ((n_indep_vars < 2) && is_depvar (0, depvar))
+ if ((n_indep_vars < 1) && is_depvar (0, depvar))
{
/*
There is only one independent variable, and it is the same
struct moments_var *mom)
{
int n_data;
- struct ccase c;
+ struct ccase *c;
size_t i;
assert (vars != NULL);
cat_stored_values_create (vars[i]);
n_data = 0;
- for (; casereader_read (input, &c); case_destroy (&c))
+ for (; (c = casereader_read (input)) != NULL; case_unref (c))
{
/*
The second condition ensures the program will run even if
*/
for (i = 0; i < n_vars; i++)
{
- const union value *val = case_data (&c, vars[i]);
+ const union value *val = case_data (c, vars[i]);
if (var_is_alpha (vars[i]))
cat_value_update (vars[i], val);
else
pspp_coeff_init (c->coeff, 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 (struct casereader *input, struct cmd_regression *cmd,
struct dataset *ds, pspp_linreg_cache **models)
size_t i;
int n_indep = 0;
int k;
- struct ccase c;
+ struct ccase *c;
const struct variable **indep_vars;
struct design_matrix *X;
struct moments_var *mom;
assert (models != NULL);
- if (!casereader_peek (input, 0, &c))
+ c = casereader_peek (input, 0);
+ if (c == NULL)
{
casereader_destroy (input);
return true;
}
- output_split_file_values (ds, &c);
- case_destroy (&c);
+ output_split_file_values (ds, c);
+ case_unref (c);
if (!v_variables)
{
const struct variable *dep_var;
struct casereader *reader;
casenumber row;
- struct ccase c;
+ struct ccase *c;
size_t n_data; /* Number of valid cases. */
dep_var = cmd->v_dependent[k];
n_indep = identify_indep_vars (indep_vars, dep_var);
reader = casereader_clone (input);
reader = casereader_create_filter_missing (reader, indep_vars, n_indep,
- MV_ANY, NULL);
+ MV_ANY, NULL, NULL);
reader = casereader_create_filter_missing (reader, &dep_var, 1,
- MV_ANY, NULL);
+ MV_ANY, NULL, NULL);
n_data = prepare_categories (casereader_clone (reader),
indep_vars, n_indep, mom);
{
lopts.get_indep_mean_std[i] = 1;
}
- models[k] = pspp_linreg_cache_alloc (X->m->size1, X->m->size2);
+ models[k] = pspp_linreg_cache_alloc (dep_var, (const struct variable **) indep_vars,
+ X->m->size1, X->m->size2);
models[k]->depvar = dep_var;
/*
For large data sets, use QR decomposition.
The second pass fills the design matrix.
*/
reader = casereader_create_counter (reader, &row, -1);
- for (; casereader_read (reader, &c); case_destroy (&c))
+ for (; (c = casereader_read (reader)) != NULL; case_unref (c))
{
for (i = 0; i < n_indep; ++i)
{
const struct variable *v = indep_vars[i];
- const union value *val = case_data (&c, v);
+ const union value *val = case_data (c, v);
if (var_is_alpha (v))
design_matrix_set_categorical (X, row, v, val);
else
design_matrix_set_numeric (X, row, v, val);
}
- gsl_vector_set (Y, row, case_num (&c, dep_var));
+ gsl_vector_set (Y, row, case_num (c, dep_var));
}
/*
Now that we know the number of coefficients, allocate space
/*
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);
+ pspp_linreg ((const gsl_vector *) Y, X, &lopts, models[k]);
if (!taint_has_tainted_successor (casereader_get_taint (input)))
{