File Coverage

blib/lib/Statistics/Cook.pm
Criterion Covered Total %
statement 34 124 27.4
branch 4 20 20.0
condition 0 8 0.0
subroutine 9 18 50.0
pod 10 10 100.0
total 57 180 31.6


line stmt bran cond sub pod time code
1             package Statistics::Cook;
2 3     3   51667 use Modern::Perl;
  3         6  
  3         24  
3 3     3   2429 use Data::Dumper;
  3         12998  
  3         174  
4 3     3   19 use List::Util qw/sum/;
  3         6  
  3         277  
5 3     3   15 use Carp;
  3         5  
  3         157  
6 3     3   2548 use Moo;
  3         57191  
  3         17  
7 3     3   7645 use Types::Standard qw/Str Num Int ArrayRef/;
  3         207368  
  3         30  
8              
9             our $VERSION = '0.0.6'; # VERSION
10             # ABSTRACT: Statistics::Cook - calculate cook distance of Least squares line fit
11              
12              
13              
14             has x => (
15             is => 'rw',
16             isa => ArrayRef,
17             lazy => 1,
18             default => sub { [] },
19             trigger => 1,
20             );
21              
22              
23             has y => (
24             is => 'rw',
25             isa => ArrayRef,
26             default => sub { [] },
27             lazy => 1,
28             trigger => 1,
29             );
30              
31              
32             has weight => (
33             is => 'rw',
34             isa => ArrayRef
35             );
36              
37              
38             has slope => (
39             is => 'rw',
40             isa => Num,
41             );
42              
43              
44             has intercept=> (
45             is => 'rw',
46             isa => Num,
47             );
48              
49              
50             has regress_done => (
51             is => 'rw',
52             isa => Int,
53             default => 0,
54             lazy => 1,
55             );
56              
57              
58              
59             sub _trigger_x {
60 0     0   0 shift->regress_done(0);
61             }
62              
63             sub _trigger_y {
64 0     0   0 shift->regress_done(0);
65             }
66              
67              
68             sub regress {
69 0     0 1 0 my $self = shift;
70 0         0 my ($x, $y) = ($self->x, $self->y);
71 0 0 0     0 confess "have not got data or x y length is not same" unless(@$x and @$y and @$x == @$y);
      0        
72 0         0 my $sums = $self->computeSums;
73 0         0 my $sqdevx = $sums->{xx} - $sums->{x} ** 2 / scalar(@$x);
74 0 0       0 if ($sqdevx != 0) {
75 0         0 my $sqdevy = $sums->{yy} - $sums->{y} ** 2 / scalar(@$y);
76 0         0 my $sqdevxy = $sums->{xy} - $sums->{x} * $sums->{y} / scalar(@$x);
77 0         0 my $slope = $sqdevxy / $sqdevx;
78 0         0 my $intercept = ($sums->{y} - $slope * $sums->{x}) / @$x;
79 0         0 $self->slope($slope);
80 0         0 $self->intercept( $intercept);
81 0         0 $self->regress_done(1);
82 0         0 return ($intercept, $slope);
83             } else {
84 0         0 confess "Can't fit line when x values are all equal";
85             }
86             }
87              
88              
89             sub computeSums {
90 0     0 1 0 my $self = shift;
91 0         0 my @x = @{$self->x};
  0         0  
92 0         0 my @y = @{$self->y};
  0         0  
93 0         0 my ($sums, @weights);
94 0 0       0 if (defined (my $weight = $self->weight)) {
95 0 0       0 confess "weights does not have same length with x" unless (@$weight == @x);
96 0         0 @weights = @$weight;
97             } else {
98 0         0 @weights = (1) x scalar(@x);
99             }
100 0         0 for my $i (0..$#x) {
101 0         0 my $w = $weights[$i];
102 0         0 $sums->{x} += $w * $x[$i];
103 0         0 $sums->{y} += $w * $y[$i];
104 0         0 $sums->{xx} += $w * $x[$i] ** 2;
105 0         0 $sums->{yy} += $w * $y[$i] ** 2;
106 0         0 $sums->{xy} += $w * $x[$i] * $y[$i];
107             }
108 0         0 return $sums;
109             }
110              
111              
112             sub coefficients {
113 0     0 1 0 my $self = shift;
114 0 0       0 if ($self->regress_done) {
115 0         0 return ($self->intercept, $self->slope);
116             } else {
117 0         0 return $self->regress;
118             }
119             }
120              
121              
122             sub fitted {
123 0     0 1 0 my $self = shift;
124 0 0       0 if ($self->regress_done) {
125 0         0 return map {$self->intercept + $self->slope * $_ } @{$self->x};
  0         0  
  0         0  
126             } else {
127 0         0 my ($a, $b) = $self->regress;
128 0         0 return map {$a + $b * $_} @{$self->x};
  0         0  
  0         0  
129             }
130             }
131              
132              
133             sub residuals {
134 0     0 1 0 my $self = shift;
135 0         0 my @y = @{$self->y};
  0         0  
136 0         0 my @yf = $self->fitted;
137 0         0 return map { $y[$_] - $yf[$_] } 0..$#y;
  0         0  
138             }
139              
140              
141             sub cooks_distance {
142 0     0 1 0 my ($self, @cooks) = shift;
143 0         0 my @yr = $self->residuals;
144 0         0 my @y = @{$self->y};
  0         0  
145 0         0 my @x = @{$self->x};
  0         0  
146 0         0 my $statis = Statistics::Cook->new();
147 0         0 for my $i (0..$#y) {
148 0         0 my @xi = @x;
149 0         0 my @yi = @y;
150 0         0 splice(@xi, $i, 1);
151 0         0 splice(@yi, $i, 1);
152 0         0 $statis->x(\@xi);
153 0         0 $statis->y(\@yi);
154 0         0 my ($a, $b) = $statis->coefficients;
155 0         0 my @yf_new = map {$a + $b * $_ } @x;
  0         0  
156 0         0 my @yf = $self->fitted;
157 0         0 my ($sum1, $sum2) = (0, 0);
158 0         0 for my $j (0..$#yf) {
159 0         0 $sum1 += ($yf[$j] - $yf_new[$j]) ** 2;
160 0         0 $sum2 += $yr[$j] ** 2;
161             }
162 0         0 my $cook = $sum1 * (@y - 2) / $sum2 / 2;
163 0         0 push @cooks, $cook;
164             }
165 0         0 return @cooks;
166             }
167              
168              
169             sub N {
170 0     0 1 0 my ($self, $num, $N) = @_;
171 0   0     0 $N ||= 50;
172 0         0 my @nums = sort { $b <=> $a } @$num;
  0         0  
173 0         0 my $sum = sum(@nums);
174 0         0 my $tmp = 0;
175 0         0 for my $i (0..$#nums) {
176 0         0 $tmp += $nums[$i];
177 0 0       0 return ($nums[$i], $i+1) if ($tmp > $sum * $N / 100);
178             }
179             }
180              
181              
182             sub mean {
183 2     2 1 4 my $self = shift;
184 2 50       7 my @arr = ref $_[0] eq 'ARRAY' ? @{$_[0]} : @_;
  0         0  
185 2         4 my $sum = 0;
186 2         7 $sum += $_ for @arr;
187 2         7 return $sum / @arr;
188             }
189              
190              
191             sub var {
192 2     2 1 4 my $self = shift;
193 2 50       9 my @arr = ref $_[0] eq 'ARRAY' ? @{$_[0]} : @_;
  0         0  
194 2         6 my $m = $self->mean(@arr);
195 2         4 my $sum = 0;
196 2         20 $sum += ($_ - $m) ** 2 for (@arr);
197 2         30 return $sum / $#arr;
198             }
199              
200              
201             sub sd {
202 2     2 1 10370 my $self = shift;
203 2 100       8 my @arr = ref $_[0] eq 'ARRAY' ? @{$_[0]} : @_;
  1         3  
204 2         8 my $var = $self->var(@arr);
205 2         56 return sqrt($var);
206             }
207              
208             1;
209              
210             __END__