matlab: praticle state simulation

Question

Lets say I want to simulate a particle state, which can be normal (0) or excited (1) in given frame. The particle is in excited state f % of time. If the particle is in excited state, it lasts for ~L frames (with poisson distribution). I want to simulate that state for N time points. So the input is for example:

N = 1000;
f = 0.3;
L = 5;

and the result will be something like

state(1:N) = [0 0 1 1 1 1 0 0 0 0 1 1 1 0 0 0 0 0 0 0 0 0 1 1 1 1 1 0 0 0 0 0 ... and so on]

with sum(state)/N close to 0.3

How to do that? Thanks!

I don't really know what you mean by that. What I exactly want to do, is to simulate diffusion behavior of a particle with two different diffusion coefficients and defined fraction of faster and slower component (f) and some kind of lifetime in one state or another. I wanted to simulate state first (in this case two, but possibly more) and then just simulate displacement and coordinates, depending on state (faster or slower...). I don't know if it's the best way to do that, but that was the first one in my mind :) — Art, Mar 25 '12 at 08:53
@NoamN.Kremen As f=0.3, and the length of state 1 is 5. The length of state zero on average should be about 17 (5/0.3), so the change of flipping from 0 to 1 is 0.06. Edit: Not sure if this statement is completely true. — Bernhard, Mar 25 '12 at 10:00

score 2 · Answer 1 · answered Mar 25 '12 at 10:07

The average length of the excited state is 5. The average length of the normal state, should thus be around 12 to obtain.

The strategy can be something like this.

Start in state 0
Draw a random number a from a Poisson distribution with mean L*(1-f)/f
Fill the state array with a zeroes
Draw a random number b from a Poission distribution with mean L
Fill the state array witb b ones.
Repeat

Another option would be to think in terms of switching probabilities, where the 0->1 and 1->0 probabilities are unequal.

Serg · Accepted Answer · 2012-03-25T22:45:35.363

%% parameters
f = 0.3; % probability of state 1
L1 = 5;  % average time in state 1
N = 1e4;
s0 = 1; % init. state
%% run simulation
L0 = L1 * (1 / f - 1); % average time state 0 lasts
p01 = 1 / L0; % probability to switch from 0 to 1
p10 = 1 / L1; % probability to switch from 1 to 0
p00 = 1 - p01;
p11 = 1 - p10;
sm = [p00, p01; p10, p11];  % build stochastic matrix (state machine)
bins = [0, 1]; % possible states
states = zeros(N, 1);
assert(all(sum(sm, 2) == 1), 'not a stochastic matrix');
smc = cumsum(sm, 2); % cummulative matrix
xi = find(bins == s0);
for k = 1 : N
    yi = find(smc(xi, :) > rand, 1, 'first');
    states(k) = bins(yi);
    xi = yi;
end
%% check result
ds = [states(1); diff(states)];
idx_begin = find(ds == 1 & states == 1);
idx_end = find(ds == -1 & states == 0);
if idx_end(end) < idx_begin(end)
    idx_end = [idx_end; N + 1];
end
df = idx_end - idx_begin;
fprintf('prob(state = 1) = %g; avg. time(state = 1) = %g\n', sum(states) / N, mean(df));

matlab: praticle state simulation

2 Answers2