Symbols and bar thickness in BAR graphs

Hi,

I am using the below code to generate bar graphs. I am not able to get the symbols as typed in the latex commands with subscripts. The symbols appear as they are typed in the Matlab code and I do not see the subscripts. Can someone please suggest how to make the symbols with subscripts appear on the y axis. Can someone suggest how to reduce the bar thickness.

%%%%%%%%%%%%%%%%%

tic % All values are for DDT and p=0.04, t=90 days clear all clear memory Nh=33898857; Nc=21571585; H=7933615; % beta=4.6; % fixed number cattle per cattle shed p=0.093; % p = percent increase in natural sand fly death rate t=90;% model is run with time after spray as 90 days % Type of Analysis % analysis=1; % For Uncertainity analysis analysis=2; % For Sensitivity analysis

if analysis==1 dim=1e5; % dim=3; end

if analysis==2 dim=1e5; % dim=1e4; % dim=5; end

model='Sand flies killed'; % examine the effect of Sand flies killed with different ...

    % response from my model is the optimal values of x and y 
    x_rlz=zeros(dim,1); % a zero matrix called "R0 realizations" of dim rows and reali columns
    y_rlz=zeros(dim,1);

     % 6 input parameters for getting the above response
    beta_all=zeros(dim,1);
    ah_all=zeros(dim,1);
    muv_all=zeros(dim,1);
    Ct0_all=zeros(dim,1); % DDT
%     Ct0_all=zeros(dim,reali); % Deltamethrin
    b1_all=zeros(dim,1);
    b2_all=zeros(dim,1);

    beta=zeros(dim,1);
    ah=zeros(dim,1);
    muv=zeros(dim,1);
    Ct0=zeros(dim,1);
    b1=zeros(dim,1);
    b2=zeros(dim,1);

    beta_vec=zeros(1,6);
    ah_vec=zeros(1,6);
    muv_vec=zeros(1,6);
    Ct0_vec=zeros(1,6);
    b1_vec=zeros(1,6);
    b2_vec=zeros(1,6);

    x_vec=zeros(1,6);
    y_vec=zeros(1,6);

    tag=zeros(dim,1);

      name1='exponential';
      name2='beta';
      name3='gamma'; % syntax of random function, A: shape parameter  B: scale parameter
      name4='uniform';
      name5='normal';
      name6='gev'; % generalized extreme value distribution
      name7='logn'; % Log normal distribution

% I/P parameter, ah has Normal Distribution, mean is 0.1792 and variance is 0.00396 par_ah=0.1792; % var_mu=0.00396344; std_ah=0.028; ah=random(name5,par_ah,std_ah,dim,1); % VECTOR

% THIS BELOW CODE MAKES IT A TRUNCATED NORMAL distribution l=0; % to the letter small "l" number zero is assigned for k1=1:dim if ah(k1) < (0.0001) l=l+1; tag(l)=k1; % assigns the row index of negative values to vector tag end end

%%%%%%added lines to check while min(ah)<0 for k1=1:dim if ah(k1) < 0 %(0.0001) ah(k1)=random(name5,par_ah,std_ah,1,1); end end end

ah_all(:,1)=ah; % to the j th column of array muv, the column vector muv is assigneds

% Below Loop introduced by Kaushik on 9_14_2012 to check if negative values do occur in the column vector ah, even after executing the % above two loops if min(ah_all)<0 fprintf(2,['One of ah values is negative ','\n']); end

      ah_vec(1,1)=mean(ah);
      ah_vec(1,2)=median(ah);
      ah_vec(1,3)=std(ah);
      ah_vec(1,4)=var(ah);
      ah_vec(1,5)=min(ah);
      ah_vec(1,6)=max(ah);

      ac=1-ah;       %  VECTOR
       Q=(Nh.*ah)./(ah.*Nh+ac*Nc);           %  VECTOR

       % I/P parameter muv also has a normal distribution just as ah above
%       par_muv = 1/(0.47*30);  std_muv = 1/(0.42*30); %         var_mu=0.00396344;
%       ShapeParameter_muv= 23.611; ScaleParameter_muv=0.00321387;
%       muv=random(name5,par_muv,std_muv,dim,1);                 %  VECTOR
% low_muv=1/(0.47000000000001*30);high_muv=1/(0.46999999999999*30); % from
% first paper of srinivasan
low_muv=0.07499999999999;high_muv=0.07511111111111;
muv=random(name4,low_muv,high_muv,dim,1);  
%      muv=random(name3,ShapeParameter_muv,ScaleParameter_muv,dim,1);                 %  VECTOR

% % % % % % % % To have all non-negative values of muv we choose again values which % % % % % % % % are negative % % % % % % % l=0; % to the letter small "l" number zero is assigned % % % % % % % for k1=1:dim % % % % % % % if muv(k1) < 0 %(0.0001) % % % % % % % % k1 % % % % % % % % mu(k1) % % % % % % % l=l+1; % % % % % % % tag(l)=k1; % % % % % % % end % % % % % % % end % % % % % % % % % % % % % % while min(muv)<0 % % % % % % % for k1=1:dim % % % % % % % if muv(k1) < 0 %(0.0001) % % % % % % % % muv(k1)=random(name3,par_muv,std_muv,1,1); % % % % % % % muv(k1)=random(name5,par_muv,std_muv,1,1); % % % % % % % end % % % % % % % end % % % % % % % end

      muv_all(:,1)=muv;      % to the j th column of array muv, the column vector muv is assigneds
if min(muv_all)<0
     fprintf(2,['One of muv values is negative  ','\n']); 
end

      muv_vec(1,1)=mean(muv);
      muv_vec(1,2)=median(muv);
      muv_vec(1,3)=std(muv);
      muv_vec(1,4)=var(muv);
      muv_vec(1,5)=min(muv);
      muv_vec(1,6)=max(muv);

 % Uniform distribution
%         low_Ct0= 0.4941;  high_Ct0=0.5858;                
%         Ct0=random(name4,low_Ct0,high_Ct0,dim,1); %  assuming uniform distribution
%         par_Ct0=(low_Ct0+high_Ct0)/2; var_Ct0 = (high_Ct0-low_Ct0)^2/12; std_Ct0 = sqrt(var_Ct0);
       Mean_Ct0=0.54; Std_Ct0=0.0265;
       Ct0=random(name5,Mean_Ct0,Std_Ct0,dim,1);

    Ct0_all(:,1)=Ct0;

      Ct0_vec(1,1)=mean(Ct0);
      Ct0_vec(1,2)=median(Ct0);
      Ct0_vec(1,3)=std(Ct0);
      Ct0_vec(1,4)=var(Ct0);
      Ct0_vec(1,5)=min(Ct0);
      Ct0_vec(1,6)=max(Ct0);

      % To have all non-negative values of muv we choose again values which
      % are negative
      l=0;            % to the letter small "l" number zero is assigned
     for k1=1:dim
         if Ct0(k1) < 0  %(0.0001)
    %         k1
    %         mu(k1)
             l=l+1;
             tag(l)=k1;
         end
     end

     while min(Ct0)<0
        for k1=1:dim 
         if Ct0(k1) < 0  %(0.0001)
            Ct0(k1)=random(name5,Mean_Ct0,Std_Ct0,1,1);    
         end 
        end
     end

