%% Physics 233, MATLAB Exercises for Hughes and Hase: Chapter 5
% MATLAB activities to complement _Measurements and their Uncertainties_,
% Chapter 5, "Data visualization and reduction." Section numbers in the
% book correspond to those listed below.
%
% Author: J. S. Dodge, 2016

%% 5.1 Producing a good graph
%% 5.1.1 The independent and dependent variables
% Consider the data tabulated in Prob. (5.3), taken from an experiment to
% verify Ohm's law.

current = (10:10:90);   % (uA)
voltage = [0.98, 1.98, 2.98, 3.97, 4.95, 5.95, 6.93, 7.93, 8.91];   % (mV)

%%
% The problem states that the voltmeter precision is 0.01 mV, and that the
% uncertainty in the current is negligible. This suggests we use the
% current as the independent variable. Plot data points using markers.
% MATLAB provides several choices; see the |plot| documentation. Use black
% circles (|ko|) here. 

plot(current, voltage, 'ko');
xlabel('Current (\muA)')
ylabel('Voltage (mV)')

%% 5.1.2 Linearizing the data
% Linearizing a graph can be a very effective way to demonstrate the
% validity of certain functional relationships between two quantities.
% Consider first the example given in the text, the period $T$ of a
% pendulum as a function of its length $L$. Use the random number generator
% to simulate period measurements with a constant standard error $\alpha_T
% = 0.1$ s, for pendulum lengths ranging from 10 cm to 100 cm. Use |T| for
% the model prediction, and |Tm| for the simulated measurements. 

L = (0.1:0.1:1.0);  % m
g = 9.8;            % m/s^2
T = 2*pi*sqrt(L/g); % s

rand('seed',0)
%rng('default')
alphaT = 0.03;
Tm = T + alphaT*randn(size(T));

plot(L,Tm,'ko');
xlabel('Length (m)')
ylabel('Period (s)')

%%
% Now plot |Tm.^2| versus |L| to linearize the graph. One advantage of a
% linear plot is that you can assess its quality visually just by
% positioning a straight edge through the data---try this with both of
% these plots, and notice that the original plot _looks_ fairly linear, but
% that a linear fit doesn't pass through the origin. This is obviously
% wrong: even if you didn't know the functional relationship between the
% period and the pendulum length, the period clearly must be zero when the
% length is zero. In contrast, the plot of |Tm.^2| versus |L| is consistent
% with a linear fit that passes through the origin.

plot(L,Tm.^2,'ko')
xlabel('Length (m)')
ylabel('(Period)^2 (s^2)')

%%
% Yet another way to linearize the same data is to use a
% <http://www.mathworks.com/help/matlab/ref/loglog.html |loglog|> plot.
% These are useful for identifying a _power-law_ relationship between
% two quantities, i.e.,
%
% $$ y = Ax^\alpha, $$
%
% where $A$ and $\alpha$ are unknown. Taking the logarithm of both sides
% yields a _linear_ relatonship between $\log y$ and $\log x$:
%
% $$ \log y = log A + \alpha\log x. $$
%
% Use |loglog| to demonstrate this with the |L|, |Tm| data. Show the data
% as black circles (|ko|), and add a red line (|r-|) to show the model
% behavior, |T|, and a blue dashed line (|b--|) to show |sqrt(L)|, to
% demonstrate the advantage of using |loglog| to identify a power-law
% functional dependence. The data fall approximately on a straight line,
% with a _slope_ similar to the plot of |sqrt(L)| versus |L|.

loglog(L, Tm, 'ko', L, T, 'r-', L, sqrt(L), 'b--');
xlabel('Length (m)')
ylabel('Period (s)')

%%
% Yet another common linearization method uses the
% <http://www.mathworks.com/help/matlab/ref/semilogx.html |semilogx|> and
% <http://www.mathworks.com/help/matlab/ref/semilogy.html |semilogy|>
% functions, which is useful to demonstrate an _exponential_ functional
% dependence. Consider, for example, the voltage across a 47 nF capacitor
% as it discharges from an initial value of 1 V across a 100 kOhm resistor.
% Add noise with a 10 microvolt amplitude using the random number
% generator.

x = linspace(0,100);    % s
R = 100e3;              % Ohm
C = 47e-9;              % Farad
tau = R*C*1000;         % ms
V0 = 1;                 % V
alphaV = 10e-6;         % V

V = V0*exp(-x/tau);
Vm = V + alphaV*randn(size(V));

%%
% First plot the data using linear scales.

plot(x,Vm,'ko',x,V,'r-')
xlabel('Time (ms)')
ylabel('Voltage (V)')

%%
% Now plot it on a semilog scale.

semilogy(x,Vm,'ko',x,V,'r-')
xlabel('Time (ms)')
ylabel('Voltage (V)')

%%
% The linear plot shows little variation for |t| greater than 20 ms, while
% the semilog plot allows this data to be visualized much more effectively.
% The model function appears as a straight line on the |semilogy| plot,
% while the data follows a straight line for |t| less than about 50 ms. As
% |t| increases above 50 ms, the noise begins to dominate and the voltage
% fluctuates randomly with a 10 microvolt width about zero. About half of
% these fluctuations have negative values that must be excluded from the
% plot because their logarithm is undefined; MATLAB issues a |warning:
% negative data ignored| statement to let you know. Most of the remaining
% positive values are greater than 1 microvolt: with 10 microvolt Gaussian
% noise centered at zero, it is more than ten times more likely that the
% voltage will be greater than one microvolt than between one microvolt and
% zero.


