From a784d04546d9585118a94877b5dee9c8ab714c55 Mon Sep 17 00:00:00 2001 From: Jingwen Yao Date: Thu, 17 Jun 2021 16:42:44 -0700 Subject: [PATCH] 06/17/2021 Add fitting parameter confidence intervals (Jingwen) --- Fitting/utebrain_multiscan_model_fit.m | 17 +++++++++++++++-- Fitting/utebrain_t1_model_fit.m | 19 +++++++++++++++++-- 2 files changed, 32 insertions(+), 4 deletions(-) diff --git a/Fitting/utebrain_multiscan_model_fit.m b/Fitting/utebrain_multiscan_model_fit.m index 2c84ad4..3566da6 100644 --- 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 index 87119b7..589d7ed 100644 --- 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);