Live Script edama_13
This Live Script supports Chapter 13 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.
% edapy_13_01: exemplary errors related to redefining variables
fprintf('the variable pi: %.4f (before)\n', pi);
pi=2;
fprintf('the variable pi: %.4f (after)\n', pi);
zeros = 3;
x = zeros(3,1);
% edapy_13_02: time arithmetic
% the clear undos the destruction of the last cell
clear all
fprintf('seconds DT between 2008/2/11 4:4:1 and 2008/2/11 4:4:0\n');
DT = 86400*(datenum(2008,2,11,4,4,1)- ...
datenum(2008,2,11,4,4,0));
fprintf('DT: %.2f\n',DT);
DT = 86400*(datenum(1975,1,1,0,0,1)-datenum(1974,12,31,23,59,59));
fprintf('leap second checker\n');
fprintf('seconds DT between 1975/1/1 0:0:1 and 1974/12/31 23:59:59\n');
fprintf('right answer is 3, calculation gives DT: %.2f\n',DT);
% edapy_13_03: reading idiosynchratic text file
clear all;
% since size of file unknaon, allocate very big arrays
NBIG=1000000;
t = zeros(NBIG,1);
d = zeros(NBIG,1);
% names of months
moname = {'Jan','Feb','Mar','Apr','May','Jun','Jul','Aug','Sep','Oct','Nov','Dec'};
% tab
tab = sprintf('\t');
% open file
filename='../Data/brf_raw.txt';
fid = fopen(filename,"r");
N=0; % line number counter
while 1 % loop over lines in file
myline = fgetl(fid); % read line
if ~ischar(myline) % check for end of file
fclose(fid);
break;
end
N=N+1; % increent line counter
% now process the line
myline = replace(myline,tab," "); % tab to space
Ns = length(myline);
while 1 % replace double spaces with a single space
myline = replace(myline," "," ");
Nsnew = length(myline);
if( Nsnew==Ns)
break;
end
end
s = strsplit(myline); % split into parts, result is celstr
da = str2num(char(s(2))); % decode day
yr = str2num(char(s(4))); % decode year
j= ismember( moname, s(3) ); % look-up month name in table
mo = find(j==1,1); % month in 1-12 range
s2 = strsplit(char(s(1)),"-"); % split time into start, end
s2a = char(s2(1)); % start piece
s2b = char(s2(2)); % end piece
hr1 = str2num(s2a(1:2)); % start hour
hr2 = str2num(s2b(1:2)); % end hour
mn1 = str2num(s2a(3:4)); % start minute
mn2 = str2num(s2b(3:4)); % end minute
sc1 = 0; % start second
sc2 = 0; % end second
% detect end time of 2400 and turn it into 23:00:00
if( (hr2==24) && (mn2==0) )
hr2 = 23;
mn2 = 59;
sc2 = 59;
end
if(N==1) % refernce time is start of year of first sampple
tref = datenum( yr, 1, 1, 0, 0, 0 );
end
t1=datenum(yr,mo,da,hr1,mn1,sc1); % start datenum
t2=datenum(yr,mo,da,hr2,mn2,sc2); % end datenum
t(N)=0.5*(t2+t1) - tref; % average time - reference time
va = str2num(char(s(5))); % decode data value
d(N) = va; % data value
end
% truncate to actual size
t = t(1:N);
d = d(1:N);
figure(1);
clf;
set(gca,'LineWidth',2);
set(gca,'FontSize',14);
hold on;
plot(t,d,'k-','LineWidth', 2);
xlabel('time (days)');
ylabel('temperature (degC)');
% edama_13_04, example of use of eda_draw() function
clear all;
% define an examplary matrix, G
N=40;
G=[1:N]'*[1:N];
% deinine a exemplacy vector, m
m = [1:N]';
% multiply
d = G*m;
% draw
eda_draw(d, 'caption d', '=', G, 'caption G', m, 'caption m');
% additional text
fprintf('Fig. The equation Gm=d, depicted using eda_dra\n');
% edama_13_05: simple examples of a function
% see function definitions in the last dell of this LiveScript
% function called as follows
radius=2;
area = areaofcircle(radius);
fprintf('circle:\n');
fprintf('radius=%.4f area=%.4f\n',radius,area);
fprintf('\n');
% function called as follows
a=2;
b=4;
[circ, area] = CandAofrectangle(a,b);
fprintf('rectangle:\n');
fprintf('length=%.4f width=%.4f circumference=%.4f area=%.4f\n',a,b,circ,area);
fprintf('\n');
% eda13_06: examples of reorganizing a matrix
% define a matrix
A=[ [1,2,3,4]; [5,6,7,8]; [9,10,11,12] ];
[N,M] = size(A);
fprintf('a test matrix A\n');
A
% note that a matrix can be referred to using only a single
% index, in which case, it acts like a vector constructed
% by concatenating the columns of the matrix, one on top
% of the other
fprintf('A(4) is A(1,2)\n');
A(4)
% create a column-vector, c, containing the elements of A,
% organized column-wise.
c = A(:);
fprintf('transpose of a vector formed from A by method 1\n');
c'
% here's another way to do it
% create a column-vector, c, containing the elements of A,
% organized column-wise.
c = reshape(A,M*N,1);
fprintf('transpose of a vector formed from A by method 2\n');
c'
% row vector
r = reshape(A,1,M*N);
fprintf('row vector formed from A\n');
c'
% the matrix reformed from the vector
B = reshape(c,N,M);
fprintf('the vector reformed from the vector\n');
B
% when you 'unravel' an array with elements A[i,j] into
% column vector with elements b[k,0], you sometimes need
% to determine the correspondence beteen (i,j) and k.
% and vice versa:
i= 2;
j= 1;
fprintf('for a (%d,%d) matrix M[i,j] unraveled into a vector c[k]\n',N,M);
% corresponding index into the array b
k = sub2ind(size(A),i,j);
fprintf('(i,j) of (%d,%d) corresponds to k of %d\n', i,j,k );
fprintf('A(i,j)=%.0f and c[k]=%.0f', A(i,j),b(k) );
[ii,jj] = ind2sub(size(A),k)
fprintf('k of %d corresponds to (i,j) of (%d,%d) \n', k, i,j );
% edama_12_07: function Phi used in Lagrange multiplier figure
% standard set-up of (x,y) axes
N=51;
M=51;
Dx=0.01;
Dy=0.01;
x=Dx*[0:N-1]';
y=Dy*[0:M-1]';
% a long-wavelngth, 2D sinusoid
Phi=sin(x)*sin(y)';
eda_draw(' ', Phi);
% edama_12_08: Supports figure 13.3
% covariance correspondig to second derivative
% set constants
s = 0.1; % scale factor
g2 = 1.0; % amplitude
% standard set-up of x-axis
N2 = 100;
N = 2*N2+1;
Dx = 1.0;
x = Dx*[-N2:N2]';
xmin = min(x);
xmax = max(x);
% autocovariance, with |x| = |xi - xj|
c = g2 * (1.0+s*abs(x)) .* exp(-s*abs(x));
% Gaussian autocovariance for comparison
cg = g2*exp(-0.5*s*s*(x.^2));
% plot c(x) vs. x
figure(1);
clf;
set(gca,'LineWidth', 2);
set(gca,'FontSize',14);
hold on;
axis([xmin, xmax, -0.1, 1.1])
plot(x,c,'k-','LineWidth',3);
plot(x,cg,'k:','LineWidth',2);
plot(x,0.0*c,'k-','LineWidth',1);
xlabel('x');
ylabel('c(x)');
% edama_12_09, computing ln det Cm
% model parameter axis, t, with sampling Dt
M=64; % number of model parameters
Dt=1.0; % time sanpling
tmax=Dt*(M-1); % maximum time
t=Dt*[0:M-1]'; % time
fprintf('M: %d\n', M);
% covariance matrix, of cos(w0t) exp(-st) form
n=6; % number of oscillations on (0,tmax) interal
w0 = n*2*pi/tmax; % angular frequency of oscilation
tmid=tmax/2.0;
s = Dt/60.0; % decay rate
gamma = 1.0;
gamma2 = gamma^2; % variance
% build covariances via outer product
Tm = abs(t*ones(1,M)-ones(M,1)*t');
Cm = gamma2*cos(w0*Tm).*exp(-s*Tm);
eda_draw('caption Cm',Cm);
% Brute force method
D = det(Cm);
logdetCm_A = log(D);
fprintf('Brute force method, log det Cm = %e\n', logdetCm_A);
% Eigenvalue method
LAMBDA = eig(Cm);
logdetCm_B = sum(log(LAMBDA));
fprintf('Eigenvalue method, log det Cm = %e\n', logdetCm_B);
% Cholesky method
R = chol(Cm);
logdetCm_C = 2.0*sum(log(diag(R)));
fprintf('Cholesky method, log det Cm = %e\n', logdetCm_C);
fprintf('\n');
fprintf("brute force works fine for M=%d, but fails for M=230 and above\n",M);
function a = areaofcircle(r)
% computes area, a, of circle of radius, r.
a = pi * (r^2);
return
end
function [c,a] = CandAofrectangle(l, w)
% computes circumference, c, and area, a, of
% a rectangle of length, l, and width, w.
c = 2*(l+w);
a = l*w;
return
end