From 7f52920b0bc861e22e841b7d3a0fed4d6c1bab5e Mon Sep 17 00:00:00 2001
From: Jan Grewe <jan.grewe@g-node.org>
Date: Mon, 26 Oct 2015 18:35:00 +0100
Subject: [PATCH] solutions for fitting exercises

---
 statistics/code/boltzmann.m                  |  7 ++++
 statistics/code/create_linear_data.m         |  7 ++++
 statistics/code/estimate_regression.m        |  7 ++++
 statistics/code/exponential.m                |  5 +++
 statistics/code/lsq_gradient_sigmoid.m       |  9 ++++
 statistics/code/lsq_sigmoid_error.m          |  8 ++++
 statistics/code/sigmoidal_gradient_descent.m | 44 ++++++++++++++++++++
 7 files changed, 87 insertions(+)
 create mode 100644 statistics/code/boltzmann.m
 create mode 100644 statistics/code/create_linear_data.m
 create mode 100644 statistics/code/estimate_regression.m
 create mode 100644 statistics/code/exponential.m
 create mode 100644 statistics/code/lsq_gradient_sigmoid.m
 create mode 100644 statistics/code/lsq_sigmoid_error.m
 create mode 100644 statistics/code/sigmoidal_gradient_descent.m

diff --git a/statistics/code/boltzmann.m b/statistics/code/boltzmann.m
new file mode 100644
index 0000000..8190460
--- /dev/null
+++ b/statistics/code/boltzmann.m
@@ -0,0 +1,7 @@
+function y = boltzmann(parameter, x)
+% parameter 1: alpha
+% parameter 2: k
+% parameter 3: x_0
+% parameter 4: y_0
+
+y = (parameter(1) ./ (1 + exp(-parameter(2) .* (x - parameter(3))))) + parameter(4);
\ No newline at end of file
diff --git a/statistics/code/create_linear_data.m b/statistics/code/create_linear_data.m
new file mode 100644
index 0000000..f929f23
--- /dev/null
+++ b/statistics/code/create_linear_data.m
@@ -0,0 +1,7 @@
+function y = create_linear_data(x)
+
+m = 2.5;
+n = -0.35;
+d = 2.5;
+
+y = x .* m + n + randn(size(x)) .* d;
diff --git a/statistics/code/estimate_regression.m b/statistics/code/estimate_regression.m
new file mode 100644
index 0000000..fd3d2c4
--- /dev/null
+++ b/statistics/code/estimate_regression.m
@@ -0,0 +1,7 @@
+function param = estimate_regression(x,y,p_0)
+
+objective_function = @(p)lsq_error(p, x, y);
+param = fminunc(objective_function, p_0);
+disp(param)
+param1 = fminsearch(objective_function, p_0);
+disp(param1)
\ No newline at end of file
diff --git a/statistics/code/exponential.m b/statistics/code/exponential.m
new file mode 100644
index 0000000..1bb2a04
--- /dev/null
+++ b/statistics/code/exponential.m
@@ -0,0 +1,5 @@
+function y = exponential(parameter, x)
+% Function implements an exponential function with two parameters
+% controlling the amplitude and the time constant.
+
+y = parameter(1) .* exp(x./parameter(2));
\ No newline at end of file
diff --git a/statistics/code/lsq_gradient_sigmoid.m b/statistics/code/lsq_gradient_sigmoid.m
new file mode 100644
index 0000000..05e9721
--- /dev/null
+++ b/statistics/code/lsq_gradient_sigmoid.m
@@ -0,0 +1,9 @@
+function gradient = lsq_gradient_sigmoid(parameter, x, y)
+h = 1e-6;
+
+gradient = zeros(size(parameter));
+for i = 1:length(parameter)
+    parameter_h = parameter;
+    parameter_h(i) = parameter_h(i) + h;
+    gradient(i) = (lsq_sigmoid_error(parameter_h, x, y) - lsq_sigmoid_error(parameter, x, y)) / h;
+end
\ No newline at end of file
diff --git a/statistics/code/lsq_sigmoid_error.m b/statistics/code/lsq_sigmoid_error.m
new file mode 100644
index 0000000..5f05f21
--- /dev/null
+++ b/statistics/code/lsq_sigmoid_error.m
@@ -0,0 +1,8 @@
+function error = lsq_sigmoid_error(parameter, x, y)
+% p(1) the amplitude
+% p(2) the slope
+% p(3) the x-shift
+% p(4) the y-shift
+
+y_est = parameter(1)./(1+ exp(-parameter(2) .* (x - parameter(3)))) + parameter(4);
+error = mean((y_est - y).^2);
\ No newline at end of file
diff --git a/statistics/code/sigmoidal_gradient_descent.m b/statistics/code/sigmoidal_gradient_descent.m
new file mode 100644
index 0000000..9763f5a
--- /dev/null
+++ b/statistics/code/sigmoidal_gradient_descent.m
@@ -0,0 +1,44 @@
+
+
+%% fit the sigmoid
+
+clear 
+close all
+
+load('iv_curve.mat')
+
+figure()
+plot(voltage, current, 'o')
+xlabel('voltate [mV]')
+ylabel('current [pA]')
+
+% amplitude, slope, x-shift, y-shift
+%parameter = [10 0.25 -50, 2.5];
+parameter = [20 0.5 -50, 2.5];
+
+eps = 0.1;
+% do the descent
+gradient = [];
+steps = 0;
+error = [];
+
+while isempty(gradient) || norm(gradient) > 0.01
+    steps = steps + 1; 
+    gradient = lsq_gradient_sigmoid(parameter, voltage, current);
+    error(steps) = lsq_sigmoid_error(parameter, voltage, current);
+    parameter = parameter - eps .* gradient; 
+end
+plot(1:steps, error)
+
+disp('gradient descent done!')
+disp(strcat('final position: ', num2str(parameter)))
+disp(strcat('final error: ', num2str(error(end))))
+
+%% use fminsearch
+parameter = [10 0.5 -50, 2.5];
+
+objective_function = @(p)lsq_sigmoid_error(p, voltage, current);
+param = fminunc(objective_function, parameter);
+disp(param)
+param1 = fminsearch(objective_function, parameter);
+disp(param1)