function nregress1

global xd yd

% read in data for traffic speed
% density (veh/km), speed (km/hr)
load amst.txt
density=amst(:,1);
speed=amst(:,2);
xd=density;
yd=speed;

% find v1 and v2, and then determine resulting curve
v2=max(yd); v1=-5;
sol = fminsearch(@E, [v1 v2],optimset('MaxFunEvals',1000000,'MaxIter',1000000,'TolX',1e-6));
v1=sol(1)
v2=sol(2)
x=linspace(min(xd),max(xd),100);
for i=1:100
    y(i)=f(x(i),v1,v2);
end

% plot data and regression curve
clf
% get(gcf)
set(gcf,'Position', [25 1115 658 230])
co = [0 0 1;
      0 0.5 0;
      1 0 0;
      0 0.75 0.75;
      0.75 0 0.75;
      0.75 0.75 0;
      0.25 0.25 0.25];
set(groot,'defaultAxesColorOrder',co)

hold on
plot(xd,yd,'o','LineWidth',1.5,'MarkerSize',6,'MarkerEdgeColor','g')
box on
grid on
xlabel('Density (cars/km)')
ylabel('Speed (km/hr)')
set(gca,'FontSize',16,'FontWeight','bold') 

pause

plot(x,y,'r','LineWidth',1.6)
legend({' Data',' Nonlinear Fit'},'Location','NorthEast','FontSize',16,'FontWeight','bold')

pause

% find V1 and V2
nd=length(xd);
for i=1:nd
    XD(i)=xd(i)^4;
    YD(i)=1/yd(i);
end
SX=sum(XD); 
SX2=dot(XD,XD);
SY=sum(YD); 
SXY=dot(XD,YD);
D=nd*SX2-SX^2;
V1=(SX2*SY-SX*SXY)/D;
V2=(-SX*SY+nd*SXY)/D;
vv1=1/V1
vv2=V2/V1
for i=1:100
    yy(i)=f(x(i),vv1,vv2);
end
plot(x,yy,'b','LineWidth',1.5)
legend({' Data',' Nonlinear Fit',' Transformed Fit'},'Location','NorthEast','FontSize',16,'FontWeight','bold')

pause

%%%%%%  find min using relative error
A=0;  B=0;  D=0; B1=0;  B2=0;
for i=1:nd
    A=A+1/YD(i)^2;
    B=B+XD(i)/YD(i)^2;
    D=D+XD(i)^2/YD(i)^2;
    B1=B1+1/YD(i);
    B2=B2+XD(i)/YD(i);
end
VV=[[A B];[B D]]\[B1; B2];
vvv1=1/VV(1)
vvv2=vvv1*VV(2)
for i=1:100
    yyy(i)=f(x(i),vvv1,vvv2);
end
plot(x,yyy,'k','LineWidth',1.5)
legend({' Data',' Nonlinear Fit',' Transformed Fit',' Scaled Transformed Fit'},'Location','NorthEast','FontSize',16,'FontWeight','bold')




% error function
function s=E(v)
global xd yd
v1=v(1); v2=v(2);
s=0;
for i=1:length(xd)
	s = s + (f(xd(i),v1,v2) - yd(i))^2;
end

% model function
function g=f(x,v1,v2)
g=v1/(1+v2*x^4);