0

在 MATLAB 中，我为一个简单的生死过程编写了随机模拟算法 (Gillespie)，并通过循环使用得到hold on了一个绘图。for

PStoch每个值我有 100 个Qp值，因为我为每个Qp值运行了 100 次模拟。在下图中很难看到这些值，因为它们都聚集在一起。

如何将绘图中的数据保存在矩阵中，以便之后对它们进行一些计算？具体来说，我需要一个大小为 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

剧情在这里：在此处输入图像描述

感谢帮助。

4

1 回答 1

2

我在下面的代码中添加了三行。这假设您说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 目前无疑会警告您。

于 2015-12-15T06:54:40.720 回答