Upload files to 'hsmm'

hsmm/calcBumps.m

@@ -0,0 +1,31 @@
+function bumps = calcBumps(data)
+% it puts on each sample the correlation of that sample and the previous
+% five samples with a Bump morphology on time domain.
+% @ data: 'zscores' [nxch] a normalizad PCA from all trials concatenated
+% This function has been modified to accpet 125 Hz sampling frequency
+
+% Bump normalization sum(P) = 1.294.  
+template=[0.3090    0.8090    1.0000    0.8090    0.3090]';
+%template = [0.309016994374947;0.728968627421412;0.968583161128631;0.968583161128631;0.728968627421411;0.309016994374948];
+template=template/sum(template.^2);
+
+width = size(data,2); % number of EEG channels or PCA components
+nsamples = size(data,1); % total samples
+bumps = zeros(size(data));
+
+for j = 1:width % PCA components
+    temp=zeros(nsamples,5);
+    temp(:,1)=data(:,j);
+    for i = 2:5
+        temp(:,i)=cat(1,temp(2:end,i-1),0); % concatante along 1D
+        % puts the a component in a [nsamples X length(bump)] matrix shifted.
+        % each column is a copy of the first one but shifted one sample
+        % upwards
+    end
+    bumps(:,j)=temp*template;
+end
+
+bumps(3:end,:)=bumps(1:end-2,:); %old
+bumps([1,2,nsamples-1,nsamples],:)=0; %old
+% bumps(4:end,:)=bumps(1:end-3,:); % the edges are centered (changed)
+% bumps([1,2,3,nsamples-2,nsamples-1,nsamples],:)=0; % the leakage from the edges gets zeros (changed)

hsmm/calcEEG50CondshAllD.m

@@ -0,0 +1,89 @@
+function [likelihood eventprobs forward backward]=calcEEG50CondshAllD(bumps,conds,mags,params,x, y,map)
+spacing=5;
+offset = 2;
+%offset is how soon the first peak can be or how late the last
+%space is the distance between peaks
+nconds=max(conds);
+nstates=size(mags,2);
+lens=y-x+1;
+maxDur  = max(lens);
+ndims=size(bumps,2);
+nsamples=size(bumps,1);
+ntrials=length(x);
+gains=zeros(nsamples,nstates);
+for i = 1:ndims
+    %gains=gains+(2*bumps(:,i)*mags(i,:)-repmat(mags(i,:).^2,nsamples,1))/(2/5);
+    gains=gains+bumps(:,i)*mags(i,:)-repmat(mags(i,:).^2,nsamples,1)/2;
+end
+gains=exp(gains);
+probs=zeros(maxDur,ntrials,nstates);
+probsB=zeros(maxDur,ntrials,nstates);
+for i = 1:ntrials
+    probs(offset+1:y(i)-x(i)+1-offset,i,:)=gains(x(i)+offset:y(i)-offset,:);
+    for j = 1:nstates
+        probsB(offset+1:y(i)-x(i)+1-offset,i,j)=wrev(gains(x(i)+offset:y(i)-offset,nstates+1-j));
+    end
+end
+LP=zeros(maxDur,nstates+1,nconds);
+for j = 1:nstates+1
+    for i = 1:nconds
+        if j==nstates+1 || map(i,j)~=0 
+            % do gamma PDF for last stage and requested stage/condition (map)
+            LP(:,j,i)=gammaEEG(params(j,1,i),params(j,2,i),maxDur);
+        end
+    end
+end 
+BLP=LP(:,nstates+1:-1:1,:);
+forward = zeros(maxDur,ntrials,nstates);
+backward = zeros(maxDur,ntrials,nstates);
+for j = 1:nconds  % iterate conditions
+    % it does forward & backward by conditions. If a stage=S in a condition is not
+    % requested, it assign to stage=S forward and backward from stage=S-1
+    nj = sum(conds==j);  % number of trials of condition
+    lensT=lens(conds==j);
+    forwardT=zeros(maxDur,nj,nstates);
+    forwardBT=zeros(maxDur,nj,nstates);
+    backwardT=zeros(maxDur,nj,nstates);
+    forwardBT(offset+1:maxDur,:,1)=repmat(BLP(1:maxDur-offset,1,j),1,nj);  % reversed gamma PDF
+    forwardT(offset+1:maxDur,:,1)=repmat(LP(1:maxDur-offset,1,j),1,nj).*probs(offset+1:maxDur,conds==j,1); % Gamma PDF * Gains of condition
+    probsT=probs(:,conds==j,:);
+    probsBT=probsB(:,conds==j,:);
+    for i = 2:nstates   % iterate bumps
+        if map(j,i)>0  % if condition/current stage estimation requested
+            next=cat(1,zeros(spacing,1),LP(1:maxDur-spacing,i,j)); % gamma PDF (shifted by offset)
+            for k = 1:nj % iterate trials of condition
+                temp=conv(forwardT(:,k,i-1),next); % convolution with previous stage ~~ conv(gammaPDF(i-1)*gains(i-1), gammaPDF(i))
+                forwardT(:,k,i)=temp(1:maxDur); 
+            end
+            forwardT(:,:,i)=forwardT(:,:,i).*probsT(:,:,i);  % forwardT ~~ conv(gammaPDF(i-1)*gains(i-1), gammaPDF(i)) * gains(i)
+        else % if not requested: assign to current stage forwardT of previous stage/state
+            forwardT(:,:,i)=forwardT(:,:,i-1);
+        end
+        if map(j,nstates+2-i)>0 % if condition/last minus current stage estimation requested 
+            nextB=cat(1,zeros(spacing,1),BLP(1:maxDur-spacing,i,j));
+            addB=forwardBT(:,:,i-1).*probsBT(:,:,i-1); % add ~~ gammaPDF_b(i-1)*gains_b(i-1)     ; _b = reverted back to forth
+            for k=1:nj % iterate trials of condition
+                temp=conv(addB(:,k),nextB); % temp ~~ conv(gammaPDF_b(i-1)*gains_b(i-1), gammaPDF_b(i))
+                forwardBT(:,k,i)=temp(1:maxDur);  
+            end
+        else % if not requested: assign to current stage forwardBT of previous stage/state
+            forwardBT(:,:,i)=forwardBT(:,:,i-1);  
+        end
+    end
+    forwardBT=forwardBT(:,:,nstates:-1:1);
+    for k=1:nj
+        for i = 1:nstates
+            backwardT(1:lensT(k),k,i)=wrev(forwardBT(1:lensT(k),k,i));  % revers forwardBT by trial and stage
+        end
+    end
+    backwardT(1:offset,:,:)=0;
+    backward(:,conds==j,:)=backwardT;
+    forward(:,conds==j,:)=forwardT;
+end
+temp=forward.*backward;
+likelihood=sum(log(sum(temp(:,:,1))));
+eventprobs=temp./repmat(sum(temp),[maxDur,1,1]);
+
+
+
+

