User Tools

Site Tools


analysis:cosmo2014:module4

Differences

This shows you the differences between two versions of the page.

Link to this comparison view

Both sides previous revision Previous revision
Next revision
Previous revision
analysis:cosmo2014:module4 [2014/08/13 18:34]
mvdm
analysis:cosmo2014:module4 [2018/07/07 10:19] (current)
Line 13: Line 13:
 In the previous module, we applied the decoder to each time bin independently,​ using a flat spatial prior. In effect, this assumes that the place representation can move around arbitrarily from one time step to the next. Clearly, however, a rat cannot move around arbitrarily but instead moves subject to smoothness and continuity constraints! We can use this domain knowledge to improve the performance of our decoder. In the previous module, we applied the decoder to each time bin independently,​ using a flat spatial prior. In effect, this assumes that the place representation can move around arbitrarily from one time step to the next. Clearly, however, a rat cannot move around arbitrarily but instead moves subject to smoothness and continuity constraints! We can use this domain knowledge to improve the performance of our decoder.
  
-Our approach will be similar to Kalman filtering, in that we can construct a model of the rat's movement. We can then use this model to generate a prediction of the rat's position $P(\hat{x}_t)|P(\hat{x}_{t-1})$. The hat is meant to indicate that these are all estimates since the decoder does not have access to the rat's true position $x$.+Our approach will be similar to Kalman filtering, in that we can construct a model of the rat's movement. We can then use this model to generate a prediction of the rat's position ​($P(\hat{\mathbf{x}}_t)|P(\hat{\mathbf{x}}_{t-1}))$) which can function as a prior in our Bayesian decoder. The hat is meant to indicate that these are all estimatessince the decoder does not have access to the rat's true position $x$
 + 
 +To do this, first we need to estimate $P(\hat{\mathbf{x}}_t)|P(\hat{\mathbf{x}}_{t-1}))$ from the data. 
 + 
 +==== Estimating movement kernels ==== 
 + 
 +First, start with a clean slate: 
 + 
 +<code matlab>​ 
 +%% load the data 
 +clear all; pack 
 +cd('​C:​\Users\mvdm\Dropbox\teaching\CoSMo2014\R042-2013-08-18'​);​ % isidro 
 + 
 +load(FindFile('​*vt.mat'​));​ % from position_sandbox.m 
 +load(FindFile('​*times.mat'​));​ 
 +load(FindFile('​*CoordL.mat'​));​  
 + 
 +cfg = []; 
 +cfg.load_questionable_cells = 1; 
 +S = LoadSpikes(cfg);​ 
 + 
 +cfg = []; 
 +cfg.fc = {'​R042-2013-08-18-CSC03a.ncs'​};​ 
 +csc = LoadCSC(cfg);​ 
 +</​code>​ 
 + 
 +Next, we resample the position data at our desired binsize that matches the decoding (50ms in this case): 
 + 
 +<code matlab>​ 
 +%% resample position data at desired bin size 
 +binsize = 0.05; 
 + 
 +pos_ups = tsd; 
 +pos_ups.tvec = pos.tvec(1):​binsize:​pos.tvec(end);​ 
 +pos_ups.data(1,:​) = interp1(pos.tvec,​pos.data(1,:​),​pos_ups.tvec,'​linear'​);​ 
 +pos_ups.data(2,:​) = interp1(pos.tvec,​pos.data(2,:​),​pos_ups.tvec,'​linear'​);​ 
 +</​code>​ 
 + 
 +As before, some housekeeping to take care of: 
 + 
 +<code matlab>​ 
 +%% restrict to trials 
 +ENC_S = restrict(S,​run_start,​run_end);​ 
 +ENC_pos = restrict(pos,​run_start,​run_end);​ 
 +ENC_pos_ups = restrict(pos_ups,​run_start,​run_end);​ 
 + 
 +% check for empties and remove 
 +keep = ~cellfun(@isempty,​ENC_S.t);​ 
 +ENC_S.t = ENC_S.t(keep);​ 
 +ENC_S.label = ENC_S.label(keep);​ 
 + 
 +S.t = S.t(keep);​ 
 +S.label = S.label(keep);​ 
 +</​code>​ 
 + 
 +Now, we can obtain the distribution of position changes between adjacent time bins: 
 + 
 +<code matlab>​ 
 +clear diff_mat; 
 +diff_mat(:,​2) = diff(ENC_pos_ups.data(1,:​));​ 
 +diff_mat(:,​1) = diff(ENC_pos_ups.data(2,:​));​ 
 + 
 +dx = 0.5; xl = 10.25; 
 +x_edges = -xl:dx:xl; x_centers = x_edges(1:​end-1)+dx./​2;​ 
 +y_edges = -xl:dx:xl; y_centers = y_edges(1:​end-1)+dx./​2;​ 
 + 
 +[diff_hist,​~,​~,​pos_idx] = histcn(diff_mat,​y_edges,​x_edges);​ 
 +diff_hist = diff_hist./​sum(diff_hist(:​));​ % normalize to 1 total area 
 + 
 +subplot(221);​ 
 +imagesc(y_centers,​x_centers,​log(diff_hist));​ shading flat; colorbar 
 +hold on; 
 +plot([x_edges(1) x_edges(end)],​[0 0],'​k--'​);​ 
 +plot([0 0],​[y_edges(1) y_edges(end)],'​k--'​);​ 
 +</​code>​ 
 + 
 +{{ :​analysis:​cosmo2014:​2d_deltaspace_fine.png?​600 |}} 
 + 
 +This shows an overall tendency for the rat to travel right, and to a lesser extent up and down, which makes sense given the shape of the maze and the fact that the rat starts on the left side of the central stem. 
 + 
 +Let's construct the same at the same spatial resolution as our tuning curves: 
 + 
 +<code matlab>​ 
 +dx = 10; dx = 10; 
 +xl = 105; 
 +x_edges = -xl:dx:xl; x_centers = x_edges(1:​end-1)+dx./​2;​ 
 +y_edges = -xl:dx:xl; y_centers = y_edges(1:​end-1)+dx./​2;​ 
 + 
 +[diff_hist,​~,​~,​pos_idx] = histcn(diff_mat,​y_edges,​x_edges);​ 
 + 
 +diff_hist = diff_hist./​sum(diff_hist(:​));​ % normalize to 1 total area 
 +</​code>​ 
 + 
 +The above kernel corresponds to our model of $P(\hat{\mathbf{x}}_t)|P(\hat{\mathbf{x}}_{t-1}))$. By recursively applying this kernel we can obtain $P(\hat{\mathbf{x}}_t)|P(\hat{\mathbf{x}}_{t-n}))$:​ 
 + 
 +<code matlab>​ 
 +%% plot a few more at different speeds 
 +x1 = diff_hist;​ 
 +subplot(221);​ 
 +imagesc(y_centers,​x_centers,​log(x1));​ shading flat; colorbar 
 +hold on; 
 +plot([x_edges(1) x_edges(end)],​[0 0],'​k--'​);​ 
 +plot([0 0],​[y_edges(1) y_edges(end)],'​k--'​);​ 
 +title('​1x'​);​ 
 + 
 +x7 = diff_hist;​ 
 +for iI = 2:7 
 +    x7 = conv2(x7,​diff_hist,'​same'​);​ 
 +end 
 +subplot(222);​ 
 +imagesc(y_centers,​x_centers,​log(x7));​ shading flat; colorbar 
 +hold on; 
 +plot([x_edges(1) x_edges(end)],​[0 0],'​k--'​);​ 
 +plot([0 0],​[y_edges(1) y_edges(end)],'​k--'​);​ 
 +title('​7x'​);​ 
 + 
 +x15 = diff_hist;​ 
 +for iI = 2:15 
 +    x15 = conv2(x15,​diff_hist,'​same'​);​ 
 +end 
 +subplot(223);​ 
 +imagesc(y_centers,​x_centers,​log(x15));​ shading flat; colorbar 
 +hold on; 
 +plot([x_edges(1) x_edges(end)],​[0 0],'​k--'​);​ 
 +plot([0 0],​[y_edges(1) y_edges(end)],'​k--'​);​ 
 +title('​15x'​);​ 
 + 
 +x99 = diff_hist;​ 
 +for iI = 2:99 
 +    x99 = conv2(x99,​diff_hist,'​same'​);​ 
 +end 
 +subplot(224);​ 
 +imagesc(y_centers,​x_centers,​log(x99));​ shading flat; colorbar 
 +hold on; 
 +plot([x_edges(1) x_edges(end)],​[0 0],'​k--'​);​ 
 +plot([0 0],​[y_edges(1) y_edges(end)],'​k--'​);​ 
 +title('​99x'​);​ 
 +</​code>​ 
 + 
 +This gives: 
 + 
 +{{ :​analysis:​cosmo2014:​2d_deltaspace_all.png?​900 |}} 
 + 
 +These kernels are the expected probability distributions for different movement speeds, with 1x being the distribution actually observed in our data set. We are now ready to implement an improved version of our decoding algorithm that can use these kernels. 
 + 
 +==== Building the two-step decoder ==== 
 + 
 +As before, we first need to construct our tuning curves and Q-matrix: 
 + 
 +<code matlab>​ 
 +kernel = gausskernel([4 4],2); % 2-D gaussian for smoothing: 30 points in each direction, SD of 8 bins 
 + 
 +SET_xmin = 80; SET_ymin = 0; SET_xmax = 660; SET_ymax = 520; 
 +SET_xBinSz = 10; SET_yBinSz = 10; 
 + 
 +x_edges = SET_xmin:​SET_xBinSz:​SET_xmax;​ 
 +y_edges = SET_ymin:​SET_yBinSz:​SET_ymax;​ 
 + 
 +clear pos_mat; 
 +pos_mat(:,​2) = getd(ENC_pos,'​x'​);​ pos_mat(:,​1) = getd(ENC_pos,'​y'​);​ 
 + 
 +[occ_hist,​~,​~,​pos_idx] = histcn(pos_mat,​y_edges,​x_edges);​ 
 +no_occ_idx = find(occ_hist == 0); 
 + 
 +occ_hist = conv2(occ_hist,​kernel,'​same'​);​ 
 +occ_hist(no_occ_idx) = NaN; 
 + 
 +occ_hist = occ_hist .* (1/30); % convert to seconds using video frame rate 
 + 
 +subplot(221);​ 
 +pcolor(occ_hist);​ shading flat; axis off; colorbar 
 +title('​occupancy'​);​ 
 + 
 +
 +iC = 7; 
 +spk_x = interp1(pos.tvec,​getd(pos,'​x'​),​S.t{iC},'​linear'​);​ 
 +spk_y = interp1(pos.tvec,​getd(pos,'​y'​),​S.t{iC},'​linear'​);​ 
 + 
 +clear spk_mat; 
 +spk_mat(:,​2) = spk_x; spk_mat(:,​1) = spk_y; 
 + 
 +spk_hist = histcn(spk_mat,​y_edges,​x_edges);​ 
 +spk_hist = conv2(spk_hist,​kernel,'​same'​);​ 
 + 
 +spk_hist(no_occ_idx) = NaN; 
 + 
 +subplot(222) 
 +pcolor(spk_hist);​ shading flat; axis off; colorbar 
 +title('​spikes'​);​ 
 + 
 +
 +tc = spk_hist./​occ_hist;​ 
 + 
 +subplot(223) 
 +pcolor(tc); shading flat; axis off; colorbar 
 +title('​rate map'​);​ 
 +</​code>​ 
 + 
 +Next, the tuning curves: 
 + 
 +<code matlab>​ 
 +%% do it for all cells 
 +clear tc 
 +nCells = length(ENC_S.t);​ 
 +for iC = 1:nCells 
 +    spk_x = interp1(ENC_pos.tvec,​getd(ENC_pos,'​x'​),​ENC_S.t{iC},'​linear'​);​ 
 +    spk_y = interp1(ENC_pos.tvec,​getd(ENC_pos,'​y'​),​ENC_S.t{iC},'​linear'​);​ 
 +  
 +    clear spk_mat; 
 +    spk_mat(:,​2) = spk_x; spk_mat(:,​1) = spk_y; 
 +    spk_hist = histcn(spk_mat,​y_edges,​x_edges);​ 
 +    spk_hist = conv2(spk_hist,​kernel,'​same'​);​ 
 +     
 +    spk_hist(no_occ_idx) = NaN; 
 +  
 +    tc = spk_hist./​occ_hist;​ 
 +  
 +    all_tc{iC} = tc; 
 +  
 +end 
 + 
 +%% rehape tuning curves for decoding 
 +clear tc 
 +nBins = numel(occ_hist);​ 
 +nCells = length(ENC_S.t);​ 
 +for iC = nCells:​-1:​1 
 +    tc(:,:,iC) = all_tc{iC};​ 
 +end 
 +tc = reshape(tc,​[size(tc,​1)*size(tc,​2) size(tc,​3)]);​ 
 +occUniform = repmat(1/​nBins,​[nBins 1]); 
 + 
 +</​code>​ 
 + 
 +The Q-matrix: 
 + 
 +<code matlab>​ 
 +%% make Q-mat 
 +clear Q; 
 +binsize = 0.05; 
 + 
 +% assemble tvecs 
 +Q_tvec_edges = run_start(1):​binsize:​run_end(end);​ 
 +Q_tvec_centers = Q_tvec_edges(1:​end-1)+binsize/​2;​ 
 + 
 +for iC = length(ENC_S.t):​-1:​1 
 +  
 +    spk_t = ENC_S.t{iC};​ 
 +    Q(iC,:) = histc(spk_t,​Q_tvec_edges);​ 
 +    Q(iC,end-1) = Q(iC,​end-1)+Q(iC,​end);​ 
 +  
 +end 
 +Q = Q(:,​1:​end-1);​ 
 + 
 +Q_tsd = tsd(Q_tvec_centers,​Q);​ 
 +Q_tsd = restrict(Q_tsd,​run_start,​run_end);​ 
 + 
 +Q_tvec_centers = Q_tsd.tvec;​ 
 +Q = Q_tsd.data;​ 
 +nActiveNeurons = sum(Q > 0); 
 +</​code>​ 
 + 
 +The above steps have all been identical to those for the flat-prior decoding in the previous module. However, the actual decoding algorithm is now different:​ 
 + 
 +<code matlab>​ 
 +p1 = nan(length(Q_tvec_centers),​nBins);​ 
 +goodOccInd = find(occ_hist > 0); 
 + 
 +SET_nxBins = length(x_edges)-1;​ SET_nyBins = length(y_edges)-1;​ 
 + 
 +for iT = 1:​size(p1,​1) 
 + 
 +    L = ones(nBins,​1);​ L = L./​length(L);​ % start with uniform likelihood 
 +     
 +    if nActiveNeurons(iT) ~= 0 
 +        % obtain likelihood p(x|s_t) 
 +        for iB = 1:​length(goodOccInd);​%1:​nBins 
 +            iB = goodOccInd(iB);​ 
 +            tempProd = nansum(log(tc(iB,:​)'​.^Q(:,​iT)));​ 
 +            tempSum = exp(-binsize*nansum(tc(iB,:​),​2));​ 
 +            L(iB) = exp(tempProd)*tempSum*occUniform(iB);​ 
 +        end 
 +         
 +    end 
 +     
 +    Lu = nan(SET_nyBins,​SET_nxBins);​ % unwrap L into space for convolution 
 +    Lu(goodOccInd) = L(goodOccInd);​ 
 +     
 +    Lu = Lu./​nansum(Lu(:​));​ % normalize to 1 
 +     
 +    if iT == 1 
 +         
 +        % first sample, no movement prior 
 +        p1(iT,:) = Lu(:); 
 +        continue; 
 +         
 +    else 
 +         
 +        % obtain movement prior p(x|x_{t-1}) 
 +        temp = reshape(p1(iT-1,:​),​[SET_nyBins SET_nxBins]);​ % could optimize by carrying over to next iteration 
 +        pr = zeros(SET_nyBins,​SET_nxBins);​ 
 +        pr(goodOccInd) = temp(goodOccInd);​ 
 +        pr = conv2(pr,​x1,'​same'​);​ 
 +        pr = pr + eps; 
 +         
 +        % posterior 
 +        pst = log(Lu)+log(pr);​ 
 +        pst = exp(pst); 
 +        pst = pst./​nansum(pst(:​));​  
 + 
 +        p1(iT,:) = pst(:); 
 +         
 +    end 
 +     
 +end 
 +</​code>​ 
 + 
 +The part that computes the likelihood $P(x|s)$ is identical to the one-step case, but now there is the addition of applying our movement kernel to the previous time step's output to make the prior prediction. 
 + 
 +☛ Plot the decoding error by trial and overall as in the previous module. How does the overall error for this two-step 50ms bin size decoding compare to the version without the prior? 
 + 
 +You should find that the current implementation is much more accurate. Before declaring victory, however, let's inspect the frame-by-frame visualization -- again you can run the code from the last module (that was used to make the movie). 
 + 
 +You should find that the decoding tracks the rat's location quite well, but there is evidence for a systematic deviation: the decoding often lags behind the rat's true location. 
 + 
 +☛ We used the '​x1'​ kernel to generate our spatial prior. It is possible that the rat often runs faster than this. Try the '​x7'​ kernel, comparing the decoding error and your observations on the dynamics with the '​x1'​ version. What do you notice? What about the '​x99'​ kernel? 
 + 
 +An important technical point when using priors is implicit in the following line of the code: 
 + 
 +<code matlab>​ 
 +pr = pr + eps; 
 +</​code>​ 
 + 
 +☛ Comment out this line from the decoder and run it again. You should see that things look normal for a while but then fail catastrophically. Why? 
 + 
 +Raw decoding accuracy is an important goal in the brain-machine interface field. However, if we are interested in characterizing the dynamics of neural activity, we need to be open to the possibility that internal representations can change in ways that are not tied to actual location. We will explore this idea in the next module.
analysis/cosmo2014/module4.1407969241.txt.gz · Last modified: 2018/07/07 10:19 (external edit)