# n Accuracy - Code

# To measure the accuracy of the estimation proceedure as I vary "n".


#set argument values
alpha <- 2.25
beta <- 3
lambda <- 1

#set n values
nvalues<-c(10,25,50,75,100,125,150,175,200,225,250,300,350,400,450,500,550,600,650,700,750,800,850,900
           ,950,1000,1100,1200,1300,1400,1500,1600,1700,1800,1900,2000,2500)

# clear the error vector
error=c()

# run for loop
for(i in nvalues){
  n<-i
  newerror<-errors(n, alpha, beta, lambda) # this "errors" function is in the functions tab
  error=c(error,newerror)
}
error # run this line, then copy and paste the error vector below

error=abs(error) 
error=error*100
plot(nvalues,error, type="l", ylim=c(0,150), main="Estimate Error vs Number of Arrivals (Varying Parameters)", 
     ylab="Percentage Error", xlab="Number of Arrivals")









############################################################################################
### Here are the values found so that I dont need to run that again  (removing set.seed) ###
############################################################################################
# these errors have been divided by 3

#set n values
nvalues<-c(10,25,50,75,100,125,150,175,200,225,250,300,350,400,450,500,550,600,650,700,750,800,850,900
           ,950,1000,1100,1200,1300,1400,1500,1600,1700,1800,1900,2000,2500)

# errors found for n up to 2500 (alpha <- 2.25, beta <- 3, lambda <- 1)
errors1=c(1.68698232, 0.62856779, 1.32470589, 0.79422817, 0.47623088, 0.29251429, 0.44912733, 0.25767379, 0.09273400, 0.24485118,
          0.59353323, 0.17411433, 0.07988639, 0.12799148, 0.14915892, 0.30260157, 0.25954783, 0.08576282, 0.09311230, 0.08635029,
          0.09990123, 0.09546072, 0.10540401, 0.09576042, 0.07610974, 0.20534800, 0.13639272, 0.01663841, 0.06330755, 0.10013126,
          0.03386863, 0.03341786, 0.09036264, 0.09619811, 0.03834620, 0.08999917, 0.04373520)

# errors found for n up to 2500 (alpha <- 2.25, beta <- 3, lambda <- 1)
errors2=c(0.28700873, 0.31931982, 1.11168007, 0.55880354, 0.29107740, 0.39735057, 0.18916372, 0.31269181, 0.19350461, 0.22768275,
          0.15001536, 0.04600497, 0.06242419, 0.15346328, 0.14244700, 0.19974377, 0.06852106, 0.13774865, 0.07009237, 0.08544778,
          0.05653525, 0.02516168, 0.02903747, 0.09465042, 0.16136880, 0.06426391, 0.12490454, 0.03857611, 0.02041611, 0.12846649,
          0.04750295, 0.05442422, 0.02822789, 0.06269866, 0.03659746, 0.07481723, 0.05119580)

# errors found for n up to 2500 (alpha <- 2.25, beta <- 3, lambda <- 1)
errors3=c(1.13548919, 0.55295011, 0.16104273, 0.39699594, 0.30116435, 0.15453258, 0.34843850, 0.15353964, 0.11028405, 0.06959256,
          0.03156990, 0.02566006, 0.10442374, 0.22269423, 0.08942004, 0.08687360, 0.08062569, 0.04522704, 0.08873554, 0.19971647,
          0.07088254, 0.11217833, 0.04109870, 0.04203557, 0.06085740, 0.02774021, 0.05433981, 0.05600135, 0.15772356, 0.07913931,
          0.03940372, 0.04349799, 0.03432131, 0.09613920, 0.10691901, 0.08289716, 0.05410117)

# errors found for n up to 2500 (alpha <- 2.25, beta <- 3, lambda <- 1)
errors4=c(0.36711109, 0.10538777, 0.17510842, 0.81340404, 0.54424788, 0.24677990, 0.13023401, 0.76422307, 0.21239338, 0.12439860,
          0.33159212, 0.17638941, 0.14528097, 0.14167297, 0.07602200, 0.02695593, 0.11227598, 0.12515564, 0.14764396, 0.03653284,
          0.01549124, 0.11586081, 0.16132640, 0.03708148, 0.07289278, 0.09959860, 0.02310074, 0.07331059, 0.06910577, 0.05023177,
          0.06635710, 0.07408741, 0.17978943, 0.08985802, 0.08787736, 0.09176981, 0.04349867)

# errors found for n up to 2500 (alpha <- 2.25, beta <- 3, lambda <- 1)
errors5=c(1.47348683, 3.53335467, 0.34194399, 0.20683933, 0.93896112, 0.13014722, 0.16311658, 0.15031503, 0.12805141, 0.19842525,
          0.29592702, 0.16968282, 0.14552593, 0.09636042, 0.18481554, 0.08174549, 0.12946416, 0.08062132, 0.19514862, 0.11470684,
          0.12400592, 0.14568271, 0.21715374, 0.11944951, 0.10396573, 0.05406316, 0.12634583, 0.11938980, 0.08836416, 0.10588617,
          0.04425322, 0.08287835, 0.03085518, 0.04588232, 0.06265193, 0.09968219, 0.02240977)

# errors found for n up to 2500 (alpha <- 2.25, beta <- 3, lambda <- 1)
errors6=c(0.82007816, 0.11530079, 0.29256931, 0.31567964, 0.02707893, 0.18118903, 0.30006903, 0.14992453, 0.20065110, 0.13590271,
          0.11349497, 0.32877474, 0.08568590, 0.28015696, 0.15233986, 0.11992335, 0.07058811, 0.05230925, 0.13600071, 0.05399470,
          0.15942050, 0.08560285, 0.09854435, 0.04279388, 0.12798431, 0.04333639, 0.09418721, 0.08806445, 0.01401458, 0.06716407,
          0.07938661, 0.03434735, 0.04649978, 0.06241612, 0.02526109, 0.08909531, 0.01984788)

# errors found for n up to 2500 (alpha <- 2.25, beta <- 3, lambda <- 1)
errors7=c(2.23515513, 0.55065039, 1.72541667, 0.23093006, 0.35955511, 0.25564816, 0.03867841, 0.16030675, 0.23296199, 0.19039690,
          0.12846701, 0.27805078, 0.16169784, 0.19796703, 0.07911744, 0.15919477, 0.30559940, 0.16401727, 0.02245569, 0.06520802,
          0.14381769, 0.25535772, 0.14900726, 0.17896214, 0.06791195, 0.07510559, 0.14897200, 0.10077359, 0.13023300, 0.07700237,
          0.12876516, 0.11946685, 0.07331343, 0.04245584, 0.05399121, 0.01878355, 0.10271367)

# errors found for n up to 2500 (alpha <- 2.25, beta <- 3, lambda <- 1)
errors8=c(1.54283000, 1.74987291, 0.44735487, 0.44681199, 0.35197040, 0.32960340, 0.28298543, 0.23025654, 0.04805144, 0.14689511,
          0.29817065, 0.11190638, 0.12717153, 0.04744580, 0.10239420, 0.13131867, 0.01724246, 0.01614890, 0.08930352, 0.07268733,
          0.07935311, 0.10753303, 0.06136514, 0.03178560, 0.02559235, 0.03276829, 0.07515130, 0.06145899, 0.02561290, 0.13199597,
          0.06530630, 0.01991813, 0.08613773, 0.03776365, 0.05820599, 0.06699449, 0.08284698)




# Plot all 
plot(nvalues, errors1*100, type="l", ylim=c(0,150), main="Estimate Error vs Number of Arrivals", 
     ylab="Percentage Error", xlab="Number of Arrivals", yaxt='n')
lines(nvalues, errors2*100, col="blue")
lines(nvalues, errors3*100, col="red")
lines(nvalues, errors4*100, col="green")
lines(nvalues, errors5*100, col="#CC33FF")
lines(nvalues, errors6*100, col="orange")
lines(nvalues, errors7*100, col="#0033FF")
lines(nvalues, errors8*100, col="#993300")
# to add my prefered y axis
g<-c(0,50,100,150)
axis(2, at = g, labels=c("0%","50%","100%","150%"), las=2) 










