[regression] finished exercise

regression/code/gradientDescent.m

@@ -13,7 +13,7 @@ function [p, ps, mses] = gradientDescent(x, y, func, p0, epsilon, threshold)
 %            mses: vector with the corresponding mean squared errors
 
   p = p0;
-  gradient = ones(1, length(p0)) * 1000.0;
+  gradient = ones(size(p0)) * 1000.0;
   ps = [];
   mses = [];
   while norm(gradient) > threshold
@@ -29,13 +29,12 @@ function mse = meanSquaredError(x, y, func, p)
 end
 
 function gradmse = meanSquaredGradient(x, y, func, p)
-  gradmse = zeros(size(p, 1), size(p, 2));
+  gradmse = zeros(size(p));
   h = 1e-7;               % stepsize for derivatives
+  ph = eye(length(p))*h;  % ... for each parameter
   mse = meanSquaredError(x, y, func, p);
   for i = 1:length(p)     % for each coordinate ...
-    pi = p;
+    msepi = meanSquaredError(x, y, func, p + ph(:,i));
-    pi(i) = pi(i) + h;    %   displace i-th parameter
-    msepi = meanSquaredError(x, y, func, pi);
     gradmse(i) = (msepi - mse)/h;
   end
 end

regression/exercises/exp2errorsurface.m

@@ -0,0 +1,15 @@
+expdecaydata;                               % generate data
+close all;
+
+cs = 0.5:0.02:1.5;                          % range of c values
+taus = 4.0:0.2:16.0;                        % range of tau values
+msesurf = zeros(length(taus), length(cs));  % mean squared error
+for i = 1:length(taus)
+    for k = 1:length(cs)
+        msesurf(i, k) = mean((voltage - ...
+                              exp2func(time, [cs(k), taus(i)])).^2);
+    end
+end
+%surface(cs, taus, mses)
+%contour(cs, taus, mses)
+contourf(cs, taus, log10(msesurf))

regression/exercises/exp2fit.m

@@ -0,0 +1,15 @@
+expdecaydata;
+close all;
+
+p0 = [0.5, 6.0];
+pest = lsqcurvefit(@(p, x) exp2func(x, p), p0, time, voltage);
+fprintf('c=%4.1fmV, tau=%4.1fms\n', pest(1), pest(2));
+
+hold on;
+plot(time, voltage)
+plot(time, exp2func(time, pest), 'LineWidth', 2)
+hold off;
+xlabel('time [ms]')
+ylabel('voltage [mv]')
+savefigpdf(gcf, 'exp2fit.pdf', 15, 8.5);
+

regression/exercises/exp2func.m

@@ -0,0 +1,4 @@
+function x = exp2func(t, p)
+% return the exponential function x = c*e^{-t/tau}
+  x = p(1)*exp(-t./p(2));
+end

regression/exercises/exp2plot.m

@@ -0,0 +1,27 @@
+...
+    
+p0 = [0.5, 6.0];
+[pest, ps, mses] = gradientdescent(time, voltage, @exp2func, p0, 1.0, 0.00001);
+fprintf('c=%4.1fmV, tau=%4.1fms\n', pest(1), pest(2));
+
+subplot(1, 2, 1);
+hold on;
+%surface(cs, taus, msesurf);
+%plot3(ps(1,:), ps(2,:), mses, '.r');
+%view(80, 35);
+contourf(cs, taus, log10(msesurf));
+plot(ps(1,:), ps(2,:), '.r');
+hold off;
+xlabel('c [mV]')
+ylabel('tau [ms]')
+zlabel('mse [mV^2]')
+
+subplot(1, 2, 2);
+hold on;
+plot(time, voltage)
+plot(time, exp2func(time, pest), 'LineWidth', 2)
+hold off;
+xlabel('time [ms]')
+ylabel('voltage [mv]')
+savefigpdf(gcf, 'exp2plot.pdf', 15, 8.5);
+

regression/exercises/expdecaydata.m

@@ -2,7 +2,7 @@ tau = 10.0;                      % membrane time constant in ms
 dt = 0.05;                       % sampling interval in ms
 noisesd = 0.05;                  % measurement noise in mV
 
-time = 0.0:dt:5*tau;             % time vector
+time = 0.0:dt:6*tau;             % time vector
 voltage = expdecay(time, tau);   % exponential decay
 voltage = voltage + noisesd*randn(1, length(voltage));  % plus noise
 

regression/exercises/expdecayplot.m

@@ -8,7 +8,7 @@ thresh = 0.00001;
 subplot(2, 2, 1);                         % top left panel
 hold on;
 plot(taus, '-o');
