DEA estima la medida de eficiencia de una DMU específica por comparación relativa de su desempeño con las otras. DEA asume que hay NDMUs (j = 1,2,...,N). El desempeño de cada DMU es caracterizada por su proceso de producción de m entradas (xlj para l = 1,...,m) para producir S salidas (yrj para r = 1,..., S).

Uno de los métodos clásicos para tratar la incertidumbre en DEA es el CCP (por sus siglas en inglés "chance-constrained programming"), es ampliamente tratado en Land et al. [12, 13, 14, 15]. El modelo (1) es un CCP extendido del CCR (Chames et al. [11]) en su forma primal.

El modelo (1) está basado en los multiplicadores u y v los cuales son variables que hacen que se encuentre la máxima , si una DMU es ineficiente no es porque se hayan asignado mal sus multiplicadores sino porque en verdad es ineficiente.

Se puede consultar en Cooper et al. [16] cómo el modelo (1) se transforma en su equivalente determinístico (2).

En modelo (2), 𝐾 _{(1−𝛼𝑗 )}=ϕ−1(1 - 𝛼_𝑗 ) y ϕ⁻¹ es la función fractile asociada con la distribución normal estándar. 𝛴_𝑗 es la matriz de covarianza y las 𝜂_𝑗 son llamadas “splitting variables”.

Sin embargo, para la aplicación del generador propuesto a DEA estocástico Bayesiano, en el presente artículo la postura parte del modelo CCP dual definido en Subhash [17]:

En el modelo (3), ϕ es la medida de eficiencia o score, los xj son los valores de la variable de entrada en la DMU j, yj son los valores de la variable de salida en la DMU j; los λj son valores que hacen que se maximice ϕ y son valores que ayudan a construir la frontera de eficiencia; xo y yo son los valores de entrada y salida respectivamente para la DMU observada; N simboliza el número de DMUs, así que este modelo se corre N veces para hallar las eficiencias relativas de todas las DMUs.

Este modelo se transforma a su equivalente determinístico a través de la creación de una nueva variable aleatoria u usando ecuación (4) la cual se distribuye normalmente con parámetros dados en ecuaciones (5) y (6). Dicha transformación se puede consultar en Subhash [17].

Continuando con la transformación y haciendo reemplazos correspondientes se llega al modelo (7).

De esta forma el modelo CCP dado en (3), que es dual, se transforma en su equivalente determinístico (9).

El gradiente y la derivada direccional son los fundamentos matemáticos y base para la propuesta del generador bivariado. Se parte del hecho de que la derivada direccional de una distribución de probabilidad acumulada F(x,y) nos indica la velocidad a la que se acumula probabilidad en un punto (x,y) en una dirección dada por un vector unitario U. Para el caso univariado, si en una función de densidad f (x) un valor x¡ tiene mayor densidad de probabilidad que un valor X2 entonces se cumple que . Es claro que si se toma un incremento Δx = x² − x¹ y se calcula la pendiente estaríamos observando la velocidad a la que se acumula probabilidad en un intervalo Δx. Por ejemplo en la distribución normal, si se parte de la media μ, y formamos un intervalo (μ,μ+ Δx) se espera más tasa de acumulación de probabilidad en este intervalo que en otro intervalo, digamos (P95, P95 + Δx). Dicho de otra manera, , siendo P95 el percentil 95. Este concepto es llevado al caso bivariado y se propone un algoritmo a partir de esta idea. Dicho algoritmo será capaz de generar números aleatorios que sigan una distribución de probabilidad bivariada. Si f es una función en dos variables x e y, y las derivadas parciales f_x y f_y existen, entonces el gradiente de f denotado ∇f está definido por la ecuación (16).

Y la derivada direccional se puede escribir como el producto punto entre un vector unitario y el vector gradiente.

donde es la derivada parcial de la función 𝑓 (𝑥,𝑦) y describe la variación de f en la dirección positiva del eje es la derivada parcial de la función f (x,y) y describe la variación de f en la dirección positiva del eje y; U = cosθi + sinθj representa el vector unitario que forma un ángulo de θ radianes con la parte positiva del eje x; y 𝐷𝑈 𝑓 (𝑥,𝑦) =𝑓_𝑥 (𝑥,𝑦)𝑐𝑜𝑠𝜃+ 𝑓_𝑦 (𝑥,𝑦)𝑠𝑖𝑛𝜃 representa la derivada direccional de f en la dirección U. Una derivada direccional de una función diferenciable puede obtenerse mediante el producto punto del gradiente y un vector unitario de la dirección deseada, ecuación (17). El valor máximo de DU f (x0, y0) ocurre cuando cosα = 1 o cuando α = 0 en la ecuación (18).