%% 5.1.3 Appropriate scales for the axes
%%
% MATLAB plotting functions compute the axis limits automatically from the
% data range, but the user can override these choices. In the
% current-voltage plot above, for example, we can use
% <http://www.mathworks.com/help/matlab/ref/ylim.html |ylim|> to readjust
% the vertical scale so that the data lies closer to the diagonal. This
% allows the viewer to identify the proportionality relationship more
% easily, and it is also more appealing visually.

plot(current, voltage, 'ko');
xlabel('Current (\muA)')
ylabel('Voltage (mV)')

ylim([0 10])

%%
% We can also rescale the the pendulum period to include the origin; this
% makes it clear that |T| has a nonlinear relationship with |L|. The |ylim|
% function requires a two-element vector input, |[ymin, ymax]|. To set
% |ymin| without changing |ymax|, get the limits by calling |ylim| without
% any argument, changing the first element, then calling |ylim| again with
% the modified vector as an argument.

plot(L,Tm,'ko');
xlabel('Length (m)')
ylabel('Period (s)')

yl = ylim;
disp(yl)

yl(1) = 0;
ylim(yl)

%%
% We can also do this using
% <http://www.mathworks.com/help/matlab/graphics-objects.html graphics
% objects>, which were implemented in MATLAB in version 2014b.

plot(L,Tm,'ko')
xlabel('Length (m)')
ylabel('Period (s)')

Ax = gca;
disp(Ax)
Ax.YLim(1) = 0;

%%
% The command |gca| stands for _get current axes_, and it returns the
% handle to the current _axes object_. These objects have standardized
% <http://www.mathworks.com/help/matlab/ref/axes-properties.html
% _properties_>, like |YLim|, that can be accessed by appending the
% property name after the object name with a dot (|.|) in between. Notice
% that two of these axes properties are |XScale| and |YScale|, and they are
% both set to |linear|. We can also set them both to |log| to reproduce the
% behavior of the |loglog| plot function.

Ax.XScale = 'log';
Ax.YScale = 'log';

%%
% Axes are only one of several types of graphics objects, which are
% organized hierarchically, as shown
% <http://www.mathworks.com/help/matlab/graphics-objects.html here>. The
% <http://www.mathworks.com/help/matlab/ref/groot.html graphics root
% object> usually refers to your computer monitor, and the |groot| function
% provides access to its
% <http://www.mathworks.com/help/matlab/ref/root-properties.html
% properties>. All _figure objects_ are _children_ of |groot|, which in
% turn is the _parent_ object to the figures. The function
% <http://www.mathworks.com/help/matlab/ref/gcf.html |gcf|> returns the
% handle of the current figure object, which provides access to its
% <http://www.mathworks.com/help/matlab/ref/figure-properties.html
% properties>.  As we've already seen, the function
% <http://www.mathworks.com/help/matlab/ref/gca.html |gca|> returns the
% handle of the current _axes object_, which provides access to its
% <http://www.mathworks.com/help/matlab/ref/axes-properties.html
% properties>. Below the axes object level, MATLAB defines a variety of
% _chart objects_ that are designed for specific types of plots. At the
% lowest level of the hierarchy, _primitive objects_ include lines, text,
% and other elements of each plot. 
% If no figure or axes objects yet exist, the |plot|
% function creates a _figure object_, which contains an _axes_ object,
% which in turn contains the _line_ objects that comprise the plot. Take
% some time to familiarize yourself with the properties of each of these
% object types; see
% <http://www.mathworks.com/help/matlab/graphics-object-properties.html
% here> for a complete list.
%
% Return now to the exponential decay plot and examine its components. Note
% that |plot| returns an array of two _line_ objects, one for the data and
% one for the model. Use the appropriate index for |Plt| to access these
% properties individually.

Root = groot;
disp('Root:')
disp(Root);

Fig = figure;
disp('Fig:')
disp(Fig)

Ax = axes;
disp('Ax:')
disp(Ax)

Plt = plot(x,Vm,'ko',x,V,'r-');
XLab = xlabel('Time (ms)');
YLab = ylabel('Voltage (V)');

disp('Plt:')
disp(Plt)

disp('Plt(1):')
disp(Plt(1))

disp('Plt(2):')
disp(Plt(2))

disp('XLab:')
disp(XLab)

disp('YLab:')
disp(YLab)

%%
% Now we can modify the plot by changing the graphics object properties
% directly. For example, we can convert the linear plot to a semilog plot
% by changing the |YScale| axes property; we also can change the appearance
% of the text, markers, and lines.

Ax.YScale = 'log';

Ax.FontSize = 18;
Plt(1).Marker = 's';
Plt(1).MarkerSize = 6;
Plt(1).Color = 'b';

Plt(2).LineWidth = 2;
Plt(2).LineStyle = '--';
Plt(2).Color = 'k';

%% 5.1.4 Labelling the axes
% Your turn: make a plot of the current through an ideal diode as a
% function of applied voltage.

V = (0:0.05:0.6);   % V

kB = 1.3806e-23;    % J/K
qe = 1.6022e-19;    % C
Tr = 296;           % Room temperature, in degrees Kelvin

VT = kB*Tr/qe;      % Thermal voltage, in Volts
Is = 1e-12;         % Saturation current, in Amperes

I = Is*(exp(V/VT) -1);

