Loading required libraries
library(MASS)
library(scoringTools)
library(rpart)
Mean vectors and variances for each class
d <- 2
mu0 <- array(0, c(1, d))
mu1 <- array(1, c(1, d))
sigma0 <- diag(1, d, d)
sigma1 <- diag(1, d, d)
Data generation
m <- 10000
set.seed(21)
y <- rbinom(m, 1, 0.5)
data <- array(0, c(m, d + 1))
x <- array(0, c(m, d))
x[y == 0, ] <- mvrnorm(n = sum(y == 0), mu0, sigma0)
x[y == 1, ] <- mvrnorm(n = sum(y == 1), mu1, sigma1)
data <- as.matrix(cbind.data.frame(y = y, x = x))
rm(x)
rm(y)
train <- as.data.frame(data[1:(m / 2), ])
test <- as.data.frame(data[(m / 2 + 1):m, ])
train[, "y"] <- as.factor(train[, "y"])
Oracle
modele_complet_arbre <- rpart(y ~ ., data = train, method = "class")
modele_complet_reglog <- glm(y ~ ., data = train, family = binomial(link = "logit"))
Loop over cut-off values
list_gini_arbre <- list()
list_gini_reglog <- list()
for (i in seq(0.2, 0.7, 0.05)) {
ind_refuses_arbre <- predict(modele_complet_arbre, train)[, 1] < i
ind_refuses_reglog <- predict(modele_complet_reglog, train, type = "response") < i
train_partiel_arbre <- train[!ind_refuses_arbre, ]
train_partiel_reglog <- train[!ind_refuses_reglog, ]
modele_partiel_arbre <- rpart(y ~ ., data = train_partiel_arbre, method = "class")
modele_partiel_reglog <- glm(y ~ ., data = train_partiel_reglog, family = binomial(link = "logit"))
list_gini_arbre <- append(list_gini_arbre, normalizedGini(test[, "y"], predict(modele_partiel_arbre, test)[, 2]))
list_gini_reglog <- append(list_gini_reglog, normalizedGini(test[, "y"], predict(modele_partiel_reglog, test, type = "response")))
}
Figures
plot(seq(0.2, 0.7, 0.05),
list_gini_arbre,
ylim = c(0.35, 0.7),
xlab = "Rejection rate",
ylab = "Gini",
col = "red"
)
lines(seq(0.2, 0.7, 0.05), list_gini_reglog, ylim = c(0.35, 0.7), type = "p", col = "blue")
legend(1, 0.45,
pch = c(1, 1), lty = c(1, 1),
col = c("red", "blue"),
legend = c("Decision tree", "Logistic regression"),
cex = 1
)
![plot of chunk plot](data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAfgAAAH4CAIAAAApSmgoAAAACXBIWXMAAAsSAAALEgHS3X78AAAgAElEQVR4nO3de1yUZf7/8Q8wIqaAoySeqMQCQ9mfhxwEzPGAuXlEkfXQ4VG5pIGW1qOjrm3YYvXVbdtaSdss89vJUyFa4LaRfjVNNM+IGuKAhRAoAjoIyvz+GNcDGdrGfU1e83r+dd/3zPi+7sK3N9fcBw+HwyEAAH15unoAAABjUfQAoDmKHgA0R9EDgOYoegDQHEUPAJqj6AFAcxQ9AGiOogcAzVH0AKA5ih4ANEfRA4DmKHoA0BxFDwCao+gBQHMUPQBojqIHAM1R9ACgOYoeADRH0QOA5ih6ANAcRQ8AmqPoAUBzFD0AaI6iBwDNUfQAoDmKHgA0R9EDgOYoegDQHEUPAJqj6AFAcxQ9AGiOogcAzVH0AKA5ih4ANEfRA4DmKHoA0BxFDwCao+gBQHMUPQBojqIHAM1R9ACgOYoeADRH0QOA5ih6ANAcRQ8AmqPoAUBzFD0AaI6iBwDNUfQAoDmKHgA0R9EDgOYoegDQHEUPAJqj6AFAcxQ9AGiOogcAzZlcPYBfprS0NCsry9WjAIBG5unpOXLkyCZNmhjxh19nRf/ll1+mpaVZrVZXDwQAGtN7773XvXv3zp07G/GHX2dFLyLR0dEPP/ywq0cBAI1p69atxv3hqufoKysrHQ6H4lAAcGcqit5ut6ekpISGhvr4+Pj5+Xl7e4eEhCQnJ9fU1ChIBwA3p6Lok5KSsrOzFy9eXFRUVFtbW1JSsmTJkpycnMTERAXpAODmVMzRp6enFxQUNGvWzLlqNpsjIyMtFktwcLCCdABwcyqO6IOCgjIzM+tt3LBhQ0BAgIJ0AHBzKo7oU1NT4+PjZ8+eHRYW5uvrW1VVlZubW1pampaWpiAdANyciqKPiIjIz8/Pysqy2WxlZWVmszkhIcFqtXp5eSlIBwA3p6LoN23aZLFYYmJiFGQBAOpRMUfft2/fiIgIQy8HAAD8HEUXTC1cuPCZZ5655557cnNz1SQCAJwUFX1YWNiXX345atSo+Pj4QYMGLV68+OjRo2qiAcDNKb3XzR/+8IexY8d+/vnnS5YsefTRR81mc2Fh4c+9+Ysvvli+fHm9jXv27OnSpYvBwwQArai+qZmnp+ewYcOGDRtWW1vb8Kx9r169zGZzvY3z58+vrKw0coAAoBsVRT9v3jxvb+96G5s0aRIdHd3Ap8xmc69eveptDAwMLCoqauTxAYDWVMzRP/HEEz+9m355eflPL5cFADQ6lz1KMCcnJzY21lXpAOA+XFb0UVFRdrvdVekA4D548AgAaI4HjwCA5njwCABojgePAIDmePAIAGiOB48AgOZ48AgAaE7RvW68vLx48AgAuITLLpgCAKhB0QOA5ih6ANAcRQ8AmqPoAUBzFD0AaI6iBwDNUfQAoDmKHgA0R9EDgOYoegDQHEUPAJqj6AFAcxQ9AGiOogcAzVH0AKA5RQ8eUc9ulw8+kIIC6d5dYmPFw0NF6JEjsmmTmM0SEyPe3ioS6+pk2zYpL5fevcVsVpF47pwsWya7dsktt8h990nz5ipCP/5YVq4UEYmPl/h4FYm5uTJvnths0qOHPPWUKHi8cW2tpKbKl1+K2SxTpkhEhOGJIpKXJ++8IxUVMnCgxMaqSKypkXfeOf/zk5Cg6Id2715JSxMPDxk1Srp2VZEoIvv3i80m4eHSoYOixAboeURfUSF33SVnz0pMjOzeLePHi8NheOi770pSktjtsnOnDBokx44Znnj8uAweLB99JNu2yejRouARvHV1EhcnhYUyerQ0bSpDhsjJk4aHvvKKbN4sb7whb7wh//d/Mn++4Yl5eZKQII8+Kp9/LnfdJWPGyOnThoc+8ICYTPLuu/Lss5KcLFlZhidu3y4PPSSDBsnDD8uWLfL444Ynnj0ro0ZJXZ1MmyadO8uwYVJWZnjoqlXyzDNisUjv3vLkk/LJJ4Yn1tXJ/ffLSy/Jzp0yZYq8/LLhiVfnuK5Mnz593LhxV33b3LmOZcsurj71lOOLLwwclcPhOHnS0a+f4+zZ86tbtzoeesjYRIfD8cgjjq++Or9cXe2IjnacOmVsYkaG47nnLq6uWOGYO9fYRIfD0bevo67u/HJdnaNvX8MTn3vOsX79xdX58x2rVhmbWFDguPTnuqzMMXy4sYkOhyM21vH99xdXhw93FBcbm7h2rePPf764unq1Y84cYxMdDkd0tMNuP798+rSKn59//tMxb97F1QkTHDt2XP1TkyZN+u677wwakp5H9Lm5YrFcXI2IkJwcYxNzciQqSi487bx3bzl82NhEEdm/X6zW88tNm0p0tOzda2ziwYPSq9fF1TvukNxcYxPPnJHmzS/OvHl4SPPmcuaMsaFFRRIUdHH1ppvkhx+MTfz+e+nU6eJqq1ZSVWVsooiUl0v79hdXw8MlL8/YxLw8CQ9XmnjmjLRoIT4+51ebNVPx8/P11zJ69MXV0aPl66+NTbwqPYs+JER27Li4+u23EhpqbGLHjpKff3G1rEx8fY1NFJFmzS6rgx9+kMBAYxNDQ2XbtourW7fK7bcbm9i0qdTUyI8/nl/98UeprZWmTY0NjYiQtWsvrq5ZI336GJvYtat8/bWcO3d+detW6dzZ2EQRadtWDhw4v+xwyDffSJcuxiY6d/OCTZukWzdjE5s2lVOnLjZ7dbWcOmX4z0+bNpfN3Cr4i3lVen4Z+8gjMny4VFVJWJhkZcnevZKcbGxix47SpIm8/rrce68UF8uMGfLMM8YmisjkyfLHP8qCBeLvL8uXy+nTcvPNxiYOHixvvikpKTJggOzdK0uXypo1xiaKyNy5Mnq0xMeLwyErVsirrxqeOGmSTJggOTnStausXy+3337Z7zFG8PWVSZNk+HAZPVpKSiQzU1asMDZRRObMkfvukylTpHVr+fBDGTrU8K9GBw6UxYvluefEapV9+2TNGhU/P48+KmPHnv8G4q9/lenTDU984AGZPFnefls6d5avv5YVK+Tzzw0PbZiHQ8HXlI1nxowZRUVFH3300VXfWVkp774rR45Ir14ybtzFSRXj1NbKokXy739Lq1YyebL07m14ooisXStvvinV1RIdLU8+qeIcmLo6WblSdu6UTp3knnukWTPDE0Xk5EnJyhIPD+nfX/z9VSSKyPbtkp8v/+//yW23KUosKJBNm6RVKxk4UJo0UZF44oSsXCnl5TJokPTooSJRRDIyZMcO6dRJ4uIU7ebu3ee/gx0z5rK5I+Ps3CmvvCLHjsntt8vMmZdNkf2cP/7xj88++2xnY36V07boAeA6YmjRu2CO/vjx46dOnVKfCwDuSUXRjxgxorq6WkTy8vIsFkvr1q39/f2HDh1aXFysIB0A3JyKol+zZs25c+dEZPLkyVFRUTabrbCwsFu3blOnTlWQDgBuTulZN9u3b8/IyDCZTCIyd+7c9tfyDQUA4NdRNEfvnKUJDg4uKChwbjl06JC/svMnAMCNqSj6fv36Wa1WX1/fw4cPT5s2TUQ2btxotVqdywAAQ6mYulm/fr2IVFdXHzlypKKiQkRMJtPChQtj1dwuDwDcm7o5eh8fny5duojIiRMnunbt2sfoi8oBACKiZurGZrPde++9hYWFNpstKiqqdevWLVu2HDhw4NGjRxWkA4CbU1H0EyZM6Ny5c7t27R599NG+fft+//33R48e7d2795QpUxSkA4CbUzF1k5ub+9VXX5lMpt27d7///vstWrQQkblz5wa6/JZuAOAGVBR9nz59Vq5cOWHChKioqHXr1o0ZM0ZE1q1b17FjxwY+9e233/7rX/+qt3Hbtm2clAkAv4iKon/rrbfuvvvu+fPnt23bNi4urmfPnnV1dUVFRatXr27gU/7+/sHBwT/d6LzeCgBwjVSUZocOHXbt2rV+/fpDhw7169fP39//lltuGTRoUMOV3blz55/eyO3rr78uKioycrAAoBtFR8ceHh79+/fv37//hS3l5eXffPPNkCFD1AwAANyWyx4lmJOTwwVTAKCAy4o+KirKbre7Kh0A3Ifqoq+srLy+nmkFANc7FUVvt9tTUlJCQ0N9fHz8/Py8vb1DQkKSk5NramoUpAOAm1NR9ElJSdnZ2YsXLy4qKqqtrS0pKVmyZElOTk5iYqKCdABwcyrOuklPTy8oKGjWrJlz1Ww2R0ZGWiyWn54mDwBodCqO6IOCgjIzM+tt3LBhQ0BAgIJ0AHBzKo7oU1NT4+PjZ8+eHRYW5uvrW1VVlZubW1pampaWpiAdANyciqKPiIjIz8/Pysqy2WxlZWVmszkhIcFqtXp5eSlIBwA3p+jKWC8vr5iYGDVZAIBLueyCKQCAGhQ9AGiOogcAzVH0AKA5ih4ANEfRA4DmKHoA0BxFDwCao+gBQHMUPQBojqIHAM1R9ACgOYoeADRH0QOA5ih6ANAcRQ8AmqPoAUBzFD0AaI6iBwDNUfQAoDmKHgA0R9EDgOYoegDQHEUPAJpTXfSVlZUOh0NxKAC4MxVFb7fbU1JSQkNDfXx8/Pz8vL29Q0JCkpOTa2pqFKQDgJtTUfRJSUnZ2dmLFy8uKiqqra0tKSlZsmRJTk5OYmKignQAcHMmBRnp6ekFBQXNmjVzrprN5sjISIvFEhwcrCAdANyciiP6oKCgzMzMehs3bNgQEBCgIB0A3JyKI/rU1NT4+PjZs2eHhYX5+vpWVVXl5uaWlpampaUpSAcAN6ei6CMiIvLz87Oysmw2W1lZmdlsTkhIsFqtXl5eCtIBwM2pKHoR8fLyiomJcS6fOHHCZDLR8gCghoo5epvNdu+99xYWFtpstqioqNatW7ds2XLgwIFHjx5VkA4Abk5F0U+YMKFz587t2rV79NFH+/bt+/333x89erR3795TpkxRkA4Abk7F1E1ubu5XX31lMpl27979/vvvt2jRQkTmzp0bGBioIB0A3JyKou/Tp8/KlSsnTJgQFRW1bt26MWPGiMi6des6duzYwKcKCwsPHDhQb2NBQUFdXZ2BYwUA7ago+rfeeuvuu++eP39+27Zt4+LievbsWVdXV1RUtHr16gY+VVRUtH379nobS0pK/Pz8jBwsAOhGRdF36NBh165d69evP3ToUL9+/fz9/W+55ZZBgwaZTA2lWywWi8VSb+OxY8eKioqMHCwA6EbR6ZUeHh79+/fv37//hS3l5eXffPPNkCFD1AwAANyWy+5Hn5OTExsb66p0AHAfLiv6qKgou93uqnQAcB88eAQANMeDRwBAczx4BAA0x4NHAEBzPHgEADTHg0cAQHM8eAQANOeCB48AAFRy2QVTAAA1KHoA0BxFDwCao+gBQHMUPQBojqIHAM1R9ACgOYoeADRH0QOA5ih6ANAcRQ8AmqPoAUBzFD0AaI6iBwDNUfQAoDmKHgA0R9EDgOYUPWHKNXbulIICCQ+XTp0UJZaUSHa2tGwpffoID0oE8Nug6RH9uXMyfrwsWCCHD8sTT0hysorQFStkwgTJzZW0NImJkbIyFaEAcDWaHtG//bZERspjj4mITJ8u998v334rPXsamFhVJX/7m2RlSZMmIiL/938ya5akphqYCADXRtMj+i1bZOTIi6sjRsiWLcYm7tsnffueb3kRufNOyc01NhEAro2mRd+2rXz//cXVo0elbVuliRUV0ry5sYkAcG00LfqHHpLnnpP9+6WuTtavl08+kSFDjE28+Wax22XpUqmtleJiefBBSUoyNhEAro2mc/S33iqpqfLyy1JQIN27y/LlKo6v33tPXn1Vhg0TPz9JTJRBgwxPBIBr4IKiP378eNOmTZsb3bxdu8q77xobUc8NN8jMmTJzptJQALgaFVM3I0aMqK6uFpG8vDyLxdK6dWt/f/+hQ4cWFxcrSAcAN6ei6NesWXPu3DkRmTx5clRUlM1mKyws7Nat29SpUxWkA4CbUzp1s3379oyMDJPJJCJz585t3769ynQAcE+KzrpxztIEBwcXFBQ4txw6dMjf319NOgC4MxVF369fP6vV6uvre/jw4WnTponIxo0brVarcxkAYCgVUzfr168Xkerq6iNHjlRUVIiIyWRauHBhbGysgnQAcHPq5uh9fHy6dOkiIsePHw8PD+/Tp4+yaABwZ5xeCQCa4/RKANAcp1cCgOYUFX1xcXFwcLDz9Mrg4GC5htMrq6qqSkpK6m08efKk85cDAMA1UlH0ztMry8vLTSbTtGnT1q5du3Hjxri4uFmzZjXwqU2bNq1atarexm3btrU1+obDAKCX3+7plUOGDBnyk3sLz5gxo6ioyLihAoB+XHB6pYj06dOnvLw8MzPzp1UOAGhcLnvwSE5ODhdMAYACLiv6qKgou93uqnQAcB+qi76ystLhcCgOBQB3pqLo7XZ7SkpKaGioj4+Pn5+ft7d3SEhIcnJyTU2NgnQAcHMqij4pKSk7O3vx4sVFRUW1tbUlJSVLlizJyclJTExUkA4Abk7FWTfp6ekFBQXNmjVzrprN5sjISIvF4rxyCgBgKBVH9EFBQZmZmfU2btiwISAgQEE6ALg5FUf0qamp8fHxs2fPDgsL8/X1raqqys3NLS0tTUtLU5AOAG5ORdFHRETk5+dnZWXZbLaysjKz2ZyQkGC1Wr28vBSkA4CbU3RlrJeXV0xMjJosAMClXHbBFABADYoeADR32dRN7969//GPfyQlJf30fdnZ2aqGBABoTJcV/RtvvNGlS5c333zTVaMBADS6y4o+IiJCRHr16uWiwQAAGt8Vzrr58ssvn3/++fLy8ks37tmzR9WQAACN6QpF/9xzzw0fPnzMmDGc5w4AGrhC0R8+fHjmzJkeHh7qRwMAaHRXOL0yJibm9ddfr6qqUj8aAECju8IRfUFBwYcffvjkk0926tTpwnH9/v371Q4MANA4rlD0CxcuVD8OAIBBrlD0Xbt2VT8OAIBBuDIWADTHlbEAoLn6V8bW1tbecMMNt99++1/+8pezZ886t3fs2JHLZQHgOnXZ6ZXHjh3r0aPHa6+9JiIvvvji6dOnV69e/fbbb4eEhLhoeACAX+uyop85c2Z0dPSFqZuXX345Ozt70KBB3333nSvGBgBoBJcVfUZGxpQpUy572dMzOTk5NTVV7agAAI3msqIvLy8PDQ11Lu/atcu5EBAQcPDgQdXjAgA0ksuKPiQkZMuWLReWnQtZWVldunRRPS5co/R0GTZMYmJk9mzhrhUAruSyon/88ccTExMvnZH/+uuvH3nkkalTpyofGK7BqlXy0Ufy/vuybp2Eh8t997l6QAB+iy47vfK+++47depUz549f/e73914443ffffdkSNHnn/++XvvvddV40ND/vlPWb5cmjcXEYmPlzVr5MgRueUWF48KwG9M/VsgTJkyJTY29ptvvsnPz3/wwQcjIyNvvPFGl4wMV2e3n295p7ZtpaSEogdQzxXuddO2bdtRo0apHwp+sa5d5csvZeBAEZHqatm4Uf78ZxcPCcBvzxWK3lCVlZUtWrTgqSaN48UXJT5ePv1UWreWf/9bnn1WmjVz9ZgA/OZc4cEjjc5ut6ekpISGhvr4+Pj5+Xl7e4eEhCQnJ9fU1ChI11nLlrJunTz4oNx5p6Sny/Dhrh4QgN8iFUWflJSUnZ29ePHioqKi2trakpKSJUuW5OTkJCYmKkjXnIeH9OghAweKv7+rhwLgN0rF1E16enpBQUGz/8wqmM3myMhIi8USHBysIB2N7ORJee012bFDOnWS6dPlpptUhH76qSxfLh4e8oc/yMiRKhIBjag4og8KCsrMzKy3ccOGDQEBAQrS0Ziqq2XUKPnd72TxYomLk/Hj5ehRw0NffVXWrZNXXpGXXpLPPpPXXjM8EdCLiiP61NTU+Pj42bNnh4WF+fr6VlVV5ebmlpaWpqWlKUhHY8rIkKFDJTZWRCQ6WpKT5a235IUXjA1dtUrWrxdPTxGRBQvEapXHHjM2EdCLiqKPiIjIz8/Pysqy2WxlZWVmszkhIcFqtXp5eSlIR2M6elQ6d7642rmzfPCBsYlnzkjz5udbXkQ8PaV5czlzRpo2NTYX0Iii0yu9vLxiYmKcy8ePH2/atCktf13q1UuWLJG4uPOr69aJ0U+kadpUzpyR48elVSsRkbIyqamh5YFfRMUc/YgRI6qrq0UkLy/PYrG0bt3a399/6NChxcXFCtLRmCIjxcdHHnpI3ntPnn5a1q6Vhx82PPTFF2XUKElNlQULJDZW/vIXwxMBvago+jVr1pw7d05EJk+eHBUVZbPZCgsLu3Xrxr3Srkt/+5tMmyZ1dXL33ZKWJk2aGJ4YHS2ffCKtW0tAgHz6qURGGp4I6EXplbHbt2/PyMgwmUwiMnfu3Pbt26tMR6Pp0UN69FCaGBAgf/iD0kRAIyqO6EXEOUsTHBxcUFDg3HLo0CF/rvEBAOOpOKLv16+f1WotLy83mUzTpk1bu3btxo0b4+LiZs2a1cCnVq5ceeHptRccPHjw1ltvNXKwAKAbFUW/fv16Eamurj5y5EhFRYWImEymhQsXxjpPx/4ZcXFxcRfO7viPGTNmFBUVGTdUANCPujl6Hx+fC48k7NOnT3l5eWZm5pAhQ5QNAADck6I5+p/Kyclp+IgeANAoXFb0UVFRdrvdVekA4D5UF31lZaXD4VAcCgDujAePAIDmePAIAGiOB48AgOZ48AgAaI4HjwCA5njwCABozgUPHgEAqOSyC6YAAGpQ9ACgOYoeADRH0QOA5ih6ANAcRQ8AmqPoAUBzFD0AaI6iBwDNUfQAoDmKHgA0R9EDgOYoegDQHEUPAJqj6AFAcxQ9AGiOogcAzVH0AKA5ih4ANEfRA4DmKHoA0BxFDwCao+gBQHMUPQBojqIHAM25oOiPHz9+6tQp9bkA4J5UFP2IESOqq6tFJC8vz2KxtG7d2t/ff+jQocXFxQrSAcDNqSj6NWvWnDt3TkQmT54cFRVls9kKCwu7des2depUBekA4OZMKsO2b9+ekZFhMplEZO7cue3bt1eZDgDuSdEcvXOWJjg4uKCgwLnl0KFD/v7+atIBwJ2pKPp+/fpZrVZfX9/Dhw9PmzZNRDZu3Gi1Wp3LAABDqZi6Wb9+vYhUV1cfOXKkoqJCREwm08KFC2NjYxWkA4CbUzdH7+Pj06VLF+eyw+Gg5QFADddcMHXXXXe5JBcA3JCKop8+fXqzy1VVVTkXFKQDgJtTMXWTkpJit9v37duXmprapk0bEQkODj58+LCCaACAiqK/4YYbFi5cuHr16okTJ7700kvDhg3z9PQMDAxs+FMrV65888036208ePDgrbfeathIAUBD6r6MHTlypMViefDBB9PT06/l/XFxcXFxcfU2zpgxo6ioyIDRAYC2lH4Z27Zt288++yw8PJxTbgBAGdVn3Xh4eCQlJS1durS8vDwzM1NxOgC4IZfdjz4nJ4fjegBQwGVFHxUVZbfbXZUOAO5DddFXVlY6HA7FoQDgzlQUvd1uT0lJCQ0N9fHx8fPz8/b2DgkJSU5OrqmpUZAOAG5ORdEnJSVlZ2cvXry4qKiotra2pKRkyZIlOTk5iYmJCtIBwM2pOI8+PT29oKDgwg0PzGZzZGSkxWIJDg5WkA4Abk7FEX1QUNBPz6TcsGFDQECAgnQAcHMqjuhTU1Pj4+Nnz54dFhbm6+tbVVWVm5tbWlqalpamIB0A3JyKoo+IiMjPz8/KyrLZbGVlZWazOSEhwWq1enl5KUgHADen6F43Xl5eMTExarIAAJdy2QVTAAA1KHoA0BxFDwCao+gBQHMUPQBojqIHAM1R9ACgOYoeADRH0QOA5ih6ANAcRQ8AmqPoAUBzFD0AaI6iBwDNUfQAoDmKHgA0R9EDgOYoegDQHEUPAJqj6AFAcxQ9AGiOogcAzVH0AKA5ih4ANKe66CsrKx0Oh+JQAHBnKorebrenpKSEhob6+Pj4+fl5e3uHhIQkJyfX1NQoSAcAN6ei6JOSkrKzsxcvXlxUVFRbW1tSUrJkyZKcnJzExEQF6QDg5kwKMtLT0wsKCpo1a+ZcNZvNkZGRFoslODhYQToAuDkVR/RBQUGZmZn1Nm7YsCEgIEBBOgC4ORVH9KmpqfHx8bNnzw4LC/P19a2qqsrNzS0tLU1LS1OQDgBuTkXRR0RE5OfnZ2Vl2Wy2srIys9mckJBgtVq9vLwUpAOAm1NR9CLi5eUVExNzYXXz5s20PACo4ZoLpu666y6X5AKAG1JR9NOnT292uaqqKueCgnQAcHMqpm5SUlLsdvu+fftSU1PbtGkjIsHBwYcPH1YQDQBQUfQ33HDDwoULV69ePXHixJdeemnYsGGenp6BgYENf2rlypVvvvlmvY0HDx687bbbDBspAGhI0ZexIjJy5EiLxfLggw+mp6dfy/vj4uLi4uLqbZwxY0ZRUZEBowMAbSn9MrZt27afffZZeHh4bGysylwAcGeqz7rx8PBISkpaunRpeXn5Ty+XBQA0Opfdjz4nJ4fjegBQwGVFHxUVZbfbXZUOwL1UV0t1tdJEh0NOnlSa+PPUfRnrVFlZ2aJFCw8PD8W5ANzUjz/KlCly6pSISIsW8uabouB2in//u3z4oQQGSmmpzJwpd99teGKDePAIAK1NnSqPPy4ZGZKRIY89JlOnGp74+eeyd69s2iSffirr1slf/yo2m+GhDeLBIwD0dfasHD8u0dHnV++8U8rK5OxZY0PXrJHp08XTU0Tkhhvk/vvl3/82NvFqePAIAH39dJbY4bjCxsbl5SV1dRdXz50TV9/DkQePANCXl5e0ayeff35+de1a6dDB8NodPVrmzhXn1PTx47JkiQwebGzi1fDgEQBae/11mT5dXnlFRCQ4WF5/3fDEAQOksFAGDRJvbxGROXOkfXvDQxvEg0cAaM3fX955RxwOkSvN5Bjk/vvl/vulru78TL2ruebBIwCglEtO6f5ttLyoP48euD6cPi0ffCA2m3TvLqNHK/obu369ZGVJy5YycaK0aaMisbhYPvhAKvH4GjYAABDzSURBVCpk4EC5804ViXV1snKl7N4tt9wiEycKD6VQ4rfyDw7wG1JRIUOGSF2d/P73sn+/jBt3/hd/Q73wgnzwgQwYIB06SGys5OYanrh/v4wZI0FB0r+/vPeezJljeKLDIfHxcuCA/P73cu6cDBkilZWGh4IjeuAKFiyQxx6TsWNFRKKj5Zln5N//FkPnHktLZfNmycg4v9qzpzz3nHz8sYGJIvL887J0qTjPcrZaZcgQOX5cWrUyMPGLLyQkRGbNEhGJjpaWLSU1VZ56ysBEiAhH9MAV5OZK794XVy0W2b/f2MRDh6Rnz4urnTtLSYmxiSLy449y6bUsvXrJgQPGJtb7DxsRITk5xiZCRCh64ApCQuTbby+ubt8uoaHGJt56q+zceXH1yBEV92Np3fqyS/N37BCjH98WEiI7dlxcVfAfFiLC1A1wBYmJMmyYVFRIWJh89ZXk5Bg+f33jjdK9u0ybJvfcIz/+KC+9JIsWGZsoIs8/LxMnyrPPSkCALF0qd9xh+L8ugwfLggUyb55YrbJvn/zzn7JmjbGJEBGO6IEraNlSMjOlqkqWLZOOHWXFChVn3aSkyMiR8umnsn+/LF8uXbsanhgeLsuWyb598umnMmaMii9jPT1l1Spp316WLZPTpyUzU1q2NDwUHNEDV9aihSQlqQ4dPFj1tfIdOsjTTytN9PKSiRNl4kSloW6PI3oA0BxFDwCao+gBQHMUPQBojqIHAM1R9ACgOYoeADRH0QOA5ih6ANAcRQ8AmqPoAUBzFD0AaI6iBwDNUfQAoDmKHgA0R9EDgOZcUPTHjx8/deqU+lwAcE8qin7EiBHV1dUikpeXZ7FYWrdu7e/vP3To0OLiYgXpAODmVBT9mjVrzp07JyKTJ0+Oioqy2WyFhYXdunWbOnWqgnQAcHNKnxm7ffv2jIwMk8kkInPnzm3fvr3KdABwT4rm6J2zNMHBwQUFBc4thw4d8vf3V5MOAO5MRdH369fParX6+voePnx42rRpIrJx40ar1epcBgAYSsXUzfr160Wkurr6yJEjFRUVImIymRYuXBgbG6sgHQDcnLo5eh8fny5dujiX77jjDudMPQDAaCqmboqLi//4xz8OGDBg0aJF3377bUhIiI+PT79+/Ww2m4J0AHBzKop+0qRJJ06cSExM/Oijj0aMGDFnzpzDhw93796d0ysBQAEV8yebNm3Kz89v2bLljTfeOGPGjHHjxonIiy++GBwcrCAdANyciqL39fU9duxYy5Ytf/e7382fP9+50Wazmc3mBj71xRdfLF++vN7GzZs3d+jQwaiBAoCOVBT9M888079//3vvvXfevHkDBw4UkUWLFv3tb39LSEho4FO9evX66b8EQUFBTZs2NXCsAKAdFUWfmJhotVrz8vIubKmsrHzxxRfHjBnTwKfMZnOvXr3qbczLyystLTVklACgKUXnOHbt2rVr164XVp944ony8vLMzMwhQ4aoGQAAuC2X3Y8+JyeHC6YAQAGXFX1UVJTdbndVOgC4D9VFX1lZ6XA4FIcCgDtTUfR2uz0lJSU0NNTHx8fPz8/b2zskJCQ5ObmmpkZBOgC4ORVFn5SUlJ2dvXjx4qKiotra2pKSkiVLluTk5CQmJipIBwA3p+Ksm/T09IKCgmbNmjlXzWZzZGSkxWLhylgAUEDFEX1QUFBmZma9jRs2bAgICFCQDgBuTsURfWpqanx8/OzZs8PCwnx9fauqqnJzc0tLS9PS0hSkA4CbU1H0ERER+fn5WVlZNputrKzMbDYnJCRYrVYvLy8F6QDg5hRdGevl5RUTE6MmCwBwqevsMU/+/v4pKSmffPLJNb5/8+bN2t8Era6urqamxsfHx9UDMVZtba2INGnSxNUDMVZ1dbW3t7enp8uuZFTj9OnTN9xwg6tHYbi6uro77rjjGt9cWFh44YyVRueh9+VLAwYMyMrKcvUojLVnz5633nrr73//u6sHYqx3331XRB544AEXj8Ngjz76aEJCQnh4uKsHYix3+Ispv6Xd1PzAAQBA0QOA5ih6ANAcRQ8AmqPoAUBz19nplb+UyaT5DoqIp6en9mfjiYibXF7nJv833eEvpvyWdlPz0ysrKyt9fX1dPQpjORyOU6dOtWjRwtUDMZbzptbe3t6uHoixqqqqmjdv7uHh4eqBGMsd/mLKb2k3NS96AID+vyQCgJuj6AFAcxQ9AGiOogcAzVH0AKA5ih4ANEfRA4Dm9Cn6rVu39urV68Ybb5w0aZLdbq/36ueff969e3d/f/++ffvu3bvXJSNsFA3vZmpq6s0332w2m0eMGFFcXOySETaKhnfT6fTp09OnT1c8sEbU8D5arVaP/xgxYoRLRtgoGt7NH374YejQoa1atYqOjrbZbC4ZYaNoYDfnz5/vcblNmzapHp9DC2fPng0KClq1atXp06dHjRr1l7/85dJXCwsLmzdvnpmZWVlZ+eSTT4aFhblqnL9Sw7u5d+/ewMDAnJycEydOjBkzZtKkSa4a56/U8G5e8PDDD/v4+CgeW2O56j527NgxLy+vqqqqqqrKbre7ZJC/XsO7WVdXFx4e/v777585c+ZPf/rTPffc46px/koN72ZNTU3Vf2zdurVHjx41NTWKR6hJ0f/rX/8KDw93Lm/cuDEkJOTSVz/88MO77rrLuVxeXu7h4XHixAnVQ2wMDe9mWlraE0884VxetmxZ//79VY+vkTS8m05paWnR0dHXb9E3vI92u71FixauGFcja3g3s7KyevXq5Vyura0tKChQPb5Gci0/sU4DBgz45ptvVI3rIk2mbmw2W1hYmHM5LCys3u+AI0aM+PDDD53LO3bs6NSpU8uWLVUPsTE0vJsjR46cN2/eiRMnvvrqq9TU1PHjx7tijI2g4d0UkWPHjs2cOTM1NVX50BpNw/t45MgRk8l0xx13mM3mwYMHHz582BVjbAQN7+ahQ4fatWt3zz33BAcHjx071nHd3o7lqj+xTsuWLWvXrp3FYlE4tPM0KfqysrILNw/y9fU9c+ZMZWXlhVebN2/eqlUrEVm9evWECRP++te/umaUv1rDu+m0Y8eOGTNmHDhwoHfv3soH2Dga3k2Hw/Hggw++9NJLbdq0cdEAG0HD+1hRUREVFfXxxx+XlJSEhYVdv/9mN7ybJ06cWLNmzejRo3ft2hUSEqLrbjrV1dXNmTNn1qxZykcnos1tis1mc1VVlXO5srLSZDLVu5vjyZMnp0yZsnPnzhUrVkRHR7tijI3gqrspIgMHDtyxY8fHH388fvz4gwcPKh9jI2h4NxcsWHDzzTcPGzbsuv62ueF9tFgsa9eudS4/8cQTN9988/Hjx50HK9eXhnfT39+/R48eY8eOFZEZM2Y4fx81m82uGeuvcC1/MT///PM2bdrcfvvtykcnos0RfXBw8IEDB5zLubm5nTp1uvRGr7W1tUOGDGnVqtXu3buv35aXq+3m/PnzFy1a5FyOjIw8fPjw2bNnXTDKX63h3dyyZcv777/v6+sbHBxcXV3t6+u7efNmF430v3fVfdy4caNzuUmTJl5eXk2bNnXBKH+1hnfzlltuuXS6xsPD47dzA/dfpOHddFqyZMnEiROVD+0/1H8tYISzZ8926NDhyy+/PHv27Lhx41544QXn9k2bNp04cWL58uW9evWqvoRrR/tfa3g3V65cGRYWtmfPnpMnT06bNu3OO+907Wj/aw3v5oW3HTt27Pr9Mrbhffzkk0/atGmzf//+c+fOPfXUU7///e9dO9r/WsO7WVNTExgYmJGRce7cuaeffnrAgAGuHe1/7ao/sTU1NS1atCgsLHTVCDUpeofDsW3bth49etx0000PPfTQmTNnnBt9fHwyMjKeeeaZev+8lZeXu3a0/7UGdtPhcLz88sudO3du1arV0KFDjxw54tKR/ioN76bTdV30jqvt42uvvXbrrbe2adMmLi7u2LFjLh3pr9Lwbm7ZsqV79+4BAQFDhw51YQ/+eg3v5qZNmzp06ODC4fHgEQDQnCZz9ACAn0PRA4DmKHoA0BxFDwCao+gBQHMUPQBojqIHAM1R9ACgOYoeADRH0QOA5ih6ANAcRQ8AmqPoAUBzFD0AaI6iBwDNUfQAoDmKHgA0R9EDgOYoegDQHEWP60lQUJDJZGrSpInJZPLz8xs2bFhRUdHPvXnfvn3h4eG/6M8fNmyYiHzzzTe9e/f+tWO9WgqgDA8Hx/UkKChoxYoVERERInLixImhQ4eGhYW9/fbbV3zzqVOn9uzZ06dPn2v/85s0aVJbW1tRUZGbm2uxWH7NUM+cOdO0adMGUn7NHw78IhzR43plNpvHjh1bUFDgXF2/fn337t3btWs3bty448ePi8iRI0cSEhJ+7lUR+eijj2677ba2bdtOmTKlurp65MiRZ8+eDQ8P37dvX1JSkvM9K1asCA0NvfHGG+Pi4kpLS0Vk3759ffv2veeeewIDAyMjI7Oysi4d1a5duwYMGDBr1qzo6GgR+fvf/96pUyc/P78+ffocOHBARC6k1NTUXHFUQKOj6HG9ys/P//TTT0eOHCkiZWVlsbGxb7zxxtGjR2+66aYHHnjg0nde8dX9+/dPnz49PT199+7du3bt+sc//rF69WqTybRnzx5vb2/nB/Py8h5++OH//d///f777wMCAh577DHn9k2bNlkslmPHjg0dOnT27Nn1Bpadne3v779582abzfb000+np6f/+OOPPXr0ePXVV0XkQkplZWUDYwYakwO4fnTs2NHb29vHx8c5KzJ48GDn9rfffvvC8tGjRz09PSsqKvbu3dutW7efe/XPf/7z448/7ty4b9++TZs2ORwOk8nkcDi2bdt2xx13OByOefPmTZo0yfme4uLipk2b1tXV7d2719vbu7Ky0vlBZ8QFO3fuNJvN586dczgcZ86cKS4udjgcZWVl06ZNGz9+vPM9zpQrjsqI/2iAydX/0AC/zIYNG5xz9AcPHoyLi1uwYEFiYmJhYeGOHTsufPUaFhZWVVV14SNXfPXo0aOXbrli1rFjx2666Sbncps2bTw8PJwTLJ06dWrRooWIeHl5/fRT7dq18/T0FBFPT89XXnll2bJlbdu2bdGiRWBg4KVvu+KofH19/7v/LEADKHpcr0JCQkaNGrV3714RCQwMHD169KJFi0Tk7NmzW7dubdeu3YVZ7yu+GhgY+MMPPzjfsH379vz8/LFjx9aLaNu2bW5urnO5rKysrq7ObDYfO3bMw8OjgYE5W15Eli9fnpWVtWPHjtatWy9duvSzzz679G1XHNWv+08CXBlz9LiOtWrVqrKyUkSGDRuWlpa2ZcuWmpqaF154YdasWZe+7YqvjhkzZunSpXl5eWVlZdOnT7/wpe6ZM2cufHDkyJGrVq3auXPn2bNn//SnP8XGxl4o8WtRVlbWsmVLPz+/kydPLlq0qKam5sJLZ86caXjMQCOi6HEdu+2227744ovKysqgoKC33377oYceCgwM3LJlyzvvvHPp2674as+ePefMmTN48OAuXbrcdtttU6dOFZHhw4ffeuutFxr5tttuS01NHTduXPv27Y8dO/bGG2/8ouHdd999np6eQUFBw4cPf+SRR77++utPPvnkQkpgYGADYwYaEefRQ1v79u0bP378nj17XD0QwMU4ooe29u/f7/zKFHBzFD30tHPnzv/5n/+ZOXOmqwcCuB5TNwCgOY7oAUBzFD0AaI6iBwDNUfQAoDmKHgA0R9EDgOYoegDQHEUPAJqj6AFAcxQ9AGiOogcAzVH0AKA5ih4ANEfRA4DmKHoA0Nz/ByY9pP/Goh6IAAAAAElFTkSuQmCC)