lsq.cpp

orange源码 数据挖掘技术

			pos1 = start;
			pos2 = pos;
			total = zero;
			for( k = row+1; k <= col-1; ++k) {
				pos2 = pos2 + nreq - k;
				total = total - r[pos1] * rinv[pos2];
				pos1 = pos1 + 1;
			}
			rinv[pos] = total - r[pos1];
			pos = pos - 1;
		}
	}
}


void lsq::partial_corr(int in, double *cormat, int dimc, double *ycorr, int& ifault) {
/*
!     Replaces subroutines PCORR and COR of:
!     ALGORITHM AS274  APPL. STATIST. (1992) VOL.41, NO. 2

  !     Calculate partial correlations after the variables in rows
  !     1, 2, ..., IN have been forced into the regression.
  !     If IN = 1, and the first row of R represents a constant in the
  !     model, then the usual simple correlations are returned.
  !	    If IN = 0, the value returned in array CORMAT for the correlation
  !     of variables Xi & Xj is:
  !       sum ( Xi.Xj ) / sqrt ( sum (Xi^2) . sum (Xj^2) )
  !     On return, array CORMAT contains the upper triangle of the matrix of
  !     partial correlations stored by rows, excluding the 1's on the diagonal.
  !     e.g. if IN = 2, the consecutive elements returned are:
  !     (3,4) (3,5) ... (3,ncol), (4,5) (4,6) ... (4,ncol), etc.
  !     Array YCORR stores the partial correlations with the Y-variable
  !     starting with YCORR(IN+1) = partial correlation with the variable in
  !     position (IN+1).
  !
  !--------------------------------------------------------------------------
	*/
	
	
	int base_pos, pos, row, col, col1, col2, pos1, pos2,k;
	double *rms, sumxx, sumxy, sumyy, *work;
	
	rms = new double[ncol+1];
	work = new double[ncol+1];
	
	
	ifault = 0;
	if (in < 0 || in > ncol-1) ifault = ifault + 4;
	if (dimc < (ncol-in)*(ncol-in-1)/2) ifault = ifault + 8;
	if (ifault != 0) return;
	
	//!     Base position for calculating positions of elements in row (IN+1) of R.
	
	base_pos = in*ncol - (in+1)*(in+2)/2;
	
	
	//!     Calculate 1/RMS of elements in columns from IN to (ncol-1).
	
	if (d[in+1] > zero) 
		rms[in+1] = one / sqrt(d[in+1]);
	for (col = in+2;col<= ncol; ++col) {
		pos = base_pos + col;
		sumxx = d[col];
		for(row = in+1; row <= col-1; ++row) {
			sumxx = sumxx + d[row] * r[pos]*r[pos];
			pos = pos + ncol - row - 1;
		}
		if (sumxx > zero) 
			rms[col] = one / sqrt(sumxx);
		else {
			rms[col] = zero;
			ifault = -col;
			
		}
	}
	
	//    Calculate 1/RMS for the Y-variable
	
	sumyy = sserr;
	for(row = in+1;row <= ncol;++row) {
		sumyy = sumyy + d[row]* rhs[row]*rhs[row];
	}
	if (sumyy > zero) sumyy = one / sqrt(sumyy);
	
	/*
	!     Calculate sums of cross-products.
	!     These are obtained by taking dot products of pairs of columns of R,
	!     but with the product for each row multiplied by the row multiplier
	!     in array D.
	*/
	pos = 1;
	for (col1 = in+1; col1 <= ncol; ++col1) {
		sumxy = zero;
		for (k = col1+1; k <= ncol; ++k) {
			work[k] = zero;
		}
		pos1 = base_pos + col1;
		for(row = in+1; row <= col1-1; ++row) {
			pos2 = pos1 + 1;
			for(col2 = col1+1; col2 <= ncol; ++col2) {
				work[col2] = work[col2] + d[row] * r[pos1] * r[pos2];	  
				pos2 = pos2 + 1;
			}
			sumxy = sumxy + d[row] * r[pos1] * rhs[row];
			pos1 = pos1 + ncol - row - 1;
		}
		
		//!     Row COL1 has an implicit 1 as its first element (in column COL1)
		
		pos2 = pos1 + 1;
		for (col2 = col1+1; col2<= ncol; ++col2) {
			work[col2] = work[col2] + d[col1] * r[pos2];
			pos2 = pos2 + 1;
			cormat[pos] = work[col2] * rms[col1] * rms[col2];
			pos = pos + 1;
		}
		sumxy = sumxy + d[col1] * rhs[col1];
		ycorr[col1] = sumxy * rms[col1] * sumyy;
	}
	
	for (k = 1; k <= in; ++k) {
		ycorr[k] = zero;
	}
	
	return;
}




