Live Script edama_11
This Live Script supports Chapter 11 of Environmental Data Analysis with MATLAB or PYTHON, Third Edition, by WIlliam Menke
The script naming system has the form edama_CC_SS, where CC is the chapter number and SS is the script number.
% edama_11_01: example of a Taylor Series
% using the function
% 1/(1-x) = 1 + x + x^2 + x^3 + ...
clear all;
N = 101;
tmin = -1.2;
tmax = 0.8;
Dt = (tmax-tmin)/(N-1);
t = tmin+Dt*[0:N-1]';
figure(1);
clf;
set(gca,'LineWidth',2);
hold on;
axis([tmin,tmax,0,2]);
hold on;
xlabel('t');
ylabel('y(t)');
y = 1./(1-t);
plot(t,y,'-','LineWidth',7,'Color',[0.8,0.8,0.8]);
t0=0;
y0=1./(1-t0);
plot(t0,y0,'ko','LineWidth',2);
y1 = 1+t;
plot(t,y1,'k:','LineWidth',2);
y2 = 1+t+t.^2;
plot(t,y2,'k--','LineWidth',1.5);
y3 = 1+t+t.^2+t.^3;
plot(t,y3,'k-','LineWidth',3);
% edama_11_02: small number approximation for distances on a sphere
clear all;
RTOD = 180/pi;
DTOR = pi/180;
lat1 = 30*DTOR;
lon1 = 30*DTOR;
N = 101;
lat2min = 30*DTOR;
lat2max = 40*DTOR;
lat2 = lat2min + (lat2max-lat2min)*[1:N-1]'/(N-1);
lon2 = lat2;
latbar = (lat1 + lat2)/2;
Dlon = lon2 - lon1;
Dlat = lat2 - lat1;
cr = sin(lat1)*sin(lat2) + cos(lat1)*cos(lat2).*cos(Dlon);
r = acos(cr);
coslatbar2 = cos(latbar).^2;
Dlat2 = (Dlat).^2;
Dlon2 = (Dlon).^2;
t1 = Dlat2 + (Dlon).^2 .* coslatbar2;
r1 = sqrt( t1 );
t2 = Dlat2.*Dlon2.*(0.25-(coslatbar2/6));
r2 = sqrt( t1-t2 );
figure(1);
clf;
subplot(2,1,1);
set(gca,'LineWidth',2);
hold on;
axis( [30, 40, 0, 15 ] );
xlabel('latitude (deg)');
ylabel('r (deg)');
plot( RTOD*lat2, RTOD*r, '-', 'LineWidth', 5, 'Color', [0.8,0.8,0.8] );
plot( RTOD*lat2, RTOD*r1, 'k:', 'LineWidth', 2 );
plot( RTOD*lat2, RTOD*r2, 'k-', 'LineWidth', 2 );
subplot(2,1,2);
set(gca,'LineWidth',2);
hold on;
axis( [30, 40, 0, 0.02 ] );
xlabel('latitude (deg)');
ylabel('error (deg)');
plot( RTOD*lat2, RTOD*(r1-r), 'k:', 'LineWidth', 2 );
plot( RTOD*lat2, RTOD*(r2-r), 'k-', 'LineWidth', 2 );
% edama_11_03: variance using a small number approximation
% variance of T = 1/f = 2 pi / w computed by:
% small number approximation; and
% numerical simulation
clear all;
% abbreiation for 2 pi
tpi = 2*pi;
% mean and variance of Normally-disstributed variable w
wbar = 1.3*tpi;
sigmaw = 0.2;
% standard set-up for wb-axis used in histogram
Nbins = 201;
wbmin = 0;
wbmax = 20;
Dwb = (wbmax-wbmin)/(Nbins-1);
wb = wbmin + Dwb*[0:Nbins-1]';
% Normal pdf pw1(wb) with mean wbar and sqrt(variance) sigmaw
pw1 = (1/(sqrt(tpi)*sigmaw))*exp(-((wb-wbar).^2)/(2*(sigmaw^2)));
% realizations of a random variable
% with mean wbar and sqrt(variance) sigmaw
N = 10000; % number of realizations
wr = random('Normal',wbar,sigmaw,N,1);
% histogram of the realizatios
pw2 = hist(wr,wb)';
pw2 = pw2 / (Dwb*sum(pw2)); % normalize to empirial pdf pw2
% plot pw1 and pw2
clf;
subplot(2,1,1);
set(gca,'LineWidth',2);
hold on;
axis( [0, 10, 0, 3] );
xlabel('w (radians/s)');
ylabel('p(w)');
plot( wb, pw1, '-', 'LineWidth', 4, 'Color', [0.8,0.8,0.8] );
plot( wb, pw2, 'k:', 'LineWidth', 2 );
% standard setup of a Tb-axis used in histogram
Tbmin = 0;
Tbmax = 2;
DTb = (Tbmax-Tbmin)/(Nbins-1);
Tb = Tbmin + DTb*[0:Nbins-1]';
% mean and sqrt(variance) of T
Tbar = tpi/wbar;
sigmaT = sigmaw*tpi/(wbar^2);
% Normal pdf pT1(Tb) with mean Tbar and sqrt(variance) sigmaT
pT1 = (1/(sqrt(tpi)*sigmaT))*exp(-((Tb-Tbar).^2)/(2*(sigmaT^2)));
Tr = tpi./wr; % realizations of T computed from w
pT2 = hist(Tr,Tb)'; % histogram
pT2 = pT2 / (DTb*sum(pT2)); % normalize to empirical pdf pT2
% plot pT1 and pT2
subplot(2,1,2);
set(gca,'LineWidth',2);
hold on;
axis( [0, 2, 0, 25] );
xlabel('T (s)');
ylabel('p(T)');
plot( Tb, pT1, '-', 'LineWidth', 4, 'Color', [0.8,0.8,0.8] );
plot( Tb, pT2, 'k:', 'LineWidth', 2 );
% eda11_04: non-linear least squares with Newton's Method
% fit of a sinusoid of unknown (but approximately annual)
% frequency to Black Rock Forest temperature data
clear all;
% load data from file
D=load('../Data/brf_temp.txt');
% break out time t and temperature d data
t=D(:,1);
dobs=D(:,2);
Ns=size(D);
% clean the data by deleting outliers and zero (= no data) values
k = find( (dobs>-35) & (dobs<40) & (dobs~=0.0) );
t = D(k,1);
dobs = D(k,2);
N = length(dobs);
% some statistics
dmean = mean(dobs);
amp = 0.5*(max(dobs)-min(dobs));
% Annual period, frequency, angular frequency
Ta = 365;
fa = 1/Ta;
wa = 2*pi*fa;
% nonlinear theory and its derivative
% note model parameters are scaled to be order unity
% d = m1*dmean + amp*m2*sin(m4 wa t) + amp*m3*cos(m4 wa t)
% dd/dm1 = dmean
% dd/dm2 = amp*sin(m4 wa t)
% dd/dm3 = amp*cos(m4 wa t)
% dd/dm4 = amp*m2*wa*t*cos(m4 wa t) - amp*m3*wa*t*sin(m4 wa t)
% inital guesses of model
M=4;
m = zeros(M,1);
m(1) = 1; % average level is mean
m(2) = 0.01; % phase with max in summer
m(3) = -0.99; % so only -cos() term
m(4) = 1; % annual frequency
% data predicted by initial guess
s = amp*sin(m(4)*wa*t);
c = amp*cos(m(4)*wa*t);
dpre = dmean*m(1) + m(2)*s + m(3)*c;
% plot predicted data
figure(1);
clf;
set(gca,'LineWidth',2);
hold on;
axis([0,5000,-30,50]);
plot(t,dobs,'k.','LineWidth',2);
hold on;
plot(t,dpre,'-','LineWidth',2,'Color',[0.7,0.7,0.7]);
xlabel('time (days)');
ylabel('temperature (deg C)');
% plot observed data
figure(2);
clf;
set(gca,'LineWidth',2);
hold on;
axis([0,5000,-30,50]);
plot(t,dobs,'k.','LineWidth',2);
hold on;
xlabel('time (days)');
ylabel('temperature (deg C)');
% Newton's method
s = amp*sin(m(4)*wa*t);
c = amp*cos(m(4)*wa*t);
dpre = dmean*m(1) + m(2)*s + m(3)*c; % predicted data
dd = dobs - dpre; % deviation of data from current prediction
E = dd'*dd;
MAXITER = 10; % maximum number of iterations
for iter=[1:MAXITER] % interatively improve solution
Eprev = E;
% data kernel contains derivatives
G = [ dmean*ones(N,1), s, c, m(2)*wa*t.*c-m(3)*wa*t.*s ];
dm = (G'*G)\(G'*dd); % model perturbation dm by least sqiars
m = m + dm; % model update
s = amp*sin(m(4)*wa*t);
c = amp*cos(m(4)*wa*t);
dpre = m(1)*dmean + m(2)*s + m(3)*c; % predicted data
dd = dobs - dpre; % deviation of data from current prediction
E = dd'*dd; % total error
if( (abs(E-Eprev)/E) < 1.0e-6 )
break % terminate loop when errror fails to decrease
end
end
figure(1);
plot(t,dpre,'-','LineWidth',2,'Color',[0.9,0.9,0.9]);
figure(2);
plot(t,dpre,'-','LineWidth',2,'Color',[0.9,0.9,0.9]);
% undo scaling to arrive at final unknowns
ffinal = (wa*m(4))/(2*pi);
Tfinal = (2*pi)/(wa*m(4));
% plot final solution
figure(3);
clf;
set(gca,'LineWidth',2);
hold on;
axis([320,420,0,1]);
plot( [Tfinal, Tfinal]', [0, 0.75]', 'k-', 'LineWidth', 3 );
xlabel('period (days)');
ylabel('amplitude spectral density');
% covariance
% f = f0 + df and T = T0 + dT and T0 = 1 / f0
% T = 1 / (f0 + df ) = (1/f0) (1+df/f0)^-1
% = (1/f0) (1-df/f0) = (1/f0) - df/(f0^2)
% dT = (-1/(f0^2)) df
% var(dT) = (-1/(f0^2))^2 var(df)
% data kernel after final iteration
G = [ dmean*ones(N,1), s, c, m(2)*wa*t.*c-m(3)*wa*t.*s ];
sigma2d = E/(N-M); % posterior variance
GTGI = inv(G'*G); % unit covariance of m
sigma2m4 = sigma2d*GTGI(4,4); % variance of m(4)
sigma2f = (wa^2)*sigma2m4/((2*pi)^2); % variane of f
sigma2T = sigma2f / (ffinal^4); % variance of T
sigmaT = sqrt(sigma2T); % sqrt(variance) of T
DT95 = 2*sigmaT; % 95 percent confidence interval
plot( [Tfinal-DT95, Tfinal+DT95]', [0.7, 0.7]', 'k-', 'LineWidth', 3 );
fprintf('iterations %d\n', iter);
fprintf('posterior sqrt(sigma2d) %f\n', sqrt(sigma2d));
fprintf('Estimated period %.3f +/- %.3f days\n', Tfinal, DT95 );
% compare with result of FFT
% standard time & frequency definitions
Dt = t(2)-t(1);
tstart = t(1);
tend = t(N);
N2 = floor((tend-tstart)/Dt)+1;
N2 = 2*floor(N2/2); % insure even
t2 = tstart + Dt*[0:N2-1]';
fny = 1/(2*Dt);
Nf = N2/2+1;
Df = fny / (N2/2);
f = Df*[0:Nf]';
% interpolate data to even sampling
% but note that data gaps filled in with stright lines
d2 = interp1(t,dobs,t2,'linear','extrap');
% demean, take FFT and compute amplitude spectral density
d2 = d2-mean(d2);
d2tilda = Dt*fft(d2);
d2tilda = d2tilda(1:Nf);
asd = abs(d2tilda);
maxasd = max(asd);
plot( 1./f(2:Nf), 0.9*asd(2:Nf)/maxasd, 'k-', 'LineWidth', 2 );
plot( 1./f(2:Nf), 0.9*asd(2:Nf)/maxasd, 'ko', 'LineWidth', 2 );
% edama_11_05: non-linear least squares by gradient descent
% fit of a sinusoid of unknown (but approximately annual)
% frequency to Black Rock Forest temperature data
% note: Bill Menke patched and error in this code 09/20/2020
clear all;
% load Black Rock Forest temperature data
D=load('../Data/brf_temp.txt');
% break out time t and observed discharge dobs
t=D(:,1);
dobs=D(:,2);
Ns=size(D);
% clean the data by deleting outliers and zero (= no data) values
k = find( (dobs>-35) & (dobs<40) & (dobs~=0.0) );
t = D(k,1);
dobs = D(k,2);
N = length(dobs);
% some statistics
dmean = mean(dobs);
amp = 0.5*(max(dobs)-min(dobs));
% Annual period, frequency, angular frequency
Ta = 365;
fa = 1/Ta;
wa = 2*pi*fa;
% nonlinear theory and its derivative
% d = m1*dmean + m2*amp*sin(m4 wa t) + m3*amp*cos(m4 wa t)
% dd/dm1 = dmean
% dd/dm2 = amp*sin(m4 wa t)
% dd/dm3 = amp*cos(m4 wa t)
% dd/dm4 = amp*m2*wa*t*cos(m4 wa t) - amp*m3*wa*t*sin(m4 wa t)
% inital guesses of model parameters
M=4;
m = zeros(M,1);
m(1) = 1; % average level is mean
m(2) = 0; % phase with max in summer
m(3) = -0.75; % so only -cos() term
m(4) = 1; % fraction of annual frequency wa
Lm = sqrt(m'*m); % length of model
% prdited data based on starting model
s = amp*sin(m(4)*wa*t);
c = amp*cos(m(4)*wa*t);
dpre = dmean*m(1) + m(2)*s + m(3)*c;
figure(1);
clf;
set(gca,'LineWidth',2);
hold on;
axis([0,5000,-30,50]);
plot(t,dobs,'k.','LineWidth',2);
hold on;
% plot(t,dpre,'-','LineWidth',2,'Color',[0.7,0.7,0.7]);
xlabel('time (days)');
ylabel('temperature (deg C)');
% error and its gradient at the trial solution
mgo = m;
s = amp*sin(mgo(4)*wa*t);
c = amp*cos(mgo(4)*wa*t);
dprego = mgo(1)*dmean + mgo(2)*s + mgo(3)*c;
ego = dobs - dprego;
Ego = ego'*ego;
Estart = Ego;
Go = [ amp*ones(N,1), s, c, mgo(2)*wa*t.*c-mgo(3)*wa*t.*s ];
dEdmo = -2*Go'*ego;
% gradient descent parameters
alpha = 1.0;
c1 = 0.0001;
tau = 0.5;
MAXITER = 4000; % maximum numbe of iterations
for iter = [1:MAXITER] % interative improvement of solution
% downslope unit vector
v = -dEdmo / sqrt(dEdmo'*dEdmo);
% backstep
for kk=[1:100]
mg = mgo+alpha*v; % update the trail solution
s = amp*sin(mg(4)*wa*t);
c = amp*cos(mg(4)*wa*t);
dpreg = dmean*mg(1) + mg(2)*s + mg(3)*c; % predicted data
eg = dobs - dpreg; % individual errors
Eg = eg'*eg; % total error
% data kernel
G = [ dmean*ones(N,1), s, c, mg(2)*wa*t.*c-mg(3)*wa*t.*s ];
dEdm = -2*G'*eg; % gradient
% trminate backstep by Armijo's rule
c1 = 1.0e-4;
if( (Eg<=(Ego + c1*alpha*v'*dEdmo)) )
break;
end
% if the error went up, decrease the step size
alpha = tau*alpha;
end
% change in solution (with respect to initial guess)
Dmg = sqrt( (mg-mgo)'*(mg-mgo) );
% update
mgo=mg;
dprego = dpreg;
ego = eg;
Ego = Eg;
Go = G;
dEdmo = dEdm;
% terminate loop if solution hasn't changd much
if( (Dmg/Lm) < 1.0e-4 )
break;
end
end
m = mg; % final solution
dpre = dpreg; % final prediction of data
sigma2d = Ego/(N-M);
plot(t,dpre,'-','LineWidth',2,'Color',[0.9,0.9,0.9]);
% undo scalings
ffinal = wa*m(4)/(2*pi);
Tfinal = (2*pi)/(wa*m(4));
fprintf('Iterations: %d\n', iter );
fprintf('posterior sqrt(sigma2d) %f\n', sqrt(sigma2d));
fprintf('mean: %f amp: %f %f Dm: %f\n', m(1)*dmean, m(2)*amp, m(3)*amp, Dmg );
fprintf('Estimated period %.3f error ratio %.3f\n', Tfinal, Eg/Estart );
% edama_11_06: stochastic gradient method
% least squares fit of sinusoid of unknown frequencyk
% to Black Rock Forest temperature dataset
% The data are subsampled on most each iterations,
% % but full dataset is used ecvery 100 iterations.
% A fixed step size alpha is used to adcance the
% solution anti-parallel to the gradient.
% The solution is considered good enough
% when its posterior variance is reduced to a
% specififed value.
clear all;
% load Black Rock Forest temperature data
D=load('../Data/brf_temp.txt');
% break out time t and discharge dobs
t=D(:,1);
dobs=D(:,2);
Ns=size(D);
% clean the data by deleting outliers and zero (= no data) values
k = find( (dobs>-35) & (dobs<40) & (dobs~=0.0) );
t = D(k,1);
dobs = D(k,2);
N = length(dobs);
% subsampling information
decimation_factor = 1000;
Ns = floor(N/decimation_factor);
% some statistics
dmean = mean(dobs);
amp = 0.5*(max(dobs)-min(dobs));
% the solution considered good enough when
% its posterior variance is smaller than
low_enough_variance = 35.0;
% Annual period, frequency, angular frequency
Ta = 365.0;
fa = 1/Ta;
wa = 2*pi*fa;
% the model parameters m in this problem have been
% normalized so that we expect them to be of order
% 1. The DC offset is parameterized as a fraction
% of the mean of the data, dmean*m(1). The amplitudes
% of the sine and cosine are parameterized as a fraction
% of the overall amplitude of the data, m(2)*amp and
% m(3)*amp. The frequency is parameterized as a
% fraction of the annual frequecy, m(4)*wa.
% The resulting nonlinear theory and its derivative
% d = m1*dmean + m2*amp*sin(m4 wa t) + m3*amp*cos(m4 wa t)
% dd/dm1 = dmean
% dd/dm2 = amp*sin(m4 wa t)
% dd/dm3 = amp*cos(m4 wa t)
% dd/dm4 = amp*m2*wa*t*cos(m4 wa t) - amp*m3*wa*t*sin(m4 wa t)
% inital guesses of model parameters (trial solution)
M=4;
m = zeros(M,1);
m(1) = 1; % average level is close to the mean
m(2) = 0; % phase with max in summer, so no sine term
m(3) = -0.8; % cosine term is negative (hot in summer) and
% probably less than 1 sine data contains noise
m(4) = 1; % frequency close to annual frequency
% data predicted by trial solution
s = amp*sin(m(4)*wa*t);
c = amp*cos(m(4)*wa*t);
dpre = dmean*m(1) + m(2)*s + m(3)*c;
% plot observations
figure(1);
clf;
set(gca,'LineWidth',2);
hold on;
axis([0,5000,-30,50]);
plot(t,dobs,'k.','LineWidth',2);
hold on;
% plot(t,dpre,'-','LineWidth',2,'Color',[0.7,0.7,0.7]);
xlabel('time (days)');
ylabel('temperature (deg C)');
% trial solution
mg = m;
% the step size alpha is fixed
alpha = 1e-4;
% number of iterations is expected range of model
% parameters divided by step size
range = 1.0;
MAXITER = floor(1.5*range/alpha);
for iter = [1:MAXITER]
% subsample data
ns = unidrnd(N,Ns,1);
ts = t(ns);
dobss = dobs(ns);
% compute predicted data, the error and the
% gradient of the error of the subsampled data
s = amp*sin(mg(4)*wa*ts);
c = amp*cos(mg(4)*wa*ts);
dpreg = dmean*mg(1) + mg(2)*s + mg(3)*c;
eg = dobss - dpreg;
Eg = eg'*eg;
if( iter==1 )
Estart = Eg;
end
G = [ dmean*ones(Ns,1), s, c, m(2)*wa*ts.*c-m(3)*wa*ts.*s ];
dEdm = -2*G'*eg;
% posterior variance of data
sigma2d = Eg/(Ns-M);
% compute just the error for the full dataset
% but only every 100 iterations
if( mod(iter-1, 100) == 0 )
sfull = amp*sin(mg(4)*wa*t);
cfull = amp*cos(mg(4)*wa*t);
dpregfull = dmean*mg(1) + mg(2)*sfull + mg(3)*cfull;
egfull = dobs - dpregfull;
Egfull = egfull'*egfull;
% posterior variance of data
sigma2dfull = Egfull/(N-M);
if( iter==1 )
Estartfull = Egfull;
end
% stop when the veriance is low enough
if( sigma2dfull <= low_enough_variance )
break;
end
end
% downslope unit vector
v = -dEdm / sqrt(dEdm'*dEdm);
% update solution
mg = mg+alpha*v;
end
% final solution
m = mg;
s = amp*sin(m(4)*wa*t);
c = amp*cos(m(4)*wa*t);
% data predicted by final solution
dpre = dmean*m(1) + m(2)*s + m(3)*c;
plot(t,dpre,'-','LineWidth',2,'Color',[0.9,0.9,0.9]);
% undo scaling
ffinal = wa*m(4)/(2*pi);
Tfinal = (2*pi)/(wa*m(4));
fprintf('Iterations: %d\n', iter );
fprintf('mean: %f amp: %f %f\n', m(1)*dmean, m(2)*amp, m(3)*amp );
fprintf('Estimated period %.3f error ratio %.3f\n', Tfinal, Egfull/Estartfull );
% eda11_07: table lookup during a grid search
% example of a grid search that uses a table
% lookup of the function
% f(z) = sin(pi*sin(pi*z*sin(pi*(1-z^2))))
clear all;
% true values of the model parameters
mtrue = [0.5142, 0.7524]';
% the data are a function of time
N = 10001;
tmin = 0;
tmax = 1;
Dt = (tmax-tmin)/(N-1);
t = tmin+Dt*[0:N-1]';
% subsamples indices for plotting purposes
Nsub = 20;
isub = [1:Nsub:N]';
% the true data obey the formula
% dtrue = f(z) with z = m1 + m2*t
z = mtrue(1)+mtrue(2)*t;
dtrue = sin(pi*sin(pi*z.*sin(pi*(1-z.^2))));
% synthetic estimated data is true data plus random noise
sigmad = 0.1;
dobs = dtrue + random('Normal',0,sigmad,N,1);
etrue = dobs - dtrue;
Etrue = etrue'*etrue;
fprintf('true: m %f %f E %f\n', mtrue(1), mtrue(2), Etrue );
% plot the true and observed data
figure(1);
clf;
set(gca,'LineWidth',2);
hold on;
axis([tmin,tmax,-1.2,1.2]);
plot(t(isub),dobs(isub),'k.','LineWidth',2);
plot(t(isub),dtrue(isub),'-','LineWidth',7,'Color',[0.9,0.9,0.9]);
xlabel('t');
ylabel('d');
% range of m1 for grid search
M1 = 301;
M1=101;
m1min = 0;
m1max = 1;
Dm1 = (m1max-m1min)/(M1-1);
m1 = m1min+Dm1*[0:M1-1]';
% range of m2 for grid search
M2 = 301;
M2=101;
m2min = 0;
m2max = 1;
Dm2 = (m2max-m2min)/(M1-1);
m2 = m2min+Dm2*[0:M2-1]';
% standard grid search that computes f(z) exactly
% mote use of tic/toc pair to time exectiion
tic;
first = 1;
for m1trial = m1'
for m2trial = m2'
z = m1trial+m2trial*t;
dtrial = sin(pi*sin(pi*z.*sin(pi*(1-z.^2))));
etrial = dobs - dtrial;
Etrial = etrial'*etrial;
if( first==1 )
Ebest1 = Etrial;
m1best1 = m1trial;
m2best1 = m2trial;
first=0;
elseif (Etrial < Ebest1 )
Ebest1 = Etrial;
m1best1 = m1trial;
m2best1 = m2trial;
end
end
end
duration1=toc;
% print and plot results
fprintf('exact: m %f %f E %f time %f\n', m1best1, m2best1, Ebest1, duration1 );
z = m1best1+m2best1*t;
dest = sin(pi*sin(pi*z.*sin(pi*(1-z.^2))));
plot(t(isub),dest(isub),'-','LineWidth',4,'Color',[0.6,0.6,0.6]);
% grid search with table lookup of the funtion f(z)
tic;
% compute the table
K = 1001;
zmintable = -0.1;
zmaxtable = 2.1;
Dztable = (zmaxtable-zmintable)/(K-1);
ztable = zmintable+Dztable*[0:K-1]';
dtable=sin(pi*sin(pi*ztable.*sin(pi*(1-ztable.^2))));
% grid search with table lookup
first = 1;
for m1trial = m1'
for m2trial = m2'
z = m1trial+m2trial*t;
n = floor((z-zmintable)/Dztable)+1; % table index
dtrial = dtable(n); % table lookup
etrial = dobs - dtrial;
Etrial = etrial'*etrial;
if( first==1 )
Ebest2 = Etrial;
m1best2 = m1trial;
m2best2 = m2trial;
first=0;
elseif (Etrial < Ebest2 )
Ebest2 = Etrial;
m1best2 = m1trial;
m2best2 = m2trial;
end
end
end
duration2=toc;
% print and plot estimated solution
fprintf('lookup: m %f %f E %f time %f\n', m1best2, m2best2, Ebest2, duration2 );
fprintf('improvement in speed: %f\n', duration1/duration2 );
z = m1best2+m2best2*t;
dest = sin(pi*sin(pi*z.*sin(pi*(1-z.^2))));
plot(t(isub),dest(isub),'k-','LineWidth',2);
% edama_11_08: sigmoid function
% plot sigma(x) = (1+exp(-wx))^-1
% for a variety of w's. Note that
% dsigma/dx = w exp(-w) (1+exp(-wx))^-2
% so the slope at x=0 is s=w/4
clear all;
s = [100, 3, 1, 0.1]'; %slopes
Ns = length(s);
% standard x axis setup
Nx = 501;
xmin = -2;
xmax = 2;
x = xmin + (xmax-xmin)*[0:Nx-1]'/(Nx-1);
figure(1);
clf;
% loop over plots with different w's
for i = [1:Ns]
w=4*s(i);
sigma = (1+exp(-w*x)).^(-1);
subplot( Ns, 1, i );
set(gca,'LineWidth',2);
hold on;
axis([xmin,xmax,-1.5,1.5]);
plot(x,sigma,'k-','LineWidth',3);
xlabel('x');
ylabel('sigma(x)');
end
% eda11_09: neural net that approximates a step function
clear all;
% Step function definitions
Xc = 2; % center of the step function
h = 1; % height of the step function
slope = 200.0; % slope of the step function
% Neural Net definitions
Lmax = 3; % number of layers;
Nmax = 2; % maximum number of neurons in any layer
% ititialize the network. This step just
% creates zeroed arrays of the right size
[N,w,b,a] = eda_net_init(Lmax,Nmax);
% set up th network
N = [1, 1, 1]'; % one neuron in each layer
% weights
w(1,1,2) = 4*slope; % (neuron 1 on layer 2) fr (neuron 1 on previous layer)
w(1,1,3) = h; % (neuron 1 on layer 3) fr (neuron 1 on previous layer)
% biases
b(1,2) = -4*slope*Xc; % (neuron 1 on layer 2)
% plot this network
figure(1);
clf;
eda_net_view(N,w,b);
% standard setup of x-axis
Nx= 100;
xmin = 0;
xmax = 4;
x = xmin + (xmax-xmin)* [0:Nx-1]'/(Nx-1);
% evaluate the network for each x to get m(x)
m = zeros(Nx,1);
for kx = [1:Nx] % loop of x's
a(1:N(1),1) = [x(kx)]; % activity of (neuron 1 on layer 1) is x
a = eda_net_eval( N,w,b,a); % evaluate network
m(kx) = a(1:N(3),3); % activity of (neuron 1 on layer 3) is m(x)
end
% plot
figure(2);
clf;
subplot(2,1,1);
set(gca,'LineWidth',2);
hold on;
axis xy;
axis( [xmin xmax -1 2] );
plot(x,m,'k-','LineWidth',2);
title('step function; slope = 999');
xlabel('a_1_1');
ylabel('a_3_1(a_1_1)');
title(sprintf('step function; slope = %.0f',slope));
%----------------- do it again for slope = 1 ------------
% set up network
slope = 1.0;
w(1,1,2) = 4*slope;
w(1,1,3) = h;
b(1,2) = -4*slope*Xc;
% evaluate network for each x to met m2(x)
m2= zeros(Nx,1);
for kx = [1:Nx]
a(1:N(1),1) = [x(kx)];
a = eda_net_eval( N,w,b,a);
m2(kx) = a(1:N(3),3);
end
% plot
subplot(2,1,2);
set(gca,'LineWidth',2);
hold on;
axis xy;
axis( [xmin xmax -1 2] );
plot(x,m2,'k-','LineWidth',2);
title(sprintf('step function; slope = %.0f',slope));
xlabel('a_1_1');
ylabel('a_3_1(a_1_1)');
% edama_11_10: neural net approximating a 1D tower function
clear all;
% Tower function definitions
Xc = 2.0; % center of the tower function
Dx = 1.0; % width of the tower function
h = 1.5; % height of the tower function
slope = 200.0; % slope of the step function
% Neural Net definitions
Lmax = 3; % number of layer
Nmax = 2; % maximum number of neurons in any layer
% ititialize the network. This step just
% creates zeroed arrays of the right size
[N,w,b,a] = eda_net_init(Lmax,Nmax);
% implement a three-layer Neural net
N = [1, 2, 1]'; % number of neurons in eachc layer
% weights
w12 = 4*slope; % (neuron 1 on layer 2) fr (neuron 1 on previous layer)
w22 = 4*slope; % (neuron 1 on layer 2) fr (neuron 1 on previous layer)
w(1:N(2),1,2) = [w12;w22];
w13 = h; % (neuron 1 on layer 3) fr (neuron 1 on previous layer)
w23 = -h; % (neuron 1 on layer 3) fr (neuron 2 on previous layer)
w(1,1:N(2),3) = [w13;w23];
% biases
b12 = -w12 * ( Xc - Dx/2); % (neuron 1 on layer 2)
b22 = -w22 * ( Xc + Dx/2); % (neuron 2 on layer 2)
b(1:N(2),2) = [b12;b22];
% plot the net
figure(1);
clf;
eda_net_view(N,w,b);
% standard setup of x-axis
Nx= 101;
xmin = -1;
xmax = 5;
x = xmin + (xmax-xmin)* [0:Nx-1]'/(Nx-1);
% evaluate neural net for all x's to get m(x)
m = zeros(Nx,1);
for kx = [1:Nx] % loop over x
a(1:N(1),1) = [x(kx)]; % activity of (neuron 1 on layer 1) is x
a = eda_net_eval( N,w,b,a); % evaluate network
m(kx) = a(1:N(3),3); % m(x) is activity of (neuron 1 on layer 3)
end
% plot the result
figure(2);
clf;
subplot(2,1,1);
set(gca,'LineWidth',2);
hold on;
axis xy;
axis( [xmin xmax -1 2] );
plot(x,m,'k-','LineWidth',2);
title(sprintf('tower function; slope = %.0f',slope));
xlabel('a_1_1');
ylabel('a_3_1(a_1_1)');
%----------------- do it again for slope = 1 ------------
% set network parametres
slope = 1.0;
w12 = 4*slope;
w22 = 4*slope;
w(1:N(2),1,2) = [w12;w22];
b12 = -w12 * ( Xc - Dx/2);
b22 = -w22 * ( Xc + Dx/2);
b(1:N(2),2) = [b12;b22];
% evaluate for each x to get m2(x)
m2= zeros(Nx,1);
for kx = [1:Nx]
a(1:N(1),1) = [x(kx)];
a = eda_net_eval( N,w,b,a);
m2(kx) = a(1:N(3),3);
end
% plot
subplot(2,1,2);
set(gca,'LineWidth',2);
hold on;
axis xy;
axis( [xmin xmax -1 2] );
plot(x,m2,'k-','LineWidth',2);
title(sprintf('tower function; slope = %.0f',slope));
xlabel('a_1_1');
ylabel('a_3_1(a_1_1)');
% edama_11_11: neural net approximating a parabola y(x)
% by superimposing 1D towers
clear all;
% standard setup of an axis xc for centers of towers
Nc= 6; % number of centers in the approximation
xcmin = 0;
xcmax = 1;
xc = xcmin + (xcmax-xcmin)* [0:Nc-1]'/(Nc-1);
% yc: true value of the parabola
yc = 0.2 + 1.75*xc.^2;
% eight of the towers
h = yc;
% all the towers have the same width
Dx = ones(Nc,1)*(xcmax-xcmin)/(Nc-1);
% all the towers have the same slope
slope = 250;
% ---Neural Net definitions---
Lmax = 3; % number of layers;
Nmax = Nc * 2; % maximum number of neurons in any layer
% initialize network to zeroes
[N,w,b,a] = eda_net_init(Lmax,Nmax);
% create network of Nc 1D towers
[N,w,b] = eda_net_1dtower( Nc, slope, xc, Dx, h, N, w, b );
% plot
figure(1);
clf;
eda_net_view(N,w,b);
% ---Picking the points to plot---
% standard setup of an x-axis
Nx= 100;
xmin = xcmin - 0.1;
xmax = xcmax + 0.1;
x = xmin + (xmax-xmin)* [0:Nx-1]'/(Nx-1);
% evaluate netowrk at each x to get m(x)
m= zeros(Nx,1);
for kx = [1:Nx] % loop over x's
a(1:N(1),1) = [x(kx)]; % activity of (neuron 1 in layer 1) is x
a = eda_net_eval( N,w,b,a); % evaluate network
m(kx) = a(1:N(3),3); % m(x) is activity of (neuron 1 in layer 3)
end
% plot m(x)
figure(2);
clf;
subplot(2,1,1);
set(gca,'LineWidth',2);
hold on;
axis xy;
axis( [xmin xmax -1 2] );
plot(x,m,'k-','LineWidth',2);
title(sprintf('approximation with slope = %.2f',slope));
xlabel('a_1_1=x');
ylabel('y(x)=a_3_1(a_1_1)');
% true function
mtrue=zeros(Nx,1);
mtrue = 0.2 + 1.75*x.^2;
% plot true function
subplot(2,1,2);
set(gca,'LineWidth',2);
hold on;
axis xy;
axis( [xmin xmax -1 2] );
plot(x,mtrue,'k-','LineWidth',2);
title('true function');
xlabel('x');
ylabel('y(x)');
% edama_11_12: neural net approximating a 2D step function
clear all;
% this scripts plots an intermediate value z13 of the
% network when PLOTZ13=1. When PLOTZ13=0, the newtaork
% is evaluated using eda_net_eval(). However, when
% PLOTZ13=1, in order to access z13, the network is evaluated
% within this scipt.
PLOTZ13 = 1;
% ---Step function definitions---
Nc = 1; % the number of centers
Xc = [1.5]; % the center x-value of the step function
Yc = [1.5]; % the center y-value of the step function
Dx = [1]; % x-width of the step function
Dy = [1]; % y-width of the step function
h = [1]; % height of the step function
slope = 200; % slope of the step function
% ---Neural Net definitions---
Lmax = 4; % number of layers;
Nmax = 4*Nc; % maximum number of neurons in any layer
% initialize net to all zeroes
[N,w,b,a] = eda_net_init(Lmax,Nmax);
% initialize net to a seigle 2D tower
[N,w,b] = eda_net_2dtower( Nc, slope, Xc, Yc, Dx, Dy, h, N, w, b );
N = [2, Nc*4, Nc, 1]'; % reset N to correct number of neurons
% z: a(1:N(i),i) is a column-vector that gives the
% activitites of the N(i) neurons on the i-th layer
z = zeros( Nmax, Lmax );
% plot the network
figure(5);
clf;
eda_net_view(N,w,b);
% standard setup of x-axis and y-axis
Nx = 50;
Ny = 50;
xmin = 0;
ymin = 0;
xmax = 3.1;
ymax = 3.1;
x = xmin + (xmax-xmin)* [0:Nx-1]'/(Nx-1);
y = ymin + (ymax-ymin)* [0:Ny-1]'/(Ny-1);
if( PLOTZ13 )
% evaluate the network at each (x,y) to get m(x,y)
m= zeros(Nx,Ny);
% also save the input of to the last layer
zprev = zeros(Nx,Ny);
% ---Calculating the Activities---
for kx = [1:Nx]
for ky = [1:Ny]
a(1:N(1),1) = [x(kx);y(ky)];
for i=[2:Lmax]
z(1:N(i),i) = w(1:N(i),1:N(i-1),i) * a(1:N(i-1),i-1) + b(1:N(i),i);
if i < Lmax
a(1:N(i),i) = 1.0./(1.0+exp(-z(1:N(i),i)));
else
a(1:N(i),i) = z(1:N(i),i);
end
end
zprev(kx,ky) = z(1,3);
m(kx,ky) = a(1:N(3),4);
end
end
hold on;
eda_draw(zprev,'caption z(1,3)',' ',m,'caption a(1,4)');
else
% evaluate the network at each (x,y) to get m(x,y)
m= zeros(Nx,Ny);
% also save the input of to the last layer-
for kx = [1:Nx]
for ky = [1:Ny]
% activity of neuron 1 on layer 1) is x
% activity of neuron 2 on layer 1) is y
a(1:N(1),1) = [x(kx);x(ky)];
a = eda_net_eval( N,w,b,a); % evaluate network
m(kx,ky) = a(1:N(3),4); % f(x,y) is activity of (neuron 1 on layer 4)
end
end
hold on;
eda_draw(m,'caption a(1,4)');
end
% edama_11_13: least squares fit of single 2D tower to a Gaussian
% standard set-up of (x,y) axes
Nx = 21;
Ny = 21;
xmax = 5;
ymax = 5;
xmin = 0;
ymin = 0;
x = xmin + (xmax-xmin)* [0:Nx-1]'/(Nx-1);
y = ymin + (ymax-ymin)* [0:Ny-1]'/(Ny-1);
% (x,y) index arrays for easy access
kofij=zeros(Nx,Ny); % linear index k of a particular (i,j)
iofk=zeros(Nx*Ny,1); % row index of a particular k
jofk=zeros(Nx*Ny,1); % col index of a particular k
Nxy=0;
for i=[1:Nx]
for j=[1:Ny]
Nxy=Nxy+1;
kofij(i,j)=Nxy;
iofk(Nxy)=i;
jofk(Nxy)=j;
end
end
% Gaussian evaluated at (x,y) grid is true data array dobs
dobs = zeros(Nx,Ny); % data array
dobsvec = zeros(Nxy,1); % data array in vector form
amp = 5.0; % amplitude
xbar = 2.4; % mean in x
ybar = 2.6; % mean in y
sigmax = 1.2; % sqrt(variance) in x
sigmay = 0.8; % sqrt(variance) in y
cx = 1/(2*sigmax*sigmax); % shorthand
cy = 1/(2*sigmay*sigmay); % shorthand
% evauate Gaussian at (x,y) grid points
for i=[1:Nx]
for j=[1:Ny]
dobs(i,j) = amp*exp(-cx*((x(i)-xbar)^2) - cy*((y(j)-ybar)^2) );
dobsvec(kofij(i,j))=dobs(i,j); % map to vector
end
end
% build the test network for a single 2D tower
Nc = 1; % One tower
Xc = [2.5]; % Xc: the center of the tower
Yc = [2.5];
Dx = [2]; % Dx: the width of the tower
Dy = [2];
h = [1]; % h: the height of the tower
slope = 5/4; % slope: slope of the step function
% ---Neural Net definitions---
% Nmax: number of layers;
Lmax = 4;
% Lmax: maximum number of neurons in any layer
Nmax = Nc * 4;
% initialize neural net
[N,w,b,a] = eda_net_init(Lmax,Nmax);
% daLmaxdb(i,j,k): the change in the activity of the i-th neuron in layer Lmax
% due to a change in the bias of the of the j-th neuron in the layer k
daLmaxdb = zeros(Nmax,Nmax,Lmax); % initialize
% daLmaxdw(i,j,k,l): the change in the activity of the i-th neuron in layer Lmax
% due to a change in the weight connecting the j-th neuron in the layer l
% with the k-th neuron of layer l-1
daLmaxdw = zeros(Nmax,Nmax,Nmax,Lmax); % initialize
% create a net correspnding to a single 2D tower
[N,w,b] = eda_net_2dtower( Nc, slope, Xc, Yc, Dx, Dy, h, N, w, b );
figure(2);
clf;
eda_net_view(N,w,b);
% create a vector-to-net map that facilitates the updating
% of network biases and weights using linearized least squares
[ Nb,Nw,layer_of_bias,neuron_of_bias,index_of_bias,layer_of_weight, ...
r_neuron_of_weight,l_neuron_of_weight,index_of_weight ] = eda_net_map(N,w,b);
% linearized data kernel
M = Nb+Nw; % biases + weigthts
G = zeros(Nxy,M); % initialize
% predicted data
dpre = zeros(Nx,Ny); % initialize
dprevec = zeros(Nxy,1); % initialize
% initial predicted data
d0 = zeros(Nx,Ny); % initialize
d0vec = zeros(Nxy,1); % initialize
Niter=20;
for iter=[1:Niter] % iterative improvement using Newton's method
% loop over (x,y)'s of grid
for kx=[1:Nx]
for ky=[1:Ny]
% left activities are current (x,y)
a(1:N(1),1) = [x(kx);x(ky)];
% evaluate the network
a = eda_net_eval( N,w,b,a);
% appoximate Gaussian is the right activity
dpre(kx,ky) = a(1,4);
dprevec(kofij(kx,ky))=dpre(kx,ky);
% compute derivatives with respect to biases and weights
[daLmaxdw,daLmaxdb] = eda_net_deriv( N,w,b,a );
% first rows of G are derivatives with respect to biases
for ib=[1:Nb]
nob=neuron_of_bias(ib); % use net-vector mapping
lob=layer_of_bias(ib); % use net-bector mapping
% set row of G
G(kofij(kx,ky), ib ) = daLmaxdb(1, nob, lob );
end
% subsequent rows of G are derivatives with respect to weights
for iw=[1:Nw]
rnow = r_neuron_of_weight(iw);
lnow = l_neuron_of_weight(iw);
low = layer_of_weight(iw);
G(kofij(kx,ky), Nb+iw ) = daLmaxdw( 1, rnow, lnow, low );
end
end % next (x,y)
end
% deviation in data
Dd = dobsvec - dprevec;
if( iter==1 ) % save error of first iteration
d0 = dpre;
d0vec = dobsvec;
E0 = Dd'*Dd;
end
% implement some damping to improve stability
% but relax it after a few iterations
if( iter<5 )
epsi=1.0;
elseif (iter<10)
epsi=0.1;
else
epsi=0.001;
end
% damped least squares estimate of deviation in model
Dm = (G'*G+epsi*eye(M,M))\(G'*Dd);
% update biases
for ib=[1:Nb]
% use net-vector mapping
nob=neuron_of_bias(ib);
lob=layer_of_bias(ib);
b(nob, lob) = b(nob, lob) + Dm(ib);
end
% update weights
for iw=[1:Nw]
% use net-vector mapping
rnow = r_neuron_of_weight(iw);
lnow = l_neuron_of_weight(iw);
low = layer_of_weight(iw);
w(rnow, lnow, low) = w(rnow, lnow, low) + Dm( Nb+iw );
end
end % next iteration
% evaluate the Gaussian by looping over (x,y) grid
for kx=[1:Nx]
for ky=[1:Ny]
% left activities are current (x,y)
a(1:N(1),1) = [x(kx);x(ky)];
% evaluate net
a = eda_net_eval( N,w,b,a);
% approximate Gaussian is right activity
dpre(kx,ky) = a(1,4); % final prediction
dprevec(kofij(kx,ky))=dpre(kx,ky); % vector veresion
end
end
eda_draw( dobs, 'caption d-obs', d0, 'caption d-start',dpre, 'caption d-pre');
% final error
Dd = dobsvec - dprevec;
error = Dd'*Dd;
fprintf('starting error %.2f\n', E0 );
fprintf('ending error %.2f\n', error );
% edama_11_14 neural net approximating the linear function y=cx+d
% a network that implements the linear function y=cx+d
% for small w112 and w113 so that the approximation
% is reasonably accurate in the range -1<x<+1
% (when c, d, are both or order 1)
c = 0.5; % true slope
d = 2.0; % true intercept
% standard set-up of x axis
Nx = 101;
xmin = -1;
xmax = 1;
x = xmin + (xmax-xmin)*[0:Nx-1]'/(Nx-1);
% the function we are trying o approximate
ytrue = c*x+d;
% ---Neural Net definitions---
% Nmax: number of layers;
Lmax = 4;
% Lmax: maximum number of neurons in any layer
Nmax = 1;
[N,w,b,a] = eda_net_init(Lmax,Nmax);
N = [1, 1, 1, 1]';
% these two slopes are arbitrary, as long as they are small
slope1 = 0.01;
slope2 = 0.02;
w(1,1,2) = 4*slope1;
w(1,1,3) = 4*slope2;
w2w3 = w(1,1,2)*w(1,1,3);
b(1,2) = 0;
% the details:
% near x=0, sigma(wx) = (w/4)*x + (1/2)
% z2 = w2*x + b2 = w2*x + 0
% a2 = sigma(z2) = w2/4+(1/2)
% z3 = w3*a2 + b3 = w3*(w2/4+(1/2)+b3 = (w2*w3/4)*x + w3/2 + b3
% so with choice b3 = -w3/2
% z3 = (w2*w3/4)*x
% a3 = sigma(z3) = (w2*w3/16)*x + (1/2)
% z4 = w4*a3+b4
% so with choices w4=16c/(w2*w3) and b4=d-8c(w2*w3)
% a4 = c*x + d
% set weights and biases
w(1,1,4) = 16*c/w2w3;
b(1,2) = 0;
b(1,3) = -w(1,1,3)/2;
b(1,4) = d-8*c/w2w3;
figure(1);
clf;
eda_net_view(N,w,b);
% evaluate linear function for each x
yest = zeros(Nx,1);
for kx=[1:Nx] % loop over points on x-axis
% left activity is current x
a(1,1) = x(kx);
% evaluate network
a = eda_net_eval( N,w,b,a);
% estimated function y(x) is top activity
yest(kx)=a(1,4);
end
% plot the results
figure(2);
clf;
set(gca,'LineWidth',2);
hold on;
axis( [xmin, xmax, d+c*xmin, d+c*xmax] );
plot( x, ytrue, '-', 'LineWidth', 7, 'Color', [0.8,0.8,0.8]);
plot( x, yest, 'k-', 'LineWidth', 2);
xlabel('x');
ylabel('y(x)');
% edama_11_15: network to match non-linear response function
% the data d(t) (equal) F[x(t)] (about equal) h(t) (convolved with) x(t)
% the process starts with network that implements the linear relation
% d(t) (about equal) h(t) (convolved with) x(t)
% and then uses iterative least squares to improve the network
% so that it implements d(t) (equal) F[x(t)]
% least squares parameters
epsi1 = 1; % initial damping
epsi2 = 0.05; % final damping
Niter=100; % number of iterations
% load data from file
D = load('../Data/nonlinearresponse.txt');
% break out into separate variables
t = D(:,1); % time in days
% discharge in cubic meters per second
dobs = D(:,2)/35.3146;
% precipitation in millimeters per day
x = D(:,3)*25.4;
% in this example, we use only a Nx day portion of the
% dataset, mainly so that when we plot it, storm events
% are clearly visible
Nx=501;
istart = 1;
t = t(istart:istart+Nx-1);
x = x(istart:istart+Nx-1);
dobs = dobs(istart:istart+Nx-1);
tmin = min(t);
tmax = max(t);
xnorm = 1/max(x);
dobsnorm = 1/max(dobs);
x = x * xnorm;
dobs = dobs * dobsnorm;
% The parameters for the filter h(t) have been chosen
% by eye to lead to a reasonable prediction of the data
Nc = 10;
c0 = 2;
To = 3;
f = exp(-[0:Nc-1]'/To);
f = [f(1)/1000;f(1:end-1)];
h = c0*f/sum(f);
% now implement this filter with a neural net
% ---Neural Net definitions---
% Nmax: number of layers;
Lmax = 4;
% Lmax: maximum number of neurons in any layer
Nmax = Nc ;
[N,w,b,a] = eda_net_init(Lmax,Nmax); % initialize the net
% linear neural net that implements the filter h
slope1 = 0.01;
slope2 = 0.02;
[N,w,b] = eda_net_linear( Nc, slope1, slope2, 0.0, h, N, w, b );
% Phase 1: Improve only biases and weights in the linear filter
figure(1);
clf;
eda_net_view(N,w,b);
% the neural net will now be adjusted via non-linear least squares
% to improve its prediction. We first create a vector-net mapping
% of the net's biases and weights to facilitate the process
[ Nb,Nw,layer_of_bias,neuron_of_bias,index_of_bias,layer_of_weight, ...
r_neuron_of_weight,l_neuron_of_weight,index_of_weight ] = eda_net_map(N,w,b);
% linearized data kernel
M = Nb+Nw; % improve both biases and weights
G = zeros(Nx,M); % initialize
% predicted data
dpre = zeros(Nx,1); % initialize
% initial predicted data
d0 = zeros(Nx,1);
for iter=[1:Niter] % iterative improvemy via Newton's Method
% loop over x's
for kx=[1:Nx]
% loop of windows of data upon which filter h acts
for j=[1:N(1)]
k=kx-j+1;
% left activities are window of x's
% except if window does beyong available x's
if( k>0 )
a(j,1) = x(k);
else
a(j,1) = 0;
end
end
a = eda_net_eval( N,w,b,a); % evaluate net
% predicted data is right activity
dpre(kx,1) = a(1,4);
% derivatives with respect to weigts and biases
[daLmaxdw,daLmaxdb] = eda_net_deriv( N,w,b,a );
% build linearized data kernel G
% first rows of G have derivatives with respect to biases
for ib=[1:Nb] % loop over biases
% use vector-net mapping
nob=neuron_of_bias(ib);
lob=layer_of_bias(ib);
G(kx, ib ) = daLmaxdb(1, nob, lob );
end
% subsequent rows of G have derivatives with respect to weights
for iw=[1:Nw] % loop over weights
% use vector-net mapping
rnow = r_neuron_of_weight(iw);
lnow = l_neuron_of_weight(iw);
low = layer_of_weight(iw);
% use vector-net mapping
G(kx, Nb+iw ) = daLmaxdw( 1, rnow, lnow, low );
end
end % next x
% deviation in data
Dd = dobs - dpre;
if( iter==1 )
d0 = dpre; % save prediction
E0 = Dd'*Dd; % save error
end
% implement some damping to improve stability
% but relax it after a few iterations
if( iter==50 )
epsi1=epsi1 / 10;
end
Dm = (G'*G+epsi1*eye(M,M))\(G'*Dd);
% add in bias perturbations
for ib=[1:Nb]
nob=neuron_of_bias(ib);
lob=layer_of_bias(ib);
b(nob, lob) = b(nob, lob) + Dm(ib);
end
% add in weights perturbations
for iw=[1:Nw]
rnow = r_neuron_of_weight(iw);
lnow = l_neuron_of_weight(iw);
low = layer_of_weight(iw);
w(rnow, lnow, low) = w(rnow, lnow, low) + Dm( Nb+iw );
end
end % next iteration
for kx=[1:Nx]
% left activities are window of x's
% except if window does beyong available x's
for j=[1:N(1)] % loop over x's
k=kx-j+1;
if( k>0 )
a(j,1) = x(k);
else
a(j,1) = 0;
end
end
a = eda_net_eval( N,w,b,a); % evaluate net
% right activity is predicted data
dpre(kx) = a(1,4);
end
Dd = dobs - dpre;
Ei = Dd'*Dd;
% the discharge and precipitation data
figure(2);
clf;
subplot(2,1,1);
hold on;
set(gca,'LineWidth',3);
axis( [tmin, tmax, -1, 1 ] );
plot( t, x, 'k-', 'LineWidth', 3 );
plot( [tmin,tmax]', [0,0]', 'k-', 'LineWidth', 2 );
xlabel('time (days)');
ylabel('precipitation (mm/day)');
subplot(2,1,2);
set(gca,'LineWidth',3);
hold on;
axis( [tmin, tmax, -1, 1] );
plot( t, dobs, '-', 'LineWidth', 5, 'Color', [0.8,0.8,0.8] );
plot( [tmin,tmax]', [0,0]', 'k-', 'LineWidth', 2 );
xlabel('time (days)');
ylabel('discharge (m3/s)');
plot( t, d0, 'k--', 'LineWidth', 2 );
fprintf('starting error %.2f\n', E0 );
figure(3);
clf;
subplot(2,1,1);
hold on;
set(gca,'LineWidth',3);
axis( [tmin, tmax, -1, 1 ] );
plot( t, x, 'k-', 'LineWidth', 3 );
plot( [tmin,tmax]', [0,0]', 'k-', 'LineWidth', 2 );
xlabel('time (days)');
ylabel('precipitation (mm/day)');
subplot(2,1,2);
set(gca,'LineWidth',3);
hold on;
axis( [tmin, tmax, -1, 1] );
plot( t, dobs, '-', 'LineWidth', 5, 'Color', [0.8,0.8,0.8] );
plot( [tmin,tmax]', [0,0]', 'k-', 'LineWidth', 2 );
xlabel('time (days)');
ylabel('discharge (m3/s)');
plot( t, dpre, 'k--', 'LineWidth', 2 );
fprintf('intermediate error %.2f with %d unknowns\n', Ei, Nb+Nw );
% Phase 2: Add links in network that make it more complicated than
% a linear filter, iitially setting these weights to small random
% numbers, and the improve this network
% add links
for i=[1]
for j=[1:N(1)]
if( w(i,j,2) == 0 )
w(i,j,2) = random('Normal',0,0.0001,1,1);
end
end
end
for i=[1]
for j=[1:N(1)]
if( w(i,j,3) == 0 )
w(i,j,3) = random('Normal',0,0.0001,1,1);
end
end
end
figure(4);
clf;
eda_net_view(N,w,b);
% the neural net will now be adjusted via non-linear least squares
% to improve its prediction. We first create a vector-net mapping
% of the net's biases and weights to facilitate the process
[ Nb,Nw,layer_of_bias,neuron_of_bias,index_of_bias,layer_of_weight, ...
r_neuron_of_weight,l_neuron_of_weight,index_of_weight ] = eda_net_map(N,w,b);
% linearized data kernel
M = Nb+Nw;
G = zeros(Nx,M); % initialize
% predicted data
dpre = zeros(Nx,1); % initialize
% initial predicted data
d0 = zeros(Nx,1); % initialize
for iter=[1:Niter] % iteratively improve via Newton's Method
% loop over x's
for kx=[1:Nx]
% loop of windows of data upon which filter h acts
for j=[1:N(1)]
k=kx-j+1;
% left activities are window of x's
% except if window does beyong available x's
if( k>0 )
a(j,1) = x(k);
else
a(j,1) = 0;
end
end
a = eda_net_eval( N,w,b,a);
% right activity is predicted data
dpre(kx,1) = a(1,4);
% compute derivatives with respect to weights and biases
[daLmaxdw,daLmaxdb] = eda_net_deriv( N,w,b,a );
% build data kernel
% first rows are derivatives with respect to biases
for ib=[1:Nb]
% use vector-net mapping
nob=neuron_of_bias(ib);
lob=layer_of_bias(ib);
G(kx, ib ) = daLmaxdb(1, nob, lob );
end
% subsequent rows are derivatives with respect to weights
for iw=[1:Nw]
% use vector-net mapping
rnow = r_neuron_of_weight(iw);
lnow = l_neuron_of_weight(iw);
low = layer_of_weight(iw);
G(kx, Nb+iw ) = daLmaxdw( 1, rnow, lnow, low );
end
end % next x
Dd = dobs - dpre;
% implement some damping to improve stability
% but relax it after a few iterations
if( iter==50 )
epsi2=epsi2/10;
end
% damped least squares
Dm = (G'*G+epsi2*eye(M,M))\(G'*Dd);
% add in bias perturbations
for ib=[1:Nb]
nob=neuron_of_bias(ib);
lob=layer_of_bias(ib);
b(nob, lob) = b(nob, lob) + Dm(ib);
end
% add in weights perturbations
for iw=[1:Nw]
rnow = r_neuron_of_weight(iw);
lnow = l_neuron_of_weight(iw);
low = layer_of_weight(iw);
w(rnow, lnow, low) = w(rnow, lnow, low) + Dm( Nb+iw );
end
end % next iteration
% loop over (x,y)
for kx=[1:Nx]
% activities
for j=[1:N(1)]
k=kx-j+1;
if( k>0 )
a(j,1) = x(k);
else
a(j,1) = 0;
end
end
a = eda_net_eval( N,w,b,a);
dpre(kx) = a(1,4);
end
Dd = dobs - dpre;
error = Dd'*Dd;
fprintf('final error %.2f with %d unknowns\n', error, Nb+Nw );
% break out into separate variables, metric conversion
% time in days
t = D(:,1);
% discharge in cubic meters per second
dobs = D(:,2)/35.3146;
% precipitation in millimeters per day
x = D(:,3)*25.4;
Nx=501;
istart = 499;
t = t(istart:istart+Nx-1);
x = x(istart:istart+Nx-1);
dobs = dobs(istart:istart+Nx-1);
tmin = min(t);
tmax = max(t);
x = x * xnorm;
dobs = dobs * dobsnorm;
% one last evaluation of predicted data
for kx=[1:Nx]
% activities
for j=[1:N(1)]
k=kx-j+1;
if( k>0 )
a(j,1) = x(k);
else
a(j,1) = 0;
end
end
a = eda_net_eval( N,w,b,a);
dpre(kx) = a(1,4);
end
figure(5);
clf;
subplot(2,1,1);
hold on;
set(gca,'LineWidth',3);
axis( [tmin, tmax, -1, 1 ] );
plot( t, x, 'k-', 'LineWidth', 3 );
plot( [tmin,tmax]', [0,0]', 'k-', 'LineWidth', 2 );
xlabel('time (days)');
ylabel('precipitation (mm/day)');
subplot(2,1,2);
set(gca,'LineWidth',3);
hold on;
axis( [tmin, tmax, -1, 1] );
plot( t, dobs, '-', 'LineWidth', 5, 'Color', [0.8,0.8,0.8] );
plot( [tmin,tmax]', [0,0]', 'k-', 'LineWidth', 2 );
xlabel('time (days)');
ylabel('discharge (m3/s)');
plot( t, dpre, 'k--', 'LineWidth', 2 );
% edama_11_16, make time sereis x(t) and d(t)
% used by edama_11_15, by solving the
% nonlinear differential equation in the text
% by Euler's method
clear all;
% Note: in order to prevent the overwriting of the
% '../Data/nonlinearresponse.txt file of the standard
% distribution, this script writes the file
% '../Data/mynonlinearresponse.txt, instead.
% stardard setup of a time axis, t
N=1001;
Dt = 1;
t = Dt*[0:N-1]';
% main variables
x = zeros(N,1); % precipitation
y = zeros(N,1); % watershed volume
d = zeros(N,1); % discharge
% precipitation events are randomly placed spikes
Np = floor(N/10);
p = random('unid',N,Np,1); % uniformly distributed times
cexp = 0.5;
x(p) = random('exp',cexp,Np,1); % exponentially-distributed amplitudes
% plot precipitation, x(t)
figure(1);
clf;
subplot(3,1,1);
set(gca,'LineWidth',2);
hold on;
axis( [0, N, 0, 2 ] );
plot(t,x,'k-','LineWidth',2);
xlabel('time,t');
ylabel('precipitation, x(t)');
% quick and dirty integration of the differential
% equation by Euler's method. Since y(t) is only
% exemplary, there is no special reason for the
% solution to be especially accurate
cf = 0.2;
for i=[2:N] % loop over times
d(i) = cf*(y(i-1)^2);
dydt = x(i) - d(i);
y(i) = y(i-1) + dydt * Dt;
if( y(i) < 0 ) % negative volumes not allowed
y(i)=0;
end
end
% plot discharge
subplot(3,1,2);
set(gca,'LineWidth',2);
hold on;
axis( [0, N, 0, 1.1*max(d) ] );
plot(t,d,'k-','LineWidth',2);
xlabel('time,t');
ylabel('discharge, d(t)');
% plot volume
subplot(3,1,3);
set(gca,'LineWidth',2);
hold on;
axis( [0, N, 0, 1.1*max(y) ] );
plot(t,y,'k-','LineWidth',2);
xlabel('time,t');
ylabel('volume, y(t)');
dlmwrite('../Data/mynonlinearresponse.txt', [t,d,x], 'delimiter', '\t' );