/*Copyright (C) 1992, 1995 by Thomas Glen Smith.  All Rights Reserved.*/
/* regress APL2 V1.0.0 *************************************************
* Called from matinv and mdivide when the right argument is non-square.*
* Simple linear regression will be performed with the left argument    *
* containing values for the dependent variable y, and the right        *
* argument containing values for the independent variable x.  Each     *
* column of the right argument contains the value of x to some         *
* exponent, e.g. the first column may contain x to the 0th power,      *
* the second column x to the 1st power, etc.  The result will be a     *
* vector of coefficients providing the best fit in the least squares   *
* sense.  Given -                                                      *
* X = matrix pointed to by rite.                                       *
* Y = matrix pointed to by left. if left is NIL, use identity matrix.  *
* BETA = desired result.                                               *
* In matrix notation, beta is produced as follows:                     *
*   XPX = X'*X cross products matrix.                                  *
*   XPXI = INV(XPX) inverse cross products matrix.                     *
*   XPY = X'*Y sum(of sy) and sum(of xy).                              *
*   BETA = XPXI * XPY.                                                 *
***********************************************************************/
#define INCLUDES APLCB+APLMEM
#include "includes.h"
double *regress(left,rite,lrows,lcols,rrows,rcols)
Aplcb left,rite;
int lrows,lcols,rrows,rcols;
{
	Getcb; Regrest;
	extern int aplerr;
	Aplcb xpx,xpy;
	double *beta;

	beta = malloc(rcols*lcols*sizeof(double));
	if (aplerr) return(NULL); /* out of memory */
	xpx = getcb(NULL,rcols*rcols,APLNUMB+APLTEMP,2,NULL);
	*(xpx->apldim) = rcols;
	*(xpx->apldim+1) = rcols;
	xpy = getcb(NULL,rcols*lcols,APLNUMB+APLTEMP,2,NULL);
	*(xpy->apldim) = rcols;
	*(xpy->apldim+1) = lcols;
	return(regrest(left,rite,lrows,lcols,rrows,rcols,beta,xpx,xpy));
}