1 /* PSPP - a program for statistical analysis.
2 Copyright (C) 2011 Free Software Foundation, Inc.
3
4 This program is free software: you can redistribute it and/or modify
5 it under the terms of the GNU General Public License as published by
6 the Free Software Foundation, either version 3 of the License, or
7 (at your option) any later version.
8
9 This program is distributed in the hope that it will be useful,
10 but WITHOUT ANY WARRANTY; without even the implied warranty of
11 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 GNU General Public License for more details.
13
14 You should have received a copy of the GNU General Public License
15 along with this program. If not, see <http://www.gnu.org/licenses/>.
16 */
17
18 #include <config.h>
19
20
21 #include "t-test.h"
22
23 #include <math.h>
24 #include <gsl/gsl_cdf.h>
25
26 #include "data/variable.h"
27 #include "data/format.h"
MCObjectStreamer(MCContext & Context,std::unique_ptr<MCAsmBackend> TAB,std::unique_ptr<MCObjectWriter> OW,std::unique_ptr<MCCodeEmitter> Emitter)28 #include "data/casereader.h"
29 #include "data/dictionary.h"
30 #include "libpspp/hash-functions.h"
31 #include "libpspp/hmapx.h"
32 #include "math/moments.h"
33
34 #include "output/pivot-table.h"
35
36 #include "gettext.h"
37 #define N_(msgid) msgid
38 #define _(msgid) gettext (msgid)
39
40
41 struct per_var_stats
42 {
43 const struct variable *var;
44
45 /* N, Mean, Variance */
46 struct moments *mom;
47
48 /* Sum of the differences */
49 double sum_diff;
50 };
addPendingLabel(MCSymbol * S)51
52
53 struct one_samp
54 {
55 struct per_var_stats *stats;
56 size_t n_stats;
57 double testval;
58 };
59
60
61 static void
62 one_sample_test (const struct tt *tt, const struct one_samp *os)
63 {
64 struct pivot_table *table = pivot_table_create (N_("One-Sample Test"));
65 pivot_table_set_weight_var (table, tt->wv);
66
67 struct pivot_dimension *statistics = pivot_dimension_create (
68 table, PIVOT_AXIS_COLUMN, N_("Statistics"));
69 struct pivot_category *group = pivot_category_create_group__ (
70 statistics->root, pivot_value_new_user_text_nocopy (
flushPendingLabels(MCFragment * F,uint64_t FOffset)71 xasprintf (_("Test Value = %.*g"), DBL_DIG + 1, os->testval)));
72 pivot_category_create_leaves (
73 group,
74 N_("t"), PIVOT_RC_OTHER,
75 N_("df"), PIVOT_RC_COUNT,
76 N_("Sig. (2-tailed)"), PIVOT_RC_SIGNIFICANCE,
77 N_("Mean Difference"), PIVOT_RC_OTHER);
78 struct pivot_category *subgroup = pivot_category_create_group__ (
79 group, pivot_value_new_user_text_nocopy (
80 xasprintf (_("%g%% Confidence Interval of the Difference"),
81 tt->confidence * 100.0)));
82 pivot_category_create_leaves (subgroup,
83 N_("Lower"), PIVOT_RC_OTHER,
84 N_("Upper"), PIVOT_RC_OTHER);
85
86 struct pivot_dimension *dep_vars = pivot_dimension_create (
87 table, PIVOT_AXIS_ROW, N_("Dependent Variables"));
88
89 for (size_t i = 0; i < os->n_stats; i++)
90 {
91 const struct per_var_stats *per_var_stats = &os->stats[i];
flushPendingLabels()92 const struct moments *m = per_var_stats->mom;
93
94 int dep_var_idx = pivot_category_create_leaf (
95 dep_vars->root, pivot_value_new_variable (per_var_stats->var));
96
97 double cc, mean, sigma;
98 moments_calculate (m, &cc, &mean, &sigma, NULL, NULL);
99 double tval = (mean - os->testval) * sqrt (cc / sigma);
100 double mean_diff = per_var_stats->sum_diff / cc;
101 double se_mean = sqrt (sigma / cc);
102 double df = cc - 1.0;
103 double p = gsl_cdf_tdist_P (tval, df);
104 double q = gsl_cdf_tdist_Q (tval, df);
105 double sig = 2.0 * (tval > 0 ? q : p);
106 double tval_qinv = gsl_cdf_tdist_Qinv ((1.0 - tt->confidence) / 2.0, df);
107 double lower = mean_diff - tval_qinv * se_mean;
108 double upper = mean_diff + tval_qinv * se_mean;
109
110 double entries[] = { tval, df, sig, mean_diff, lower, upper };
111 for (size_t j = 0; j < sizeof entries / sizeof *entries; j++)
112 pivot_table_put2 (table, j, dep_var_idx,
113 pivot_value_new_number (entries[j]));
resolvePendingFixups()114 }
115
116 pivot_table_submit (table);
117 }
118
119 static void
120 one_sample_summary (const struct tt *tt, const struct one_samp *os)
121 {
122 struct pivot_table *table = pivot_table_create (N_("One-Sample Statistics"));
123 pivot_table_set_weight_var (table, tt->wv);
124
125 pivot_dimension_create (table, PIVOT_AXIS_COLUMN, N_("Statistics"),
126 N_("N"), PIVOT_RC_COUNT,
127 N_("Mean"), PIVOT_RC_OTHER,
128 N_("Std. Deviation"), PIVOT_RC_OTHER,
129 N_("S.E. Mean"), PIVOT_RC_OTHER);
absoluteSymbolDiff(const MCSymbol * Hi,const MCSymbol * Lo)130
131 struct pivot_dimension *variables = pivot_dimension_create (
132 table, PIVOT_AXIS_ROW, N_("Variables"));
133
134 for (size_t i = 0; i < os->n_stats; i++)
135 {
136 const struct per_var_stats *per_var_stats = &os->stats[i];
137 const struct moments *m = per_var_stats->mom;
138
139 int var_idx = pivot_category_create_leaf (
emitAbsoluteSymbolDiff(const MCSymbol * Hi,const MCSymbol * Lo,unsigned Size)140 variables->root, pivot_value_new_variable (per_var_stats->var));
141
142 double cc, mean, sigma;
143 moments_calculate (m, &cc, &mean, &sigma, NULL, NULL);
144
145 double entries[] = { cc, mean, sqrt (sigma), sqrt (sigma / cc) };
146 for (size_t j = 0; j < sizeof entries / sizeof *entries; j++)
147 pivot_table_put2 (table, j, var_idx,
148 pivot_value_new_number (entries[j]));
emitAbsoluteSymbolDiffAsULEB128(const MCSymbol * Hi,const MCSymbol * Lo)149 }
150
151 pivot_table_submit (table);
152 }
153
154 void
155 one_sample_run (const struct tt *tt, double testval, struct casereader *reader)
156 {
reset()157 struct one_samp os;
158 os.testval = testval;
159 os.stats = xcalloc (tt->n_vars, sizeof *os.stats);
160 os.n_stats = tt->n_vars;
161 for (size_t i = 0; i < tt->n_vars; ++i)
162 {
163 struct per_var_stats *per_var_stats = &os.stats[i];
164 per_var_stats->var = tt->vars[i];
165 per_var_stats->mom = moments_create (MOMENT_VARIANCE);
166 }
167
emitFrames(MCAsmBackend * MAB)168 struct casereader *r = casereader_clone (reader);
169 struct ccase *c;
170 for (; (c = casereader_read (r)); case_unref (c))
171 {
172 double w = dict_get_case_weight (tt->dict, c, NULL);
173 for (size_t i = 0; i < os.n_stats; i++)
174 {
175 const struct per_var_stats *per_var_stats = &os.stats[i];
176 const struct variable *var = per_var_stats->var;
177 const union value *val = case_data (c, var);
178 if (var_is_value_missing (var, val, tt->exclude))
getCurrentFragment() const179 continue;
180
181 moments_pass_one (per_var_stats->mom, val->f, w);
182 }
183 }
184 casereader_destroy (r);
185
186 r = reader;
187 for (; (c = casereader_read (r)); case_unref (c))
canReuseDataFragment(const MCDataFragment & F,const MCAssembler & Assembler,const MCSubtargetInfo * STI)188 {
189 double w = dict_get_case_weight (tt->dict, c, NULL);
190 for (size_t i = 0; i < os.n_stats; i++)
191 {
192 struct per_var_stats *per_var_stats = &os.stats[i];
193 const struct variable *var = per_var_stats->var;
194 const union value *val = case_data (c, var);
195 if (var_is_value_missing (var, val, tt->exclude))
196 continue;
197
198 moments_pass_two (per_var_stats->mom, val->f, w);
199 per_var_stats->sum_diff += w * (val->f - os.testval);
200 }
201 }
202 casereader_destroy (r);
getOrCreateDataFragment(const MCSubtargetInfo * STI)203
204 one_sample_summary (tt, &os);
205 one_sample_test (tt, &os);
206
207 for (size_t i = 0; i < os.n_stats; i++)
208 {
209 const struct per_var_stats *per_var_stats = &os.stats[i];
210 moments_destroy (per_var_stats->mom);
211 }
visitUsedSymbol(const MCSymbol & Sym)212 free (os.stats);
213 }
214
215