Skip to content

Commit

Permalink
updating
Browse files Browse the repository at this point in the history
  • Loading branch information
scott committed Mar 24, 2017
1 parent 5252665 commit 2119019
Show file tree
Hide file tree
Showing 5 changed files with 288 additions and 19 deletions.
234 changes: 234 additions & 0 deletions quantile.m
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,234 @@
function [q,N] = quantile(X,p,dim,method,weights)
%QUANTILE Quantiles of a sample via various methods.
%
% Q = IOSR.STATISTICS.QUANTILE(X,P) returns quantiles of the values in X.
% P is a scalar or a vector of cumulative probability values. When X is
% a vector, Q is the same size as P, and Q(i) contains the P(i)-th
% quantile. When X is a matrix, the i-th row of Q contains the P(i)-th
% quantiles of each column of X. For N-D arrays,
% IOSR.STATISTICS.QUANTILE operates along the first non-singleton
% dimension.
%
% Q = IOSR.STATISTICS.QUANTILE(X,P,DIM) calculates quantiles along
% dimension DIM. The DIM'th dimension of Q has length LENGTH(P).
%
% Q = IOSR.STATISTICS.QUANTILE(X,P,DIM,METHOD) calculates quantiles using
% one of the methods described in http://en.wikipedia.org/wiki/Quantile.
% The method are designated 'R-1'...'R-9'; the default is R-8 as
% described in http://bit.ly/1kX4NcT, whereas Matlab uses 'R-5'.
%
% Q = IOSR.STATISTICS.QUANTILE(X,P,[],METHOD) uses the specified METHOD,
% but calculates quantiles along the first non-singleton dimension.
%
% Q = IOSR.STATISTICS.QUANTILE(X,P,[],METHOD,WEIGHTS) and
% IOSR.STATISTICS.QUANTILE(X,P,[],[],WEIGHTS) uses the array WEIGHTS to
% weight the values in X when calculating quantiles. If no weighting is
% specified, the method determines the real-valued index in to the data
% that is used to calculate the P(i)-th quantile. When a weighting array
% WEIGHTS is specified (WEIGHTS should be the same size as X), this index
% is mapped to the cumulative weights (the weights are scaled to sum to
% N(i) - see below), and a new weighted index is returned (using linear
% interpolation) for the point where the cumulative weights equal the
% unweighted index. The weighted index is used to calculate the P(i)-th
% quantile. If the values in WEIGHTS are equal, then the weighted and
% unweighted index (and correpsonding quantile) are identical. The
% default method R-8 is used if METHOD is specified as an empty array
% ([]).
%
% [Q,N] = IOSR.STATISTICS.QUANTILE(...) returns an array that is the same
% size as Q such that N(i) is the number of points used to calculate
% Q(i).
%
% Further reading
%
% Hyndman, R.J.; Fan, Y. (November 1996). "Sample Quantiles in
% Statistical Packages". The American Statistician 50 (4): 361-365.
% Frigge, Michael; Hoaglin, David C.; Iglewicz, Boris (February 1989).
% "Some Implementations of the Boxplot". The American Statistician 43
% (1): 50-54.
%
% See also QUANTILE.

% Copyright 2016 University of Surrey.

%% Check input and make default assignments

assert(isnumeric(X),'X must be a numeric');
assert(isvector(p) & isnumeric(p),'P must be a numeric vector');
assert(all(p>=0 & p<=1),'Values in P must be in the interval [0,1].')

if nargin<2
error('Not enough input arguments.')
end

dims = size(X);
if nargin<3 || isempty(dim)
dim = find(dims>1,1,'first'); % default dim
else % validate input
assert(isnumeric(dim) | isempty(dim),'DIM must be an integer or empty');
assert(isint(dim) | isempty(dim),'DIM must be an integer or empty');
assert(dim>0,'DIM must be greater than 0')
end

if nargin<4
method = 'r-8'; % default method
else % validate input
if isempty(method)
method = 'r-8'; % default method
else
assert(ischar(method),'METHOD must be a character array')
end
end

if nargin<5
weights = [];
else
assert(isequal(size(X),size(weights)) || isempty(weights),'WEIGHTS must be the same size as X');
end

%% choose method

% See http://en.wikipedia.org/wiki/Quantile#Estimating_the_quantiles_of_a_population

