MultipleLinearRegression.cpp

#include "StdAfx.h"
#include "MultipleLinearRegression.h"

CMultipleLinearRegression::CMultipleLinearRegression(LR_CountT nRegressors)
{
	m_pSumXn = NULL;
	ReInit(nRegressors);
}
void CMultipleLinearRegression::ReInit(LR_CountT nRegressors)
{
	LR_CountT nSize;
	LR_CountT nSumXnCnt;
	LR_CountT nSumXnYCnt;
	LR_CountT nSumXnXmCnt;
	LR_CountT nResultsCnt;
	if( nRegressors < 1 )
	{
		nRegressors = 1;
	}
	nSumXnCnt = nRegressors + 1;
	nSumXnYCnt = nSumXnCnt;
	if( (nRegressors&1) == 0 )
	{
		nSumXnXmCnt = (nRegressors+1) * (nRegressors/2);
	}
	else
	{
		nSumXnXmCnt = (nRegressors+2) * (nRegressors/2) + 1; // ((n+1)*(n\2)) + ((n\2) + 1)
	}
	nResultsCnt = nSumXnCnt;
	nSize = nSumXnCnt + nSumXnYCnt + nSumXnXmCnt + nResultsCnt;
	do
	{
		if( m_pSumXn != NULL )
		{
			if( m_nRegressors == nRegressors )
			{
				break;
			}
			delete []m_pSumXn;
		}
		m_pSumXn = new LR_DataT[nSize];
		m_pSumXnY = &m_pSumXn[nSumXnCnt];
		m_pSumXnXm = &m_pSumXnY[nSumXnYCnt];
		m_pResults = &m_pSumXnXm[nSumXnXmCnt];
	}while(0);
	m_nRegressors = nRegressors;
	
	// for( LR_CountT i = 0 ; i < nSize ; i++ )
	// {
		// m_pSumXn[i] = 0.0;
	// }
	memset(m_pSumXn , 0 , nSize*sizeof(LR_DataT));

	m_uRefreshType = 0;
}

CMultipleLinearRegression::~CMultipleLinearRegression()
{
	delete []m_pSumXn;
	m_pSumXn = NULL;
}

BOOL CMultipleLinearRegression::AddXY(LR_DataT* pX , LR_DataT y)
{
	LR_CountT i , j , k , n;
	LR_DataT x;
	m_pSumXn[0] += 1.0;
	m_pSumXnY[0] += y;
	for( i = 0 , k = 0 , n = m_nRegressors ; i < n ; i++ )
	{
		x = pX[i];
		m_pSumXn[i+1] += x;
		m_pSumXnY[i+1] += x * y;
		m_pSumXnXm[k++] += x * x;
		for( j = i + 1 ; j < n ; j++ )
		{
			m_pSumXnXm[k++] += x * pX[j];
		}
	}
	m_uRefreshType = 0;
	return TRUE;
}
void CMultipleLinearRegression::SetRefresh()
{
	m_uRefreshType = 0;
}

LR_DataT* CMultipleLinearRegression::GetResults()
{
	if( m_uRefreshType != 1 )
	{
		if( !Calculate() )
		{
			return NULL;
		}
	}
	return m_pResults;
}

BOOL CMultipleLinearRegression::IsValid()
{
	return (m_pSumXn[0] >= (LR_DataT)(m_nRegressors + 1));
}

BOOL CMultipleLinearRegression::Calculate()
{
	LR_CountT i , j , k;
	LR_CountT nRows;
	LR_CountT nColumns;
	LR_CountT iMaxIdx;
	LR_DataT* pMatrix;
	LR_DataT** ppIndex;
	LR_DataT* pMaxRow;
	LR_DataT Val;
	
	nRows = m_nRegressors + 1;
	nColumns = m_nRegressors + 2;
	
	pMatrix = new LR_DataT[nColumns * nRows];
	ppIndex = new LR_DataT*[nRows];
	for( i = 0 ; i < nRows ; i++ )
	{
		ppIndex[i] = &pMatrix[i*nColumns];
	}
	m_uRefreshType = 0;

	/*
	E:Sigma
	num    Ex1      Ex2      Ex3      Ex4      Ey
	Ex1    Ex1*x1   Ex1*x2   Ex1*x3   Ex1*x4   Ex1*y
	Ex2    Ex2*x1   Ex2*x2   Ex2*x3   Ex2*x4   Ex2*y
	Ex3    Ex3*x1   Ex3*x2   Ex3*x3   Ex3*x4   Ex3*y
	Ex4    Ex4*x1   Ex4*x2   Ex4*x3   Ex4*x4   Ex4*y
	*/
	// Filling polynomial matrix
	pMatrix[0] = m_pSumXn[0];
	for( i = 1 ; i < nRows ; i++ )
	{
		pMatrix[i*nColumns] = pMatrix[i] = m_pSumXn[i];
	}
	for( i = 0 ; i < nRows ; i++ )
	{
		pMatrix[(i+1)*nColumns-1] = m_pSumXnY[i];
	}
	for( i = 1 , k = 0 ; i < nRows ; i++ )
	{
		pMatrix[(i*nColumns) + i] = m_pSumXnXm[k++];
		for( j = i + 1 ; j < nRows ; j++ )
		{
			pMatrix[(i*nColumns) + j] = pMatrix[(j*nColumns) + i] = m_pSumXnXm[k++];
		}
	}
	
	/*
	// Output raw data
	CString str;
	for( i = 0 , k = (nColumns * nRows) ; i < k ; i++ )
	{
		str.Format(_T("%g  ") , (double)pMatrix[i]);
		TRACE(str);
	}
	*/

	// Gauss elimination
	for( i = 0 ; i < nRows ; i++ )
	{
		Val = LR_ABS(ppIndex[i][i]);
		iMaxIdx = i;
		for( j = i+1 ; j < nRows ; j++ )
		{
			if( Val < LR_ABS(ppIndex[j][i]) )
			{
				Val = LR_ABS(ppIndex[j][i]);
				iMaxIdx = j;
			}
		}
		if( iMaxIdx != i )
		{
			pMaxRow = ppIndex[i];
			ppIndex[i] = ppIndex[iMaxIdx];
			ppIndex[iMaxIdx] = pMaxRow;
		}
		if( Val <= LR_MinValue )
		{
			goto PLR_Calculate_ErrorEnd;
		}
		for( j = i+1 ; j < nRows ; j++ )
		{
			Val = ppIndex[j][i] / ppIndex[i][i];
			for( k = i+1 ; k <= nRows ; k++ )
			{
				ppIndex[j][k] -= ppIndex[i][k] * Val;
			}
			ppIndex[j][i] = 0.0;
		}
	}

	// Polynomial solving
	for( i = nRows ; i != 0 ; )
	{
		i--;
		m_pResults[i] = ppIndex[i][nRows] / ppIndex[i][i]; // nRows : nColumns - 1 (Y)
		if( LR_ABS(ppIndex[i][i]) <= LR_MinValue ) // Save extreme values : m_pResults[i]
		{
			goto PLR_Calculate_ErrorEnd;
		}
		for( j = i ; j != 0 ; )
		{
			j--;
			ppIndex[j][nRows] -= ppIndex[j][i] * m_pResults[i]; // nRows : nColumns - 1 (Y)
			ppIndex[j][i] = 0.0;
		}
	}
	m_uRefreshType = 1;
PLR_Calculate_ErrorEnd:
	delete []pMatrix;
	delete ppIndex;
	return m_uRefreshType != 0;
}

