line |
stmt |
bran |
cond |
sub |
pod |
time |
code |
1
|
|
|
|
|
|
|
|
2
|
|
|
|
|
|
|
# |
3
|
|
|
|
|
|
|
# GENERATED WITH PDL::PP! Don't modify! |
4
|
|
|
|
|
|
|
# |
5
|
|
|
|
|
|
|
package PDL::Stats::Basic; |
6
|
|
|
|
|
|
|
|
7
|
|
|
|
|
|
|
@EXPORT_OK = qw( binomial_test rtable which_id PDL::PP stdv PDL::PP stdv_unbiased PDL::PP var PDL::PP var_unbiased PDL::PP se PDL::PP ss PDL::PP skew PDL::PP skew_unbiased PDL::PP kurt PDL::PP kurt_unbiased PDL::PP cov PDL::PP cov_table PDL::PP corr PDL::PP corr_table PDL::PP t_corr PDL::PP n_pair PDL::PP corr_dev PDL::PP t_test PDL::PP t_test_nev PDL::PP t_test_paired ); |
8
|
|
|
|
|
|
|
%EXPORT_TAGS = (Func=>[@EXPORT_OK]); |
9
|
|
|
|
|
|
|
|
10
|
5
|
|
|
5
|
|
147103
|
use PDL::Core; |
|
5
|
|
|
|
|
72629
|
|
|
5
|
|
|
|
|
31
|
|
11
|
5
|
|
|
5
|
|
1104
|
use PDL::Exporter; |
|
5
|
|
|
|
|
5
|
|
|
5
|
|
|
|
|
21
|
|
12
|
5
|
|
|
5
|
|
87
|
use DynaLoader; |
|
5
|
|
|
|
|
5
|
|
|
5
|
|
|
|
|
195
|
|
13
|
|
|
|
|
|
|
|
14
|
|
|
|
|
|
|
|
15
|
|
|
|
|
|
|
|
16
|
|
|
|
|
|
|
|
17
|
|
|
|
|
|
|
@ISA = ( 'PDL::Exporter','DynaLoader' ); |
18
|
|
|
|
|
|
|
push @PDL::Core::PP, __PACKAGE__; |
19
|
|
|
|
|
|
|
bootstrap PDL::Stats::Basic ; |
20
|
|
|
|
|
|
|
|
21
|
|
|
|
|
|
|
|
22
|
|
|
|
|
|
|
|
23
|
|
|
|
|
|
|
|
24
|
|
|
|
|
|
|
|
25
|
5
|
|
|
5
|
|
680
|
use PDL::LiteF; |
|
5
|
|
|
|
|
903
|
|
|
5
|
|
|
|
|
27
|
|
26
|
5
|
|
|
5
|
|
143093
|
use PDL::NiceSlice; |
|
5
|
|
|
|
|
7
|
|
|
5
|
|
|
|
|
25
|
|
27
|
5
|
|
|
5
|
|
56244
|
use Carp; |
|
5
|
|
|
|
|
6
|
|
|
5
|
|
|
|
|
7419
|
|
28
|
|
|
|
|
|
|
|
29
|
|
|
|
|
|
|
$PDL::onlinedoc->scan(__FILE__) if $PDL::onlinedoc; |
30
|
|
|
|
|
|
|
|
31
|
|
|
|
|
|
|
eval { require PDL::GSL::CDF; }; |
32
|
|
|
|
|
|
|
my $CDF = 1 if !$@; |
33
|
|
|
|
|
|
|
|
34
|
|
|
|
|
|
|
=head1 NAME |
35
|
|
|
|
|
|
|
|
36
|
|
|
|
|
|
|
PDL::Stats::Basic -- basic statistics and related utilities such as standard deviation, Pearson correlation, and t-tests. |
37
|
|
|
|
|
|
|
|
38
|
|
|
|
|
|
|
=head1 DESCRIPTION |
39
|
|
|
|
|
|
|
|
40
|
|
|
|
|
|
|
The terms FUNCTIONS and METHODS are arbitrarily used to refer to methods that are threadable and methods that are NOT threadable, respectively. |
41
|
|
|
|
|
|
|
|
42
|
|
|
|
|
|
|
Does not have mean or median function here. see SEE ALSO. |
43
|
|
|
|
|
|
|
|
44
|
|
|
|
|
|
|
=head1 SYNOPSIS |
45
|
|
|
|
|
|
|
|
46
|
|
|
|
|
|
|
use PDL::LiteF; |
47
|
|
|
|
|
|
|
use PDL::NiceSlice; |
48
|
|
|
|
|
|
|
use PDL::Stats::Basic; |
49
|
|
|
|
|
|
|
|
50
|
|
|
|
|
|
|
my $stdv = $data->stdv; |
51
|
|
|
|
|
|
|
|
52
|
|
|
|
|
|
|
or |
53
|
|
|
|
|
|
|
|
54
|
|
|
|
|
|
|
my $stdv = stdv( $data ); |
55
|
|
|
|
|
|
|
|
56
|
|
|
|
|
|
|
=cut |
57
|
|
|
|
|
|
|
|
58
|
|
|
|
|
|
|
|
59
|
|
|
|
|
|
|
|
60
|
|
|
|
|
|
|
|
61
|
|
|
|
|
|
|
|
62
|
|
|
|
|
|
|
|
63
|
|
|
|
|
|
|
|
64
|
|
|
|
|
|
|
=head1 FUNCTIONS |
65
|
|
|
|
|
|
|
|
66
|
|
|
|
|
|
|
|
67
|
|
|
|
|
|
|
|
68
|
|
|
|
|
|
|
=cut |
69
|
|
|
|
|
|
|
|
70
|
|
|
|
|
|
|
|
71
|
|
|
|
|
|
|
|
72
|
|
|
|
|
|
|
|
73
|
|
|
|
|
|
|
|
74
|
|
|
|
|
|
|
|
75
|
|
|
|
|
|
|
=head2 stdv |
76
|
|
|
|
|
|
|
|
77
|
|
|
|
|
|
|
=for sig |
78
|
|
|
|
|
|
|
|
79
|
|
|
|
|
|
|
Signature: (a(n); float+ [o]b()) |
80
|
|
|
|
|
|
|
|
81
|
|
|
|
|
|
|
|
82
|
|
|
|
|
|
|
|
83
|
|
|
|
|
|
|
=for ref |
84
|
|
|
|
|
|
|
|
85
|
|
|
|
|
|
|
Sample standard deviation. |
86
|
|
|
|
|
|
|
|
87
|
|
|
|
|
|
|
=cut |
88
|
|
|
|
|
|
|
|
89
|
|
|
|
|
|
|
|
90
|
|
|
|
|
|
|
=for bad |
91
|
|
|
|
|
|
|
|
92
|
|
|
|
|
|
|
stdv processes bad values. |
93
|
|
|
|
|
|
|
It will set the bad-value flag of all output piddles if the flag is set for any of the input piddles. |
94
|
|
|
|
|
|
|
|
95
|
|
|
|
|
|
|
|
96
|
|
|
|
|
|
|
=cut |
97
|
|
|
|
|
|
|
|
98
|
|
|
|
|
|
|
|
99
|
|
|
|
|
|
|
|
100
|
|
|
|
|
|
|
|
101
|
|
|
|
|
|
|
|
102
|
|
|
|
|
|
|
|
103
|
|
|
|
|
|
|
*stdv = \&PDL::stdv; |
104
|
|
|
|
|
|
|
|
105
|
|
|
|
|
|
|
|
106
|
|
|
|
|
|
|
|
107
|
|
|
|
|
|
|
|
108
|
|
|
|
|
|
|
|
109
|
|
|
|
|
|
|
=head2 stdv_unbiased |
110
|
|
|
|
|
|
|
|
111
|
|
|
|
|
|
|
=for sig |
112
|
|
|
|
|
|
|
|
113
|
|
|
|
|
|
|
Signature: (a(n); float+ [o]b()) |
114
|
|
|
|
|
|
|
|
115
|
|
|
|
|
|
|
|
116
|
|
|
|
|
|
|
|
117
|
|
|
|
|
|
|
=for ref |
118
|
|
|
|
|
|
|
|
119
|
|
|
|
|
|
|
Unbiased estimate of population standard deviation. |
120
|
|
|
|
|
|
|
|
121
|
|
|
|
|
|
|
=cut |
122
|
|
|
|
|
|
|
|
123
|
|
|
|
|
|
|
|
124
|
|
|
|
|
|
|
=for bad |
125
|
|
|
|
|
|
|
|
126
|
|
|
|
|
|
|
stdv_unbiased processes bad values. |
127
|
|
|
|
|
|
|
It will set the bad-value flag of all output piddles if the flag is set for any of the input piddles. |
128
|
|
|
|
|
|
|
|
129
|
|
|
|
|
|
|
|
130
|
|
|
|
|
|
|
=cut |
131
|
|
|
|
|
|
|
|
132
|
|
|
|
|
|
|
|
133
|
|
|
|
|
|
|
|
134
|
|
|
|
|
|
|
|
135
|
|
|
|
|
|
|
|
136
|
|
|
|
|
|
|
|
137
|
|
|
|
|
|
|
*stdv_unbiased = \&PDL::stdv_unbiased; |
138
|
|
|
|
|
|
|
|
139
|
|
|
|
|
|
|
|
140
|
|
|
|
|
|
|
|
141
|
|
|
|
|
|
|
|
142
|
|
|
|
|
|
|
|
143
|
|
|
|
|
|
|
=head2 var |
144
|
|
|
|
|
|
|
|
145
|
|
|
|
|
|
|
=for sig |
146
|
|
|
|
|
|
|
|
147
|
|
|
|
|
|
|
Signature: (a(n); float+ [o]b()) |
148
|
|
|
|
|
|
|
|
149
|
|
|
|
|
|
|
|
150
|
|
|
|
|
|
|
|
151
|
|
|
|
|
|
|
=for ref |
152
|
|
|
|
|
|
|
|
153
|
|
|
|
|
|
|
Sample variance. |
154
|
|
|
|
|
|
|
|
155
|
|
|
|
|
|
|
=cut |
156
|
|
|
|
|
|
|
|
157
|
|
|
|
|
|
|
|
158
|
|
|
|
|
|
|
=for bad |
159
|
|
|
|
|
|
|
|
160
|
|
|
|
|
|
|
var processes bad values. |
161
|
|
|
|
|
|
|
It will set the bad-value flag of all output piddles if the flag is set for any of the input piddles. |
162
|
|
|
|
|
|
|
|
163
|
|
|
|
|
|
|
|
164
|
|
|
|
|
|
|
=cut |
165
|
|
|
|
|
|
|
|
166
|
|
|
|
|
|
|
|
167
|
|
|
|
|
|
|
|
168
|
|
|
|
|
|
|
|
169
|
|
|
|
|
|
|
|
170
|
|
|
|
|
|
|
|
171
|
|
|
|
|
|
|
*var = \&PDL::var; |
172
|
|
|
|
|
|
|
|
173
|
|
|
|
|
|
|
|
174
|
|
|
|
|
|
|
|
175
|
|
|
|
|
|
|
|
176
|
|
|
|
|
|
|
|
177
|
|
|
|
|
|
|
=head2 var_unbiased |
178
|
|
|
|
|
|
|
|
179
|
|
|
|
|
|
|
=for sig |
180
|
|
|
|
|
|
|
|
181
|
|
|
|
|
|
|
Signature: (a(n); float+ [o]b()) |
182
|
|
|
|
|
|
|
|
183
|
|
|
|
|
|
|
|
184
|
|
|
|
|
|
|
|
185
|
|
|
|
|
|
|
=for ref |
186
|
|
|
|
|
|
|
|
187
|
|
|
|
|
|
|
Unbiased estimate of population variance. |
188
|
|
|
|
|
|
|
|
189
|
|
|
|
|
|
|
=cut |
190
|
|
|
|
|
|
|
|
191
|
|
|
|
|
|
|
|
192
|
|
|
|
|
|
|
=for bad |
193
|
|
|
|
|
|
|
|
194
|
|
|
|
|
|
|
var_unbiased processes bad values. |
195
|
|
|
|
|
|
|
It will set the bad-value flag of all output piddles if the flag is set for any of the input piddles. |
196
|
|
|
|
|
|
|
|
197
|
|
|
|
|
|
|
|
198
|
|
|
|
|
|
|
=cut |
199
|
|
|
|
|
|
|
|
200
|
|
|
|
|
|
|
|
201
|
|
|
|
|
|
|
|
202
|
|
|
|
|
|
|
|
203
|
|
|
|
|
|
|
|
204
|
|
|
|
|
|
|
|
205
|
|
|
|
|
|
|
*var_unbiased = \&PDL::var_unbiased; |
206
|
|
|
|
|
|
|
|
207
|
|
|
|
|
|
|
|
208
|
|
|
|
|
|
|
|
209
|
|
|
|
|
|
|
|
210
|
|
|
|
|
|
|
|
211
|
|
|
|
|
|
|
=head2 se |
212
|
|
|
|
|
|
|
|
213
|
|
|
|
|
|
|
=for sig |
214
|
|
|
|
|
|
|
|
215
|
|
|
|
|
|
|
Signature: (a(n); float+ [o]b()) |
216
|
|
|
|
|
|
|
|
217
|
|
|
|
|
|
|
|
218
|
|
|
|
|
|
|
|
219
|
|
|
|
|
|
|
=for ref |
220
|
|
|
|
|
|
|
|
221
|
|
|
|
|
|
|
Standard error of the mean. Useful for calculating confidence intervals. |
222
|
|
|
|
|
|
|
|
223
|
|
|
|
|
|
|
=for usage |
224
|
|
|
|
|
|
|
|
225
|
|
|
|
|
|
|
# 95% confidence interval for samples with large N |
226
|
|
|
|
|
|
|
|
227
|
|
|
|
|
|
|
$ci_95_upper = $data->average + 1.96 * $data->se; |
228
|
|
|
|
|
|
|
$ci_95_lower = $data->average - 1.96 * $data->se; |
229
|
|
|
|
|
|
|
|
230
|
|
|
|
|
|
|
|
231
|
|
|
|
|
|
|
|
232
|
|
|
|
|
|
|
|
233
|
|
|
|
|
|
|
=for bad |
234
|
|
|
|
|
|
|
|
235
|
|
|
|
|
|
|
se processes bad values. |
236
|
|
|
|
|
|
|
It will set the bad-value flag of all output piddles if the flag is set for any of the input piddles. |
237
|
|
|
|
|
|
|
|
238
|
|
|
|
|
|
|
|
239
|
|
|
|
|
|
|
=cut |
240
|
|
|
|
|
|
|
|
241
|
|
|
|
|
|
|
|
242
|
|
|
|
|
|
|
|
243
|
|
|
|
|
|
|
|
244
|
|
|
|
|
|
|
|
245
|
|
|
|
|
|
|
|
246
|
|
|
|
|
|
|
*se = \&PDL::se; |
247
|
|
|
|
|
|
|
|
248
|
|
|
|
|
|
|
|
249
|
|
|
|
|
|
|
|
250
|
|
|
|
|
|
|
|
251
|
|
|
|
|
|
|
|
252
|
|
|
|
|
|
|
=head2 ss |
253
|
|
|
|
|
|
|
|
254
|
|
|
|
|
|
|
=for sig |
255
|
|
|
|
|
|
|
|
256
|
|
|
|
|
|
|
Signature: (a(n); float+ [o]b()) |
257
|
|
|
|
|
|
|
|
258
|
|
|
|
|
|
|
|
259
|
|
|
|
|
|
|
|
260
|
|
|
|
|
|
|
=for ref |
261
|
|
|
|
|
|
|
|
262
|
|
|
|
|
|
|
Sum of squared deviations from the mean. |
263
|
|
|
|
|
|
|
|
264
|
|
|
|
|
|
|
=cut |
265
|
|
|
|
|
|
|
|
266
|
|
|
|
|
|
|
|
267
|
|
|
|
|
|
|
=for bad |
268
|
|
|
|
|
|
|
|
269
|
|
|
|
|
|
|
ss processes bad values. |
270
|
|
|
|
|
|
|
It will set the bad-value flag of all output piddles if the flag is set for any of the input piddles. |
271
|
|
|
|
|
|
|
|
272
|
|
|
|
|
|
|
|
273
|
|
|
|
|
|
|
=cut |
274
|
|
|
|
|
|
|
|
275
|
|
|
|
|
|
|
|
276
|
|
|
|
|
|
|
|
277
|
|
|
|
|
|
|
|
278
|
|
|
|
|
|
|
|
279
|
|
|
|
|
|
|
|
280
|
|
|
|
|
|
|
*ss = \&PDL::ss; |
281
|
|
|
|
|
|
|
|
282
|
|
|
|
|
|
|
|
283
|
|
|
|
|
|
|
|
284
|
|
|
|
|
|
|
|
285
|
|
|
|
|
|
|
|
286
|
|
|
|
|
|
|
=head2 skew |
287
|
|
|
|
|
|
|
|
288
|
|
|
|
|
|
|
=for sig |
289
|
|
|
|
|
|
|
|
290
|
|
|
|
|
|
|
Signature: (a(n); float+ [o]b()) |
291
|
|
|
|
|
|
|
|
292
|
|
|
|
|
|
|
|
293
|
|
|
|
|
|
|
|
294
|
|
|
|
|
|
|
=for ref |
295
|
|
|
|
|
|
|
|
296
|
|
|
|
|
|
|
Sample skewness, measure of asymmetry in data. skewness == 0 for normal distribution. |
297
|
|
|
|
|
|
|
|
298
|
|
|
|
|
|
|
=cut |
299
|
|
|
|
|
|
|
|
300
|
|
|
|
|
|
|
|
301
|
|
|
|
|
|
|
=for bad |
302
|
|
|
|
|
|
|
|
303
|
|
|
|
|
|
|
skew processes bad values. |
304
|
|
|
|
|
|
|
It will set the bad-value flag of all output piddles if the flag is set for any of the input piddles. |
305
|
|
|
|
|
|
|
|
306
|
|
|
|
|
|
|
|
307
|
|
|
|
|
|
|
=cut |
308
|
|
|
|
|
|
|
|
309
|
|
|
|
|
|
|
|
310
|
|
|
|
|
|
|
|
311
|
|
|
|
|
|
|
|
312
|
|
|
|
|
|
|
|
313
|
|
|
|
|
|
|
|
314
|
|
|
|
|
|
|
*skew = \&PDL::skew; |
315
|
|
|
|
|
|
|
|
316
|
|
|
|
|
|
|
|
317
|
|
|
|
|
|
|
|
318
|
|
|
|
|
|
|
|
319
|
|
|
|
|
|
|
|
320
|
|
|
|
|
|
|
=head2 skew_unbiased |
321
|
|
|
|
|
|
|
|
322
|
|
|
|
|
|
|
=for sig |
323
|
|
|
|
|
|
|
|
324
|
|
|
|
|
|
|
Signature: (a(n); float+ [o]b()) |
325
|
|
|
|
|
|
|
|
326
|
|
|
|
|
|
|
|
327
|
|
|
|
|
|
|
|
328
|
|
|
|
|
|
|
=for ref |
329
|
|
|
|
|
|
|
|
330
|
|
|
|
|
|
|
Unbiased estimate of population skewness. This is the number in GNumeric Descriptive Statistics. |
331
|
|
|
|
|
|
|
|
332
|
|
|
|
|
|
|
=cut |
333
|
|
|
|
|
|
|
|
334
|
|
|
|
|
|
|
|
335
|
|
|
|
|
|
|
=for bad |
336
|
|
|
|
|
|
|
|
337
|
|
|
|
|
|
|
skew_unbiased processes bad values. |
338
|
|
|
|
|
|
|
It will set the bad-value flag of all output piddles if the flag is set for any of the input piddles. |
339
|
|
|
|
|
|
|
|
340
|
|
|
|
|
|
|
|
341
|
|
|
|
|
|
|
=cut |
342
|
|
|
|
|
|
|
|
343
|
|
|
|
|
|
|
|
344
|
|
|
|
|
|
|
|
345
|
|
|
|
|
|
|
|
346
|
|
|
|
|
|
|
|
347
|
|
|
|
|
|
|
|
348
|
|
|
|
|
|
|
*skew_unbiased = \&PDL::skew_unbiased; |
349
|
|
|
|
|
|
|
|
350
|
|
|
|
|
|
|
|
351
|
|
|
|
|
|
|
|
352
|
|
|
|
|
|
|
|
353
|
|
|
|
|
|
|
|
354
|
|
|
|
|
|
|
=head2 kurt |
355
|
|
|
|
|
|
|
|
356
|
|
|
|
|
|
|
=for sig |
357
|
|
|
|
|
|
|
|
358
|
|
|
|
|
|
|
Signature: (a(n); float+ [o]b()) |
359
|
|
|
|
|
|
|
|
360
|
|
|
|
|
|
|
|
361
|
|
|
|
|
|
|
|
362
|
|
|
|
|
|
|
=for ref |
363
|
|
|
|
|
|
|
|
364
|
|
|
|
|
|
|
Sample kurtosis, measure of "peakedness" of data. kurtosis == 0 for normal distribution. |
365
|
|
|
|
|
|
|
|
366
|
|
|
|
|
|
|
=cut |
367
|
|
|
|
|
|
|
|
368
|
|
|
|
|
|
|
|
369
|
|
|
|
|
|
|
=for bad |
370
|
|
|
|
|
|
|
|
371
|
|
|
|
|
|
|
kurt processes bad values. |
372
|
|
|
|
|
|
|
It will set the bad-value flag of all output piddles if the flag is set for any of the input piddles. |
373
|
|
|
|
|
|
|
|
374
|
|
|
|
|
|
|
|
375
|
|
|
|
|
|
|
=cut |
376
|
|
|
|
|
|
|
|
377
|
|
|
|
|
|
|
|
378
|
|
|
|
|
|
|
|
379
|
|
|
|
|
|
|
|
380
|
|
|
|
|
|
|
|
381
|
|
|
|
|
|
|
|
382
|
|
|
|
|
|
|
*kurt = \&PDL::kurt; |
383
|
|
|
|
|
|
|
|
384
|
|
|
|
|
|
|
|
385
|
|
|
|
|
|
|
|
386
|
|
|
|
|
|
|
|
387
|
|
|
|
|
|
|
|
388
|
|
|
|
|
|
|
=head2 kurt_unbiased |
389
|
|
|
|
|
|
|
|
390
|
|
|
|
|
|
|
=for sig |
391
|
|
|
|
|
|
|
|
392
|
|
|
|
|
|
|
Signature: (a(n); float+ [o]b()) |
393
|
|
|
|
|
|
|
|
394
|
|
|
|
|
|
|
|
395
|
|
|
|
|
|
|
|
396
|
|
|
|
|
|
|
=for ref |
397
|
|
|
|
|
|
|
|
398
|
|
|
|
|
|
|
Unbiased estimate of population kurtosis. This is the number in GNumeric Descriptive Statistics. |
399
|
|
|
|
|
|
|
|
400
|
|
|
|
|
|
|
=cut |
401
|
|
|
|
|
|
|
|
402
|
|
|
|
|
|
|
|
403
|
|
|
|
|
|
|
=for bad |
404
|
|
|
|
|
|
|
|
405
|
|
|
|
|
|
|
kurt_unbiased processes bad values. |
406
|
|
|
|
|
|
|
It will set the bad-value flag of all output piddles if the flag is set for any of the input piddles. |
407
|
|
|
|
|
|
|
|
408
|
|
|
|
|
|
|
|
409
|
|
|
|
|
|
|
=cut |
410
|
|
|
|
|
|
|
|
411
|
|
|
|
|
|
|
|
412
|
|
|
|
|
|
|
|
413
|
|
|
|
|
|
|
|
414
|
|
|
|
|
|
|
|
415
|
|
|
|
|
|
|
|
416
|
|
|
|
|
|
|
*kurt_unbiased = \&PDL::kurt_unbiased; |
417
|
|
|
|
|
|
|
|
418
|
|
|
|
|
|
|
|
419
|
|
|
|
|
|
|
|
420
|
|
|
|
|
|
|
|
421
|
|
|
|
|
|
|
|
422
|
|
|
|
|
|
|
=head2 cov |
423
|
|
|
|
|
|
|
|
424
|
|
|
|
|
|
|
=for sig |
425
|
|
|
|
|
|
|
|
426
|
|
|
|
|
|
|
Signature: (a(n); b(n); float+ [o]c()) |
427
|
|
|
|
|
|
|
|
428
|
|
|
|
|
|
|
|
429
|
|
|
|
|
|
|
|
430
|
|
|
|
|
|
|
=for ref |
431
|
|
|
|
|
|
|
|
432
|
|
|
|
|
|
|
Sample covariance. see B for ways to call |
433
|
|
|
|
|
|
|
|
434
|
|
|
|
|
|
|
=cut |
435
|
|
|
|
|
|
|
|
436
|
|
|
|
|
|
|
|
437
|
|
|
|
|
|
|
=for bad |
438
|
|
|
|
|
|
|
|
439
|
|
|
|
|
|
|
cov processes bad values. |
440
|
|
|
|
|
|
|
It will set the bad-value flag of all output piddles if the flag is set for any of the input piddles. |
441
|
|
|
|
|
|
|
|
442
|
|
|
|
|
|
|
|
443
|
|
|
|
|
|
|
=cut |
444
|
|
|
|
|
|
|
|
445
|
|
|
|
|
|
|
|
446
|
|
|
|
|
|
|
|
447
|
|
|
|
|
|
|
|
448
|
|
|
|
|
|
|
|
449
|
|
|
|
|
|
|
|
450
|
|
|
|
|
|
|
*cov = \&PDL::cov; |
451
|
|
|
|
|
|
|
|
452
|
|
|
|
|
|
|
|
453
|
|
|
|
|
|
|
|
454
|
|
|
|
|
|
|
|
455
|
|
|
|
|
|
|
|
456
|
|
|
|
|
|
|
=head2 cov_table |
457
|
|
|
|
|
|
|
|
458
|
|
|
|
|
|
|
=for sig |
459
|
|
|
|
|
|
|
|
460
|
|
|
|
|
|
|
Signature: (a(n,m); float+ [o]c(m,m)) |
461
|
|
|
|
|
|
|
|
462
|
|
|
|
|
|
|
|
463
|
|
|
|
|
|
|
|
464
|
|
|
|
|
|
|
=for ref |
465
|
|
|
|
|
|
|
|
466
|
|
|
|
|
|
|
Square covariance table. Gives the same result as threading using B but it calculates only half the square, hence much faster. And it is easier to use with higher dimension pdls. |
467
|
|
|
|
|
|
|
|
468
|
|
|
|
|
|
|
=for usage |
469
|
|
|
|
|
|
|
|
470
|
|
|
|
|
|
|
Usage: |
471
|
|
|
|
|
|
|
|
472
|
|
|
|
|
|
|
# 5 obs x 3 var, 2 such data tables |
473
|
|
|
|
|
|
|
|
474
|
|
|
|
|
|
|
perldl> $a = random 5, 3, 2 |
475
|
|
|
|
|
|
|
|
476
|
|
|
|
|
|
|
perldl> p $cov = $a->cov_table |
477
|
|
|
|
|
|
|
[ |
478
|
|
|
|
|
|
|
[ |
479
|
|
|
|
|
|
|
[ 8.9636438 -1.8624472 -1.2416588] |
480
|
|
|
|
|
|
|
[-1.8624472 14.341514 -1.4245366] |
481
|
|
|
|
|
|
|
[-1.2416588 -1.4245366 9.8690655] |
482
|
|
|
|
|
|
|
] |
483
|
|
|
|
|
|
|
[ |
484
|
|
|
|
|
|
|
[ 10.32644 -0.31311789 -0.95643674] |
485
|
|
|
|
|
|
|
[-0.31311789 15.051779 -7.2759577] |
486
|
|
|
|
|
|
|
[-0.95643674 -7.2759577 5.4465141] |
487
|
|
|
|
|
|
|
] |
488
|
|
|
|
|
|
|
] |
489
|
|
|
|
|
|
|
# diagonal elements of the cov table are the variances |
490
|
|
|
|
|
|
|
perldl> p $a->var |
491
|
|
|
|
|
|
|
[ |
492
|
|
|
|
|
|
|
[ 8.9636438 14.341514 9.8690655] |
493
|
|
|
|
|
|
|
[ 10.32644 15.051779 5.4465141] |
494
|
|
|
|
|
|
|
] |
495
|
|
|
|
|
|
|
|
496
|
|
|
|
|
|
|
for the same cov matrix table using B, |
497
|
|
|
|
|
|
|
|
498
|
|
|
|
|
|
|
perldl> p $a->dummy(2)->cov($a->dummy(1)) |
499
|
|
|
|
|
|
|
|
500
|
|
|
|
|
|
|
|
501
|
|
|
|
|
|
|
|
502
|
|
|
|
|
|
|
|
503
|
|
|
|
|
|
|
=for bad |
504
|
|
|
|
|
|
|
|
505
|
|
|
|
|
|
|
cov_table processes bad values. |
506
|
|
|
|
|
|
|
It will set the bad-value flag of all output piddles if the flag is set for any of the input piddles. |
507
|
|
|
|
|
|
|
|
508
|
|
|
|
|
|
|
|
509
|
|
|
|
|
|
|
=cut |
510
|
|
|
|
|
|
|
|
511
|
|
|
|
|
|
|
|
512
|
|
|
|
|
|
|
|
513
|
|
|
|
|
|
|
|
514
|
|
|
|
|
|
|
|
515
|
|
|
|
|
|
|
|
516
|
|
|
|
|
|
|
*cov_table = \&PDL::cov_table; |
517
|
|
|
|
|
|
|
|
518
|
|
|
|
|
|
|
|
519
|
|
|
|
|
|
|
|
520
|
|
|
|
|
|
|
|
521
|
|
|
|
|
|
|
|
522
|
|
|
|
|
|
|
=head2 corr |
523
|
|
|
|
|
|
|
|
524
|
|
|
|
|
|
|
=for sig |
525
|
|
|
|
|
|
|
|
526
|
|
|
|
|
|
|
Signature: (a(n); b(n); float+ [o]c()) |
527
|
|
|
|
|
|
|
|
528
|
|
|
|
|
|
|
|
529
|
|
|
|
|
|
|
|
530
|
|
|
|
|
|
|
=for ref |
531
|
|
|
|
|
|
|
|
532
|
|
|
|
|
|
|
Pearson correlation coefficient. r = cov(X,Y) / (stdv(X) * stdv(Y)). |
533
|
|
|
|
|
|
|
|
534
|
|
|
|
|
|
|
=for usage |
535
|
|
|
|
|
|
|
|
536
|
|
|
|
|
|
|
Usage: |
537
|
|
|
|
|
|
|
|
538
|
|
|
|
|
|
|
perldl> $a = random 5, 3 |
539
|
|
|
|
|
|
|
perldl> $b = sequence 5,3 |
540
|
|
|
|
|
|
|
perldl> p $a->corr($b) |
541
|
|
|
|
|
|
|
|
542
|
|
|
|
|
|
|
[0.20934208 0.30949881 0.26713007] |
543
|
|
|
|
|
|
|
|
544
|
|
|
|
|
|
|
for square corr table |
545
|
|
|
|
|
|
|
|
546
|
|
|
|
|
|
|
perldl> p $a->corr($a->dummy(1)) |
547
|
|
|
|
|
|
|
|
548
|
|
|
|
|
|
|
[ |
549
|
|
|
|
|
|
|
[ 1 -0.41995259 -0.029301192] |
550
|
|
|
|
|
|
|
[ -0.41995259 1 -0.61927619] |
551
|
|
|
|
|
|
|
[-0.029301192 -0.61927619 1] |
552
|
|
|
|
|
|
|
] |
553
|
|
|
|
|
|
|
|
554
|
|
|
|
|
|
|
but it is easier and faster to use B. |
555
|
|
|
|
|
|
|
|
556
|
|
|
|
|
|
|
=cut |
557
|
|
|
|
|
|
|
|
558
|
|
|
|
|
|
|
|
559
|
|
|
|
|
|
|
=for bad |
560
|
|
|
|
|
|
|
|
561
|
|
|
|
|
|
|
corr processes bad values. |
562
|
|
|
|
|
|
|
It will set the bad-value flag of all output piddles if the flag is set for any of the input piddles. |
563
|
|
|
|
|
|
|
|
564
|
|
|
|
|
|
|
|
565
|
|
|
|
|
|
|
=cut |
566
|
|
|
|
|
|
|
|
567
|
|
|
|
|
|
|
|
568
|
|
|
|
|
|
|
|
569
|
|
|
|
|
|
|
|
570
|
|
|
|
|
|
|
|
571
|
|
|
|
|
|
|
|
572
|
|
|
|
|
|
|
*corr = \&PDL::corr; |
573
|
|
|
|
|
|
|
|
574
|
|
|
|
|
|
|
|
575
|
|
|
|
|
|
|
|
576
|
|
|
|
|
|
|
|
577
|
|
|
|
|
|
|
|
578
|
|
|
|
|
|
|
=head2 corr_table |
579
|
|
|
|
|
|
|
|
580
|
|
|
|
|
|
|
=for sig |
581
|
|
|
|
|
|
|
|
582
|
|
|
|
|
|
|
Signature: (a(n,m); float+ [o]c(m,m)) |
583
|
|
|
|
|
|
|
|
584
|
|
|
|
|
|
|
|
585
|
|
|
|
|
|
|
|
586
|
|
|
|
|
|
|
=for ref |
587
|
|
|
|
|
|
|
|
588
|
|
|
|
|
|
|
Square Pearson correlation table. Gives the same result as threading using B but it calculates only half the square, hence much faster. And it is easier to use with higher dimension pdls. |
589
|
|
|
|
|
|
|
|
590
|
|
|
|
|
|
|
=for usage |
591
|
|
|
|
|
|
|
|
592
|
|
|
|
|
|
|
Usage: |
593
|
|
|
|
|
|
|
|
594
|
|
|
|
|
|
|
# 5 obs x 3 var, 2 such data tables |
595
|
|
|
|
|
|
|
|
596
|
|
|
|
|
|
|
perldl> $a = random 5, 3, 2 |
597
|
|
|
|
|
|
|
|
598
|
|
|
|
|
|
|
perldl> p $a->corr_table |
599
|
|
|
|
|
|
|
[ |
600
|
|
|
|
|
|
|
[ |
601
|
|
|
|
|
|
|
[ 1 -0.69835951 -0.18549048] |
602
|
|
|
|
|
|
|
[-0.69835951 1 0.72481605] |
603
|
|
|
|
|
|
|
[-0.18549048 0.72481605 1] |
604
|
|
|
|
|
|
|
] |
605
|
|
|
|
|
|
|
[ |
606
|
|
|
|
|
|
|
[ 1 0.82722569 -0.71779883] |
607
|
|
|
|
|
|
|
[ 0.82722569 1 -0.63938828] |
608
|
|
|
|
|
|
|
[-0.71779883 -0.63938828 1] |
609
|
|
|
|
|
|
|
] |
610
|
|
|
|
|
|
|
] |
611
|
|
|
|
|
|
|
|
612
|
|
|
|
|
|
|
for the same result using B, |
613
|
|
|
|
|
|
|
|
614
|
|
|
|
|
|
|
perldl> p $a->dummy(2)->corr($a->dummy(1)) |
615
|
|
|
|
|
|
|
|
616
|
|
|
|
|
|
|
This is also how to use B and B with such a table. |
617
|
|
|
|
|
|
|
|
618
|
|
|
|
|
|
|
|
619
|
|
|
|
|
|
|
|
620
|
|
|
|
|
|
|
|
621
|
|
|
|
|
|
|
=for bad |
622
|
|
|
|
|
|
|
|
623
|
|
|
|
|
|
|
corr_table processes bad values. |
624
|
|
|
|
|
|
|
It will set the bad-value flag of all output piddles if the flag is set for any of the input piddles. |
625
|
|
|
|
|
|
|
|
626
|
|
|
|
|
|
|
|
627
|
|
|
|
|
|
|
=cut |
628
|
|
|
|
|
|
|
|
629
|
|
|
|
|
|
|
|
630
|
|
|
|
|
|
|
|
631
|
|
|
|
|
|
|
|
632
|
|
|
|
|
|
|
|
633
|
|
|
|
|
|
|
|
634
|
|
|
|
|
|
|
*corr_table = \&PDL::corr_table; |
635
|
|
|
|
|
|
|
|
636
|
|
|
|
|
|
|
|
637
|
|
|
|
|
|
|
|
638
|
|
|
|
|
|
|
|
639
|
|
|
|
|
|
|
|
640
|
|
|
|
|
|
|
=head2 t_corr |
641
|
|
|
|
|
|
|
|
642
|
|
|
|
|
|
|
=for sig |
643
|
|
|
|
|
|
|
|
644
|
|
|
|
|
|
|
Signature: (r(); n(); [o]t()) |
645
|
|
|
|
|
|
|
|
646
|
|
|
|
|
|
|
|
647
|
|
|
|
|
|
|
|
648
|
|
|
|
|
|
|
=for usage |
649
|
|
|
|
|
|
|
|
650
|
|
|
|
|
|
|
$corr = $data->corr( $data->dummy(1) ); |
651
|
|
|
|
|
|
|
$n = $data->n_pair( $data->dummy(1) ); |
652
|
|
|
|
|
|
|
$t_corr = $corr->t_corr( $n ); |
653
|
|
|
|
|
|
|
|
654
|
|
|
|
|
|
|
use PDL::GSL::CDF; |
655
|
|
|
|
|
|
|
|
656
|
|
|
|
|
|
|
$p_2tail = 2 * (1 - gsl_cdf_tdist_P( $t_corr->abs, $n-2 )); |
657
|
|
|
|
|
|
|
|
658
|
|
|
|
|
|
|
=for ref |
659
|
|
|
|
|
|
|
|
660
|
|
|
|
|
|
|
t significance test for Pearson correlations. |
661
|
|
|
|
|
|
|
|
662
|
|
|
|
|
|
|
=cut |
663
|
|
|
|
|
|
|
|
664
|
|
|
|
|
|
|
|
665
|
|
|
|
|
|
|
=for bad |
666
|
|
|
|
|
|
|
|
667
|
|
|
|
|
|
|
t_corr processes bad values. |
668
|
|
|
|
|
|
|
It will set the bad-value flag of all output piddles if the flag is set for any of the input piddles. |
669
|
|
|
|
|
|
|
|
670
|
|
|
|
|
|
|
|
671
|
|
|
|
|
|
|
=cut |
672
|
|
|
|
|
|
|
|
673
|
|
|
|
|
|
|
|
674
|
|
|
|
|
|
|
|
675
|
|
|
|
|
|
|
|
676
|
|
|
|
|
|
|
|
677
|
|
|
|
|
|
|
|
678
|
|
|
|
|
|
|
*t_corr = \&PDL::t_corr; |
679
|
|
|
|
|
|
|
|
680
|
|
|
|
|
|
|
|
681
|
|
|
|
|
|
|
|
682
|
|
|
|
|
|
|
|
683
|
|
|
|
|
|
|
|
684
|
|
|
|
|
|
|
=head2 n_pair |
685
|
|
|
|
|
|
|
|
686
|
|
|
|
|
|
|
=for sig |
687
|
|
|
|
|
|
|
|
688
|
|
|
|
|
|
|
Signature: (a(n); b(n); indx [o]c()) |
689
|
|
|
|
|
|
|
|
690
|
|
|
|
|
|
|
|
691
|
|
|
|
|
|
|
|
692
|
|
|
|
|
|
|
=for ref |
693
|
|
|
|
|
|
|
|
694
|
|
|
|
|
|
|
Returns the number of good pairs between 2 lists. Useful with B (esp. when bad values are involved) |
695
|
|
|
|
|
|
|
|
696
|
|
|
|
|
|
|
=cut |
697
|
|
|
|
|
|
|
|
698
|
|
|
|
|
|
|
|
699
|
|
|
|
|
|
|
=for bad |
700
|
|
|
|
|
|
|
|
701
|
|
|
|
|
|
|
n_pair processes bad values. |
702
|
|
|
|
|
|
|
It will set the bad-value flag of all output piddles if the flag is set for any of the input piddles. |
703
|
|
|
|
|
|
|
|
704
|
|
|
|
|
|
|
|
705
|
|
|
|
|
|
|
=cut |
706
|
|
|
|
|
|
|
|
707
|
|
|
|
|
|
|
|
708
|
|
|
|
|
|
|
|
709
|
|
|
|
|
|
|
|
710
|
|
|
|
|
|
|
|
711
|
|
|
|
|
|
|
|
712
|
|
|
|
|
|
|
*n_pair = \&PDL::n_pair; |
713
|
|
|
|
|
|
|
|
714
|
|
|
|
|
|
|
|
715
|
|
|
|
|
|
|
|
716
|
|
|
|
|
|
|
|
717
|
|
|
|
|
|
|
|
718
|
|
|
|
|
|
|
=head2 corr_dev |
719
|
|
|
|
|
|
|
|
720
|
|
|
|
|
|
|
=for sig |
721
|
|
|
|
|
|
|
|
722
|
|
|
|
|
|
|
Signature: (a(n); b(n); float+ [o]c()) |
723
|
|
|
|
|
|
|
|
724
|
|
|
|
|
|
|
|
725
|
|
|
|
|
|
|
|
726
|
|
|
|
|
|
|
=for usage |
727
|
|
|
|
|
|
|
|
728
|
|
|
|
|
|
|
$corr = $a->dev_m->corr_dev($b->dev_m); |
729
|
|
|
|
|
|
|
|
730
|
|
|
|
|
|
|
=for ref |
731
|
|
|
|
|
|
|
|
732
|
|
|
|
|
|
|
Calculates correlations from B vals. Seems faster than doing B from original vals when data pdl is big |
733
|
|
|
|
|
|
|
|
734
|
|
|
|
|
|
|
=cut |
735
|
|
|
|
|
|
|
|
736
|
|
|
|
|
|
|
|
737
|
|
|
|
|
|
|
=for bad |
738
|
|
|
|
|
|
|
|
739
|
|
|
|
|
|
|
corr_dev processes bad values. |
740
|
|
|
|
|
|
|
It will set the bad-value flag of all output piddles if the flag is set for any of the input piddles. |
741
|
|
|
|
|
|
|
|
742
|
|
|
|
|
|
|
|
743
|
|
|
|
|
|
|
=cut |
744
|
|
|
|
|
|
|
|
745
|
|
|
|
|
|
|
|
746
|
|
|
|
|
|
|
|
747
|
|
|
|
|
|
|
|
748
|
|
|
|
|
|
|
|
749
|
|
|
|
|
|
|
|
750
|
|
|
|
|
|
|
*corr_dev = \&PDL::corr_dev; |
751
|
|
|
|
|
|
|
|
752
|
|
|
|
|
|
|
|
753
|
|
|
|
|
|
|
|
754
|
|
|
|
|
|
|
|
755
|
|
|
|
|
|
|
|
756
|
|
|
|
|
|
|
=head2 t_test |
757
|
|
|
|
|
|
|
|
758
|
|
|
|
|
|
|
=for sig |
759
|
|
|
|
|
|
|
|
760
|
|
|
|
|
|
|
Signature: (a(n); b(m); float+ [o]t(); [o]d()) |
761
|
|
|
|
|
|
|
|
762
|
|
|
|
|
|
|
|
763
|
|
|
|
|
|
|
|
764
|
|
|
|
|
|
|
=for usage |
765
|
|
|
|
|
|
|
|
766
|
|
|
|
|
|
|
my ($t, $df) = t_test( $pdl1, $pdl2 ); |
767
|
|
|
|
|
|
|
|
768
|
|
|
|
|
|
|
use PDL::GSL::CDF; |
769
|
|
|
|
|
|
|
|
770
|
|
|
|
|
|
|
my $p_2tail = 2 * (1 - gsl_cdf_tdist_P( $t->abs, $df )); |
771
|
|
|
|
|
|
|
|
772
|
|
|
|
|
|
|
=for ref |
773
|
|
|
|
|
|
|
|
774
|
|
|
|
|
|
|
Independent sample t-test, assuming equal var. |
775
|
|
|
|
|
|
|
|
776
|
|
|
|
|
|
|
=cut |
777
|
|
|
|
|
|
|
|
778
|
|
|
|
|
|
|
|
779
|
|
|
|
|
|
|
=for bad |
780
|
|
|
|
|
|
|
|
781
|
|
|
|
|
|
|
t_test processes bad values. |
782
|
|
|
|
|
|
|
It will set the bad-value flag of all output piddles if the flag is set for any of the input piddles. |
783
|
|
|
|
|
|
|
|
784
|
|
|
|
|
|
|
|
785
|
|
|
|
|
|
|
=cut |
786
|
|
|
|
|
|
|
|
787
|
|
|
|
|
|
|
|
788
|
|
|
|
|
|
|
|
789
|
|
|
|
|
|
|
|
790
|
|
|
|
|
|
|
|
791
|
|
|
|
|
|
|
|
792
|
|
|
|
|
|
|
*t_test = \&PDL::t_test; |
793
|
|
|
|
|
|
|
|
794
|
|
|
|
|
|
|
|
795
|
|
|
|
|
|
|
|
796
|
|
|
|
|
|
|
|
797
|
|
|
|
|
|
|
|
798
|
|
|
|
|
|
|
=head2 t_test_nev |
799
|
|
|
|
|
|
|
|
800
|
|
|
|
|
|
|
=for sig |
801
|
|
|
|
|
|
|
|
802
|
|
|
|
|
|
|
Signature: (a(n); b(m); float+ [o]t(); [o]d()) |
803
|
|
|
|
|
|
|
|
804
|
|
|
|
|
|
|
|
805
|
|
|
|
|
|
|
|
806
|
|
|
|
|
|
|
=for ref |
807
|
|
|
|
|
|
|
|
808
|
|
|
|
|
|
|
Independent sample t-test, NOT assuming equal var. ie Welch two sample t test. Df follows Welch-Satterthwaite equation instead of Satterthwaite (1946, as cited by Hays, 1994, 5th ed.). It matches GNumeric, which matches R. |
809
|
|
|
|
|
|
|
|
810
|
|
|
|
|
|
|
=for usage |
811
|
|
|
|
|
|
|
|
812
|
|
|
|
|
|
|
my ($t, $df) = $pdl1->t_test( $pdl2 ); |
813
|
|
|
|
|
|
|
|
814
|
|
|
|
|
|
|
=cut |
815
|
|
|
|
|
|
|
|
816
|
|
|
|
|
|
|
|
817
|
|
|
|
|
|
|
=for bad |
818
|
|
|
|
|
|
|
|
819
|
|
|
|
|
|
|
t_test_nev processes bad values. |
820
|
|
|
|
|
|
|
It will set the bad-value flag of all output piddles if the flag is set for any of the input piddles. |
821
|
|
|
|
|
|
|
|
822
|
|
|
|
|
|
|
|
823
|
|
|
|
|
|
|
=cut |
824
|
|
|
|
|
|
|
|
825
|
|
|
|
|
|
|
|
826
|
|
|
|
|
|
|
|
827
|
|
|
|
|
|
|
|
828
|
|
|
|
|
|
|
|
829
|
|
|
|
|
|
|
|
830
|
|
|
|
|
|
|
*t_test_nev = \&PDL::t_test_nev; |
831
|
|
|
|
|
|
|
|
832
|
|
|
|
|
|
|
|
833
|
|
|
|
|
|
|
|
834
|
|
|
|
|
|
|
|
835
|
|
|
|
|
|
|
|
836
|
|
|
|
|
|
|
=head2 t_test_paired |
837
|
|
|
|
|
|
|
|
838
|
|
|
|
|
|
|
=for sig |
839
|
|
|
|
|
|
|
|
840
|
|
|
|
|
|
|
Signature: (a(n); b(n); float+ [o]t(); [o]d()) |
841
|
|
|
|
|
|
|
|
842
|
|
|
|
|
|
|
|
843
|
|
|
|
|
|
|
|
844
|
|
|
|
|
|
|
=for ref |
845
|
|
|
|
|
|
|
|
846
|
|
|
|
|
|
|
Paired sample t-test. |
847
|
|
|
|
|
|
|
|
848
|
|
|
|
|
|
|
=cut |
849
|
|
|
|
|
|
|
|
850
|
|
|
|
|
|
|
|
851
|
|
|
|
|
|
|
=for bad |
852
|
|
|
|
|
|
|
|
853
|
|
|
|
|
|
|
t_test_paired processes bad values. |
854
|
|
|
|
|
|
|
It will set the bad-value flag of all output piddles if the flag is set for any of the input piddles. |
855
|
|
|
|
|
|
|
|
856
|
|
|
|
|
|
|
|
857
|
|
|
|
|
|
|
=cut |
858
|
|
|
|
|
|
|
|
859
|
|
|
|
|
|
|
|
860
|
|
|
|
|
|
|
|
861
|
|
|
|
|
|
|
|
862
|
|
|
|
|
|
|
|
863
|
|
|
|
|
|
|
|
864
|
|
|
|
|
|
|
*t_test_paired = \&PDL::t_test_paired; |
865
|
|
|
|
|
|
|
|
866
|
|
|
|
|
|
|
|
867
|
|
|
|
|
|
|
|
868
|
|
|
|
|
|
|
|
869
|
|
|
|
|
|
|
=head2 binomial_test |
870
|
|
|
|
|
|
|
|
871
|
|
|
|
|
|
|
=for Sig |
872
|
|
|
|
|
|
|
|
873
|
|
|
|
|
|
|
Signature: (x(); n(); p_expected(); [o]p()) |
874
|
|
|
|
|
|
|
|
875
|
|
|
|
|
|
|
=for ref |
876
|
|
|
|
|
|
|
|
877
|
|
|
|
|
|
|
Binomial test. One-tailed significance test for two-outcome distribution. Given the number of successes, the number of trials, and the expected probability of success, returns the probability of getting this many or more successes. |
878
|
|
|
|
|
|
|
|
879
|
|
|
|
|
|
|
This function does NOT currently support bad value in the number of successes. |
880
|
|
|
|
|
|
|
|
881
|
|
|
|
|
|
|
=for usage |
882
|
|
|
|
|
|
|
|
883
|
|
|
|
|
|
|
Usage: |
884
|
|
|
|
|
|
|
|
885
|
|
|
|
|
|
|
# assume a fair coin, ie. 0.5 probablity of getting heads |
886
|
|
|
|
|
|
|
# test whether getting 8 heads out of 10 coin flips is unusual |
887
|
|
|
|
|
|
|
|
888
|
|
|
|
|
|
|
my $p = binomial_test( 8, 10, 0.5 ); # 0.0107421875. Yes it is unusual. |
889
|
|
|
|
|
|
|
|
890
|
|
|
|
|
|
|
=cut |
891
|
|
|
|
|
|
|
|
892
|
|
|
|
|
|
|
*binomial_test = \&PDL::binomial_test; |
893
|
|
|
|
|
|
|
sub PDL::binomial_test { |
894
|
0
|
|
|
0
|
0
|
0
|
my ($x, $n, $P) = @_; |
895
|
|
|
|
|
|
|
|
896
|
0
|
0
|
|
|
|
0
|
carp 'Please install PDL::GSL::CDF.' unless $CDF; |
897
|
0
|
0
|
|
|
|
0
|
carp 'This function does NOT currently support bad value in the number of successes.' if $x->badflag(); |
898
|
|
|
|
|
|
|
|
899
|
0
|
|
|
|
|
0
|
my $pdlx = pdl($x); |
900
|
0
|
|
|
|
|
0
|
$pdlx->badflag(1); |
901
|
0
|
|
|
|
|
0
|
$pdlx = $pdlx->setvaltobad(0); |
902
|
|
|
|
|
|
|
|
903
|
0
|
|
|
|
|
0
|
my $p = 1 - PDL::GSL::CDF::gsl_cdf_binomial_P( $pdlx - 1, $P, $n ); |
904
|
0
|
|
|
|
|
0
|
$p = $p->setbadtoval(1); |
905
|
0
|
|
|
|
|
0
|
$p->badflag(0); |
906
|
|
|
|
|
|
|
|
907
|
0
|
|
|
|
|
0
|
return $p; |
908
|
|
|
|
|
|
|
} |
909
|
|
|
|
|
|
|
|
910
|
|
|
|
|
|
|
|
911
|
|
|
|
|
|
|
=head1 METHODS |
912
|
|
|
|
|
|
|
|
913
|
|
|
|
|
|
|
=head2 rtable |
914
|
|
|
|
|
|
|
|
915
|
|
|
|
|
|
|
=for ref |
916
|
|
|
|
|
|
|
|
917
|
|
|
|
|
|
|
Reads either file or file handle*. Returns observation x variable pdl and var and obs ids if specified. Ids in perl @ ref to allow for non-numeric ids. Other non-numeric entries are treated as missing, which are filled with $opt{MISSN} then set to BAD*. Can specify num of data rows to read from top but not arbitrary range. |
918
|
|
|
|
|
|
|
|
919
|
|
|
|
|
|
|
*If passed handle, it will not be closed here. |
920
|
|
|
|
|
|
|
|
921
|
|
|
|
|
|
|
*PDL::Bad::setvaltobad only works consistently with the default TYPE double before PDL-2.4.4_04. |
922
|
|
|
|
|
|
|
|
923
|
|
|
|
|
|
|
=for options |
924
|
|
|
|
|
|
|
|
925
|
|
|
|
|
|
|
Default options (case insensitive): |
926
|
|
|
|
|
|
|
|
927
|
|
|
|
|
|
|
V => 1, # verbose. prints simple status |
928
|
|
|
|
|
|
|
TYPE => double, |
929
|
|
|
|
|
|
|
C_ID => 1, # boolean. file has col id. |
930
|
|
|
|
|
|
|
R_ID => 1, # boolean. file has row id. |
931
|
|
|
|
|
|
|
R_VAR => 0, # boolean. set to 1 if var in rows |
932
|
|
|
|
|
|
|
SEP => "\t", # can take regex qr// |
933
|
|
|
|
|
|
|
MISSN => -999, # this value treated as missing and set to BAD |
934
|
|
|
|
|
|
|
NROW => '', # set to read specified num of data rows |
935
|
|
|
|
|
|
|
|
936
|
|
|
|
|
|
|
=for usage |
937
|
|
|
|
|
|
|
|
938
|
|
|
|
|
|
|
Usage: |
939
|
|
|
|
|
|
|
|
940
|
|
|
|
|
|
|
Sample file diet.txt: |
941
|
|
|
|
|
|
|
|
942
|
|
|
|
|
|
|
uid height weight diet |
943
|
|
|
|
|
|
|
akw 72 320 1 |
944
|
|
|
|
|
|
|
bcm 68 268 1 |
945
|
|
|
|
|
|
|
clq 67 180 2 |
946
|
|
|
|
|
|
|
dwm 70 200 2 |
947
|
|
|
|
|
|
|
|
948
|
|
|
|
|
|
|
($data, $idv, $ido) = rtable 'diet.txt'; |
949
|
|
|
|
|
|
|
|
950
|
|
|
|
|
|
|
# By default prints out data info and @$idv index and element |
951
|
|
|
|
|
|
|
|
952
|
|
|
|
|
|
|
reading diet.txt for data and id... OK. |
953
|
|
|
|
|
|
|
data table as PDL dim o x v: PDL: Double D [4,3] |
954
|
|
|
|
|
|
|
0 height |
955
|
|
|
|
|
|
|
1 weight |
956
|
|
|
|
|
|
|
2 diet |
957
|
|
|
|
|
|
|
|
958
|
|
|
|
|
|
|
Another way of using it, |
959
|
|
|
|
|
|
|
|
960
|
|
|
|
|
|
|
$data = rtable( \*STDIN, {TYPE=>long} ); |
961
|
|
|
|
|
|
|
|
962
|
|
|
|
|
|
|
=cut |
963
|
|
|
|
|
|
|
|
964
|
|
|
|
|
|
|
sub rtable { |
965
|
|
|
|
|
|
|
# returns obs x var data matrix and var and obs ids |
966
|
3
|
|
|
3
|
1
|
17858
|
my ($src, $opt) = @_; |
967
|
|
|
|
|
|
|
|
968
|
3
|
|
|
|
|
4
|
my $fh_in; |
969
|
3
|
50
|
33
|
|
|
29
|
if ($src =~ /STDIN/ or ref $src eq 'GLOB') { $fh_in = $src } |
|
3
|
|
|
|
|
5
|
|
970
|
0
|
0
|
|
|
|
0
|
else { open $fh_in, $src or croak "$!" } |
971
|
|
|
|
|
|
|
|
972
|
3
|
|
|
|
|
11
|
my %opt = ( V => 1, |
973
|
|
|
|
|
|
|
TYPE => double, |
974
|
|
|
|
|
|
|
C_ID => 1, |
975
|
|
|
|
|
|
|
R_ID => 1, |
976
|
|
|
|
|
|
|
R_VAR => 0, |
977
|
|
|
|
|
|
|
SEP => "\t", |
978
|
|
|
|
|
|
|
MISSN => -999, |
979
|
|
|
|
|
|
|
NROW => '', |
980
|
|
|
|
|
|
|
); |
981
|
3
|
|
33
|
|
|
64
|
$opt and $opt{uc $_} = $opt->{$_} for (keys %$opt); |
982
|
3
|
50
|
|
|
|
11
|
$opt{V} and print STDERR "reading $src for data and id... "; |
983
|
|
|
|
|
|
|
|
984
|
3
|
|
|
|
|
5
|
local $PDL::undefval = $opt{MISSN}; |
985
|
|
|
|
|
|
|
|
986
|
3
|
|
|
|
|
15
|
my $id_c = []; # match declaration of $id_r for return purpose |
987
|
3
|
50
|
|
|
|
6
|
if ($opt{C_ID}) { |
988
|
3
|
|
|
|
|
15
|
chomp( $id_c = <$fh_in> ); |
989
|
3
|
|
|
|
|
42
|
my @entries = split $opt{SEP}, $id_c; |
990
|
3
|
50
|
|
|
|
9
|
$opt{R_ID} and shift @entries; |
991
|
3
|
|
|
|
|
7
|
$id_c = \@entries; |
992
|
|
|
|
|
|
|
} |
993
|
|
|
|
|
|
|
|
994
|
3
|
|
|
|
|
14
|
my ($c_row, $id_r, $data, @data) = (0, [], PDL->null, ); |
995
|
3
|
|
|
|
|
55
|
while (<$fh_in>) { |
996
|
132
|
|
|
|
|
97
|
chomp; |
997
|
132
|
|
|
|
|
272
|
my @entries = split /$opt{SEP}/, $_, -1; |
998
|
|
|
|
|
|
|
|
999
|
132
|
50
|
|
|
|
202
|
$opt{R_ID} and push @$id_r, shift @entries; |
1000
|
|
|
|
|
|
|
|
1001
|
|
|
|
|
|
|
# rudimentary check for numeric entry |
1002
|
132
|
100
|
66
|
|
|
131
|
for (@entries) { $_ = $opt{MISSN} unless defined $_ and /\d\b/ } |
|
669
|
|
|
|
|
1875
|
|
1003
|
|
|
|
|
|
|
|
1004
|
132
|
|
|
|
|
216
|
push @data, pdl( $opt{TYPE}, \@entries ); |
1005
|
132
|
|
|
|
|
1615
|
$c_row ++; |
1006
|
|
|
|
|
|
|
last |
1007
|
132
|
50
|
33
|
|
|
408
|
if $opt{NROW} and $c_row == $opt{NROW}; |
1008
|
|
|
|
|
|
|
} |
1009
|
|
|
|
|
|
|
# not explicitly closing $fh_in here in case it's passed from outside |
1010
|
|
|
|
|
|
|
# $fh_in will close by going out of scope if opened here. |
1011
|
|
|
|
|
|
|
|
1012
|
3
|
|
|
|
|
12
|
$data = pdl $opt{TYPE}, @data; |
1013
|
3
|
|
|
|
|
160
|
@data = (); |
1014
|
|
|
|
|
|
|
# rid of last col unless there is data there |
1015
|
|
|
|
|
|
|
$data = $data(0:$data->getdim(0)-2, )->sever |
1016
|
3
|
100
|
|
|
|
54
|
unless ( nelem $data(-1, )->where($data(-1, ) != $opt{MISSN}) ); |
1017
|
|
|
|
|
|
|
|
1018
|
3
|
|
|
|
|
453
|
my ($idv, $ido) = ($id_r, $id_c); |
1019
|
|
|
|
|
|
|
# var in columns instead of rows |
1020
|
3
|
50
|
|
|
|
22
|
$opt{R_VAR} == 0 |
1021
|
|
|
|
|
|
|
and ($data, $idv, $ido) = ($data->inplace->transpose, $id_c, $id_r); |
1022
|
|
|
|
|
|
|
|
1023
|
3
|
50
|
|
|
|
71
|
if ($opt{V}) { |
1024
|
0
|
|
|
|
|
0
|
print STDERR "OK.\ndata table as PDL dim o x v: " . $data->info . "\n"; |
1025
|
0
|
|
0
|
|
|
0
|
$idv and print STDERR "$_\t$$idv[$_]\n" for (0..$#$idv); |
1026
|
|
|
|
|
|
|
} |
1027
|
|
|
|
|
|
|
|
1028
|
3
|
|
|
|
|
56
|
$data = $data->setvaltobad( $opt{MISSN} ); |
1029
|
3
|
|
|
|
|
19
|
$data->check_badflag; |
1030
|
3
|
50
|
|
|
|
184
|
return wantarray? (@$idv? ($data, $idv, $ido) : ($data, $ido)) : $data; |
|
|
50
|
|
|
|
|
|
1031
|
|
|
|
|
|
|
} |
1032
|
|
|
|
|
|
|
|
1033
|
|
|
|
|
|
|
=head2 group_by |
1034
|
|
|
|
|
|
|
|
1035
|
|
|
|
|
|
|
Returns pdl reshaped according to the specified factor variable. Most useful when used in conjunction with other threading calculations such as average, stdv, etc. When the factor variable contains unequal number of cases in each level, the returned pdl is padded with bad values to fit the level with the most number of cases. This allows the subsequent calculation (average, stdv, etc) to return the correct results for each level. |
1036
|
|
|
|
|
|
|
|
1037
|
|
|
|
|
|
|
Usage: |
1038
|
|
|
|
|
|
|
|
1039
|
|
|
|
|
|
|
# simple case with 1d pdl and equal number of n in each level of the factor |
1040
|
|
|
|
|
|
|
|
1041
|
|
|
|
|
|
|
pdl> p $a = sequence 10 |
1042
|
|
|
|
|
|
|
[0 1 2 3 4 5 6 7 8 9] |
1043
|
|
|
|
|
|
|
|
1044
|
|
|
|
|
|
|
pdl> p $factor = $a > 4 |
1045
|
|
|
|
|
|
|
[0 0 0 0 0 1 1 1 1 1] |
1046
|
|
|
|
|
|
|
|
1047
|
|
|
|
|
|
|
pdl> p $a->group_by( $factor )->average |
1048
|
|
|
|
|
|
|
[2 7] |
1049
|
|
|
|
|
|
|
|
1050
|
|
|
|
|
|
|
# more complex case with threading and unequal number of n across levels in the factor |
1051
|
|
|
|
|
|
|
|
1052
|
|
|
|
|
|
|
pdl> p $a = sequence 10,2 |
1053
|
|
|
|
|
|
|
[ |
1054
|
|
|
|
|
|
|
[ 0 1 2 3 4 5 6 7 8 9] |
1055
|
|
|
|
|
|
|
[10 11 12 13 14 15 16 17 18 19] |
1056
|
|
|
|
|
|
|
] |
1057
|
|
|
|
|
|
|
|
1058
|
|
|
|
|
|
|
pdl> p $factor = qsort $a( ,0) % 3 |
1059
|
|
|
|
|
|
|
[ |
1060
|
|
|
|
|
|
|
[0 0 0 0 1 1 1 2 2 2] |
1061
|
|
|
|
|
|
|
] |
1062
|
|
|
|
|
|
|
|
1063
|
|
|
|
|
|
|
pdl> p $a->group_by( $factor ) |
1064
|
|
|
|
|
|
|
[ |
1065
|
|
|
|
|
|
|
[ |
1066
|
|
|
|
|
|
|
[ 0 1 2 3] |
1067
|
|
|
|
|
|
|
[10 11 12 13] |
1068
|
|
|
|
|
|
|
] |
1069
|
|
|
|
|
|
|
[ |
1070
|
|
|
|
|
|
|
[ 4 5 6 BAD] |
1071
|
|
|
|
|
|
|
[ 14 15 16 BAD] |
1072
|
|
|
|
|
|
|
] |
1073
|
|
|
|
|
|
|
[ |
1074
|
|
|
|
|
|
|
[ 7 8 9 BAD] |
1075
|
|
|
|
|
|
|
[ 17 18 19 BAD] |
1076
|
|
|
|
|
|
|
] |
1077
|
|
|
|
|
|
|
] |
1078
|
|
|
|
|
|
|
ARRAY(0xa2a4e40) |
1079
|
|
|
|
|
|
|
|
1080
|
|
|
|
|
|
|
# group_by supports perl factors, multiple factors |
1081
|
|
|
|
|
|
|
# returns factor labels in addition to pdl in array context |
1082
|
|
|
|
|
|
|
|
1083
|
|
|
|
|
|
|
pdl> p $a = sequence 12 |
1084
|
|
|
|
|
|
|
[0 1 2 3 4 5 6 7 8 9 10 11] |
1085
|
|
|
|
|
|
|
|
1086
|
|
|
|
|
|
|
pdl> $odd_even = [qw( e o e o e o e o e o e o )] |
1087
|
|
|
|
|
|
|
|
1088
|
|
|
|
|
|
|
pdl> $magnitude = [qw( l l l l l l h h h h h h )] |
1089
|
|
|
|
|
|
|
|
1090
|
|
|
|
|
|
|
pdl> ($a_grouped, $label) = $a->group_by( $odd_even, $magnitude ) |
1091
|
|
|
|
|
|
|
|
1092
|
|
|
|
|
|
|
pdl> p $a_grouped |
1093
|
|
|
|
|
|
|
[ |
1094
|
|
|
|
|
|
|
[ |
1095
|
|
|
|
|
|
|
[0 2 4] |
1096
|
|
|
|
|
|
|
[1 3 5] |
1097
|
|
|
|
|
|
|
] |
1098
|
|
|
|
|
|
|
[ |
1099
|
|
|
|
|
|
|
[ 6 8 10] |
1100
|
|
|
|
|
|
|
[ 7 9 11] |
1101
|
|
|
|
|
|
|
] |
1102
|
|
|
|
|
|
|
] |
1103
|
|
|
|
|
|
|
|
1104
|
|
|
|
|
|
|
pdl> p Dumper $label |
1105
|
|
|
|
|
|
|
$VAR1 = [ |
1106
|
|
|
|
|
|
|
[ |
1107
|
|
|
|
|
|
|
'e_l', |
1108
|
|
|
|
|
|
|
'o_l' |
1109
|
|
|
|
|
|
|
], |
1110
|
|
|
|
|
|
|
[ |
1111
|
|
|
|
|
|
|
'e_h', |
1112
|
|
|
|
|
|
|
'o_h' |
1113
|
|
|
|
|
|
|
] |
1114
|
|
|
|
|
|
|
]; |
1115
|
|
|
|
|
|
|
|
1116
|
|
|
|
|
|
|
|
1117
|
|
|
|
|
|
|
=cut |
1118
|
|
|
|
|
|
|
|
1119
|
|
|
|
|
|
|
*group_by = \&PDL::group_by; |
1120
|
|
|
|
|
|
|
sub PDL::group_by { |
1121
|
3
|
|
|
3
|
0
|
4269
|
my $p = shift; |
1122
|
3
|
|
|
|
|
4
|
my @factors = @_; |
1123
|
|
|
|
|
|
|
|
1124
|
3
|
100
|
|
|
|
8
|
if ( @factors == 1 ) { |
1125
|
2
|
|
|
|
|
2
|
my $factor = $factors[0]; |
1126
|
2
|
|
|
|
|
1
|
my $label; |
1127
|
2
|
50
|
|
|
|
5
|
if (ref $factor eq 'ARRAY') { |
1128
|
0
|
|
|
|
|
0
|
$label = _ordered_uniq($factor); |
1129
|
0
|
|
|
|
|
0
|
$factor = _array_to_pdl($factor); |
1130
|
|
|
|
|
|
|
} else { |
1131
|
2
|
|
|
|
|
6
|
my $perl_factor = [$factor->list]; |
1132
|
2
|
|
|
|
|
34
|
$label = _ordered_uniq($perl_factor); |
1133
|
|
|
|
|
|
|
} |
1134
|
|
|
|
|
|
|
|
1135
|
2
|
|
|
|
|
4
|
my $p_reshaped = _group_by_single_factor( $p, $factor ); |
1136
|
|
|
|
|
|
|
|
1137
|
2
|
100
|
|
|
|
29
|
return wantarray? ($p_reshaped, $label) : $p_reshaped; |
1138
|
|
|
|
|
|
|
} |
1139
|
|
|
|
|
|
|
|
1140
|
|
|
|
|
|
|
# make sure all are arrays instead of pdls |
1141
|
1
|
50
|
|
|
|
3
|
@factors = map { ref($_) eq 'PDL'? [$_->list] : $_ } @factors; |
|
2
|
|
|
|
|
7
|
|
1142
|
|
|
|
|
|
|
|
1143
|
1
|
|
|
|
|
2
|
my (@cells); |
1144
|
1
|
|
|
|
|
3
|
for my $ele (0 .. $#{$factors[0]}) { |
|
1
|
|
|
|
|
3
|
|
1145
|
10
|
|
|
|
|
8
|
my $c = join '_', map { $_->[$ele] } @factors; |
|
20
|
|
|
|
|
19
|
|
1146
|
10
|
|
|
|
|
12
|
push @cells, $c; |
1147
|
|
|
|
|
|
|
} |
1148
|
|
|
|
|
|
|
# get uniq cell labels (ref List::MoreUtils::uniq) |
1149
|
1
|
|
|
|
|
1
|
my %seen; |
1150
|
1
|
|
|
|
|
1
|
my @uniq_cells = grep {! $seen{$_}++ } @cells; |
|
10
|
|
|
|
|
14
|
|
1151
|
|
|
|
|
|
|
|
1152
|
1
|
|
|
|
|
3
|
my $flat_factor = _array_to_pdl( \@cells ); |
1153
|
|
|
|
|
|
|
|
1154
|
1
|
|
|
|
|
2
|
my $p_reshaped = _group_by_single_factor( $p, $flat_factor ); |
1155
|
|
|
|
|
|
|
|
1156
|
|
|
|
|
|
|
# get levels of each factor and reshape accordingly |
1157
|
1
|
|
|
|
|
2
|
my @levels; |
1158
|
1
|
|
|
|
|
2
|
for (@factors) { |
1159
|
2
|
|
|
|
|
2
|
my %uniq; |
1160
|
2
|
|
|
|
|
5
|
@uniq{ @$_ } = (); |
1161
|
2
|
|
|
|
|
4
|
push @levels, scalar keys %uniq; |
1162
|
|
|
|
|
|
|
} |
1163
|
|
|
|
|
|
|
|
1164
|
1
|
|
|
|
|
5
|
$p_reshaped = $p_reshaped->reshape( $p_reshaped->dim(0), @levels )->sever; |
1165
|
|
|
|
|
|
|
|
1166
|
|
|
|
|
|
|
# make labels for the returned data structure matching pdl structure |
1167
|
1
|
|
|
|
|
29
|
my @labels; |
1168
|
1
|
50
|
|
|
|
4
|
if (wantarray) { |
1169
|
1
|
|
|
|
|
2
|
for my $ifactor (0 .. $#levels) { |
1170
|
2
|
|
|
|
|
2
|
my @factor_label; |
1171
|
2
|
|
|
|
|
4
|
for my $ilevel (0 .. $levels[$ifactor]-1) { |
1172
|
4
|
|
|
|
|
4
|
my $i = $ifactor * $levels[$ifactor] + $ilevel; |
1173
|
4
|
|
|
|
|
5
|
push @factor_label, $uniq_cells[$i]; |
1174
|
|
|
|
|
|
|
} |
1175
|
2
|
|
|
|
|
3
|
push @labels, \@factor_label; |
1176
|
|
|
|
|
|
|
} |
1177
|
|
|
|
|
|
|
} |
1178
|
|
|
|
|
|
|
|
1179
|
1
|
50
|
|
|
|
8
|
return wantarray? ($p_reshaped, \@labels) : $p_reshaped; |
1180
|
|
|
|
|
|
|
} |
1181
|
|
|
|
|
|
|
|
1182
|
|
|
|
|
|
|
# get uniq cell labels (ref List::MoreUtils::uniq) |
1183
|
|
|
|
|
|
|
sub _ordered_uniq { |
1184
|
2
|
|
|
2
|
|
2
|
my $arr = shift; |
1185
|
|
|
|
|
|
|
|
1186
|
2
|
|
|
|
|
2
|
my %seen; |
1187
|
2
|
|
|
|
|
2
|
my @uniq = grep { ! $seen{$_}++ } @$arr; |
|
20
|
|
|
|
|
42
|
|
1188
|
|
|
|
|
|
|
|
1189
|
2
|
|
|
|
|
5
|
return \@uniq; |
1190
|
|
|
|
|
|
|
} |
1191
|
|
|
|
|
|
|
|
1192
|
|
|
|
|
|
|
sub _group_by_single_factor { |
1193
|
3
|
|
|
3
|
|
3
|
my $p = shift; |
1194
|
3
|
|
|
|
|
3
|
my $factor = shift; |
1195
|
|
|
|
|
|
|
|
1196
|
3
|
|
|
|
|
7
|
$factor = $factor->squeeze; |
1197
|
3
|
50
|
|
|
|
93
|
die "Currently support only 1d factor pdl." |
1198
|
|
|
|
|
|
|
if $factor->ndims > 1; |
1199
|
|
|
|
|
|
|
|
1200
|
3
|
50
|
|
|
|
11
|
die "Data pdl and factor pdl do not match!" |
1201
|
|
|
|
|
|
|
unless $factor->dim(0) == $p->dim(0); |
1202
|
|
|
|
|
|
|
|
1203
|
|
|
|
|
|
|
# get active dim that will be split according to factor and dims to thread over |
1204
|
3
|
|
|
|
|
6
|
my @p_threaddims = $p->dims; |
1205
|
3
|
|
|
|
|
34
|
my $p_dim0 = shift @p_threaddims; |
1206
|
|
|
|
|
|
|
|
1207
|
3
|
|
|
|
|
9
|
my $uniq = $factor->uniq; |
1208
|
|
|
|
|
|
|
|
1209
|
3
|
|
|
|
|
438
|
my @uniq_ns; |
1210
|
3
|
|
|
|
|
7
|
for ($uniq->list) { |
1211
|
9
|
|
|
|
|
205
|
push @uniq_ns, which( $factor == $_ )->nelem; |
1212
|
|
|
|
|
|
|
} |
1213
|
|
|
|
|
|
|
|
1214
|
|
|
|
|
|
|
# get number of n's in each group, find the biggest, fit output pdl to this |
1215
|
3
|
|
|
|
|
57
|
my $uniq_ns = pdl \@uniq_ns; |
1216
|
3
|
|
|
|
|
42
|
my $max = pdl(\@uniq_ns)->max; |
1217
|
|
|
|
|
|
|
|
1218
|
3
|
|
|
|
|
113
|
my $badvalue = int($p->max + 1); |
1219
|
3
|
|
|
|
|
73
|
my $p_tmp = ones($max, @p_threaddims, $uniq->nelem) * $badvalue; |
1220
|
3
|
|
|
|
|
160
|
for (0 .. $#uniq_ns) { |
1221
|
9
|
|
|
|
|
879
|
my $i = which $factor == $uniq($_); |
1222
|
9
|
|
|
|
|
299
|
$p_tmp->dice_axis(-1,$_)->squeeze->(0:$uniq_ns[$_]-1, ) .= $p($i, ); |
1223
|
|
|
|
|
|
|
} |
1224
|
|
|
|
|
|
|
|
1225
|
3
|
|
|
|
|
367
|
$p_tmp->badflag(1); |
1226
|
3
|
|
|
|
|
29
|
return $p_tmp->setvaltobad($badvalue); |
1227
|
|
|
|
|
|
|
} |
1228
|
|
|
|
|
|
|
|
1229
|
|
|
|
|
|
|
=head2 which_id |
1230
|
|
|
|
|
|
|
|
1231
|
|
|
|
|
|
|
=for ref |
1232
|
|
|
|
|
|
|
|
1233
|
|
|
|
|
|
|
Lookup specified var (obs) ids in $idv ($ido) (see B) and return indices in $idv ($ido) as pdl if found. The indices are ordered by the specified subset. Useful for selecting data by var (obs) id. |
1234
|
|
|
|
|
|
|
|
1235
|
|
|
|
|
|
|
=for usage |
1236
|
|
|
|
|
|
|
|
1237
|
|
|
|
|
|
|
my $ind = which_id $ido, ['smith', 'summers', 'tesla']; |
1238
|
|
|
|
|
|
|
|
1239
|
|
|
|
|
|
|
my $data_subset = $data( $ind, ); |
1240
|
|
|
|
|
|
|
|
1241
|
|
|
|
|
|
|
# take advantage of perl pattern matching |
1242
|
|
|
|
|
|
|
# e.g. use data from people whose last name starts with s |
1243
|
|
|
|
|
|
|
|
1244
|
|
|
|
|
|
|
my $i = which_id $ido, [ grep { /^s/ } @$ido ]; |
1245
|
|
|
|
|
|
|
|
1246
|
|
|
|
|
|
|
my $data_s = $data($i, ); |
1247
|
|
|
|
|
|
|
|
1248
|
|
|
|
|
|
|
=cut |
1249
|
|
|
|
|
|
|
|
1250
|
|
|
|
|
|
|
sub which_id { |
1251
|
38
|
|
|
38
|
1
|
39
|
my ($id, $id_s) = @_; |
1252
|
|
|
|
|
|
|
|
1253
|
38
|
|
|
|
|
28
|
my %ind; |
1254
|
38
|
|
|
|
|
107
|
@ind{ @$id } = ( 0 .. $#$id ); |
1255
|
|
|
|
|
|
|
|
1256
|
38
|
|
|
|
|
32
|
my @ind_select; |
1257
|
38
|
|
|
|
|
54
|
for (@$id_s) { |
1258
|
65
|
50
|
|
|
|
143
|
defined( $ind{$_} ) and push @ind_select, $ind{$_}; |
1259
|
|
|
|
|
|
|
} |
1260
|
38
|
|
|
|
|
63
|
return pdl @ind_select; |
1261
|
|
|
|
|
|
|
} |
1262
|
|
|
|
|
|
|
|
1263
|
|
|
|
|
|
|
sub _array_to_pdl { |
1264
|
133
|
|
|
133
|
|
2055
|
my ($var_ref) = @_; |
1265
|
|
|
|
|
|
|
|
1266
|
133
|
|
|
|
|
96
|
my (%level, $l); |
1267
|
133
|
|
|
|
|
109
|
$l = 0; |
1268
|
133
|
|
|
|
|
179
|
for (@$var_ref) { |
1269
|
4684
|
100
|
66
|
|
|
18100
|
if (defined($_) and $_ ne '' and $_ ne 'BAD') { |
|
|
|
100
|
|
|
|
|
1270
|
|
|
|
|
|
|
$level{$_} = $l ++ |
1271
|
4676
|
100
|
|
|
|
6785
|
if !exists $level{$_}; |
1272
|
|
|
|
|
|
|
} |
1273
|
|
|
|
|
|
|
} |
1274
|
|
|
|
|
|
|
|
1275
|
133
|
100
|
100
|
|
|
146
|
my $pdl = pdl( map { (defined($_) and $_ ne '' and $_ ne 'BAD')? $level{$_} : -1 } @$var_ref ); |
|
4684
|
|
|
|
|
17979
|
|
1276
|
133
|
|
|
|
|
3282
|
$pdl = $pdl->setvaltobad(-1); |
1277
|
133
|
|
|
|
|
836
|
$pdl->check_badflag; |
1278
|
|
|
|
|
|
|
|
1279
|
133
|
100
|
|
|
|
4538
|
return wantarray? ($pdl, \%level) : $pdl; |
1280
|
|
|
|
|
|
|
} |
1281
|
|
|
|
|
|
|
|
1282
|
|
|
|
|
|
|
|
1283
|
|
|
|
|
|
|
=head1 SEE ALSO |
1284
|
|
|
|
|
|
|
|
1285
|
|
|
|
|
|
|
PDL::Basic (hist for frequency counts) |
1286
|
|
|
|
|
|
|
|
1287
|
|
|
|
|
|
|
PDL::Ufunc (sum, avg, median, min, max, etc.) |
1288
|
|
|
|
|
|
|
|
1289
|
|
|
|
|
|
|
PDL::GSL::CDF (various cumulative distribution functions) |
1290
|
|
|
|
|
|
|
|
1291
|
|
|
|
|
|
|
=head1 REFERENCES |
1292
|
|
|
|
|
|
|
|
1293
|
|
|
|
|
|
|
Hays, W.L. (1994). Statistics (5th ed.). Fort Worth, TX: Harcourt Brace College Publishers. |
1294
|
|
|
|
|
|
|
|
1295
|
|
|
|
|
|
|
=head1 AUTHOR |
1296
|
|
|
|
|
|
|
|
1297
|
|
|
|
|
|
|
Copyright (C) 2009 Maggie J. Xiong |
1298
|
|
|
|
|
|
|
|
1299
|
|
|
|
|
|
|
All rights reserved. There is no warranty. You are allowed to redistribute this software / documentation as described in the file COPYING in the PDL distribution. |
1300
|
|
|
|
|
|
|
|
1301
|
|
|
|
|
|
|
=cut |
1302
|
|
|
|
|
|
|
|
1303
|
|
|
|
|
|
|
|
1304
|
|
|
|
|
|
|
|
1305
|
|
|
|
|
|
|
; |
1306
|
|
|
|
|
|
|
|
1307
|
|
|
|
|
|
|
|
1308
|
|
|
|
|
|
|
|
1309
|
|
|
|
|
|
|
# Exit with OK status |
1310
|
|
|
|
|
|
|
|
1311
|
|
|
|
|
|
|
1; |
1312
|
|
|
|
|
|
|
|
1313
|
|
|
|
|
|
|
|