########### get the average of each error vector and plot that ###########
avgerror = c()
for (i in 1:length(nvalues)){
  avgerror[i] = (errors1[i]+errors2[i]+errors3[i]+errors4[i]+errors5[i]+errors6[i]+errors7[i]+errors8[i])/8
}
avgerror
plot(nvalues, avgerror*100, type="l", ylim=c(0,150), main="Average Estimate Error vs Number of Arrivals (Multiple Simulations)", 
     ylab="Percentage Error (Average for all parameters)", xlab="Number of Arrivals", yaxt='n', col="red")
g<-c(0,50,100,150)
axis(2, at = g, labels=c("0%","50%","100%","150%"), las=2) 

# plot on a log scale
plot(log(nvalues), avgerror*100, type="l", ylim=c(0,150), main="Average Estimate Error vs Number of Arrivals (log scaled)", 
     ylab="Percentage Error (Average for each parameter)", xlab="Number of Arrivals", yaxt='n', xaxt='n')
g<-c(0,50,100,150) # y axis
axis(2, at = g, labels=c("0%","50%","100%","150%"), las=2) 
axis(1, at = log(nvalues), labels=nvalues, las=2) 





######    vary each parameter one-by-one (and all together) by an     ##########
######    order of magnitude and see how that affects the error plot  ##########
######    plot it on the avg error plot
alpha <- 2.25
beta <- 30
lambda <- 1

# clear the error vector
error=c()
# run for loop
for(i in nvalues){
  n<-i
  newerror<-errors(n, alpha, beta, lambda) # this "errors" function is in the functions tab
  error=c(error,newerror)
}
error # paste this below

# varying UP by an order of magnitude

# errors found for n up to 2500 (alpha <- 22.5, beta <- 30, lambda <- 10) changing all up (alpha wont work up)
errors1=c(1.87080162, 1.51714335, 0.83306281, 0.77792315, 0.16234540, 0.06102634, 0.11283971, 0.22021508, 0.08863498, 0.27723599,
          0.09551576, 0.02583214, 0.24964262, 0.05177590, 0.17387340, 0.11010308, 0.16098389, 0.10534065, 0.18754867, 0.13843700,
          0.13331066, 0.08837246, 0.02727609, 0.07582609, 0.09732341, 0.07952910, 0.02462238, 0.02436249, 0.04332119, 0.08369200,
          0.09114491, 0.03557326, 0.02398421, 0.18723121, 0.02380288, 0.11439677, 0.04479409) 

# errors found for n up to 2500 (alpha <- 2.25, beta <- 30, lambda <- 1) changing beta
errors2=c(1.21636200, 0.83139254, 1.38799350, 0.65238838, 0.68018296, 0.10705205, 0.42805037, 0.09925392, 0.18136082, 0.55261941,
          0.39796663, 0.28391955, 0.11479223, 0.16832919, 0.21316417, 0.19712212, 0.27247244, 0.53899342, 0.12319267, 0.25781520,
          0.10727385, 0.14959995, 0.15583202, 0.05085551, 0.10091016, 0.07242381, 0.12807135, 0.14552988, 0.05465267, 0.18136294,
          0.10567308, 0.23239753, 0.03695703, 0.06272712, 0.18515092, 0.08034624, 0.02552389)

# errors found for n up to 2500 (alpha <- 2.25, beta <- 3, lambda <- 10) changing lambda
errors3=c(2.36399495, 0.84403892, 0.51302914, 4.33934347, 1.21033657, 0.13019513, 0.25538850, 0.69737945, 0.69392144, 0.14943617,
          1.02787903, 0.18285121, 0.33838880, 0.43913551, 0.51770468, 0.32434083, 0.27153834, 0.05683072, 0.20143372, 0.76123397,
          0.07001602, 0.16961454, 0.53575571, 0.49990709, 0.11678775, 0.17321222, 0.19133916, 0.31409615, 0.18272530, 0.33776692,
          0.37955037, 0.32332712, 0.08223859, 0.15283952, 0.02394233, 0.12190118, 0.13189039)

