#!/bin/sh 

# $Id: linear-reg.awk,v 1.1 2005/04/21 16:02:14 rajo Exp $

# Author: Lubomir Host 'rajo' <host8 AT kepler.fmph.uniba.sk>
# Description: Linear regresion with awk

# Input data format:
# <x_i> <y_i> [delta_x] [delta_y]

# Run this with awk:
exec awk '
function linreg_uy()
{
	# a,b must be known!
	a = linreg_a();
	b = linreg_b();

	if (N < 3) {
		print "Warning: N <= 2 !!  u(y) = ??" > "/dev/stderr"
		return 0;
	}
	uy = (sumy2 - 2*a*sumxy + 2*b*(a*sumx - sumy) + a*a*sumx2 + b*b) / (N - 2);
	return uy;
}

function linreg_a()
{
	a = (N * sumxy - sumx * sumy) / (N * sumx2 - sumx * sumx);
	return a;
}

function linreg_b()
{
	b = (sumy - a * sumx) / N
	return b;
}

function linreg_da()
{
	da = N * uy * uy / (N * sumx2 - sumx * sumx)
	return da;
}

function linreg_db()
{
	db = sumx2 * uy * uy / (N * sumx2 - sumx * sumx)
	return db;
}

function linreg_r()
{
	#r = (-1 * sumx) / sqrt(N * sumx2)
	r = (N*sumxy - sumx*sumy) / sqrt((N*sumx2 - sumx*sumx) * (N*sumy2 - sumy*sumy))
	return r;
}

function quadrreg_A()
{
	return A;
}

function quadrreg_B()
{
	return B;
}

function quadrreg_C(){
	return C;
}

BEGIN{
	N  = 0;
	x  = 0;  y = 0;
	dx = 0; dy = 0;
	sumx   = 0; sumy   = 0;
	sumdx  = 0; sumdy  = 0;
	sumxy  = 0;
	sumx2  = 0; sumy2  = 0;
	sumxdx = 0; sumdxy = 0;
	sumxdy = 0;
	a  = 0;  b = 0;
	da = 0; db = 0;
	uy = 0;  r = 0;
	A  = 0;  B = 0; C = 0;
}

// {
	if (NF < 4)
		dy = 0;
	else
		dy = $4;
	if (NF < 3)
		dx = 0;
	else
		dx = $4;
	if (NF < 2) {
		print "Data format error: Ignoring line " NR > "/dev/stderr"
		next;
	}
	N += 1;
	x = $1; y = $2;
	sumx += x; sumy += y;
	sumdx += dx; sumdy += dy;
	sumxy += x * y;
	sumx2y += x * x * y;
	sumx2 += x * x;
	sumy2 += y * y;
	sumx3 += x * x * x;
	#sumy3 += y * y * y;
	sumx4 += x * x * x * x;
	#sumy4 += y * y * y * y;
	sumxdx = x * dx;
	sumdxy = y * dx;
	sumxdy = x * dy;
}

END {
	print "# ----------------------"
	print "# N         = " N;
	print "# Sum(x)    = " sumx;
	print "# Sum(y)    = " sumy;
	print "# Sum(dx)   = " sumdx;
	print "# Sum(dy)   = " sumdy;
	print "# Sum(x*y)  = " sumxy;
	print "# Sum(x*x)  = " sumx2;
	print "# Sum(y*y)  = " sumy2;
	print "# Sum(x*dx) = " sumxdx;
	print "# Sum(x*dy) = " sumxdy;
	print "# Sum(dx*y) = " sumdxy;
	print "# <x>       = " sumx / N;
	print "# <y>       = " sumy / N;
	print "# u(y)      = " linreg_uy();
	print "# y         = a*x + b"
	print "# a         = " a;
	print "# b         = " b;
	print "# da        = " linreg_da();
	print "# db        = " linreg_db();
	print "# r         = " linreg_r();

	print "# Sum(x*x*y)= " sumx2y;
	print "# Sum(x*x*x)= " sumx3;
	print "# Sum(x^4)  = " sumx4;
	print "# Y         = A*x^2 + B*x + C"
	det  = N    * (sumx2*sumx4  - sumx3*sumx3 ) - sumx * (sumx *sumx4  - sumx3*sumx2 ) + sumx2 * (sumx *sumx3  - sumx2*sumx2 )
	det1 = sumy * (sumx2*sumx4  - sumx3*sumx3 ) - sumx * (sumxy*sumx4  - sumx3*sumx2y) + sumx2 * (sumxy*sumx3  - sumx2*sumx2y)
	det2 = N    * (sumxy*sumx4  - sumx3*sumx2y) - sumy * (sumx *sumx4  - sumx3*sumx2 ) + sumx2 * (sumx *sumx2y - sumxy*sumx2 )
	det3 = N    * (sumx2*sumx2y - sumxy*sumx3 ) - sumx * (sumx *sumx2y - sumxy*sumx2 ) + sumy  * (sumx *sumx3  - sumx2*sumx2 )
	print "# det()     = " det ;
	print "# det1      = " det1 ;
	print "# det2      = " det2 ;
	print "# det3      = " det3 ;
	print "# A         = " det3 / det ; #quadrreg_A();
	print "# B         = " det2 / det ; #quadrreg_B();
	print "# C         = " det1 / det ; #quadrreg_C();
	
}

' $* # end of awk script