Merge branch 'develop'

jdiedrichsen · Apr 27, 2018 · 8ee1caa · 8ee1caa
2 parents d11ebe6 + 3de3f76
commit 8ee1caa
Show file tree

Hide file tree

Showing 10 changed files with 321 additions and 185 deletions.
diff --git a/documentation/.DS_Store b/documentation/.DS_Store
diff --git a/documentation/pcm_toolbox_manual.md b/documentation/pcm_toolbox_manual.md
@@ -72,10 +72,12 @@ These functions are higher level functions that perform fitting and crossvalidat
 | 	Function 			       | Comment  
 |:-----------------------------|:-----------------------------
 | `pcm_fitModelIndivid`        | Fits G and noise parameter to individual data
+| `pcm_fitModelIndividCrossval`| Within-subject crossvalidation of models
 | `pcm_fitModelGroup`          | Fit common G to all subjects
 | `pcm_fitModelGroupCrossval`  | Between-subject crossvalidation of models  
 | `pcm_setUpFit`         |  Generally prepares models and data for fitting 
-| `pcm_componentPosterior`     | Inference on components of a Model family 
+| `pcm_knockModels`     | Inference on components using simple knock-in knock-out likelihoods
+| `pcm_componentPosterior`     | Inference on components by model averaging 
 
 ### Visualization functions
 | 	Function 			       | Comment  
@@ -124,7 +126,6 @@ These functions are higher level functions that perform fitting and crossvalidat
 |  `pcm_NR_diag`			| Newton Raphson for diagonalized component models
 |  `pcm_NR_free`			| Newton Raphson for a free model 
 |  `pcm_EM`				| Expectation-Maximization 
-| `pcm_fitModelIndividCrossval`| Within-subject crossvalidation of models 
 
 
 
@@ -413,7 +414,23 @@ r_model1    = (cov12./sqrt(var1.*var2))';
 
 ## Fitting to individual data sets with cross-validation across partitions 
 
-Crossvalidation within subject ([@fig:Fig3]a) is the standard for encoding models and can also be applied to PCM-models.  Note that this part of the toolbox is currently under development. 
+Crossvalidation within subject ([@fig:Fig3]a) is the standard for encoding models and can also be applied to PCM-models.  Note that this part of the toolbox is still under development.
+
+```
+function [T,D,theta_hat]=pcm_fitModelIndividCrossval(Y,M,partitionVec,conditionVec,varargin);
+```
+
+The use is very similar `pcm_fitModelIndivid`, except for two optional input parameters: 
+```
+'crossvalScheme': Crossvalidation scheme on the different partitions 
+                  'leaveOneOut': Leave one partition at a time out 
+                  'leaveTwoOut': Leave two consecutive partitions out
+                  'oddeven': Split half by odd and even runs 
+'evaluation':  So far implemented are only the standard encoding-model 
+               style evalutation criteria.
+                  'R2': Crossvalidated R^2 
+                  'R' : Correlation between observed and predicted  
+```
 
 ##  Fitting to group data sets
 

diff --git a/documentation/pcm_toolbox_manual.pdf b/documentation/pcm_toolbox_manual.pdf
diff --git a/documentation/readme.md b/documentation/readme.md
@@ -11,49 +11,75 @@ What the toolbox does *not* provide are functions to extract the required data f
 
 The functions provided in the toolbox can be categorized into different categories:  
 
-## Basic likelihood and optimization
+### Basic likelihood and optimization
 
 These are the functions that perform the core statistical functions of the toolbox.  
 
 | 	Function 			    | Comment  
