Symbols and bar thickness in BAR graphs

8 views (last 30 days)
Anant
Anant on 29 Oct 2012
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 %%%%%%%%%%%%%%%%%%

Answers (0)

Categories

Find more on Nonlinear Mixed-Effects Modeling in Help Center and File Exchange

Products

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!