#======================================================================
#                    P E A R S N . P L 
#                    doc: Wed Mar 24 11:23:29 1999
#                    dlm: Mon Jul 24 15:02:14 2006
#                    (c) 1999 A.M. Thurnherr
#                    uE-Info: 26 0 NIL 0 0 72 66 2 4 NIL ofnI
#======================================================================

# HISTORY:
#	Mar 24, 1999: - created from NR c-version
#	Mar 26, 1999: - cosmetic changes
#	May 23, 1999: - allowed for N==2
#	Oct 04, 1999: - changed from specfuns to funs
#	Oct 13, 1999: - BUG: had to change TINY from 1e-20 to 1e-19
#	Nov 11, 1999: - BUG: had to change TINY from to 1e-16
#	Dec 11, 2001: - BUG: NaNs had not been handled correctly
#	Dec 12, 2001: - BUG: croak() had been used (this produces
#						 a pipe #ERROR# output which is wrong
#						 if pearsn is called within eval as in [fit])
#	Jan  9, 2006: - removed @antsFlagged

require "$ANTS/libfuns.pl";

{											# static scope

my($TINY) = 1e-16;							# for complete correlation

# get correlation coefficient (retval); N (ref); significance level at which 
# null hypothesis of zero correlation is disproved (ref; small value
# indicates significant correlation); and Fisher's z (ref). Missing refs
# indicate values are not returned. Adapted for ANTS.

sub pearsn(@)
{
	my($xfnr,$yfnr,$NR,$pR,$zR) = @_;
	my($n,$r);
	my($j);
	my($yt,$xt,$t,$df);
	my($syy,$sxy,$sxx,$ay,$ax);

	for ($j=0; $j<=$#ants_; $j++) {
		next unless (numberp($ants_[$j][$xfnr]) && numberp($ants_[$j][$yfnr]));
		$n++;
		$ax += $ants_[$j][$xfnr];
		$ay += $ants_[$j][$yfnr];
	}
	die("$0 (pearsn.pl): no data\n") unless ($n >= 2);
	$ax /= $n;
	$ay /= $n;
	for ($j=0; $j<=$#ants_; $j++) {
		next unless (numberp($ants_[$j][$xfnr]) && numberp($ants_[$j][$yfnr]));
		$xt = $ants_[$j][$xfnr] - $ax;
		$yt = $ants_[$j][$yfnr] - $ay;
		$sxx += $xt * $xt;
		$syy += $yt * $yt;
		$sxy += $xt * $yt;
	}
	$r   = $sxy/(sqrt($sxx * $syy) + $TINY);
	$df  = $n - 2;
	$t   = $r * sqrt($df/((1-$r+$TINY) * (1+$r+$TINY)));
	$$NR = $n if (defined($NR));
	$$pR = ($n == 2) ? 0 : &betai(0.5*$df,0.5,$df/($df+$t*$t))
		if (defined($pR));
	$$zR = 0.5 * log((1+$r+$TINY) / (1-$r+$TINY))
		if (defined($zR));
	return $r;
}

}													# end of static scope