-|:--------------------------|:-----------------------------
-|	`pcm_likelihoodIndivid` | Likelihood of a single data set under a model
-|   `pcm_likelihoodGroup`   | Likelihood of a group data set under a model
-|   `pcm_NR`				| Newton Raphson optimisation 
-|   `pcm_NR_diag`			| Newton Raphson for diagonalized models (faster)
-|   `pcm_NR_free`			| Newton Raphson for a free model 
-|   `pcm_EM`				| Expectation-Maximization 
-|   `pcm_minimize`			| Conjugate gradient descent 
-
-## Model Evaluation
+|:-----------------------------------|:-----------------------------
+|  `pcm_likelihood` |  Likelihood of a single data set under a model
+|  `pcm_likelihoodIndivid` | pcm_likelihood with optional random or fixed block effect
+|  `pcm_likelihoodGroup`   | Likelihood of a group data set under a model
+|  `pcm_NR`				| Newton Raphson optimisation 
+|  `pcm_minimize`			| Conjugate gradient descent 
+
+### Model Evaluation
 These functions are higher level functions that perform fitting and crossvalidation of either individual data set or group data sets.  
 
 | 	Function 			       | Comment  
 |:-----------------------------|:-----------------------------
 | `pcm_fitModelIndivid`        | Fits G and noise parameter to individual data
-| `pcm_fitModelIndividCrossval`| Within-subject crossvalidation of models 
-| `pcm_fitModelGroup`          | Fit common G to all subjects, using individual noise and scale parameter
+| `pcm_fitModelIndividCrossval`| Within-subject crossvalidation of models
+| `pcm_fitModelGroup`          | Fit common G to all subjects
 | `pcm_fitModelGroupCrossval`  | Between-subject crossvalidation of models  
+| `pcm_setUpFit`         |  Generally prepares models and data for fitting 
+| `pcm_knockModels`     | Inference on components using simple knock-in knock-out likelihoods
+| `pcm_componentPosterior`     | Inference on components by model averaging 
 
+### Visualization functions
+| 	Function 			       | Comment  
+|:-----------------------------|:-----------------------------
+| `pcm_classicalMDS`           | Multidimensional scaling  
+| `pcm_estimateU`         | Estimates voxel-patterns under model M
+| `pcm_estimateW`         | Estimates  voxel-feature weights under  model M
+| `pcm_plotModelLikelihood`    | Displays marginal likelihood for different models
 
-## Utility functions
+### Model building 
 | 	Function 			       | Comment  
 |:-----------------------------|:-----------------------------
-| `pcm_checkderiv`             | Checks derivate of a nonlinear model 
-| `pcm_estGcrossval`           | Cross-validated estimate of G 
-| `pcm_indicatorMatrix`        | Generates indicator matrices 
-| `pcm_vararginoptions`		   | Handles optional input arguments
+| `pcm_constructModelFamily`     | Makes a family of models from different components
+| `pcm_buildModelFromFeatures`     | Makes a component model from featuresets
+
 
-## Visualization functions
+### Utility functions
 | 	Function 			       | Comment  
 |:-----------------------------|:-----------------------------
-| `pcm_classicalMDS`           | Multidimensional scaling  
-| `pcm_plotModelLikelihood`    | Displays marginal likelihood (results from `pcm_fitModel` functions)
-| `pcm_plotFittedG`            | Displays second-moment matrices with fitted parameters
+| `pcm_addModelComp`                  | Adds a model component to model M  
+| `pcm_blockdiag`                  | Makes a blockdiagonal matrix  
+| `pcm_calculateG`                  | Determines G for models   
+| `pcm_checkderiv`                  | Checks derivate of a nonlinear model 
+| `pcm_diagonalize`                  | Decomposes G into A * A'  
+| `pcm_estGcrossval`              | Cross-validated estimate of G 
+| `pcm_generateData`                  | Generates data simulations from model M  
+| `pcm_getStartingval`               | Provides a starting value estimate for model M  
+| `pcm_indicatorMatrix`       | Generates indicator matrices 
+| `pcm_vararginoptions`	    | Deals with variable input options 
+| `pcm_getUserOptions`	    | Deals with variable input options (struct-based, faster) 
+| `pcm_makePD`    | Moves a matrix to closest semi-postive definite Matrix 
+| `pcm_optimalAlgorithm`    | Recommendation for the best optimisation algorithm
+| `pcm_prepFreeModel`         |  Sets up free model (freechol)
+
 
-## Recipes 
+### Recipes 
 | 	Function 			       | Comment  
 |:-----------------------------|:-----------------------------
-| `pcm_recipe_finger`          | Example of fixed ad component models 
-| `pcm_recipe_correlation`     | Example of feature model  
+| `pcm_recipe_finger`          | Example of a fixed and component models 
+| `pcm_recipe_correlation`     | Example of feature model   
 | `pcm_recipe_nonlinear`       | Example for non-linear model   