%% 5.1.5 Adding data points and error bars to graphs
% To include errorbars, use the
% <http://www.mathworks.com/help/matlab/ref/errorbar.html |errorbar|>
% function. Returning to the pendulum example, the error vector |errT| must
% be the same size as the primary data vectors |L| and |Tm|, so we use
% <http://www.mathworks.com/help/matlab/ref/ones.html |ones|> to make an
% array of ones that has the same size as |T|, and multiply this vector by
% the constant uncertainty |alphaT|.

errT = alphaT*ones(size(Tm));

errorbar(L,Tm,errT,'ko');
xlabel('Length (m)')
ylabel('Period (s)')

%%
% Now plot |Tm.^2| versus |L| to linearize the graph. Remember that we also
% need to propagate the uncertainty in |Tm| to produce the correct
% uncertainty in |Tm.^2|; we'll use the calculus-based method of Sec.
% 4.1.2.

TmSq = Tm.^2;
alphaTmSq = 2*alphaT.*Tm;

errorbar(L,TmSq,alphaTmSq,'ko')
xlabel('Length (m)')
ylabel('(Period)^2 (s^2)')

%% 5.1.6 Adding a fit or trend line
% MATLAB provides a variety of ways to fit a model to data. One is the
% Curve Fitting Tool, a graphical user interface (GUI) designed for
% interactive fitting. It relies on the
% <http://www.mathworks.com/help/curvefit/fit.html |fit|> function of the
% <http://www.mathworks.com/help/curvefit/index.html Curve Fitting
% Toolbox>. Use this to fit the power-law model $y = a x^b$ to the pendulum
% data. See the
% <http://www.mathworks.com/help/curvefit/list-of-library-models-for-curve-and-surface-fitting.html
% documentation> for a list of models that are already included in the
% Curve Fitting Toolbox. You can also define your own models; see the
% documentation for
% <http://www.mathworks.com/help/curvefit/fit.html#examples examples>.
%
% The |fit| function requires a minimum of three inputs: the independent
% variable, |x|, the dependent variable |y|, and the model |fitType|. The
% independent and dependent variables must be column vectors, but |L| and
% |Tm| are currently row vectors; we can use the colon operator (|:|) to
% convert the data to column vector format.

disp('L is a row vector:')
disp(L)
L = L(:);
disp('Now L is a column vector:')
disp(L)

Tm = Tm(:);

[fitModel,gof,output] = fit(L,Tm,'power1');

disp('fitModel:')
disp(fitModel)

disp('gof:')
disp(gof)

disp('output:')
disp(output)

%%
% Review the
% <http://www.mathworks.com/help/curvefit/fit.html#zmw57dd0e32281
% documentation> for detailed descriptions of each of these output
% variables. Right now, we just want to plot the fit with the data. The
% first output of |fit| is |fitModel|, which is a _cfit object_ that
% contains the model and the best-fit parameters, among other things. The
% syntax |fitModel(L)| treats the best-fit model as a function on the
% independent variable |L|.

Lmodel = linspace(1e-3,1);
Lmodel = Lmodel(:);
Tmodel = fitModel(Lmodel);

errorbar(L,Tm,errT,'ko');

hold on
plot(Lmodel,Tmodel,'r-')
hold off

xlabel('Length (m)')
ylabel('Period (s)')

%% 5.1.7 Adding a title or caption
% In your notebook, use titles for figures. In a report, _never_ use titles
% for figures. That's what the caption is for. To put a title on a figure
% for your notebook, use the |title| function. The conventions around
% legends are less strict: they are usually useful for notebook figures,
% and it is usually acceptable, if not recommended, in reports.
%
% Try to choose an informative title; a common error is to just restate the
% axis labels, like 'Period versus pendulum length.'

title('Pendulum period is proportional to sqrt(L)')
legend('Simulated measurements','Model: T = A*sqrt(L)',...
    'Location','SouthEast')

%% 5.2 Using a graph to see trends in the data
% Reproduce Fig. 5.5. The text gives the fit parameters, and we can
% estimate the noise level from the tail in Fig. 5.5(b).

tau(1) = 8.67;
tau(2) = 67;
a(1) = 6.5e7;
a(2) = 0.5e7;

alphaN = 1e5;

alphaTau(1) = 0.04;
alphaTau(2) = 4;

nPoints = 200;
x = linspace(0, 60, nPoints);
x = x(:);

n = a(1)*exp(-x/tau(1)) + a(2)*exp(-x/tau(2));
nm = n + alphaN*randn(size(n));

figure
plot(x,nm/1e6,'ko','MarkerFaceColor','k','MarkerSize',6)
xlabel('Time (seconds)')
ylabel('Number of atoms/10^6')
box off
text(0.95,0.95,'(a)',...
    'Units','Normalized')

figure
semilogy(x,nm/1e6,'ko','MarkerFaceColor','k','MarkerSize',6)
xlabel('Time (seconds)')
ylabel('Number of atoms/10^6')
box off
text(0.95,0.95,'(b)',...
    'Units','Normalized')

%%
% To show the trend lines in Fig. 5.5(c), fit separately to the data for
% |t < 20| and |t > 40|. Use
% <http://www.mathworks.com/help/matlab/math/matrix-indexing.html#bq7egb6-1
% logical indexing> to select data in these two ranges. To see how logical
% indexing works, consider the following expressions.

seq = (1:10);
disp('seq')
disp(seq)

disp('seq > 5')
disp(seq > 5)

disp('seq < 5')
disp(seq < 5)

disp('seq(seq > 5)')
disp(seq(seq > 5))

disp('seq(seq < 5)')
disp(seq(seq < 5))

