Skip to content

Commit 46336b5

Browse files
mohawk2mohawk2
committedNov 22, 2024·
use Test::PDL in tests for Basic and GLM
·
0.8550.850
1 parent 602de55 commit 46336b5

File tree

3 files changed

+255
-377
lines changed

3 files changed

+255
-377
lines changed
 

‎Makefile.PL

+1
Original file line numberDiff line numberDiff line change
@@ -32,6 +32,7 @@ WriteMakefile(
3232
},
3333
TEST_REQUIRES => {
3434
'Test::More' => '0.88', # done_testing
35+
'Test::PDL' => '0.21',
3536
},
3637
$got_PDL ? () : (DIR => []), # just write MYMETA if no PDL
3738
dist => { PREOP => 'gsl-config --version && $(PERL) -MPDL::Core::Dev -e pdlpp_mkgen $(DISTVNAME)' },

‎t/stats_basic.t

+85-114
Original file line numberDiff line numberDiff line change
@@ -5,120 +5,113 @@ use Test::More;
55
use PDL::LiteF;
66
use PDL::NiceSlice;
77
use PDL::Stats::Basic;
8-
9-
sub tapprox {
10-
my($a,$b, $eps) = @_;
11-
$eps ||= 1e-6;
12-
my $diff = abs($a-$b);
13-
# use max to make it perl scalar
14-
ref $diff eq 'PDL' and $diff = $diff->max;
15-
return $diff < $eps;
16-
}
8+
use Test::PDL;
179

1810
my $a = sequence 5;
19-
20-
is( tapprox( $a->stdv, 1.4142135623731 ), 1, "standard deviation of $a");
21-
is( tapprox( $a->stdv_unbiased, 1.58113883008419 ), 1, "unbiased standard deviation of $a");
22-
is( tapprox( $a->var, 2 ), 1, "variance of $a");
23-
is( tapprox( $a->var_unbiased, 2.5 ), 1, "unbiased variance of $a");
24-
is( tapprox( $a->se, 0.707106781186548 ), 1, "standard error of $a");
25-
is( tapprox( $a->ss, 10 ), 1, "sum of squared deviations from the mean of $a");
26-
is( tapprox( $a->skew, 0 ), 1, "sample skewness of $a");
27-
is( tapprox( $a->skew_unbiased, 0 ), 1, "unbiased sample skewness of $a");
28-
is( tapprox( $a->kurt, -1.3 ), 1, "sample kurtosis of $a");
29-
is( tapprox( $a->kurt_unbiased, -1.2 ), 1, "unbiased sample kurtosis of $a");
30-
31-
{
32-
ok(tapprox($_->ss, (($_ - $_->avg)**2)->sum), "ss for $_") for
33-
pdl('[1 1 1 1 2 3 4 4 4 4 4 4]'),
34-
pdl('[1 2 2 2 3 3 3 3 4 4 5 5]'),
35-
pdl('[1 1 1 2 2 3 3 4 4 5 5 5]');
36-
}
11+
is_pdl $a->stdv, pdl( 1.4142135623731 ), "standard deviation of $a";
12+
is_pdl $a->stdv_unbiased, pdl( 1.58113883008419 ), "unbiased standard deviation of $a";
13+
is_pdl $a->var, pdl( 2 ), "variance of $a";
14+
is_pdl $a->var_unbiased, pdl( 2.5 ), "unbiased variance of $a";
15+
is_pdl $a->se, pdl( 0.707106781186548 ), "standard error of $a";
16+
is_pdl $a->ss, pdl( 10 ), "sum of squared deviations from the mean of $a";
17+
is_pdl $a->skew, pdl( 0 ), "sample skewness of $a";
18+
is_pdl $a->skew_unbiased, pdl( 0 ), "unbiased sample skewness of $a";
19+
is_pdl $a->kurt, pdl( -1.3 ), "sample kurtosis of $a";
20+
is_pdl $a->kurt_unbiased, pdl( -1.2 ), "unbiased sample kurtosis of $a";
21+
22+
is_pdl $_->ss, (($_ - $_->avg)**2)->sumover, "ss for $_" for
23+
pdl('[1 1 1 1 2 3 4 4 4 4 4 4]'),
24+
pdl('[1 2 2 2 3 3 3 3 4 4 5 5]'),
25+
pdl('[1 1 1 2 2 3 3 4 4 5 5 5]');
3726

3827
my $a_bad = sequence 6;
3928
$a_bad->setbadat(-1);
40-
41-
is( tapprox( $a_bad->stdv, 1.4142135623731 ), 1, "standard deviation of $a_bad");
42-
is( tapprox( $a_bad->stdv_unbiased, 1.58113883008419 ), 1, "unbiased standard deviation of $a_bad");
43-
is( tapprox( $a_bad->var, 2 ), 1, "variance of $a_bad");
44-
is( tapprox( $a_bad->var_unbiased, 2.5 ), 1, "unbiased variance of $a_bad");
45-
is( tapprox( $a_bad->se, 0.707106781186548 ), 1, "standard error of $a_bad");
46-
is( tapprox( $a_bad->ss, 10 ), 1, "sum of squared deviations from the mean of $a_bad");
47-
is( tapprox( $a_bad->skew, 0 ), 1, "sample skewness of $a_bad");
48-
is( tapprox( $a_bad->skew_unbiased, 0 ), 1, "unbiased sample skewness of $a_bad");
49-
is( tapprox( $a_bad->kurt, -1.3 ), 1, "sample kurtosis of $a_bad");
50-
is( tapprox( $a_bad->kurt_unbiased, -1.2 ), 1, "unbiased sample kurtosis of $a_bad");
51-
52-
my $b = sequence 5;
53-
$b %= 2;
54-
$b = qsort $b;
55-
56-
is( tapprox( $a->cov($b), 0.6 ), 1, "sample covariance of $a and $b" );
57-
is( tapprox( $a->corr($b), 0.866025403784439 ), 1, "Pearson correlation coefficient of $a and $b");
58-
is( tapprox( $a->n_pair($b), 5 ), 1, "Number of good pairs between $a and $b");
59-
is( tapprox( $a->corr($b)->t_corr( 5 ), 3 ), 1, "t significance test of Pearson correlation coefficient of $a and $b");
60-
is( tapprox( $a->corr_dev($b), 0.903696114115064 ), 1, "correlation calculated from dev_m values of $a and $b");
61-
62-
my $b_bad = sequence 6;
63-
$b_bad = qsort( $b_bad % 2 );
64-
$b_bad->setbadat(0);
65-
66-
is( tapprox( $a_bad->cov($b_bad), 0.5 ), 1, "sample covariance with bad data of $a_bad and $b_bad");
67-
is( tapprox( $a_bad->corr($b_bad), 0.894427190999916 ), 1, "Pearson correlation coefficient with bad data of $a_bad and $b_bad");
68-
is( tapprox( $a_bad->n_pair($b_bad), 4 ), 1, "Number of good pairs between $a_bad and $b_bad with bad values taken into account");
69-
is( tapprox( $a_bad->corr($b_bad)->t_corr( 4 ), 2.82842712474619 ), 1, "t signifiance test of Pearson correlation coefficient with bad data of $a_bad and $b_bad");
70-
is( tapprox( $a_bad->corr_dev($b_bad), 0.903696114115064 ), 1, "correlation calculated from dev_m values with bad data of $a_bad and $b_bad");
29+
is_pdl $a_bad->stdv, pdl( 1.4142135623731 ), "standard deviation of $a_bad";
30+
is_pdl $a_bad->stdv_unbiased, pdl( 1.58113883008419 ), "unbiased standard deviation of $a_bad";
31+
is_pdl $a_bad->var, pdl( 2 ), "variance of $a_bad";
32+
is_pdl $a_bad->var_unbiased, pdl( 2.5 ), "unbiased variance of $a_bad";
33+
is_pdl $a_bad->se, pdl( 0.707106781186548 ), "standard error of $a_bad";
34+
is_pdl $a_bad->ss, pdl( 10 ), "sum of squared deviations from the mean of $a_bad";
35+
is_pdl $a_bad->skew, pdl( 0 ), "sample skewness of $a_bad";
36+
is_pdl $a_bad->skew_unbiased, pdl( 0 ), "unbiased sample skewness of $a_bad";
37+
is_pdl $a_bad->kurt, pdl( -1.3 ), "sample kurtosis of $a_bad";
38+
is_pdl $a_bad->kurt_unbiased, pdl( -1.2 ), "unbiased sample kurtosis of $a_bad";
39+
40+
my $b = pdl '0 0 0 1 1';
41+
is_pdl $a->cov($b), pdl( 0.6 ), "sample covariance of $a and $b";
42+
is_pdl $a->corr($b), pdl( 0.866025403784439 ), "Pearson correlation coefficient of $a and $b";
43+
is_pdl $a->n_pair($b), indx( 5 ), "Number of good pairs between $a and $b";
44+
is_pdl $a->corr($b)->t_corr( 5 ), pdl( 3 ), "t significance test of Pearson correlation coefficient of $a and $b";
45+
is_pdl $a->corr_dev($b), pdl( 0.903696114115064 ), "correlation calculated from dev_m values of $a and $b";
46+
47+
my $b_bad = pdl 'BAD 0 0 1 1 1';
48+
is_pdl $a_bad->cov($b_bad), pdl( 0.5 ), "sample covariance with bad data of $a_bad and $b_bad";
49+
is_pdl $a_bad->corr($b_bad), pdl( 0.894427190999916 ), "Pearson correlation coefficient with bad data of $a_bad and $b_bad";
50+
is_pdl $a_bad->n_pair($b_bad), indx( 4 ), "Number of good pairs between $a_bad and $b_bad with bad values taken into account";
51+
is_pdl $a_bad->corr($b_bad)->t_corr( 4 ), pdl( 2.82842712474619 ), "t signifiance test of Pearson correlation coefficient with bad data of $a_bad and $b_bad";
52+
is_pdl $a_bad->corr_dev($b_bad), pdl( 0.903696114115064 ), "correlation calculated from dev_m values with bad data of $a_bad and $b_bad";
7153