El algoritmo consiste en conseguir direcciones aleatorias volviendo a U un vector aleatorio. Podemos generar aleatoriamente puntos (x,y) con x ∈ X y y ∈ Y y calcular las derivadas direccionales máximas aplicando ‖∇F(x,y)‖ y obtener mediante procesos iterativos dicha derivada direccional máxima. Una vez obtenida DU máx se puede encontrar la derivada direccional en un punto aleatorio (x,y) con la dirección del vector aleatorio U y hacer la relación ´ la cual es menor a 1 y x e y son valores aleatorios generados por ecuaciones (19) y (20), donde los respectivos Ri son variables aleatorias uniformes en el intervalo (0,1).

En las ecuaciones (19) y (20) a y b son, respectivamente, el límite inferior y superior de la variable aleatoria X; c y d son, respectivamente, el límite inferior y superior de la variable aleatoria Y. Un vector aleatorio U será generado con ecuaciones (21), (22) y (23).

● Paso 1. Encontrar la función de distribución mediante

● Paso 2. Calcular DU F(x,y) máxima usando DU F(x,y) = ‖∇𝐹(x,y)‖y ecuaciones (19) y (20).

● Paso 3. Hacer G(x,y,R3) = DUF(x,y) = U · ∇F(x,y)

● Paso 4. Hacer

● Paso 5. En paso 4 sustituir y = c + (d − c) ∗ R5 y despejar x, donde x = H (R3, R4, R5)

● Paso 6. Generar R3, R4 y R5, sustituir en x = H (R3, R4, R5) y obtener la pareja (x, y) = (H (R3, R4, R5), c + (d − c) ∗ R5)

Al repetir este algoritmo se obtienen parejas aleatorias (x, y) que siguen una función de densidad de probabilidad f (x, y)

Supóngase que se van a generar números aleatorios que sigan la función dada por la ecuación (24).

A continuación, se presentan lo pasos resueltos después de los desarrollos algebraicos.

• Paso 1.

• Paso 2.

• Paso 3. Hacer G (x,y,R3)

• Paso 4. Hacer R₄

• Paso 5.

• Paso 6. Generar R₃, R₄ y R₅, sustituir en x = H (R₃, R₄, R₅) y obtener la pareja (x, y) = (H (R₃, R₄, R₅), c + (d − c) ∗ R₅). Con este paso se generaron 3586 parejas de números aleatorios. La forma como se distribuyeron se encuentra en tabla 1.

Tabla 1
Frecuencias observadas en la generación bivariada de la función de la ecuación (24)

Fuente: elaboración propia.

Las integrales dobles se consiguieron para revisar la frecuencia esperada a la luz de ecuación (24) y se usaron junto con las frecuencias de la taba 1 para aplicar el estadístico el cual dio un resultado de 26,296 indica que 21.7831 que comparado con Xq₀₅₁₆ los números generados se distribuyen como (24).

En esta subsección se presenta una adaptación del generador mostrado en [18]. Debido a que en muchas ocasiones es difícil obtener la integral doble en el paso 1 del algoritmo y tampoco resulta factible despejar x en forma explícita en el paso 5, se ha extendido el generador propuesto al caso discreto. Esta discretización del algoritmo propuesto se ha aplicado a una distribución normal bivariada.

Simulamos 100 vectores bivariados aleatorios que sigan la distribución X ∼ N(μ,Σ), donde μ = (3,9658,7,7132) y .

Paso 1 y 2. Estos pasos se resolvieron de la siguiente forma. Lo que se hizo fue elegir valores posibles x₁ e X₂ que pertenezcan al dominio de la función. Se eligió un paso pequeño de 0.1, con este construimos una matriz de puntos bivariados, como una malla de puntos en una cuadrícula. Mediante una estructura de programación "for anidado" se fueron encontrando los volúmenes bajo la curva estableciendo las alturas de acuerdo con la función bivariada. Los volúmenes se acumularon para encontrar una función de distribución aproximada. Los gradientes se construyeron de acuerdo con el paso = 0,1 y a la nueva función acumulada discretizada. La tabla 2 presenta en las columnas 2 y 3 x₁ y X₂ respectivamente.

Tabla 2
Discretización del algoritmo propuesto

Fuente: elaboración propia.

En dicha tabla se presentan los registros 1500 a 1510 solamente, ya que esta tiene 6222 registros. La columna 4 corresponde a las probabilidades acumuladas o los volúmenes bajo la curva acumulados.

Paso 3, 4, 5 y 6. Se continúa construyendo la tabla 2. Las columnas 5 y 6 corresponde al componente gradiente de x₁ y componente gradiente en x₂ respectivamente. La columna 7 es la norma del gradiente construido de las anteriores dos columnas. La columna 8 es la derivada direccional. La columna 9 el cociente entre derivada direccional y derivada direccional máxima. Una vez se tiene la tabla 2 lo que sigue es ordenar la tabla completa por la columna 9, de esta forma es fácil generar un número aleatorio que recorra la tabla hasta que este sea menor a algún valor de esta última columna ordenada de menor a mayor, de acuerdo con este resultado se selecciona el vector (x₁ , x₂) respectivo.

