"

From 2012.igem.org

---

Yeast Growth Model

Code

data1=[1,1.3625,2.875,5.25];
data1=[1,2.6279,8.3721,12.6279];
data1=[1,2.0299,6.8955,9.5671];
data1=[1,1.12777,1.68888,4.044444,6.5555];
data1=[1,0.9339,1.2264,2.0377,4.7169];
data1=[1,1.5,2.2857,6.3571,14.2857];
data1=[1,1.1666,1.51851,3.85185,8.87037];
data1=[1,1.1961,1.7058,5.0980,12.4509];
data1=[1,1.0888888,0.95555,1.0222222];
data1=[1,0.87234,0.936170,1.361702,3.510638];
data1=[1,1.48863,2.272727,9.363636];
data1=[1,1.441860,2.465116,7.2558139];
data1=[1,1.52439,2.17073,9.097560];
data1=[1,1.04109,1.232876,1.6712328];
data1=[1,0.70967,0.854837,1.0322580,2.9516162];
data1=[1,0.93617,1.021276,1.7446808,3.8297872];
data1=[1,2.342857,4.342857,14.8857142];
data1=[1,1.8625,3.425,12.4];
data1=[1,2.4848,4.3030303,15.45454545];
data1=[1,1.603773,4.603773];
data1=[1,1.895844,6.541666];
data1=[1,2.702127,8.936170];
data1=[1,1.42857,4.5];
data1=[1,1.4791,3.83333];
data1=[1,3.15,11.25];
data1=[1,2.56,9.12];
data1=[1,2.4814,10.1851];
data1=[1,2.393939,9.363636];
data1=[1,2,8,17.25];
data1=[1,2.16666,9,16.1666];
data1=[1,2.1428,9.5714,24];


t=[0,80,270,450];
t=[0,70,145,340,520];
t=[0,70,145,340];
t=[0,80,210];
t=[0,100,250,335];
t=[1,2.5,14.16666,34.1666]



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

[k, RN] = fminsearch(@(k) fun(k(1),data1)/k(2) - 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,data1)/k(2);
    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;
            % update waitbar
            %waitbar((k-1)/n)
            % save current acceptance rate for post-processing
            lh(s)=acc/(nacc+acc)*100;
            s=s+1;
            %p = [kcur sigma];
            %plotETOHPROM
            %drawnow
        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))])