function Model_1
rand('state',sum(clock));
ILILAB = load('ILILAB.csv');
ILILAB(1,:)=1;
ILILAB = (ILILAB')/sum(ILILAB);

cum_ILILAB = cumsum(ILILAB);

data = load('grouped_data.csv');
data(:,15) = data(:,33);
data(:,6) = data(:,32);
data = data(:,1:15);


para1 = [1 2.5 1 0.05 0 0 0 0 0 0 0 0 0 0 0];


mcmc_n = 11000;

[y1 LL ] = mcmc(data,mcmc_n,para1,cum_ILILAB);


startpt =1001;
endpt = 11000;


csvwrite('mcmc_result_model1.csv',y1);

end

%%%%%%%%%%%%%%%%%%%%%%%%
function result = F_serial_density(x,para)

result = (wblcdf(x+1,para(2),para(1)) - wblcdf(1,para(2),para(1)))/(1-wblcdf(1,para(2),para(1)));
result(result < 0) = 0;
end

%%%%%%%%%%%%%%%%%%%%%
% likelihood function

function result = loglikelihood1(data,para1,cum_ILILAB)

result = -1000000;

if sum(serial_density(1:14,para1(1:2)))>=0.8;
ncolb = max(data(:,3));
ncola = size(data,2)/(ncolb+1);

% the first part of the likelihood

data1 = data(data(:,4) == 1 & data(:,2) ~= 0,:);

a=(para1(1,4))*sum(exp(para1(1,8)*(data1(:,(1:ncolb)*ncola+8)<=18)).* ...
    exp(para1(1,9)*(data1(:,(1:ncolb)*ncola+2)==0)).* ...
    exp(para1(1,10)*(data1(:,(1:ncolb)*ncola+14)==1)).* ...    
    exp(para1(1,11)*(data1(:,(1:ncolb)*ncola+15)==2)).* ...     
    exp(para1(1,12)*(data1(:,(1:ncolb)*ncola+15)==3)).* ...  
    (data1(:,(1:ncolb)*ncola+4)==1).* ...
    serial_density( bsxfun(@minus,data1(:,5),data1(:,(1:ncolb)*ncola+5)),para1(1:2))  ...
    ,2);

c=exp(para1(1,7)*(data1(:,10)==1)+para1(1,6)*(data1(:,8) > 50)+para1(1,5)*(data1(:,8) <= 18)+...
    para1(1,13)*(data1(:,6)==0)+para1(1,14)*(data1(:,6)==2)+para1(1,15)*(data1(:,6)==3) );

aa = a.*c;

bb = para1(1,3)*(cum_ILILAB(data1(:,5)-1)-cum_ILILAB(data1(:,5)-2))'.*c;

L1 = sum(log(1-exp(-aa - bb)));

% the second part of the likelihood
% for all observation.

data2 = data(data(:,2) ~=0,:);

c=exp(para1(1,7)*(data2(:,10)==1)+para1(1,6)*(data2(:,8) > 50)+para1(1,5)*(data2(:,8) <= 18) + ...
    para1(1,13)*(data2(:,6)==0)+para1(1,14)*(data2(:,6)==2)+para1(1,15)*(data2(:,6)==3));

temp_se = F_serial_density( bsxfun(@minus,data2(:,5)-1,data2(:,(1:ncolb)*ncola+5)).*(data2(:,(1:ncolb)*ncola+4)==1),para1(1:2));
temp_se(:,1) = temp_se(:,1) - F_serial_density( data2(:,7),para1(1:2));

a = para1(1,4)*sum( exp(para1(1,8)*(data2(:,(1:ncolb)*ncola+8)<=18)).* ...
    exp(para1(1,9)*(data2(:,(1:ncolb)*ncola+2)==0)).* ...
    exp(para1(1,10)*(data2(:,(1:ncolb)*ncola+14)==1)).* ... 
    exp(para1(1,11)*(data2(:,(1:ncolb)*ncola+15)==2)).* ... 
    exp(para1(1,12)*(data2(:,(1:ncolb)*ncola+15)==3)).* ...    
    (data2(:,(1:ncolb)*ncola+4)==1).*  ...
    temp_se,2);

aa = a.*c;

bb = para1(1,3)*(cum_ILILAB(data2(:,5)-1)-cum_ILILAB(data2(:,11)+data2(:,7)))'.*c;

L2 = - sum(aa) - sum(bb);

result = L1 + L2 ;

end;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function result = logprior1(para)
lik = log(unifpdf(para(1,1),0,1000)) + log(unifpdf(para(1,2),0,1000)) + log(unifpdf(para(1,3),0,1000)) + ...
      log(unifpdf(para(1,4),0,1000)) + log(unifpdf(para(1,5),-1000,1000)) + log(unifpdf(para(1,6),-1000,1000)) + ...
      log(unifpdf(para(1,7),-1000,1000)) + log(unifpdf(para(1,8),-1000,1000)) + log(unifpdf(para(1,9),-0.001,0.001)) + ...
      log(unifpdf(para(1,10),-1000,1000)) + log(unifpdf(para(1,11),-1000,1000)) + log(unifpdf(para(1,12),-1000,1000)) + ...
      log(unifpdf(para(1,13),-1000,1000)) + log(unifpdf(para(1,14),-1000,1000)) + log(unifpdf(para(1,15),-1000,1000));
result = lik;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% main mcmc function

function [result1 LL ] = mcmc(data,mcmc_n,para1,cum_ILILAB)
posterior_para1 = zeros(mcmc_n,length(para1));
LL = zeros(mcmc_n,1);

sigma_para1 = [0.4 0.7 1.8 0.0397 0.88 2.3 1.1 0.8 1 1 1 1 2.1 1 1.3];
posterior_para1(1,:) = para1;



% reshape the data in the big format
data_reshape = reshape_data(data);
LL(1,1) = loglikelihood1(data_reshape,posterior_para1(1,:),cum_ILILAB);

temploglik1 = LL(1,1)+logprior1(posterior_para1(1,:));        

        
for j = 2:mcmc_n;
 


data_reshape = reshape_data(data);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % main mcmc
    % propose new value to paramter first
    proposedvalue1 = normrnd(posterior_para1(j-1,:),sigma_para1,1,length(para1));  
  
    pvalue1 = [proposedvalue1(1,1:2) posterior_para1(j-1,3:length(para1))];
    new = loglikelihood1(data_reshape,pvalue1,cum_ILILAB)+logprior1(pvalue1);
    loglikratio = new-temploglik1;    
    if (loglikratio-log(unifrnd(0,1,1,1))>0);
    posterior_para1(j,1:2)=pvalue1(1,1:2);
    temploglik1 = new;
    else
    posterior_para1(j,1:2)=posterior_para1(j-1,1:2);    
    end;
	
  
    for k = 3:length(para1);
    pvalue1 = [posterior_para1(j,1:(k-1)) proposedvalue1(1,k) posterior_para1(j-1,(k+1):length(para1))];
    new = loglikelihood1(data_reshape,pvalue1,cum_ILILAB)+logprior1(pvalue1);
    loglikratio = new-temploglik1;    
    if (loglikratio-log(unifrnd(0,1,1,1))>0);
    posterior_para1(j,k)=pvalue1(1,k);
    temploglik1 = new;
    else
    posterior_para1(j,k)=posterior_para1(j-1,k);    
    end;
    end;
    
    LL(j,1) = temploglik1-logprior1(posterior_para1(j,:));
	

 
end;

result1 = posterior_para1;
end
%%%%%%%%%%%%%%%%%%%%%%%%
% metropolis-hastings method
function result = metropolis1(data,old_para,new_para,cum_ILILAB)
result = old_para;
L_new = loglikelihood1(data,new_para,cum_ILILAB);
L_old = loglikelihood1(data,old_para,cum_ILILAB); 
    lu = log(unifrnd(0,1,1,1));
    if lu < L_new + logprior1(new_para) - L_old - logprior1(old_para);
        result = new_para;
    end;
end

%%%%%%%%%%%%%%%%%%%%%%
% summary of parameter
function result = para_summary(para,k,l)
[nrow ncol ] = size(para);
mat = zeros(ncol,4);
for i = 1:ncol;
mat(i,:) = [ quantile(para(:,i),[0.5 0.025 0.975]) sum( para(1:(length(para(:,i))-1),i) ~= para(2:length(para(:,i)),i))/(length(para)-1)];   
subplot(k,l,i),plot(para(:,i));
subplot(k,l,i+ncol),hist(para(:,i));
end;
result = mat;
end
%%%%%%%%%%%%%%%%%%%%%%%%
% function to reshape the data so that for every observation, their
% household contact are in the same row

function result = reshape_data(data)

max_member = max(data(:,2));

data1 = [ data ones(size(data,1),(max_member+1)*size(data,2))*-1 ];

for i = 0:max_member;
    data_temp = data(data(:,2) == i,:);
    idx = [];
    idx(cumsum([ 1 (data_temp(:,3)-1)'])) = 1;
    data1(data1(:,3)>=i+1 & data1(:,2) ~= i,size(data,2)*(i+1)+1:(i+2)*size(data,2)) = data_temp(cumsum(idx(1:length(idx)-1)),:);
end;



result = data1;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% function to update the infection time


%%%% parameter
%%% 1: alpha 2:gamma

function result = serial_density(x,para)

x(x == 0) = -1;

result = (wblcdf(x+1,para(2),para(1))-wblcdf(x,para(2),para(1)))/(1-wblcdf(1,para(2),para(1)));



end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%