7254
my ($t, $df) = $a->t_test($b);
73-
is( tapprox( $t, 2.1380899352994 ), 1, "t-test between $a and $b - 't' output");
74-
is( tapprox( $df, 8 ), 1, "t-test between $a and $b - 'df' output");
55+
is_pdl $t, pdl( 2.1380899352994 ), "t-test between $a and $b - 't' output";
56+
is_pdl $df, pdl( 8 ), "t-test between $a and $b - 'df' output";
7557

7658
($t, $df) = $a->t_test_nev($b);
77-
is( tapprox( $t, 2.1380899352994 ), 1, "t-test with non-equal variance between $a and $b - 't' output");
78-
is( tapprox( $df, 4.94637223974763 ), 1, "t-test with non-equal variance between $a and $b - 'df' output");
59+
is_pdl $t, pdl( 2.1380899352994 ), "t-test with non-equal variance between $a and $b - 't' output";
60+
is_pdl $df, pdl( 4.94637223974763 ), "t-test with non-equal variance between $a and $b - 'df' output";
7961

8062
($t, $df) = $a->t_test_paired($b);
81-
is( tapprox( $t, 3.13785816221094 ), 1, "paired sample t-test between $a and $b - 't' output");
82-
is( tapprox( $df, 4 ), 1, "paired sample t-test between $a and $b - 'df' output");
63+
is_pdl $t, pdl( 3.13785816221094 ), "paired sample t-test between $a and $b - 't' output";
64+
is_pdl $df, pdl( 4 ), "paired sample t-test between $a and $b - 'df' output";
8365

8466
($t, $df) = $a_bad->t_test($b_bad);
85-
is( tapprox( $t, 1.87082869338697 ), 1, "t-test with bad values between $a_bad and $b_bad - 't' output");
86-
is( tapprox( $df, 8 ), 1, "t-test with bad values between $a_bad and $b_bad - 'd' output");
67+
is_pdl $t, pdl( 1.87082869338697 ), "t-test with bad values between $a_bad and $b_bad - 't' output";
68+
is_pdl $df, pdl( 8 ), "t-test with bad values between $a_bad and $b_bad - 'd' output";
8769

8870
($t, $df) = $a_bad->t_test_nev($b_bad);
89-
is( tapprox( $t, 1.87082869338697 ), 1, "t-test with non-equal variance with bad values between $a_bad and $b_bad - 't' output");
90-
is( tapprox( $df, 4.94637223974763 ), 1, "t-test with non-equal variance with bad values between $a_bad and $b_bad - 'df' output");
71+
is_pdl $t, pdl( 1.87082869338697 ), "t-test with non-equal variance with bad values between $a_bad and $b_bad - 't' output";
72+
is_pdl $df, pdl( 4.94637223974763 ), "t-test with non-equal variance with bad values between $a_bad and $b_bad - 'df' output";
9173

9274
($t, $df) = $a_bad->t_test_paired($b_bad);
93-
is( tapprox( $t, 4.89897948556636 ), 1, "paired sample t-test with bad values between $a_bad and $b_bad - 't' output");
94-
is( tapprox( $df, 3 ), 1, "paired sample t-test with bad values between $a_bad and $b_bad - 'df' output");
75+
is_pdl $t, pdl( 4.89897948556636 ), "paired sample t-test with bad values between $a_bad and $b_bad - 't' output";
76+
is_pdl $df, pdl( 3 ), "paired sample t-test with bad values between $a_bad and $b_bad - 'df' output";
9577

9678
{
9779
my ($data, $idv, $ido) = rtable(\*DATA, {V=>0});
98-
is( tapprox( sum(pdl($data->dims) - pdl(14, 5)), 0 ), 1, 'rtable data dim' );
99-
is( tapprox( $data->sum / $data->nbad, 1.70731707317073 ), 1, 'rtable bad elem' );
80+
is_pdl $data, pdl '
81+
[ 5 BAD BAD 2 BAD 5 BAD 9 4 BAD BAD BAD 5 BAD]
82+
[ 7 BAD 3 7 0 BAD 0 8 BAD 0 3 0 BAD 0]
83+
[BAD BAD BAD BAD BAD 1 BAD 1 BAD BAD BAD BAD 1 BAD]
84+
[BAD BAD BAD BAD BAD 0 BAD 5 BAD BAD BAD BAD 0 BAD]
85+
[BAD BAD 0 BAD 2 BAD 0 BAD BAD 0 0 2 BAD 0]
86+
';
10087
}
10188

10289
{
103-
my $a = random 10, 3;
104-
is( tapprox( sum($a->cov_table - $a->cov($a->dummy(1))), 0 ), 1, 'cov_table' );
105-
90+
my $a = pdl '
91+
0.045 0.682 0.290 0.024 0.598 0.321 0.772 0.375 0.237 0.811;
92+
0.356 0.094 0.925 0.139 0.701 0.849 0.689 0.109 0.240 0.847;
93+
0.822 0.492 0.351 0.860 0.400 0.243 0.313 0.011 0.437 0.480
94+
';
95+
is_pdl $a->cov_table, $a->cov($a->dummy(1)), 'cov_table';
10696
$a->setbadat(4,0);
107-
is( tapprox( sum($a->cov_table - $a->cov($a->dummy(1))), 0 ), 1, 'cov_table bad val' );
97+
is_pdl $a->cov_table, $a->cov($a->dummy(1)), 'cov_table bad val';
10898
}
10999

