function ret = lecture4B_sample_RVMean_Var_Entropy(task, with_replace_in, k_in)
% TO RUN THIS: uncomment the top section-Setup for the simulation you want to run (only one section setup should be active at a time) 
global N;
global D;
global DATA;
global SAMPLED;
global SAMPLEDCOUNT;
global DIMSAM;
global TITLETEXT;
global WITH_REPLACE_OR_NOT;
global FIGNOW
global XLIM
global LASTK

LASTK =k_in;
% h =  findobj('type','figure');
% FIGNOW = length(h)+1;

XLIM=-1;
if (	strcmp(with_replace_in, 'with_replace'))
	WITH_REPLACE_OR_NOT=1;
else
	WITH_REPLACE_OR_NOT=0;
end;
%%%%%%%%% ONE DIE ROLL
if (strcmp (task,'one die roll uniform'))
	D=1; N=6;
	PROBABS = ones(6,1)*1/6;
	DATA = [PROBABS  [1:6]']
	DIMSAM=2;
	TITLETEXT= 'ONE DIE ROLLS';
	FIGNOW =1;
	XLIM=[1 6];
end;
%%%%%%%%%% ONE DIE ROLL NONUNIFORM-6
if (strcmp (task,'one die roll nonuniform-6'))
	D=1; N=6;
	PROBABS = [0.1 0.1 0.1 0.1 0.1 0.5]';
	DATA = [PROBABS  [1:6]']
	DIMSAM=2;
	TITLETEXT= 'ONE DIE ROLL NONUNIFORM-6';
	FIGNOW =1;
	XLIM=[1 6];
end;
%%%%%%%%%% ONE DIE ROLL EXTREME-1-6
if (strcmp (task,'one die roll nonuniform-1-6'))
	D=1; N=6;
	PROBABS = [0.3 0.1 0.1 0.1 0.1 0.3]';
	DATA = [PROBABS  [1:6]']
	DIMSAM=2;
	TITLETEXT= 'ONE DIE ROLL EXTREME-1-6';
	FIGNOW =2;
	XLIM=[1 6];
end;
% %%%%%%%%%% ONE DIE ROLL MIDDLE-3-4
if (strcmp (task,'one die roll nonuniform-3-4'))
	D=1; N=6;
	PROBABS = [0.05 0.05 0.4 0.4 0.05 0.05]';
	DATA = [PROBABS  [1:6]']
	DIMSAM=2;
	TITLETEXT= 'ONE DIE ROLL MIDDLE-3-4';
	FIGNOW =3;
	 XLIM=[1 6];
end;
%%%%%%%%%% TWO DICE ROLL
if (strcmp (task,'two dice roll uniform'))
	D=1; N=36;
	PROBABS = ones(N,1)*1/N;
	DATA = [PROBABS  [2:7 3:8 4:9 5:10 6:11 7:12]']
	DIMSAM=2;
	TITLETEXT= 'TWO DIE ROLLS';
	FIGNOW =1;
	 XLIM=[2 12];
end;
 
 % %%%%%%%%% 1000 CAR PRICES
 if (strcmp (task,'1000 car prices'))
	 D=1;
	 N=1000;
	 PROBABS = ones(N,1)*1/N;
	 DATA = [PROBABS  normrnd(14000,1000,N,1)];
	 DIMSAM=2;
	 TITLETEXT= '1000 CAR PRICES';
	 FIGNOW =1;
	 XLIM=[0 30000];
end;
%%%%%%%%%% COIN TOSS
if (strcmp (task,'coin toss fair'))
	D=1; N=2;
	PROBABS = ones(N,1)*1/N;
	DATA = [PROBABS  [0:1]']
	DIMSAM=2;
	TITLETEXT= 'COIN TOSS';
	FIGNOW =1;
	 XLIM=[0 1];
 end;
%%%%%%%%%% COIN TOSS NONUNIFORM
if (strcmp (task,'coin toss heads'))
	D=1; N=2;
	PROBABS = [0.3 0.7]';
	DATA = [PROBABS  [0:1]']
	DIMSAM=2;
	TITLETEXT= 'COIN TOSS';
	FIGNOW =2;
	 XLIM=[0 1];
end;
%%%%%%%%% 600 uniform random integers 1-100
if (strcmp (task,'600 uniform integers'))
	D=1;
	N=600;
	PROBABS = ones(N,1)*1/N;
	DATA = [PROBABS  randi(100,N,D)]
	DIMSAM=2;
	TITLETEXT= '600 RANDOM INTEGERS UPTO 100';
	FIGNOW =1
	XLIM=[1 100];
end;
%%%%%%%%%% 17F MIDTERM GRADES
if (strcmp (task,'midterm 17F'))
	grades = importdata('17F_midterm_grades.txt')
	[N, D]=size(grades)
	PROBABS = ones(N,1)*1/N;
	DATA = [PROBABS  grades]
	DIMSAM=2;
	TITLETEXT= '17F MIDTERM GRADES';
	FIGNOW =2;
	 XLIM=[0 100];
end;
%%%%%%%%%% 600 INTEGERS CLOSE-MEAN-76
if (strcmp (task,'600 integers close-mean-76'))
	D=1;
	N=600;
	PROBABS = ones(N,1)*1/N;
	DATA = [PROBABS  normrnd(76,4,N,1)];
	DIMSAM=2;
	TITLETEXT= '600 INTEGERS CLOSE-MEAN-76';
	FIGNOW =1;
	XLIM=[0 100];
end;

%%%%%%%%% 600 INTEGERS EXTREME
if (strcmp (task,'600 integers extremes'))
	D=1;
	N=600;
	PROBABS = ones(N,1)*1/N;
	FB=2;
	DATA = [PROBABS  round(100*random('Beta',0.76*FB,0.24*FB,N,1)) ];
	DIMSAM=2;
	TITLETEXT= '600 INTEGERS EXTREME';
	FIGNOW =1;
	XLIM=[0 100];
end;



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
SAMPLEDCOUNT = zeros(N,1);
SAMPLED = [];

for T=1:1000
	if (WITH_REPLACE_OR_NOT==1) 
		[iseq, seq] = sample_with_replacement(LASTK);
	else
		[iseq, seq] = sample_without_replacement(LASTK);
	end;
	seq_sampled = seq	
	plot_sampled();
	pause
end;
return;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [isamp,samp] = sample_with_replacement(k)
global N;
global D;
global DATA;
global SAMPLED;
global SAMPLEDCOUNT;
global DIMSAM;

isamp = zeros(1,k);
samp = -ones(1,k);
probscumul = cumsum(DATA(:,1));
for t=1:k
	v = rand();
	p = 0; q=N; m = round((p+q)/2);
	while(m<q)
		if(v<=probscumul(m)) q=m; else p=m; end;	
		m = round((p+q)/2);	
	end;
	isamp(t)=q;
	samp(t) = DATA(q,DIMSAM);
	SAMPLED = [SAMPLED ; DATA(q,DIMSAM)];
	SAMPLEDCOUNT(q) = SAMPLEDCOUNT(q) + 1;
end;
	 
return;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [isamp,samp] = sample_without_replacement(k)
global N;
global D;
global DATA;
global SAMPLED;
global SAMPLEDCOUNT;
global DIMSAM;
global TITLETEXT;

% available = find(SAMPLEDCOUNT==0);
% permavailable = available (randperm(length(available)));
% isam = permavailable(1: min(k,length(permavailable)));
isamp = [];
samp = [];

for t=1:k
	availind= find(SAMPLEDCOUNT==0);
	if (length(availind)<1) break; end;
	probscumul = cumsum(DATA(availind,1) / sum(DATA(availind,1)) );
	v = rand();
	p = 0; q=length(availind); m = round((p+q)/2);
	while(m<q)
		if(v<=probscumul(m)) q=m; else p=m; end;	
		m = round((p+q)/2);	
	end;
	q = availind(q)
	isamp(t)=q;
	samp(t) = DATA(q,DIMSAM);
	SAMPLED = [SAMPLED ; DATA(q,DIMSAM)];
	SAMPLEDCOUNT(q) = SAMPLEDCOUNT(q) + 1;
end;
	 
return;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function []= plot_sampled()
global N;
global D;
global DATA;
global SAMPLED;
global SAMPLEDCOUNT;
global TITLETEXT;
global WITH_REPLACE_OR_NOT;
global FIGNOW;
global XLIM;
global LASTK;

figure(FIGNOW); clf;
% histogram(SAMPLED(SAMPLED>0))

% histogram(SAMPLED,20);
[C1,BINE] = histcounts (SAMPLED,20);
C2 = histcounts (SAMPLED(1:end-LASTK),BINE);
C3 = histcounts (SAMPLED(end-LASTK+1:end),BINE);
BIN= (BINE(1:end-1) + BINE(2:end) )/2;
% C2
% C3
% bar(BIN, [C2 ; C3]','stacked' );
 bar(BIN,C1);
grid on; hold on;
ax=gca; 
fy = max  (1,  max(ax.YLim) /10);
bar(BIN,C3*fy,'g');

% ax.YLim = [0 max(50,max(SAMPLEDCOUNT))];
ax.XLim = [0 max(SAMPLED)];
if (max(XLIM)>0) ax.XLim = XLIM; end;
% xtickformat('%.d')
ax.XAxis.Exponent = 0;

MU = mean(SAMPLED);
SSTD = std(SAMPLED);
CS = sum(SAMPLEDCOUNT);
% plot(MU, 0 ,'r*');
line([MU MU],ax.YLim/10, 'Color','red','linewidth',4 )
if (WITH_REPLACE_OR_NOT) w_rep_ornot_text='W/ REPLACE'; else w_rep_ornot_text='WITHOUT REPLACE';end; 
set(gca, 'FontSize', 16);
 set(gca,'units','normalized','position',[0.06 0.12 0.89 0.8]);
title(sprintf('%s  sampling histogram, %s   samplesize=%d   mean=%.3f   std=%.3f ',TITLETEXT,w_rep_ornot_text, CS,MU,SSTD),'fontsize',18 );
ylabel('count (frequency in bucket) ','fontsize',18,'FontWeight','bold');
xlabel('values sampled','fontsize',18,'FontWeight','bold');
return;