Las figuras 1 y 2 muestran las pruebas Mardia y Henze-Zirkler de bondad de ajuste normal multivariadas corridas en el Rcomander [19]. Un diagrama Q-Q fue construido en Matlab en la figura 3, basado en la ecuación (25), donde se puede mostrar que si los Yi siguen una distribución normal multivariada, ui en ecuación (26) tiene distribución Beta con α = p/2 y β = (n − p − 1)/2 con p el número de variables n las observaciones. La prueba de bondad de ajuste de ui a una Beta se ve en la figura 4 que es una salida del software Infostat.

Figura 1
Prueba de normalidad MardiaTest.

Figura 2
Prueba de normalidad Henze-Zirkler.

Figura 3
Diagrama Q-Q para los datos simulados.

Figura 4
Prueba de bondad para los u_i.

Una vez logrado que el generador bivariado tanto discreto como continuo pasen las pruebas de bondad de ajuste, se puede poner en práctica al utilizarlo en un modelo DEA estocástico Bayesiano. Para correr dicho DEA Bayesiano se han tenido en cuenta varios aspectos. Primero se usaron datos de DANE Colombia (ver siguiente subsección "descripción de los datos"). Segundo se transformó un DEA estocástico expresándolo en forma de conjuntos convexos para garantizar que las soluciones de un problema no lineal son óptimos globales. Tercero se corrió el modelo estocástico transformado utilizando la distribución a posteriori de las variables de salida. Cuarto se obtuvieron, gracias a los aspectos anteriores, los determinados interevalos de confianza, las distribuciones de probabilidad de las eficiencias y las comparaciones con DEA clásico. En las siguientes subsecciones se amplía cada uno de los anteriores aspectos.

Se han usado datos del DANE Colombia, en un estudio sobre la investigación de educación formal en Colombia. Este estudio está enfocado a los niveles de preescolar, básica primaria, básica secundaria y educación media [20]. En dicho estudio se cuenta con una base de datos llamada indicadores de eficiencia 2004-2014 (es la versión más actual encontrada que está organizada por departamentos) que contiene la información de porcentajes de alumnos aprobados, reprobados, desertores y transferidos de todos los municipios de Colombia.

Para probar el generador aplicado al DEA estocástico Bayesiano, se seleccionaron 4 variables medidas en el año 2014, dos de entrada y dos de salida. Como variables de entradas se definieron el porcentaje de estudiantes reprobados y porcentaje de estudiantes desertores. Como variables de salida el porcentaje de estudiantes aprobados y porcentaje de transferidos. Como DMUs se tuvieron en cuenta los 19 municipios de mayor tamaño de Colombia. Con estas 4 variables se desea crear una medida comparativa de resultados para los municipios. La lógica es que las salidas mientras más altas mejor, es decir mientras más estudiantes aprobados y transferidos mejor. Las entradas mientras menores mejor, así que mejor será que haya menor reprobados y menos desertores para la educación en Colombia. Se tienen en cuenta las siguientes definiciones de acuerdo con el Ministerio de Educación Nacional [21]:

• Reprobados. Son los alumnos que durante el período electivo no logran cumplir con todos los requisitos académicos para ingresar al siguiente período escolar.

• Desertores. Son los estudiantes que abandonan el sistema escolar durante el período académico, o que finalizan el período lectivo, pero no regresan en el período siguiente.

• Aprobados. Son alumnos que durante el periodo lectivo cumplen con todos los requisitos académicos para ingresar al siguiente período escolar.

• Transferidos. Son los alumnos que se encontraban matriculados en un establecimiento educativo y que cambiaron de institución educativa durante el año escolar, por lo tanto, siguen siendo cubiertos por el sistema educativo.

En [22] se puede consultar el problema de optimización cono compuesto como sigue:

donde C bebe satisfacer los siguientes requerimientos. Sea x^t ∈ 𝑅^𝑛𝑡 , t = 1,...,k vectores incluidos en las variables de decisión, entonces se define: C := { 𝑥 ∈ 𝑅^𝑛∶𝑥^𝑡 ∈ 𝐶_𝑡 ,𝑡 = 1,2,...,𝑘} o donde Ct debe tener una de las siguientes formas:

• Conjunto R:

𝐶𝑡= {𝑥 ∈ 𝑅^𝑛𝑡}

• Cono cuadrático:

• Cono cuadrático rotado:

