0

MATLAB で、確率シミュレーション アルゴリズム (Gillespie) を単純な生死プロセス用にコーディングhold onし、forループ内で使用してプロットを取得しました。

値ごとに 100 のシミュレーションを実行したため、値PStochごとに 100 の値があります。値はすべて一緒にクラスター化されているため、以下のプロットでは見にくいです。QpQp

後で計算できるように、プロットからのデータをマトリックスに保存するにはどうすればよいですか? PStoch具体的には、すべての値が各値に対応するサイズ 100 x 100 のマトリックスが必要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

以下のコードに 3 行追加しました。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 に答える