%% 
clear all; 
clc;
K = 6; 
Lamda = 1; 
% K = 1./(Lambda-t); 
tt = Lamda -(1/K); 
mmt= Lamda./(Lamda-tt);
Num = 10000;
Random = rand(1,Num)'; 
% Random = u;

Gx1 = zeros(Num,1);
tic 
%% classical estimate 
xx1(1:Num) = -log(Random(1:Num));
% xx1(1:Num) = -log(u(1:Num));
% X1 = -log(Z);
Gx1(xx1(1:Num)> K) = 1;
Gx1(xx1(1:Num)<=K) = 0;
ThetaHat(1:Num,1) = 1./(1:Num)' .* cumsum(Gx1(:));
cclassical = toc; 
% thetahat(i) 
plot(ThetaHat) 
% find(ThetaHat' ~=thetahat)
%% tilted sampling 
xx(1:Num) = -K .* log(Random(1:Num));
% xx(1:Num) = -K .* log(Random(1:Num));
Gx(xx(1:Num)>= K) = exp(-tt.*xx(xx(1:Num)>= K));
Gx(xx(1:Num)<  K) = 0;
ThetaTilda(1:Num,1) = mmt./(1:Num)' .* cumsum(Gx(:));
tilted = toc;
% thetatilda(i) ;
plot(ThetaTilda) 
% find(ThetaTilda' ~=thetatilda)
% classical
% tilted
%%

%%

%% Commands ¡V Part 4(b) 
clear all; clc; 
% aa = 6; 
aa(1:4) = 2*(1:4) ;
Num = 10000;
Gx1 = zeros(Num,4);
Lamda = 1; 

Random = rand(1,Num)'; 
for i =1:4
tt = Lamda -(1/aa(i)); 
mmt= Lamda./(Lamda-tt);
%% classical estimate 
% xx1(1:Num) = -log(Random(1:Num));
xx1(1:Num) = -log(Random(1:Num));
% X1 = -log(Z);
Gx1(xx1(1:Num)>aa(i),aa(i)) = 1;
Gx1(xx1(1:Num)<=aa(i),aa(i)) = 0;
ThetaHat(1:Num,aa(i)) = 1./(1:Num)' .* cumsum(Gx1(:,aa(i)));
plot(ThetaHat) 
%% Tilted Sampling 
xx(1:Num) = -aa(i) .* log(Random(1:Num));
Gx(xx(1:Num)>= aa(i),aa(i)) = exp(-tt.*xx(xx(1:Num)>= aa(i)));
Gx(xx(1:Num)<  aa(i),aa(i)) = 0;
ThetaTilda(1:Num,aa(i)) = mmt./(1:Num)' .* cumsum(Gx(:,aa(i)));
plot(ThetaTilda)
%%

end
plot(ThetaHat(:,aa(1:4))) ;
 
display(ThetaHat(Num,aa(1:4))) ;
  
plot(ThetaTilda(:,aa(1:4))) 

display(ThetaTilda(Num,aa(1:4))); 
 
display(ThetaHat(Num,aa(1:4))); 