Basado en los conceptos anteriores, se transformó el modelo (9) en un problema convexo usando un cono cuadrático y restricciones lineales para garantizar de que, a pesar de tener un problema no lineal, se puedan lograr soluciones óptimas globales. Se puede cambiar el primer conjunto de restricciones de la salida creando nuevas variables y obteniendo conos cuadráticos. Entonces el primer conjunto de restricciones de (9) queda así: que se puede expresar como:

Reemplazando el lado izquierdo por una variable x₁ como el cono cuadrático tenemos:

Si por ejemplo la DMU observada es la 19, entonces se pueden hacer los siguientes reemplazamientos:

Ahora, para la DMU 19, se puede expresar (28) así:

que es un cono cuadrático convexo. Las restricciones del modelo se completan con las siguientes restricciones lineales que parten de (29).

Las variables x1, x2,... , x20 son todas expresiones de tipo lineal. Un procedimiento muy similar se realiza para las restricciones de la segunda salida. De esta forma se demuestra que el problema es de programación no lineal con restricciones convexas.

Se programó en el software Matlab el generador bivariado discretizado propuesto en el presente artículo para simular valores que siguieran la distribución a posteriori (15) cuya marginal para μy es de Student multivariada la cual está definida por ecuaciones (32) y (33) de acuerdo con siguiente definición:

Sea Y un vector aleatorio continuo K × 1. Que μ sea un vector K × 1, V una matriz simétrica definida positiva y n ∈ ℜⅇ++. Se dice que, Y tiene distribución Student multivariada con media μ, matriz de escala V y n grados de libertad si su función de densidad de probabilidad conjunta es:

Donde:

Para las matrices de varianzas y covarianzas se usó la simulación de Wishart de Matlab. Al obtener las medias y varianzas necesarias que provienen de distribuciones a posteriori, se alimenta el modelo DEA estocástico (9).

Para evaluar la incertidumbre de los parámetros de la distribución a posteriori, se tuvo en cuenta que, como era de esperarse, la parte derecha de la expresión (15) es normal multivariada ya que la marginal de μy es Student multivariada. Se hizo la prueba de multinormalidad de Henze-Zirkler a los datos obtenidos del generador discretizado que sigue la distribución a posteriori en expresión (15) y dicha prueba fue superada de acuerdo con la figura 5.

Figura 5
Prueba de normalidad Henze-Zirkler a posteriori.

De acuerdo con la figura 5 se puede obtener que P ((84,09 < μY 1< 93,29)|Y ) = 0,90 y P ((0,27<μY2< 9,13)|Y ) = 0,90

En las subsecciones anteriores se logró obtener p(θ|Z), lo que permite llegar al cálculo de las eficiencias p(r|Z). De esta forma se ha podido llegar al cálculo de las eficiencias de los principales municipios de Colombia teniendo en cuenta las medidas del DANE Colombia. Se usaron rutinas del programa Mosek [22], específicamente el toolbox Mosekopt, para correrse en Matlab. Estas rutinas usadas son diseñanadas por Mosek para resolver un problema de optimización cónica a través de de conos convexos. Para cada municipio (DMU) se corrió el modelo estocástico Bayesiano 30 veces. La tabla 3 muestra los puntajes de DEA clásico orientado a las salidas con retorno a escala contante (CRS por sus siglas en inglés) comparado con del DEA Bayesiano estocástico. De este último se presenta también la media, desviación estándar y los intervalos de confianza del 95% construidos con el estadístico t student (usando el teorema del límite central), donde L.I significa límite inferior y L.S el límite superior. También se presentan el ranking alcanzado por DEA clásico y Bayesiano.

Tabla 3
Resultados empíricos

Fuente: elaboración propia.

Las pruebas de bondad de ajuste para tres distribuciones de probabilidad adecuadas de acuerdo con los histogramas que presentaron las eficiencias de los municipios son presentadas en la tabla 4.

Tabla 4
Pruebas de bondad de ajuste Kolmogorov

Fuente: elaboración propia.

La figura 6 muestra una de las corridas que se usó en la DMU 19 usando software MOSEK. En la figura 7 se muestra una de las corridas hechas, para la DMU 19, se puede apreciar que la solución es óptima debido a que el problema primal y dual llegaron al mismo resultado.

Figura 6
Salida Mosek.

Figura 7
Solución óptima usando software Mosek.

Las pruebas de bondad de ajuste para tres distribuciones de probabilidad adecuadas de acuerdo con los histogramas que presentaron las eficiencias de los municipios son presentadas en la tabla 4.

La figura 6 muestra una de las corridas que se usó en la DMU 19 usando software MOSEK. En la figura 7 se muestra una de las corridas hechas, para la DMU 19, se puede apreciar que la solución es óptima debido a que el problema primal y dual llegaron al mismo resultado.