simulation help - Programmers Heaven

Howdy, Stranger!

It looks like you're new here. If you want to get involved, click one of these buttons!

simulation help

Posts: 2Member
Hello,

I am trying to generate a stochastic run of 100 simulations and plot the result of each simulation in a histogram. Thus far, I have everything correct except that I can only get 1 simulation each time in my histogram.

Here is my code. Any suggestions for getting all 100 values in it?

n=100; % number of generations to run
x=zeros(1,n+1); %Vector to keep track of time/given generation for potential fixation
t=0:n; % time indices
allele=.5; % Initial frequency of focal alleles
Nlow=2; %Low population size
Nhigh=10; %High population size
mu=[2,10];%Initial distribution
histo2e=zeros(1,n+1);
%2 state Markov Chain; Nlow_pop_size to Nhigh_pop_size
P=[[.2,.8];[.4,.6]];

x(1)=rando(mu); % generate first x value (time 0, not time 1)

for i=1:n %Cycle through 100 generations

%Generate value based on transition probability Matrix

x(i+1) = rando(P(x(i),);

gencount=zeros(1,n);
N=x(i+1);%Variable to keep track of varying population size for each simulation

%Initialize frequency of focal allele to 0.5
init=1;
thestate=init;%Initial state begins at 0.5
m=0; %Keep track number of times given state is reached
maxgen=100; % Random cutoff for W-F simulation
traj=zeros(1,maxgen);%A useful vector

while (m < maxgen && thestate > 0 && thestate < 2*N)
m=m+1; %Bump counter
%Sample from binomial scheme W-F here
thestate=binornd(2*N,thestate/(2*N));
traj(m) = thestate;

if ((thestate==0 || thestate==2*N) && m < maxgen) % fixation has occured: fill the rest of traj with the appropriate value
traj((m+1):maxgen) = repmat(thestate,1,maxgen-m);

end; %if

end;%while

twin=max(m);

gencount(i)=twin;

continue;

end

hist(gencount(i)); My 100 generated values here!