% drifr.m
% Markov chain model for genetic drift
%                       2/22/03     SG


clear all
close all

warning off

tic

% PARAMETERS

T = 100;        % number of generations to run
N =  10;        % number of diploid organisms
p0 = 0.5;       % initial frequency of A



% TRANSITION MATRIX: Prob(i->j)

for i = 0:2*N
    for j = 0:2*N
        pi = i/(2*N);
        qi = 1-pi;
        P(i+1,j+1) = nchoosek(2*N,j)*pi^j*qi^(2*N-j);
    end
end
        


% INITIAL CONDITIONS

n = round(1+2*N*p0);            % initial number of alleles A
                                % extra 1 because no i=0 in P
x(1,:) = zeros(1,2*N+1);       
x(1,n) = 1;                     % initially there are n alleles A

% MAIN PROGRAMM

for t = 2:T
    
    x(t,:) = x(1,:)*P^t;
    
    % GRAPHICS
    
    bar(0:2*N, x(t,:))
    xlabel('number of alleles A')
    ylabel('probability')
    axis([-1 2*N+1 0 1.1*max(x(t,:))])
    title(['Genetic Drift:   t = ', num2str(t)])
    pause(0.1)
    
    
end

toc