+
+### Currently not maintained / under developement
+| 	Function 			       | Comment  
+|:-----------------------------|:-----------------------------
+|  `pcm_NR_diag`			| Newton Raphson for diagonalized component models
+|  `pcm_NR_free`			| Newton Raphson for a free model 
+|  `pcm_EM`				| Expectation-Maximization 
diff --git a/pcm_NR.m b/pcm_NR.m
@@ -14,6 +14,7 @@
 %   'HessReg':      Starting regulariser on the Hessian matrix (default starts at 1 
 %                   then being adjusted in decibels)
 %   'regularization': 'L': Levenberg 'LM': Levenberg-Marquardt 
+%   'verbose':     0: No feedback, 1:Important warnings, 2:regularisation 
 % 
 % OUTPUT:
 %   theta : Variance coefficients 
@@ -39,9 +40,12 @@
 OPT=pcm_getUserOptions(varargin,OPT,{'HessReg','thres','numIter','verbose','regularization'});
 
 % Set warning to error, so it can be caught
-warning('error','MATLAB:nearlySingularMatrix');
-warning('error','MATLAB:singularMatrix'); 
-warning('error','MATLAB:illConditionedMatrix'); 
+CATCHEXP = {'MATLAB:nearlySingularMatrix','MATLAB:singularMatrix',...
+    'MATLAB:illConditionedMatrix','MATLAB:posdef',...
+    'MATLAB:nearlySingularMatrix','MATLAB:eig:matrixWithNaNInf'}; 
+for i=1:numel(CATCHEXP)
+    warning('error',CATCHEXP{i});
+end; 
 
 % Initialize Interations
 %--------------------------------------------------------------------------
@@ -70,15 +74,19 @@
                 dtheta   =  dFdhh\dFdh;
                 theta = theta - dtheta;
                 break; 
-            catch % Likely matrix close to singluar
+            catch ME% Likely matrix close to singluar
+                if any(strcmp(ME.identifier,CATCHEXP))
                 OPT.HessReg = OPT.HessReg*10;
                 if (OPT.HessReg>100000) 
                     if (OPT.verbose)
-                        fprintf('Cant regularise second derivative.. Giving up\n');
+                        warning('Cant regularise second derivative.. Giving up\n');
                     end;
                     break;  % Give up  
                 end; 
                 dFdhh = dFdhh_old + diag(diag(dFdhh_old))*OPT.HessReg;
+                else 
+                    ME.rethrow; 
+                end; 
             end; 
         end;
     end;
@@ -89,7 +97,7 @@
     try
         [nl(k),dFdh,dFdhh]=likefcn(theta);
     catch ME  % Catch errors based on invalid parameter settings
-        if any(strcmp(ME.identifier,{'MATLAB:posdef','MATLAB:nearlySingularMatrix','MATLAB:eig:matrixWithNaNInf','MATLAB:singularMatrix'}))
+        if any(strcmp(ME.identifier,CATCHEXP))
             if (k==1) 
                 error('bad initial values for theta'); 
             else 
@@ -103,7 +111,7 @@
     %----------------------------------------------------------------------
     if (k>1 & (nl(k)-nl(k-1))>eps)      % If not....
         OPT.HessReg = OPT.HessReg*10;   % Increase regularisation
-        if (OPT.verbose)
+        if (OPT.verbose==2)
             fprintf('Step back. Regularisation %2.3f\n',OPT.HessReg(1));
         end;
         theta = thetaH(:,k-1);

diff --git a/pcm_componentPosterior.m b/pcm_componentPosterior.m
@@ -1,6 +1,6 @@
 function [postProb,logBayes]=pcm_componentPosterior(likelihood,compI,varargin); 
 % function p=pcm_componentPosterior(likelihood,compI); 
-% Calculates a posterior proability for the presence of each 
+% Calculates a posterior probability for the presence of each 
 % component from a family of models 
 % INPUT: 
 %   likelihood: (numSubj x numModels) estimated marginal log-likelihood of each model. This can be 