110100
{
111-
my $a = random 10, 3;
112-
is( tapprox( sum(abs($a->corr_table - $a->corr($a->dummy(1)))), 0 ), 1, "Square Pearson correlation table");
113-
101+
my $a = pdl '
102+
0.045 0.682 0.290 0.024 0.598 0.321 0.772 0.375 0.237 0.811;
103+
0.356 0.094 0.925 0.139 0.701 0.849 0.689 0.109 0.240 0.847;
104+
0.822 0.492 0.351 0.860 0.400 0.243 0.313 0.011 0.437 0.480
105+
';
106+
is_pdl $a->corr_table, $a->corr($a->dummy(1)), "Square Pearson correlation table";
114107
$a->setbadat(4,0);
115-
is( tapprox( sum(abs($a->corr_table - $a->corr($a->dummy(1)))), 0 ), 1, "Square Pearson correlation table with bad data");
108+
is_pdl $a->corr_table, $a->corr($a->dummy(1)), "Square Pearson correlation table with bad data";
116109
}
117110

118111
{
119112
my $a = pdl([0,1,2,3,4], [0,0,0,0,0]);
120113
$a = $a->setvaltobad(0);
121-
is( $a->stdv->nbad, 1, "Bad value input to stdv makes the stdv itself bad");
114+
ok $a->stdv->nbad, "Bad value input to stdv makes the stdv itself bad";
122115
}
123116

124117
SKIP: {
@@ -127,67 +120,45 @@ SKIP: {
127120
my $x = pdl(1, 2);
128121
my $n = pdl(2, 10);
129122
my $p = .5;
130-
131-
my $a = pdl qw[ 0.75 0.9892578125 ];
132-
133-
is (tapprox( sum(abs(binomial_test( $x,$n,$p ) - $a)) ,0), 1, 'binomial_test');
123+
is_pdl binomial_test( $x,$n,$p ), pdl(0.75, 0.9892578125), 'binomial_test';
134124
}
135125

136126
{
137127
my $a = sequence 10, 2;
138128
my $factor = sequence(10) > 4;
139129
my $ans = pdl( [[0..4], [10..14]], [[5..9], [15..19]] );
140-
141130
my ($a_, $l) = $a->group_by($factor);
142-
is( tapprox( sum(abs($a_ - $ans)), 0 ), 1, 'group_by single factor equal n' );
143-
is_deeply( $l, [0, 1], 'group_by single factor label');
131+
is_pdl $a_, $ans, 'group_by single factor equal n';
132+
is_deeply $l, [0, 1], 'group_by single factor label';
144133

145134
$a = sequence 10,2;
146135
$factor = qsort sequence(10) % 3;
147136
$ans = pdl( [1.5, 11.5], [5, 15], [8, 18] );
148-
149-
is( tapprox( sum(abs($a->group_by($factor)->average - $ans)), 0 ), 1, 'group_by single factor unequal n' );
137+
is_pdl $a->group_by($factor)->average, $ans, 'group_by single factor unequal n';
150138

151139
$a = sequence 10;
152140
my @factors = ( [qw( a a a a b b b b b b )], [qw(0 1 0 1 0 1 0 1 0 1)] );
153-
$ans = pdl(
154-
[
155-
[0,2,-1],
156-
[1,3,-1],
157-
],
158-
[
159-
[4,6,8],
160-
[5,7,9],
161-
]
162-
);
163-
$ans->badflag(1);
164-
$ans = $ans->setvaltobad(-1);
165-
141+
$ans = pdl '[ 0 2 BAD; 1 3 BAD ], [ 4 6 8; 5 7 9 ]';
166142
($a_, $l) = $a->group_by( @factors );
167-
is(tapprox(sum(abs($a_ - $ans)), 0), 1, 'group_by multiple factors') or diag($a_, $ans);
168-
is_deeply($l, [[qw(a_0 a_1)], [qw( b_0 b_1 )]], 'group_by multiple factors label');
143+
is_pdl $a_, $ans, 'group_by multiple factors';
144+
is_deeply $l, [[qw(a_0 a_1)], [qw( b_0 b_1 )]], 'group_by multiple factors label';
169145
}
170146

171147

172148
{
173149
my @a = qw(a a b b c c);
174150
my $a = PDL::Stats::Basic::_array_to_pdl( \@a );
175151
my $ans = pdl( 0,0,1,1,2,2 );
176-
is( tapprox( sum(abs($a - $ans)), 0 ), 1, '_array_to_pdl' );
152+
is_pdl $a, $ans, '_array_to_pdl';
177153

178154
$a[-1] = undef;
179155
my $a_bad = PDL::Stats::Basic::_array_to_pdl( \@a );
180-
my $ans_bad = pdl( 0,0,1,1,2,2 );
181-
$ans_bad = $ans_bad->setbadat(-1);
182-
183-
like( $a_bad(-1)->isbad(), qr/1/, '_array_to_pdl with missing value undef' );
184-
is( tapprox( sum(abs($a_bad - $ans_bad)), 0 ), 1, '_array_to_pdl with missing value undef correctly coded' );
156+
my $ans_bad = pdl '0 0 1 1 2 BAD';
157+
is_pdl $a_bad, $ans_bad, '_array_to_pdl with missing value undef correctly coded';
185158

186159
$a[-1] = 'BAD';
187160
$a_bad = PDL::Stats::Basic::_array_to_pdl( \@a );
188-
189-
like( $a_bad(-1)->isbad(), qr/1/, '_array_to_pdl with missing value BAD' );
190-
is( tapprox( sum(abs($a_bad - $ans_bad)), 0 ), 1, '_array_to_pdl with missing value BAD correctly coded' );
161+
is_pdl $a_bad, $ans_bad, '_array_to_pdl with missing value BAD correctly coded';
191162
}
192163

193164
done_testing();

‎t/stats_glm.t

+169-263
Original file line numberDiff line numberDiff line change
@@ -5,47 +5,30 @@ use PDL::Stats::Basic;
55
use PDL::Stats::GLM;
66
use PDL::LiteF;
77
use PDL::NiceSlice;
8+
use Test::PDL;
89

9-
sub tapprox {
10-
my($a,$b, $eps) = @_;
11-
$eps ||= 1e-6;
12-
my $diff = abs($a-$b);
13-
# use max to make it perl scalar
14-
ref $diff eq 'PDL' and $diff = $diff->max;
15-
return $diff < $eps;
16-
}
17-
18-
my $a = sequence 5;
19-
my $b = pdl(0, 0, 0, 1, 1);
20-
21-
is( t_fill_m(), 1, "fill_m replaces bad values with sample mean");
22-
sub t_fill_m {
23-
my $aa = sequence 5;
24-
$aa = $aa->setvaltobad(0);
25-
tapprox( $aa->fill_m->sum, 12.5 );
26-
}
10+
is_pdl pdl('BAD 1 2 3 4')->fill_m, pdl('2.5 1 2 3 4'), "fill_m replaces bad values with sample mean";
2711

28-
is( t_fill_rand(), 1, "fill_rand replaces bad values with random sample of good values from same variable");
29-
sub t_fill_rand {
30-
my $aa = sequence 5;
31-
$aa = $aa->setvaltobad(0);
32-
my $stdv = $aa->fill_rand->stdv;
33-
tapprox( $stdv, 1.01980390271856 ) || tapprox( $stdv, 1.16619037896906 );
12+
{
13+
my $stdv = pdl('BAD 1 2 3 4')->fill_rand->stdv;
14+
ok PDL::Core::approx( $stdv, 1.01980390271856 ) || PDL::Core::approx( $stdv, 1.16619037896906 ), "fill_rand replaces bad values with random sample of good values from same variable";
3415
}
3516

36-
ok tapprox( $a->dev_m->avg, 0 ), "dev_m replaces values with deviations from the mean on $a";
37-
ok tapprox( $a->stddz->avg, 0 ), "stddz standardizes data on $a";
17+
my $a = sequence 5;
18+
is_pdl $a->dev_m, pdl('-2 -1 0 1 2'), "dev_m replaces values with deviations from the mean on $a";
19+
is_pdl $a->stddz, pdl('-1.41421356237309 -0.707106781186547 0 0.707106781186547 1.41421356237309'), "stddz standardizes data on $a";
3820

39-
ok tapprox( $a->sse($b), 18), "sse gives sum of squared errors between actual and predicted values between $a and $b";
40-
ok tapprox( $a->mse($b), 3.6), "mse gives mean of squared errors between actual and predicted values between $a and $b";
41-
ok tapprox( $a->rmse($b), 1.89736659610103 ), "rmse gives root mean squared error, ie. stdv around predicted value between $a and $b";
21+
my $b = pdl(0, 0, 0, 1, 1);
22+
is_pdl $a->sse($b), pdl(18), "sse gives sum of squared errors between actual and predicted values between $a and $b";
23+
is_pdl $a->mse($b), pdl(3.6), "mse gives mean of squared errors between actual and predicted values between $a and $b";
24+
is_pdl $a->rmse($b), pdl(1.89736659610103), "rmse gives root mean squared error, ie. stdv around predicted value between $a and $b";
4225

43-
ok tapprox( $b->glue(1,ones(5))->pred_logistic(pdl(1,2))->sum, 4.54753948757851 ), "pred_logistic calculates predicted probability value for logistic regression";
26+
is_pdl $b->glue(1,ones(5))->pred_logistic(pdl(1,2)), pdl('0.880797077977882 0.880797077977882 0.880797077977882 0.952574126822433 0.952574126822433'), "pred_logistic calculates predicted probability value for logistic regression";
4427

4528
my $y = pdl(0, 1, 0, 1, 0);
46-
ok tapprox( $y->d0(), 6.73011667009256 ), 'd0';
47-
ok tapprox( $y->dm( ones(5) * .5 ), 6.93147180559945 ), 'dm';
48-
ok tapprox( sum($y->dvrs(ones(5) * .5) ** 2), 6.93147180559945 ), 'dvrs';
29+
is_pdl $y->d0(), pdl( 6.73011667009256 ), 'd0';
30+
is_pdl $y->dm( ones(5) * .5 ), pdl( 6.93147180559945 ), 'dm';
31+
is_pdl $y->dvrs(ones(5) * .5) ** 2, pdl('1.38629436111989 1.38629436111989 1.38629436111989 1.38629436111989 1.38629436111989'), 'dvrs';
4932

5033
{
5134
my $a = pdl(ushort, [0,0,1,0,1], [0,0,0,1,1] );
@@ -67,17 +50,15 @@ ok tapprox( sum($y->dvrs(ones(5) * .5) ** 2), 6.93147180559945 ), 'dvrs';
6750
[qw( 0.0071428571 0.035714286 -0.057142857)],
6851
],
6952
);
70-
ok tapprox( sum( abs($m{R2} - $rsq) ), 0 ), 'ols_t R2';
71-
ok tapprox( sum( abs($m{b} - $coeff) ), 0 ), 'ols_t b';
53+
is_pdl $m{R2}, $rsq, 'ols_t R2';
54+
is_pdl $m{b}, $coeff, 'ols_t b';
7255

7356
my %m0 = $a->ols_t(sequence(5), {CONST=>0});
7457
my $b0 = pdl ([ 0.2 ], [ 0.23333333 ]);
75-
76-
ok tapprox( sum( abs($m0{b} - $b0) ), 0 ), 'ols_t, const=>0';
58+
is_pdl $m0{b}, $b0, 'ols_t, const=>0';
7759
}
7860