LR_DataT* CMultipleLinearRegression::GetResultsInterceptZero()
{
	if( m_uRefreshType != 2 )
	{
		if( !CalculateInterceptZero() )
		{
			return NULL;
		}
	}
	return m_pResults;
}
BOOL CMultipleLinearRegression::CalculateInterceptZero()
{
	LR_CountT i , j , k;
	LR_CountT nRows;
	LR_CountT nColumns;
	LR_CountT iMaxIdx;
	LR_DataT* pMatrix;
	LR_DataT** ppIndex;
	LR_DataT* pMaxRow;
	LR_DataT Val;
	
	nRows = m_nRegressors;
	nColumns = m_nRegressors + 1;
	
	pMatrix = new LR_DataT[nColumns * nRows];
	ppIndex = new LR_DataT*[nRows];
	for( i = 0 ; i < nRows ; i++ )
	{
		ppIndex[i] = &pMatrix[i*nColumns];
	}
	m_uRefreshType = 0;

	/*
	Ex1*x1   Ex1*x2   Ex1*x3   Ex1*x4   Ex1*y
	Ex2*x1   Ex2*x2   Ex2*x3   Ex2*x4   Ex2*y
	Ex3*x1   Ex3*x2   Ex3*x3   Ex3*x4   Ex3*y
	Ex4*x1   Ex4*x2   Ex4*x3   Ex4*x4   Ex4*y
	*/
	// Filling polynomial matrix
	for( i = 0 ; i < nRows ; i++ )
	{
		pMatrix[(i+1)*nColumns-1] = m_pSumXnY[i+1];
	}
	for( i = 0 , k = 0 ; i < nRows ; i++ )
	{
		pMatrix[(i*nColumns) + i] = m_pSumXnXm[k++];
		for( j = i + 1 ; j < nRows ; j++ )
		{
			pMatrix[(i*nColumns) + j] = pMatrix[(j*nColumns) + i] = m_pSumXnXm[k++];
		}
	}

	// Gauss elimination
	for( i = 0 ; i < nRows ; i++ )
	{
		Val = LR_ABS(ppIndex[i][i]);
		iMaxIdx = i;
		for( j = i+1 ; j < nRows ; j++ )
		{
			if( Val < LR_ABS(ppIndex[j][i]) )
			{
				Val = LR_ABS(ppIndex[j][i]);
				iMaxIdx = j;
			}
		}
		if( iMaxIdx != i )
		{
			pMaxRow = ppIndex[i];
			ppIndex[i] = ppIndex[iMaxIdx];
			ppIndex[iMaxIdx] = pMaxRow;
		}
		if( Val <= LR_MinValue )
		{
			goto PLR_Calculate_ErrorEnd;
		}
		for( j = i+1 ; j < nRows ; j++ )
		{
			Val = ppIndex[j][i] / ppIndex[i][i];
			for( k = i+1 ; k <= nRows ; k++ )
			{
				ppIndex[j][k] -= ppIndex[i][k] * Val;
			}
			ppIndex[j][i] = 0.0;
		}
	}

	// Polynomial solving
	for( i = nRows ; i != 0 ; )
	{
		i--;
		m_pResults[i] = ppIndex[i][nRows] / ppIndex[i][i]; // nRows : nColumns - 1 (Y)
		if( LR_ABS(ppIndex[i][i]) <= LR_MinValue ) // Save extreme values : m_pResults[i]
		{
			goto PLR_Calculate_ErrorEnd;
		}
		for( j = i ; j != 0 ; )
		{
			j--;
			ppIndex[j][nRows] -= ppIndex[j][i] * m_pResults[i]; // nRows : nColumns - 1 (Y)
			ppIndex[j][i] = 0.0;
		}
	}
	m_uRefreshType = 2;
PLR_Calculate_ErrorEnd:
	delete []pMatrix;
	delete ppIndex;
	return m_uRefreshType != 0;
}