function [e]=errors(koeff,z)
	X=[0:0.4:2]*%pi;
	rand("seed",1); 
	noise=.15*(rand(1,length(X))-0.5);
	Y1=sin(X)+noise;
endfunction

function [e]=errors(koeff,z)
  [m,n] = size(koeff);
  e=z(2);
  for i=1:1:m
    e=e-koeff(i)*z(1).^(i-1);
  end	
endfunction

function [g]=derrors(koeff,z)
  [m,n] = size(koeff);
  g=koeff;
  for i=1:1:m
    g(i)=-z(1).^(i-1);
  end	
endfunction

function [koeff]=generalfitalgo(X,Y1,t,lambda)
	measdata=[X;Y1]; // measurement data as a matrix
	if ( ~exists('lambda') ) then
  		lambda = 0
	end

	//koeff0=[1;2;3]; // This is a starting guess
	koeff0=rand(t+1,1);
	// And here the command datafit:
	[koeff,SSE]=datafit(errors,measdata,koeff0); //SSE=SumSquaredErrors
endfunction

function [koeff]=ridgefitalgo(X,Y1,t,lambda)
	[m1,n1] = size(X);
	if ( ~exists('lambda') ) then
  		lambda = 0
	end

	phi=zeros(n1,t+1);
	for i=1:1:n1
		for j=1:1:t+1
    			phi(i,j)=X(i)^(j-1);
		end
	end

	//Matrix inverse gives inaccurate results for large values of "t" 
	koeff=(phi'*phi + lambda*eye(t+1,t+1))^(-1)*(phi'*Y1'); 	
endfunction

function [koeff]=specificfitalgo(X,Y1,t)
	[m1,n1] = size(X);
	
	phi=zeros(n1,t+1);
	for i=1:1:n1
		for j=1:1:t+1
    			phi(i,j)=X(i)^(j-1);
		end
	end

	//Matrix inverse gives inaccurate results for large values of "t" 
	//koeff=(phi'*phi)\(phi'*Y1');	
	//[koeff,ker]=linsolve((phi'*phi),-(phi'*Y1'));
	
	koeff=lsq(phi,Y1');
	
	//Experimental -- did not work so far
	//rols=ols2(Y1',phi); koeff=rols('beta');
endfunction

function [koeff,sizex] = findfit(t,lambda,ridge)

	close;
	//exec('grocer/loader.sce',-1)
	if ( ~exists('lambda') ) then
		lambda = 0
	end
	[X,Y1]=getRegressionData();
        if (ridge == 1) then
		koeff=ridgefitalgo(X,Y1,t,lambda);
	else
		koeff=specificfitalgo(X,Y1,t);
		//koeff=generalfitalgo(X,Y1,t,lambda);
	end
	

	// Now the result is in koeff and the fitted polynomial becomes 
	Xplot=[0:0.02:2]*%pi
	[m1,n1] = size(Xplot);
	y=ones(m1,n1)*koeff(1);
	for i=2:1:t+1
	    y=y+koeff(i)*Xplot^(i-1);
	end
	//plot2d(X',[Y1' y'], style=[-9 2], leg="estimated@data");
	plot(Xplot',sin(Xplot'),"g",X',Y1',"b+",Xplot',y',"r");
	hl=legend(['cos(t)';'cos(2*t)';'cos(3*t)']);
	koeff;
	meanwt=mean(koeff);
	varwt=variance(koeff);
	print(%io(2),meanwt,varwt);
	[dummy,sizex]=size(X);	
endfunction

function [X,Y1] = getRegressionData()
	X=[0:0.4:2]*%pi;
	rand("seed",1); 
	noise=.15*(rand(1,length(X))-0.5);
	Y1=sin(X)+noise;
endfunction

function [rmsvalue]=rms(koeff,measdata)
  [m,n] = size(koeff);
  [dummy,trainsize]=size(measdata);
  e=measdata(2,:);
  for i=1:1:m
    e=e-koeff(i)*measdata(1,:).^(i-1);
  end	
  rmsvalue=sqrt(sum((e.*e))/trainsize);
endfunction

///////////////////////////////////////////////////////////////////

function [trainRMS] = plotRMSAgainstDegree(numpoints)
	close;
	//exec('grocer/loader.sce',-1)
	if ( ~exists('lambda') ) then
		lambda = 0
	end
	trainX=[0:0.4:2]*%pi;
	rand("seed",1); 
	noise=.15*(rand(1,length(trainX))-0.5);
	trainY1=sin(trainX)+noise;
	trainRMS=getRMSAgainstDegree(numpoints,trainX,trainY1);

	testX=[1:0.4:3]*%pi; 
	noise=.15*(rand(1,length(testX))-0.5);
	testY1=sin(testX)+noise;
	testRMS=getRMSAgainstDegree(numpoints,testX,testY1);

	plot2d([1:1:numpoints]',[trainRMS' testRMS'], style=[4 5], leg="test@train");
endfunction

function [rmsvalues] = getRMSAgainstDegree(numpoints,X,Y1)
	exec('grocer/loader.sce',-1)
	measdata=[X;Y1]; //measurement data as a matrix
	rmsvalues=zeros(1,numpoints);
	for t=1:1:numpoints
		koeff=specificfitalgo(X,Y1,t);
		rmsvalues(t)=rms(koeff,measdata);
	end
endfunction

/////////////////////////////////////////////////////

function [trainRMS,testRMS] = plotRMSAgainstLambda(t,numlambda)
	close;
	if ( ~exists('lambda') ) then
		lambda = 0
	end
	base=-40	;
	trainX=[0:0.4:2]*%pi;
	rand("seed",1); 
	noise=.15*(rand(1,length(trainX))-0.5);
	trainY1=sin(trainX)+noise;
	trainRMS=getRMSAgainstLambda(numlambda,trainX,trainY1,t);
	
	testX=[1:0.4:3]*%pi; 
	noise=.15*(rand(1,length(testX))-0.5);
	testY1=sin(testX)+noise;
	testRMS=getRMSAgainstLambda(numlambda,testX,testY1,t);

	plot2d([base:1:(base-1+numlambda)]',[trainRMS' testRMS'], style=[4 5], leg="test@train");
endfunction

function [rmsvalues] = getRMSAgainstLambda(numlambda,X,Y1,t)
	measdata=[X;Y1]; //measurement data as a matrix
	rmsvalues=zeros(1,numlambda);
	for lnlambda=1:1:numlambda
		koeff=ridgefitalgo(X,Y1,t,%e^(base+lnlambda));
		rmsvalues(lnlambda)=rms(koeff,measdata);
	end
endfunction
////////////////////////////////////////////////////////