%%
% The expression |seq > 5| evaluates the inequality for each
% element of |seq| to assign a |true (1)| or |false (0)| value to the
% corresponding element in a boolean array of the same size. The expression
% |seq(seq > 5)|, then, returns only those elements of |seq| for which |seq
% > 5| is |true|.
%
% Now use it this to fit to |t < 20| and |t > 40|.

initialDecay = fit(x(x < 20),log(nm(x < 20)/1e6),'poly1');
laterDecay = fit(x(x > 40),log(nm(x > 40)/1e6),'poly1');

figure
plot(x,log(nm/1e6),'ko','MarkerSize',6)
hold on
plot(initialDecay)
plot(laterDecay)
hold off
legend off
ylim([0.5 4.5])
xlabel('Time (seconds)')
ylabel('ln(No. of atoms/10^6)')
box off
text(0.95,0.95,'(c)',...
    'Units','Normalized')

%% 5.2.1 Adding a linear trend line
% Determine the best linear fit to the current-voltage data of Sec. 5.1.1,
% using Eqs. (5.1) and (5.2) to get the intercept and slope; then, get 
% parameter uncertainties from Eqs. (5.3) and (5.4), and relate them to 
% the confidence bounds that MATLAB gives.
%
% The current (in µA) is the independent variable, which corresponds to $x$
% in Eqs. (5.1) and (5.2). The voltage (in mV) is the dependent variable,
% which corresponds to $y$. Use the colon operator to put the data in
% column-vector format, so we can use it with the |fit| function.

current = current(:);
voltage = voltage(:);

N = length(current);
Delta = N*sum(current.^2) - (sum(current))^2;

cFit = (sum(current.^2)*sum(voltage) - ...
    sum(current)*sum(current.*voltage))/Delta;

mFit = (N*sum(current.*voltage) - sum(current)*sum(voltage))/Delta;

alphaCU = sqrt(sum((voltage - mFit*current - cFit).^2)/(N-2));

alphac = alphaCU*sqrt(sum(current.^2)/Delta);
alpham = alphaCU*sqrt(N/Delta);

fprintf('voltage = mFit*current + cFit, fit with Eqs. (5.1) and (5.2)\n')
fprintf('Best-fit parameters +/- one-sigma uncertainty:\n')
fprintf('mFit:\t %.4g +/- %.4g\n',mFit,alpham)
fprintf('cFit:\t %.4g +/- %.4g\n',cFit,alphac)

%%
% Compare this with the result from the |fit| function. Plot the fit and
% show the residuals separately.

[cvFit,gof,output] = fit(current,voltage,'poly1');
disp(cvFit)

figure;
subplot(3,1,[1 2])
plot(cvFit,current,voltage)
legend('Location','SouthEast')
xlabel('')
set(gca,'XTickLabel',{})
ylabel('Voltage (mV)')
subplot(3,1,3)
stem(current,voltage-cvFit(current))
xlabel('Current (\muA)')
ylabel('Residuals (mV)')

%%
% The |cfit| object displays the parameter uncertainties in terms of 95%
% "confidence bounds" instead of the one-standard-deviation uncertainty
% values $\alpha_c$ and $\alpha_m$. To examine the relationship between
% them, use the <http://www.mathworks.com/help/curvefit/coeffvalues.html
% |coeffvalues|> and <http://www.mathworks.com/help/curvefit/confint.html
% |confint|> functions to get the best-fit parameters and their confidence
% bounds. Next, find the distance between the best-fit parameter estimates
% and their confidence bounds, in units of $\alpha_m$ and $\alpha_c$.

betaFit = coeffvalues(cvFit);
ciFit = confint(cvFit);

disp('betaFit:')
disp(betaFit)
disp('ciFit:')
disp(ciFit)

ciErrm = (ciFit(:,1)-betaFit(1))/alpham;
ciErrc = (ciFit(:,2)-betaFit(2))/alphac;

disp('ciErrm:')
disp(ciErrm)
disp('ciErrc:')
disp(ciErrc)

