/*******************************************************************************************
Linear and Multiple Regression 

Author: Juan Pablo Romero Bernal - Grupo Linux Universidad Distrital
e-mail: jromerobernal AT gmail DOT com
Date: May 2006
License: GPL


Purpose: Estimate the regression (lineal and multiple) coefficients and others
related values (as variance, true intervals, etc.) using matrix theory.

	
********************************************************************************************/

fpprec:10$

/**
	Coefficients Linear Regression Estimation (CRE). Take a matrix as argument, being 
	this matrix 2xN. Return the a and b values.  
	
*/


crle(D):=block
	([n:length(D),
	 M:copymatrix(D)],	
	t1:sum((M[i][1])*(M[i][2]),i,1,n),
	t2:sum(M[i][1],i,1,n),
	t3:sum(M[i][2],i,1,n),
	t4:sum((M[i][1])^2,i,1,n),
	b:(n*t1-(t2)*(t3))/(n*t4-(t2)^2),
	a:(t3-b*t2)/(n),
	display(b),
	display(a),
	cx:makelist(M[i][1],i,1,n),
	cy:makelist(M[i][2],i,1,n),
	plot2d(a+b*x,[x,M[1][1],M[n][2]],[gnuplot_preamble,"set title 'Regression Line'"]),
	plot2d([discrete,cx,cy],[x,M[1][1],M[n][2]],[gnuplot_curve_styles,["with points"]],[gnuplot_preamble,"set title 'Dispersion Diagram'"]),
	plot2d([a+b*x,[discrete,cx,cy]],[x,M[1][1],M[n][2]],[gnuplot_preamble,"set title 'Regression Line vs Dispersion Diagram'"],[gnuplot_curve_styles,["with points"]])
)$

/**
	Only Dispersion Diagram

*/

dd(D):=block
	([n:length(D),
	 M:copymatrix(D)],
	 cx:makelist(M[i][1],i,1,n),
	 cy:makelist(M[i][2],i,1,n),
	 plot2d([discrete,cx,cy],[gnuplot_curve_styles,["with points"]],[gnuplot_preamble,"set title 'Dispersion Diagram'"])
)$


/**

	Coefficients Multiple Regression Estimation. Take a matrix as argument. The 
	matrix is NxN.  Return the vector with regression coefficients and variance estimation. 	

*/

crme(A):=block
	([n:length(A)],
	M:copymatrix(A),
	k:length(transpose(M)),
	for i:1 step 1 thru k do
		(
		if(i = 1) then
			a[i]:sum(M[t][1],t,1,n)
		else			
			a[i]:sum(M[t][1]*M[t][i],t,1,n)		
		),
	y1:matrix(makelist(a[i],i,1,k)),
	g:transpose(y1),
	X:ident(k),
	for i:1 step 1 thru k do (
		for j:i step 1 thru k do 
			(
			if(j = 1 and  i = 1) then
				X[i][j]:n
			elseif(i = 1) then (
				X[i][j]:sum(M[t][j],t,1,n),
				X[j][i]:sum(M[t][j],t,1,n))
			else (
				X[i][j]:sum(M[t][j]*M[t][i],t,1,n),
				X[j][i]:sum(M[t][j]*M[t][i],t,1,n))
			)
	),
	X1:X^^-1,
	b:X1.g,
	display(b),
	k1:sum((M[t][1])^2,t,1,n),
	k2:(sum(M[t][1],t,1,n))^2,
	k3:sum(b[f]*g[f],f,1,k),
	SST:k1-((k2)/(n)),
	SSR:k3-((k2)/(n)),
	est:(SST-SSR)/(n-k),
	display(est)
)$