function X=mdsim(P,n,start,option)

% mdsim(P,n,start) simulerar markovkedja upp till tid 1,
% P övergångsmatris, n antal tidpunkter, start=starttillstånd.
% Om starttillstånd ej anges ansätts det till 1.
% Tillståndsrum {1,2,...,N} om P NxN-matris.
% Första figuren visar en simulering.
% Den andra figuren visar ett stolpdiagram med relativa tider i tillstånden.
% Om en entydig stationär fördelning existerar, läggs denna in i
% stolpdiagrammet. Anropet 
% 
% X=mdsim(P,n,start) ger dessutom en vektor med de simulerade
% tillstånden. Anropet 
%
% mdsim(P,n,start,pauslgd) där pauslgd är icke-negativt tal ger
% simulering med pausen i figur 1 satt till pauslgd (tiondels sekunder,
% standardvärde 1). Anropet 
%
% X=mdsim(P,n,start,pauslgd) ger dessutom vektor med de simulerade
% tillstånden. Anropet 
%
% X=mdsim(P,n,start,fig) där fig ett negativt tal, ger vektor med
% de simulerade tillstånden, men ej några figurer.

% Jan Enger, KTH,  ht 2001

[N, N1]=size(P); 

if nargin==4 & option>=0
  pauslgd=option;
else
  paudlgd=1;    
end

% koll om P övergångsmatris
if any(P(:)<-10^(-12)) | max(abs(sum(P')-ones(size(P(1,:)))))>10^(-12)|N~=N1
  disp('       P ej övergångsmatris');
  return;
end

if n~=round(n) | n<0
  disp('   felaktigt antal tidpunkter');
  return
end

if nargin>=3
  if start<0 | start >N | start~=round(start)
    disp('   felaktigt startvärde');
    return
  end
   x(1)=start;
else
   x(1)=1;
end

for k=2:n
  % söker nytt tillstånd att hoppa till
  p=P(x(k-1),:); % sannfkn för hoppet
  F=cumsum(p); % fördelningsfkn för hoppet
  u=rand(1); 
  x(k)=min(find(u<F)); % nytt slumpmässigt tillstånd
end


if nargout==1
  X=x;
end

if nargin==4 & option<0
  return
end

figure
pos=get(gcf,'posit');
set(gcf,'posit',pos+[0 0 200 0]); % drar ut plotten
plot([1 2],[x(1) x(1)])
axis([0 n+1 0 N+1])
c=get(gca,'ch');
set(gca,'ytick',1:N)

title('Simulering markovkedja');
xlabel('tid')
ylabel('tillstånd')

hold on

for j=2:n
   yg=get(c,'Yda');
   xg=get(c,'Xda');
   set(c,'Xda',[xg j j+1],'Yda',[yg x(j) x(j)]);
   if pauslgd >0
     pause(pauslgd/10);
   end
end

plot(x,'.')

for j=1:N
   y(j)=sum(x==j);
end

hold off

figure
set(gcf,'posit',pos+[-200 -200 0 0]) 

title('Andel besök i tillstånd')
xlabel('tillstånd')
ylabel('relativ frekvens')
hold on

r=size(P);
if rank(P-eye(r))~=r(1)-1; % ingen entydig stationär fördelning
  bar(y/(n+1),.35);
  set(gca,'xtick',1:N)  
else
  bar(-0.175+(1:N),y/(n+1),.35);
% beräkning av stationär fördelning
  p0=[zeros(size(P(1,:))) 1];
  P=[P-eye(r) ones(size(P(:,1)))];
  z=p0/P;

  bar(0.175+(1:N),z,0.35,'r');
  set(gca,'xtick',1:N)  
  title('Andel besök i tillstånd (blått), stationär fördelning (rött)')
  ylabel('relativ frekvens, sannolikhet')
end

grid

hold off
  