Extracting data from a plot within a for loop

Question

In MATLAB, I have coded the stochastic simulation algorithm (Gillespie) for a simple birth-death process, and have gotten a plot by using hold on in a for loop.

I have 100 values of PStoch for each Qp value, because I have ran 100 simulations for each Qp value. The values are hard to see in the plot below because they are all clustered together.

How can I save the data from the plot in a matrix, so I can do some calculations on them afterward? Specifically, I need a matrix of size 100 x 100 with all PStoch values corresponding to each Qp value.

My code is below:

rng('shuffle')

%% Pre-defined variables
Qpvec = logspace(-2,1,100);
len = length(Qpvec);
delta = 1e-3;
P0vec = Qpvec./delta;
V = [1,-1];
tmax = 10000;

%% Begin simulation
figure(1)
for k = 1:len
    t0 = 0;
    tspan = t0;
    Qp = Qpvec(k);
    P0 = P0vec(k);
    Pstoch = P0;
    while t0 < tmax && length(Pstoch) < 100
        a = [Qp, delta*P0];
        tau = -log(rand)/sum(a);
        t0 = t0 + tau;
        asum = cumsum(a)/sum(a);
        chosen_reaction = find(rand < asum,1);
        if chosen_reaction == 1;
            P0 = P0 + V(:,1);
        else
            P0 = P0 + V(:,2);
        end
        tspan = [tspan,t0];
        Pstoch = [Pstoch;P0];
    end
    plot(Qp,Pstoch)
    hold on
    axis([0 max(Qp) 0 max(Pstoch)])
end

The plot is here:

Thanks for help.

Make `Qp` a vector and `Pstoch` a matrix, and index everything into those, rather than scalars. — David, Dec 15 '15 at 03:22
Can you elaborate? I understand making vectors and matrices, but not 'indexing everything'. — abscissa, Dec 15 '15 at 03:31
Actually, you don't even need `Qp`, just use `Qpvec` in the plot. `Pstoch` is at more 100 columns, and `len` rows, so make it `zeros(len,100)`, then put elements in it like `Pstoch(k,i)`, you just have to work out `i`, just count the number of iterations of your while loop. — David, Dec 15 '15 at 03:58

score 2 · Accepted Answer · answered Dec 15 '15 at 06:54

I've added three lines in the code below. This assumes you are right in saying that Pstoch always has 100 elements (or less than 100):

Pstoch_M = zeros(len, 100)

for

    k = 1:len
    t0 = 0;
    tspan = t0;
    Qp = Qpvec(k);
    P0 = P0vec(k);

    Pstoch = zeros(100,1);
    counter = 1;    

    Pstoch(1) = P0;
    while t0 < tmax && counter < 100 %// length(Pstoch) < 100
        a = [Qp, delta*P0];
        tau = -log(rand)/sum(a);
        t0 = t0 + tau;
        asum = cumsum(a)/sum(a);
        chosen_reaction = find(rand < asum,1);
        if chosen_reaction == 1;
            P0 = P0 + V(:,1);
        else
            P0 = P0 + V(:,2);
        end
        counter = counter + 1;
        tspan = [tspan,t0];
        Pstoch(counter) P0;;
    end

    Pstock_M(:,k) = Pstoch;

    plot(Qp,Pstoch)
    hold on
    axis([0 max(Qp) 0 max(Pstoch)])
end

Note the pre-allocation added for Pstoch should make your code considerably faster. You should do the same for tspan. Growing variables inside loops in MATLAB is extremely inefficient, as m-lint is no doubt currently warning you.

Extracting data from a plot within a for loop

1 Answers1