plot(nvalues, avgerror*100, type="l", ylim=c(0,300), main="Estimate Error vs Number of Arrivals (varying parameters)", 
     ylab="Percentage Error (Average for each parameter)", xlab="Number of Arrivals", yaxt='n')
g<-c(0,50,100,150,200,250,300)
axis(2, at = g, labels=c("0%","50%","100%","150%","200%","250%","300%"), las=2) 
lines(nvalues, errors1*100, col="blue")
lines(nvalues, errors2*100, col="red")
lines(nvalues, errors3*100, col="green")
# legend
legend("topright",legend=c("Average Error","Varying all","Varying beta","Varying lambda"),
       text.col=c("black","blue","red","green"),pch=c(16,15),col=c("black","blue","red","green"))










# varying DOWN by an order of magnitude


# errors found for n up to 2500 (alpha <- 0.225, beta <- 0.3, lambda <- 0.1) changing all down
errors0=c(1.78379895, 0.41999919, 0.22310964, 0.28748660, 0.18328143, 0.12791260, 0.13458584, 0.14731262, 0.46557044, 0.03316837,
          0.19090959, 0.26239480, 0.04793602, 0.05219560, 0.13361894, 0.04684903, 0.09681905, 0.03492994, 0.03054034, 0.02735415,
          0.14409236, 0.03980513, 0.03666408, 0.07463000, 0.04672714, 0.10935067, 0.01272532, 0.09748947, 0.03140024, 0.06118681,
          0.01680504, 0.12986691, 0.09236004, 0.02121779, 0.06246884, 0.05917334, 0.06520968)

# errors found for n up to 2500 (alpha <- 0.225, beta <- 3, lambda <- 1) changing alpha
errors1=c(6.400637e+00, 5.613634e+03, 1.531801e+00, 4.291560e-01, 3.645569e-01, 6.130213e-01, 4.640072e-01, 1.767041e+00, 2.083731e+00,
          4.485721e-01, 6.398394e-01, 3.290168e-01, 7.191689e-02, 1.275403e+00, 1.460020e-01, 4.570586e-01, 4.098739e-01, 1.679925e-01,
          8.645303e-01, 3.656308e-01, 1.114110e-01, 3.050284e-01, 1.649805e-01, 1.118920e+00, 2.270381e-01, 1.886771e-01, 4.652195e+00,
          1.182026e+00, 1.355609e-01, 2.079126e-01, 1.427573e-01, 7.586800e-01, 2.139195e-01, 2.858915e-01, 1.280236e-01, 1.891759e-01,
          2.353148e-02)

# errors found for n up to 2500 (alpha <- 2.25, beta <- 0.3, lambda <- 1) changing beta
errors2=c(13.6667567,  4.3936574,  2.4959747,  1.9743257,  3.9459114,  1.0964327,  1.9282714,  0.8524422,  2.5792072,  0.7997612,
          2.5207133,  2.4957801,  2.8411803,  1.9216438,  1.8984646,  1.6018619,  4.7405027,  0.7216055,  2.0838040,  1.7166778,
          0.6594127,  1.8988483,  1.4188028,  2.0513108,  0.6892958,  1.8047611,  2.3168749,  0.7614604,  0.2014355,  0.4737996,
          1.4013568,  1.8753846,  1.1810689,  0.8312702,  0.5767367,  2.7148810,  2.1231275)

# errors found for n up to 2500 (alpha <- 2.25, beta <- 3, lambda <- 0.1) changing lambda
errors3=c(1.67875685, 0.22673024, 0.09369349, 0.23601014, 0.25481094, 0.10883754, 0.13027954, 0.15000341, 0.12834472, 0.15053676,
          0.12216558, 0.21737891, 0.11437908, 0.07141566, 0.02962836, 0.13317320, 0.03253026, 0.05910694, 0.09212685, 0.12342043,
          0.08420388, 0.02341176, 0.12694428, 0.03669339, 0.05941969, 0.02080171, 0.03383746, 0.02327558, 0.12324168, 0.06940207,
          0.07545020, 0.07686385, 0.01197320, 0.05991222, 0.04772412, 0.01421757, 0.01776045)


