clc
clear all
close all
for x=1:1000:50000
syms Rab Rca Rbc
a = 1;
b = 5.45*(10^4);
y_gam11(x) = gamcdf(x,a,b);
y_gam1(x)=1-y_gam11(x);
a = 0.01*7.06;
b = 1.2*5.45*(10^4);
y_gam22(x) = gamcdf(x,a,b);
y_gam2(x)=1-y_gam22(x);
a = 0.075*7.06;
b = 1.5*5.45*(10^4);
y_gam33(x) = gamcdf(x,a,b);
y_gam3(x)=1-y_gam33(x);
gammaab(x)=y_gam1(x).*y_gam2(x);
gammabc(x)=y_gam2(x).*y_gam3(x);
gammaca(x)=y_gam3(x).*y_gam1(x);
eq1(x)=( 1-( (1-Rab)*(1-(Rca*Rbc)) ) )==gammaab(x);
eq2(x)=( 1-( (1-Rbc)*(1-(Rab*Rca)) ) )==gammabc(x);
eq3(x)=( 1-( (1-Rca)*(1-(Rab*Rbc)) ) )==gammaca(x);
sol=solve([eq1(x),eq2(x),eq3(x)], [Rab, Rbc, Rca]);
if Rab>0 && Rab<1
RabSol{x} = sol.Rab;
end
if Rbc>0 && Rbc<1
RbcSol{x} = sol.Rbc;
end
if Rca>0 && Rca<1
RcaSol{x} = sol.Rca;
end
end
plot(x,RabSol)
plot(x,RbcSol)
plot(x,RcaSol)