void lsq::vmove(int from, int to, int& ifault) {
/*
!     ALGORITHM AS274 APPL. STATIST. (1992) VOL.41, NO. 2
!     Move variable from position FROM to position TO in an
!     orthogonal reduction produced by AS75.1.
!
!--------------------------------------------------------------------------
	*/
	
	double d1, d2, x, d1new, d2new, cbar, sbar, y;
	int m, first, last, inc, m1, m2, mp1, col, pos, row,k;
	
	//!     Check input parameters
	
	ifault = 0;
	if (from < 1 || from > ncol) ifault = ifault + 4;
	if (to < 1 || to > ncol) ifault = ifault + 8;
	if (ifault != 0) return;
	
	if (from == to) {
		return;
	}
	
	if (!rss_set) ss();
	
	if (from < to) {
		first = from;
		last = to;
		inc = 1;
	} else {
		first = from - 1;
		last = to-1;
		inc = -1;
	}
	
	for (m = first; m != last; m += inc) { 
		
		//!     Find addresses of first elements of R in rows M and (M+1).
		
		m1 = (m-1)*(ncol+ncol-m)/2 + 1;
		m2 = m1 + ncol - m;
		mp1 = m + 1;
		d1 = d[m];
		d2 = d[mp1];
		
		//!     Special cases.
		
		if ((d1 < vsmall && d2 < vsmall)) {
			goto l40;
		}
		x = r[m1];
		if (fabs(x) * sqrt(d1) < tol[mp1]) {
			x = zero;
		}
		if (d1 < vsmall || fabs(x) < vsmall) {
			d[m] = d2;
			d[mp1] = d1;
			r[m1] = zero;
			for(col = m+2;col<= ncol; ++col) {
				m1 = m1 + 1;
				x = r[m1];
				r[m1] = r[m2];
				r[m2] = x;
				m2 = m2 + 1;
			}
			x = rhs[m];
			rhs[m] = rhs[mp1];
			rhs[mp1] = x;
		} else if (d2 < vsmall) {
			d[m] = d1 * x*x;
			r[m1] = one / x;
			for (k = m1+1; k <=m1+ncol-m-1; ++k) {
				r[k] /= x;
			}
			rhs[m] = rhs[m] / x;
		} else {
			l40:
			
			//!
			//!     Planar rotation in regular case.
			//!
			d1new = d2 + d1*x*x;
			cbar = d2 / d1new;
			sbar = x * d1 / d1new;
			d2new = d1 * cbar;
			d[m] = d1new;
			d[mp1] = d2new;
			r[m1] = sbar;
			for(col = m+2; col <= ncol; ++col) {
				m1 = m1 + 1;
				y = r[m1];
				r[m1] = cbar*r[m2] + sbar*y;
				r[m2] = y - x*r[m2];
				m2 = m2 + 1;
			}
			y = rhs[m];
			rhs[m] = cbar*rhs[mp1] + sbar*y;
			rhs[mp1] = y - x*rhs[mp1];
		}
		//!     Swap columns M and (M+1) down to row (M-1).
		
		if (m != 1) {
			pos = m;
			for (row = 1; row<= m-1; ++row) {
				x = r[pos];
				r[pos] = r[pos-1];
				r[pos-1] = x;
				pos = pos + ncol - row - 1;
			}
			
			//!     Adjust variable order (VORDER), the tolerances (TOL) and
			//!     the vector of residual sums of squares (RSS).
			
			m1 = vorder[m];
			vorder[m] = vorder[mp1];
			vorder[mp1] = m1;
			//printf("%i %i %f %f\n",m,mp1,tol[m],tol[mp1]);
			x = tol[m];
			tol[m] = tol[mp1];
			tol[mp1] = x;
			rss[m] = rss[mp1] + d[mp1] * rhs[mp1]*rhs[mp1];
		}
	}
	
	return;
}