hsmm/gammaEEG.m

@@ -0,0 +1,8 @@
+function d = gammaEEG(a,b,maxlength)
+% returns gamma probability density function normalized to one
+d=zeros(maxlength,1);
+for i = 1:maxlength
+    d(i)=gampdf(i-.5,a,b);
+end
+d = d/sum(d);
+

hsmm/gammaParamsCondshD.m

@@ -0,0 +1,42 @@
+function [params events temp] = gammaParamsCondshD(eventprobs,lens,shape,conds,map,maxDur)
+% @ lens: mean trial duration per condition [N Conditions]
+offset1=3;
+spacing=5;
+nconds=max(conds);
+map=cat(2,map,ones(nconds,1)); % add one column maps[conditions, N stages]
+nstates=size(eventprobs,3);  % number of bumps
+temp=zeros(maxDur,nstates,nconds);
+params=zeros(nstates+1,2,nconds);
+for i = 1:nconds
+    temp(:,:,i)=reshape(mean(eventprobs(:,conds==i,:),2),maxDur,nstates);
+    % temp[samples,bumps,condition] mean accros trials per sample, like an ERP
+end
+maxobs=(maxDur+1)-(2*(offset1-1)+(nstates-1)*spacing); % max duration of trial without bumps' length and offset
+events=zeros(nstates+2,nconds);
+for j = 1:nconds
+    offset=-2;
+    for i = 1:nstates
+        if map(j,i)~=0 % if condition/current stage estimation requested
+            offset=offset+spacing;
+            temp(1:maxobs,i,j)=temp(offset:offset+maxobs-1,i,j);
+            % global expected location of a bump
+            events(i+1,j)=[1:maxobs]*temp(1:maxobs,i,j);
+        else % if not requested: assign to next stage 'events' of current stage/state
+            events(i+1,j)=events(i,j);
+        end
+    end
+    events(i+2,j)=lens(j)-offset-1; % mean trial duration per condition - offest - 1
+end
+temp(maxobs+1:end,:,:)=[];
+for j = 1:nstates+1 % iterate stages
+    % get the index of requested stages to compute the duration of the flats.
+    uniques=unique(map(:,j)); % take uinque requested map/condition values in bump
+    uniques=uniques(uniques>0); % discard zeros (not requested)
+    for i = 1:length(uniques) % iteratee unique values
+        subs=find(map(:,j)==i); % index of unque value in a stage
+        % Gamma scale : duration of a flat (event(j+1) - event(j)) divided by shape
+        params(j,:,subs)=repmat([shape,mean(events(j+1,subs)-events(j,subs))/shape],[1,1,length(subs)]);
+    end
+end
+params(1,2,:)=params(1,2,:)-.5/shape;
+params(2:nstates,2,:)=params(2:nstates,2,:)+.5/shape;

