forked from mgschiavon/GZM_InducibleTF
-
Notifications
You must be signed in to change notification settings - Fork 0
/
FN_FitMRW.m
90 lines (84 loc) · 2.97 KB
/
FN_FitMRW.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
%% [GZM] Inducible Transcription Factors %%
% ------------------------------------------- %
% FUNCTION: Perform MCMC parameter estimation %
% Created by Mariana Gómez-Schiavon
% August 2019
% FN_FitMRW : Find the set of parameters that best fit the data.
%
% [] = FN_FitMRW(X,H,p,M,D,s,fp,ExID,printAll)
% X : Vector (array) of TF concentration
% H : Vector of inducer (hormone) concentration
% M : Transcriptional model to consider
% ('Mechanistic','HillxBasal','SimpleHill')
% p : Structure with the kinetic parameters
% D : Measured output (data) matrix [length(H) x length(X)]
% s : Random number generator seed
% f : Structure (array) with information to fit parameters
% .par : Parameter to fit (e.g. 'mY')
% .cov : Covariance to calculate parameter perturbations (e.g. 1e-3)
% .lim : Range of acceptable values (e.g. [0,1])
% I : Number of iterations
% ExID : Code for output file name
% printAll : Flag for printing full random walk
%
% OUTPUT bestP : Array of the best set of parameters
% minE : Error of the best set of parameters
%
% See also FN_SS_SimpleHill.m
% See also FN_SS_HillxBasal.m
% See also FN_SS_Mechanistic.m
% See also FN_SS_Allosteric.m
% See also FN_FitError.m
function [bestP,minE] = FN_FitMRW(X,H,p,M,D,s,f,I,ExID,printAll)
mrw.s = s;
mrw.f = f;
mrw.P = zeros(I,length(f)); % OUTPUT: Parameters.
mrw.e = zeros(I,1); % OUTPUT: Error function values.
Mstep = 2;
% (1) Define random number generator:
rng(s,'twister');
r.tL = rand(I,1); % To evaluate proposal acceptance.
r.Pe = zeros(I,length(f)); % Parameter perturbations.
for i = 1:length(f)
r.Pe(:,i) = mvnrnd(zeros(I,1),f(i).cov);
% (2) Initialize system:
mrw.P(1,i) = 10.^((rand()*(log10(f(i).lim(2)) ...
- log10(f(i).lim(1)))) ...
+ log10(f(i).lim(1)));
p.(f(i).par) = mrw.P(1,i);
end
mrw.e(1) = FN_FitError(X*p.nM,H,p,M,D*p.nM);
% (3) Iterate:
for j = 2:I
mrw.P(j,:) = mrw.P(j-1,:);
mrw.e(j) = mrw.e(j-1,:);
% Alternative parameter set:
for i = 1:length(f)
p.(f(i).par) = min(max(f(i).lim(1),mrw.P(j,i)*(Mstep^r.Pe(j,i))),f(i).lim(2));
end
% Generate proposal:
myE = FN_FitError(X*p.nM,H,p,M,D*p.nM);
% If proposal is accepted, update system:
if(r.tL(j) < exp(mrw.e(j)-myE))
for i = 1:length(f)
mrw.P(j,i) = p.(f(i).par);
end
mrw.e(j) = myE;
end
% Save progress:
if(printAll && (mod(j,10000)==0))
j0 = j + 1
save('TEMP_MRW.mat','mrw','j0','r','p');
end
end
clear j i myE
% (4) Save:
if(printAll)
save(cat(2,'MRW_',ExID,'_s',num2str(s),'.mat'),'mrw','p');
delete('TEMP_MRW.mat');
end
[a b] = min(mrw.e);
bestP = mrw.P(b,:);
minE = a;
clear a b
end