%%
% After reviewing the results in Table 3.1, you might have expected the 95%
% confidence bound on $p = \{m, c\}$ to be $\hat{p} \pm 1.96\alpha_p$, not
% $\hat{p} \pm 2.36\alpha_p$ as shown above (we use $\hat{p}$ to denote the
% best-fit value of $p$). The difference appears because Table 3.1 assumes
% that the standard deviation $\sigma$ _is known independently_, whereas
% our values for $\alpha_m$ and $\alpha_c$ are _derived from the fit_.
%
% To clarify this distinction, let $\hat{p}(X)$ be the best-fit estimate of
% the parameter $p$, given a set of data $X$. It's important to note here
% that $\hat{p}(X)$ is a random variable, since it is a function of data
% that contains random fluctuations. Now define $L(X)$ to be the lower
% confidence bound and $U(X)$ the upper confidence bound on $p$, where
% $L(X)$ and $U(X)$ are equidistant from $\hat{p}(X)$ and enclose 95% of
% the probability distribution for $p$, which we may either know
% independently or infer from $X$:
% 
% $$ P(L(X) < p < U(X)) = 0.95 $$
%
% Here also, $L(X)$ and $U(X)$ are _functions of the data_, so we should
% expect them to fluctuate randomly from one measurement to the next.
% If we know _a priori_ that $\hat{p}(X)$ will be normally distributed
% with a standard deviation $\sigma_p$, then the 95% confidence interval
% $(L(X),U(X))$ is just the usual $(\hat{p}(X) - 1.96\sigma_p, \hat{p}(X) +
% 1.96\sigma_p)$ from Table 3.1, since the fluctuations in $\hat{p}(X)$
% don't change the shape of the inferred normal distribution. If we _also_
% have to estimate $\sigma_p$ from the data, however, then the uncertainty
% in the width _will_ change the inferred distribution, since we
% effectively need to consider a _range_ of normal distributions, all with
% different widths. The resulting distribution is known as the _Student's t
% distribution_, and its shape depends on how many _degrees of freedom_ are
% available to determine the width: in the current-voltage fit, we have
% nine data points to estimate two parameters, so the number of degrees of
% freedom is seven (9-2 = 7). As the number of degrees of freedom
% increases, our ability to estimate the distribution width $\sigma$
% improves, and the Student's t distribution approaches the normal
% distribution with $\mu = 1, \sigma = 1$. Use the
% <http://www.mathworks.com/help/stats/tpdf.html |tpdf|> function to plot
% some examples, and the <http://www.mathworks.com/help/stats/tinv.html
% |tinv|> function to compute the 95% confidence bounds on $p$ in terms of
% the standard error $\alpha_p$. Use the
% <http://www.mathworks.com/help/matlab/creating_plots/specify-legend-descriptions-during-line-creation.html
% |DisplayName|> to assign a legend label to each distribution function as
% it's plotted, then use
% <http://www.mathworks.com/help/matlab/ref/legend.html |legend('show')|>
% to show them all at the end.

p = linspace(-5,5);
figure;
plot(p,normpdf(p),...
    'DisplayName','Normal')
xlabel('p')
ylabel('Probability distribution function')
hold on

alpha = 0.05;
fprintf(...
    'For alpha = %.2f, the confidence level is %.2f%%, and \n', ...
    alpha, (1- alpha)*100)
fprintf(...
    'the confidence bounds are set at P(t(L,DOF) < %.3f) and P(t(U,DOF) < %.3f).\n', ...
    alpha/2, 1- alpha/2)

fprintf('\n')
fprintf('DOF\tL\tU\n')
dof = [1 3 7 10 100];
for iDOF = dof
    fprintf(...
        '%d\t%.4f\t%.4f\n',...
        iDOF, tinv(alpha/2,iDOF), tinv(1-alpha/2,iDOF))
    plot(p,tpdf(p,iDOF),...
        'DisplayName',sprintf('t, %d DOF', iDOF))
end

hold off
legend('show')
%%
% Let's summarize. When we fit a linear model to data, the best-fit
% parameters are given by Eqs. (5.1) and (5.2). They are functions of the
% data, so we expect them to fluctuate from one set of equivalent
% measurements to another. Eqs. (5.3) and (5.4) give estimates for the
% standard deviation for these fluctuations, assuming normally distributed
% noise in the $y$ measurements with a standard deviation $\alpha_{CU}$
% (Hughes and Hase refer to this as the *common uncertainty*). We may know
% this independently, but if we don't, then Eq. (5.6) provides an estimate
% _under the assumption that the fit limited by statistical uncertainty_.
% (As you know by know, this assumption is often violated in real
% experiments!) Under this assumption, the confidence interval for $p =
% \{m, c\}$ is $(\hat{p} - t\alpha_p, \hat{p} + t\alpha_p)$, where the
% prefactor $t$ is a function of the confidence level (95% in our example)
% and the number of degrees of freedom that we used to determine the
% |voltage| uncertainty (9-2=7 in our example). If we know $\alpha_{CU}$
% independently, then the confidence bounds are the usual ones for a normal
% distribution with standard deviation $\alpha_p$. When we estimate
% $\alpha_{CU}$ from the fit, the confidence bounds will be wider.

%%
% The |fit| function has an optional output |gof| that includes the common
% uncertainty (in the field |rmse|) and the number of degrees of freedom
% for the fit (in the field |dfe|). The field |sse| is the sum squared
% error that is minimized in the fit. The other fields, |rsquare| and
% |adjrsquare|, describe how much of the overall variation in the data can
% be accounted for with the model---this can be important in biological or
% social sciences, where the experimenter may not have any control over the
% variables, but it's usually not of much interest in experimental physics,
% where we can purposely vary the independent parameter to see how the
% dependent parameter depends on it.

disp('Goodness of fit output structure:')
disp(gof)

sseAlt = sum((voltage-cvFit(current)).^2);
numParameters = length(coeffvalues(cvFit));
numPoints = length(voltage);
dfeAlt =  numPoints - numParameters;
rmseAlt = sqrt(sseAlt/dfeAlt);

fprintf('Explicitly calculated goodness of fit parameters:\n')
fprintf('Sum-squared error:\t%.4e\n',sseAlt)
fprintf('Degrees of freedom:\t%d\n',dfeAlt)
fprintf('Root-mean-square error:\t%.4f\n',rmseAlt)

