% Figure 1.6
x=linspace(1,9,21);
y=x+0.5*sin(5*x);
net=newrb(x,y);
xf=linspace(0,10,101);
yf=xf+0.5*sin(5*xf );
ysim=sim(net,xf );
plot(xf,yf ); hold on;
plot(xf,ysim,'ro');
%-----------------------------------------------------------------------------
%
% Figure 2.1
noise=randn(1,30);
x=1:1:30; y=x+noise;
[p,s]=polyfit(x,y,1);
yfit=polyval(p,x);
plot(x,y,'+',x,x,'r',x,yfit,'b');
legend('Data','True function','Fitted function');
RMS_data=sqrt(sum((x-y).^2)/30)
RMS_fit=sqrt(sum((x-yfit).^2)/30)
mean(noise)
%-----------------------------------------------------------------------------
%
% Figure 2.2
noise=0.3*randn(1,31); noise=noise-mean(noise);
x=linspace(0,3,31); ytrue=x.^2; y=ytrue+noise;
[p,s]=polyfit(x,y,1);
yfit=polyval(p,x);
plot(x,y,'+',x,ytrue,'r',x,yfit,'b');
legend('Data','True function','Fitted function');
RMS_data=sqrt(sum((ytrue - y).^2)/31)
noise_std=std(noise)
RMS_fit=sqrt(sum((ytrue - yfit).^2)/31)
RMS_disc=sqrt(sum((y - yfit).^2)/31)
%-----------------------------------------------------------------------------
%
% Example 2.2
% Quadratic PRS for profit-per-can
x=[1.8 3.6; 2.4 4.8; 3.0 6.0; 1.5 4.6; 2.1 5.8; 2.7 7.0; 1.4 5.6; 2.0 6.8; 2.6 8.0];
p=[5.9 13.1 22.0 4.7 12.0 20.9 4.9 12.5 21.4]';
D=x(:,1); H=x(:,2);
ptrue=0.2361*pi*D.^2.*H-3.858E-4*pi^2*D.^4.*H.^2-0.1*pi*(0.5*D.^2+D.*H);
%
X=[ones(9,1) D H D.^2 D.*H H.^2];
A=X'*X;
Btrue=X'*ptrue;
B=X'*p;
btrue=A\Btrue;
b=A\B;
%
d=linspace(1.5,3.5,10);
h=linspace(3.0,7.0,10);
[D,H]=meshgrid(d,h);
P=b(1)+b(2)*D+b(3)*H+b(4)*D.^2+b(5)*D.*H+b(6)*H.^2;
Ptrue=btrue(1)+btrue(2)*D+btrue(3)*H+btrue(4)*D.^2+btrue(5)*D.*H+btrue(6)*H.^2;
surf(D, H, P); hold on; surf(D, H, Ptrue);
%-----------------------------------------------------------------------------
%
% Example 2.3
x=[0 2.5 5 7.5 10]';
y=[1.5326 1.7303 5.3714 7.2744 11.1174]';
xp=0:0.5:10; ytrue=xp;
[b,s]=polyfit(x,y,4);
yfit=polyval(b,xp);
plot(x,y,'+',xp,ytrue,'r',xp,yfit,'b');
legend('Samples','True function','Fitted function');
%-----------------------------------------------------------------------------
%
% Example 2.4
rng default;
x=linspace(0,10,20);
ytrue=x.^2;
noise=randn(1,20); noise=noise-mean(noise);
y=ytrue+noise;
bL=polyfit(x,y,1);
yL=polyval(bL,x);
sigmaL=sqrt((y-yL)*(y-yL)'/(20-2))
bQ=polyfit(x,y,2);
yQ=polyval(bQ,x);
sigmaQ=sqrt((y-yQ)*(y-yQ)'/(20-3))
plot(x,y,'+',x,x.^2,'r',x,yL,'b',x,yQ,'k');
legend('Samples','True function','Linear PRS','Quadratic PRS');
figure(2);
plot(x,y-yL,'+',x,y-yQ,'o');
legend('e_L','e_Q');
%-----------------------------------------------------------------------------
%
% Example 2.6
rng default;
x=linspace(0,30,31); noise=randn(1,31);
y1true=x; y1=y1true + noise;
y2true=0.1*x; y2=y2true + noise;
b1=polyfit(x,y1,1);
b2=polyfit(x,y2,1);
y1p=polyval(b1,x);
y2p=polyval(b2,x);
eva1= sum(abs(y1p-y1true))/31;
eva2= sum(abs(y2p-y2true))/31;
meany1=mean(y1);
meany2=mean(y2);
SSy1= sum((y1-meany1).^2);
SSy2=sum((y2-meany2).^2);
SSr1=sum((y1p-meany1).^2);
SSr2=sum((y2p-meany2).^2);
R1=SSr1/SSy1;
R2=SSr2/SSy2;
figure(1)
plot(x,y1,'+',x,y1true,'b',x,y1p,'r');
legend('Samples','True function','Linear PRS');
figure(2)
plot(x,y2,'+',x,y2true,'b',x,y2p,'r');
legend('Samples','True function','Linear PRS');
%-----------------------------------------------------------------------------
%
% Example 2.7
rng default;
x=[1:30]'; noise=randn(30,1);y=x+noise;
X=[ones(30,1) x];
b=regress(y,X)
yfit=b(1)+b(2)*x;
error=y-yfit;
sigma=sqrt(error'*error/28)
mean(noise)
std(noise)
M=X'*X; E=X*inv(M)*X';
d=diag(E)
ep=error./(1-d);
PRESS=ep'*ep
ePRESS=sqrt(PRESS/29)
%-----------------------------------------------------------------------------
%
% Example 2.10
x=[-1 -1 1]; y=[-1 1 -1];
[X,Y]=meshgrid(-1:.1:1, -1:.1:1);
Z=sqrt(.5*(1+X+Y+X.^2+Y.^2+X.*Y));
v=linspace(0.6,1.8,7)
scatter(x,y,'filled');
grid on; hold on
[C,h]=contour(X,Y,Z,v);
clabel(C,h)
%-----------------------------------------------------------------------------
%
% Example 2.13
rng default; % Control random number sequence
x=linspace(-5,5,50)'; % Equally spaces sample locations
ytrue=5*x.^3 - x.^2 + x; % True function
noise=100*randn(50,1); noise=noise-mean(noise); % Unbiased random noise
y=ytrue+noise; % Generate samples
%X=[ones(50,1) x];
X=[ones(50,1) x x.^2 x.^3]; % Design matrix
%X=[ones(50,1) x x.^2 x.^3 x.^4 x.^5 x.^6];
%X=[ones(50,1) x x.^2 x.^3 x.^4 x.^5 x.^6 x.^7 x.^8 x.^9 x.^10];
[b,bint,r,rint,stats]=regress(y,X); % Fitting PRS
yfit=X*b; % Prediction at sample points
s=sqrt(r'*r/(50-size(X,2)); % Standard error of noise
sy=s.*diag(sqrt(X*inv(X'*X)*X')); %Standard error of prediction
plot(x,y,'+',x,yfit,'b',x,yfit-2*sy,'k',x,yfit+2*sy,'k'); % Plotting
%-----------------------------------------------------------------------------