1c0b746e5SOllivier Robert;# 2c0b746e5SOllivier Robert;# lr.pl,v 3.1 1993/07/06 01:09:08 jbj Exp 3c0b746e5SOllivier Robert;# 4c0b746e5SOllivier Robert;# 5c0b746e5SOllivier Robert;# Linear Regression Package for perl 6c0b746e5SOllivier Robert;# to be 'required' from perl 7c0b746e5SOllivier Robert;# 8c0b746e5SOllivier Robert;# Copyright (c) 1992 9c0b746e5SOllivier Robert;# Frank Kardel, Rainer Pruy 10c0b746e5SOllivier Robert;# Friedrich-Alexander Universitaet Erlangen-Nuernberg 11c0b746e5SOllivier Robert;# 12224ba2bdSOllivier Robert;# Copyright (c) 1997 by 13224ba2bdSOllivier Robert;# Ulrich Windl <Ulrich.Windl@rz.uni-regensburg.de> 14224ba2bdSOllivier Robert;# (Converted to a PERL 5.004 package) 15c0b746e5SOllivier Robert;# 16c0b746e5SOllivier Robert;############################################################# 17c0b746e5SOllivier Robert 18224ba2bdSOllivier Robertpackage lr; 19224ba2bdSOllivier Robert 20c0b746e5SOllivier Robert## 21c0b746e5SOllivier Robert## y = A + Bx 22c0b746e5SOllivier Robert## 23c0b746e5SOllivier Robert## B = (n * Sum(xy) - Sum(x) * Sum(y)) / (n * Sum(x^2) - Sum(x)^2) 24c0b746e5SOllivier Robert## 25c0b746e5SOllivier Robert## A = (Sum(y) - B * Sum(x)) / n 26c0b746e5SOllivier Robert## 27c0b746e5SOllivier Robert 28c0b746e5SOllivier Robert## 29c0b746e5SOllivier Robert## interface 30c0b746e5SOllivier Robert## 31224ba2bdSOllivier Robert;# init(tag); initialize data set for tag 32224ba2bdSOllivier Robert;# sample(x, y, tag); enter sample 33224ba2bdSOllivier Robert;# Y(x, tag); compute y for given x 34224ba2bdSOllivier Robert;# X(y, tag); compute x for given y 35224ba2bdSOllivier Robert;# r(tag); regression coefficient 36224ba2bdSOllivier Robert;# cov(tag); covariance 37224ba2bdSOllivier Robert;# A(tag); 38224ba2bdSOllivier Robert;# B(tag); 39224ba2bdSOllivier Robert;# sigma(tag); standard deviation 40224ba2bdSOllivier Robert;# mean(tag); 41c0b746e5SOllivier Robert######################### 42c0b746e5SOllivier Robert 43224ba2bdSOllivier Robertsub init 44224ba2bdSOllivier Robert{ 45224ba2bdSOllivier Robert my $self = shift; 46c0b746e5SOllivier Robert 47224ba2bdSOllivier Robert $self->{n} = 0; 48224ba2bdSOllivier Robert $self->{sx} = 0.0; 49224ba2bdSOllivier Robert $self->{sx2} = 0.0; 50224ba2bdSOllivier Robert $self->{sxy} = 0.0; 51224ba2bdSOllivier Robert $self->{sy} = 0.0; 52224ba2bdSOllivier Robert $self->{sy2} = 0.0; 53c0b746e5SOllivier Robert} 54c0b746e5SOllivier Robert 55ea906c41SOllivier Robertsub sample($$) 56c0b746e5SOllivier Robert{ 57224ba2bdSOllivier Robert my $self = shift; 58224ba2bdSOllivier Robert my($_x, $_y) = @_; 59c0b746e5SOllivier Robert 60224ba2bdSOllivier Robert ++($self->{n}); 61224ba2bdSOllivier Robert $self->{sx} += $_x; 62224ba2bdSOllivier Robert $self->{sy} += $_y; 63224ba2bdSOllivier Robert $self->{sxy} += $_x * $_y; 64224ba2bdSOllivier Robert $self->{sx2} += $_x**2; 65224ba2bdSOllivier Robert $self->{sy2} += $_y**2; 66c0b746e5SOllivier Robert} 67c0b746e5SOllivier Robert 68ea906c41SOllivier Robertsub B() 69c0b746e5SOllivier Robert{ 70224ba2bdSOllivier Robert my $self = shift; 71c0b746e5SOllivier Robert 72224ba2bdSOllivier Robert return 1 unless ($self->{n} * $self->{sx2} - $self->{sx}**2); 73224ba2bdSOllivier Robert return ($self->{n} * $self->{sxy} - $self->{sx} * $self->{sy}) 74224ba2bdSOllivier Robert / ($self->{n} * $self->{sx2} - $self->{sx}**2); 75c0b746e5SOllivier Robert} 76c0b746e5SOllivier Robert 77ea906c41SOllivier Robertsub A() 78c0b746e5SOllivier Robert{ 79224ba2bdSOllivier Robert my $self = shift; 80c0b746e5SOllivier Robert 81ea906c41SOllivier Robert return ($self->{sy} - B() * $self->{sx}) / $self->{n}; 82c0b746e5SOllivier Robert} 83c0b746e5SOllivier Robert 84ea906c41SOllivier Robertsub Y() 85c0b746e5SOllivier Robert{ 86224ba2bdSOllivier Robert my $self = shift; 87c0b746e5SOllivier Robert 88ea906c41SOllivier Robert return A() + B() * $_[$[]; 89c0b746e5SOllivier Robert} 90c0b746e5SOllivier Robert 91ea906c41SOllivier Robertsub X() 92c0b746e5SOllivier Robert{ 93224ba2bdSOllivier Robert my $self = shift; 94c0b746e5SOllivier Robert 95ea906c41SOllivier Robert return ($_[$[] - A()) / B(); 96c0b746e5SOllivier Robert} 97c0b746e5SOllivier Robert 98ea906c41SOllivier Robertsub r() 99c0b746e5SOllivier Robert{ 100224ba2bdSOllivier Robert my $self = shift; 101c0b746e5SOllivier Robert 102224ba2bdSOllivier Robert my $s = ($self->{n} * $self->{sx2} - $self->{sx}**2) 103224ba2bdSOllivier Robert * ($self->{n} * $self->{sy2} - $self->{sy}**2); 104c0b746e5SOllivier Robert 105c0b746e5SOllivier Robert return 1 unless $s; 106c0b746e5SOllivier Robert 107224ba2bdSOllivier Robert return ($self->{n} * $self->{sxy} - $self->{sx} * $self->{sy}) / sqrt($s); 108c0b746e5SOllivier Robert} 109c0b746e5SOllivier Robert 110ea906c41SOllivier Robertsub cov() 111c0b746e5SOllivier Robert{ 112224ba2bdSOllivier Robert my $self = shift; 113c0b746e5SOllivier Robert 114224ba2bdSOllivier Robert return ($self->{sxy} - $self->{sx} * $self->{sy} / $self->{n}) 115224ba2bdSOllivier Robert / ($self->{n} - 1); 116c0b746e5SOllivier Robert} 117c0b746e5SOllivier Robert 118ea906c41SOllivier Robertsub sigma() 119c0b746e5SOllivier Robert{ 120224ba2bdSOllivier Robert my $self = shift; 121c0b746e5SOllivier Robert 122224ba2bdSOllivier Robert return 0 if $self->{n} <= 1; 123224ba2bdSOllivier Robert return sqrt(($self->{sy2} - ($self->{sy} * $self->{sy}) / $self->{n}) 124224ba2bdSOllivier Robert / ($self->{n})); 125c0b746e5SOllivier Robert} 126c0b746e5SOllivier Robert 127ea906c41SOllivier Robertsub mean() 128c0b746e5SOllivier Robert{ 129224ba2bdSOllivier Robert my $self = shift; 130c0b746e5SOllivier Robert 131224ba2bdSOllivier Robert return 0 if $self->{n} <= 0; 132224ba2bdSOllivier Robert return $self->{sy} / $self->{n}; 133c0b746e5SOllivier Robert} 134c0b746e5SOllivier Robert 135224ba2bdSOllivier Robertsub new 136224ba2bdSOllivier Robert{ 137224ba2bdSOllivier Robert my $class = shift; 138224ba2bdSOllivier Robert my $self = { 139224ba2bdSOllivier Robert (n => undef, 140224ba2bdSOllivier Robert sx => undef, 141224ba2bdSOllivier Robert sx2 => undef, 142224ba2bdSOllivier Robert sxy => undef, 143224ba2bdSOllivier Robert sy => undef, 144224ba2bdSOllivier Robert sy2 => undef) 145224ba2bdSOllivier Robert }; 146224ba2bdSOllivier Robert bless $self, $class; 147224ba2bdSOllivier Robert init($self); 148224ba2bdSOllivier Robert return $self; 149224ba2bdSOllivier Robert} 150c0b746e5SOllivier Robert 151c0b746e5SOllivier Robert1; 152