hsmm/gammaParamsN.m

@@ -0,0 +1,18 @@
+function params = gammaParamsN(eventprobs,lens,shape)
+%4 samples are not accounted for in the flats
+width=4;
+N=size(eventprobs,1);
+% given that the shape is fixed the calculation of the maximum likelihood
+% scales becomes simple.  One just calculates the means expected lengths of
+% the flats and divides by the shape=
+nstates = size(eventprobs,3);
+averagepos=cat(2,sum(repmat([1:N]',1,nstates).*reshape(mean(eventprobs,2),N,nstates)),mean(lens));
+averagepos=averagepos-(2+[0:width:(nstates-1)*width,(nstates-1)*width+1]);
+flats=averagepos-cat(2,0,averagepos(1:end-1));
+params=zeros(nstates+1,2);
+params(:,1)=shape;
+params(:,2)=flats'/shape;
+%correct flats between bumps for the fact that the gamma is calculated at
+%midpoint
+params(:,2)=params(:,2)- .5/shape;
+

hsmm/hsmmEEG50CondMapAllDNew.m

@@ -0,0 +1,62 @@
+function [lkh1 mags1 params1 eventprobs1]=hsmmEEG50CondAllDNew(bumps,conds,mags,params,thresh,x,y,map,shape)
+% @ conds: [N trials] condition per trial as integres from 1 to N conditions 
+% @ map: [conditions, stages?]
+
+lens    = y - x + 1;
+maxDur  = max(lens);
+
+% bumps=calcBumps(bumps); % do it outside if commented
+nconds=max(conds);
+if length(size(params))==2;
+    params=repmat(params,[1,1,nconds]); % params [N stages, (shape,scale), N conds]
+end
+lkh1=-Inf;
+nstates=size(mags,2);
+ncases=length(x);
+ndims=size(bumps,2);
+means=zeros(maxDur,ncases,ndims);
+lengths=zeros(nconds,1);
+for i = 1:nconds
+    lengths(i)=mean(lens(conds==i)); % mean condition duration. It's used to estimate Gamma parameteres 
+end
+for i = 1:ncases
+    means(1:lens(i),i,:)=bumps(x(i):y(i),:); % mean amplitude
+end
+[lkh, eventprobs,~,~]=calcEEG50CondshAllD(bumps,conds,mags,params,x,y,map);
+%lkh
+if thresh==0
+    lkh1=lkh;
+    mags1 = mags;
+    params1 =params;
+    eventprobs1=eventprobs;
+else
+    while (lkh-lkh1)>thresh
+        lkh1=lkh;
+        mags1 = mags;
+        params1 =params;
+        eventprobs1=eventprobs;
+        for i = 1:nstates
+            a=find(map(conds,i)>0);
+            for j = 1:ndims
+                mags(j,i)=mean(sum(eventprobs(:,a,i).*means(:,a,j)));
+            end
+        end
+        params=gammaParamsCondshD(eventprobs,lengths,shape,conds,map,maxDur);
+        for i = 1:nstates+1
+            for j = 1:nconds
+                if prod(params(i,:,j))< 5
+                    params(i,:,j)=params1(i,:,j);
+                end
+            end
+        end
+        [lkh, eventprobs]=calcEEG50CondshAllD(bumps,conds,mags,params,x,y,map);
+%         reshape(params(:,2,:),nstates+1,nconds)'
+        %lkh-lkh1
+        %lkh
+    end
+    eventprobs1 = single(eventprobs1);
+    mags1 = single(mags1);
+    params1 = single(params1);
+end
+    
+