@@ -10,7 +10,8 @@
 %               of the components for each fitted (pcm_constructModelFamily)
 % VARARGIN: 
 %  prior:       Prior probability of each model
-%               - A numModelsx1 vector  
+%               - A numModels x 1 vector of prior probabilities for the models 
+%               - A numComponent x 1 vector or prior probabilities for each of the model components                
 %               - A scalar is interpreted as the base-probability for each component to be present  
 %               - If empty, it defaults to a flat prior over models 
 % OUTPUT:      
@@ -23,15 +24,29 @@
 pcm_vararginoptions(varargin,{'prior'}); 
 
 [numModels,numComp]=size(compI); 
+if (size(likelihood,2)~=numModels)
+    error('log-likelihood needs to be a numSubj x numModels matrix'); 
+end; 
 
+% Generate a log-prior for each MODEL
 if (isempty(prior))
     logPrior = log(ones(1,numModels)./numModels);
 elseif isscalar(prior) 
     logPrior = sum(compI.*log(prior)+(1-compI).*log(1-prior),2)'; 
+elseif (numel(prior)==numComp) 
+    logPrior = zeros(1,numModels); 
+    for i=1:numComp 
+        logPrior=logPrior+compI(:,i)'.*log(prior(i))+(1-compI(:,i))'.*log(1-prior(i)); 
+    end; 
+elseif (numel(prior)==numModels) 
+    logPrior = log(prior); 
+else 
+    error('Prior must either have the length numModels, numComponents, 1 or 0'); 
 end; 
 
 logJoint = bsxfun(@plus,likelihood,logPrior);
 logJoint = bsxfun(@minus,logJoint,max(logJoint,[],2))+50;   % Scale the log-joint, so it doesn't give errors with exp
+logLike  = bsxfun(@minus,likelihood,max(likelihood,[],2))+50;
 
 for i=1:numComp 
     indx = compI(:,i)==1; 
@@ -40,10 +55,10 @@
    % calculate logBayes factor in favor of component 
    if sum(indx)==length(indx) || sum(indx)==0
        logPosterior(1:length(logJoint),i)=logPrior(i);
-       logBayes(1:length(logJoint),i)=NaN;
+       logBayes(1:length(logLike),i)=NaN;
    else
        logPosterior(:,i)=log(sum(exp(logJoint(:,indx)),2))-log(sum(exp(logJoint),2));
-       logBayes(:,i)=log(sum(exp(logJoint(:,indx)),2))-log(sum(exp(logJoint(:,~indx)),2));
+       logBayes(:,i)=log(sum(exp(logLike(:,indx)),2))-log(sum(exp(logLike(:,~indx)),2));
    end
 end; 
 

diff --git a/pcm_fitModelIndivid.m b/pcm_fitModelIndivid.m
@@ -56,9 +56,10 @@
 %   'MaxIteration': Number of max minimization iterations. Default is 1000.
 % 
 %   'verbose':      Optional flag to show display message in the command
-%                   line (e.g., elapsed time). Default is 1.
+%                   line. Default is 1. Setting to 2 gives more detailed
+%                   feedback on pcm_NR
 % 
-%   'S':             Optional specific covariance structureof the noise 
+%   'S':             Optional specific covariance structure of the noise 
 % 
 %   'fitAlgorithm': Either 'NR' or 'minimize' - provides over-write for
 %                   model specific algorithms 
@@ -160,7 +161,7 @@
                 T.likelihood(s,m) =  -fX(end);  %invert the sign 
             case 'NR' 
                 fcn = @(x) pcm_likelihoodIndivid(x,YY{s},M{m},Z{s},X{s},P(s),OPT);
-                [theta_hat{m}(:,s),T.likelihood(s,m),T.iterations(s,m),T.reg(s,m)]=pcm_NR(theta0{m}(:,s),fcn,'verbose',verbose==2); 
+                [theta_hat{m}(:,s),T.likelihood(s,m),T.iterations(s,m),T.reg(s,m)]=pcm_NR(theta0{m}(:,s),fcn,'verbose',verbose); 
         end; 
 
         G_pred{m}(:,:,s)  =  pcm_calculateG(M{m},theta_hat{m}(1:M{m}.numGparams,s));