(* :Title: GLib *)

(* :Author: Robert G. Buice,Jr *)

(* :Summary:
	This package provides a set of functions that
	are commonly required in the analysis of spectral
	data files.  These routines were developed by the
	Analytical Spectroscopy Group at the University
	of Kentucky, College of Pharmacy.
*)

(* :Context: GLib` *)

(* :Package Version: 0.2.4 *)

(* :Copyright: Copyright 1994,1995 by Robert G. Buice, Jr *)

(* :History: *)

(* :Keywords: *)

(* :Source: None. *)

(* :Warning: *)

(* :Mathematica Version: 2.0 *)

(* :Limitation: None known. *)

(* :Discussion: *)

Needs["Statistics`ContinuousDistributions`"]
Needs["Statistics`DataManipulation`"]
Needs["Statistics`DescriptiveStatistics`"]
Needs["Statistics`LinearRegression`"]
Needs["Statistics`InverseStatisticalFunctions`"]
Needs["Statistics`DiscreteDistributions`"]

GLib::usage =
"GLib is a package that serves as a chemometric toolbox for analysis
of real data sets."

MSCorrect::usage =
"newspec = MSCorrect[tnspec]. The scatter corrected spectra
 (newspec) are returned from the input spectra (tnspec)."

FixMSCorrect::usage = 
"newspec = FixMSCorrect[tnspec, mnspec].  The scatter corrected
spectra (newspec) are returned from the input spectra (tnspec).
The correction is calculated from a user supplied reference 
spectrum (mnspec) with the same number of wavelengths as tnspec."

Replica::usage =
"{btrain,cnter} = Replica[tnspec, reps]. A List is returned which contains
first the replicates and second the center of the replicates."

CenterSpec::usage = 
"newspec = CenterSpec[tnspec].  Returns the average or center 
(newspec) of an input set of spectra (tnspec)."

SOBeast::usage = 
"SOBeast[btrain_List,ctn_List,testspec_List].  Returns the 
correlation between two groups in hyperspace.  (testspec) is the
test set and (btrain) is the set of replicates of the calibration 
set, having center (ctn).  The groups must be of the same size."

Limits::usage = 
"confid = Limits[tnspec].  The confidence limit (confid) on the
calibration set (tnspec) is returned."

Runmed::usage =
"outspec = Runmed[spec, n].  A running median filter is applied
to the input spectra (spec) in windows of (n) points to porduce the 
filtered spectra(outpec).  This routine is to be used on one dimensional
arrays only (a single spectra)."

Runmeds::usage =
"outspec = Runmed[spec, n].  A running median filter is applied
to the input spectra (spec) in windows of (n) points to porduce the 
filtered spectra(ouspec).  This routine operates on a two dimensional
array (multiple spectra)."

CorrelRows::usage =
"correlation_matrix = CorrelRows[spectra].  This routine
calculates the correlations between the rows of a set of spectra
and returns them to correlation_matrix."

CorrelCols::usage =
"correlation_matrix = CorrelCols[spectra]  This routine
calculates the correlations between the columns of a set of 
spectra and returns them to correlation_matrix."

QBeast::usage =
"QBeast[tnspec_, btrain_, newspec_,cnter_,radfrac_,sensitiv_]  This routine 
calculates the beast standard deviation between the group
tnspec and the single spectra (newspec)."

StdCols::usage =
"StdCols[tnspec] returns the standard deviation calculated with
n rather than n-1."

MatchIndex::usage =
"MatchIndex[spec1_, spec2_] returns the match index between the 
two spectra spec1 and spec2. Reed and Wong, Applied Spectroscopy,
Vol. 20, 320-325."

SEP::usage =
"SEP[res] returns the standard error of performance using 
cross-validation residuals from a fit."

Mahalanobis::usage =
"Mahalanobis[tnspec_List, spec_List] calculates the distance in SD's
from a single sample (spec) to a homogeneous calibration set (tnspec).
If the number of wavelengths exceeds the number of samples significant
errors are introduced."

Kubelka::usage =
"Kubleka[spec_List] takes 1/R data in (spec) and converts it to 
z (Kubelka-Munk) scores via (1-R)^2/(2R).  H. Martens, Applied
Spectroscopy, 1988, 42(5), 722-728."

Spec3M::usage =
"Spec3M[spec_List]  returns the average of every three rows of the 2D
array spec."

XValidQ::usage =
"XValidQ[tnspec_List] returns the t-statistics for leave-one-out
cross-validation of the calibration set (tnspec) using the quantile
beast metric."

FSOBeast::usage =
"FSOBeast[btrain_List,ntn_Integer,ctn_List,testspec_List].  Returns the 
correlation between two groups in hyperspace.  (testspec) is the
test set and (btrain) is the set of replicates of the calibration 
set, having center (ctn). FSOBeast does not require that the groups be
of the same size."

BigPComp::usage =
"BigPComp[spectra_List].  Returns a list of three elements.  The first is
the principal components and the second is the loadings matrix and
the third is the zscored input spectra."

OLSQ::usage =
"OLSQ[indep_List, dep_List].  Returns a list of two elements.  The first
is the predicted values, the second is the residuals.  This routine is 
designed to emulate the Speakeasy OLSQ command which is a quick way of 
doing a regression of a set of spectra against concentrations."


Begin["Private`"]

MSCorrect[tnspec_List]:= Module[{norows,nocols,temp,a,b,mnspec,i,
							coeff,x,output},
	norows = (Dimensions[tnspec])[[1]];
	nocols = (Dimensions[tnspec])[[2]];
	temp = Array[a, {2,nocols}];
	outspec = Array[b, {norows, nocols}];
	mnspec = (Plus @@ tnspec)/(Dimensions[tnspec])[[1]];
	temp[[1]] = mnspec;
	For[i = 1, i < norows + 1, ++i,
		temp[[2]] = tnspec[[i]];
		output = Regress[Transpose[temp], {1,x}, x, 
			     OutputControl->NoPrint, 
			     OutputList->BestFitCoefficients];
		coeff = BestFitCoefficients /. output;
		Print[coeff];
		outspec[[i]] = (tnspec[[i]] - coeff[[1]])/coeff[[2]];
	];
	Return[outspec]]
	
FixMSCorrect[tnspec_List, mnspec_List]:= Module[{norows,nocols,
	temp,a,b,i,output,coeff,x},
	norows = (Dimensions[tnspec])[[1]];
	nocols = (Dimensions[tnspec])[[2]];
	temp = Array[a, {2,nocols}];
	outspec = Array[b, {norows, nocols}];
	temp[[1]] = mnspec;
	For[i = 1, i < norows + 1, ++i,
		temp[[2]] = tnspec[[i]];
		output = Regress[Transpose[temp], {1,x}, x, 
			     OutputControl->NoPrint, 
			     OutputList->BestFitCoefficients];
		coeff = BestFitCoefficients /. output;
		Print[coeff];
		outspec[[i]] = (tnspec[[i]] - coeff[[1]])/coeff[[2]];
	];
	Return[outspec]]
	
Replica[tnspec_List,reps_Integer] := 
	Module[{norows,nocols,picks,a,i,bsamp,x,c,btrain,cnter},
	norows = (Dimensions[tnspec])[[1]];
	Print["Replica output:"];
	Print["---------------"];
	Print["Return[[1]] = the replicates"];
	Print["Return[[2]] = the center of the replicates"];
	nocols = (Dimensions[tnspec])[[2]];
	btrain = Array[a, {reps, nocols}];
	picks = Table[Random[Integer,{1,norows}], {reps},{norows}];
	For[i = 1, i <= reps, ++i,
		bsamp = Table[tnspec[[(picks[[i]])[[x]]]], {x, 1, norows}];
		btrain[[i]] = (Plus @@ bsamp)/((Dimensions[bsamp])[[1]])
	];
	cnter = ((Plus @@ btrain) /((Dimensions[btrain])[[1]]));
	outs = Array[c,{2}];
	outs[[1]] = btrain;
	outs[[2]] = cnter;
	Return[outs]]
	
CenterSpec[tnspec_List]:= (Plus @@ tnspec) /(Dimensions[tnspec])[[1]];

	
SOBeast[btrain_List,ctn_List,testspec_List] :=
	Module[{btest,ctest,m,p,st,sx,trim,tecdf,ttcdf,temp,output},
	btest = Replica[testspec,100];
	m = (Dimensions[btrain])[[1]]; m *= 2;
	p = 0.1;
	st = Table[N[Sqrt[Plus @@ ((btrain[[x]]-ctn)^2)]], {x,m/2}];
	sx = Table[N[Sqrt[Plus @@ ((btest[[x]]-ctn)^2)]], {x,m/2}];
	trim = Range[(m*p)+1,m-(m*p)];
	tecdf = Sort[Join[st,sx]];
	ttcdf = Sort[Join[st,st]];
	ecdf = Table[tecdf[[trim[[x]]]],{x,Length[trim]}];
	tcdf = Table[ttcdf[[trim[[x]]]],{x,Length[trim]}];
	temp = Array[a, {2,Length[trim]}];
	temp[[1]]=tcdf; temp[[2]]=ecdf;
	temp = Transpose[temp];
	ListPlot[temp,PlotRange->All,Frame->True,FrameLabel->{"ecdf","tcdf"},
										PlotLabel->"QQ Plot",Axes->False];
	output = Regress[temp,{1,x},x,
							OutputControl->NoPrint,
							OutputList->RSquared];
	rsqr = RSquared /. output;
	Print[rsqr];
	Return[rsqr]]
	
FSOBeast[btrain_List,ntn_Integer,ctn_List,testspec_List] :=
	Module[{btest,ctest,m,p,st,sx,trim,tecdf,ttcdf,temp,output},
	btest = Replica[testspec,100];
	m = (Dimensions[btrain])[[1]]; m *= 2;
	p = 0.1;
	ntest = (Dimensions[testspec])[[1]];
	st = Table[N[Sqrt[Plus @@ (((btrain[[x]]-ctn)(Sqrt[ntn]))^2)]], {x,m/2}];
	sx = Table[N[Sqrt[Plus @@ (((btest[[x]]-ctn)(Sqrt[ntest]))^2)]], {x,m/2}];
	trim = Range[(m*p)+1,m-(m*p)];
	tecdf = Sort[Join[st,sx]];
	ttcdf = Sort[Join[st,st]];
	ecdf = Table[tecdf[[trim[[x]]]],{x,Length[trim]}];
	tcdf = Table[ttcdf[[trim[[x]]]],{x,Length[trim]}];
	temp = Array[a, {2,Length[trim]}];
	temp[[1]]=tcdf; temp[[2]]=ecdf;
	temp = Transpose[temp];
	ListPlot[temp];
	output = Regress[temp,{1,x},x,
							OutputControl->NoPrint,
							OutputList->RSquared];
	rsqr = RSquared /. output;
	Print[rsqr];
	Return[rsqr]]
	
Limits[training_List] := 
	Module[{norows,nocols,btrain,ncors,cors,ctr,picks,check,
			rc,mn,sd},
	norows = (Dimensions[training])[[1]];
	nocols = (Dimensions[training])[[2]];
	btrain = Replica[training, 100];
	ncors = 10;
	Print["ncors = ",ncors];
	cors = Array[a, {ncors}];
	For[ctr = 1, ctr <= ncors, ++ctr,
		picks = Table[Random[Integer,{1,norows}], {norows}];
		check = Table[training[[picks[[x]]]], {x,norows}];
		Dimensions[check];
		rc = SOBeast[btrain,CenterSpec[btrain],check];
		cors[[ctr]] = rc;
	];
	cors = Sort[cors];
	mn = Mean[cors]; Print[mn];
	sd = StandardDeviation[cors]; Print[sd];
	confid = mn - (2.0537sd); Print[confid];
	Return[confid]]
	
Runmed[spec_List, n_Integer] := 
	Module[{mshift,shifted,a,b,i,j,chop,rchop},
	mshift = Array[b, {Length[spec]}];
	shifted = Array[a, {n, Length[spec]}];
	shifted[[1]] = spec;
	For[i = 2, i < n + 1, ++i,
		shifted[[i]] = RotateRight[spec, i - 1];
		Do[shifted[[i,y]] = 0, {y,1,i-1}];	
	];
	For[j = 1, j < Length[spec] + 1, ++j,
		mshift[[j]] = Median[Column[shifted, j]];
	];
	chop = Floor[n/2];
	outspec = RotateLeft[mshift,chop];
	Do[outspec[[-z]]= 0, {z,1,chop}];
	rchop = Range[(Length[spec]-chop),Length[spec]];
	Do[outspec[[rchop[[x]]]] = spec[[rchop[[x]]]], 
				{x,1,Length[rchop]}];
	Return[outspec]]
	
Runmeds[ins_List,n_Integer] := 
	Module[{r,c,holdit,a,i},
	r = (Dimensions[ins])[[1]];
	c = (Dimensions[ins])[[2]];
	outs = Array[a, {r,c}];
	For[i = 1, i < r +1, ++i,
		holdit = Runmed[ins[[i]],n];
		outs[[i]] = holdit;
	];
	Return[outs]]
	
CorrelRows[spect_List] := 
	Module[{zscore,spec,junk,cmat,x},
	zscore = Function[x, (x - Mean[x])/StandardDeviation[x]];
	spec = N[Map[zscore, spect]];
	junk = spec.(Transpose[spec]);
	cmat = junk/junk[[1,1]];
	Return[cmat]]
	
CorrelCols[spect_List] := 
	Module[{zscore,spec,junk,x},
	zscore = Function[x, (x - Mean[x])/StandardDeviation[x]];
	spec = Transpose[spect];
	spec = N[Map[zscore, spec]];
	junk = spec.(Transpose[spec]);
	cmat = junk/junk[[1,1]];
	Return[cmat]]
	
QBeast[tnspec_, btrain_, newspec_,cnter_,radfrac_:0.35,sensitiv_:0.01] := 
	Module[{b,so2,sor,s2r,sub,area,radial,project,x,qdist,qr,radii,
			lindex,uindex,sd,sds,alpha,za,tcenter,i,cso2,csor,cs2r,csub,
			carea,cradial,cproject,n,md,fdist,index,zelement,zo,lowind,
			upind,lowlim,uplim,euc,fac,erd},
	Print["radfrac: ", radfrac];
	Print["sensitiv:" , sensitiv];
	b = (Dimensions[btrain])[[1]];
	so2 = Sqrt[Plus @@ ((newspec - cnter)^2)];
	sor = Table[Sqrt[Plus @@ ((btrain[[x]] - cnter)^2)], {x,1,b}];
	s2r = Table[Sqrt[Plus @@ ((btrain[[x]] - newspec)^2)], {x,1,b}];
	sub = (so2 + sor + s2r) / 2;
	area = Sqrt[sub (sub - so2) (sub - sor) (sub - s2r)];
	radial = (2 area)/so2;
	project = Sqrt[sor^2 - radial^2];
	Do[
		If[	((so2^2 + sor[[x]]^2) < s2r[[x]]^2),
			project[[x]] *= -1],
	{x,1,Length[s2r]}];
	qdist = project;
	qr = Sort[radial]; 
	radii = qr[[Floor[radfrac*b]]];
	Do[
		If[radial[[x]] > radii, qdist[[x]] = 0],
	{x,1,Length[radial]}];
	qdist = Cases[qdist, x_ /; x != 0];
	qdist = Sort[qdist];
	lindex = Floor[0.16 Length[qdist]];
	uindex = Floor[0.84 Length[qdist]];
	If[Length[qdist] < 50, Print["** Need more replicates in hypercylinder **"]]; 
	If[Length[qdist] < 50, Return[]];
	sd = StandardDeviation[qdist] Sqrt[(Dimensions[tnspec])[[1]]];
	sds = N[Sqrt[Plus @@ ((cnter - newspec)^2)]/sd];
	Print["sds: ", sds];
	alpha = CDF[NormalDistribution[0,1], -1];
	za = N[Sqrt[2] InverseErf[-1 + 2alpha]];
	tcenter = Array[j,{(Dimensions[tnspec])[[2]]}];
	For[i = 1, i < (Dimensions[tnspec])[[2]] + 1, ++i,
		tcenter[[i]] = Median[Column[tnspec,i]];
	];
	Label[adjust];
	cso2 = so2;
	csor = Sqrt[Plus @@ ((tcenter - cnter)^2)];
	cs2r = N[Sqrt[Plus @@ ((tcenter - newspec)^2)]];	
	csub = (cso2 + csor + cs2r)/2;
	carea = N[Sqrt[csub (csub-cso2) (csub - csor) (csub - cs2r)]];
	cradial = (2carea)/cso2;
	cproject = N[Sqrt[csor^2 - cradial^2]];
	If[((so2^2 + csor^2) > cs2r^2), cproject *= -1];
	n = Length[qdist];
	If[IntegerQ[n/2],
		md = (qdist[[n/2]] + qdist[[n/2 + 1]])/2,
		md = qdist[[n/2 + 0.5]]
	];
	cproject *= sensitiv;
	cproject += md;
	fdist = qdist - cproject;
	index = Range[Length[fdist]];
	If[cproject > Max[qdist], zelement = Length[qdist] - 1];
	If[cproject > Max[qdist], Goto[setcor]];
	If[cproject < Min[qdist], zelement =1];
	If[cproject < Min[qdist], Goto[setcor]];
	zelement = (Flatten[Position[Abs[fdist], Min[Abs[fdist]]]])[[1]];
	Label[setcor];
	zo = N[(Sqrt[2] InverseErf[-1 + 2 (zelement/Length[qdist])])];
	If[(Abs[2zo] > Abs[za]), Print["** decrease skew sensitivity **"]];
	If[(Abs[2zo] > Abs[za]), Return[]];
	lowind = Floor[CDF[NormalDistribution[0,1], (2zo+za)] Length[qdist]];
	upind = Floor[CDF[NormalDistribution[0,1], (2zo-za)] Length[qdist]];
	If[lowind < 2, Print["** warning too few replicates **"]];
	If[upind > Length[qdist] -2, Print["** warning too few replicates **"]];
	If[lowind < 1, lowind =1];
	If[upind > Length[qdist], upind = Length[qdist]];
	lowlim = qdist[[lowind]];
	uplim = qdist[[upind]];
	euc = N[Sqrt[Plus @@ ((cnter - newspec)^2)]];
	fac = Abs[N[Sqrt[2] InverseErf[-1 + 2alpha]]];
	erd = N[Sqrt[(Dimensions[tnspec])[[1]]]];
	If[(Abs[2zo] > Abs[fac]), Print["** skew correction exceeds replicates **"]];
	sdskew = N[euc/((uplim/fac)erd)];
	Print["sdskew: ", sdskew];
	Return[sdskew]]
	
StdCols[useit_List] := Module[{s},
	useit = Transpose[useit];
	dosd = Function[xxx, StandardDeviation[xxx]];
    s = N[Map[dosd,useit]];
	Return[s]]
	
MatchIndex[spec1_List, spec2_List] := 
 	Plus @@ (spec1 spec2) / Sqrt[(Plus @@ (spec1^2))(Plus @@ (spec2^2))]
	
SEP[res_] := Sqrt[(Plus @@ (res^2)) / (Length[res] - 1)]
	
Mahalanobis[tnspec_List, spec_List] := 
	Module[{gmean,t,m2,d,temp},
	gmean = (Plus @@ (tnspec)) / (Dimensions[tnspec])[[1]];
	t = spec - gmean;
	temp = Transpose[tnspec];
	m2 = Inverse[temp.(Transpose[temp])];
	d = Sqrt[t.m2.t];
	Print["d = ",d];
	Return[d]]
	
Kubelka[spec_List] := 
	Module[{r,newspec},
	r = 10^(-spec);
	newspec = ((1-r)^2)/(2r);
	Return[newspec]]
	
Spec3M[spec_List] := 	
	Module[{r,c,a,taker,takec,i,temp,new},
	r = (Dimensions[spec])[[1]];
	c = (Dimensions[spec])[[2]];
	If[!IntegerQ[r/3], 
		Print["Invalid array size!"];
		Return[]];
	new = Array[a, {r/3,c}];
	taker = {1,2,3};
	takec = Range[c];
	For[i = 1, i <= r/3, ++i,
		temp = spec[[taker,takec]];
		new[[i]] = (Plus @@ temp)/3;
		taker += 3;
	];
	Return[new]]
	
XValidQ[tnspec_] := 
	Module[{nsamp,dis,i,trn,bdiststuff,bdist,ct,dd,avdis,df,tfunc,temp},
	nsamp = (Dimensions[tnspec])[[1]];
	dis = Array[a, {nsamp}];
	For[i = 1, i < nsamp + 1, ++i,
		Print[i];
		trn = Delete[tnspec, i];
		Print["replicating..."];
		bdiststuff = Replica[trn, 200];
		bdist = bdiststuff[[1]];
		ct = bdiststuff[[2]];
		Print["doing qb"];
		dd = QBeast[trn,bdist,tnspec[[i]],ct];
		dis[[i]] = dd;
	];
	avdis = Mean[dis];
	Print["avdis = ", avdis];
	Print[Min[dis],Max[dis]];
	Print["Range = ", Max[dis] - Min[dis]];
	Print["t-corrected 3 SD limit ="];
	df = nsamp - 2;
	If[df > 30, df = 30];
	Print[dis];
	tfunc = Function[x,1 - CDF[StudentTDistribution[df],x]];
	temp = Map[tfunc,dis];
	Print[temp];
	Print["Elements failing xvalidation:"];
	Print[Position[temp, x_ /; x < 0.0026]];
	Return[dis]]
	
BigPComp[spectra_] := Module[{tspec,zscore,zspec,cmat,evec,eval,sqrteval,
						dump,loadings,rz,cz,pri,cumvar,teval,roots,
						rootnum,propvar,junk,dtable,take,ctr},
	Print["Copyright 1995 by Robert G. Buice, Jr"];
	Print["Returns a list with the first item being the pc's"];
	Print["and the second being the loadings matrix and the "];
	Print["third being the z-scored input spectra."];
	tspec = Transpose[spectra];
	zscore = Function[x, (x - Mean[x])/StandardDeviation[x]];
	zspec = N[Map[zscore, tspec]];
	cmat = CorrelCols[spectra];
	evec = Eigenvectors[cmat];
	evec = Transpose[Reverse[evec]];
	eval = Eigenvalues[cmat];
	eval = Chop[eval];
	eval = Reverse[eval];
	eval = Cases[eval, x_ /; x>0];
	sqrteval = Sqrt[eval];
	sqrteval = DiagonalMatrix[sqrteval];
	If[((Dimensions[evec])[[2]] - (Dimensions[sqrteval])[[2]]) == 0, Goto[nodump]];
		dump = ((Dimensions[evec])[[2]] - (Dimensions[sqrteval])[[2]]);
		Print["evec has lost column numbers:"];
		Print [dump];
		evec = Transpose[evec];
		evec = Drop[evec,dump];
		evec = Transpose[evec];
	Label[nodump];
	loadings = evec.sqrteval;
	loadings = Transpose[Reverse[Transpose[loadings]]];
	rz = (Dimensions[zspec])[[1]];
	cz = (Dimensions[zspec])[[2]];
	pri = PseudoInverse[loadings].zspec;
	pri = Transpose[pri];
	cumvar = Array[c, {Length[eval]}];
	teval = Reverse[eval];
	For[x = Length[teval], x >= 1, --x,
		cumvar[[x]] = teval[[x]];
		For[y = x-1, y >= 1, --y,
			cumvar[[x]] += teval[[y]];
		];
	];
	cumvar = cumvar/(Dimensions[spectra])[[2]];
	roots = Reverse[eval];
	rootnum = Range[1,Length[eval]];
	propvar = Reverse[eval]/(Dimensions[spectra])[[2]];
	take = rootnum;
	ctr = -1;
	Do[
		If[cumvar[[x]] == 1, 
	 		ctr += 1;
	 		Print["ctr:", ctr];
	 		col = x - ctr;
	 		Print["Dropped column:", col];
			take = Drop[take, {col}]
		],
	{x,1,Length[cumvar]}];
	cumvar = cumvar[[take]];
	roots = roots[[take]];
	rootnum = rootnum[[take]];
	propvar = propvar[[take]];
	junk = {rootnum,roots,propvar,cumvar};
	junk = Transpose[junk];
	dtable = TableForm[junk,TableHeadings->{None,{"Rootnum","Roots",
			"Propvar","Cumvar"}}];
	FontForm[dtable, {"Courier-Bold", 14}];
	Print[dtable];
	pri = pri[[Range[(Dimensions[pri])[[1]]],take]];
	outs = Array[f, {3}];
	outs[[1]] = pri;
	outs[[2]] = loadings;
	outs[[3]] = Transpose[zspec];
	Return[outs]]
	
	OLSQ[indep_List, dep_List] := Module[{ncols,nrows,usedat,tindep,
					xxx,varstouse,mymodel,myvars,output,predict,resids,
					ptable,rsqr,arsqr},
	Print["Copyright 1996 Robert G. Buice, Jr."];
	Print["Returns a list with two elements:"];
	Print["The first is the predicted results."];
	Print["The second is the fit residuals"];
	Print["indep is:", Dimensions[indep]];
	Print["dep is:", Dimensions[dep]];
	ncols = (Dimensions[indep])[[2]];
	Print["ncols:"ncols];
	nrows = (Dimensions[indep])[[1]];
	Print["nrows:" nrows];
	usedat = Array[ppp, {nrows, ncols + 1}];
	Print["usedat is:", Dimensions[usedat]];
	usedat = Transpose[usedat];
	tindep = Transpose[indep];
	For[xxx = 1, xxx <= ncols, xxx++,
		usedat[[xxx]] = tindep[[xxx]];
	];
	usedat[[ncols + 1]] = dep;
	usedat = Transpose[usedat];
	varstouse = {1,c1,c2,c3,c4,c5,c6,c7,c8,c9,c10,
				 c11,c12,c13,c14,c15,c16,c17,c18,c19,c20,
				 c21,c22,c23,c24,c25,c26,c27,c28,c29,c30,
				 c30,c31,c32,c33,c34,c35,c36,c37,c38,c39,
				 c40,c41,c42,c43,c44,c45,c46,c47,c48,c49,
				 c50,c51,c52,c53,c54,c55,c56,c57,c58,c59,
				 c60,c61,c62,c63,c64,c65,c66,c67,c68,c69,
				 c70,c71,c72,c73,c74,c75,c76,c77,c78,c79,
				 c80,c81,c82,c83,c84,c85,c86,c87,c88,c89,
				 c90,c91,c92,c93,c94,c95,c96,c97,c98,c99,
				 c100};
	mymodel = varstouse[[Range[ncols + 1]]];
	myvars = Drop[mymodel, 1];
	output = Regress[usedat, mymodel, myvars, 
		OutputList->{PredictedResponse, FitResiduals,
			         ParameterTable, RSquared, 
			         AdjustedRSquared}];
	predict = PredictedResponse /. output;
	resids = FitResiduals /. output;
	ptable = ParameterTable /. output;
	Print[ptable];
	rsqr = RSquared /. output;
	arsqr = AdjustedRSquared /. output;
	Print["RSquared:", rsqr];
	Print["AdjustedRSquared:", arsqr];
	o11 = Array[rrr, {2}];
	o11[[1]] = predict;
	o11[[2]] = resids;
	Return[o11];
	]

	End[ ]