% if min(Ct0_all)<low_Ct0 % fprintf(2,['Lowest sampled Ct0 value is lesser than low_Ct0 ','\n']); % end % if min(Ct0_all)>high_Ct0 % fprintf(2,['Highest sampled Ct0 value is greater than high_Ct0 ','\n']); % end

      %%%%%%
% % %       % Uniform distribution
% % %         low_b1 = 0.0086794; high_b1= 0.04332; b1=random(name4,low_b1,high_b1,dim,1);           % 
% % %         par_b1=0.026; %par_b1=(low_b1+high_b1)/2;

% K_b1=0.741502;% shape parameter % mu_b1=0.00664696; % ?: location parameter % sigma_b1=0.00351183; % ?: scale parameter % b1=random(name6,K_b1,mu_b1,sigma_b1,dim,1); % assuming generalized extreme value distribution % Gamma distribution % ShapeParameter_b1 = 2.30441; ScaleParameter_b1=0.00527972; % b1=random(name3,ShapeParameter_b1,ScaleParameter_b1,dim,1);

% % % % % % % Uniform Distribution % % % % % % low_b1 = 0; high_b1= 0.032; b1=random(name4,low_b1,high_b1,dim,1); % % % % % % % par_b1=0.012; %par_b1=(low_b1+high_b1)/2;

% standard uniform distribution obtained by setting both parameters of beta % distribution to 1 StandardUniform_b1=random(name2,1,1,dim,1); for i=1:dim if (0 < StandardUniform_b1(i)) && (StandardUniform_b1(i) <= 1/6) b1(i)=0.004; elseif (1/6 < StandardUniform_b1(i)) && (StandardUniform_b1(i) <= 2/6) b1(i)=0.006; elseif (2/6 < StandardUniform_b1(i)) && (StandardUniform_b1(i) <= 3/6) b1(i)=0.007; elseif (3/6 < StandardUniform_b1(i)) && (StandardUniform_b1(i)<= 4/6) b1(i)=0.009; elseif (4/6 < StandardUniform_b1(i)) && (StandardUniform_b1(i)<= 5/6) b1(i)=0.019; elseif (5/6 < StandardUniform_b1(i)) && (StandardUniform_b1(i)<= 6/6) b1(i)=0.028; end end b1_all(:,1)= b1;

      b1_vec(1,1)=mean(b1);
      b1_vec(1,2)=median(b1);
      b1_vec(1,3)=std(b1);
      b1_vec(1,4)=var(b1);
      b1_vec(1,5)=min(b1);
      b1_vec(1,6)=max(b1);
      %%%%%%

% if min(b1_all)<low_b1 % fprintf(2,['Lowest sampled b1 value is lesser than low_b1 ','\n']); % end % if min(b1_all)>high_b1 % fprintf(2,['Highest sampled b1 value is greater than high_b1 ','\n']); % end

       %%%%%%
% % % %       % Uniform distribution
% % % %         low_b2 = 0.026358; high_b2= 0.095641; b2=random(name4,low_b2,high_b2,dim,1);  
% % % %         par_b2=0.061; %par_b2=(low_b2+high_b2)/2; %var_b2 = (high_b2-low_b2)^2/12;
% % % %         std_b2 =0.02; %sqrt(var_b2);
%         mu_b2 = -4.64085; sigma_b2 = 0.742362;
%         b2=random(name7,mu_b2,sigma_b2,dim,1);        % taking lognormal distribution

% ShapeParameter_b2= 3.03234; ScaleParameter_b2= 0.026712; % taking Gamma distribution % b2=random(name3,ShapeParameter_b2,ScaleParameter_b2,dim,1);

% % % % % % % % Uniform Distrbution % % % % % % % low_b2 = 0.002; high_b2= 0.21; b2=random(name4,low_b2,high_b2,dim,1); % % % % % % % par_b2=0.081; %par_b2=(low_b2+high_b2)/2; %var_b2 =(high_b2-low_b2)^2/12;

% standard uniform distribution obtained by setting both parameters of beta % distribution to 1 StandardUniform_b2=random(name2,1,1,dim,1); for i=1:dim if (0 < StandardUniform_b2(i)) && (StandardUniform_b2(i) <= 1/6) b2(i)=0.0280; elseif (1/6 < StandardUniform_b2(i)) && (StandardUniform_b2(i) <= 2/6) b2(i)=0.0440; elseif (2/6 < StandardUniform_b2(i)) && (StandardUniform_b2(i) <= 3/6) b2(i)=0.0670; elseif (3/6 < StandardUniform_b2(i)) && (StandardUniform_b2(i)<= 4/6) b2(i)=0.0720; elseif (4/6 < StandardUniform_b2(i)) && (StandardUniform_b2(i)<= 5/6) b2(i)=0.0910; elseif (5/6 < StandardUniform_b2(i)) && (StandardUniform_b2(i)<= 6/6) b2(i)=0.1840; end end

b2_all(:,1)= b2;

      b2_vec(1,1)=mean(b2);
      b2_vec(1,2)=median(b2);
      b2_vec(1,3)=std(b2);
      b2_vec(1,4)=var(b2);
      b2_vec(1,5)=min(b2);
      b2_vec(1,6)=max(b2);

% if min(b2_all)<low_b2 % fprintf(2,['Lowest sampled b2 value is lesser than low_b2 ','\n']); % end % if min(b2_all)>high_b2 % fprintf(2,['Highest sampled b2 value is greater than high_b2 ','\n']); % end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % I/P parameter beta, has gamma Distribution, mean is 4.6 and s.d. is 2.6 % % % % % % par_beta=4.6; std_beta=2.6; % var_mu=0.00396344; % % % % % % % beta=0.0001+random(name1,par_beta,[dim,1]); %Anuj said add 0.0001 % % % % % % ShapeParameter_beta=8.1385; ScaleParameter_beta=0.5652; low_beta=4.59999999999999; high_beta=4.600000000000010; % % % % % % beta=0.0001+random(name3,ShapeParameter_beta,ScaleParameter_beta,[dim,1]); beta=random(name4,low_beta,high_beta,dim,1); beta_all(:,1)=beta;

