Commit 01143eec authored by Antonio Politi's avatar Antonio Politi

Include computation of CI

parent a42475b1
......@@ -52,10 +52,13 @@ classdef absDelayedMaturation < handle
% Run fit of model
% ifit: index of parameter(s) to fit
% par: parameter values
function [parf, norm] = fitModel(MO, ifit, par)
function [parf, norm] = fitModel(MO, ifit, par, options)
if nargin < 4
options = optimset('Display', 'iter');
end
lb = MO.LB(ifit);
hb = MO.HB(ifit);
options = optimset('Display', 'iter');
[parf, norm] = lsqnonlin(@MO.distMO, par(ifit), lb, hb, ...
options, ifit, par);
end
......@@ -143,7 +146,7 @@ classdef absDelayedMaturation < handle
parfit = par;
parfit(ifit) = parfit(ifit)+ fact*parfit(ifit).*(1-2*rand(1,length(ifit)));
parfit = (parfit > MO.HB).*(MO.HB) + (parfit <= MO.LB).*(MO.LB) + (parfit <= MO.HB).*(parfit > MO.LB).*parfit;
[parf, norm] = MO.fitModel(ifit, parfit);
[parf, norm] = MO.fitModel(ifit, parfit, optimset('Display', 'iter'));
parfit(ifit) = parf;
for i=1:length(parfit)-1
fprintf(fid, '%.3e\t', parfit(i));
......@@ -157,6 +160,91 @@ classdef absDelayedMaturation < handle
fclose all;
end
%%
% par: parameter values
% ifit: indexes of parameters to fit
% iC: index of paramater to systematically vary
% bound: Boundary where to search
% step: number of steps
% outfile: output file
%%
function findCI(MO, par, ifit, iC, step, outfile, bounds)
NrDataPts = size(MO.dataF(1).time,1)*2;
if nargin < 7
%% find boundary of parameter values
options = optimset('Display', 'iter', 'MaxFunEvals', 100)
parloc = par;
[parf, norm] = MO.fitModel( ifit, parloc);
parL = par(iC)/20;
parH = par(iC);
parloc(iC) = parL;
[parf, normL] = MO.fitModel(ifit, parloc, options);
logproL = MO.LogPro(normL, norm, NrDataPts);
if (abs(logproL) > chi2inv(0.95, 1))
while abs(parL - parH) > par(iC)/5
parL2 = (parH + parL)/2;
parloc(iC) = parL2;
[parf, normL] = MO.fitModel(ifit, parloc, options);
logproL = MO.LogPro(normL, norm, NrDataPts)
if abs(logproL) < chi2inv(0.95, 1)
parH = parL2;
else
parL = parL2;
end
end
else
lowPar = parL;
end
lowPar = parL;
parL = par(iC);
parH = 20*par(iC);
parloc(iC) = parH;
[parf, normL] = MO.fitModel(ifit, parloc, options);
logproL = MO.LogPro(normL, norm, NrDataPts);
if abs(logproL) > chi2inv(0.95, 1)
while abs(parL - parH) > par(iC)/5
parH2 = (parH + parL)/2;
parloc(iC) = parH2;
[parf, normL] = MO.fitModel(ifit, parloc, options);
logproL = MO.LogPro(normL, norm, NrDataPts)
if abs(logproL) < chi2inv(0.95, 1)
parL = parH2;
else
parH = parH2;
end
end
else
highPar = parH;
end
highPar = parH;
else
lowPar = bounds(1);
highPar = bounds(2);
end
lowPow = floor(log10(lowPar/par(iC))*10)/10;
highPow = ceil(log10(highPar/par(iC))*10)/10;
parRange = par(iC)*10.^[lowPow*[step:-1:1]/step highPow*[0:step]/step];
%try to find lower border
for i = 1:length(parRange)
parloc = par;
parloc(iC) = parRange(i);
MO.runFit(parloc, ifit, i, 0, outfile);
end
pCI = load(fullfile(MO.outdir, outfile));
%%
CI = pCI(:,iC).*(MO.LogPro(pCI(:,end-1), min(pCI(:,end-1)), NrDataPts)<chi2inv(0.95,1));
CITM = (1./pCI(:,3) + pCI(:,4)).*(MO.LogPro(pCI(:,end-1), min(pCI(:,end-1)), NrDataPts)<chi2inv(0.95,1))
CI = [min(CI(CI>0)) max(CI(CI>0))]
CITM = [min(CITM(CITM>0)) max(CITM(CITM>0))]
end
function logPro = LogPro(MO, norm1,norm2, nrpts)
logPro = nrpts*log(norm1/nrpts) - nrpts*log(norm2/nrpts);
end
function par = getBestfit(MO, filename)
try
parV = load(filename);
......
......@@ -23,7 +23,7 @@ classdef delayedMaturation < absDelayedMaturation
end
MO.parnames = 'ts v kM tauM kd P01 M01 P02 M02';
MO.LB = [0, 0, 0.01, 10, 0, 0, 0,0,0]; %Low boundary values
MO.LB = [0, 0, 0.001, 0, 0, 0, 0,0,0]; %Low boundary values
MO.HB = [15, 1, 10, 70, 1, 10000, 10000, 10000, 10000]; %High boundary values
MO.name = mfilename('class')
MO.N = 3;
......
%%
%fit just Mature pore
clear;
clear all;
fclose all
MO = delayedMaturationExponentialV('24h', 0,0)
MO.useIntermediate = 0;
tsv = [5 7 10 12 15];
nrFitvec = [5 5 32 5 5];
MO.outdir = fullfile(MO.outdir, 'delay')
%%
for its = 3 %1:length(tsv)
%%
ts = tsv(its);
nrfits = nrFitvec(its);
At0 = MO.getSurf(ts);
P01 = 0*4.493*At0;
M01 = 4.79*At0;
P02 = 0*3.935*At0;
M02 = 9.048*At0;
%%
MO.useIntermediate = 1;
ts = tsv(its);
At0 = MO.getSurf(ts);
P01 = 0*4.493*At0;
M01 = 4.79*At0;
P02 = 0*3.935*At0;
M02 = 9.048*At0;
par = [ts, 8.993, 2, 0.015, 2, 20, 4.2e-04, P01, M01, P02, M02];
ifit = [2:3 5:6 9 11];
parfile = ['delayedMaturationExponentialV_noPrecursor_ts' num2str(ts) '_v0.015'];
[dist, distI, distM] = afitround(MO, par, ifit, 35, parfile)
end
%%
iC = 6;
ifit = [2:3 5 9 11];
parf = MO.getBestfit(fullfile(MO.outdir, [parfile '.txt']));
MO.findCI(parf,ifit, iC, 25, [parfile '_tauMIC.txt'], [7.066e-01 5.613e+01]);
%%
iC = 5;
ifit = [2:3 6 9 11];
parf = MO.getBestfit(fullfile(MO.outdir, [parfile '.txt']))
MO.findCI(parf,ifit, iC, 25, [parfile '_kMIC.txt'], [1.579e-02 2.503e+00])
%%
%fit just Mature pore
clear;
clear all;
fclose all
MO = delayedMaturation('24h', 0,0)
MO.outdir = fullfile(MO.outdir, 'delay')
MO.useIntermediate = 0;
tsv = [5 7 10 12 15];
nrFitvec = [10 10 10 10 10];
%%
for its = 3 %1:length(tsv)
ts = tsv(its);
MO.useIntermediate = 0;
nrfits = nrFitvec(its);
At0 = MO.getSurf(ts);
P01 = 0*4.493*At0;
M01 = 4.79*At0;
P02 = 0*3.935*At0;
M02 = 9.048*At0;
%%
P01 = 4.493*At0;
P02 = 3.935*At0;
MO.useIntermediate = 1;
par = [ts, 0.015, 1, 40, 4.2e-4, P01, M01, P02, M02];
ifit = [3:4 6:9];
parfile = ['delayedMaturation_v0.015_ts' num2str(ts)];
[dist, distI, distM] = afitround(MO, par, ifit, nrfits, parfile)
end
%%
iC = 4;
ifit = [3 6:9];
parf = MO.getBestfit(fullfile(MO.outdir, [parfile '.txt']));
MO.findCI(parf,ifit, iC, 25, [parfile '_tauMIC.txt'], [7.066e-01 5.613e+01]);
%%
iC = 3;
ifit = [4 6:9];
parf = MO.getBestfit(fullfile(MO.outdir, [parfile '.txt']))
MO.findCI(parf,ifit, iC, 25, [parfile '_kMIC.txt'], [1.579e-02 2.503e+00])
%%
MO.useIntermediate = 1;
ts = tsv(its);
nrfits = nrFitvec(its);
At0 = MO.getSurf(ts);
P01 = 0*4.493*At0;
M01 = 4.79*At0;
P02 = 0*3.935*At0;
M02 = 9.048*At0;
par = [ts, 0.014, 1.136, 40, 4.2e-04, P01, M01, P02, M02];
ifit = [2:4 7 9];
parfile = ['delayedMaturation_noPrecursor_ts' num2str(ts)];
afitround(MO, par, ifit, nrfits, parfile);
%%
iC = 4;
ifit = [3 7 9];
parf = MO.getBestfit(fullfile(MO.outdir, [parfile '.txt']));
MO.findCI(parf,ifit, iC, 25, [parfile '_tauMIC.txt'], [7.257e-01 3.637e+01]);
%%
iC = 3;
ifit = [4 7 9];
parf = MO.getBestfit(fullfile(MO.outdir, [parfile '.txt']))
MO.findCI(parf,ifit, iC, 25, [parfile '_kMIC.txt'], [0.03 2] ) %no changes up to a value of 20
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment