updated functional dir - nishant

This commit is contained in:
tekhnofiend
2007-11-16 23:49:48 +00:00
parent f1ac910379
commit 7fdb4a4a5c
8 changed files with 389 additions and 23 deletions
+13
View File
@@ -0,0 +1,13 @@
function [filtered_data] = CAR(data)
%CAR Apply CAR filter to EEG data
%
% [filtered_data] = CAR(DATA)
% (time,numChannels,numEpochs)
% Nishant Mehta 2007
disp('applying CAR filter');
for z = 1:28
filtered_data(:,z,:) = data(:,z,:) - (sum(data,2) / 28);
end
@@ -0,0 +1,25 @@
% find_opt_subrotation(X)
function [theta_star, rotator_star] = find_opt_subrotation(X);
thetas = linspace(-pi/4, pi/4, 200);
max_h_sum = -Inf;
for theta = thetas
rotator =[cos(theta) -sin(theta); sin(theta) cos(theta)];
rotated_X = rotator * X;
h_sum = get_vasicek_entropy_estimate(rotated_X(1,:)) + ...
get_vasicek_entropy_estimate(rotated_X(2,:));
if h_sum > max_h_sum
max_h_sum = h_sum;
theta_star = theta;
end
end
rotator_star = ...
[cos(theta_star) -sin(theta_star); ...
sin(theta_star) cos(theta_star)];
theta_star
@@ -0,0 +1,47 @@
% find_opt_unmixing_matrix(X) - find the unmixing matrix of W that
% best separates the components of X
% X is a p x n matrix
function [Y_pos,Y_neg,W_pos,W_neg] = find_opt_unmixing_matrix(X);
p = size(X, 1);
current_X = X;
total_rotator = eye(p);
for epoch = 1:1
for i = 1:p
for j = i+1:p
disp(sprintf('[i,j] = [%d,%d]', i, j));
subspace = [current_X(i,:) ; current_X(j,:)];
[theta_star, rotator_star] = find_opt_subrotation(subspace);
new_rotator=eye(p);
new_rotator(i,i)=cos(theta_star);
new_rotator(i,j)=-sin(theta_star);
new_rotator(j,i)=sin(theta_star);
new_rotator(j,j)=cos(theta_star);
total_rotator = new_rotator * total_rotator;
current_X = total_rotator * X;
end
end
end
Y_pos = current_X;
neg_total_rotator=eye(p);
neg_total_rotator(1,1)=cos(theta_star + pi);
neg_total_rotator(1,2)=-sin(theta_star + pi);
neg_total_rotator(2,1)=sin(theta_star + pi);
neg_total_rotator(2,2)=cos(theta_star + pi);
Y_neg = neg_total_rotator * X;
W_pos = total_rotator;
W_neg = neg_total_rotator;
+71 -20
View File
@@ -15,8 +15,11 @@ rand('state', sum(100*clock))
% sample from laplacian
clear;
D = 2;
N = 10000;
p = 30;
mu = 0;
sigma = 1;
@@ -58,22 +61,22 @@ l_x = l_x - repmat(mean(l_x')', 1, N);
% generate b-spline basis curves
t = linspace(0,1,1000);
mybasis = create_bspline_basis([0 1], 30, 4);
mybasis = create_bspline_basis([0 1], p, 4);
basis_curves = eval_basis(t, mybasis);
load e1e2;
e = [e1(t); e2(t)]';
load s1s2;
s = [s1(t); s2(t)]';
data = e * l_x;
data = s * l_x;
myfd_data = data2fd(data, t, mybasis);
coef = getcoef(myfd_data);
%data1 = basis_curves * coef(:,1);
pca_results = pca_fd(myfd_data, 30);
pca_results = pca_fd(myfd_data, p);
pc_coef = getcoef(pca_results.harmfd);
pc_curves = basis_curves * pc_coef;
pc_scores = pca_results.harmscr;
@@ -83,31 +86,75 @@ pc_scores = pca_results.harmscr;
% encode our source functions e1 and e2 using the pc basis
for i=1:30
e1_weights(i) = ...
diff(ppval(fnint(spline(t, e1(t) .* pc_curves(:,i)')), ...
for i=1:p
s1_weights(i) = ...
diff(ppval(fnint(spline(t, s1(t) .* pc_curves(:,i)')), ...
[0 1]));
e2_weights(i) = ...
diff(ppval(fnint(spline(t, e2(t) .* pc_curves(:,i)')), ...
s2_weights(i) = ...
diff(ppval(fnint(spline(t, s2(t) .* pc_curves(:,i)')), ...
[0 1]));
end
for i=1:N
e1_scores(i) = dot(e1_weights, pc_scores(i,:));
e2_scores(i) = dot(e2_weights, pc_scores(i,:));
s1_scores(i) = dot(s1_weights, pc_scores(i,:));
s2_scores(i) = dot(s2_weights, pc_scores(i,:));
end
p_small = 2;
sub_pc_coef = pc_coef(:,1:p_small);
E = pc_scores(:,1:p_small)';
[Y_pos,Y_neg,W_pos,W_neg] = find_opt_unmixing_matrix(E);
for i=1:p_small
h_E(i) = get_vasicek_entropy_estimate(E(i,:));
h_Y_pos(i) = get_vasicek_entropy_estimate(Y_pos(i,:));
h_Y_neg(i) = get_vasicek_entropy_estimate(Y_neg(i,:));
end
ic_coef_pos = (W_pos * sub_pc_coef')';
ic_coef_neg = (W_neg * sub_pc_coef')';
sub_pc_curves = basis_curves * sub_pc_coef;
ic_curves_pos = basis_curves * ic_coef_pos;
ic_curves_neg = basis_curves * ic_coef_neg;
figure(1);
clf;
hold on;
plot(s, 'b');
plot(sub_pc_curves, 'r');
plot(ic_curves_pos, 'g');
plot(ic_curves_neg, 'c');
% using sub_pc_coef', recover the data
% now we want to find a matrix W that unmixes well
% let f be some candidate solution
let the data be encoded as your mother
f1_weights = rand(30,1);
f1_weights = f1_weights / norm(f1_weights);
%f1_weights = rand(p,1);
%f1_weights = f1_weights / norm(f1_weights);
f1 = pc_curves * f1_weights;
%f1 = pc_curves * f1_weights;
@@ -121,11 +168,15 @@ f1 = pc_curves * f1_weights;
% curves with splines, then use quadrature to evaluate the
% integral in [0,1]
f1_scores = zeros(N,1);
%f1_scores = zeros(N,1);
%for i=1:N
% f1_scores(i) = dot(f1_weights, pc_scores(i,:));
%end
for i=1:N
f1_scores(i) = dot(f1_weights, pc_scores(i,:));
end
+205
View File
@@ -0,0 +1,205 @@
% Given
% distribution F
% observations X (dims N x T)
% construct random vector w by sampling from F
% sort w descending
% attribute highest values of w to observations with maximal variance
% initialize random number generator
rand('state', sum(100*clock))
% sample from laplacian
clear;
D = 2;
N = 10000;
p = 30;
mu = 0;
sigma = 1;
b = sigma/2;
% laplacian
clear l_x l_px;
for i=1:D
l_x(i,:) = laplacinv(rand(N, 1), mu, b);
l_px(i,:) = laplacpdf(l_x(i,:), mu, b);
end
% center the sampling distribution
l_x = l_x - repmat(mean(l_x')', 1, N);
% gaussian
%clear g_x g_px;
%for i=1:D
% g_x(i,:) = norminv(rand(N, 1), mu, sigma);
% g_px(i,:) = normpdf(g_x(i,:), mu, sigma);
%end
% center the sampling distribution
%g_x = g_x - repmat(mean(g_x')', 1, N);
% mix the source signals
%a = rand(D,D);
%X = a * l_x;
%shuffling = shuffle(1:N);
%X_shuffled = X(:,shuffling);
% generate b-spline basis curves
t = linspace(0,1,1000);
mybasis = create_bspline_basis([0 1], p, 4);
basis_curves = eval_basis(t, mybasis);
load s1s2;
s = [s1(t); s2(t)]';
data = s * l_x;
myfd_data = data2fd(data, t, mybasis);
coef = getcoef(myfd_data);
%data1 = basis_curves * coef(:,1);
pca_results = pca_fd(myfd_data, p);
pc_coef = getcoef(pca_results.harmfd);
pc_curves = basis_curves * pc_coef;
pc_scores = pca_results.harmscr;
% encode our source functions e1 and e2 using the pc basis
for i=1:p
s1_weights(i) = ...
diff(ppval(fnint(spline(t, s1(t) .* pc_curves(:,i)')), ...
[0 1]));
s2_weights(i) = ...
diff(ppval(fnint(spline(t, s2(t) .* pc_curves(:,i)')), ...
[0 1]));
end
for i=1:N
s1_scores(i) = dot(s1_weights, pc_scores(i,:));
s2_scores(i) = dot(s2_weights, pc_scores(i,:));
end
p_small = 2;
sub_pc_coef = pc_coef(:,1:p_small);
E = pc_scores(:,1:p_small)';
[Y_pos,Y_neg,W_pos,W_neg] = find_opt_unmixing_matrix(E);
for i=1:p_small
h_E(i) = get_vasicek_entropy_estimate(E(i,:));
h_Y_pos(i) = get_vasicek_entropy_estimate(Y_pos(i,:));
h_Y_neg(i) = get_vasicek_entropy_estimate(Y_neg(i,:));
end
ic_coef_pos = (W_pos * sub_pc_coef')';
ic_coef_neg = (W_neg * sub_pc_coef')';
sub_pc_curves = basis_curves * sub_pc_coef;
ic_curves_pos = basis_curves * ic_coef_pos;
ic_curves_neg = basis_curves * ic_coef_neg;
figure(1);
clf;
hold on;
plot(s, 'b');
plot(sub_pc_curves, 'r');
plot(ic_curves_pos, 'g');
plot(ic_curves_neg, 'c');
% using sub_pc_coef', recover the data
% now we want to find a matrix W that unmixes well
% let f be some candidate solution
%f1_weights = rand(p,1);
%f1_weights = f1_weights / norm(f1_weights);
%f1 = pc_curves * f1_weights;
% we evaluate some f by considering projections P of the data onto f
% define the l2 norm for functional space:
% given some vector a and another vector b, we dot multiply the
% two vectors at the specified values, then approximate the
% curves with splines, then use quadrature to evaluate the
% integral in [0,1]
%f1_scores = zeros(N,1);
%for i=1:N
% f1_scores(i) = dot(f1_weights, pc_scores(i,:));
%end
%given the f1_scores, what to do now?
% objective function
% min sigma H(X_i)
% X i
%for a given input variable X, we seek to minimize the sum of the ...
% entropies of the marginal distributions we consider the sum
% of the entropies of the marginal distributions
% in the case of one dimension, we are given a set of scalar values
% - we can study the distribution of these values
% in the case of two dimensions, we are given a set of 2-vector
% values
% we want to know the entropy of this distribution
% the m spacing estimator studies the spacing between the sample
% points
@@ -0,0 +1,23 @@
% get_vasicek_entropy_estimate() - estimate entropy using vasicek's
% entropy estimator
function h = get_vasicek_entropy_estimate(X);
n = length(X);
m = round(sqrt(n)/2);
Z = sort(X, 'ascend');
sum_logs = 0;
for i = 1:n
if (i + m) > n
sum_logs = sum_logs + log(Z(n) - Z(i-m));
elseif (i-m) < 1
sum_logs = sum_logs + log(Z(i+m) - Z(1));
else
sum_logs = sum_logs + log(Z(i+m) - Z(i-m));
end
end
h = (sum_logs/n) + log(n/(2*m));
+1 -1
View File
@@ -5,4 +5,4 @@ if size(f) == size(g')
g = g';
end
norm = quad(@ppval, domain(1), domain(end), [], [], spline(domain, f .* g));
norm = diff(ppval(fnint(spline(domain, f .* g)), [domain(1) domain(end)]));
+4 -2
View File
@@ -5,7 +5,9 @@ c4 = 18;
mybasis = create_bspline_basis([0 .5], 30, 4);
load ../bci_comp_2003_dataSet_IV/sp1s_aa_1000Hz;
load ~/brains/bci_comp_2003_dataSet_IV/sp1s_aa_1000Hz;
x_train = CAR(x_train);
ground_truth = y_train;
@@ -71,4 +73,4 @@ right_correct = ...
abs(score_test(rights_test, pc) - lmean);
[sum(left_correct) / length(left_correct)
sum(right_correct) / length(right_correct)]
sum(right_correct) / length(right_correct)]