matlab: симуляция состояния статьи - PullRequest
Новая
       89

matlab: симуляция состояния статьи

2 голосов
/ 25 марта 2012

Допустим, я хочу смоделировать состояние частицы, которое может быть нормальным (0) или возбужденным (1) в данном кадре. Частица находится в возбужденном состоянии f% времени. Если частица находится в возбужденном состоянии, она длится ~ L кадров (с распределением Пуассона). Я хочу смоделировать это состояние для N временных точек. Таким образом, ввод, например: 

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

и результат будет что-то вроде 

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]

с суммой (состоянием) / N, близкой к 0,3

Как это сделать? Спасибо!

Ответы [ 2 ]

2 голосов
/ 25 марта 2012 
%% 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));
2 голосов
/ 25 марта 2012

Средняя длина возбужденного состояния равна 5. Средняя длина нормального состояния, таким образом, должна составлять около 12 для получения.

Стратегия может быть примерно такой.

  • Запуск в состоянии 0
  • Нарисуйте случайное число a из распределения Пуассона со средним L*(1-f)/f
  • Заполнить массив состояний a нулями
  • Нарисуйте случайное число b из распределения Poission со средним L
  • Заполнить массив состояний b единицами.
  • Повторите

Другой вариант - думать о терминах вероятностей переключения, где вероятности 0-> 1 и 1-> 0 неравны.

...