Извлечение данных из графика в цикле for

В MATLAB я закодировал алгоритм стохастического моделирования (Гиллеспи) для простого процесса рождения-смерти и получил график, используя hold on в цикле for.

У меня есть 100 значений PStoch для каждого значения Qp, потому что я провел 100 симуляций для каждого значения Qp. Значения трудно увидеть на графике ниже, потому что все они сгруппированы вместе.

Как сохранить данные с графика в матрицу, чтобы впоследствии выполнить над ними некоторые вычисления? В частности, мне нужна матрица размером 100 x 100 со всеми значениями PStoch, соответствующими каждому значению Qp.

Мой код ниже:

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

Сюжет здесь: введите здесь описание изображения

Спасибо за помощь.


person abscissa    schedule 15.12.2015    source источник
comment
Сделайте Qp вектором и Pstoch матрицей и индексируйте все в них, а не в скаляры.   -  person David    schedule 15.12.2015
comment
Можете ли вы уточнить? Я понимаю создание векторов и матриц, но не «индексировать все».   -  person abscissa    schedule 15.12.2015
comment
На самом деле вам даже не нужно Qp, просто используйте Qpvec в сюжете. Pstoch содержит более 100 столбцов и len строк, поэтому сделайте его zeros(len,100), затем поместите в него такие элементы, как Pstoch(k,i), вам просто нужно отработать i, просто подсчитайте количество итераций вашего цикла while.   -  person David    schedule 15.12.2015


Ответы (1)


Я добавил три строки в приведенный ниже код. Это предполагает, что вы правы, говоря, что Pstoch всегда имеет 100 элементов (или меньше 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

Обратите внимание, что предварительное выделение, добавленное для Pstoch, должно значительно ускорить ваш код. Вы должны сделать то же самое для tspan. Выращивание переменных внутри циклов в MATLAB крайне неэффективно, о чем m-lint, несомненно, в настоящее время предупреждает вас.

person Dan    schedule 15.12.2015