switch lower(method)
case 'r-1'
min_con = @(N,p)(p==0);
max_con = @(N,p)(false);
h = @(N,p)((N*p)+.5);
Qp = @(x,h)(x(ceil(h-.5)));
case 'r-2'
min_con = @(N,p)(p==0);
max_con = @(N,p)(p==1);
h = @(N,p)((N*p)+.5);
Qp = @(x,h)((x(ceil(h-.5))+x(floor(h+.5)))/2);
case 'r-3'
min_con = @(N,p)(p<=(.5/N));
max_con = @(N,p)(false);
h = @(N,p)(N*p);
Qp = @(x,h)(x(round(h)));
case 'r-4'
min_con = @(N,p)(p<(1/N));
max_con = @(N,p)(p==1);
h = @(N,p)(N*p);
Qp = @(x,h)(x(floor(h)) + ((h-floor(h))*(x(floor(h)+1)-x(floor(h)))));
case 'r-5'
min_con = @(N,p)(p<(.5/N));
max_con = @(N,p)(p>=((N-.5)/N));
h = @(N,p)((N*p)+.5);
Qp = @(x,h)(x(floor(h)) + ((h-floor(h))*(x(floor(h)+1)-x(floor(h)))));
case 'r-6'
min_con = @(N,p)(p<(1/(N+1)));
max_con = @(N,p)(p>=(N/(N+1)));
h = @(N,p)((N+1)*p);
Qp = @(x,h)(x(floor(h)) + ((h-floor(h))*(x(floor(h)+1)-x(floor(h)))));
case 'r-7'
min_con = @(N,p)(false);
max_con = @(N,p)(p==1);
h = @(N,p)(((N-1)*p)+1);
Qp = @(x,h)(x(floor(h)) + ((h-floor(h))*(x(floor(h)+1)-x(floor(h)))));
case 'r-8'
min_con = @(N,p)(p<((2/3)/(N+(1/3))));
max_con = @(N,p)(p>=((N-(1/3))/(N+(1/3))));
h = @(N,p)(((N+(1/3))*p)+(1/3));
Qp = @(x,h)(x(floor(h)) + ((h-floor(h))*(x(floor(h)+1)-x(floor(h)))));
case 'r-9'
min_con = @(N,p)(p<((5/8)/(N+.25)));
max_con = @(N,p)(p>=((N-(3/8))/(N+.25)));
h = @(N,p)(((N+.25)*p)+(3/8));
Qp = @(x,h)(x(floor(h)) + ((h-floor(h))*(x(floor(h)+1)-x(floor(h)))));
otherwise
error(['Method ''' method ''' does not exist'])
end

%% calculate quartiles

% reshape data so function works down columns
order = mod(dim-1:dim+length(dims)-2,length(dims))+1;
dims_shift = dims(order);
x = rearrange(X,order,[dims_shift(1) prod(dims_shift(2:end))]);
if ~isempty(weights)
weights = rearrange(weights,order,[dims_shift(1) prod(dims_shift(2:end))]);
cumwfunc = @accumulateWeights;
wfunc = @weightedIndex;
else
cumwfunc = @(~,~,~,N) 1:N;
wfunc = @(x,~) x;
end

% pre-allocate q
q = zeros([length(p) prod(dims_shift(2:end))]);
N = zeros([length(p) prod(dims_shift(2:end))]);
for m = 1:length(p)
for n = 1:numel(q)/length(p)
[xSorted,ind] = sort(x(~isnan(x(:,n)),n)); % sort
N(m,n) = length(xSorted); % sample size
k = cumwfunc(weights,ind,n,N(m,n));
switch N(m,n)
case 0
q(m,n) = NaN;
case 1
q(m,n) = xSorted;
otherwise
if min_con(N(m,n),p(m)) % at lower limit
q(m,n) = xSorted(1);
elseif max_con(N(m,n),p(m)) % at upper limit
q(m,n) = xSorted(N(m,n));
else % everything else
huw = h(N(m,n),p(m)); % unweighted index
hw = wfunc(huw,k);
q(m,n) = Qp(xSorted,hw);
end
end
end
end

% restore dims of q to equate to those of input
q = irearrange(q,order,[length(p) dims_shift(2:end)]);
N = irearrange(N,order,[length(p) dims_shift(2:end)]);

% if q is a vector, make same shape as p
if numel(p)==numel(q)
q=reshape(q,size(p));
N=reshape(N,size(p));
end

end

function cumweights = accumulateWeights(weights,ind,n,N)
%ACCUMULATEWEIGHTS accumulate the weights

wSorted = weights(ind,n); % sort weights
wSorted = wSorted*N/sum(wSorted); % normalize weights to sum to N
cumweights = cumsum(wSorted); % cumulative weights

end

function hw = weightedIndex(huw, cumweights)
%WEIGHTEDINDEX calculate index from cumulative weights

ii = find(sign(cumweights-huw)<0,1,'last');
jj = find(sign(cumweights-huw)>0,1,'first');
if isempty(ii) || isempty(jj)
hw = huw;
else
hw = ii + (huw-cumweights(ii))/(cumweights(jj)-cumweights(ii)); % weighted index
end

end

function y = isint(x)
%ISINT check if input is whole number
y = x==round(x);
end

function y = rearrange(x,order,shape)
%REARRANGE reshape and permute to make target dim column
y = permute(x,order);
y = reshape(y,shape);
end

function y = irearrange(x,order,shape)
%IREARRANGE reshape and permute to original size
y = reshape(x,shape);
y = ipermute(y,order);
end
Binary file modified svca4_PlotBPvoiGui.fig
View file
Binary file not shown.
14 changes: 4 additions & 10 deletions svca4_extractRefGui.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -46,20 +46,13 @@

% --- Executes just before svca4_extractRefGui is made visible.
function svca4_extractRefGui_OpeningFcn(hObject, eventdata, handles, varargin)
% This function has no output args, see OutputFcn.
% hObject handle to figure
% eventdata reserved - to be defined in a future version of MATLAB
% handles structure with handles and user data (see GUIDATA)
% varargin command line arguments to svca4_extractRefGui (see VARARGIN)

% Choose default command line output for svca4_extractRefGui
handles.output = hObject;

% Update handles structure
guidata(hObject, handles);

% UIWAIT makes svca4_extractRefGui wait for user response (see UIRESUME)
% uiwait(handles.figure1);
global svca4
handles.listsubs.String = svca4.PET_list;
handles.listsubs.Max = length(svca4.PET_list);
Expand All @@ -70,7 +63,6 @@ function svca4_extractRefGui_OpeningFcn(hObject, eventdata, handles, varargin)
% Get default command line output from handles structure
varargout{1} = handles.output;


% --- Executes on selection change in listsubs.
function listsubs_Callback(hObject, eventdata, handles)
global svca4
Expand Down Expand Up @@ -104,7 +96,7 @@ function extract_Callback(hObject, eventdata, handles)
GRAYtac = zeros(1,length(svca4.PET_standardDurations));

for s = 1:length(handles.listsubs.Value)
ifeedback=str2num(handles.iteration.String); % this can be changed later
ifeedback=str2num(handles.iteration.String);
q = str2num(handles.it_q.String);

% load the brain mask
Expand Down Expand Up @@ -181,7 +173,7 @@ function extract_Callback(hObject, eventdata, handles)
fprintf('BLOOD selection : %d voxels remaining\n', numel(indGRAY));
indGRAY = intersect(indGRAY, find(TSPO.*MASK<GRAY_q(4)));
fprintf('TSPO selection : %d voxels remaining\n', numel(indGRAY));
%

%----------------------------------------------
% - Calculate SVCA reference curve
%----------------------------------------------
Expand Down Expand Up @@ -258,6 +250,8 @@ function extract_Callback(hObject, eventdata, handles)
figure;
plot(svca4.PET_standardEndTimes,mean(GRAYtac))
elseif handles.save_mean.Value == 1 && handles.remCereb.Value == 1
myGRAY_TAC = [svca4.PET_standardStartTimes svca4.PET_standardEndTimes mean(GRAYtac)'];

if ifeedback == 0
fname = sprintf('%s/TACs/mean_noCB_svcaRef_G%.2fW%.2fB%.2fT%.2f.txt', svca4.outputPath, svca4.quantiles);
else fname = sprintf('%s/TACs/mean_noCB_svcaRef_q%d_it%.2d_G%.2fW%.2fB%.2fT%.2f.txt', svca4.outputPath, q*100,ifeedback,svca4.quantiles);
Expand Down
49 changes: 49 additions & 0 deletions svca4_lobeBP.m
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,49 @@
function [BP] = svca4_lobeBP(bpTable)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% gives the BP for each lobe, i.e., the average BP from all regions within
% the lobe.
% regions in lobes taken from: https://surfer.nmr.mgh.harvard.edu/fswiki/CorticalParcellation
% Desikan-Killiany ROIs
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

frontal.name = {'ctx_lh_superiorfrontal' 'ctx_rh_superiorfrontal'...
'ctx_lh_rostralmiddlefrontal' 'ctx_rh_rostralmiddlefrontal'...
'ctx_lh_caudalmiddlefrontal' 'ctx_rh_caudalmiddlefrontal'...
'ctx_lh_parsopercularis' 'ctx_rh_parsopercularis'...
'ctx_lh_parstriangularis' 'ctx_rh_parstriangularis'...
'ctx_lh_parsorbitalis' 'ctx_rh_parsorbitalis'...
'ctx_lh_lateralorbitofrontal' 'ctx_rh_lateralorbitofrontal'...
'ctx_lh_medialorbitofrontal' 'ctx_rh_medialorbitofrontal'...
'ctx_lh_precentral' 'ctx_rh_precentral'...
'ctx_lh_paracentral' 'ctx_rh_paracentral'...
'ctx_lh_frontalpole' 'ctx_rh_frontalpole'};
frontal.num = [69 104 68 103 45 80 59 94 61 96 60 95 53 88 55 90 65 100 58 93 73 108];
BP.frontal = mean(bpTable{:,frontal.num},2);

parietal.name = {'ctx_lh_superiorparietal' 'ctx_rh_superiorparietal'...
'ctx_lh_inferiorparietal' 'ctx_rh_inferiorparietal'...
'ctx_lh_supramarginal' 'ctx_rh_supramarginal'...
'ctx_lh_postcentral' 'ctx_rh_postcentral'...
'ctx_lh_precuneus' 'ctx_rh_precuneus'};
parietal.num = [70 105 49 84 72 107 63 98 66 101];
BP.parietal = mean(bpTable{:,parietal.num},2);


temporal.name = {'ctx_lh_superiortemporal' 'ctx_rh_superiortemporal'...
'ctx_lh_middletemporal' 'ctx_rh_middletemporal'...
'ctx_lh_inferiortemporal' 'ctx_rh_inferiortemporal'...
'ctx_lh_bankssts' 'ctx_rh_bankssts'...
'ctx_lh_fusiform' 'ctx_rh_fusiform'...
'ctx_lh_transversetemporal' 'ctx_rh_transversetemporal'...
'ctx_lh_entorhinal' 'ctx_rh_entorhinal'...
'ctx_lh_temporalpole' 'ctx_rh_temporalpole'...
'ctx_lh_parahippocampal' 'ctx_rh_parahippocampal'};
temporal.num = [71 106 56 91 50 85 43 78 48 83 75 110 47 82 74 109 57 92];
BP.temporal = mean(bpTable{:,temporal.num},2);

occipital.name = {'ctx_lh_lateraloccipital' 'ctx_rh_lateraloccipital'...
'ctx_lh_lingual' 'ctx_rh_lingual'...
'ctx_lh_cuneus' 'ctx_rh_cuneus'...
'ctx_lh_pericalcarine' 'ctx_rh_pericalcarine'};
occipital.num = [52 87 54 89 46 81 62 97];
BP.occipital = mean(bpTable{:,occipital.num},2);
10 changes: 1 addition & 9 deletions svca4_mainGui.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -89,9 +89,7 @@ function load_o_Callback(hObject, eventdata, handles)

% --- Executes on button press in create_svca4.
function create_svca4_Callback(hObject, eventdata, handles)
% hObject handle to create_svca4 (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB
% handles structure with handles and user data (see GUIDATA)

svca4_createSVCAGui

% --- Executes on button press in show_svca4.
Expand All @@ -112,8 +110,6 @@ function time_shift_Callback(hObject, eventdata, handles)

% --- Executes on button press in freesurfer.
function freesurfer_Callback(hObject, eventdata, handles)
% NB: Placeholder: use this button to open a gui that collects the
% SUBJECTS_DIR, the subject and runs a recon-all
global svca4
if exist('svca4','var')
svca4_freesurferGui(svca4)
Expand All @@ -139,7 +135,6 @@ function plot_class_Callback(hObject, eventdata, handles)
else warndlg('Please load a TAC_TABLE and o structure first')
end


% --- Executes on button press in calc_svca4.
function calc_svca4_Callback(hObject, eventdata, handles)
global svca4
Expand Down Expand Up @@ -274,7 +269,6 @@ function compareTacs_Callback(hObject, eventdata, handles)
svca4_compareTacsGui(svca4)
end


% --- Executes on button press in bp_voi.
function bp_voi_Callback(hObject, eventdata, handles)
% This call will create a table containing the average BP value for
Expand Down Expand Up @@ -356,11 +350,9 @@ function copy_header_Callback(hObject, eventdata, handles)
end
end


% --- Executes on button press in corr_bp.
function corr_bp_Callback(hObject, eventdata, handles)


% --- Executes on button press in plot_BP.
function plot_BP_Callback(hObject, eventdata, handles)
global svca4
Expand Down

0 comments on commit 2119019

Please sign in to comment.