function [X,T1]=mcsim(Q,n,start,option)
% mcsim(Q,n,start) simulerar markovkedja upp till tid 1,
% Q intensitetssmatris, n antal hopp, start=starttillstånd
% Diagonalen i Q kan anges godtyckligt, programmet räknar ut den korrekta.
% Om starttillstand ej anges ansätts det 1.
% Tillståndsrum {1,2,...,N} om Q NxN-matris.
% Första plotten visar en simulering.
% Om absorberande tillstånd finns, anges förväntad tid till absorption samt
% observerad tid till absorption om simuleringen nått ett sådant tillstånd.
% Om inget absorberande tillstånd finns visar den andra plotten ett
% stolpdiagram med relativa tider i tillstånden. 
% Om en entydig stationär fördelning existerar, läggs denna in i
% detta stolpdiagram. Anropet 
%
% [X,T]=mcsim(Q,n,start) ger vektor med de simulerade tillstånden i
% X och hopptiderna i T. Anropet
%
% mcsim(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, T]=mksim(P,n,start,pauslgd) ger dessutom vektor med de simulerade 
% tillstånden och hopptiderna. Anropet 
%
% [X, T]=mcsim(Q,n,start,fig) där fig är negativt tal ger vektorer med
% de simulerade tillstånden och hopptiderna, men ej några figurer.


% Jan Enger, KTH,  ht 2001

[N, N1]=size(Q); 

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

Q=Q-diag(diag(Q));

if any(Q(:)<-10^(-12)) | N~=N1
  disp('       Q ej intensitetsmatris');
  return;
end

d=sum(Q');

abs=0;

if any(d==0) 
  abs=1;
end

d1=d;
d1((d==0))=1+d(d==0);
abst=find(d==0);

N=length(d);

P=Q./(d1'*ones(1,N)); 			% hoppmatris
for j=abst
  P(j,j)=1;
end

Q=Q-diag(d);

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

if d(x(1))==0
   disp('   start i absorberande tillstånd');
   return
end

t(1)=-log(rand(1))/d(x(1));

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
  if d(x(k))==0
    abs=2;
    t(k)=0.1*sum(t);
    n=k;
    break
  end 
  t(k)=-log(rand(1))/d(x(k)); 		% tid i tillståndet
  
end

T=cumsum(t);

if nargout>=1
  X=x;
  T1=T;
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([0 t(1)],[x(1) x(1)])
axis([0 1.1*sum(t) 0 N+1])
c=get(gca,'ch');
set(gca,'ytick',1:N)

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

hold on

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

plot([0 T(1:n-1)],x ,'.')

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

switch abs

case 0

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


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


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),y1/sum(t),.35);
  set(gca,'xtick',1:N)
% beräkning av stationär fördelning
a=[zeros(1,N) 1];
Q1=[Q ones(N,1)];
z=a/Q1;

  bar(0.175+(1:N),z,0.35,'r');
  title('Relativ tid i tillstånd (blått), stationär fördelning (rött)')
  ylabel('relativ tid, sannolikhet')
end

grid

hold off

case 1

text(T(n)*.5,0.2*N,'Processen har ej nått något absorberande tillstånd')

case 2

plot([T(n-1) T(n-1)],[0 .5*x(n)]);

abstext=sprintf('observerad tid till absorbtion %3.3g',T(n-1));

text(0.5*T(n),0.1*N,abstext);


a=find(d==0);
Q(a,:)=[];
Q(:,a)=[];

ta=-inv(Q)*ones(N-length(a),1);

d(find(d~=0))=ta;

ta=d(start);

abstext=sprintf('förväntad tid till absorbtion %3.3g',ta);

text(0.5*T(n),0.2*N, abstext);


end


