06/17/2021 Add fitting parameter confidence intervals (Jingwen)

LarsonLab · Jun 17, 2021 · a784d04 · a784d04
1 parent cfe873d
commit a784d04
Show file tree

Hide file tree

Showing 2 changed files with 32 additions and 4 deletions.
diff --git a/Fitting/utebrain_multiscan_model_fit.m b/Fitting/utebrain_multiscan_model_fit.m
@@ -121,8 +121,11 @@
 %opts = optimset('MaxIter', 2000, 'MaxFunEvals', 2000, 'TolFun', 10^(-7), 'TolX', 10^(-7));
 lsq_opts = optimset('Display','none','MaxIter', 500, 'MaxFunEvals', 500);
 
-[X,resnorm,residual,exitflag] = lsqnonlin(@model_diff, X0, lb, ub, lsq_opts);
-%exitflag
+[X,resnorm,residual,~,~,~,J] = lsqnonlin(@model_diff, X0, lb, ub, lsq_opts);
+
+CI = nlparci(X,residual,'jacobian',J);
+X_lb = CI(:,1);
+X_ub = CI(:,2);
 
 fit_result = struct('rho',{}, 'T2',{}, 'df', {}, 'phi',{});
 % T2s = X(4*[0:general_opts.num_components-1] + 2);
@@ -134,6 +137,16 @@
     fit_result(n).T2 = X(component_offset + num_scans+1);
     fit_result(n).df = X(component_offset + num_scans+2);
     fit_result(n).phi = X(component_offset + num_scans+2+ [1:num_scans]);
+    % lower bound
+    fit_result(n).rho_lb(1:num_scans) = X_lb(component_offset + [1:num_scans])*Snorm;
+    fit_result(n).T2_lb = X_lb(component_offset + num_scans+1);
+    fit_result(n).df_lb = X_lb(component_offset + num_scans+2);
+    fit_result(n).phi_lb = X_lb(component_offset + num_scans+2+ [1:num_scans]);
+    % upper bound
+    fit_result(n).rho_ub(1:num_scans) = X_ub(component_offset + [1:num_scans])*Snorm;
+    fit_result(n).T2_ub = X_ub(component_offset + num_scans+1);
+    fit_result(n).df_ub = X_ub(component_offset + num_scans+2);
+    fit_result(n).phi_ub = X_ub(component_offset + num_scans+2+ [1:num_scans]);
 end
 
 rmse =sqrt(resnorm);

diff --git a/Fitting/utebrain_t1_model_fit.m b/Fitting/utebrain_t1_model_fit.m
@@ -133,8 +133,11 @@
 lsq_opts = optimset('Display','none','MaxIter', 500, 'MaxFunEvals', 500);
 
 
-[X,resnorm,residual,exitflag] = lsqnonlin(@model_diff, X0, lb, ub, lsq_opts);
-%exitflag
+[X,resnorm,residual,~,~,~,J] = lsqnonlin(@model_diff, X0, lb, ub, lsq_opts);
+
+CI = nlparci(X,residual,'jacobian',J);
+X_lb = CI(:,1);
+X_ub = CI(:,2);
 
 fit_result = struct('rho',{}, 'T2',{}, 'df', {}, 'phi',{}, 'T1', {});
 % T2s = X(4*[0:general_opts.num_components-1] + 2);
@@ -146,6 +149,18 @@
     fit_result(n).df = X(component_offset + 3);
     fit_result(n).phi = X(component_offset + 3+ [1:num_scans]);
     fit_result(n).T1 = X(component_offset + 1*num_scans+4);
+    % lower bound
+    fit_result(n).rho_lb = X_lb(component_offset + 1)*Snorm;
+    fit_result(n).T2_lb = X_lb(component_offset + 2);
+    fit_result(n).df_lb = X_lb(component_offset + 3);
+    fit_result(n).phi_lb = X_lb(component_offset + 3+ [1:num_scans]);
+    fit_result(n).T1_lb = X_lb(component_offset + 1*num_scans+4);
+    % upper bound
+    fit_result(n).rho_ub = X_ub(component_offset + 1)*Snorm;
+    fit_result(n).T2_ub = X_ub(component_offset + 2);
+    fit_result(n).df_ub = X_ub(component_offset + 3);
+    fit_result(n).phi_ub = X_ub(component_offset + 3+ [1:num_scans]);
+    fit_result(n).T1_ub = X_ub(component_offset + 1*num_scans+4);
 end
 
 rmse =sqrt(resnorm);