-plot([1, length(taus)], [tau, tau], 'k'); % line indicating true tau value
+plot([1, length(taus)], [tau, tau], 'k'); % line indicating true tau
 hold off;
 xlabel('Iteration');
 ylabel('tau');
@@ -19,11 +19,9 @@ ylabel('MSE');
 subplot(1, 2, 2);                         % right panel
 hold on;
 % generate x-values for plottig the fit:
-tt = 0.0:0.01:max(time);
+plot(time, voltage);                      % plot original data
-xx = expdecay(tt, tauest);
+plot(time, expdecay(time, tauest), 'LineWidth', 2); % plot fit
-plot(time, voltage, '.');                 % plot original data
-plot(tt, xx, '-r');                       % plot fit
 xlabel('Time [ms]');
 ylabel('Voltage [mV]');
 legend('data', 'fit', 'location', 'northeast');
-pause
+savefigpdf(gcf, 'expdecayplot.pdf', 15, 10);

regression/exercises/gradientdescent-1.tex

@@ -34,24 +34,24 @@
     \part Implement (and document!) the exponential function 
     \begin{equation}
       \label{expfunc}
-      x(t) = e^{-t/\tau}
+      x(t) = e^{-t/\tau} \quad , \qquad \tau \in \reZ
     \end{equation}
     with the membrane time constant $\tau$ as a matlab function
-    \code{expdecay(t, tau)} that takes as arguments a vector of
+    \code{expdecay(t, tau)} that takes as arguments a vector of time
-    time points and the membrane time constant. The function returns
+    points and the value of the membrane time constant. The function
-    \eqnref{expfunc} computed for each time point as a vector.
+    returns \eqnref{expfunc} computed for each time point as a vector.
     \begin{solution}
       \lstinputlisting{expdecay.m}
     \end{solution}
 
     \part Let's first generate the data. Set the membrane time
     constant to 10\,ms. Generate a time vector with sample times
-    between zero and five times the membrane time constant and a
+    between zero and six times the membrane time constant and a
     sampling interval of 0.05\,ms. Then compute a vector containing
     the corresponding measurements of the membrane potential using the
     \code{expdecay()} function and adding some measurement noise with
-    a standard deviation of 0.05\.mV (\code{randn()} function). Also
+    a standard deviation of 0.05\.mV (\code{randn()} function).
-    plot the data.
+    Plot the data.
     \begin{solution}
       \lstinputlisting{expdecaydata.m}
     \end{solution}
@@ -62,9 +62,9 @@
     for the estimation of the membrane time constant, the $\epsilon$
     factor, and the threshold for the length of the gradient where to
     terminate the algorithm. The function should return the estimated
-    membrane time constant at the minimum of the mean squared error, a
+    membrane time constant at the minimum of the mean squared error,
-    vector with the time constants, and a vector with the mean squared
+    and a vector with the time constants as well as a vector with the
-    errors for each step of the algorithm.
+    mean squared errors for each step of the algorithm.
     \begin{solution}
       \lstinputlisting{expdecaydescent.m}
     \end{solution}
@@ -78,26 +78,90 @@
     for each iteration step, (ii) the mean squared error for each
     iteration step, and (iii) the measured data and the fitted
     exponential function.
+    \newsolutionpage
     \begin{solution}
       \lstinputlisting{expdecayplot.m}
+      \includegraphics{expdecayplot}
     \end{solution}
 
   \end{parts}
 
+  \continuepage
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
   \question \qt{Read sections 8.6 -- 8.8 of chapter 8 on ``optimization
     and gradient descent!}\vspace{-3ex}
 
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
-  \question \qt{Fitting the full exponential function}
+  \question \qt{Fitting an exponential function with two parameters}
+  Usually we do not know the exact value of the voltage step. So let's
+  make this a paramter, too:
+  \begin{equation}
+    \label{exp2func}
+    x(t) = c e^{-t/\tau} \quad , \qquad c, \tau \in \reZ
+  \end{equation}
+  
   \begin{parts}