%%
% Hughes and Hase mention that Eqs. (5.1) and (5.2) follow from the _method
% of least squares_ and refer the reader to references for details. We'll
% review the derivation here, using an elegant matrix notation that can be
% implemented easily in MATLAB.
%
% Mathematically, fitting a line to data involves solving an overdetermined
% <http://www.mathworks.com/help/matlab/math/systems-of-linear-equations.html
% system of linear equations>. To see why, write the model equation $y = mx
% + b$ for each pair $\{x_k, y_k\}$ in the data;
%
% $$ y_1 = m x_1 + b $$
%
% $$ y_2 = m x_2 + b $$
%
% $$ \qquad \vdots$$
%
% $$ y_N = m x_N + b $$
%
% This is a system of $N$ linear equations in two unknowns, which will be
% overdetermined for case when more than two of these equations are
% linearly independent.
%
% Now, rewrite this in matrix form: let $\mathbf{X}$ be a column vector of
% the independent variable data $\{x_k, k = 1,...,N\}$, let $\mathbf{Y}$ be
% a column vector of the dependent variable data $\{y_k, k = 1,...,N\}$,
% and identify the intercept and slope of the linear model with elements of
% a _parameter vector_ $\mathbf{\beta} = [c\ m]^\mathrm{T}$. Then
%
% $$ \mathbf{Y} = \mathbf{A\beta}, $$
%
% where $\mathbf{A} \equiv [\mathbf{X}^0\ \mathbf{X}^1]$ is the |N x 2|
% <https://en.wikipedia.org/wiki/Vandermonde_matrix Vandermonde matrix> of
% $\mathbf{X}$; the exponential notation $\mathbf{X}^n$ applies the
% exponential operation to each element in the vector, so $\mathbf{X}^0$ is
% just a |N x 1| column vector of ones.
%
% Apply this to our current-voltage example, with |c = 0.03| and |m =0.1|.

X = current;
Y = voltage;
c = 0.03;
m = 0.1;
beta0 = [c; m];

A = [X.^0 X.^1];
disp('A:')
disp(A)
disp('beta0:')
disp(beta0)
disp('A*beta0:')
disp(A*beta0)

%%
% When we have $N = 2$ independent equations, then the matrix $\mathbf{A}$
% is square and invertible, and we have $\mathbf{\beta} =
% \mathbf{A}^{-1}\mathbf{Y}$.

Alim = A(1:2,:);
Ylim = Y(1:2,:);

disp('Alim:')
disp(Alim)
disp('Ylim:')
disp(Ylim)
disp('inv(Alim):')
disp(inv(Alim))
disp('inv(Alim)*Ylim:')
disp(inv(Alim)*Ylim)

%%
% When we have $N > 2$ independent equations, then $\mathbf{A}$ is no
% longer square, and in general the matrix equation $\mathbf{Y} =
% \mathbf{A\beta}$ does not have a solution; in geometric terms, the image
% of $\mathbf{\beta}$ in $\mathbf{A}$ may not include $\mathbf{Y}$. But we
% can still find an _approximate_ solution that minimizes the distance
% between $\mathbf{Y}$ and $\mathbf{A\beta}$, or $||\mathbf{Y} -
% \mathbf{A}\mathbf{\beta}||$. We can then define the
% <https://en.wikipedia.org/wiki/Moore?Penrose_pseudoinverse
% _pseudoinverse_> of $\mathbf{A}$ to be the matrix $\mathbf{A}^+$ that
% transforms $\mathbf{Y}$ into that solution: $\mathbf{\beta} =
% \mathbf{A}^+\mathbf{Y}$. In MATLAB, the function
% <http://www.mathworks.com/help/matlab/ref/pinv.html |pinv|> performs the
% pseudoinverse operation.

betaPinv = pinv(A)*Y;

disp('Pseudoinverse, pinv(A):')
disp(pinv(A))
disp('Approximate solution, betaPinv:')
disp(betaPinv)
disp('Residual vector, Y - A*betaPinv:')
disp(Y-A*betaPinv)
disp('Norm of the residual vector, ||Y - A*betaPinv||:')
disp(norm(Y-A*betaPinv))
disp('Mean of the residual vector elements, mean(Y - A*betaPinv):')
disp(mean(Y-A*betaPinv))

%%
% To derive an explicit expression for the pseudoinverse matrix, minimize
% $\chi^2 \equiv ||\mathbf{Y} - \mathbf{A}\mathbf{\beta}||$ with respect to
% $\mathbf{\beta}$. Express the norm as $||\mathbf{v}|| =
% \mathbf{v}^\mathrm{T}\mathbf{v}$:
%
% $$ \chi^2 = ||\mathbf{Y} - \mathbf{A}\mathbf{\beta}|| =
% \mathbf{Y}^\mathrm{T}\mathbf{Y} -
% \mathbf{Y}^\mathrm{T}\mathbf{A}\mathbf{\beta} - 
% \mathbf{\beta}^\mathrm{T}\mathbf{A}^\mathrm{T}\mathbf{Y} +
% \mathbf{\beta}^\mathrm{T}\mathbf{A}^\mathrm{T}\mathbf{A}\mathbf{\beta} $$
%
% Use the identity
%
% $$ \frac{d(\mathbf{v}^\mathrm{T}\mathbf{A})}{d\mathbf{v}} =
% \mathbf{A}^\mathrm{T} $$
%
% to obtain the derivative,
%
% $$ \frac{d\chi^2}{d\mathbf{\beta}} = -2\mathbf{Y}^\mathrm{T}\mathbf{A} +
% 2\mathbf{\beta}^\mathrm{T}\mathbf{A}^\mathrm{T}\mathbf{A},$$
%
% which we set equal to zero to find the minimum when
%
% $$ \mathbf{\beta} =
% (\mathbf{A}^\mathrm{T}\mathbf{A})^{-1}\mathbf{A}^\mathrm{T}\mathbf{Y}. $$
%
% The pseudoinverse is then
%
% $$ \mathbf{A}^+ \equiv
% (\mathbf{A}^\mathrm{T}\mathbf{A})^{-1}\mathbf{A}^\mathrm{T}. $$
%
% Confirm this for our example.