79-
ok tapprox( t_ols(), 0 ), 'ols';
80-
sub t_ols {
61+
{
8162
my $a = sequence 5;
8263
my $b = pdl(0,0,0,1,1);
8364
my %m = $a->ols($b, {plot=>0});
@@ -94,13 +75,11 @@ sub t_ols {
9475
test_stats_cmp(\%m, \%a);
9576
}
9677

97-
ok tapprox( t_ols_bad(), 0 ), 'ols with bad value';
98-
sub t_ols_bad {
99-
my $a = sequence 6;
78+
{
79+
my $a = pdl '0 1 2 3 4 BAD';
10080
my $b = pdl(0,0,0,1,1,1);
101-
$a->setbadat(5);
10281
my %m = $a->ols($b, {plot=>0});
103-
is( $b->sumover, 3, "ols with bad value didn't change caller value" );
82+
is_pdl $b, pdl(0,0,0,1,1,1), "ols with bad value didn't change caller value";
10483
ok $a->check_badflag, "ols with bad value didn't remove caller bad flag";
10584
my %a = (
10685
F => 9,
@@ -115,8 +94,7 @@ sub t_ols_bad {
11594
test_stats_cmp(\%m, \%a);
11695
}
11796

118-
ok tapprox( t_r2_change(), 0 ), 'r2_change';
119-
sub t_r2_change {
97+
{
12098
my $a = sequence 5, 2;
12199
my $b = pdl(0,0,0,1,1);
122100
my $c = pdl(0,0,2,2,2);
@@ -138,16 +116,16 @@ R2_change => pdl(.15, .15),
138116

139117
my %p = $a->pca({CORR=>1, PLOT=>0});
140118
my %a = (
141-
eigenvalue => pdl( qw( 2.786684 0.18473727 0.028578689) ),
119+
eigenvalue => float( qw( 2.786684 0.18473727 0.028578689) ),
142120
# loadings in R
143-
eigenvector => [pdl(
121+
eigenvector => [float(
144122
# v1 v2 v3
145123
[qw( 0.58518141 0.58668657 0.55978709)], # comp1
146124
[qw( -0.41537629 -0.37601061 0.82829859)], # comp2
147125
[qw( -0.69643754 0.71722722 -0.023661276)], # comp3
148126
), \&PDL::abs],
149127

150-
loadings => [pdl(
128+
loadings => [float(
151129
[qw( 0.97686463 0.97937725 0.93447296)],
152130
[qw( -0.17853319 -0.1616134 0.35601163)],
153131
[qw( -0.11773439 0.12124893 -0.0039999937)],
@@ -159,8 +137,8 @@ pct_var => pdl( qw(0.92889468 0.06157909 0.0095262297) ),
159137

160138
%p = $a->pca({CORR=>0, PLOT=>0});
161139
%a = (
162-
eigenvalue => [pdl(qw[ 22.0561695 1.581758022 0.202065959 ]), \&PDL::abs],
163-
eigenvector => [pdl(
140+
eigenvalue => [float(qw[ 22.0561695 1.581758022 0.202065959 ]), \&PDL::abs],
141+
eigenvector => [float(
164142
[qw(-0.511688 -0.595281 -0.619528)],
165143
[qw( 0.413568 0.461388 -0.78491)],
166144
[qw( 0.753085 -0.657846 0.0101023)],
@@ -177,87 +155,67 @@ pct_var => pdl( qw[0.925175 0.0663489 0.00847592] ),
177155
test_stats_cmp(\%p, \%a, 1e-4);
178156
}
179157

180-
ok tapprox( t_pca_sorti(), 0 ), "pca_sorti - principal component analysis output sorted to find which vars a component is best represented";
181-
sub t_pca_sorti {
182-
my $a = sequence 10, 5;
183-
$a = lvalue_assign_detour( $a, which($a % 7 == 0), 0 );
184-
158+
{
159+
# pca_sorti - principal component analysis output sorted to find which vars a component is best represented
160+
my $a = pdl '
161+
0 1 2 3 4 5 6 0 8 9; 10 11 12 13 0 15 16 17 18 19;
162+
20 0 22 23 24 25 26 27 0 29; 30 31 32 33 34 0 36 37 38 39;
163+
40 41 0 43 44 45 46 47 48 0
164+
';
185165
my %m = $a->pca({PLOT=>0});
186-
187166
my ($iv, $ic) = $m{loadings}->pca_sorti;
188-
189-
return sum($iv - pdl(qw(4 1 0 2 3))) + sum($ic - pdl(qw( 0 1 2 )));
167+
is_pdl $iv, indx(qw(4 1 0 2 3));
168+
is_pdl $ic, pdl(qw( 0 1 2 ));
190169
}
191170

192171
SKIP: {
193172
eval { require PDL::Fit::LM; };
194173
skip 'no PDL::Fit::LM', 1 if $@;
195-
196-
ok tapprox( t_logistic(), 0 ), 'logistic';
197-
198-
my $y = pdl( 0, 0, 0, 1, 1 );
199-
my $x = pdl(2, 3, 5, 5, 5);
200-
my %m = $y->logistic( $x, {COV=>1} );
201-
isnt $m{cov}, undef, 'get cov from logistic if ask';
202-
};
203-
sub t_logistic {
204174
my $y = pdl( 0, 0, 0, 1, 1 );
205175
my $x = pdl(2, 3, 5, 5, 5);
206176
my %m = $y->logistic( $x );
207177
my $y_pred = $x->glue(1, ones(5))->pred_logistic( $m{b} );
208178
my $y_pred_ans
209179
= pdl qw(7.2364053e-07 0.00010154254 0.66666667 0.66666667 0.66666667);
210-
return sum( $y_pred - $y_pred_ans, $m{Dm_chisq} - 2.91082711764867 );
211-
}
212-
213-
my $a_bad = sequence 6;
214-
$a_bad->setbadat(-1);
215-
my $b_bad = pdl(0, 0, 0, 0, 1, 1);
216-
$b_bad->setbadat(0);
217-
218-
ok tapprox( $a_bad->dev_m->avg, 0 ), "dev_m with bad values $a_bad";
219-
ok tapprox( $a_bad->stddz->avg, 0 ), "stdz with bad values $a_bad";
220-
221-
ok tapprox( $a_bad->sse($b_bad), 23), "sse with bad values between $a_bad and $b_bad";
222-
ok tapprox( $a_bad->mse($b_bad), 5.75), "mse with badvalues between $a_bad and $b_bad";
223-
ok tapprox( $a_bad->rmse($b_bad), 2.39791576165636 ), "rmse with bad values between $a_bad and $b_bad";
224-
225-
ok tapprox( $b_bad->glue(1,ones(6))->pred_logistic(pdl(1,2))->sum, 4.54753948757851 ), "pred_logistic with bad values";
180+
is_pdl $y_pred, $y_pred_ans;
181+
is_pdl $m{Dm_chisq}, pdl 2.91082711764867;
182+
%m = $y->logistic( $x, {COV=>1} );
183+
isnt $m{cov}, undef, 'get cov from logistic if ask';
184+
};
226185

227-
ok tapprox( $b_bad->d0(), 6.73011667009256 ), "null deviance with bad values on $b_bad";
228-
ok tapprox( $b_bad->dm( ones(6) * .5 ), 6.93147180559945 ), "model deviance with bad values on $b_bad";
229-
ok tapprox( sum($b_bad->dvrs(ones(6) * .5) ** 2), 6.93147180559945 ), "deviance residual with bad values on $b_bad";
186+
my $a_bad = pdl '0 1 2 3 4 BAD';
187+
my $b_bad = pdl 'BAD 0 0 0 1 1';
188+
is_pdl $a_bad->dev_m, pdl( '-2 -1 0 1 2 BAD' ), "dev_m with bad values $a_bad";
189+
is_pdl $a_bad->stddz, pdl( '-1.41421356237309 -0.707106781186547 0 0.707106781186547 1.41421356237309 BAD' ), "stdz with bad values $a_bad";
190+
is_pdl $a_bad->sse($b_bad), pdl(23), "sse with bad values between $a_bad and $b_bad";
191+
is_pdl $a_bad->mse($b_bad), pdl(5.75), "mse with badvalues between $a_bad and $b_bad";
192+
is_pdl $a_bad->rmse($b_bad), pdl( 2.39791576165636 ), "rmse with bad values between $a_bad and $b_bad";
193+
is_pdl $b_bad->glue(1,ones(6))->pred_logistic(pdl(1,2)), pdl( 'BAD 0.880797077977882 0.880797077977882 0.880797077977882 0.952574126822433 0.952574126822433' ), "pred_logistic with bad values";
194+
is_pdl $b_bad->d0(), pdl( 6.73011667009256 ), "null deviance with bad values on $b_bad";
195+
is_pdl $b_bad->dm( ones(6) * .5 ), pdl( 6.93147180559945 ), "model deviance with bad values on $b_bad";
196+
is_pdl $b_bad->dvrs(ones(6) * .5), pdl( 'BAD -1.17741002251547 -1.17741002251547 -1.17741002251547 1.17741002251547 1.17741002251547' ), "deviance residual with bad values on $b_bad";
230197

231198
{
232199
eval { effect_code(['a']) };
233200
isnt $@, '', 'effect_code with only one value dies';
234-
my @a = qw( a a a b b b b c c BAD );
235-
my $a = effect_code(\@a);
236-
my $ans = pdl [
237-
[qw( 1 1 1 0 0 0 0 -1 -1 -99 )],
238-
[qw( 0 0 0 1 1 1 1 -1 -1 -99 )]
239-
];
240-
$ans = $ans->setvaltobad(-99);
241-
is( sum(abs(which($a->isbad) - pdl(9,19))), 0, 'effect_code got bad value' );
242-
ok tapprox( sum(abs($a - $ans)), 0 ), 'effect_code coded with bad value';
201+
my $a = scalar effect_code([qw(a a a b b b b c c BAD)]);
202+
is_pdl $a, pdl('1 1 1 0 0 0 0 -1 -1 BAD; 0 0 0 1 1 1 1 -1 -1 BAD'), 'effect_code coded with bad value';
243203
}
244204

245-
ok tapprox( t_effect_code_w(), 0 ), 'effect_code_w';
246-
sub t_effect_code_w {
205+
{
247206
eval { effect_code_w(['a']) };
248207
isnt $@, '', 'effect_code_w with only one value dies';
249-
my @a = qw( a a a b b b b c c c );
250-
my $a = effect_code_w(\@a);
251-
return sum($a->sumover - pdl byte, (0, 0));
208+
is_pdl scalar effect_code_w([qw(a a a b b b b c c c)]), pdl '
209+
1 1 1 0 0 0 0 -1 -1 -1; 0 0 0 1 1 1 1 -1.3333333 -1.3333333 -1.3333333
210+
';
252211
}
253212

254-
ok tapprox( t_anova(), 0 ), 'anova_3w';
255-
sub t_anova {
213+
{ # anova 3 way
256214
my $d = sequence 60;
257215
my @a = map {$a = $_; map { $a } 0..14 } qw(a b c d);
258216
my $b = $d % 3;
259217
my $c = $d % 2;
260-
$d = lvalue_assign_detour( $d, 20, 10 );
218+
$d->set( 20, 10 );
261219
my %m = $d->anova(\@a, $b, $c, {IVNM=>[qw(A B C)], plot=>0});
262220
$m{'# A ~ B ~ C # m'} = $m{'# A ~ B ~ C # m'}->(,2,)->squeeze;
263221
test_stats_cmp(\%m, {
@@ -267,8 +225,7 @@ sub t_anova {
267225
});
268226
}
269227

270-
ok tapprox( t_anova_1way(), 0 ), 'anova_1w';
271-
sub t_anova_1way {
228+
{ # anova 1 way
272229
my $d = pdl qw( 3 2 1 5 2 1 5 3 1 4 1 2 3 5 5 );
273230
my $a = qsort sequence(15) % 3;
274231
my %m = $d->anova($a, {plot=>0});
@@ -280,10 +237,9 @@ sub t_anova_1way {
280237
});
281238
}
282239

283-
ok tapprox( t_anova_bad_dv(), 0 ), 'anova_3w bad dv';
284-
sub t_anova_bad_dv {
240+
{ # anova_3w bad dv
285241
my $d = sequence 60;
286-
$d = lvalue_assign_detour( $d, 20, 10 );
242+
$d->set( 20, 10 );
287243
$d->setbadat(1);
288244
$d->setbadat(10);
289245
my @a = map {$a = $_; map { $a } 0..14 } qw(a b c d);
@@ -299,14 +255,13 @@ sub t_anova_bad_dv {
299255
});
300256
}
301257

302-
ok tapprox( t_anova_bad_dv_iv(), 0 ), 'anova_3w bad dv iv';
303-
sub t_anova_bad_dv_iv {
258+
{ # anova_3w bad dv iv
304259
my $d = sequence 63;
305260
my @a = map {$a = $_; map { $a } 0..14 } qw(a b c d);
306261
push @a, undef, qw( b c );
307262
my $b = $d % 3;
308263
my $c = $d % 2;
309-
$d = lvalue_assign_detour( $d, 20, 10 );
264+
$d->set( 20, 10 );
310265
$d->setbadat(62);
311266
$b->setbadat(61);
312267
my %m = $d->anova(\@a, $b, $c, {IVNM=>[qw(A B C)], plot=>0});
@@ -318,19 +273,10 @@ sub t_anova_bad_dv_iv {
318273
});
319274
}
320275

321-
{
322-
my $a = pdl([0,1,2,3,4], [0,0,0,0,0]);
323-
$a = $a->setvaltobad(0);
324-
is( $a->fill_m->setvaltobad(0)->nbad, 5, 'fill_m nan to bad');
325-
}
276+
is_pdl pdl('BAD 1 2 3 4; BAD BAD BAD BAD BAD')->fill_m, pdl('2.5 1 2 3 4; 0 0 0 0 0'), 'fill_m nan to bad';
277+
is_pdl pdl([1,1,1], [2,2,2])->stddz, zeroes(3,2), 'stddz nan vs bad';
326278

327279
{
328-
my $a = pdl([1,1,1], [2,2,2]);
329-
is( which($a->stddz == 0)->nelem, 6, 'stddz nan vs bad');
330-
}
331-
332-
ok tapprox( t_anova_rptd_basic(), 0 ), 'anova_rptd_basic';
333-
sub t_anova_rptd_basic {
334280
# data from https://www.youtube.com/watch?v=Fh73dAOMm9M
335281
# Person,Before,After 2 weeks,After 4 weeks
336282
# P1,102,97,95
@@ -378,8 +324,7 @@ sub t_anova_rptd_basic {
378324
});
379325
}
380326

381-
ok tapprox( t_anova_rptd_1way(), 0 ), 'anova_rptd_1w';
382-
sub t_anova_rptd_1way {
327+
{ # anova_rptd_1w
383328
my $d = pdl qw( 3 2 1 5 2 1 5 3 1 4 1 2 3 5 5 );
384329
my $s = sequence(5)->dummy(1,3)->flat;
385330
my $a = qsort sequence(15) % 3;
@@ -392,14 +337,13 @@ sub t_anova_rptd_1way {
392337
});
393338
}
394339

395-
ok tapprox( t_anova_rptd_2way_bad_dv(), 0 ), 'anova_rptd_2w bad dv';
396340
my %anova_bad_a = (
397341
'| a | F' => 0.351351351351351,
398342
'| a | ms' => 0.722222222222222,
399343
'| a ~ b | F' => 5.25,
400344
'# a ~ b # m' => pdl(qw( 3 1.3333333 3.3333333 3.3333333 3.6666667 2.6666667 ))->reshape(3,2),
401345
);
402-
sub t_anova_rptd_2way_bad_dv {
346+
{ # anova_rptd_2w bad dv
403347
my $d = pdl qw( 3 2 1 5 2 1 5 3 1 4 1 2 3 5 5 3 4 2 1 5 4 3 2 2);
404348
$d = $d->setbadat(5);
405349
my $s = sequence(4)->dummy(1,6)->flat;
@@ -412,8 +356,7 @@ sub t_anova_rptd_2way_bad_dv {
412356
test_stats_cmp(\%m, \%anova_bad_a);
413357
}
414358

415-
ok tapprox( t_anova_rptd_2way_bad_iv(), 0 ), 'anova_rptd_2w bad iv';
416-
sub t_anova_rptd_2way_bad_iv {
359+
{ # anova_rptd_2w bad iv
417360
my $d = pdl qw( 3 2 1 5 2 1 5 3 1 4 1 2 3 5 5 3 4 2 1 5 4 3 2 2);
418361
my $s = sequence(4)->dummy(1,6)->flat;
419362
# [0 1 2 3 0 1 2 3 0 1 2 3 0 1 2 3 0 1 2 3 0 1 2 3]
@@ -426,8 +369,7 @@ sub t_anova_rptd_2way_bad_iv {
426369
test_stats_cmp(\%m, \%anova_bad_a);
427370
}
428371

429-
ok tapprox( t_anova_rptd_3way(), 0 ), 'anova_rptd_3w';
430-
sub t_anova_rptd_3way {
372+
{ # anova_rptd_3w
431373
my $d = pdl( qw( 3 2 1 5 2 1 5 3 1 4 1 2 3 5 5 3 4 2 1 5 4 3 2 2 ),
432374
qw( 5 5 1 1 4 4 1 4 4 2 3 3 5 1 1 2 4 4 4 5 5 1 1 2 )
433375
);
@@ -447,8 +389,7 @@ sub t_anova_rptd_3way {
447389
});
448390
}
449391

450-
ok tapprox( t_anova_rptd_mixed(), 0 ), 'anova_rptd mixed';
451-
sub t_anova_rptd_mixed {
392+
{ # anova_rptd mixed
452393
my $d = pdl qw( 3 2 1 5 2 1 5 3 1 4 1 2 3 5 5 3 4 2 1 5 4 3 2 2);
453394
my $s = sequence(4)->dummy(1,6)->flat;
454395
# [0 1 2 3 0 1 2 3 0 1 2 3 0 1 2 3 0 1 2 3 0 1 2 3]
@@ -469,130 +410,112 @@ sub t_anova_rptd_mixed {
469410

470411
# Tests for mixed anova thanks to Erich Greene
471412

472-
ok tapprox( t_anova_rptd_mixed_l2ord2(), 0, ), 'anova_rptd mixed with 2 btwn-subj var levels, data grouped by subject';
473-
SKIP: {
474-
skip "yet to be fixed", 3;
475-
ok tapprox( t_anova_rptd_mixed_l2ord1(), 0, ), 'anova_rptd mixed with 2 btwn-subj var levels, data grouped by within var';
476-
ok tapprox( t_anova_rptd_mixed_l3ord1(), 0, .001 ), 'anova_rptd mixed with 3 btwn-subj var levels, data grouped by within var';
477-
ok tapprox( t_anova_rptd_mixed_l3ord2(), 0, .001 ), 'anova_rptd mixed with 3 btwn-subj var levels, data grouped by subject';
478-
}
479413
sub test_stats_cmp {
480414
local $Test::Builder::Level = $Test::Builder::Level + 1;
481415
my ($m, $ans, $eps) = @_;
482416
$eps ||= 1e-6;
483-
my $error = pdl 0;
484417
foreach (sort keys %$ans) {
485418
my $got = PDL->topdl($m->{$_});
486419
my $exp = $ans->{$_};
487420
if (ref $exp eq 'ARRAY') {
488421
($exp, my $func) = @$exp;
489422
($got, $exp) = map &$func($_), $got, $exp;
490423
}
491-
$exp = PDL->topdl($exp);
492-
$error = $error + (my $this_diff = $got - $exp);
493-
fail($_), diag "got $m->{$_}\nexpected $exp" if any($this_diff->abs > $eps);
424+
is_pdl $got, PDL->topdl($exp), {atol=>$eps, test_name=>$_};
494425
}
495-
return $error;
496-
}
497-
sub t_anova_rptd_mixed_backend {
498-
my ($d,$s,$w,$b,$ans) = @_;
499-
my %m = $d->anova_rptd($s,$w,$b,{ivnm=>['within','between'],btwn=>[1],plot=>0, v=>0});
500-
test_stats_cmp(\%m, $ans);
501-
}
502-
sub t_anova_rptd_mixed_l2_common {
503-
my ($d,$s,$w,$b) = @_;
504-
my %ans = (
505-
'| within | df' => 2,
506-
'| within || err df' => 12,
507-
'| within | ss' => .25,
508-
'| within | ms' => .125,
509-
'| within || err ss' => 23.666667,
510-
'| within || err ms' => 1.9722222,
511-
'| within | F' => 0.063380282,
512-
'| between | df' => 1,
513-
'| between || err df' => 6,
514-
'| between | ss' => 2.0416667,
515-
'| between | ms' => 2.0416667,
516-
'| between || err ss' => 16.583333,
517-
'| between || err ms' => 2.7638889,
518-
'| between | F' => 0.73869347,
519-
'| within ~ between | df' => 2,
520-
'| within ~ between | ss' => 6.0833333,
521-
'| within ~ between | ms' => 3.0416667,
522-
'| within ~ between | F' => 1.5422535,
523-
);
524-
$ans{"| within ~ between || err $_"} = $ans{"| within || err $_"} foreach qw/df ss ms/;
525-
return t_anova_rptd_mixed_backend($d,$s,$w,$b,\%ans);
526-
}
527-
sub t_anova_rptd_mixed_l3_common {
528-
my ($d,$s,$w,$b) = @_;
529-
my %ans = (
530-
'| within | df' => 2,
531-
'| within || err df' => 12,
532-
'| within | ss' => .963,
533-
'| within | ms' => .481,
534-
'| within || err ss' => 20.889,
535-
'| within || err ms' => 1.741,
536-
'| within | F' => .277,
537-
'| between | df' => 2,
538-
'| between || err df' => 6,
539-
'| between | ss' => 1.185,
540-
'| between | ms' => .593,
541-
'| between || err ss' => 13.111,
542-
'| between || err ms' => 2.185,
543-
'| between | F' => .271,
544-
'| within ~ between | df' => 4,
545-
'| within ~ between | ss' => 4.148,
546-
'| within ~ between | ms' => 1.037,
547-
'| within ~ between | F' => .596,
548-
);
549-
$ans{"| within ~ between || err $_"} = $ans{"| within || err $_"} foreach qw/df ss ms/;
550-
return t_anova_rptd_mixed_backend($d,$s,$w,$b,\%ans);
551-
}
552-
sub t_anova_rptd_mixed_l2ord1 {
553-
my $d = pdl qw( 3 2 1 5 2 1 5 3 1 4 1 2 3 5 5 3 4 2 1 5 4 3 2 2);
554-
my $s = sequence(8)->dummy(1,3)->flat;
555-
# [0 1 2 3 4 5 6 7 0 1 2 3 4 5 6 7 0 1 2 3 4 5 6 7]
556-
my $w = qsort sequence(24) % 3;
557-
# [0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2]
558-
my $b = (sequence(8) % 2)->qsort->dummy(1,3)->flat;
559-
# [0 0 0 0 1 1 1 1 0 0 0 0 1 1 1 1 0 0 0 0 1 1 1 1]
560-
return t_anova_rptd_mixed_l2_common($d,$s,$w,$b);
561426
}
562-
sub t_anova_rptd_mixed_l2ord2 {
563-
my $d = pdl qw( 3 1 4 2 4 2 1 1 1 5 2 5 2 3 4 1 5 3 5 5 2 3 3 2);
564-
my $s = qsort sequence(24) % 8;
565-
# [0 0 0 1 1 1 2 2 2 3 3 3 4 4 4 5 5 5 6 6 6 7 7 7]
566-
my $w = sequence(3)->dummy(1,8)->flat;
567-
# [0 1 2 0 1 2 0 1 2 0 1 2 0 1 2 0 1 2 0 1 2 0 1 2]
568-
my $b = qsort sequence(24) % 2;
569-
# [0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1]
570-
return t_anova_rptd_mixed_l2_common($d,$s,$w,$b);
571-
}
572-
sub t_anova_rptd_mixed_l3ord1 {
573-
my $d = pdl qw( 5 2 2 5 4 1 5 3 5 4 4 3 4 3 4 3 5 1 4 3 3 4 5 4 5 5 2 );
574-
my $s = sequence(9)->dummy(1,3)->flat;
575-
# [0 1 2 3 4 5 6 7 8 0 1 2 3 4 5 6 7 8 0 1 2 3 4 5 6 7 8]
576-
my $w = qsort sequence(27) % 3;
577-
# [0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2]
578-
my $b = (sequence(9) % 3)->qsort->dummy(1,3)->flat;
579-
# [0 0 0 1 1 1 2 2 2 0 0 0 1 1 1 2 2 2 0 0 0 1 1 1 2 2 2]
580-
return t_anova_rptd_mixed_l3_common($d,$s,$w,$b);
581-
}
582-
sub t_anova_rptd_mixed_l3ord2 {
583-
my $d = pdl qw( 5 4 4 2 4 3 2 3 3 5 4 4 4 3 5 1 4 4 5 3 5 3 5 5 5 1 2 );
584-
my $s = qsort sequence(27) % 9;
585-
# [0 0 0 1 1 1 2 2 2 3 3 3 4 4 4 5 5 5 6 6 6 7 7 7 8 8 8]
586-
my $w = sequence(3)->dummy(1,9)->flat;
587-
# [0 1 2 0 1 2 0 1 2 0 1 2 0 1 2 0 1 2 0 1 2 0 1 2 0 1 2]
588-
my $b = qsort sequence(27) % 3;
589-
# [0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2]
590-
return t_anova_rptd_mixed_l3_common($d,$s,$w,$b);
427+
my %anova_ans_l2_common = (
428+
'| within | df' => 2,
429+
'| within || err df' => 12,
430+
'| within | ss' => .25,
431+
'| within | ms' => .125,
432+
'| within || err ss' => 23.666667,
433+
'| within || err ms' => 1.9722222,
434+
'| within | F' => 0.063380282,
435+
'| between | df' => 1,
436+
'| between || err df' => 6,
437+
'| between | ss' => 2.0416667,
438+
'| between | ms' => 2.0416667,
439+
'| between || err ss' => 16.583333,
440+
'| between || err ms' => 2.7638889,
441+
'| between | F' => 0.73869347,
442+
'| within ~ between | df' => 2,
443+
'| within ~ between | ss' => 6.0833333,
444+
'| within ~ between | ms' => 3.0416667,
445+
'| within ~ between | F' => 1.5422535,
446+
);
447+
$anova_ans_l2_common{"| within ~ between || err $_"} = $anova_ans_l2_common{"| within || err $_"} foreach qw/df ss ms/;
448+
my %anova_ans_l3_common = (
449+
'| within | df' => 2,
450+
'| within || err df' => 12,
451+
'| within | ss' => .963,
452+
'| within | ms' => .481,
453+
'| within || err ss' => 20.889,
454+
'| within || err ms' => 1.741,
455+
'| within | F' => .277,
456+
'| between | df' => 2,
457+
'| between || err df' => 6,
458+
'| between | ss' => 1.185,
459+
'| between | ms' => .593,
460+
'| between || err ss' => 13.111,
461+
'| between || err ms' => 2.185,
462+
'| between | F' => .271,
463+
'| within ~ between | df' => 4,
464+
'| within ~ between | ss' => 4.148,
465+
'| within ~ between | ms' => 1.037,
466+
'| within ~ between | F' => .596,
467+
);
468+
$anova_ans_l3_common{"| within ~ between || err $_"} = $anova_ans_l3_common{"| within || err $_"} foreach qw/df ss ms/;
469+
if (0) { # FIXME
470+
# anova_rptd mixed with 2 btwn-subj var levels, data grouped by within var
471+
my $d = pdl qw( 3 2 1 5 2 1 5 3 1 4 1 2 3 5 5 3 4 2 1 5 4 3 2 2);
472+
my $s = sequence(8)->dummy(1,3)->flat;
473+
# [0 1 2 3 4 5 6 7 0 1 2 3 4 5 6 7 0 1 2 3 4 5 6 7]
474+
my $w = qsort sequence(24) % 3;
475+
# [0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2]
476+
my $b = (sequence(8) % 2)->qsort->dummy(1,3)->flat;
477+
# [0 0 0 0 1 1 1 1 0 0 0 0 1 1 1 1 0 0 0 0 1 1 1 1]
478+
my %m = $d->anova_rptd($s,$w,$b,{ivnm=>['within','between'],btwn=>[1],plot=>0, v=>0});
479+
test_stats_cmp(\%m, \%anova_ans_l2_common);
591480
}
592-
593-
594-
ok tapprox( t_anova_rptd_mixed_bad(), 0 ), 'anova_rptd mixed bad';
595-
sub t_anova_rptd_mixed_bad {
481+
{
482+
# anova_rptd mixed with 2 btwn-subj var levels, data grouped by subject
483+
my $d = pdl qw( 3 1 4 2 4 2 1 1 1 5 2 5 2 3 4 1 5 3 5 5 2 3 3 2);
484+
my $s = qsort sequence(24) % 8;
485+
# [0 0 0 1 1 1 2 2 2 3 3 3 4 4 4 5 5 5 6 6 6 7 7 7]
486+
my $w = sequence(3)->dummy(1,8)->flat;
487+
# [0 1 2 0 1 2 0 1 2 0 1 2 0 1 2 0 1 2 0 1 2 0 1 2]
488+
my $b = qsort sequence(24) % 2;
489+
# [0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1]
490+
my %m = $d->anova_rptd($s,$w,$b,{ivnm=>['within','between'],btwn=>[1],plot=>0, v=>0});
491+
test_stats_cmp(\%m, \%anova_ans_l2_common);
492+
}
493+
if (0) { # FIXME
494+
# eps=.001 anova_rptd mixed with 3 btwn-subj var levels, data grouped by within var
495+
my $d = pdl qw( 5 2 2 5 4 1 5 3 5 4 4 3 4 3 4 3 5 1 4 3 3 4 5 4 5 5 2 );
496+
my $s = sequence(9)->dummy(1,3)->flat;
497+
# [0 1 2 3 4 5 6 7 8 0 1 2 3 4 5 6 7 8 0 1 2 3 4 5 6 7 8]
498+
my $w = qsort sequence(27) % 3;
499+
# [0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2]
500+
my $b = (sequence(9) % 3)->qsort->dummy(1,3)->flat;
501+
# [0 0 0 1 1 1 2 2 2 0 0 0 1 1 1 2 2 2 0 0 0 1 1 1 2 2 2]
502+
my %m = $d->anova_rptd($s,$w,$b,{ivnm=>['within','between'],btwn=>[1],plot=>0, v=>0});
503+
test_stats_cmp(\%m, \%anova_ans_l3_common);
504+
}
505+
if (0) { # FIXME
506+
# eps=.001 anova_rptd mixed with 3 btwn-subj var levels, data grouped by subject
507+
my $d = pdl qw( 5 4 4 2 4 3 2 3 3 5 4 4 4 3 5 1 4 4 5 3 5 3 5 5 5 1 2 );
508+
my $s = qsort sequence(27) % 9;
509+
# [0 0 0 1 1 1 2 2 2 3 3 3 4 4 4 5 5 5 6 6 6 7 7 7 8 8 8]
510+
my $w = sequence(3)->dummy(1,9)->flat;
511+
# [0 1 2 0 1 2 0 1 2 0 1 2 0 1 2 0 1 2 0 1 2 0 1 2 0 1 2]
512+
my $b = qsort sequence(27) % 3;
513+
# [0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2]
514+
my %m = $d->anova_rptd($s,$w,$b,{ivnm=>['within','between'],btwn=>[1],plot=>0, v=>0});
515+
test_stats_cmp(\%m, \%anova_ans_l3_common);
516+
}
517+
518+
{ # anova_rptd mixed bad
596519
my $d = pdl qw( 3 2 1 5 2 1 5 3 1 4 1 2 3 5 5 3 4 2 1 5 4 3 2 2 1 1 1 1 );
597520
my $s = sequence(4)->dummy(1,6)->flat;
598521
# [0 1 2 3 0 1 2 3 0 1 2 3 0 1 2 3 0 1 2 3 0 1 2 3]
@@ -617,8 +540,7 @@ sub t_anova_rptd_mixed_bad {
617540
});
618541
}
619542

620-
ok tapprox( t_anova_rptd_mixed_4w(), 0 ), 'anova_rptd_mixed_4w';
621-
sub t_anova_rptd_mixed_4w {
543+
{ # anova_rptd_mixed_4w
622544
my ($data, $idv, $subj) = rtable \*DATA, {v=>0};
623545
my ($age, $aa, $beer, $wings, $dv) = $data->dog;
624546
my %m = $dv->anova_rptd( $subj, $age, $aa, $beer, $wings, { ivnm=>[qw(age aa beer wings)], btwn=>[0,1], v=>0, plot=>0 } );
@@ -637,29 +559,13 @@ sub t_anova_rptd_mixed_4w {
637559
my $a = effect_code( sequence(12) > 5 );
638560
my $b = effect_code([ map {(0, 1)} (1..6) ]);
639561
my $c = effect_code([ map {(0,0,1,1,2,2)} (1..2) ]);
640-
641-
my $ans = pdl [
642-
[qw( 1 -1 0 -0 -1 1 -1 1 -0 0 1 -1 )],
643-
[qw( 0 -0 1 -1 -1 1 -0 0 -1 1 1 -1 )]
644-
];
562+
my $ans = pdl '1 -1 0 -0 -1 1 -1 1 -0 0 1 -1; 0 -0 1 -1 -1 1 -0 0 -1 1 1 -1';
645563
my $inter = interaction_code( $a, $b, $c);
646-
647-
is(sum(abs($inter - $ans)), 0, 'interaction_code');
564+
is_pdl $inter, $ans, 'interaction_code';
648565
}
649566

650567
done_testing();
651568

652-
sub lvalue_assign_detour {
653-
my ($pdl, $index, $new_value) = @_;
654-
655-
my @arr = list $pdl;
656-
my @ind = ref($index)? list($index) : $index;
657-
$arr[$_] = $new_value
658-
for (@ind);
659-
660-
return pdl(\@arr)->reshape($pdl->dims)->sever;
661-
}
662-
663569
__DATA__
664570
subj age Apple-android beer wings recall
665571
1 0 0 0 0 5

0 commit comments

Comments
 (0)
Please sign in to comment.