% % % % % % l=0; % to the letter small "l" number zero is assigned % % % % % % for k1=1:dim % % % % % % if beta(k1)<(0.0001) % % % % % % l=l+1; % % % % % % tag(l)=k1; % % % % % % end % % % % % % end % % % % % % % % % % % % while min(beta)<0 % % % % % % for k1=1:dim % % % % % % if beta(k1) < 0 %(0.0001) % % % % % % beta(k1)=random(name3,par_beta,1); % % % % % % end % % % % % % end % % % % % % end

      beta_vec(1,1)=mean(beta);
      beta_vec(1,2)=median(beta);
      beta_vec(1,3)=std(beta);
      beta_vec(1,4)=var(beta);
      beta_vec(1,5)=min(beta);
      beta_vec(1,6)=max(beta);

if min(beta_all)<0 fprintf(2,['One of beta values is negative ','\n']); end

       Z=Nc./beta;       %  VECTOR
%%%%%%%%%%%%%%%formula for x and y in terms of the 6 independent variables

L1=(Q.*Ct0).*(exp(-b1.*t)); L2=((1-Q).*Ct0).*(exp(-b2.*t)); Ih=0.533; % for DDT K1=Ih*7933615; %K1=Ih*H; Iz for DDT is 0.533 % only K1 is a scalar K2=(Ih/2).*Z; %K2=Iz*Z;

kk1=L1/K1; kk2=L2./K2; kk3=p.*muv; kk4=L1+L2; % count=zeros(1,4);

% defining 12 dynamic arrays below count1=[];count2=[]; count3=[]; count4=[];

for i=1:dim if (kk2(i) < kk1(i)) && (kk3(i) <= L1(i)) x_rlz(i)= K1*p*muv(i)/(Nh*L1(i)); y_rlz(i)=0; count1=[count1 i]; elseif ( kk2(i) < kk1(i) ) && ( L1(i) < kk3(i) ) && ( kk3(i) <= kk4(i) ) x_rlz(i)= K1/Nh; y_rlz(i)= K2(i)*(p*muv(i)-L1(i))/(Nc*L2(i)); count2=[count2 i]; elseif ( kk1(i) < kk2(i) ) && ( kk3(i) <= L2(i) ) x_rlz(i)= 0; y_rlz(i)= K2(i)*p*muv(i)/(Nc*L2(i)); count3=[count3 i]; elseif (kk1(i) < kk2(i)) && ( L2(i) < kk3(i)) && ( kk3(i) <= kk4(i) ) x_rlz(i)= K1*(p*muv(i)-L2(i))/(Nh*L1(i)); y_rlz(i)= K2(i)/Nc; count4=[count4 i]; end end clear i % ################################################################# % Sensitivity Analysis starts here........

if analysis==2

% IS below the rank transformation of the data

% below is for x when condition 1 is satisfied, so i/p & o/p parameter % values are extracted when condition 1 is satisfied by using count1 % try these commands to understand how extraction works: % >> x=randn(10,1) % >> y=[4 5] % >> x(y)

% Start:code with count1: sort & order matrices: condition 1 satisfied [sorted_beta_1, order_beta_1] = sort(beta_all(count1,1)); % sorting the beta matrix in ascending order [sorted_ah_1, order_ah_1] = sort(ah_all(count1,1)); [sorted_muv_1, order_muv_1] = sort(muv_all(count1,1)); [sorted_Ct0_1, order_Ct0_1] = sort(Ct0_all(count1,1)); [sorted_b1_1, order_b1_1] = sort(b1_all(count1,1)); [sorted_b2_1, order_b2_1] = sort(b2_all(count1,1)); [sorted_x_1, order_x_1]=sort(x_rlz(count1,1)); % y (R0) is rank transformed [sorted_y_1, order_y_1]=sort(y_rlz(count1,1));

for i=1:size(count1,2) beta_rank_1(order_beta_1(i))=i; ah_rank_1(order_ah_1(i))=i; muv_rank_1(order_muv_1(i))=i; Ct0_rank_1(order_Ct0_1(i))=i; b1_rank_1(order_b1_1(i))=i; b2_rank_1(order_b2_1(i))=i; x_rank_1(order_x_1(i))=i; y_rank_1(order_y_1(i))=i; end % End:code using count 1: to sort & order matrices: condition 1 satisfied

% ################################################################# % Below start:i/p param ah & o/p param x1:condition 1(use count1 OR underscore 1)satisfied X_x_ah_1 = [ones(size(beta_rank_1')) beta_rank_1' muv_rank_1' Ct0_rank_1' b1_rank_1' b2_rank_1'];

[b_x_ah_1, bint_x_ah_1, r_x_ah_1]= regress(x_rank_1', X_x_ah_1); % Removes NaN data [b_ah_x_1, bint_ah_x_1, r_ah_x_1]= regress(ah_rank_1', X_x_ah_1);

[prcc_x_ah_1, p_val_x_ah_1]= corrcoef(r_ah_x_1, r_x_ah_1);

prcc(1)=prcc_x_ah_1(1,2); p_val(1)=p_val_x_ah_1(1,2); % ################################################################# % Above end:i/p param ah & o/p param x1:condition 1 (use count1)satisfied

% ################################################################# % Below start:i/p param muv & o/p param x1:condition 1 (use count1)satisfied X_x_muv_1 = [ones(size(beta_rank_1')) beta_rank_1' ah_rank_1' Ct0_rank_1' b1_rank_1' b2_rank_1'];

[b_x_muv_1, bint_x_muv_1, r_x_muv_1]= regress(x_rank_1', X_x_muv_1); % Removes NaN data [b_muv_x_1, bint_muv_x_1, r_muv_x_1]= regress(muv_rank_1', X_x_muv_1);

[prcc_x_muv_1, p_val_x_muv_1]= corrcoef(r_muv_x_1, r_x_muv_1);

prcc(2)=prcc_x_muv_1(1,2); p_val(2)=p_val_x_muv_1(1,2); % ################################################################# % Above end:i/p param muv & o/p param x1:condition 1 (use count1)satisfied

% ################################################################# % Below start:i/p param Ct0 & o/p param x1: condition 1 (use count1)satisfied X_x_Ct0_1 = [ones(size(muv_rank_1')) beta_rank_1' ah_rank_1' muv_rank_1' b1_rank_1' b2_rank_1']; % X_y_Ct0_1 = [ones(size(muv_rank_1')) beta_rank_1' ah_rank_1' muv_rank_1' b1_rank_1' b2_rank_1'];

[b_x_Ct0_1, bint_x_Ct0_1, r_x_Ct0_1]= regress(x_rank_1', X_x_Ct0_1); % Removes NaN data [b_Ct0_x_1, bint_Ct0_x_1, r_Ct0_x_1]= regress(Ct0_rank_1', X_x_Ct0_1);

[prcc_x_Ct0_1, p_val_x_Ct0_1]= corrcoef(r_Ct0_x_1, r_x_Ct0_1);