void lsq::reordr(int *list, int n, int pos1, int& ifault) {
/*
!     ALGORITHM AS274  APPL. STATIST. (1992) VOL.41, NO. 2

  !     Re-order the variables in an orthogonal reduction produced by
  !     AS75.1 so that the N variables in LIST start at position POS1,
  !     though will not necessarily be in the same order as in LIST.
	!     Any variables in VORDER before position POS1 are not moved.		*/
	int next, i, l, j;
	
	//!     Check N.
	
	ifault = 0;
	if (n < 1 || n > ncol+1-pos1) ifault = ifault + 4;
	if (ifault != 0) return;
	
	//!     Work through VORDER finding variables which are in LIST.
	
	next = pos1;
	i = pos1;
	l10:   l = vorder[i];
	for(j = 1; j <= n; ++j) {
		if (l == list[j]) goto l40;
	}
	l30:
	i = i + 1;
	if (i <= ncol) goto l10;
	
	//!     If this point is reached, one or more variables in LIST has not
	//!     been found.
	
	ifault = 8;
	return;
	
	//!     Variable L is in LIST; move it up to position NEXT if it is not
	//!     already there.
	
	l40:
	if (i > next) 
		vmove(i, next,ifault);
	next = next + 1;
	if (next < n+pos1) goto l30;
	
	return;
}



void lsq::hdiag(double *xrow, int nreq, double& hii,int& ifault) {
	int  col, row, pos;
	double total, *wk;
	
	//!     Some checks
	
	ifault = 0;
	if (nreq > ncol) ifault = ifault + 4;
	if (ifault != 0) {
		return;
	}
	
	wk = new double[ncol+1];
	
	//!     The elements of xrow.inv(R).sqrt(D) are calculated and stored
	//!     in WK.
	
	hii = zero;
	for( col = 1; col <= nreq; ++col) {
		if (sqrt(d[col]) <= tol[col]) {
			wk[col] = zero;
		} 
		else {
			pos = col - 1;
			total = xrow[col];
			for (row = 1; row <= col-1; ++row) {
				total = total - wk[row]*r[pos];
				pos = pos + ncol - row - 1;
			}
			wk[col] = total;
			hii = hii + total*total / d[col];
		}
	}
	delete[] wk;
	return;
}		


double lsq::varprd(double *x, int nreq, double& var, int& ifault) {
	
	//!     Calculate the variance of x'b where b consists of the first nreq
	//!     least-squares regression coefficients.
	
	int  row;
	double *wk, rvarprd;
	
	//!     Check input parameter values
	
	rvarprd = zero;
	ifault = 0;
	if (nreq < 1 || nreq > ncol) ifault = ifault + 4;
	if (nobs <= nreq) ifault = ifault + 8;
	if (ifault != 0) 
		return 0.0;
	
	wk = new double[nreq+1];
	
	//!     Calculate the residual variance estimate.
	
	var = sserr / (nobs - nreq);
	
	//!     Variance of x'b = var.x'(inv R)(inv D)(inv R')x
	//!     First call BKSUB2 to calculate (inv R')x by back-substitution.
	
	bksub2(x, wk, nreq);
	for( row = 1; row <= nreq; ++row) {
		if(d[row] > tol[row]) 
			rvarprd += wk[row]*wk[row] / d[row];
	}
	rvarprd = rvarprd * var;
	
	delete[] wk;
	
	return rvarprd;
}



void lsq::bksub2(double *x, double *b, int nreq) {
	
	//!     Solve x = R'b for b given x, using only the first nreq rows and
	//!     columns of R, and only the first nreq elements of R.
	
	int  pos, row, col;
	double temp;
	
	//!     Solve by back-substitution, starting from the top.
	
	for( row = 1; row<= nreq; ++row) {
		pos = row - 1;
		temp = x[row];
		for(col = 1; col <= row-1; ++col) {
			temp = temp - r[pos]*b[col];
			pos = pos + ncol - col - 1;
		}
		b[row] = temp;
	}
}