teaching/scientificComputing
Archived
3
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

linear algebra first scripts

Browse Source
This commit is contained in:
Jan Benda 2014-10-28 17:17:27 +01:00
parent 948e00a21a
commit 490c6161f8
10 changed files with 458 additions and 0 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

82
linearalgebra/simulations/covareigen.m Normal file
View File

@ -0,0 +1,82 @@
function covareigen( x, y, pm, w, h )
% computes covariance matrix from the pairs x, y
% diagonalizes covariance matrix
% x and y are column vectors
% pm: 0 - scatter of data, 1 histogram, 2 multivariate gauus from
% covariance matrix, 1 and 2 require w and h
% w and h are the width and the height (in data units) for plotting
% histograms instead of scatter
% covariance matrix:
cv = cov( [ x y ] );
% eigen values:
[ v , d] = eig( cv )
s = sign( v );
s(1,:) = s(2,:);
v = v .* s;
% plots:
subplot( 2, 2, 1 );
hold on;
if (pm > 0) & (nargin == 5)
if pm == 1
% histogram of data:
xp = -w:0.5:w;
yp = -h:0.5:h;
[n,c] = hist3([x y], { xp, yp } );
contourf( c{1}, c{2}, n' );
else
xp = -w:0.1:w;
yp = -h:0.1:h;
[xg,yg] = meshgrid( xp, yp );
xy = [ xg(:) yg(:) ];
gauss = reshape( exp(-0.5*diag(xy*inv(cv)*xy'))/sqrt((2.0*pi)^2.0*det(cv)), size( xg ) );
contourf( xp, yp, gauss )
end
colormap( 'gray' );
else
% scatter plot:
scatter( x, y, 'b', 'filled', 'MarkerEdgeColor', 'white' );
end
% plot eigenvectors:
quiver( ones( 1, 2 ).*mean( x ), ones( 1, 2 )*mean(y), v(1,:).*sqrt(diag(d))', v(2,:).*sqrt(diag(d))', 'r', 'LineWidth', 3, 'AutoScale', 'off', 'AutoScaleFactor', 1.0, 'MaxHeadSize', 0.7 )
axis( 'equal' );
hold off;
subplot( 2, 2, 3 );
hist( x, 50, 'b' );
% sort the eigenvalues:
[d,inx] = sort( diag(d), 'descend' );
subplot( 2, 2, 2 );
hold on;
% subtract means:
x = x - mean( x );
y = y - mean( y );
% project onto eigenvectors:
nx = [ x y ] * v(:,inx(1));
ny = [ x y ] * v(:,inx(2));
cv = cov( [nx, ny] )
[ v , d] = eig( cv )
if (pm > 0) & (nargin == 5)
if pm == 1
[n,c] = hist3([nx ny], { xp, yp } );
contourf( c{1}, c{2}, n' );
else
gauss = reshape( exp(-0.5*diag(xy*inv(cv)*xy'))/sqrt((2.0*pi)^2.0*det(cv)), size( xg ) );
contourf( xp, yp, gauss )
end
else
scatter( nx, ny, 'b', 'filled', 'MarkerEdgeColor', 'white' );
end
axis( 'equal' );
hold off;
subplot( 2, 2, 4 );
hist( nx, 50, 'b' );
end

45
linearalgebra/simulations/covareigen3.m Normal file
View File

@ -0,0 +1,45 @@
function covareigen3( x, y, z )
% computes covariance matrix from the triples x, y, z
% diagonalizes covariance matrix
% x, y and z are column vectors
% covariance matrix:
cv = cov( [ x y z ] );
% eigen values:
[ v , d] = eig( cv )
% plots:
subplot( 1, 2, 1 );
hold on;
% scatter plot:
view( 3 );
scatter3( x, y, z, 0.1, 'b', 'filled', 'MarkerEdgeColor', 'blue' );
xlabel( 'x' );
ylabel( 'y' );
zlabel( 'z' );
grid on;
% plot eigenvectors:
quiver3( ones( 1, 3 ).*mean( x ), ones( 1, 3 )*mean(y), ones( 1, 3 )*mean(z), v(1,:).*sqrt(diag(d))', v(2,:).*sqrt(diag(d))', v(3,:).*sqrt(diag(d))', 'r', 'LineWidth', 3, 'AutoScale', 'off', 'AutoScaleFactor', 1.0, 'MaxHeadSize', 0.7 )
%axis( 'equal' );
hold off;
% sort the eigenvalues:
[d,inx] = sort( diag(d), 'descend' );
subplot( 2, 2, 2 );
hold on;
% subtract means:
x = x - mean( x );
y = y - mean( y );
% project onto eigenvectors:
nx = [ x y z ] * v(:,inx(1));
ny = [ x y z ] * v(:,inx(2));
scatter( nx, ny, 'b', 'filled', 'MarkerEdgeColor', 'white' );
axis( 'equal' );
hold off;
end

14
linearalgebra/simulations/covareigen3examples.m Normal file
View File

@ -0,0 +1,14 @@
n = 10000;
% three distributions:
x = randn( n, 1 );
y = randn( n, 1 );
z = randn( n, 1 );
f = figure( 'Position', [ scrsz(3)/2 scrsz(4)/2 scrsz(3)/2 scrsz(4)/2 ]);
d = [ 0 4 -4 8 ];
for k = 1:4
x((k-1)*n/4+1:k*n/4) = x((k-1)*n/4+1:k*n/4) + d(k);
y((k-1)*n/4+1:k*n/4) = y((k-1)*n/4+1:k*n/4) - d(k);
z((k-1)*n/4+1:k*n/4) = z((k-1)*n/4+1:k*n/4) + d(k);
end
covareigen3( x, y, z );

35
linearalgebra/simulations/covareigenexamples.m Normal file
View File

@ -0,0 +1,35 @@
scrsz = get( 0, 'ScreenSize' );
set( 0, 'DefaultFigurePosition', [ scrsz(3)/2 scrsz(4)/2 scrsz(3)/2 scrsz(4)/2 ] );
% correlation coefficients:
n = 10000;
x = randn( n, 1 );
f = figure( 1 );
for r = 0.01:0.199:1
fprintf( 'Correlation = %g\n', r );
clf( f );
y = r*x + sqrt(1-r^2)*randn( n, 1 );
covareigen( x, y, 2, 5.0, 3.0 );
key = waitforbuttonpress;
end
return
% two distributions:
n = 10000;
x1 = randn( n/2, 1 );
y1 = randn( n/2, 1 );
x2 = randn( n/2, 1 );
y2 = randn( n/2, 1 );
f = figure( 1 );
pause( 'on' );
for d = 0:1:5
fprintf( 'Distance = %g\n', d );
clf( f );
d2 = d / sqrt( 2.0 );
x = [ x1; x2 ];
y = [ y1+d2; y2-d2 ];
scrsz = get(0,'ScreenSize');
covareigen( x, y, 1, 10.0, 7.0 );
%key = waitforbuttonpress;
pause( 1.0 );
end

74
linearalgebra/simulations/matrix2dtrafos.m Normal file
View File

@ -0,0 +1,74 @@
m = [ 1 0; 0 1];
matrixbox( m );
title( 'Identity' );
waitforbuttonpress;
m = [ 2 0; 0 1];
matrixbox( m );
title( 'Scale x' );
waitforbuttonpress;
m = [ 1 0; 0 2];
matrixbox( m );
title( 'Scale y' );
waitforbuttonpress;
m = [ 2 0; 0 2];
matrixbox( m );
title( 'Scale both' );
waitforbuttonpress;
m = [ -1 0; 0 1];
matrixbox( m );
title( 'Flip x' );
waitforbuttonpress;
m = [ 1 0; 0 -1];
matrixbox( m );
title( 'Flip y' );
waitforbuttonpress;
m = [ -1 0; 0 -1];
matrixbox( m );
title( 'Flip both' );
waitforbuttonpress;
m = [ 1 0; 1 1];
matrixbox( m );
title( 'Shear x' );
waitforbuttonpress;
m = [ 1 1; 0 1];
matrixbox( m );
title( 'Shear y' );
waitforbuttonpress;
phi = 0.1667*pi;
m = [ cos(phi) -sin(phi); sin(phi) cos(phi)];
matrixbox( m );
title( 'Rotate 30' );
waitforbuttonpress;
phi = 0.333*pi;
m = [ cos(phi) -sin(phi); sin(phi) cos(phi)];
matrixbox( m );
title( 'Rotate 60' );
waitforbuttonpress;
phi = 0.5*pi;
m = [ cos(phi) -sin(phi); sin(phi) cos(phi)];
matrixbox( m );
title( 'Rotate 90' );
waitforbuttonpress;
phi = 0.75*pi;
m = [ cos(phi) -sin(phi); sin(phi) cos(phi)];
matrixbox( m );
title( 'Rotate 135' );
waitforbuttonpress;
phi = 1.0*pi;
m = [ cos(phi) -sin(phi); sin(phi) cos(phi)];
matrixbox( m );
title( 'Rotate 180' );
waitforbuttonpress;

24
linearalgebra/simulations/matrixbox.m Normal file
View File

@ -0,0 +1,24 @@
function matrixbox( m )
% visualizes the effect of a matrix m on a set of vectors forming a box
% m: a 2x2 matrix
v = [ 0 1 1 0 0; 0 0 1 1 0];
w = m*v;
clf;
hold on;
%set(gca, 'Xtick', [], 'Ytick', [], 'box', 'off')
% axis:
plot( [-2 2], [0 0], 'k', 'LineWidth', 1 );
plot( [0 0], [-2 2], 'k', 'LineWidth', 1 );
% old box:
plot( v(1,:), v(2,:), 'k', 'LineWidth', 1.0 )
quiver( [v(1,1)], [v(2,1)], [v(1,3)], [v(2,3)], 1.0, 'k', 'LineWidth', 2.0 );
scatter( [v(1,2)], [v(2,2)], 60.0, 'filled', 'k' );
% transfomred box:
plot( w(1,:), w(2,:), 'b', 'LineWidth', 2.0 )
quiver( [w(1,1)], [w(2,1)], [w(1,3)], [w(2,3)], 1.0, 'b', 'LineWidth', 3.0 );
scatter( [w(1,2)], [w(2,2)], 100.0, 'filled', 'b' );
hold off;
xlim( [ -2 2 ] );
ylim( [ -2 2 ] );
end

22
linearalgebra/simulations/plotvector.m Normal file
View File

@ -0,0 +1,22 @@
function plotvector( v, s, sp )
% plot an arrow for the 2D-vector v
% v: the 2D vector
% s: a string passed to quiver for the color
% sp: if this string is set it defines the color of some helper lines
if nargin < 2
s = 'b';
end
quiver( [0.0], [0.0], [v(1)], [v(2)], 1.0, s );
grid on;
if nargin == 3
hh = ishold;
if ~hh
hold on;
end
plot( [0.0 v(1) v(1)], [0.0 0.0 v(2)], sp );
if ~hh
hold off;
end
end
end

24
linearalgebra/simulations/plotvector3.m Normal file
View File

@ -0,0 +1,24 @@
function plotvector3( v, s, sp )
% plot an arrow for the 3D-vector v
% v: the 3D vector
% s: a string passed to quiver for the color
% sp: if this string is set it defines the color of some helper lines
if nargin < 2
s = 'b';
end
quiver3( [0.0], [0.0], [0.0], [v(1)], [v(2)], [v(3)], 1.0, s );
view( 42, 12 );
grid on
if nargin == 3
hh = ishold;
if ~hh
hold on;
end
plot3( [v(1) v(1)], [v(2) v(2)], [0.0 v(3)], sp );
plot3( [0.0 v(1) v(1)], [v(2) v(2) 0.0], [0.0 0.0 0.0], sp );
if ~hh
hold off;
end
end
end

127
linearalgebra/simulations/spikesortingwave.m Normal file
View File

@ -0,0 +1,127 @@
% generate spikes:
n = 1000;
misi = 0.01;
isis = exprnd( misi, n, 1 );
isis = isis + 0.01;
spikes = cumsum( isis );
p = rand( size( spikes ) );
% spike waveforms:
dt = 0.0001;
x = -0.01:dt:0.01;
y1 = 2.0*exp( -(x-0.0003).^2/2.0/0.0005^2 ) - 1.4*exp( -(x-0.0005).^2/2.0/0.0005^2 );
y2 = exp( -x.^2/2.0/0.002^2 ).*cos(2.0*pi*x/0.005);
p12 = 0.5;
% voltage trace:
noise = 0.2;
t = 0:dt:spikes(end)+3.0*x(end);
v = noise*randn( 1, length( t ) );
for k = 1:length( spikes )
inx = ceil( spikes(k)/dt );
if p(k) < p12
v(inx:inx+length(x)-1) = v(inx:inx+length(x)-1) + y1;
else
v(inx:inx+length(x)-1) = v(inx:inx+length(x)-1) + y2;
end
end
figure( 1 );
clf;
plot( t(t<1.0), v(t<1.0) );
hold on;
% find peaks:
thresh = 0.7;
inx = find( v(2:end-1) > thresh & v(1:end-2)<v(2:end-1) & v(2:end-1) > v(3:end) ) + 1;
spiket = t(inx);
spikev = v(inx);
tinx = inx;
for k=1:2
inx = find( ( spiket(2:end-1)-spiket(1:end-2)>0.005 | spikev(2:end-1) > spikev(1:end-2) ) & ( spiket(3:end)-spiket(2:end-1)>0.005 | spikev(2:end-1) > spikev(3:end) ) )+1;
spiket = spiket(inx);
spikev = spikev(inx);
tinx = tinx(inx);
end
scatter( spiket(spiket<1.0), spikev(spiket<1.0), 100.0, 'r', 'filled' );
%scatter( t(tinx), v(tinx), 100.0, 'r', 'filled' );
hold off;
% spike waveform snippets:
w = ceil( 0.005/dt );
vs = [];
for k=1:length(tinx)
vs = [ vs; v(tinx(k)-w:tinx(k)+w-1) ];
end
ts = t(1:size(vs,2));
ts = ts - ts(floor(length(ts)/2));
figure( 2 );
clf;
hold on
for k=1:size(vs,1)
plot( ts, vs(k,:), '-b' );
end
hold off
% pca:
cv = cov( vs );
[ ev , ed ] = eig( cv );
[d,dinx] = sort( diag(ed), 'descend' );
% spectrum of eigenvalues:
figure( 3 );
clf;
subplot( 4, 1, 1 );
scatter( 1:length(d), d, 'b', 'filled' );
% project onto eigenvectors:
nx = vs * ev(:,dinx(1));
ny = vs * ev(:,dinx(2));
%scatter( nx, ny, 'b', 'filled', 'MarkerEdgeColor', 'white' );
%hold on;
% features:
subplot( 4, 2, 5 );
plot( 1000.0*ts, ev(:,dinx(1)), 'g', 'LineWidth', 2 );
subplot( 4, 2, 6 );
plot( 1000.0*ts, ev(:,dinx(2)), 'r', 'LineWidth', 2 );
% clustering:
%kx = kmeans( [ nx, ny ], 2 );
% nx smaller or greater a threshold:
kthresh = 1.6;
kx = ones( size( vs, 1 ), 1 );
kx(nx>kthresh) = 2;
subplot( 4, 1, 2 );
scatter( nx(kx==1), ny(kx==1), 'g', 'filled', 'MarkerEdgeColor', 'white' );
hold on;
scatter( nx(kx==2), ny(kx==2), 'r', 'filled', 'MarkerEdgeColor', 'white' );
hold off;
% show sorted spike waveforms:
subplot( 4, 2, 7 );
hold on
kinx1 = find(kx==1);
for k=1:length(kinx1)
plot( 1000.0*ts, vs(kinx1(k),:), '-g' );
end
plot( 1000.0*x, y1, '-k', 'LineWidth', 2 );
xlim( [ 1000.0*ts(1) 1000.0*ts(end) ] )
hold off
subplot( 4, 2, 8 );
hold on
kinx2 = find(kx==2);
for k=1:length(kinx2)
plot( 1000.0*ts, vs(kinx2(k),:), '-r' );
end
plot( 1000.0*x, y2, '-k', 'LineWidth', 2 );
xlim( [ 1000.0*ts(1) 1000.0*ts(end) ] )
hold off
% spike trains:
figure( 1 );
hold on;
six1 = find( spiket(kinx1)<1.0 );
scatter( spiket(kinx1(six1)), spikev(kinx1(six1)), 100.0, 'g', 'filled' );
six2 = find( spiket(kinx2)<1.0 );
scatter( spiket(kinx2(six2)), spikev(kinx2(six2)), 100.0, 'r', 'filled' );
hold off;

11
linearalgebra/simulations/spiketime.m Normal file
View File

@ -0,0 +1,11 @@
p = 0.0;
n=0;
for k = 1:length( spikes )
%a= 0.5*(spikes(k+1)+spikes(k));
%fprintf( 'set label %d "$T_{%d}$" at %g, -0.5 center\n', k+2, k, a );
%fprintf( 'set arrow %d from %g, 0.3 to %g, 0.3 heads\n', k+2, spikes(k), spikes(k+1) );
%fprintf( '%g %d\n%g %d\n\n', p, n, spikes(k), n );
fprintf( '%g %d\n', spikes(k), n+1 );
n = n+1;
p = spikes(k);
end