Skip to content

HTTPS clone URL

Subversion checkout URL

You can clone with
or
.
Download ZIP
Newer
Older
100644 142 lines (117 sloc) 4.364 kb
f3a1ee7 @hanke Imported Upstream version 8.20100721.1~dfsg.1
hanke authored
1 function [p,pc,R2] = spm_mvb_cvk(MVB,k)
2 % K-fold cross validation of a multivariate Bayesian model
3 % FORMAT [p_value,percent,R2] = spm_mvb_cvk(MVB,k)
4 %
5 % MVB - Multivariate Bayes structure
6 % k - k-fold cross-validation ('0' implies a leave-one-out scheme)
7 %
8 % p - p-value: under a null GLM
9 % percent: proportion correct (median threshold)
10 % R2 - coefficient of determination
11 %
12 % spm_mvb_cvk performs a k-fold cross-validation by trying to predict
13 % the target variable using training and test partitions on orthogonal
14 % mixtures of data (from null space of confounds)
15 %__________________________________________________________________________
16 % Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging
17
18 % Karl Friston
19 % $Id: spm_mvb_cvk.m 3806 2010-04-06 14:42:32Z ged $
20
21
22 %-partition order
23 %--------------------------------------------------------------------------
24 try
25 k;
26 catch
27 str = 'k-fold cross-validation';
28 k = spm_input(str,'!+1','b',{'2','4','8','loo'},[2 4 8 0]);
29 end
30
31 %-Get figure handles and set title
32 %--------------------------------------------------------------------------
33 Fmvb = spm_figure('GetWin','MVB');
34 spm_clf(Fmvb);
35
36 % get MVB results
37 %==========================================================================
38 try
39 MVB;
40 catch
41 mvb = spm_select(1,'mat','please select models',[],pwd,'MVB_*');
42 MVB = load(mvb(1,:));
43 MVB = MVB.MVB;
44 end
45
46 % check under null hypothesis
47 %--------------------------------------------------------------------------
48 % MVB.Y = randn(size(MVB.Y));
49
50 % whiten target and predictor (X) variables (Y) (i.e., remove correlations)
51 %--------------------------------------------------------------------------
52 K = MVB.K;
53 X = K*MVB.X;
54 Y = K*MVB.Y;
55 X0 = K*MVB.X0;
56 U = MVB.M.U;
57 Ni = length(MVB.M.F) - 1;
58
59
60 % create orthonormal projection to remove confounds
61 %--------------------------------------------------------------------------
62 Ns = length(X);
63 X0 = spm_svd(X0);
64 R = speye(Ns) - X0*X0';
65 R = spm_svd(R);
66 X = R'*X;
67 Y = R'*Y;
68 V = R'*R;
69
70 % leave-one-out
71 %--------------------------------------------------------------------------
72 if ~k, k = length(X); end
73 Ns = length(X);
74 qX = zeros(Ns,1);
75 qE = zeros(size(Y,2),k);
76 P = zeros(size(Y,2),k);
77
78 % k-fold cross-validation
79 %==========================================================================
80 for i = 1:k
81
82 % specify indices of training and test data
83 %----------------------------------------------------------------------
84 ns = floor(Ns/k);
85 test = (1:ns) + (i - 1)*ns;
86
87 % orthogonalise test and training partition
88 %----------------------------------------------------------------------
89 tran = 1:Ns;
90 tran(test) = [];
91
92 % Training
93 %======================================================================
94 M = spm_mvb(X(tran,:),Y(tran,:),[],U,[],Ni,MVB.sg);
95
96 % Test
97 %======================================================================
98 qX(test) = qX(test) + Y(test,:)*M.qE;
99
100 % record feature weights
101 %----------------------------------------------------------------------
102 qE(:,i) = M.qE;
103
104 % and posterior probabilities
105 %----------------------------------------------------------------------
106 P(:,i) = 1 - spm_Ncdf(0,abs(M.qE),M.qC);
107
108 end
109
110 % parametric inference
111 %==========================================================================
112
113 % test correlation
114 %--------------------------------------------------------------------------
115 [T df] = spm_ancova(X,V,qX,1);
116 p = 1 - spm_Tcdf(T,df(2));
117
118 % percent correct (after projection)
119 %--------------------------------------------------------------------------
120 pX = R*X;
121 qX = R*qX;
122 T = sign(pX - median(pX)) == sign(qX - median(qX));
123 pc = 100*sum(T)/length(T);
124 R2 = corrcoef(pX,qX);
125 R2 = 100*(R2(1,2)^2);
126
127 % assign in base memory
128 %--------------------------------------------------------------------------
129 MVB.p_value = p;
130 MVB.percent = pc;
131 MVB.R2 = R2;
132 MVB.cvk = struct('qX',qX,'qE',qE,'P',P);
133
134 % save results
135 %--------------------------------------------------------------------------
136 save(MVB.name,'MVB')
137 assignin('base','MVB',MVB)
138
139 % display and plot validation
140 %--------------------------------------------------------------------------
141 spm_mvb_cvk_display(MVB)
Something went wrong with that request. Please try again.