add code files and mat files

statistics/code/lsq_error.m

@@ -0,0 +1,6 @@
+function error = lsq_error(parameter, x, y)
+% parameter(1) is the slope
+% parameter(2) is the intercept
+
+f_x = x .* parameter(1) + parameter(2);
+error = mean((f_x - y).^2);

statistics/code/lsq_gradient.m

@@ -0,0 +1,7 @@
+function gradient = lsq_gradient(parameter, x, y)
+h = 1e-6;
+
+partial_m = (lsq_error([parameter(1)+h, parameter(2)],x,y) - lsq_error(parameter,x,y))/ h;
+partial_n = (lsq_error([parameter(1), parameter(2)+h],x,y) - lsq_error(parameter,x,y))/ h;
+
+gradient = [partial_m, partial_n];

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

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);

statistics/code/plot_error_surface.m

@@ -0,0 +1,112 @@
+clear
+close all
+
+%% first, plot the raw data
+load('lin_regression.mat');
+
+figure()
+plot(x,y, 'o')
+xlabel('Input')
+ylabel('Output')
+
+%% plot the error surface
+clear
+load('lin_regression.mat')
+ms = -5:0.25:5;
+ns = -30:1:30;
+
+error_surf = zeros(length(ms), length(ns));
+
+for i = 1:length(ms)
+    for j = 1:length(ns)
+        error_surf(i,j) = lsq_error([ms(i), ns(j)], x, y);
+    end
+end
+
+
+% plot the error surface
+figure()
+[N,M] = meshgrid(ns, ms);
+s = surface(M,N,error_surf);
+xlabel('slope')
+ylabel('intercept')
+zlabel('error')
+view(3)
+% rotate(s, [1 1 0], 25 )
+
+%% Plot the gradient at different points in the surface
+clear
+load('lin_regression.mat')
+
+ms = -1:0.5:5;
+ns = -10:1:10;
+
+error_surf = zeros(length(ms), length(ns));
+gradient_m = zeros(size(error_surf));
+gradient_n = zeros(size(error_surf));
+
+for i = 1:length(ms)
+    for j = 1:length(ns)
+        error_surf(i,j) = lsq_error([ms(i), ns(j)], x, y);
+        grad = lsq_gradient([ms(i), ns(j)], x, y);
+        gradient_m(i,j) = grad(1);
+        gradient_n(i,j) = grad(2);
+    end
+end
+
+figure()
+hold on
+[N, M] = meshgrid(ns, ms);
+surface(M,N, error_surf, 'FaceAlpha', 0.5);
+contour(M,N, error_surf, 50);
+quiver(M,N, gradient_m, gradient_n)
+view(3)
+xlabel('slope')
+ylabel('intercept')
+zlabel('error')
+
+%% do the gradient descent
+clear 
+close all
+
+load('lin_regression.mat')
+
+ms = -1:0.5:5;
+ns = -10:1:10;
+
+position = [-2. 10.];
+gradient = [];
+error = [];
+eps = 0.01;
+
+% claculate error surface
+error_surf = zeros(length(ms), length(ns));
+for i = 1:length(ms)
+    for j = 1:length(ns)
+        error_surf(i,j) = lsq_error([ms(i), ns(j)], x, y);
+    end
+end
+figure()
+hold on
+[N, M] = meshgrid(ns, ms);
+surface(M,N, error_surf, 'FaceAlpha', 0.5);
+view(3)
+xlabel('slope')
+ylabel('intersection')
+zlabel('error')
+
+% do the descent
+
+while isempty(gradient) || norm(gradient) > 0.1
+    gradient = lsq_gradient(position, x,y);
+    error = lsq_error(position, x, y);
+    plot3(position(1), position(2), error, 'o', 'color', 'red')
+    position = position - eps .* gradient; 
+    pause(0.25)
+end
+disp('gradient descent done!')
+disp(strcat('final position: ', num2str(position)))
+disp(strcat('final error: ', num2str(error)))
+
+
+

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)