+    \part Implement the exponential function \eqref{exp2func} as a
+    matlab function \code{exp2func(t, p)} that takes as arguments a
+    vector of time points and a vector with the two parameter values
+    $c$ and $\tau$. The function returns \eqnref{exp2func} computed for
+    each time point as a vector.
+    \begin{solution}
+      \lstinputlisting{exp2func.m}
+    \end{solution}
+
+    \part Plot the error surface, the mean squared error, for a range
+    of values for the two parameters $c$ and $\tau$. Try the
+    \code{surface()}, the \code{contour()} as well as the
+    \code{contourf()} functions to visualize the error surface for the
+    data from the previous exercise. You may also want to try to plot
+    the logarithm of the mean squared error.
+    \begin{solution}
+      \lstinputlisting{exp2errorsurface.m}
+    \end{solution}
+
+    \newsolutionpage
+    \part Implement the gradient descent algorithm for finding the
+    least squares of an arbitrary function. The function takes as
+    arguments the measured data, a vector holding the initial values
+    of the function parameters, the $\epsilon$ factor, and the
+    threshold for the length of the gradient where to terminate the
+    algorithm. The function should return the estimated parameter
+    values in a vector, and a 2D vector with the parameter values as
+    well as a vector with the mean squared errors for each step of the
+    algorithm.
+    \begin{solution}
+      \lstinputlisting{gradientdescent.m}
+    \end{solution}
+
+    \part Call the gradient descent function with the generated data.
+    Set $\epsilon=1.0$ and the threshold value initially to 0.001.
+
+    \part Plot the path of the gradient descent algorithm into the
+    error surface.
+
+    Why appears (potentially) the path of the gradient descent
+    algorithm not perpendicular to the contour lines of the error
+    surface?
+
+    \part Plot the data together with the fitted exponential function.
+
+    \part Then adapt $\epsilon$ and the threshold value to improve the
+    convergence of the algortihm to the minimum of the cost function.
+    \begin{solution}
+      \lstinputlisting{exp2plot.m}
+      \includegraphics{exp2plot}
+    \end{solution}
 
-    \part Use the functions \code{polyfit()} and \code{lsqcurvefit()}
+    \part Use the function \code{lsqcurvefit()} provided by matlab to
-    provided by matlab to find the slope and intercept of a straight
+    find $c$ and $\tau$ for the exponential function describing the
-    line that fits the data. Compare the resulting fit parameters of
+    data.  Compare the resulting fit parameters of this function with
-    those functions with the ones of your gradient descent algorithm.
+    the ones of your gradient descent algorithm.  Plot the data
+    together with the fitted exponential function.
     \begin{solution}
-      \lstinputlisting{linefit.m}
+      \lstinputlisting{exp2fit.m}
+      \includegraphics{exp2fit}\\
+      Resulting values for $c$ and $\tau$ are similar but not identical.
     \end{solution}
 
   \end{parts}

regression/exercises/gradientdescent.m

@@ -0,0 +1,39 @@
+function [p, ps, mses] = gradientdescent(t, x, func, p0, epsilon, threshold)
+% Gradient descent for leas squares fit.
+%
+% Arguments: t, vector of time points.
+%            x, vector of the corresponding measured data values.
+%            p0, vector with initial parameter values.
+%            func, function to be fit to the data. 
+%            epsilon: factor multiplying the gradient.
+%            threshold: minimum value for gradient
+%
+% Returns:   p, vector with the final estimates of the parameter values.
+%            ps: 2D vector with all the parameter values traversed.
+%            mses: vector with the corresponding mean squared errors
+  p = p0;
+  gradient = ones(size(p0)) * 1000.0;
+  ps = [];
+  mses = [];
+  while norm(gradient) > threshold
+      ps = [ps, p(:)];
+      mses = [mses, meansquarederror(t, x, func, p)];
+      gradient = msegradient(t, x, func, p);
+      p = p - epsilon * gradient; 
+  end
+end
+
+function mse = meansquarederror(t, x, func, p)
+  mse = mean((x - func(t, p)).^2);
+end
+
+function gradmse = msegradient(t, x, func, p)
+  gradmse = zeros(size(p));
+  h = 1e-3;   % stepsize for derivative
+  ph = eye(length(p))*h;  % ... for each parameter
+  mse = meansquarederror(t, x, func, p);
+  for i = 1:length(p)     % for each coordinate ...
+    msepi = meansquarederror(t, x, func, p(:) + ph(:,i));
+    gradmse(i) = (msepi - mse)/h;
+  end
+end

regression/exercises/savefigpdf.m

@@ -0,0 +1,28 @@
+function savefigpdf(fig, name, width, height)
+% Saves figure fig in pdf file name.pdf with appropriately set page size
+% and fonts
+
+  % default width:
+  if nargin < 3
+      width = 11.7;
+  end
+  % default height:
+  if nargin < 4
+      height = 9.0;
+  end
+
+  % paper:
+  set(fig, 'PaperUnits', 'centimeters');
+  set(fig, 'PaperSize', [width height]);
+  set(fig, 'PaperPosition', [0.0 0.0 width height]);
+  set(fig, 'Color', 'white')
+
+  % font:
+  set(findall(fig, 'type', 'axes'), 'FontSize', 12)
+  set(findall(fig, 'type', 'text'), 'FontSize', 12)
+
+  % save:
+  saveas(fig, name, 'pdf')
+
+end
+