clc;clear;
aa = [0.01,0.02,0.01,0.03,0.01];
bb = [0.1,0.2,0.3,0.2,0.1,0.2];
nn = [5,5,4,5,3];
kk = [2,2,2,3,1];
As=1;
for j = 1:length(aa)
a = aa(j); b = bb(j); n = nn(j); k = kk(j);
disp('=================================')
t=30;c=0.1;
Q=zeros(n-k+2);
for i=1:n-k+1
Q(i,i+1)=k*a;
end
% disp(Q)
for i=1:n-k+1
Q(i+1,i)=i*b;
end
% disp(Q)
for i=1:n-k+2
Q(i,i)=-sum( Q(i,:) );
end
disp(Q)
la=1.5*max(-diag(Q));
P=eye(n-k+2)+Q/la;
%%
error=1;
M=50;
while error>0.2
M=M+1;
%round(la*t+(la*t)^0.5*c);
error=1-poisscdf(M,la*t);
end
error
%%
res_sum=0;
for n=0:M
res_sum=res_sum+exp(-la*t)*(la*t)^n/factorial(n)*(eye(1,size(P,1))*P^n);
end
res_sum
%% t ok
% Q=Q';
% q=size(Q,1)+1;
% Q(q,1:q)=1;
% rref(Q)
%%
As=As*(1-res_sum(length(res_sum)));
end
% for i = 1:length(a)
% for j = 1:length(b)
% for l = 1:length(n)
% for m = 1:length(k)
%
% sum = calc_sum(a(i),b(j),n(l),k(m))
% end
% end
% end
% end