prcc(3)=prcc_x_Ct0_1(1,2); p_val(3)=p_val_x_Ct0_1(1,2); % ################################################################# % Above end: i/p param Ct0 & o/p param x1:condition 1 (use count1)satisfied

% ################################################################# % Below start: i/p param b1 & o/p param x1 when condition 1 (use count1)satisfied X_x_b1_1 = [ones(size(beta_rank_1')) beta_rank_1' muv_rank_1' Ct0_rank_1' ah_rank_1' b2_rank_1'];

[b_x_b1_1, bint_x_b1_1, r_x_b1_1]= regress(x_rank_1', X_x_b1_1); % Removes NaN data [b_b1_x_1, bint_b1_x_1, r_b1_x_1]= regress(b1_rank_1', X_x_b1_1);

[prcc_x_b1_1, p_val_x_b1_1]= corrcoef(r_b1_x_1, r_x_b1_1);

prcc(4)=prcc_x_b1_1(1,2); p_val(4)=p_val_x_b1_1(1,2); % ################################################################# % Above end: i/p param b1 & o/p param x1:condition 1(use count1 OR % underscore 2)satisfied

% Start code using count2 to sort & order matrices:condition 2 satisfied [sorted_beta_2, order_beta_2] = sort(beta_all(count2,1)); % sorting the beta matrix in ascending order [sorted_ah_2, order_ah_2] = sort(ah_all(count2,1)); [sorted_muv_2, order_muv_2] = sort(muv_all(count2,1)); [sorted_Ct0_2, order_Ct0_2] = sort(Ct0_all(count2,1)); [sorted_b1_2, order_b1_2] = sort(b1_all(count2,1)); [sorted_b2_2, order_b2_2] = sort(b2_all(count2,1)); [sorted_x_2, order_x_2]=sort(x_rlz(count2,1));% y (R0) is rank transformed [sorted_y_2, order_y_2]=sort(y_rlz(count2,1));

for i=1:size(count2,2) beta_rank_2(order_beta_2(i))=i; i ah_rank_2(order_ah_2(i))=i; muv_rank_2(order_muv_2(i))=i; Ct0_rank_2(order_Ct0_2(i))=i; b1_rank_2(order_b1_2(i))=i; b2_rank_2(order_b2_2(i))=i; x_rank_2(order_x_2(i))=i; y_rank_2(order_y_2(i))=i; end % End code: that uses count 2 to sort & order matrices:condition 2 satisfied

% ################################################################# % Below start: i/p param beta & o/p param y2:condition 2 (use count2)satisfied X_y_beta_2 = [ones(size(beta_rank_2')) b2_rank_2' muv_rank_2' Ct0_rank_2' ah_rank_2' b1_rank_2'];

[b_y_beta_2, bint_y_beta_2, r_y_beta_2]= regress(y_rank_2', X_y_beta_2); % Removes NaN data [b_beta_y_2, bint_beta_y_2, r_beta_y_2]= regress(beta_rank_2', X_y_beta_2);

[prcc_y_beta_2, p_val_y_beta_2]= corrcoef(r_beta_y_2, r_y_beta_2);

prcc(5)=prcc_y_beta_2(1,2); p_val(5)=p_val_y_beta_2(1,2); % ################################################################# % Above end: i/p param beta & o/p param y2:condition 2 (use count2)satisfied

% ################################################################# % Below start: i/p param ah & o/p param y2: condition 2 (use count2 OR underscore 2)satisfied X_y_ah_2 = [ones(size(beta_rank_2')) beta_rank_2' muv_rank_2' Ct0_rank_2' b1_rank_2' b2_rank_2'];

[b_y_ah_2, bint_y_ah_2, r_y_ah_2]= regress(y_rank_2', X_y_ah_2); % Removes NaN data [b_ah_y_2, bint_ah_y_2, r_ah_y_2]= regress(ah_rank_2', X_y_ah_2);

[prcc_y_ah_2, p_val_y_ah_2]= corrcoef(r_ah_y_2, r_y_ah_2);

prcc(6)=prcc_y_ah_2(1,2); p_val(6)=p_val_y_ah_2(1,2); % ################################################################# % Above end: i/p param ah & o/p param y2:condition 2 (use count2)satisfied

% ################################################################# % Below start: i/p param muv & o/p param y2: condition 2 (use count2)satisfied X_y_muv_2 = [ones(size(beta_rank_2')) beta_rank_2' ah_rank_2' Ct0_rank_2' b1_rank_2' b2_rank_2'];

[b_y_muv_2, bint_y_muv_2, r_y_muv_2]= regress(y_rank_2', X_y_muv_2); % Removes NaN data [b_muv_y_2, bint_muv_y_2, r_muv_y_2]= regress(muv_rank_2', X_y_muv_2);

[prcc_y_muv_2, p_val_y_muv_2]= corrcoef(r_muv_y_2, r_y_muv_2);

prcc(7)=prcc_y_muv_2(1,2); p_val(7)=p_val_y_muv_2(1,2); % ################################################################# % Above end: i/p param muv & o/p param y2: condition 2 (use count2)satisfied

% ################################################################# % Below start: i/p param Ct0 & o/p param y2:condition 2 (use count2)satisfied X_y_Ct0_2 = [ones(size(muv_rank_2')) beta_rank_2' ah_rank_2' muv_rank_2' b1_rank_2' b2_rank_2'];

[b_y_Ct0_2, bint_y_Ct0_2, r_y_Ct0_2]= regress(Ct0_rank_2', X_y_Ct0_2); [b_Ct0_y_2, bint_Ct0_y_2, r_Ct0_y_2]= regress(Ct0_rank_2', X_y_Ct0_2);

[prcc_y_Ct0_2, p_val_y_Ct0_2]= corrcoef(r_Ct0_y_2, r_y_Ct0_2);

prcc(8)=prcc_y_Ct0_2(1,2); p_val(8)=p_val_y_Ct0_2(1,2); % ################################################################# % Above end: i/p param Ct0 & o/p param x1:condition 1 (use count1)satisfied

% ################################################################# % Below start: i/p param b1 & o/p param y2:condition 2 (use count2)satisfied

