%
==================================================================================
%
==================================================================================
% Fitting function used on simulated data – may be used for real data
as well – in the paper:
%
---------------------------------------------------------------------------------------------------------------------------------------
% GP Keith, JFX DeSouza, X Yan, H Wang, JD Crawford (accepted
March 2009) A method for mapping
% response fields and determining intrinsic reference frames of
single-unit activity: Applied to 3-D
% head-unrestrained gaze shifts. Journal of Neuroscience Methods.
%
==================================================================================
% Code written by Gerald Keith in MATlab
%
==================================================================================
%
==================================================================================
%
==================================================================================
%
==================================================================================
% This function was called by a program that generated the target
position, gaze and head orientation,
% and neural activity associated with a set of trials – for
either simulated or neurophysiological
% experimental data
%
% The inputs to the function are, for all trials:
% Th = home-target position: a 3-D unit vector
pointing from the subject to the home target
% Ts = saccade-target position: a 3-D unit
vector pointing from the subject to the saccade/destination
%
target
% qGi, qGf = initial and final 3-D gaze orientation:
quaternions
% qHi, qHf = initial and final 3-D head orientation:
quaternions
% F_hold = the neuron firing rate
associated with the gaze shift: actually, number of spikes
%
==================================================================================
%
==================================================================================
function [q_PRESS_resid, best_kw] =
G_DoFits(Th,Ts,qGi,qGf,qHi,qHf,F_hold)
%
==================================================================================
%
==================================================================================
global SS;
global Tpos_4fit;
global GiD_4fit;
global F_4fit;
global kw_4fit;
global ThD_4fit;
global pref_dir_vec;
global gfT;
% ---------------------------------------------
% transform into reference frames
% ---------------------------------------------
% transform gaze (Gf) and target position (T) into the three
% reference frames: space(S), head(H), and eye(E)
% plus final gaze and target position in linearly-transformed eye frame
[pos_hold] = G_Transformations(Th,Ts,qGi,qGf,qHi,qHf);
% eliminate all trials where the saccade target was the home target
%
----------------------------------------------------------------------------------------
bad_ct = 0;
Tct_good = 0;
for Tct = 1:length(Th(:,1));
if (Th(Tct,1) == Ts(Tct,1) && Th(Tct,2) ==
Ts(Tct,2) && Th(Tct,3) == Ts(Tct,3))
bad_ct = bad_ct + 1;
else
Tct_good = Tct_good + 1;
pos(Tct_good,:,:) =
pos_hold(Tct,:,:);
F(Tct_good) = F_hold(Tct);
% calc home target direction
in space
%
-------------------------------------------------
% (convert the input 3D
target position unit-vector to a
% UV4, suitable for
manipulating with quaternions)
ThinS_UV4_hold = [0,
Th(Tct,:)];
% target direction in space
as 2D AV
ThinS_2D_AV_hold =
UV4_to_2D_AV(ThinS_UV4_hold);
ThD_4fit(Tct_good,:) =
ThinS_2D_AV_hold;
end;
end;
bad_ct
size(pos_hold)
size(pos)
NT = length(F)
ThD_vals = ThD_4fit;
% output is array: pos(trial(1:NT),dim(1:2),rep(1:8))
% =========================================
% =========================================
% do gain-field modifications and fits with gain fields
% =========================================
% =========================================
if (do_gf == 1)
% generate preferred-direction-vector, a vector that
maximizes difference
% in home-target position across all trials
max_var = 0;
for dir = 1:180;
[dir_vec(1) dir_vec(2)] =
pol2cart(dir * pi/180, 1);
cur_var =
sum(abs(dir_vec(1)*ThD_vals(:,1) + dir_vec(2)*ThD_vals(:,2)));
if (cur_var > max_var)
max_var = cur_var;
pref_dir_vec = dir_vec;
end;
end;
% =============
% do gain field fits
% =============
% -------------------------------------------------
% initial stuff
% -------------------------------------------------
% calc H directions as 2-D AVs
for cur_trial = 1:NT;
% get initial gaze quaternion
cur_qGi = qGi(cur_trial,:);
% calc initial
gaze-direction unit vector
UV04 = [0 1 0 0];
UV4_GiD =
rotateq(UV04,cur_qGi);
% get initial gaze-direction
2-D AV
cur_GiD =
UV4_to_2D_AV(UV4_GiD);
% save in initial
gaze-direction 2-D AV array
GiD(cur_trial,1) =
cur_GiD(1);
GiD(cur_trial,2) =
cur_GiD(2);
end;
TinE = pos(:,:,6);
GfinE = pos(:,:,3);
GiD_vals = GiD;
F_vals = F;
ThD_vals = ThD_4fit;
%
--------------------------------------------------------------------
% do gain-field fits for TARGET-position in eye frame
%
--------------------------------------------------------------------
% initial coefficient values
a0 = ones(3,1);
for kwi = 1:10;
Tpos_4fit = TinE;
GiD_4fit = GiD_vals;
F_4fit = F_vals;
kw_4fit = kwi;
[avals,res] =
fminsearch(@G_DoGainFieldFit,a0);
gf_TinE(:,kwi) = avals;
res_TinE(kwi) = res;
end;
figure(711);
title('Fit residuals');
hold on;
plot(1:10,res_TinE,'ks');
param_dim = 3;
figure(712)
for ai = 1:param_dim;
subplot(3,1,ai);
hold on;
title(strcat('Fit
coefficient',int2str(ai)));
plot(1:10,gf_TinE(ai,:),'ks');
end;
% plot gain field surface
% ------------------------------
XI = 0;
YI = 0;
ZI = 0;
clear XI;
clear YI;
clear ZI;
figure(720);
hold on;
title('Gi gain field COEFFICIENT for TinE');
minx = 5*floor((min(ThD_vals(:,1))-0.1)/5);
miny = 5*floor((min(ThD_vals(:,2))-0.1)/5);
maxx = 5*ceil((max(ThD_vals(:,1))+0.1)/5);
maxy = 5*ceil((max(ThD_vals(:,2))+0.1)/5);
tix = minx:1:maxy;
tiy = miny:1:maxy;
[XI,YI] = meshgrid(tix,tiy);
% use parameters for best fit kernel bandwidth
(without gain field) for representation TinE (6)
a = gf_TinE(:,best_kw(6));
gfT_coef = -(a(1) + a(2)*
(ThD_vals(:,1)*pref_dir_vec(1) + ThD_vals(:,2)*pref_dir_vec(2)) + ...
+ a(3)* (ThD_vals(:,1)*pref_dir_vec(1) +
ThD_vals(:,2)*pref_dir_vec(2)).^2);
stem3(ThD_vals(:,1),ThD_vals(:,2),gfT_coef,'ko');
for ict = 1:length(XI(:,1));
for jct = 1:length(XI(1,:));
ZI(ict,jct) = -(a(1) + a(2)* (XI(ict,jct)*pref_dir_vec(1) +
YI(ict,jct)*pref_dir_vec(2)) + ...
+ a(3)* (XI(ict,jct)*pref_dir_vec(1) +
YI(ict,jct)*pref_dir_vec(2)).^2);
end;
end;
mesh(XI,YI,ZI);
axis tight;
figure(721);
hold on;
title('Gi gain field for TinE');
gfT = exp(gfT_coef);
stem3(ThD_vals(:,1),ThD_vals(:,2),gfT,'ko');
for ict = 1:length(XI(:,1));
for jct = 1:length(XI(1,:));
ZI(ict,jct) = exp(-(a(1) + a(2)* (XI(ict,jct)*pref_dir_vec(1) +
YI(ict,jct)*pref_dir_vec(2)) + ...
+ a(3)* (XI(ict,jct)*pref_dir_vec(1) +
YI(ict,jct)*pref_dir_vec(2)).^2));
end;
end;
mesh(XI,YI,ZI);
axis tight;
%
-------------------------------------------------------------------------
% do gain-field fits for FINAL-GAZE-position in eye
frame
%
-------------------------------------------------------------------------
% initial coefficient values
a0 = ones(3,1);
for kwi = 1:10;
Tpos_4fit = GfinE;
GiD_4fit = GiD_vals;
F_4fit = F_vals;
kw_4fit = kwi;
[avals,res] =
fminsearch(@G_DoGainFieldFit,a0);
gf_GfinE(:,kwi) = avals;
res_GfinE(kwi) = res;
end;
figure(711);
title('Fit residual');
hold on;
plot(1:10,res_GfinE,'ko');
figure(712)
for ai = 1:param_dim;
subplot(3,1,ai);
hold on;
plot(1:10,gf_GfinE(ai,:),'ko');
end;
% plot gain field surface
% ------------------------------
figure(730);
hold on;
title('Gi gain field COEFFIENT for GfinE');
% use parameters for best fit kernel bandwidth
(without gain field) for representation GfinE (3)
a = gf_GfinE(:,best_kw(3));
gfG_coef = -(a(1) + a(2)*
(ThD_vals(:,1)*pref_dir_vec(1) + ThD_vals(:,2)*pref_dir_vec(2)) + ...
+ a(3)* (ThD_vals(:,1)*pref_dir_vec(1) +
ThD_vals(:,2)*pref_dir_vec(2)).^2);
stem3(ThD_vals(:,1),ThD_vals(:,2),gfG_coef,'ko');
for ict = 1:length(XI(:,1));
for jct = 1:length(XI(1,:));
ZI(ict,jct) = -(a(1) + a(2)* (XI(ict,jct)*pref_dir_vec(1) +
YI(ict,jct)*pref_dir_vec(2)) + ...
+ a(3)* (XI(ict,jct)*pref_dir_vec(1) +
YI(ict,jct)*pref_dir_vec(2)).^2);
end;
end;
mesh(XI,YI,ZI);
axis tight;
figure(731);
hold on;
title('Gi gain field for GfinE');
gfG = exp(gfG_coef);
stem3(ThD_vals(:,1),ThD_vals(:,2),gfG,'ko');
for ict = 1:length(XI(:,1));
for jct = 1:length(XI(1,:));
ZI(ict,jct) = exp(-(a(1) + a(2)* (XI(ict,jct)*pref_dir_vec(1) +
YI(ict,jct)*pref_dir_vec(2)) + ...
+ a(3)* (XI(ict,jct)*pref_dir_vec(1) +
YI(ict,jct)*pref_dir_vec(2)).^2));
end;
end;
mesh(XI,YI,ZI);
axis tight;
% -----------------------------------------
% do fit in each reference frame
% -----------------------------------------
% repeat for all representations
for rep = 1:8;
% repeat for all kernel
bandwidths
cur_best_kw = 0;
min_mean_q_PRESS_resid =
10000;
for kw = 1:15;
%
repeat for all trials
for
cur_trial = 1:NT;
cur_xpos = pos(cur_trial,1,rep);
cur_ypos = pos(cur_trial,2,rep);
numer = 0;
denom = 0;
for other_trial = 1:NT;
if (other_trial ~= cur_trial)
distsq = (cur_xpos - pos(other_trial,1,rep))^2 + (cur_ypos -
pos(other_trial,2,rep))^2;
sigsq = kw^2;
wt = exp(-distsq/sigsq);
% use either gfT (using Th position)
% or gfG (using Gi position)
numer = numer + F(other_trial) * (1/gfG(other_trial)) * wt;
denom = denom + wt;
end;
end;
fit_val = numer/denom;
PRESS_resid(cur_trial) = fit_val - F(cur_trial);
end;
%
calc q_PRESS_resid
RMS_PRESS_resid(rep,kw) = sqrt(dot(PRESS_resid,PRESS_resid)/NT);
med_PRESS_resid = median(PRESS_resid);
cur_q_PRESS_resid = abs(PRESS_resid - med_PRESS_resid);
% if
min then save kw
if
(mean(cur_q_PRESS_resid) < min_mean_q_PRESS_resid)
min_mean_q_PRESS_resid = mean(cur_q_PRESS_resid);
cur_best_kw = kw;
q_PRESS_resid_hold = cur_q_PRESS_resid;
end;
end;
best_kw(rep) = cur_best_kw;
for i =
1:length(q_PRESS_resid_hold);
q_PRESS_resid(i,rep) = q_PRESS_resid_hold(i);
end;
end;
% --------------------------------------
% plot all RMS PRESS residuals
% --------------------------------------
size(RMS_PRESS_resid)
G_PlotPRESSResids(RMS_PRESS_resid);
%
---------------------------------------------------------
% plot the best kw fits in all representations
%
---------------------------------------------------------
G_PlotFits(pos,F,best_kw);
%
----------------------------------------------------------------
% perform Levene tests on q_PRESS_resid values
%
----------------------------------------------------------------
p_vals = G_LeveneTests(q_PRESS_resid);
%
------------------------------------------------------------
% plot the p-values obtained from Levene test
%
------------------------------------------------------------
G_PlotPvals(p_vals);
%
--------------------------------------------------------------
% normalize the residuals relative to TinE
% (for use in combining data between neurons)
%
--------------------------------------------------------------
% normalization factor is mean q_PRESS_resid for
TinE representation
norm_factor = mean(q_PRESS_resid(:,4));
q_PRESS_resid = q_PRESS_resid / norm_factor;
end;
%
==================================================================================
%
==================================================================================
%
==================================================================================
%
==================================================================================
function [pos] = G_Transformations(Th,Ts,qGi,qGf,qHi,qHf)
% transforms target and final gaze directions into 2-D angle-vectors
% in different reference frames (eye, head, space)
%
==================================================================================
% outputs:
% pos(N,2,NR)
%
- positions for all N trials as 2D angle-vectors
%
in all NR representations
%
==================================================================================
%
==================================================================================
% get number of points
npts = length(Th(:,1));
% -----------------------------------------------------
% generate eye-in-head quaternion from
% gaze (eye-in-space) quaternion
% -----------------------------------------------------
for pt = 1:npts;
% Ei (initial eye-in-head)
qEi_hold = qprod(invq(qHi(pt,:)),qGi(pt,:));
qEi(pt,:) = qEi_hold;
% Ef (final eye-in-head)
qEf_hold = qprod(invq(qHf(pt,:)),qGf(pt,:));
qEf(pt,:) = qEf_hold;
end;
% ------------------------------------------
% convert quaternions to 3D AVs
% ------------------------------------------
for pt = 1:npts;
% Gi
GiAV_hold = q_to_AV(qGi(pt,:));
GiAV(pt,:) = GiAV_hold;
% Gf
GfAV_hold = q_to_AV(qGf(pt,:));
GfAV(pt,:) = GfAV_hold;
% Hi
HiAV_hold = q_to_AV(qHi(pt,:));
HiAV(pt,:) = HiAV_hold;
% Hf
HfAV_hold = q_to_AV(qHf(pt,:));
HfAV(pt,:) = HfAV_hold;
% Ei
EiAV_hold = q_to_AV(qEi(pt,:));
EiAV(pt,:) = EiAV_hold;
% Ef
EfAV_hold = q_to_AV(qEf(pt,:));
EfAV(pt,:) = EfAV_hold;
end;
%
--------------------------------------------------------------------------
% calculate unit vectors for INITIAL and FINAL
% 3D gaze direction in space, head, and eye coordinates
% AND THEN calculate associated 2D AV
%
--------------------------------------------------------------------------
for pt = 1:npts;
% final gaze in space UV
% ------------------------------
% (which is final eye-in-space)
GinS_UV4_hold = AV_to_UV4(GfAV(pt,:));
% final gaze in space 2D AV
GinS_2D_AV_hold = UV4_to_2D_AV(GinS_UV4_hold);
GfinS(pt,:) = GinS_2D_AV_hold;
% final gaze in head UV (rotate GinS UV4 by inverse
of Hi quat)
%
----------------------------------------------------------------------------------
GinH_UV4_hold =
rotateq(GinS_UV4_hold,invq(qHi(pt,:)));
% final gaze in head 2D AV
GinH_2D_AV_hold = UV4_to_2D_AV(GinH_UV4_hold);
GfinH(pt,:) = GinH_2D_AV_hold;
% final gaze in eye UV (rotate GinS UV4 by inverse
of Gi quat)
%
-------------------------------------------------------------------------------
GinE_UV4_hold =
rotateq(GinS_UV4_hold,invq(qGi(pt,:)));
% final gaze in head 2D AV
GinE_2D_AV_hold = UV4_to_2D_AV(GinE_UV4_hold);
GfinE(pt,:) = GinE_2D_AV_hold;
% INITIAL gaze in space UV
% ----------------------------------
% (which is initial eye-in-space)
GinS_UV4_hold = AV_to_UV4(GiAV(pt,:));
% final gaze in space 2D AV
GinS_2D_AV_hold = UV4_to_2D_AV(GinS_UV4_hold);
Gi_inS_2D_AV(pt,:) = GinS_2D_AV_hold;
end;
%
--------------------------------------------------------------------------
% calculate unit vectors for TARGET POSITION relative to
% subject in space, head, and eye coordinates
% AND THEN calculate associated 2D AV
%
--------------------------------------------------------------------------
for pt = 1:npts;
% target direction in space as UV
% -------------------------------------------
% (convert the input 3D target position unit-vector
to a
% UV4, suitable for manipulating with
quaternions)
TinS_UV4_hold = [0, Ts(pt,:)];
% target direction in space as 2D AV
TinS_2D_AV_hold = UV4_to_2D_AV(TinS_UV4_hold);
TinS(pt,:) = TinS_2D_AV_hold;
% target position in head as UV
% ----------------------------------------
% (rotate TinS UV4 by inverse of Hi quat)
TinH_UV4_hold =
rotateq(TinS_UV4_hold,invq(qHi(pt,:)));
% target position in head 2D AV
TinH_2D_AV_hold = UV4_to_2D_AV(TinH_UV4_hold);
TinH(pt,:) = TinH_2D_AV_hold;
% target position in eye as UV
% --------------------------------------
% (rotate TinS UV4 by inverse of Ei quat)
TinE_UV4_hold =
rotateq(TinS_UV4_hold,invq(qGi(pt,:)));
% gaze in head 2D AV
TinE_2D_AV_hold = UV4_to_2D_AV(TinE_UV4_hold);
TinE(pt,:) = TinE_2D_AV_hold;
end;
%
-------------------------------------------------------------------------------------
% calculate 2D AV positions for final gaze and target in eye frame
%
-------------------------------------------------------------------------------------
for pt = 1:npts;
% final gaze direction in eye frame as 2D AV
%
---------------------------------------------------------
GfLinE(pt,:) = GfinS(pt,:) - Gi_inS_2D_AV(pt,:);
% target position in eye frame as 2D AV
% ---------------------------------------------------
TLinE(pt,:) = TinS(pt,:) - Gi_inS_2D_AV(pt,:);
end;
% ----------------------------------------------------
% store 2D positions in the output array
% ----------------------------------------------------
% store final gaze in space positions
pos(:,1,1) = GfinS(:,1);
pos(:,2,1) = GfinS(:,2);
% store final gaze in head positions
pos(:,1,2) = GfinH(:,1);
pos(:,2,2) = GfinH(:,2);
% store final gaze in eye positions
pos(:,1,3) = GfinE(:,1);
pos(:,2,3) = GfinE(:,2);
% store target in space positions
pos(:,1,4) = TinS(:,1);
pos(:,2,4) = TinS(:,2);
% store target in head positions
pos(:,1,5) = TinH(:,1);
pos(:,2,5) = TinH(:,2);
% store target in eye positions
pos(:,1,6) = TinE(:,1);
pos(:,2,6) = TinE(:,2);
% store final gaze in linearly-transformed eye positions
pos(:,1,7) = GfLinE(:,1);
pos(:,2,7) = GfLinE(:,2);
% store target in linearly-transformed eye positions
pos(:,1,8) = TLinE(:,1);
pos(:,2,8) = TLinE(:,2);
%
==================================================================================
%
==================================================================================
%
===============================================================================
%
===============================================================================
function G_PlotPvals(RMS_PRESSResid)
%
===============================================================================
% plot the RMS PRESS residuals for all kw values
% for all representations
%
===============================================================================
%
===============================================================================
global SS;
figure(200);
hold on;
plot(1, 0.1,'wo');
if (SS == 1)
for rep = 1:8;
switch rep
case
1
plot(1:15,RMS_PRESSResid(rep,:),'bx:');
case
2
plot(1:15,RMS_PRESSResid(rep,:),'gx:');
case
3
plot(1:15,RMS_PRESSResid(rep,:),'rx:');
case
4
plot(1:15,RMS_PRESSResid(rep,:),'b+-','MarkerFaceColor','b');
case
5
plot(1:15,RMS_PRESSResid(rep,:),'g+-','MarkerFaceColor','g');
case
6
plot(1:15,RMS_PRESSResid(rep,:),'r+-','MarkerFaceColor','r');
case
7
plot(1:15,RMS_PRESSResid(rep,:),'bd:');
case
8
plot(1:15,RMS_PRESSResid(rep,:),'rd-','MarkerFaceColor','r');
end;
end;
else
for rep = 1:8;
switch rep
case
1
plot(1:15,RMS_PRESSResid(rep,:),'bo:');
case
2
plot(1:15,RMS_PRESSResid(rep,:),'go:');
case
3
plot(1:15,RMS_PRESSResid(rep,:),'ro:');
case
4
plot(1:15,RMS_PRESSResid(rep,:),'bo-','MarkerFaceColor','b');
case
5
plot(1:15,RMS_PRESSResid(rep,:),'go-','MarkerFaceColor','g');
case
6
plot(1:15,RMS_PRESSResid(rep,:),'ro-','MarkerFaceColor','r');
case
7
plot(1:15,RMS_PRESSResid(rep,:),'bs:');
case
8
plot(1:15,RMS_PRESSResid(rep,:),'rs-','MarkerFaceColor','r');
end;
end;
end;
%
===============================================================================
%
===============================================================================
%
===============================================================================
%
===============================================================================
function G_PlotFits(pos,F,best_kw)
%
===============================================================================
% inputs:
% pos(N,2,nr) - position values for
all N trials
%
in 2-D, and in all representations
%
F(N) - firing rate for all N
trials
% best_kw(nr) - kernel bandwidth of
best fit for all representations
%
===============================================================================
%
===============================================================================
global SS;
global gfT;
global ThD_4fit;
% group trials by ThD position
for ict = 1:length(ThD_4fit(:,1));
if (ThD_4fit(ict,2) > 1)
trial_class(ict) = 1;
else
if (ThD_4fit(ict,2) < -1)
trial_class(ict) = 3;
else
trial_class(ict) = 2;
end;
end;
end;
% max act
maxAct = 5*ceil(max(F)/5)
% number of trials
NT = length(pos(:,1,1));
% number of representations
NR = length(pos(1,1,:));
% limits of plot
xmin = 5*(-1 + floor(min(min(pos(:,1,:)))/5));
xmax = 5*(1 + ceil(max(max(pos(:,1,:)))/5));
ymin = 5*(-1 + floor(min(min(pos(:,2,:)))/5));
ymax = 5*(1 + ceil(max(max(pos(:,2,:)))/5));
% choose whether to use hybrid kernel bandwidth (=1)
% or fixed kernel bandwidth (=0)
hybrid_kernel = 1;
% -----------------------
% repeat for all representations
for rep = 1:NR;
% use gain field modification of firing activity
if (SS == 1)
% get activity at all grid
positions
%
-------------------------------------------
pt = 0;
for xpos = xmin:xmax;
for
ypos = ymin:ymax;
pt = pt + 1;
numer = 0;
denom = 0;
% default value of sigma-squared
sigsq = best_kw(rep)^2;
if (hybrid_kernel == 1)
% do the sigma selection
for trial = 1:NT;
dist_sq_array(trial) = (xpos - pos(trial,1,rep))^2 + (ypos -
pos(trial,2,rep))^2;
end;
closest_distsq = min(dist_sq_array);
if (closest_distsq > sigsq)
sigsq = closest_distsq;
end;
end;
for trial = 1:NT;
distsq = (xpos - pos(trial,1,rep))^2 + (ypos - pos(trial,2,rep))^2;
wt = exp(-distsq/sigsq);
numer = numer + F(trial) * (1/gfT(trial)) * wt;
denom = denom + wt;
end;
if (denom == 0)
denom = 1;
end;
fit(pt) = numer/denom;
x(pt) = xpos;
y(pt) = ypos;
end;
end;
% get activity at all trial
positions
%
-------------------------------------------
for at_trial = 1:NT;
numer = 0;
denom = 0;
%
default value of sigma-squared
sigsq = best_kw(rep)^2;
for
trial = 1:NT;
distsq = (pos(at_trial,1,rep) - pos(trial,1,rep))^2 +
(pos(at_trial,2,rep) - pos(trial,2,rep))^2;
wt = exp(-distsq/sigsq);
numer = numer + F(trial) * (1/gfT(trial)) * wt;
denom = denom + wt;
end;
if
(denom == 0)
denom = 1;
end;
fit_trial(at_trial) = numer/denom;
end;
% DO NOT use gain field modification of firing
activity
else
% get activity at all grid
positions
%
-------------------------------------------
pt = 0;
for xpos = xmin:xmax;
for
ypos = ymin:ymax;
pt = pt + 1;
numer = 0;
denom = 0;
% default value of sigma-squared
sigsq = best_kw(rep)^2;
if (hybrid_kernel == 1)
% do the sigma selection
for trial = 1:NT;
dist_sq_array(trial) = (xpos - pos(trial,1,rep))^2 + (ypos -
pos(trial,2,rep))^2;
end;
closest_distsq = min(dist_sq_array);
if (closest_distsq > sigsq)
sigsq = closest_distsq;
end;
end;
for trial = 1:NT;
distsq = (xpos - pos(trial,1,rep))^2 + (ypos - pos(trial,2,rep))^2;
wt = exp(-distsq/sigsq);
numer = numer + F(trial) * wt;
denom = denom + wt;
end;
if (denom == 0)
denom = 1;
end;
fit(pt) = numer/denom;
x(pt) = xpos;
y(pt) = ypos;
end;
end;
% get activity at all trial
positions
%
-------------------------------------------
for at_trial = 1:NT;
numer = 0;
denom = 0;
%
default value of sigma-squared
sigsq = best_kw(rep)^2;
for
trial = 1:NT;
distsq = (pos(at_trial,1,rep) - pos(trial,1,rep))^2 +
(pos(at_trial,2,rep) - pos(trial,2,rep))^2;
wt = exp(-distsq/sigsq);
numer = numer + F(trial) * wt;
denom = denom + wt;
end;
if
(denom == 0)
denom = 1;
end;
fit_trial(at_trial) = numer/denom;
end;
end;
if (SS == 1)
figure(500 + rep);
else
figure(rep);
end;
% plot fit surface with trial activations as black
circles
%
----------------------------------------------------------------------
hold on;
[X,Y] = meshgrid(xmin:xmax,ymin:ymax);
Z = griddata(x,y,fit,X,Y);
pcolor(X,Y,Z);
shading interp;
caxis([-0.2*maxAct, maxAct]);
axis equal;
for trial = 1:NT;
if (SS == 1)
ms =
1 + (F(trial)*(1/gfT(trial))*40/maxAct);
else
ms =
1 + (F(trial)*40/maxAct);
end;
plot(pos(trial,1,rep),pos(trial,2,rep),'ko','MarkerSize',ms);
end;
axis tight;
switch (rep)
case 1
title('final gaze in space frame');
case 2
title('final gaze in head frame');
case 3
title('final gaze in eye frame');
case 4
title('target in space frame');
case 5
title('target in head frame');
case 6
title('target in eye frame');
case 7
title('final gaze in linearly-translated eye frame');
case 8
title('target in linearly-translated eye frame');
end;
end;
%
===============================================================================
%
===============================================================================
%
===============================================================================
%
===============================================================================
function p_vals = G_LeveneTests(q_PRESS_resid)
%
===============================================================================
% actually, this is the Brown-Forsythe test
%
===============================================================================
% since the PRESS residuals in all representations have
% already been transformed into the absolute value
% relative to the median value, then the Levene test
% here simply is a one-tailed t-test
%
===============================================================================
%
===============================================================================
size(q_PRESS_resid)
min_mean = 1000;
best_rep = 0;
for rep = 1:8;
mean_q_PRESS_resid = mean(q_PRESS_resid(:,rep));
if (mean_q_PRESS_resid < min_mean)
min_mean =
mean_q_PRESS_resid;
best_rep = rep;
end;
end;
for rep = 1:8;
[h_val,hold_pval,ci_val] =
ttest2(q_PRESS_resid(:,rep),q_PRESS_resid(:,best_rep),0.05,1);
p_vals(rep) = hold_pval;
end;
%
===============================================================================
%
===============================================================================
%
===============================================================================
%
===============================================================================
% UTILITY FUNCTIONS
%
===============================================================================
%
===============================================================================
%
===============================================================================
function UV4 = AV_to_UV4(AV)
% convert AV to V
%
===============================================================================
UV40 = [0 1 0 0];
q = AV_to_q(AV);
UV4 = rotateq(UV40, q);
%
===============================================================================
%
===============================================================================
function iq = invq(q)
% calculate inverse of q, store as iq
%
===============================================================================
iq(1) = q(1);
iq(2) = -q(2);
iq(3) = -q(3);
iq(4) = -q(4);
%
===============================================================================
%
===============================================================================
function rq = qprod(pq, qq)
% calc product of pq and qq, store as rq
%
===============================================================================
rq(1) = pq(1)*qq(1) - pq(2)*qq(2) - pq(3)*qq(3) - pq(4)*qq(4);
rq(2) = pq(1)*qq(2) + pq(2)*qq(1) + pq(3)*qq(4) - pq(4)*qq(3);
rq(3) = pq(1)*qq(3) + pq(3)*qq(1) - pq(2)*qq(4) + pq(4)*qq(2);
rq(4) = pq(1)*qq(4) + pq(4)*qq(1) + pq(2)*qq(3) - pq(3)*qq(2);
%
===============================================================================
%
===============================================================================
function AV = q_to_AV(q)
% convert quaternion, q, to angle-vector, AV
%
===============================================================================
topi = pi/180;
th = 2 * acos(q(1))/topi;
if (th > 0)
AV(1) = q(2) * th/sin(th*topi/2.0);
AV(2) = -q(3) * th/sin(th*topi/2.0);
AV(3) = q(4) * th/sin(th*topi/2.0);
else
AV(1) = q(2)/2;
AV(2) = -q(3)/2;
AV(3) = q(4)/2;
end;
%
===============================================================================
%
===============================================================================
function q1 = rotateq(q, rq)
% rotate quaterion q by rq to get q1
%
===============================================================================
qp = qprod(rq, q);
irq = invq(rq);
q1 = qprod(qp, irq);
%
===============================================================================
%
===============================================================================
function AV = UV4_to_2D_AV(V)
% convert V to AV
%
===============================================================================
% quaternion q (in Listing's plane) for vector V
q(1) = sqrt((1 + V(2))/2);
q(2) =
0;
% qT
q(3) = - V(4) / sqrt(2*(1 + V(2))); % qV
q(4) = V(3) / sqrt(2*(1 + V(2))); % qH
topi = pi / 180;
% angle-vector AV for quaternion q
th = 2 * acos(q(1)) / topi;
AV_hold(1) = 0;
if (th == 0)
AV_hold(2) =
-q(3);
% AV V
AV_hold(3) =
q(4);
% AV H
else
AV_hold(2) = -q(3) * th / sin(th*topi/2.0);
AV_hold(3) = q(4) * th / sin(th*topi/2.0);
end;
% store as 2D AV
AV(1) = AV_hold(3); % H
AV(2) = AV_hold(2); % V
%
===============================================================================