% fixation.m
% computes fixation probability using from matrix P
% assuming that the two absorptions states are the first and last one
% fixations menas being absorpted at the last one
%  formula:  pi=(I-\tilde{P})^(-1)*P_{i,omega}
%
%                   2/22/03     SG


Q = P;

Q(1,:) = [];        % getting rid of the first row and column corresponding to the other state
Q(:,1) = [];

omega = Q(:,end);   % omega list prob. getting to Omega in one step
omega(end) = [];

Q(:,end) = [];      % getting rid of the last column and row which went to omega
Q(end,:) = [];

prob_of_fixation = inv( eye(length(omega))-Q)*omega;

plot(prob_of_fixation)
xlabel('initial numer of alleles A')
ylabel('probability A is fixed')

J = ones(length(omega),1);
T = inv( eye(length(omega))-Q)*J;

figure
plot(T)
xlabel('initial numer of alleles A')
ylabel('average time A is fixed or lost')