disp('pinv(A):')
disp(pinv(A))
disp('inv(transpose(A)*A)*transpose(A):')
disp(inv(transpose(A)*A)*transpose(A))

%% 5.2.2 Interpolating, extrapolating and aliasing
% Use the <http://www.mathworks.com/help/matlab/ref/interp1.html |interp1|>
% function to interpolate between two measured values. The default
% algorithm is linear interpolation, but others are available. See the
% documentation for details.
% 
% The MATLAB documentation provides a nice interpolation example that we'll
% duplicate here.

x = 0:pi/4:2*pi;
y = sin(x);
xi = 0:pi/16:2*pi;
yi = interp1(x,y,xi);
xs = 0:pi/100:2*pi;
ys = sin(xs);
plot(x,y,'ko',xi,yi,'r.-',xs,ys,'k-')
xlabel('x')
ylabel('y')
legend('Measured','Interpolated','Smooth')

%%
% You can also use |interp1| to extrapolate. As Hughes and Hase emphasize,
% proceed with caution!

x = 0:0.1:1;
y = x./(1+x);
xl = x(1:5);
yl = y(1:5);
xext = x(6:end);
yext = interp1(xl,yl,xext,'linear','extrap');
plot(xl,yl,'ko',xext,yext,'ro',x,y,'k.-')
xlabel('x')
ylabel('y')
legend('Measured','Extrapolated','Actual',...
    'Location','SouthEast')


%%
% You can also use the results of a curve fit to interpolate and
% extrapolate. Add a bit of noise first to the fit more realistic. The
% <http://www.mathworks.com/help/curvefit/list-of-library-models-for-curve-and-surface-fitting.html
% |rat11|> model has the form
%
% $$ y = \frac{p_1x + p_2}{x + q_1}, $$
%
% and requires an initial guess, which we choose to be $p_1 = q_1 = 1$,
% $p_2 = 0$. The best-fit model allows us to both interpolate and
% extrapolate, and the <http://www.mathworks.com/help/curvefit/predint.html
% |predint|> function provides confidence intervals for these estimates.
% Note that the extrapolated values have greater uncertainty.

xdat = xl + 1e-3*randn(size(xl));
ydat = yl + 1e-3*randn(size(yl));
xdat = xdat(:);
ydat = ydat(:);
fitObj = fit(xdat,ydat,'rat11','StartPoint',[1,0,1]);

disp(fitObj)

plot(xl,yl,'ko',xext,fitObj(xext),'ro',x,y,'k.-')
ci = predint(fitObj,x);
hold on
plot(x,ci,'r:')
hold off
xlabel('x')
ylabel('y')
legend('Measured','Extrapolated','Actual','Prediction bounds (95%)',...
    'Location','SouthEast')

%%
% Now consider aliasing.
t = (0:5);
t = t(:);
v = 5*ones(size(t));
ts = linspace(0,5,512);
ts = ts(:);
vs = [5*ones(size(ts)) 5*cos(2*pi*ts) 5*cos(6*pi*ts)];
plot(t,v,'ko',ts,vs,'-')
ylim([-6 6])
xlabel('Time (s)')
ylabel('Potential (V)')
legend('Measurements','5 V','5cos(2\pi t) V','5cos(6\pi t) V',...
    'Location', 'SE')

%% 5.3 Introduction to the method of least squares and maximum likelihood
% To develop some intuition about the maximum likelihood principle and the
% method of least squares, consider a simple example of estimating the mass
% of a (nominally 1 kg) weight from three independent measurements, where
% the measurement uncertainty is known in advance to be 1 g.

mMeas = 1000 + randn(3,1);
disp(mMeas)

%%
% If we _assume_ that the true mass is |xest|, then we can calculate the
% _likelihood_ that we would obtain the three measurements that we did. The
% logical direction is important here: we start by _assuming_ a model, then
% calculating the likelihood that our measurements would result from this
% model. In this example, the model assumptions are that we have a Gaussian
% distribution with a mean $\mu = \mathtt{xest}$ and standard deviation
% $\sigma = 1$. The distribution tells us the probability of each
% measurement. Try this for |xest = 999|, |1000|, and |1001|.

m = linspace(997,1003);
m = m(:);
mEst = (999:1001);
sigma = 1;
for iEst = 1:3
    figure
    plot(m, normpdf(m, mEst(iEst), sigma), 'k-',...
        mMeas, normpdf(mMeas, mEst(iEst), sigma), 'ro')
    xlabel('Mass (g)')
    ylabel('Probability density (g^{-1})')
    txt = {sprintf('\\mu = %d', mEst(iEst)), ...
        sprintf('p(m_m; \\mu) = %.2g', ...
        prod(normpdf(mMeas, mEst(iEst), sigma)))};
    text(1000,0.03,txt,'HorizontalAlignment','Center')
end

%%
% Since we also assume that the measurements are independent, the
% probability density of getting all three measurements is just the product
% of the three separate probability densities. The resuling joint
% probability density is still Gaussian, with a standard deviation of
% $1/\sqrt{3}$ because it is a product of three Gaussians with $\sigma =
% 1$. The peak of the probability density is the value of $\mu$ with the
% _maximum likelihood_ of yielding our measurements. For reference, the
% measurement values are represented as red circles along the horizontal
% axis.

p = zeros(length(m),length(mMeas));
for iMeas = 1:length(mMeas)
    p(:,iMeas) = normpdf(m,mMeas(iMeas),sigma);
