simulation help - Programmers Heaven

Howdy, Stranger!

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

Categories

Welcome to the new platform of Programmer's Heaven! We apologize for the inconvenience caused, if you visited us from a broken link of the previous version. The main reason to move to a new platform is to provide more effective and collaborative experience to you all. Please feel free to experience the new platform and use its exciting features. Contact us for any issue that you need to get clarified. We are more than happy to help you.

simulation help

fraggle_28fraggle_28 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!

Sign In or Register to comment.