X_y_b1_2 = [ones(size(beta_rank_2')) beta_rank_2' muv_rank_2' Ct0_rank_2' ah_rank_2' b2_rank_2'];

[b_y_b1_2, bint_y_b1_2, r_y_b1_2]= regress(y_rank_2', X_y_b1_2); % Removes NaN data [b_b1_y_2, bint_b1_y_2, r_b1_y_2]= regress(b1_rank_2', X_y_b1_2);

[prcc_y_b1_2, p_val_y_b1_2]= corrcoef(r_b1_y_2, r_y_b1_2);

prcc(9)=prcc_y_b1_2(1,2); p_val(9)=p_val_y_b1_2(1,2);

% ################################################################# % Above end: i/p param b1 & o/p param y2:condition 2(count2 OR underscore 2)satisfied

% ################################################################# % Below start: i/p param b2 & o/p param y2:condition 2 (use count2)satisfied X_y_b2_2 = [ones(size(beta_rank_2')) beta_rank_2' muv_rank_2' Ct0_rank_2' ah_rank_2' b1_rank_2'];

[b_y_b2_2, bint_y_b2_2, r_y_b2_2]= regress(y_rank_2', X_y_b2_2); % Removes NaN data [b_b2_y_2, bint_b2_y_2, r_b2_y_2]= regress(b2_rank_2', X_y_b2_2);

[prcc_y_b2_2, p_val_y_b2_2]= corrcoef(r_b2_y_2, r_y_b2_2);

prcc(10)=prcc_y_b2_2(1,2); p_val(10)=p_val_y_b2_2(1,2);

% ################################################################# % Above end: i/p param b2 & o/p param y2:condition 2(count2 OR underscore 2)satisfied

% Start code using count3 to sort & order matrices, extracting the rows for which condition 3 was satisfied [sorted_beta_3, order_beta_3] = sort(beta_all(count3,1)); % sorting the beta matrix in ascending order [sorted_ah_3, order_ah_3] = sort(ah_all(count3,1)); [sorted_muv_3, order_muv_3] = sort(muv_all(count3,1)); [sorted_Ct0_3, order_Ct0_3] = sort(Ct0_all(count3,1)); [sorted_b1_3, order_b1_3] = sort(b1_all(count3,1)); [sorted_b2_3, order_b2_3] = sort(b2_all(count3,1)); [sorted_x_3, order_x_3]=sort(x_rlz(count3,1));% y (R0) is rank transformed [sorted_y_3, order_y_3]=sort(y_rlz(count3,1));

size(count3,2) for i=1:size(count3,2) beta_rank_3(order_beta_3(i))=i; ah_rank_3(order_ah_3(i))=i; muv_rank_3(order_muv_3(i))=i; Ct0_rank_3(order_Ct0_3(i))=i; b1_rank_3(order_b1_3(i))=i; b2_rank_3(order_b2_3(i))=i; x_rank_3(order_x_3(i))=i; y_rank_3(order_y_3(i))=i; end % End code:using count 3 to sort & order matrices:for condition 3 % #################################################################

% ################################################################# % Below start: i/p param beta & o/p param y3:condition 3(count3)satisfied X_y_beta_3 = [ones(size(beta_rank_3')) b2_rank_3' muv_rank_3' Ct0_rank_3' ah_rank_3' b1_rank_3'];

[b_y_beta_3, bint_y_beta_3, r_y_beta_3]= regress(y_rank_3', X_y_beta_3); % Removes NaN data [b_beta_y_3, bint_beta_y_3, r_beta_y_3]= regress(beta_rank_3', X_y_beta_3);

[prcc_y_beta_3, p_val_y_beta_3]= corrcoef(r_beta_y_3, r_y_beta_3);

prcc(11)=prcc_y_beta_3(1,2); p_val(11)=p_val_y_beta_3(1,2); % ################################################################# % Above end: i/p param beta & o/p param y2:condition 3 (use count3)satisfied

% ################################################################# % Below start: i/p param ah & o/p param y3:condition 3(count3 OR underscore 3)satisfied X_y_ah_3 = [ones(size(beta_rank_3')) beta_rank_3' muv_rank_3' Ct0_rank_3' b1_rank_3' b2_rank_3'];

[b_y_ah_3, bint_y_ah_3, r_y_ah_3]= regress(y_rank_3', X_y_ah_3); % Removes NaN data [b_ah_y_3, bint_ah_y_3, r_ah_y_3]= regress(ah_rank_3', X_y_ah_3);

[prcc_y_ah_3, p_val_y_ah_3]= corrcoef(r_ah_y_3, r_y_ah_3);

prcc(12)=prcc_y_ah_3(1,2); p_val(12)=p_val_y_ah_3(1,2); % ################################################################# % Above end: i/p param ah & o/p param y3:condition 3(use count3)satisfied

% ################################################################# % Below start: i/p param muv & o/p param y3:condition 3 (count3)satisfied X_y_muv_3 = [ones(size(beta_rank_3')) beta_rank_3' ah_rank_3' Ct0_rank_3' b1_rank_3' b2_rank_3'];

[b_y_muv_3, bint_y_muv_3, r_y_muv_3]= regress(y_rank_3', X_y_muv_3); % Removes NaN data [b_muv_y_3, bint_muv_y_3, r_muv_y_3]= regress(muv_rank_3', X_y_muv_3);

[prcc_y_muv_3, p_val_y_muv_3]= corrcoef(r_muv_y_3, r_y_muv_3);

prcc(13)=prcc_y_muv_3(1,2); p_val(13)=p_val_y_muv_3(1,2); % ################################################################# % Above end: i/p param muv & o/p param y3:condition 3 (count3)satisfied

% ################################################################# % Below start: i/p param Ct0 & o/p param y3:condition 3(count3)satisfied X_y_Ct0_3 = [ones(size(muv_rank_3')) beta_rank_3' ah_rank_3' muv_rank_3' b1_rank_3' b2_rank_3'];

[b_y_Ct0_3, bint_y_Ct0_3, r_y_Ct0_3]= regress(Ct0_rank_3', X_y_Ct0_3); [b_Ct0_y_3, bint_Ct0_y_3, r_Ct0_y_3]= regress(Ct0_rank_3', X_y_Ct0_3);

[prcc_y_Ct0_3, p_val_y_Ct0_3]= corrcoef(r_Ct0_y_3, r_y_Ct0_3);

prcc(14)=prcc_y_Ct0_3(1,2); p_val(14)=p_val_y_Ct0_3(1,2); % ################################################################# % Above end: i/p param Ct0 & o/p param y2:condition 3(count3)satisfied

% ################################################################# % Below start: i/p param b1 & o/p param y3:condition 3(use count3)satisfied X_y_b1_3 = [ones(size(beta_rank_3')) beta_rank_3' muv_rank_3' Ct0_rank_3' ah_rank_3' b2_rank_3'];

[b_y_b1_3, bint_y_b1_3, r_y_b1_3]= regress(y_rank_3', X_y_b1_3); % Removes NaN data [b_b1_y_3, bint_b1_y_3, r_b1_y_3]= regress(b1_rank_3', X_y_b1_3);

[prcc_y_b1_3, p_val_y_b1_3]= corrcoef(r_b1_y_3, r_y_b1_3);

prcc(15)=prcc_y_b1_3(1,2); p_val(15)=p_val_y_b1_3(1,2); % ################################################################# % Above end: i/p param b1 & o/p param y3:condition 3(count3 OR underscore 3)satisfied

% Start of code that uses count4 as this sort and order matrices are used % for condition 4, extracting the rows for which condition 4 was satisfied [sorted_beta_4, order_beta_4] = sort(beta_all(count4,1)); % sorting the beta matrix in ascending order [sorted_ah_4, order_ah_4] = sort(ah_all(count4,1)); [sorted_muv_4, order_muv_4] = sort(muv_all(count4,1)); [sorted_Ct0_4, order_Ct0_4] = sort(Ct0_all(count4,1)); [sorted_b1_4, order_b1_4] = sort(b1_all(count4,1)); [sorted_b2_4, order_b2_4] = sort(b2_all(count4,1)); [sorted_x_4, order_x_4]=sort(x_rlz(count4,1));% y (R0) is rank transformed [sorted_y_4, order_y_4]=sort(y_rlz(count4,1));

size(count4,2) for i=1:size(count4,2) beta_rank_4(order_beta_4(i))=i; ah_rank_4(order_ah_4(i))=i; muv_rank_4(order_muv_4(i))=i; Ct0_rank_4(order_Ct0_4(i))=i; b1_rank_4(order_b1_4(i))=i; b2_rank_4(order_b2_4(i))=i; x_rank_4(order_x_4(i))=i; y_rank_4(order_y_4(i))=i; end % End of code that uses count 4 as this sort & order matrices for condition 4 % #################################################################

% ################################################################# % Below start:i/p param ah & o/p param x4:condition 4(count4/underscore 4)satisfied X_x_ah_4 = [ones(size(beta_rank_4')) beta_rank_4' muv_rank_4' Ct0_rank_4' b1_rank_4' b2_rank_4'];

[b_x_ah_4, bint_x_ah_4, r_x_ah_4]= regress(x_rank_4', X_x_ah_4); % Removes NaN data [b_ah_x_4, bint_ah_x_4, r_ah_x_4]= regress(ah_rank_4', X_x_ah_4);

[prcc_x_ah_4, p_val_x_ah_4]= corrcoef(r_ah_x_4, r_x_ah_4);

prcc(16)=prcc_x_ah_4(1,2); p_val(16)=p_val_x_ah_4(1,2); % ################################################################# % Above end: i/p param ah & o/p param x4:condition 4 (use count4)satisfied

% ################################################################# % Below start: i/p param muv & o/p param x4:condition 4 (use count4)satisfied X_x_muv_4 = [ones(size(beta_rank_4')) beta_rank_4' ah_rank_4' Ct0_rank_4' b1_rank_4' b2_rank_4'];

[b_x_muv_4, bint_x_muv_4, r_x_muv_4]= regress(x_rank_4', X_x_muv_4); % Removes NaN data [b_muv_x_4, bint_muv_x_4, r_muv_x_4]= regress(muv_rank_4', X_x_muv_4);

[prcc_x_muv_4, p_val_x_muv_4]= corrcoef(r_muv_x_4, r_x_muv_4);

prcc(17)=prcc_x_muv_4(1,2); p_val(17)=p_val_x_muv_4(1,2); % ################################################################# % Above ends for: i/p parameter muv & o/p parameter x4 when condition 4 (use count4) is % satisfied

% ################################################################# % Below start: i/p param Ct0 & o/p param x4:condition 4 (use count4)satisfied X_x_Ct0_4 = [ones(size(muv_rank_4')) beta_rank_4' ah_rank_4' muv_rank_4' b1_rank_4' b2_rank_4'];

[b_x_Ct0_4, bint_x_Ct0_4, r_x_Ct0_4]= regress(x_rank_4', X_x_Ct0_4); % Removes NaN data [b_Ct0_x_4, bint_Ct0_x_4, r_Ct0_x_4]= regress(Ct0_rank_4', X_x_Ct0_4);

[prcc_x_Ct0_4, p_val_x_Ct0_4]= corrcoef(r_Ct0_x_4, r_x_Ct0_4);

prcc(18)=prcc_x_Ct0_4(1,2); p_val(18)=p_val_x_Ct0_4(1,2); % ################################################################# % Above end: i/p param Ct0 & o/p param x4:condition 4(use count4)satisfied

% ################################################################# % Below start: i/p param b1 & o/p param x4:condition 4(use count4)satisfied X_x_b1_4 = [ones(size(beta_rank_4')) beta_rank_4' muv_rank_4' Ct0_rank_4' ah_rank_4' b2_rank_4'];

[b_x_b1_4, bint_x_b1_4, r_x_b1_4]= regress(x_rank_4', X_x_b1_4); % Removes NaN data [b_b1_x_4, bint_b1_x_4, r_b1_x_4]= regress(b1_rank_4', X_x_b1_4);

[prcc_x_b1_4, p_val_x_b1_4]= corrcoef(r_b1_x_4, r_x_b1_4);

prcc(19)=prcc_x_b1_4(1,2); p_val(19)=p_val_x_b1_4(1,2); % ################################################################# % Above end: i/p param b1 & o/p param x4:condition 4(count4/underscore 4)satisfied

% ################################################################# % Below start: i/p param b2 & o/p param x4:condition 4 (use count4)satisfied X_x_b2_4 = [ones(size(beta_rank_4')) beta_rank_4' muv_rank_4' Ct0_rank_4' ah_rank_4' b1_rank_4'];

[b_x_b2_4, bint_x_b2_4, r_x_b2_4]= regress(x_rank_4', X_x_b2_4); % Removes NaN data [b_b2_x_4, bint_b2_x_4, r_b2_x_4]= regress(b1_rank_4', X_x_b2_4);

[prcc_x_b2_4, p_val_x_b2_4]= corrcoef(r_b2_x_4, r_x_b2_4);

prcc(20)=prcc_x_b2_4(1,2); p_val(20)=p_val_x_b2_4(1,2); % ################################################################# % Above end: i/p param b2 & o/p param x4:condition 4(count4/underscore 4)satisfied

end % Sensitivity Analysis ends here........ % #################################################################

% ################################################################# % Uncertainity Analysis starts here........

if analysis==1

   % mean,median,std,iqr
     mean_r0i_x= mean(x_rlz)';
     median_r0i_x=median(x_rlz)';
     mean_r0i_y= mean(y_rlz)';
     median_r0i_y=median(y_rlz)';

     std_r0i_x=std(x_rlz)';
     iqr_r0i_x=iqr(x_rlz)';
     std_r0i_y=std(y_rlz)';
     iqr_r0i_y=iqr(y_rlz)';

for k=1:length(x_rlz(1,:)) % can write as k=1:reali ??? t1=x_rlz(:,k); % pg1r0i(k)=length(t1(find(t1<=2)))/length(t1); % INTRODUCED UNDERSCORE X AS I HAVE ONE LOOP BELOW FOR Y AS WELL pg1r0i_x(k)=length(t1(find(t1<=2.5)))/length(t1); % Use this when mean(gamma)=0.23 end % You can use a logical expression to define X. For example, find(X > 2) % returns linear indices corresponding to the entries of X that are greater than 2. pg1r0i_x=pg1r0i_x'; pg1minusr0_x=1-pg1r0i_x;

 %%mean over all realizations

for k=1:length(y_rlz(1,:)) % can write as k=1:reali ??? t1=y_rlz(:,k); % pg1r0i(k)=length(t1(find(t1<=2)))/length(t1); % INTRODUCED UNDERSCORE Y BELOW pg1r0i_y(k)=length(t1(find(t1<=2.5)))/length(t1); % Use this when mean(gamma)=0.23 end % You can use a logical expression to define X. For example, find(X > 2) % returns linear indices corresponding to the entries of X that are greater than 2. pg1r0i_y=pg1r0i_y'; pg1minusr0_y=1-pg1r0i_y;

 %%mean over all realizations

   mean_alli_x=mean([mean_r0i_x median_r0i_x iqr_r0i_x std_r0i_x  pg1r0i_x pg1minusr0_x]);
   mean_alli_y=mean([mean_r0i_y median_r0i_y iqr_r0i_y std_r0i_y  pg1r0i_y pg1minusr0_y]);
%
%standard error
%
     se1_x=std([mean_r0i_x median_r0i_x iqr_r0i_x std_r0i_x pg1r0i_x pg1minusr0_x]);     
     ser1_x=se1_x/sqrt(length(mean_r0i_x));

     se1_y=std([mean_r0i_y median_r0i_y iqr_r0i_y std_r0i_y pg1r0i_y pg1minusr0_y]);     
     ser1_y=se1_y/sqrt(length(mean_r0i_y));

     %
%variation coefficient
%
     cv1_x=se1_x./mean_alli_x;
     cv1_y=se1_y./mean_alli_y;

end

% Uncertainity analysis ends here % #################################################################

%%%% File for storing data %%%%%%%%

% % % % % % % % Below needs to be changed after understanding what my UA is

% strcat - Concatenate strings horizontally, Syntax: combinedStr = strcat(s1, s2, ..., sN) diary(strcat(model,'_unif_uncrt_txt.txt'))

     disp('============================================================================')
     disp(strcat('model:  ',model,'   Analysis (Uncertainty=1 & Sensitivity=2) = ',sprintf('%2d',analysis)))
     disp(strcat('sample size = ',sprintf('%3.0e',dim),'    p =  ',sprintf('%3.2f', p), '    Days post spray =    ', sprintf('%3.2f', t)    ))
     disp('============================================================================')

%%%% File for storing data %%%%%%%%

% ################################################################# % Uncertainity results starts here

if analysis==1

%disp(' Realization Mean(R_d) Std(R_d) Pr(1<R_d<=2) Pr(2<R_d<=3)') disp(' Realization Mean(x) Std(x) Mean(y) Std(y) ') % Use this when mean(gamma)=0.23

for i=1:1 disp(strcat(sprintf('%6d',i),sprintf('%14.3f',mean_r0i_x(i)),sprintf(' %13.3f',std_r0i_x(i)),sprintf('%14.3f',mean_r0i_y(i)),sprintf(' %13.3f',std_r0i_y(i)) )); end disp('============================================================================')

% disp(strcat('mean',sprintf('%14.4f', mean_alli_x(1)),sprintf(' %13.4f', mean_alli_x(4)), ... % sprintf(' %13.4f', mean_alli_x(5)),sprintf(' %13.4f', mean_alli_x(6)),... % sprintf('%14.4f', mean_alli_y(1)),sprintf(' %13.4f', mean_alli_y(4)), ... % sprintf(' %13.4f', mean_alli_y(5)),sprintf(' %13.4f', mean_alli_y(6)) ));

disp('============================================================================') disp(strcat('se',sprintf(' %14.4f',ser1_x(1)),sprintf(' %14.4f',ser1_y(1)) ));

disp('============================================================================') disp(strcat('cv',sprintf(' %14.4f',cv1_x(1)),sprintf(' %14.4f',cv1_y(1)) ));

end

% Uncertainity results ends here % #################################################################

% ################################################################# % Senstivity results starts here

if analysis==2

disp('===================================================================') disp('') disp(' Parameter S(x1) p-value_x1 S(y2) p-value_y2 S(y3) p-value_y3 S(x4) p-value_x4') disp('-----------------------------------') disp(strcat( 'beta ',sprintf('\t\t%s ', '--'),sprintf(' %s ', ' --- '),sprintf(' \t\t%8.4f', prcc(5)),sprintf(' %13.3e', p_val(5)), sprintf('\t\t%8.4f', prcc(11)),sprintf('\t\t%13.3e', p_val(11)),sprintf('\t\t%s', '-----'),sprintf('\t\t%s', '----'))); disp(strcat( 'ah ',sprintf('\t\t%8.4f', prcc(1)),sprintf(' %13.3e', p_val(1)),sprintf('\t\t%8.4f', prcc(6)),sprintf(' %13.3e', p_val(6)), sprintf('\t\t%8.4f', prcc(12)),sprintf(' %13.3e', p_val(12)),sprintf('\t\t%8.4f', prcc(16)),sprintf(' %13.3e', p_val(16)))); disp(strcat( 'muv ',sprintf('\t\t%8.4f', prcc(2)),sprintf(' %13.3e', p_val(2)),sprintf('\t\t%8.4f', prcc(7)),sprintf(' %13.3e', p_val(7)), sprintf('\t\t%8.4f', prcc(13)),sprintf(' %13.3e', p_val(13)),sprintf('\t\t%8.4f', prcc(17)),sprintf(' %13.3e', p_val(17)) )); disp(strcat( 'Ct0 ',sprintf('\t\t%8.4f', prcc(3)),sprintf(' %13.3e', p_val(3)),sprintf('\t\t%8.4f', prcc(8)),sprintf(' %13.3e', p_val(8)), sprintf('\t\t%8.4f', prcc(14)),sprintf(' %13.3e', p_val(14)),sprintf('\t\t%8.4f', prcc(18)),sprintf(' %13.3e', p_val(18)) )); disp(strcat( 'b1 ',sprintf('\t\t%8.4f', prcc(4)),sprintf(' %13.3e', p_val(4)),sprintf('\t\t%8.4f', prcc(9)),sprintf(' %13.3e', p_val(9)), sprintf('\t\t%8.4f', prcc(15)),sprintf(' %13.3e', p_val(15)),sprintf('\t\t%8.4f', prcc(19)),sprintf(' %13.3e', p_val(19)) )); disp(strcat( 'b2 ',sprintf('\t\t%s ', ' -----'),sprintf(' %s ', ' ------ '),sprintf(' \t\t%8.4f', prcc(10)),sprintf(' %13.3e', p_val(10)), sprintf('\t\t%s ', '------'),sprintf(' %s ', '-------'),sprintf('\t\t%8.4f', prcc(20)),sprintf(' %13.3e', p_val(20)) ));

end

% Sensitivity results ends here % #################################################################

diary off

if analysis==2

 figure('Color',[1 1 1]);   
 figure(1)
 subplot(2,2,1);
    x1=[prcc(1); prcc(3); prcc(4)];
    barh(x1);
%   fname={'\beta','a_h','\mu_v','C_{t0}','b_1','b_2'}; set(gca,'yticklabel',{fname}); 
%   InputParameters = ['\beta';'a_h';'\mu_v';'C_{t0}';'b_1';'b_2']; set(gca,'YTickLabels',InputParameters)
    set(gca,'yticklabel',{'a_h','C_{t0}','b_1'})
    title('First feasible solution Sensitivity: prcc values for x_1, y_1 is fixed ')

 subplot(2,2,2);
    y2=[prcc(6); prcc(8); prcc(9); prcc(10)];
    barh(y2);
    set(gca,'yticklabel',{'a_h','C_{t0}','b_1','b_2'})
    title('Second feasible solution Sensitivity: x_2 is fixed, prcc values for y_2 ')

 subplot(2,2,3);
    y2=[prcc(12); prcc(14); prcc(15)];
    barh(y2);
    set(gca,'yticklabel',{'a_h','C_{t0}','b_1'})
    title('Third feasible solution Sensitivity: x_3 is fixed, prcc values for y_3 ')

 subplot(2,2,4);
    y2=[prcc(16); prcc(18); prcc(19); prcc(20)];
    barh(y2);
    set(gca,'yticklabel',{'a_h','C_{t0}','b_1','b_2'})
    title('Fourth feasible solution Sensitivity: prcc values for x_4, y_4 is fixed')
end

if analysis==1

figure('Color',[1 1 1]); figure(1)

% subplot(3,3,1); subplot(2,2,1); hist(ah,100) title('Sampling Distribution of a_h, mean=','fontsize',22);

% subplot(3,3,2); % hist(muv,100) % title('Distribution of \mu_v','fontsize',22);

subplot(2,2,2);
  hist(Ct0,100)
  title('Distribution of C_{t0}','fontsize',22);

% subplot(3,3,4); subplot(2,2,3); hist(b1,100) title('Distribution of b_1','fontsize',22);

    subplot(2,2,4);
    hist(b2,100)
    title('Distribution of b_2','fontsize',22);

% % % % % subplot(3,3,6); % % % % % % [xcounts1,~] = hist(x_rlz1,100);[xcounts2,~] = hist(x_rlz2,100);[xcounts3,~] = hist(x_rlz3,100);[xcounts4,~] = hist(x_rlz4,100); % % % % % % histmat_x_rlz = [reshape(xcounts1,100,1) reshape(xcounts2,100,1) reshape(xcounts3,100,1) reshape(xcounts4,100,1)]; % % % % % % bar(histmat_x_rlz, 1.2); % % % % % hist(x_rlz,100); % % % % % hist(x_rlz) % % % % % % bar(x_rlz, 1.2); % % % % % title('Distribution of x','fontsize',22) % % % % % % % % % % subplot(3,3,7); % % % % % % [ycounts1,~] = hist(y_rlz1,100);[ycounts2,~] = hist(y_rlz2,100);[ycounts3,~] = hist(y_rlz3,100);[ycounts4,~] = hist(y_rlz4,100); % % % % % % histmat_y_rlz = [reshape(ycounts1,100,1) reshape(ycounts2,100,1) reshape(ycounts3,100,1) reshape(ycounts4,100,1)]; % % % % % % bar(histmat_y_rlz, 1.2); % % % % % % bar(y_rlz, 1.2); % % % % % hist(y_rlz,100); % % % % % % hist(y_rlz(:,reali),100) % % % % % title('Distribution of y','fontsize',22)

    % beta is taken as a constant in manuscript 1
% % % % %     subplot(3,3,1);
% % % % %   hist(beta,100)
% % % % %   title('Distribution of \beta','fontsize',22);

figure('Color',[1 1 1]);  
figure(2)
  subplot(2,2,1); 
  hist(x_rlz(count1,1))
  title('Distribution of x when first pair of conditions is satisfied','fontsize',14)

    subplot(2,2,2); 
    hist(x_rlz(count2,1))
    title('Distribution of x when second pair of conditions is satisfied','fontsize',14)

    subplot(2,2,3); 
    hist(x_rlz(count3,1))
    title('Distribution of x when third pair of conditions is satisfied','fontsize',14)

    subplot(2,2,4); 
    hist(x_rlz(count4,1))
    title('Distribution of x when fourth pair of conditions is satisfied','fontsize',14)

    figure('Color',[1 1 1]);
    figure(3)
    subplot(2,2,1); 
    hist(y_rlz(count1,1))
    title('Distribution of y when first pair of conditions is satisfied','fontsize',14)

    subplot(2,2,2); 
    hist(y_rlz(count2,1))
    title('Distribution of y when second pair of conditions is satisfied','fontsize',14)

    subplot(2,2,3); 
    hist(y_rlz(count3,1))
    title('Distribution of y when third pair of conditions is satisfied','fontsize',14)

    subplot(2,2,4); 
    hist(y_rlz(count4,1))
    title('Distribution of y when fourth pair of conditions is satisfied','fontsize',14)

end

fprintf(['Number of times condition 1,2,3 & 4 are satisfied are ',num2str(size(count1,2)),', ',num2str(size(count2,2)),', ',num2str(size(count3,2)),' & ',num2str(size(count4,2)),' times respectively','\n']); fprintf(['Number of times infeasible condition arises is ',num2str(dim-size(count1,2)-size(count2,2)-size(count3,2)-size(count4,2)),' times','\n']);

fprintf(['Means of x realizations when condition 1,2,3 & 4 are satisfied are ',num2str(mean(x_rlz(count1,1))),', ',num2str(mean(x_rlz(count2,1))),', ',num2str(mean(x_rlz(count3,1))),' & ',num2str(mean(x_rlz(count4,1))),' times respectively','\n']); fprintf(['Means of y realizations when condition 1,2,3 & 4 are satisfied are ',num2str(mean(y_rlz(count1,1))),', ',num2str(mean(y_rlz(count2,1))),', ',num2str(mean(y_rlz(count3,1))),' & ',num2str(mean(y_rlz(count4,1))),' times respectively','\n']);

toc %%%%%%%%%%%%%%%%%%