Page 89 - MATLAB Recipes for Earth Sciences

4.10 Curvilinear Regression 81

plot(meters,age,'+',meters,p_age,'g-',...
meters,p_age + 2 * delta,'r', meters,p_age - 2 * delta,'r')
grid on

We now use another synthetic data set that we generate using a quadratic
relationship between the barium content (in wt.%) down a sediment core
(in meters).

meters = 20 * rand(30,1);
barium = 1.6 * meters.^2 - 1.1 * meters + 1.2;
barium = barium + 40.* randn(length(meters),1);

plot(meters,barium,'o')
bariumcont = [meters barium];

bariumcont = sortrows(bariumcont,1);
save bariumcont.txt bariumcont -ascii
The synthetic bivariate data set can be loaded from fi le bariumcont.txt.

bariumcont = load('bariumcont.txt');

meters = bariumcont(:,1);
barium = bariumcont(:,2);
plot(meters,barium,'o')

Fitting a polynomial of degree n=2 yields a convincing regression result
compared to the linear model.

p = polyfit(meters,barium,2)
p =
1.8272 -4.8390 -1.4428
As shown above, the true values for the three coefficients are +1.6, –1.1 and

+1.2. There are some discrepancies between the true values and the coeffi -
cients estimated using polyfit. The regression curve and the error bounds
can be plotted by typing (Fig. 4.8)

plot(meters,barium,'o'), hold
plot(meters,polyval(p,meters),'r')
[p,s] = polyfit(meters,barium,2);
[p_barium,delta] = polyval(p,meters,s);

plot(meters,barium,'+',meters,p_barium,'g',meters,...
p_barium+2*delta,'r--',meters,p_barium-2*delta,'r--')
grid on
xlabel('meters'), ylabel('barium content')