plot(nvalues, avgerror*100, type="l", ylim=c(0,300), main="Estimate Error vs Number of Arrivals (varying parameters)", 
     ylab="Percentage Error (Average for each parameter)", xlab="Number of Arrivals", yaxt='n')
g<-c(0,50,100,150, 200, 250, 300)
axis(2, at = g, labels=c("0%","50%","100%","150%","200%","250%","300%"), las=2) 
lines(nvalues, errors0*100, col="orange")
lines(nvalues, errors1*100, col="blue")
lines(nvalues, errors2*100, col="red")
lines(nvalues, errors3*100, col="green")
# legend
legend("topright",legend=c("Average Error","Varying all","Varying alpha","Varying beta","Varying lambda"),
       text.col=c("black","orange","blue","red","green"),pch=c(16,15),col=c("black","orange","blue","red","green"))






# the interaction between alpha and beta is very important, 
# lets vary them together and see how this affects the grpah
errors4=c(4.976480e+04, 5.829988e+01, 1.625796e+00, 2.702212e+01, 2.101953e+00, 6.774721e-01, 6.110351e-01, 2.698169e-01,
          2.487100e-01, 7.436398e-01, 5.719675e-01, 9.801478e-01, 1.422795e-01, 2.400470e-01, 9.864547e-02, 1.883747e-01,
          1.540087e-01, 2.176920e-01, 3.015003e-01, 1.566625e-01, 3.527257e-01, 3.223519e-01, 1.661477e-01, 1.125593e-01,
          2.663581e-01, 9.779508e-02, 1.821543e-01, 2.312394e-01, 1.496678e-01, 8.287792e-02, 1.599363e-01, 3.345296e-01,
          4.119845e-01, 1.047728e-01, 1.139672e-01, 7.912029e-02, 1.999210e-02)

spareerrors4=c(4.69532789,  3.90487806,  0.62282556, 12.06193081,  0.09583163,  0.45841241, 21.33583808,  0.12325331,  1.04544311,
               0.92611720,  0.49518374,  0.16561355,  0.87703612,  0.50029156,  0.24341593,  0.20935443,  0.23063612,  0.40310554,
               0.19511563,  0.29045385,  0.22765856,  0.16817100,  0.41026248,  0.09647340,  0.10624659,  0.16167166,  0.32440687,
               0.20830769,  0.05552831,  0.14332576,  0.02484425,  0.53936553,  0.17845830,  0.11164640,  0.22699157,  0.22087208,
               0.09855071) # this was another attempt at varying alpha and beta, but i prefer the look of the first one


# plot just avg, vary all, and vary lambda
plot(nvalues, avgerror*100, type="l", ylim=c(0,300), main="Estimate Error vs Number of Arrivals (varying parameters)", 
     ylab="Percentage Error (Average for each parameter)", xlab="Number of Arrivals", yaxt='n')
g<-c(0,50,100,150, 200, 250, 300)
axis(2, at = g, labels=c("0%","50%","100%","150%","200%","250%","300%"), las=2) 
lines(nvalues, errors0*100, col="orange")
lines(nvalues, errors3*100, col="green")
lines(nvalues, errors4*100, col="#9966CC")
# legend
legend("topright",legend=c("Average Error","Varying all","Varying lambda","Varying alpha and lambda"),
       text.col=c("black","orange","green","#9966CC"),pch=c(16,15),col=c("black","orange","green","#9966CC"))










# plot just avg, vary all, and vary lambda
plot(nvalues, avgerror*100, type="l", ylim=c(0,100), main="Estimate Error vs Number of Arrivals (varying parameters)", 
     ylab="Percentage Error (Average for each parameter)", xlab="Number of Arrivals", yaxt='n')
g<-c(0,50,100,150, 200, 250, 300)
axis(2, at = g, labels=c("0%","50%","100%","150%","200%","250%","300%"), las=2) 
lines(nvalues, errors0*100, col="orange")
lines(nvalues, errors3*100, col="green")
# legend
legend("topright",legend=c("Average Error","Varying all","Varying lambda"),
       text.col=c("black","orange","green"),pch=c(16,15),col=c("black","orange","green"))