% Author: Khaled Nakhleh SEED = 32 rng(SEED) clc, clear all %%%%%%%%%%%%%% CONSTANTS %%%%%%%%%%%%%%%% P = 0.2; % trans prob from good to bad Q = 0.8; % trans prob from bad to good timesteps = 1000000; kRange = 100; RUNS = 1000; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% for run = 1:RUNS fprintf("current run: %i", run) aoi = zeros(1, timesteps); deliveries = zeros(1, timesteps); B = zeros(1, timesteps); mc = createMC(P, Q); walkIns = simulate(mc, timesteps); walkIns(walkIns == 2) = 0; % converting bad channel from 2 to 0 for x = 1:length(walkIns) if x == 1 aoi(x) = 1; if walkIns(x) == 1 deliveries(x) = x; end else if walkIns(x-1) == 1 aoi(x) = 1; deliveries(x-1) = x-1; else aoi(x) = aoi(x-1) + 1; end end end walkIns(end) = []; % removing timestep+1 extra step aoi(end) = []; walkIns'; aoi; deliveries; %%%%%%%%%%%%%%%%% EMIPIRICAL AOI %%%%%%%%%%%%%%%%%%%%%%%%%%% aoi; empiricalAOI = sum(aoi)/timesteps %%%%%%%%%%%%%%%%%%%%%%%% AOI %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% walkIns'; An = nonzeros(deliveries)'; B(1) = An(1) - 0; for i = 1:(length(An)-1) B(i+1) = An(i+1) - An(i); end B = nonzeros(B)'; BSquared = sum(B.^2)/ length(B); twoB = 2 * sum(B) / length(B); avgAOI = BSquared / twoB + 0.5 %%%%%%%%%%%%%%%%%% ESTIMATED AOI %%%%%%%%%%%%%%%%%%%%%%%%%% POverQ = P / (P + Q); ms = 1 - POverQ; ms; GVal = 0.0; for k = 1:(kRange) GVal = GVal + (Q/(P+Q))*(1 - P - Q)^(k+1-1); end GVal; secondVal = POverQ - (POverQ)^2; vs = GVal*POverQ*2 + secondVal estimatedAOIFromMsVs = (0.5*((vs/ms^2) + (1/ms) - 1)) + 1 %%%%%%%%%%%%%% MU and VAR calculations %%%%%%%%%%%%%%%%%%% activations(run) = sum(walkIns); % adding the results to save array empiricalAOIs(run) = empiricalAOI; avgAOIVals(run) = avgAOI; estimatedAOIFromMsVsVals(run) = estimatedAOIFromMsVs; tableActivations = table(walkIns); end averageMu = sum(activations) / RUNS for i = 1:(RUNS) vars(i) = (activations(i) - averageMu)^2 / timesteps; end estimatedVar = sum(vars) / (RUNS - 1) muVal = averageMu / timesteps estimatedAOIFromRuns = 0.5*((estimatedVar / muVal^2) + (1/muVal) - 1) + 1 empiricalAOIs = empiricalAOIs'; avgAOIVals = avgAOIVals'; estimatedAOIFromMsVsVals = estimatedAOIFromMsVsVals'; vars = vars'; activations = activations'; estimatedAOIFromRuns = repelem(estimatedAOIFromRuns, RUNS)'; mu = repelem(averageMu, RUNS)'; results = table(empiricalAOIs, avgAOIVals, estimatedAOIFromMsVsVals, estimatedAOIFromRuns, vars, mu, activations); Filename = sprintf('single_client_results/results_p_%.2f_q_%.2f_runs_%d_timesteps_%d.csv', P, Q, RUNS, timesteps); writetable(results,Filename) function mc = createMC(P, G) TRANS = [1-P, P ; G, 1-G]; % 1 - p is trans from state G to state G. mc = dtmc(TRANS, 'StateNames', ["Good", "Bad"]); end