Team:TU Munich/Modeling/Yeast Growth

From 2012.igem.org

(Difference between revisions)
(Code)
(Code)
Line 7: Line 7:
<pre>
<pre>
data1=[1,1.3625,2.875,5.25];
data1=[1,1.3625,2.875,5.25];
 +
t=[0,80,270,450];
t=[0,80,270,450];
-
fun = @(k,data) sum((exp(k*t)-data).^2);
+
fun = @(k,data) sum(((exp(k(1)*t)-data)/k(2)).^2);
-
[k, RN] = fminsearch(@(k) fun(k(1),data1)/k(2) - length(t)*log(1/(sqrt(2*pi*k(2)))),[0.005 1]);
+
[k, RN] = fminsearch(@(k) fun(k,data1) - length(t)*log(1/(sqrt(2*pi*k(2)))),[0.005 1]);
-
% needs to be adapted such that acceptance rate is close to23%
 
kc=0.01;
kc=0.01;
acc=0;
acc=0;
Line 36: Line 36:
     sample(s,:) = kcur;
     sample(s,:) = kcur;
     % calculate the a-posteriori probability of the new sample
     % calculate the a-posteriori probability of the new sample
-
     prob(s) = -fun(kcur,data1)/k(2);
+
     prob(s) = -fun([kcur k(2)],data1);
     if s>1
     if s>1
         % check whether we reject the sample or not. mind that the
         % check whether we reject the sample or not. mind that the
Line 48: Line 48:
             lh(s)=acc/(nacc+acc)*100;
             lh(s)=acc/(nacc+acc)*100;
             s=s+1;
             s=s+1;
-
 
         else
         else
             % count not accepted samples
             % count not accepted samples
Line 58: Line 57:
             kprev = kcur;
             kprev = kcur;
             acc=acc+1;
             acc=acc+1;
 +
            %waitbar((k-1)/n)
             lh(s)=acc/(nacc+acc)*100;
             lh(s)=acc/(nacc+acc)*100;
             s=s+1;
             s=s+1;
Line 64: Line 64:
         end   
         end   
     end
     end
-
 
     catch ME
     catch ME
         disp(ME)
         disp(ME)

Revision as of 14:01, 21 October 2012


Yeast Growth Model

Code

data1=[1,1.3625,2.875,5.25];

t=[0,80,270,450];

fun = @(k,data) sum(((exp(k(1)*t)-data)/k(2)).^2);

[k, RN] = fminsearch(@(k) fun(k,data1) - length(t)*log(1/(sqrt(2*pi*k(2)))),[0.005 1]);

kc=0.01;
acc=0;
nacc=0;
n=10000;
sample=zeros(n,1);
prob=zeros(n,1);
lh=zeros(n,1);
accepted=zeros(n,1);

kprev = k(1);
s=1;
while (s<n+1)
    if(mod(s,500)==0)
        (s)/n;
        display(['Acceptance Rate: ' num2str(acc/(nacc+acc)*100) '%'])
    end
    try
    % the step is sampled from a multivariate normal distribution and
    % scaled with kc
    kcur = kprev + kc*mvnrnd(0,0.001);
    sample(s,:) = kcur;
    % calculate the a-posteriori probability of the new sample
    prob(s) = -fun([kcur k(2)],data1);
    if s>1
        % check whether we reject the sample or not. mind that the
        % probability is log-scaled
        if log(rand(1)) < +prob(s)-prob(s-1)  
            % update current sample
            kprev = kcur;
            % count accepted samples
            acc=acc+1;
            % save current acceptance rate for post-processing
            lh(s)=acc/(nacc+acc)*100;
            s=s+1;
        else
            % count not accepted samples
            nacc=nacc+1;
        end
    else
        % we need to do something else for the first sample.
        if log(rand(1)) < prob(s)+RN
            kprev = kcur;
            acc=acc+1;
            %waitbar((k-1)/n)
            lh(s)=acc/(nacc+acc)*100;
            s=s+1;
        else
            nacc=nacc+1;
        end  
    end
    catch ME
        disp(ME)
    end
end
disp(['Doubling Time: ' num2str(log(2)/k(1)) ' STD: ' num2str(std(log(2)./sample))])