end
figure
plot(m, prod(p,2), 'k-', ...
    mMeas, zeros(size(mMeas)), 'ro')
xlabel('Mass (g)')
ylabel('Joint probability density (g^{-3})')

%%
% It is also the value of the function that minimizes the least-squares
% functional, which is proportional to the negative-log-likelihood
% function.

plot(m, -log(prod(p,2)), 'k-', ...
    mMeas, zeros(size(mMeas)), 'ro')
xlabel('Mass (g)')
ylabel('Negative log-likelihood')

%%
% Return now to the matrix formulation of the linear least-squares problem,
% where we were interested in fitting a line to a set of data $\mathbf{X},
% \mathbf{Y}$ to get the best-fit values $\mathbf{\beta} = [c\
% m]^\mathrm{T}$ for the intercept and slope. In the current example, _the
% slope is zero:_ we set $\mathbf{Y}$ to be a vector of the mass
% measurements, and let $\mathbf{X}$ be a vector that enumerates the
% trials: $x_1 = 1, x_2 =2, x_3 =3.$ The mass should not depend on when we
% measured it, so the problem that we want to solve is
%
% $$ \mathbf{Y} = \mathbf{A} c, $$
%
% where $c$ is just a constant and $\mathbf{A} = \mathbf{X}^0$ is an array
% of ones that has the same size as $\mathbf{X}$ and $\mathbf{Y}$. As we
% saw earlier, the least-squares solution to this problem is $c =
% \mathbf{A}^+\mathbf{Y}$, where $\mathbf{A}^+$ is the pseudoinverse of
% $\mathbf{A}$.

A = ones(size(mMeas));
Apinv = pinv(A);
c = pinv(A)*mMeas;
disp('A:')
disp(A)
disp('pinv(A):')
disp(Apinv)
disp('c:')
disp(c)

%%
% Notice that the inverse matrix is just $\mathbf{A}^+ = [1~1~1]/N$, where
% $N = 3$, the number of measurements, so $\mathbf{A}^+\mathbf{Y}$ is just
% the mean. This justifies our use of the mean to estimate the mass: for a
% set of measurements $\mathbf{Y}$ that are normally distributed with
% constant $\sigma$ and _unknown_ $\mu$, $\hat{\mu} = (1/N)\sum_k^N y_k$
% is the _maximum-likelihood estimate_ of $\mu$.

disp(mean(mMeas))

%% 5.3.1 Example using the method of least squares
% Reproduce Fig. 5.8. By defining |x| as a |11 x 1| column vector and
% |mTrial| as a |1 x 5| row vector, we can use array multiplication to
% compute all five trial lines simultaneously.

x = (0:10);
x = x(:);
y = 2*x;
mTrial = [1.0, 1.5, 2.0, 2.25, 2.75];
yTrial = x*mTrial;

figure
hold on
plot(x,y,'ko')
plot(x,yTrial(:,3),'k-','LineWidth',2)
plot(x,yTrial(:,[2 4]),'k--')
plot(x,yTrial(:,[1 5]),'k:')
hold off
xlabel('x')
ylabel('y')
text(9,2,'(a)');

%%
% Now compute the sum-squared error, SSE $\equiv \sum_i^N (y_i - mx_i)^2$,
% as a function of the parameter $m$. Compute $y_i - mx_i$ for each column
% of |yTrial| by using the
% <https://www.mathworks.com/help/matlab/ref/repmat.html |repmat|> function
% to make |yRep|, an |11 x 5| array with identical colums, all equal to
% |y|; then, take the square of each element; finally, use the
% <https://www.mathworks.com/help/matlab/ref/sum.html |sum|> function to
% compute the SSE for each column in an array.

yRep = repmat(y,1,5);
yResidualTrial = yRep - yTrial;
sumSqErrTrial = sum(yResidualTrial.^2);
disp(sumSqErrTrial)

%%
% Use the same procedure to plot the SSE for each value of $m$ in a |1 x
% 100| row vector |mPlot|. Show SSE for |yTrial| as black circles.
mPlot = linspace(1,3);
yPlot = x*mPlot;
yRep = repmat(y,1,size(mPlot,2));
yResidualPlot = yRep - yPlot;
sumSqErrPlot = sum(yResidualPlot.^2);

figure
hold on
plot(mTrial,sumSqErrTrial,'ko')
plot(mPlot,sumSqErrPlot,'k-')
xlabel('Gradient')
ylabel('\Sigma [y_i - y(x_i)]^2')

%%
% Show the residuals for the different trial values. Resize the figure
% window to accommodate the five stacked plots, and remove the x-axis tick
% labels from all but the bottom plot.

figure
pos = get(gcf,'Position');
pos(4) = 2*pos(4);
set(gcf,'Position',pos)
nTrial = length(mTrial);
for iTrial = 1:nTrial
subplot(nTrial,1,iTrial)
stem(x,yResidualTrial(:,iTrial))
ltxt = sprintf('m_{est} = %.2f, SSE = %.1f', ...
    mTrial(iTrial), sumSqErrTrial(iTrial));
title(ltxt)
if iTrial < nTrial
    set(gca,'XTickLabel',{})
end
end
xlabel('x')

%% 5.4 Performing a least-squares fit to a straight line

%% 5.5 Using graphs to estimate random and systematic errors

%% 5.6 Residuals

%%
% $Id: P233HH3.m 151 2016-02-02 18:38:44Z physcmpt $