Técnicas de clasificación no contexto de Big Data

Traballo Fin de Grao Técnicas de clasicación no contexto de Big Data Josefa Arán Paredes 2018/2019 UNIVERSIDADE DE SANTIAGO DE COMPOSTELA GRAO DE MATEMÁTICAS Traballo Fin de Grao Técnicas de clasicación no contexto de Big Data Josefa Arán Paredes Xullo, 2019 UNIVERSIDADE DE SANTIAGO DE COMPOSTELA Traballo proposto Área de Coñecemento: Estatística Título: Técnicas de clasicación no contexto de Big Data Breve descrición do contido Trátase de estudar as técnicas de clasicación ou análise discriminante máis importantes e de referencia, e o papel que xogan ditas técnicas no moderno contexto do Big Data. Recomendacións Guión aproximado: 1) Técnicas de clasicación lineal ou cuadrática. 2) Técnicas de clasicación non paramétrica. 3) Algunhas técnicas de clasicación adaptadas ao contexto de alta dimensión en Big Data. 4) Ilustración en bases de datos reais. Outras observacións Preténdese que o alumno dedique aproximadamente tres meses ao estudo metodolóxico das técnicas correspondentes aos puntos 1), 2) e 3) do guión. Un mes para a avaliación do software dispoñible en R e outro para o desenvolvemento da aplicación con datos simulados ou reais. iii Índice xeral Resumo viii Introdución xi 1. Introdución á análise discriminante 1 1.1. Clases, etiquetas, regras e funcións de decisión. . . . . . . . . . . . . . . . . 1 1.1.1. Fronteiras e rexións discriminantes . . . . . . . . . . . . . . . . . . . 3 1.2. Avaliación das regras e probabilidade de clasicación errónea . . . . . . . . . 4 1.2.1. O problema de Bayes . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 1.2.2. AregradeBayes............................. 7 1.2.3. Pérdida e risco de Bayes . . . . . . . . . . . . . . . . . . . . . . . . . 9 1.3. Análise discriminante no caso mostral . . . . . . . . . . . . . . . . . . . . . 12 2. Técnicas de clasicación 15 2.1. Regras discriminantes lineais . . . . . . . . . . . . . . . . . . . . . . . . . . 15 2.1.1. A regra discriminante de Fisher . . . . . . . . . . . . . . . . . . . . . 15 2.2. Clasicación baixo hipótese de normalidade . . . . . . . . . . . . . . . . . . 21 2.2.1. Dúas clases normais univariantes coa mesma varianza . . . . . . . . 21 2.2.2. Dúas ou máis clases normais multivariantes coa mesma matriz de covarianzas ................................ 22 2.2.3. Regras discriminantes cuadráticas . . . . . . . . . . . . . . . . . . . . 25 2.3. Técnicas de clasicación non paramétrica . . . . . . . . . . . . . . . . . . . 29 2.3.1. Regra dos k veciños máis cercanos . . . . . . . . . . . . . . . . . . . 30 2.3.2. Regras tipo kernel . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 3. Estimación do erro e avaliación de regras discriminantes 37 4. Consideracións sobre a alta dimensión e Big Data 45 4.1. Reducción da dimensión e regularización . . . . . . . . . . . . . . . . . . . . 48 v vi ÍNDICE XERAL 4.1.1. Análise discriminante de Fisher . . . . . . . . . . . . . . . . . . . . . 48 4.1.2. Análise discriminante regularizado . . . . . . . . . . . . . . . . . . . 49 4.2. Análise dicriminante linear regularizado disperso: sparse rLDA . . . . . . . 50 4.3. Análise discriminante en alta dimensión: HDDA . . . . . . . . . . . . . . . . 50 4.3.1. Estimación dos parámetros . . . . . . . . . . . . . . . . . . . . . . . 53 5. Ilustración sobre datos simulados e reais 55 5.1. Datossimulados ................................. 55 5.2. Datos de medidas de expresión xénica . . . . . . . . . . . . . . . . . . . . . 58 A. Scripts utilizados para os exemplos e implementación do HDDA 61 A.1.Exemplos ..................................... 61 A.2. Implementación do HDDA e aplicación a datos simulados e reais . . . . . . 68 Bibliografía 73 2 CAPÍTULO 1. INTRODUCIÓN Á ANÁLISE DISCRIMINANTE Denición 1.1. Consideramos as clases C1, ..., Cκ , determinadas polas funcións de desidade f1, ... , fκ . Dicimos que un vector aleatorio d -dimensional X pertence á clase Cν para algún ν≤κ , se X ten como densidade de probabilidade fν , é dicir, satisfai as propiedades que caracterizan Cν . Esta pertenza denotámola como X ∈ Cν ou X [ν]. Para vectores aleatorios X pertencentes a unha das κ clases, a etiqueta Y de X é unha variable aleatoria que toma valores discretos 1,..., κ , tal que Y=ν se X ∈ Cν. Consideramos " X Y# como un vector aleatorio (d+ 1) -dimensional e chamámolo vector aleatorio etiquetado. Ás clases ás veces chámaselles poboacións , e usaremos ambos termos indistintamente. Se as κ clases están caracterizadas polas súas medias µν e as súas matrices de covarianzas Σν de xeito que Cν≡(µν,Σν) , escribiremos X ∼(µν,Σν) ou X ∈ Cν. De non dicir o contrario, os vectores aleatorios X pertencentes ás clases C1, ..., Cκ teñen a mesma dimensión d . Posto que se os vectores de distintas clases tivesen dimensións diferentes ou medisen distintas variables, esta información podería empregarse para clasicalos e non precisaríamos as técnicas máis sosticadas da análise discriminante. Unha regra é o mecanismo que nos permite asignar un vector aleatorio a unha clase, matemáticamente denímola da seguinte maneira: Denición 1.2. Sexan κ clases C1, ..., Cκ e X un vector aleatorio que pertence a unha delas. Unha regra discriminante ou clasicador para X é unha aplicación r que asigna a X un número l∈ {1, ... , κ} . Escribimos r( X ) = l con 1≤l≤κ. A regra r asigna X á clase correcta ou clasica X correcamente se r( X ) = ν cando X ∈ Cν e clasica X incorrectamente noutro caso. 1.1. CLASES, ETIQUETAS, REGRAS E FUNCIÓNS DE DECISIÓN. 3 Son de especial interese as situacións nas que hai dúas clases, nas cales unha función de decisión dá unha expresión alternativa á da regra en cuestión. Denición 1.3. Sexan X un vector aleatorio que pertence a unha das dúas clases C1 ou C2 e r unha regra de discriminante para X . Unha función de decisión para X , asociada a r , é unha función real h denida do seguinte xeito h( X )>0 se r( X ) = 1 h( X )<0 se r( X ) = 2. Observación 1.4 . Unha función de decisión correspondente a unha regra non é única. Por exemplo, calquera múltiplo dunha función de decisión por un escalar positivo tamén é unha función de decisión para a mesma regra discriminante. 1.1.1. Fronteiras e rexións discriminantes Posto que unha observación non é máis que un vector aleatorio d -dimensional, construir unha regra pódese ver como denir unha partición en rexións G1, ... , Gκ disxuntas de Rd , intentando que os vectores que caian na rexión Gν sexan da clase Cν . Para un problema de dúas clases, a decisión da pertenza a unha clase pódese basar nunha función de decisión que separa en rexións. As fronteiras entre as rexións son interesantes como axuda visual cando temos dúas clases e dimensión d pequena. Denición 1.5. Sexa X un vector aleatorio pertencente a C1 ou C2 . Se r é unha regra discriminante para X e h é unha función de decisión asociada, denimos a fronteira de decisión B da regra r como o conxunto formado por tódolos vectores aleatorios X tales que h( X ) = 0 . Cando a regra é linear veremos que a fronteira de decisión é un hiperplano, posto que h é unha función linear. Para regras non lineais, téñense fronteiras máis complexas. Pódese tamén estender a denición de fronteira de decisión para máis de dúas clases, pero estas complícanse ao medrar o número de clases e deixan de ser tan útiles. Tamén podemos denir unha regra de decisión a partir da fronteira, diferenciando os vectores aleatorios que quedan a un e outro lado de B . Denición 1.6. Se temos κ clases e X un vector aletorio, para ν≤κ , a rexión discriminante Gν da regra r defínese como Gν={ X :r( X ) = ν}. 4 CAPÍTULO 1. INTRODUCIÓN Á ANÁLISE DISCRIMINANTE As rexións discriminantes Gν son disxuntas, xa que cada X se asigna a unha única clase. Podemos interpretalas como as clases denidas pola regra, isto é, as determinadas por r . Recíprocamente, dadas as rexións disxuntas, podemos denir unha regra de xeito que r( X ) = ν se X ∈Gν . Polo tanto os conceptos de regra e de rexións discriminantes defínense un ao outro, e a fronteira separa as rexións Gν . Se a regra fose perfecta, entón tódolos puntos na rexión Gν serían observacións da clase Cν e viceversa. 1.2. Avaliación das regras e probabilidade de clasicación errónea Nun caso ideal, a nosa regra asignaría a X a súa etiqueta Y . Pero isto non sempre ocorre, polo que precisamos criterios que axuden a medir a calidade das nosas regras. Para iso, é necesario desenvolver un marco teórico contra o cal contrastar as regras, xa que se só nos xásemos nos datos dos que dispoñemos poderiamos atopar unha regra que clasicase á perfección as observacións dispoñibles pero fallase ao empregala sobre novos datos. Cando dispoñemos de máis dunha regra, é importante entender cal delas se comporta mellor, e baixo que condicións. Non hai unha única medida do comportamento dunha regra nin unha regra universal que funcione mellor en tódalas situacións. Cantas máis clases temos, máis erros pode cometer unha regra á hora de clasicar, como se pode ver na seguinte táboa. Nela indícase o que pode ocorrer ao clasicar unha observación en κ clases. Cadro 1.1: Etiquetas, regras e erros de clasicación Asignación da regra 1 2 ··· κ 1 (1,1) (1,2) ··· (1, κ ) Valor da 2 (2,1) (2,2) ··· (2, κ etiqueta . . .. . .. . ..... . . κ ( κ ,1) ( κ ,2) ··· ( κ , κ ) Un dato só estaría ben clasicado cando a asignación pola regra coincide coa súa etiqueta, que serían os casos da diagonal. Así, a probabilidade de acertar de maneira azarosa 1.2. AVALIACIÓN DAS REGRAS E PROBABILIDADE DE CLASIFICACIÓN ERRÓNEA 5 sería 1 κ . No caso de ter dúas clases, asignar ao azar ten un 50% de probabilidades de acerto, pero clasicar correctamente deste xeito vólvese máis difícil cantas máis clases hai, pois a probabilidade de acerto diminúe ao medrar κ . A probabilidade de cometer un erro sería entón κ−1 κ , que tende ao 100% a medida que aumenta o número de clases κ . Cando construamos unha regra, sempre quereremos que a proporción de observacións mal clasicadas fronte ao total sexa menor que κ−1 κ , pois de non ser así, daríanos peores resultados que asignando cada observación a unha clase ao azar. Non só iso, debemos averiguar como de boa pode chegar a ser a nosa regra, ou cal será a mellor regra entre unha clase de regras especíca como poden ser as lineais. Haberá problemas para os que esa regra óptima nunca sexa mellor que decidir ao azar, é dicir, cuxo erro mínimo sexa maior que κ−1 κ para κ clases. Os criterios que permiten determinar cal é a regra óptima deben medirse tanto en termos da distribución das variables aleatorias como da mostra, co cal é preciso achar maneiras de estimar o erro cometido. Unha maneira natural de cuanticar o bo comportamento dunha regra sería contar o número de clasicacións erróneas fronte ao total de observacións que intentamos clasicar, pero serán necesarias medidas máis complexas para medir se unha regra se comporta como é desexable, que introduciremos ao longo desta sección. 1.2.1. O problema de Bayes Comezamos introducindo os conceptos para o caso de dúas clases, aínda que se poden estender facilmente a problemas con máis clases. Denición 1.7. Sexa X un vector aleatorio que pertence a unha das clases C1 ou C2 e sexa Y a etiqueta de X . Consideramos r unha regra para X , e denimos a probabilidade de clasicación errónea ou probabildade de erro de r como L(r) = P{r( X )6=Y}. Quereremos entón atopar regras para as que a probabilidade de erro sexa o máis pequena posible. Denimos a regra que minimiza a probabilidade de erro de clasicación entre tódalas posibles regras como segue: r∗( X ) = arg min r:Rd→{1,...,κ} P(r( X )6=Y). Esta regra depende da distribución do vector etiquetado " X Y# , e se a coñecemos, poderemos calcular a regra r∗ explícitamente. O máis habitual é que descoñezamos a súa distribución, e con ela tamén a regra r∗ . Ao problema de achar esta regra chamarémoslle 6 CAPÍTULO 1. INTRODUCIÓN Á ANÁLISE DISCRIMINANTE problema de Bayes . No seguinte resultado veremos unha expresión alternativa para r∗ , denominada regra de Bayes . Proposición 1.8. Sexa X un vector aleatorio pertencente a unha das clases C1 ou C2 coa súa etiqueta Y . Se p( X ) = P{Y= 1| X } é a probabilidade condicional de Y= 1 dado X e denimos a regra discriminante r∗ como r∗( X ) = (1 se p( X )>1/2 2 noutro caso . Dada r outra regra discriminante para X , entón P{r∗( X )6=Y} ≤ P{r( X )6=Y} Demostración. Considerada a regra r para X , tense P{r( X )6=Y| X }= 1 −P{r( X ) = Y| X } = 1 −(P{Y= 1,r( X )=1| X }+P{Y= 2,r( X )=2| X }) = 1 −I{r( X )=1}P{Y= 1| X }+I{r( X )=2}P{Y= 2| X } = 1 −(I{r( X )=1}p( X ) + I{r( X )=2}(1 −p( X )) onde IG é a función indicador do conxunto G. Empregando os cálculos anteriores, obtemos δ∗:= P{r∗( X )6=Y}−P{r( X )6=Y} =p( X )(I{r∗( X )=1}−I{r( X )=1}) + (1 −p( X ))(I{r∗( X )=2}−I{r( X )=2}) = (2p( X )−1)(I{r∗( X )=1}−I{r( X )=1}) posto que I{r( X )=2}= 1 −I{r( X )=1} . Finalmente, pola denición de r∗ , chegamos a que δ∗≥0 . Esta proposición proba a optimalidade da regra r∗ . A probabilidade de erro desta regra é o que denominamos erro de Bayes L∗=L(r∗) = P(r∗( X )6=Y) e é o mínimo erro que podemos esperar cometer ao clasicar. Podemos calculalo cando a distribución do vector etiquetado é coñecida. 1.2. AVALIACIÓN DAS REGRAS E PROBABILIDADE DE CLASIFICACIÓN ERRÓNEA 7 Observación 1.9 . Como temos dúas clases, dada X , P(Y= 2| X )=1−p( X ) , polo tanto P(Y= 1| X )>P(Y= 2| X )⇔P(Y= 1| X )−P(Y= 2| X )>0 ⇔p( X )−(1 −p( X )) = 2p( X )−1>0⇔p( X )>1/2. Co cal podemos reescribir r∗ como r∗( X ) = (1 se P(Y= 1| X )>P(Y= 2| X ) 2 noutro caso . É dicir, a regra de Bayes asigna unha observación X á clase á que é máis probable que pertenza, como podía suxerir a nosa intuición. 1.2.2. A regra de Bayes Introducimos un entorno probabilístico de xeito que " X Y# é un vector aleatorio etiquetado que toma valores en Rd× {1, ... , κ} , cunha certa distribución e que representa a probabilidade de atoparnos certos pares vector-etiqueta na práctica. Cando Y=ν , a distribución de X virá dada por fν . Para avaliar as regras dende un punto de vista teórico é importante denir o concepto de probabilidade do erro, no que se fai ncapé en Devroye et al. (1996), o cal quereremos acotar co obxectivo de ter un indicador da propensión dos nosos datos a ser clasicados correctamente. Nun contexto bayesiano con clases C1, ..., Cκ supoñemos que coñecemos a probabilidade de que unha observación pertenza á clase Cl e denotámola πl , é dicir, πl=P( X ∈ Cl) . As probabilidades π1, ... , πκ chámanse probabilidades a priori . Se dispoñemos dunha mostra, tomamos como probabilidades a priori a proporción de observacións que pertencen a cada clase. Incluir as probabilidades a priori pode enriquecer a clasicación cando son substancialmente distintas entre as clases. Por exemplo, se temos dúas clases e sabemos que dous tercios dos datos proveñen da primeira clase e o tercio restante da segunda, unha regla debería asignar máis probablemente unha observación á clase 1. Se dispoñemos desta información, incorporada de xeito coidadoso no problema de decisión pode mellorar as nosas regras. Denición 1.10. Sexan C1, ..., Cκ clases que se diferencian nas súas medias e matrices de covarianzas. Sexa X un vector aleatorio pertencente a unha das clases. Consideramos as probabilidades a priori π1, ... , πκ asociadas ás clases. Sexa r unha regra discriminante derivada das rexións Gν . 8 CAPÍTULO 1. INTRODUCIÓN Á ANÁLISE DISCRIMINANTE 1. A probabilidade condicional de que r asigne X á clase Cν , cando X é da clase Cl é p(ν|l) = P{r( X ) = ν| X ∈ Cl}. 2. A probabilidade a posteriori de que unha observación X pertenza á clase Cl cando a regra lle asigna o valor ν é P{ X ∈ Cl|r( X ) = ν}=P{r( X ) = ν| X ∈ Cl}πl P{r( X ) = ν}=p(ν|l)πl P{r( X ) = ν}. 3. A probabilidade de clasicación errónea ou probabildade de erro de r defínese como L(r) = P{r( X )6=Y}, e usando as probabilidades condicionais e a priori ten a forma L(r) = X ν6=l p(ν|l)πl. Posto que as rexións discriminantes e as regras discriminantes se determinan unhas a outras, a probabilidade de asignar X a Gν cando X ∈ Cl é p(ν|l) = ZGν fl( x )d x =ZIGνfl, sendo IGν a función característica da rexión ν -ésima e fl a función de densidade da clase l -ésima. Tense entón que a probabilidade de acertar ao asignar X á clase Cν é P{ X ∈ Cν|r( X ) = ν}=p(ν|ν)πν P{r( X ) = ν}, mentres que a probabilidade de clasicar correctamente unha observación da clase ν -ésima é p(ν|ν) e a de cometer un erro é 1−p(ν|ν) . O seguinte Teorema dene a regra de Bayes e da a súa expresión. Teorema 1.11. Sexan C1, ..., Cκ clases con distintas medias e π= [π1, ... , πκ] as probabilidades a priori asociadas. Sexa X un vector aleatorio dunha das κ clases e fν a función de densidade da clase Cν . Denimos as rexións Gν={ X :fν( X )πν= m´ax 1≤l≤κ[fl( X )πl]}, e consideramos rBayes a regra denida a partir das Gν de tal xeito que rBayes( X ) = ν se X ∈Gν . A regra rBayes asigna X á clase Cν en preferencia a Cl cando fν( X ) fl( X )>πl πν . 1.2. AVALIACIÓN DAS REGRAS E PROBABILIDADE DE CLASIFICACIÓN ERRÓNEA 9 Demostración. A proba séguese directamente da denición da regra rBayes , xa que se ten fν( X )πν= m´ax 1≤l≤κ[fl( X )πl]⇔fν( X )πν≥fl( X )πl∀l⇔fν( X ) fl( X )≥πl πν∀l. Chamamos a rBayes regra (discriminante) de Bayes . Empregando esta regra, a probabilidade de asignar X á clase correcta é p= κ X ν=1 P{rBayes( X ) = ν| X ∈ Cν}πν= κ X ν=1 ZIGνfνπν, e a función característica IGν cumpre que IGν( X ) = (1 se fν( X )πν≥fl( X )πl∀l 0 noutro caso . Se só tiveramos dúas clases, esta sería a regra r∗ denida na sección anterior. 1.2.3. Pérdida e risco de Bayes Cando traballamos nun contexto bayesiano, a pérdida e o risco adoitan usarse para avaliar o comportamento dun método e tamén comparalo con outros diferentes. Estas ideas pertencen ao campo da teoría de decisión e nos adaptarémolas ás regras discriminantes. Á hora de clasicar, podemos acertar, obter un resultado 'non totalmente correcto' ou chegar a unha clasicación moi incorrecta. O grao de acerto está explicado por unha función de pérdida K que asigna un coste ou pérdida a unha decisión incorrecta. Denición 1.12. Sexan C unha colección de clases C1, ..., Cκ e X un vector aleatorio dunha destas clases. Sexa r unha regra discriminante para X . 1. Unha función de pérdida K é unha aplicación que leva X e a regra r nun número non negativo chamado pérdida ou coste . Se X ∈ Cl e r( X ) = ν , entón a pérdida cl,ν producida por tomar a decisión ν cando a clase real era l é K( X ,r) = cl,ν con (cl,ν = 0 se l=ν, cl,ν >0 noutro caso . 2. A función de risco R é a pérdida esperada producida ao usar a regra r . Escribimos R(ν, r) = E[K( X ,r)], onde a esperanza se toma respecto á distribución de X , dada por fν . 10 CAPÍTULO 1. INTRODUCIÓN Á ANÁLISE DISCRIMINANTE 3. Sexan π= [π1, π2, ... , πκ] as probabilidades a priori para as clases de C . O risco de Bayes B dunha regra discriminante r respecto ás probabilidades a priori π é B(π,r) = Eπ[R(C,r)], onde a media se toma respecto π . O risco considera a pérdida en tódalas clases para unha regra especíca. Agora temos unha ferramenta para escoller entre dúas regras, eliximos a que teña o menor risco. É difícil atopar unha regra que funciona ben para tódalas clases e a información adicional proporcionada polas probabilidades a priori pode facilitar a toma de decisións. Observación 1.13 . O criterio que nós estamos tomando para decidir cal é a mellor regra, aquela que minimiza o risco de Bayes, é o criterio de Bayes. Se temos dúas clases, tense que R(l, X ) = E[K( X ,r)] = cl,1P(r( X ) = 1) + cl,2P(r( X ),2) e B(π,r) = π1R(1,r) + π2R(2,r). Polo que o criterio de Bayes supón minimizar m´ın r{π1R(1,r) + π2R(2,r)}. Este criterio non é único, podemos ter unha postura pesimista e intentar minimizar o risco máximo que poida cometerse nalgunha das clases. O criterio minimax para dúas clases consistiría en m´ın r{m´ax{R(1,r),R(2,r)}}. O tipo máis común de pérdida é a pérdida cero-un , que toma cl,ν = 0 cando l=ν e cl,ν = 1 cando non coinciden. Estos costes din soamente cando non asignamos o vector á clase á que realmente pertence. En xeral, se queremos graduar como de incorrecta foi a clasicación, os costes cl,ν e cν,l poderían ser distintos. Nós empregaremos a pérdida cero- un, respecto á cal a regra de Bayes ten unha interpretación en función das probabilidades a posteriori. Se X é da clase l -ésima, K( X ,r) = (0 se fl( X )πl≥fν( X )πν∀ν≤κ, 1 noutro caso . Co seguinte Teorema relacionaremos o erro de Bayes e a regra de Bayes denidas na sección anterior coas ideas de pérdida e risco de Bayes. 1.2. AVALIACIÓN DAS REGRAS E PROBABILIDADE DE CLASIFICACIÓN ERRÓNEA 11 Teorema 1.14. Sexa C unha colección de clases C1, ..., Cκ e π= [π1, π2, ... , πκ] o vector de probabilidades a priori asociadas a C . Sexa X un vector aleatorio pertencente a unha das clases Cν , rBayes a regra de Bayes e Gν as rexións discriminantes dadas por Gν={ X :fν( X )πν= m´ax 1≤l≤κ[fl( X )πl]}. Para a pérdida cero-un, a regra de Bayes é óptima entre todas as regras discriminantes no sentido de que ten: 1. a maior probabilidade de asignar X á clase correcta, e 2. o menor risco de Bayes para a función de pérdida cero-un. Demostración. Pola denición da función característica de Gν , tense que a función de pérdida vale cero K( X ,r)=0 ⇔ X ∈Gν⇔IGν( X )=1 . Sexa rBayes a regra de Bayes e supoñamos que existe outra regra r0 baseada na mesma función de pérdida e que a súa probabilidade de asignar X á clase correcta é maior que a de rBayes . Denotamos por p0(ν|ν) a probabilidade de que r0 clasique correctamente unha observación da clase Cν , e G0 ν son as rexións discriminantes asociadas a r0 . Se p0 é a probabilidade de clasicar correctamente X usando r0 , entón p0= κ X ν=1 p0(ν|ν)πν= κ X ν=1 ZIG0 νfνπν ≤ κ X ν=1 ZIG0 νm´ax ν{fνπν}= κ X ν=1 ZG0 ν m´ax ν{fνπν} =Zm´ax ν{fνπν}= κ X ν=1 ZGν m´ax ν{fνπν}= κ X ν=1 ZIGνfνπν = κ X ν=1 p(ν|ν)πν=p. Este cálculo contradí que r0 leve a unha probabilidade de asignar correctamente mellor, é dicir, que p0>p . Polo que rBayes é óptima. Para a segunda parte do Teorema, partimos de que, B(π,r) = κ X ν=1 πνR(ν, r) = κ X ν=1 πνZK( x ,r)fν( x )d x = = κ X ν=1 πνZRd−Gν fν= κ X ν=1 πνZRd fν−ZGν fν= 18 CAPÍTULO 2. TÉCNICAS DE CLASIFICACIÓN Así, a regra de Fisher asigna a X o número l se η|µl é a media escalar máis cercana ao escalar η| X . Empregamos as cantidades escalares por simplicidade, en vez de buscar a media µl que se achegue máis a X . Ademáis, co uso de η destámoslle dando máis peso a variables importantes de X , e reducimos o efecto daquelas variables que non contribúen moito a W−1B . Para o caso de dúas clases C1 e C2 , podemos deducir unha función de decisión h para a regra linear de Fisher rF . |η| X −η|µ1|<|η| X −η|µ2| ⇐⇒ (η| X −η|µ1)|(η| X −η|µ1)<(η| X −η|µ2)|(η| X −η|µ2) ⇐⇒ X |ηη| X −2 X |ηη|µ1+µ| 1ηη|µ1< X |ηη| X −2 X |ηη|µ2+µ| 2ηη|µ2 ⇐⇒ 2 X |ηη|(µ1−µ2)>µ| 1ηη|µ1−µ| 2ηη|µ2= (µ1+µ2)|ηη|(µ1−µ2) ⇐⇒  X −1 2(µ1+µ2)| ηη|(µ1−µ2)>0 Como supoñemos η|(µ1−µ2)>0 e non depende de X , podemos denir a función de decisión h como segue: h( X ) =  X −1 2(µ1+µ2)| η. Obtendo así que h( X )>0 se X ∈ C1 e h( X )<0 cando X ∈ C2 . Cando non se coñecen os parámetros poboacionais, empregando os datos modicaremos a regra discriminante linear de Fisher presentada para a poboación cambiando o vector de medias µν e a matriz de covarianzas Σν polas súas cantidades mostrais. Para ν≤κ , a media mostral da clase ν -ésima é Xν=1 nν nν X i=1 X [ν] i e a media das medias mostrais das clases é X =1 κ κ X ν=1 Xν. Nótese que X non fai unha media ponderada das medias de cada clase dependendo do número de observacións dentro delas. Cando tódalas clases teñen o mesmo número de observacións, X será a media mostral habitual, resultante de considerar tódolos datos xuntos. Na denición da regra, cambiaremos ¯ µ por X e cada µν por Xν . 2.1. REGRAS DISCRIMINANTES LINEAIS 19 A matriz de covarianzas mostral, empregada en lugar de Σν para cada clase vén dada por Sν=1 nν−1 nν X i=1 ( X [ν] i−Xν)( X [ν] i−Xν)|. En xeral, as regras de Fisher para a poboación e para a mostra non teñen por qué coincidir. Mostraremos no seguinte exemplo, con datos simulados que proveñen de distribucións normais, como estas regras lineais poden levarnos a resultados distintos. Exemplo 2.4. Consideramos un problema de dúas clases cuns datos simulados bidimensionais. Collemos as medias e a matriz de covarianzas seguintes µ1="0 2#,µ2="1 1# e Σ = "0,5 0 0 0,5#. Simularemos unha mostra de tamaño n= 50 , con 20 observacións da primeira clase e 30 da segunda, ambas gaussianas cos parámetros elixidos. Construimos a regra de Fisher empregando os parámetros poboacionais a partir da función de decisión h( X ) =  X −1 2(µ1+µ2)|η , del tal modo que rF( X )=1 se h( X )>0 . Ademáis, facendo isto obteremos a fronteira de decisión dada polos X tales que h( X )=0 . O autovector que precisamos é η="−√2/2 √2/2# e como un múltiplo de h por un número positivo segue sendo función de decisión para a mesma regra, podemos calcular h( X ) = ( X −1 2 "0 2#+"1 1#!)|"−1 1# =( X −"1/2 3/2#)|"−1 1#=X1−1 2X2−3 2"−1 1# =1 2−X1+X2−3 2=X2−X1−1 . A fronteira virá dada entón polos vectores ("X1 X2#:X2−X1= 1) . Por outro lado, calcularemos a regra no caso mostral, estimando o vector de medias por Xν e a matriz de covarianzas por Sν para ν= 1,2 , que nos dan X1="0,1945 1,7168#,X2="0,8979 0,7606#, S1="0,2382 0,0928 0,0928 0,6109# e S2="0,4808 −0,0263 −0,0263 0,5579 #. 20 CAPÍTULO 2. TÉCNICAS DE CLASIFICACIÓN Figura 2.1: Datos simulados de dúas clases normais e as fronteiras de decisión para a regra de Fisher poblacional e mostral ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ●● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● −1.0 0.0 0.5 1.0 1.5 2.0 −1 0 1 2 3 −1.0 0.0 0.5 1.0 1.5 2.0 −1 0 1 2 3 −1.0 0.0 0.5 1.0 1.5 2.0 −1 0 1 2 3 Observemos a Figura 2.1. Na gráca da esquerda móstranse os datos da primera clase en negro e os da segunda en vermello. A liña azul que cruza o gráco é a fronteira entre as rexións discriminantes, para a regra de Fisher empregando os datos poboacionais, mentres que a verde é a fronteira asociada á regra con datos calculados a partir da mostra simulada. Estas dúas fronteiras non coinciden, co cal tampouco han coincidir as regras que denen. Nas grácas do centro e da dereita móstranse con cruces os datos coloreados segundo a clase á que foron asignados, pola regra de Fisher para a poboación e para a mostra respectivamente, e aqueles que foron clasicados erróneamente están rodeados na cor da clase á que realmente pertencen. Vemos entón que ademáis de que as funcións de decisión para estas dúas regras non coinciden, tampouco clasican os datos do mesmo xeito. Podemos ademáis contabilizar os erros cometidos cunha das medidas introducidas para a avaliación das regras. O erro de clasicación é dun 20% para as dúas regras xa que se clasicaron mal 10 observacións, ainda que as imaxes mostran como a clasicación que efectúan é distinta. A primeira regra clasica erróneamente 3 observacións da primeira clase e 7 da segunda, mentres que a regra mostral falla en 7 datos da clase 1 e en 3 da clase 2. Hai 3 datos de cada clase que clasican mal ambas regras pois están na rexión discriminante contraria, o cal se debe a que hai certa superposición das clases. O resto de datos mal clasicados están moi preto das rexións de decisión, de feito, entre as dúas rectas, co que cada regra os asigna a un grupo distinto. 2.2. CLASIFICACIÓN BAIXO HIPÓTESE DE NORMALIDADE 21 2.2. Clasicación baixo hipótese de normalidade Canto maior coñecemento teñamos dos datos, mellores decisións poderemos tomar. Se sabemos que os datos pertencen a unha poboación normal, deberíamos incorporar este coñecemento no noso proceso de decisión. 2.2.1. Dúas clases normais univariantes coa mesma varianza Comezamos co caso máis sinxelo, cando sabemos que as nosas observacións unidimensionais proveñen de dúas clases de distribucións normais coa mesma covarianza pero distintas medias. É dicir, C1=N(µ1, σ2) e C2=N(µ2, σ2) , e supoñemos que µ1> µ2 . Para construir a nosa regra, basearémonos en que se a observación X pertence a unha das clases, asignarémola á primeira clase cando sexa máis verosímil que proveña de C1 . Deste xeito, a regra que denimos será a regra de Bayes, que minimiza a probabilidade de erro de clasicación. Para isto, empregaremos a función de verosimilitude dunha distribución normal N(µ, σ2) , que non é máis que a súa función de densidade f , dada por f(X) = 1 √2πσ2exp −1 2 (X−µ)2 σ2, que depende só do parámetro µ , pois a varianza non cambia entre as clases. A regra discriminante denirase entón como r(X)=1 se f1(X)> f2(X). Para chegar a unha expresión explícita para esta regra, desenvolvemos unha serie de desigualdades equivalentes partindo da fórmula da verosimilitude para unha observación X . f1(X)> f2(X) ⇔1 √2πσ2exp −1 2 (X−µ1)2 σ2>1 √2πσ2exp −1 2 (X−µ2)2 σ2 ⇔exp −1 2 (X−µ1)2 σ2>exp −1 2 (X−µ2)2 σ2 ⇔ − 1 2σ2(X−µ1)2>−1 2σ2(X−µ2)2 ⇔(X−µ1)2<(X−µ2)2 ⇔X2−2Xµ1+µ2 1< X2−2Xµ2+µ2 2 ⇔2X(µ1−µ2)> µ2 1−µ2 2= (µ1+µ2)(µ1−µ2) ⇔2X > µ1+µ2 22 CAPÍTULO 2. TÉCNICAS DE CLASIFICACIÓN ⇔X > µ1+µ2 2. Nas últimas liñas, usamos a hipótese de que µ1> µ2 . O resultado nal ten sentido, pois se a media da primeira clase é maior que a da segunda, cando a nosa observación excede a media de µ1 e µ2 será que está máis cerca de µ1 e será máis probable que pertenza a C1 . Esta regra pode estenderse de maneira fácil e natural ao caso d -dimensional, considerando a función de verosimilitude para unha distribución normal multivariante. 2.2.2. Dúas ou máis clases normais multivariantes coa mesma matriz de covarianzas Se é sabido que os vectores aleatorios nun problema con dúas clases veñen dunha distribución normal que comparte a matriz de covarianzas, podemos presentar o seguinte Teorema: Teorema 2.5. Sexa X un vector aleatorio gaussiano que pertence a unha das clases Cν= N(µν,Σ) , con ν= 1,2 e supoñamos que µ16=µ2 . A función de densidade da clase ν -ésima vén dada por f( X ) = (2π)−d/2det(Σ)−1/2exp −1 2( X −µν)|Σ−1( X −µν), onde det(Σ) é o determinante da matriz de covarianzas común. Sexa rnorm a regra que asigna X á clase C1 se f1( X )> f2( X ) . A función h denida por h( X ) =  X −1 2(µ1+µ2)| Σ−1(µ1−µ2) será entón unha función de decisión para a regra rnorm , e h( X )>0 se e só se f1( X )> f2( X ) . Chamamos a rnorm a regra discriminante (linear) normal , basada en clases normais coa mesma matriz de covarianzas. Normalmente, regra e función de decisión considéranse unha mesma cousa. Demostración. Cunha serie de desigualdades equivalentes chegaremos á expresión de h , partindo de que rnorm( X )=1 se f1( X )> f2( X ) . Empregaremos tamén que Σ é unha matriz simétrica. f1( X )> f2( X ) ⇔(2π)−d/2det(Σ)−1/2exp −1 2( X −µ1)|Σ−1( X −µ1)> >(2π)−d/2det(Σ)−1/2exp −1 2( X −µ2)|Σ−1( X −µ2) 2.2. CLASIFICACIÓN BAIXO HIPÓTESE DE NORMALIDADE 23 ⇔exp −1 2( X −µ1)|Σ−1( X −µ1)>exp −1 2( X −µ2)|Σ−1( X −µ2) ⇔ −1 2( X −µ1)|Σ−1( X −µ1)>−1 2( X −µ2)|Σ−1( X −µ2) ⇔( X −µ1)|Σ−1( X −µ1)<( X −µ2)|Σ−1( X −µ2) ⇔ X |Σ−1 X −µ| 1Σ−1 X − X |Σ−1µ1+µ| 1Σ−1µ1< X |Σ−1 X −µ| 1Σ−1 X − X |Σ−1µ2+µ| 2Σ−1µ2 ⇔ −2 X |Σ−1µ1+µ| 1Σ−1µ1<−2 X |Σ−1µ2+µ| 2Σ−1µ2 ⇔2 X |Σ−1(µ1−µ2) + µ| 2Σ−1µ2−µ1Σ−1µ1>0 ⇔2 X |Σ−1(µ1−µ2)−µ| 1Σ−1(µ1−µ2)−µ| 1Σ−1µ2−µ| 2Σ−1(µ1−µ2) + µ| 2Σ−1µ1>0 ⇔2 X |Σ−1(µ1−µ2)−(µ1+µ2)Σ−1(µ1−µ2)>0 ⇔ X −1 2(µ1+µ2)| Σ−1(µ1−µ2)>0. Esta última desigualdade coincide con h( X )>0 . Para denir a versión mostral desta regra, simplemente temos que cambiar as cantidades poboacionais polas correspondentes mostrais obtidas a partir de X na fórmula dada pola función de decisión. Sexa X i un dos datos de X , rnorm( X ) = (1 se  X i−1 2(X1+X2)|S−1(X1−X2)>0 2 noutro caso . Mentres que se temos unha nova observación X new , tomaremos rnorm( X new)=1 se h( X new) =  X new −1 2(X1+X2)| S−1(X1−X2)>0, onde S e Xν , ν= 1,2 se calculan a partir de X sen ter en conta X new . Observación 2.6 . Algunhas veces refírese á regra discriminante linear normal como regra de Fisher, aínda que, como en Koch (2014), as denimos de maneiras distintas. A primera baséase na función de verosimilitude mentres que a segunda emprega autovectores dunha matriz. As formas das funcións de decisión asociadas son moi parecidas, e esta comparación motiva unha clase máis xeral de funcións de decisión lineais para problemas con dúas clases. Consideramos hβ( X ) =  X −1 2(µ1+µ2)| β, onde β é un vector adecuado, á nosa elección. Serían β=η para a regra de Fisher e β= Σ−1(µ1−µ2) para a regra linear normal. 24 CAPÍTULO 2. TÉCNICAS DE CLASIFICACIÓN Denindo W a matriz de variabilidade within-class e usando o Teorema 2.2 temos que W= 2 X ν=1 Σν= 2Σ, B= 2 X ν=1 (µν−µ1+µ2 2)(µν−µ1+µ2 2)|=1 2(µ1−µ2)(µ1−µ2)| e η=W−1(µ1−µ2) kW−1(µ1−µ2)k é o autovector asociado ao máximo autovalor de W−1B . Hai certas situacións nas que non se cumplen as hipóteses do Teorema 2.5 e a regra linear normal non é a adecuada. Isto ocorre cando sabemos que a distribución dos vectores non é normal, cando as matrices de covarianzas das clases non coinciden ou cando as descoñecemos. No caso de descoñecer Σν pero poder asegurar que son a mesma para ambas clases, temos a opción de empregar a matriz de covarianzas mostral agrupada (pooled) Spool = 2 X ν=1 nν−1 n−2Sν, sendo Sν a matriz de covarianzas mostral de Cν e nν o número de observacións desta clase. Antes de facer uso da matriz agrupada, teremos que asegurarnos de que é un procedemento apropiado. A regra linear normal pode non comportarse tan ben para vectores aleatorios con distribucións non coñecidas ou non normais. Aínda así, cando non se desvían demasiado da distribución gaussiana, podemos obter bos resultados con rnorm . Determinar canto é ese demasiado é difícil e para datos de dimensión moderada poderemos levar a cabo tests de normalidade ou tirar de axudas visuais. De todos modos, o proceder correcto é aplicar varias regras e avaliar o seu comportamento. Todo o introducido neste apartado pode estenderse a máis de dúas clases, baseándose igualmente na verosimilitude. Supoñamos que un vector aleatorio X pertence a unha das clases Cν=N(µν,Σ) , con ν≤κ só con medias distintas. Trataremos de atopar o parámetro µk que maximice a verosimilitude para asignar X a Ck . Denimos as funcións de decisión preferenciais h(l,ν)( X ) =  X −µl+µν 2| Σ−1(µl−µν) para l, ν = 1,2, ... , κ e l6=ν. Nótese ademáis que h(l,ν)=−h(ν,l) . Poñamos entón hnorm( X ) = m´ax (l,ν)=1,2,...,κ h(l,ν)( X ) e rnorm( X ) = k, 2.2. CLASIFICACIÓN BAIXO HIPÓTESE DE NORMALIDADE 25 sendo k o primeiro índice do par (l, ν) no que se alcanza o máximo hnorm( X ) . Alternativamente, podemos considerar as densidades fν( X ) para ν≤κ , e denir a regra ˜ rnorm( X ) = k se fk( X ) = m´ax 1≤ν≤κfν( X ). Para a poboación, as dúas regras asignarán X á mesma clase. Ao levalas ao caso da mostra pode que den lugar a dúas regras distintas, dependendo da variabilidade dentro das clases e a elección do estimador S para a matriz de covarianzas común Σ . 2.2.3. Regras discriminantes cuadráticas Ata agora consideramos únicamente regras lineais da forma h( X ) = a| X +c , para un vector a e un escalar c que non dependen do vector aleatorio X . Na regra linear normal a linealidade é unha consecuencia de que a mesma matriz de covarianzas se usa para tódalas clases. Mentres que, se sabemos que ditas matrices son diferentes (tanto que non nos serve agrupalas en Spool ) teremos que considerar matrices distintas. Isto levaranos a unha función de decisión non linear e, consecuentemente, a unha regra non linear. Consideremos entón un vector aleatorio X que pertence a unha das clases Cν=N(µν,Σν) , para ν≤κ , que se diferencian tanto no vector de medias como na matriz de covarianzas. Partimos outra vez da regra de Bayes e a densidade f para denir a regra, pero neste caso non só dependerá do parámetro µν , senón que nos interesa θν= (µν,Σν) . Asignaremos X á clase Cl se fl( X )> fν( X ) para l6=ν. Teorema 2.7. Para ν≤κ e un vector aleatorio X pertencente a unha das clases Cν= N(µν,Σν) . A regra discriminante rquad baseada nas funcións de densidade das κ clases asigna X a Cl se k X Σlk2+ log[det(Σl)] = m´ın 1≤ν≤κk X Σνk2+ log[det(Σν)], onde X Σ= Σ−1/2( X −µ) é o vector X estandarizado, é dicir, X Σ∼Nd(0, I) . Chamaremos a rquad regra discriminante cuadrática (normal) porque é cuadrática en X . Posto que derivamos esta regra de asumir unha distribución normal, non podemos esperar que o seu comportamento para vectores non gaussianos sexa o óptimo. Demostración. Obtemos coa denición da función de verosimilitude da normal multivariante que log f( X ) = −d 2log(2π)−1 2log[det(Σ)] + ( X −µ)|Σ−1( X −µ) 26 CAPÍTULO 2. TÉCNICAS DE CLASIFICACIÓN =−1 2k X Σk2+ log[det(Σ)]+c, sendo c=−d 2log(2π) independente do parámetro θ= (µ,Σ) . Asignaremos X preferentemente a Cl antes que a Cν cando fl( X )> fν( X ) ⇔log fl( X )>log fν( X ) ⇔ −1 2k X Σlk2+ log[det(Σl)]+c > −1 2k X Σνk2+ log[det(Σν)]+c ⇔ −1 2k X Σlk2+ log[det(Σl)]>−1 2k X Σνk2+ log[det(Σν)] ⇔ k X Σlk2+ log[det(Σl)] <k X Σνk2+ log[det(Σν)] E o resultado séguese desta desigualdade, tendo en conta tódalas clases e non só un par. Para a situación de dúas clases pódese obter unha expresión explícita máis sinxela para a regra rquad . Corolario 2.8. Na mesma situación que a do Teorema 2.7 con κ= 2 , supoñemos a maiores que ambas matrices de covarianzas son de rango r . Sexa Σν= ΓνΛνΓ| ν a descomposición espectral de Σν , con λν,j o autovalor j -ésimo de Λν . A regra cuadrática rquad asigna X á clase C1 cando kΛ−1/2 1Γ| 1( X −µ1)k2−kΛ−1/2 2Γ| 2( X −µ2)k2+ r X j=1 log λ1,j λ2,j <0. Nótese que se o rango r da matriz é estrictamente menor que a súa dimensión d , entón Σ ten a descomposición espectral Σ=ΓrΛrΓ| r , onde Γr é de tamaño d×r e Λr unha matriz diagonal de tamaño r×r . Ademáis, Γr non é ortogonal, e tense unha relación máis débil, Γ| rΓr=Ir×r , pero ΓrΓ| r6=Id×d . Ás matrices que cumplen isto chámaselles r -ortogonais . Tense así, para k, m ∈Z Σk/m = ΓΛk/mΓ|, e en particular, Σ−1/2= ΓΛ−1/2Γ| . Demostración. Probamos o resultado para o caso en que r=d , i.e., Σν son invertibles. Partimos dos cálculos feitos para a demostración do Teorema 2.7, pois rquad( X ) = 1 se k X Σ1k2+ log[det(Σ1)] <k X Σ2k2+ log[det(Σ2)] Tense que k X Σk2= X | Σ X Σ= [Σ−1/2( X −µ)]|[Σ−1/2( X −µ)] = 2.2. CLASIFICACIÓN BAIXO HIPÓTESE DE NORMALIDADE 27 [ΓΛ−1/2Γ|( X −µ)]|ΓΛ−1/2Γ|( X −µ)=( X −µ)|ΓΛ−1/2Γ|ΓΛ−1/2Γ|( X −µ) = ( X −µ)|ΓΛ−1/2Λ−1/2Γ|( X −µ) = [Λ−1/2Γ|( X −µ)]|[Λ−1/2Γ|( X −µ)] = kΛ−1/2Γ|( X −µ)k2, e posto que o determinante dunha matriz é o produto dos seus autovalores, temos que det Σ = d Y j=1 λj. Gracias a estos dous cálculos, desenvolvemos a desigualdade k X Σ1k2+ log[det(Σ1)] <k X Σ2k2+ log[det(Σ2)] ⇔ k X Σ1k2−k X Σ2k2+ log[det(Σ1)] −log[det(Σ2)] <0 ⇔ kΛ−1/2 1Γ| 1( X −µ1)k2−kΛ−1/2 2Γ| 2( X −µ2)k2+ log   d Y j=1 λ1,j −log   d Y j=1 λ2,j <0 ⇔ kΛ−1/2 1Γ| 1( X −µ1)k2−kΛ−1/2 2Γ| 2( X −µ2)k2+ d X j=1 log λ1,j − d X j=1 log λ2,j <0 ⇔ kΛ−1/2 1Γ| 1( X −µ1)k2−kΛ−1/2 2Γ| 2( X −µ2)k2+ d X j=1 [log λ1,j −log λ2,j]<0 ⇔ kΛ−1/2 1Γ| 1( X −µ1)k2−kΛ−1/2 2Γ| 2( X −µ2)k2+ r X j=1 log λ1,j λ2,j <0. Exemplo 2.9. Neste exemplo imos considerar simulacións por pares. En cada unha delas, haberá dúas clases de datos tridimensionais con vectores de medias distintos. En cada par, unha simulación terá matrices de covarianzas iguais e na outra serán distintas. Deste xeito poderemos ver se supón unha mellora empregar a regra cuadrática. Consideramos os vectores e matrices seguintes. µ1=    0 0 0     ,µ2=    2 1 −0,5     ,Σ1=    1 0 0 0 0,25 0 000,25     e Σ2=    1 0,25 0,125 0,25 0,5 0 0,125 0 0,25     . En primeiro lugar, simulamos datos de dúas clases normais. Na primeira gráca da Figura 2.2 aparecen 200 observacións da clase N(µ1,Σ1) e 300 de N(µ2,Σ1) , mentres que na segunda a matriz cambia dunha clase a outra, C1=N(µ1,Σ1) e C2=N(µ2,Σ2) . Os puntos da clase C1 están en negro e os da clase C2 en vermello. Vemos que na da esquerda as nubes de puntos teñen a mesma forma, pero centradas en distintas medias e na da 34 CAPÍTULO 2. TÉCNICAS DE CLASIFICACIÓN partición de Rd sería con Ai=Qn i=1[kih, (ki+ 1)h), onde h > 0 é o tamaño dos cubos e ki∈Z . Podemos escribir a regra de forma matemática como: r( X ) = arg max 1≤ν≤κ nν X i=1 I{Yi=ν}I{ X i∈A( X )}, onde A( X ) = Aj cando X ∈Aj . A regra do histograma presenta o problema de que a regra é menos precisa cando X está preto do borde das celas que cando está no medio. Polo tanto os puntos cercanos aos bordes deberían ter menos peso na decisión da clase asignada á cela. Para solventalo podemos introducir a regra da ventana móbil , que é máis suave que a do histograma porque toma os datos a certa distancia do punto que queremos clasicar e asígnalle a etiqueta maioritaria entre estos datos. Formalmente defínese como: r( X ) = arg max 1≤ν≤κ nν X i=1 I{Yi=ν, X i∈Sx,h} onde h > 0 e Sx,h denota a bóla pechada de centro X e radio h . A regra kNN é un caso particular desta, na que h é a distancia ao k -ésimo dato máis cercano, que cambia para cada X que queiramos clasicar e depende da mostra, co que non é unha constante. Podemos facer unha regra aínda máis suave se damos máis peso a aqueles puntos máis cercanos a X que aos máis distantes. Sexa K:Rd→R unha función kernel , normalmente non negativa e monótonamente decrecente sobre os radios partindo da orixe. A regra discriminante kernel vén dada por r( X ) = arg max 1≤ν≤κ nν X i=1 I{Yi=ν}K X − X i h. Ao parámetro h chámaselle ancho de banda , e proporciona unha especie de ponderación da distancia. Canto maior é h , máis contan as observacións que non están tan cerca de X . Claramente, a regra kernel é unha xeneralización da da ventana móbil, co kernel naive K(x) = I{x∈S0,1} . Algunhas das funcións kernel máis empregadas son: Kernel gaussiano K(x) = e−kxk2 , Kernel de Cauchy K(x) = 1 1+kxkd+1 , Kernel de Epanechnikov K(x) = (1 −kxk2)I{kxk≤1} , sendo k·k a distancia euclídea. Na Figura 2.4 vese a forma destas funcións para unha dimensión. Os dous primeiros teñen como soporte todo R , mentres que o de Epanechnikov soamente é non nulo no intervalo [−1,1] . 2.3. TÉCNICAS DE CLASIFICACIÓN NON PARAMÉTRICA 35 Figura 2.4: Grácas das funcións kernel máis usuais. −3 −2 −1 0 1 2 3 0.0 0.4 0.8 Kernel gaussiano −3 −2 −1 0 1 2 3 0.0 0.4 0.8 Kernel de Cauchy −3 −2 −1 0 1 2 3 0.0 0.4 0.8 Kernel de Epanechnikov Observación 2.14 . As regras tipo kernel tamén se poden denir a partir da regra de Bayes, estimando a función de densidade de cada clase empregando métodos non paramétricos, ˆ fν( X ) = 1 nν nν X i=1 I{Yi=ν}K X − X i h, onde K é unha función kernel coma as anteriores e h é un parámetro de suavización. Esta aproximación é a que se toma en Klamelä (2014). Polo que, coñecidas as probabilidades a priori π1, ... , πκ , teriamos a expresión seguinte para a regra. r( X ) = arg max 1≤ν≤κ(πν 1 nν nν X i=1 I{Yi=ν}K X − X i h). Se empregamos as proporcións mostrais como probabilidades a priori πν=nν n , esta regra é equivalente á que denimos como regra tipo kernel. 36 CAPÍTULO 2. TÉCNICAS DE CLASIFICACIÓN Capítulo 3 Estimación do erro e avaliación de regras discriminantes Tras denir unha regra discriminante e estimala a partir da mostra, é preciso ter medidas da súa calidade. Neste capítulo deniremos algunhas medidas do erro de clasicación e ilustraremos os procedementos presentados ata agora a través dalgún exemplo. Pretendemos adiviñar Y a través da regra e a partir dos datos etiquetados, xa que normalmente non coñecemos a distribución do vector etiquetado. Polo tanto é esencial estimar a probabilidade de erro dunha regra L(r) para saber que esperar da calidade desta. É de especial interese estimar a probabilidade de erro óptima L∗ , xa que se é grande, sabemos que calquera regra empregada fará unha clasicación bastante pobre, ademáis de que comparar L(r) con L∗ dinos canto podemos chegar a mellorar a regra r . Exemplo 3.1. Consideramos que temos dúas poboacións normais tridimensionais con distintos vectores de medias pero a mesma matriz de covarianzas. Simularemos datos de situacións en que cambie a proporción de datos que pertencen a cada unha das clases e tamén a dicultade do problema de clasicación, determinada pola distancia entre as medias das clases. Co segundo punto do Teorema 1.11 poderemos calcular a regra de Bayes baixo suposicións de normalidade, considerando as probabilidades a priori ou non para así comparar o seu comportamento, entre elas e respecto ao erro de Bayes. Tomamos os seguintes parámetros para denir as clases: µ1=    0 0 0     ,µ2=    2 1 −0,5     ,Σ =     0,5 0 0 0 0,25 0 0 0 0,25     e µ3=    1 1 −0,5     , de xeito que C1=N(µ1,Σ) e C2=N(µ2,Σ) para o caso 'fácil' e o mesmo pero cambiando 37 38 CAPÍTULO 3. ESTIMACIÓN DO ERRO E AVALIACIÓN DE REGRAS DISCRIMINANTES µ2 por µ3 na segunda clase para o caso 'difícil'. Na Figura 3.1 móstranse simulacións para o caso fácil e o difícil, con 200 datos da clase C1 en negro e 100 da clase C2 en vermello. Figura 3.1: Datos simulados de dous problemas de distinta dicultade con dúas clases normais. −3 −2 −1 0 1 2 3 4 −2.0 −1.5 −1.0 −0.5 0.0 0.5 1.0 1.5 −2−1 0 1 2 ● ● ● ●● ● ● ●● ● ● ● ● ● ● ●● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ●● ●● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ●● ● ● ●● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ●● ● ● ● ●● ● ●● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ●● ●●● ● ● ● ● ●● ● ● ●● ●● ● ● ● ●● ● ● ● ● ● ●● ● ● ● ● ● ● ●● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ●●● ●● ● ● ● ● ● ● ● ● ● ●● ● ●● ●● ● ● ●● −3 −2 −1 0 1 2 3 4 −2.0 −1.5 −1.0 −0.5 0.0 0.5 1.0 1.5 −2−1 0 1 2 ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ●● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●● ● ● ● ●● ● ● ● ● ●● ● ● ● ● ● ●● ● ● ● ● ● ● ● ●● ● ● ●● ●●● ● ● ● ● ●● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●●● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ●● ●● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● Podemos calcular o erro de Bayes para estas dúas situacións, no suposto de coñecer as probabilidades a priori (2/3 dos datos son da primeira clase e 1/3 da segunda) e tamén de que nos falte esta información. Como notación, r∗ representa a regra sen usar π e rBayes a que sí as emprega. Faremos as contas só para o problema 'fácil' coa regra r∗ xa que as do resto de casos son análogas. A regra óptima é lineal, como vimos no Teorema 2.5, polo tanto L∗=P(r∗( X )6=Y) = P(r∗( X ) = 1|Y= 2)π2+P(r∗( X )=2|Y= 1)π1 =P(h( X )>0|Y= 2)π2+P(h( X )<0|Y= 1)π1, sendo h( X ) = −4X1−4X2+ 2X3+ 6,5. Debido a que a matriz de covarianzas Σ é diagonal, tense que as compoñentes Xi do vector aleatorio son normais unidimensionais independentes. Co cal, cando X ∈ C1 , h( X )∈ N(6,5,13) e cando pertence á segunda clase, h( X )∈ N(−6,5,13) . Os cálculos polo tanto darán que L∗=1 3P(Z > 6,5 √13) + 2 3P(Z < −6,5 √13 =P(Z > 6,5 √13) = 0,0357. 39 En cambio, para o caso difícil, tense que L∗= 0,1217 , é dicir, espérase que clasiquemos mal o 12% das observacións. O erro de Bayes podémolo calcular cando coñecemos os parámetros poboacionais, pero de non ser así, debe estimarse a partir das regras lineais estimadas. Para ver o comportamento das regras con e sen probabilidades a priori, repetiremos a simulación dos datos e cálculo da regra 100 veces. Así podemos calcular a media e tamén a desviación típica da porcentaxe de observacións mal clasicadas. Os resultados móstranse na seguinte táboa, para todos os casos considerados e cambiando o tamaño das mostras das poboacións, n1 e n2 . Cadro 3.1: Erro de clasicación medio e desviación típica (entre parénteses) para as dúas regras e os catro conxuntos de datos simulados para o caso fácil e o difícil, en tanto por cento. Tamaño da mostra Caso fácil Caso difícil n1n2r∗rBayes r∗rBayes 30 10 2.6 2.175 7.425 6.150 (2.796) (2.427) (4.955) (3.915) 300 100 3.6675 3.11 9.06 7.2425 (0.983) (0.779) (1.61) (1.257) 20 30 2.8 2.72 8.86 8.7 (2.412) (2.503) (4.151) (4.113) 200 300 3.422 3.306 9.174 8.896 (0.799) (0.769) (1.219) (1.146) Observando os resultados do cadro, os erros cometidos no caso difícil son maiores para todas as regras e tamaños mostrais. Ademáis, canto menos se parecen as probabilidades a priori (que neste caso son a proporción de datos de cada clase), a vantaxa que supón empregar rBayes é maior respecto a usar r∗ , sobre todo nos problemas difíciles. Podemos salientar tamén que para as mostras de maior tamaño n cométese un porcentaxe maior de erros de clasicación (aínda que a desviación típica é máis pequena) que para o mesmo problema pero con menos observacións. No exemplo supuxemos que as probabilidades a priori son correctas, pero de non selo poden empeorar a toma de decisións. 40 CAPÍTULO 3. ESTIMACIÓN DO ERRO E AVALIACIÓN DE REGRAS DISCRIMINANTES Posto que en xeral quen constrúe a regra non coñece a distribución dos datos, son necesarios métodos de estimación ou cotas da probabilidade de erro que non dependan da distribución de X . Mediremos a calidade dunha regra construida a partir de X= [X[1]X[2] ···X[κ]] coa probabilidade de erro condicional P(r( X ;X)6=Y|X) . O problema disto é que avaliar a regra sobre os mesmos datos usados para construíla, lévanos a un sesgo optimista. Igual que debemos comprobar o funcionamento dunha regra con distintas medidas do erro, tamén debemos achar as limitacións das distintas medidas do erro empregándoas sobre varias regras. Denimos a continuación un par de medidas de erro de clasicación, estimadores naturais da probabilidade de erro que se calculan a partir dos datos. Denición 3.2. Sexan X= [ X 1 X 2··· X n] datos etiquetados e r unha regra. Un factor de costes c= [c1···cn] asociado a r defínese como ci(= 0 se r clasicou X i correctamente >0 se a clasicou incorretamente . O erro de clasicación mis da regra r e factor de costes c é o erro dado por mis =1 n n X i=1 ci. Cando X i é clasicado incorrectamente, adoitamos usar ci= 1 . Un erro de clasicación con factor de costes 0 e 1 é natural pois, sen importar o número de clases que haxa, o que fai é contar o número de observacións que foron mal clasicadas. De todos modos, outros factores de costes, positivos e non constantes, son útiles. Por exemplo, cando queremos distinguir entre falsos positivos e falsos negativos nunha clasi- cación binaria, como sería identicar clientes morosos. Pode que nos interese penalizar máis aqueles erros nos que se predí que un cliente pagará cando non o vai facer que a situación en que un cliente non moroso se clasique como moroso, pois a primeira pode poñer nun risco maior á empresa. Outro exemplo sería que as clases indiquen o grao de severidade dunha característica. Se temos un cancro en estado IV, que xa se estendeu a partes distantes do corpo, podemos dar pesos maiores ao erro de clasicalo nun estado 0, no que só hai presenza dalgunhas células anormais, que ao de clasicalo como estado II, que xa hai cancro presente e estendido a tecidos cercanos. Así reexaríamos que o primeiro erro é máis grave que o segundo. A fase de aprendizaxe consiste na construción dunha regra empregando os datos etiquetados X0 e a fase de predicción en aplicar dita regra a Xnew para predicir as respostas Y new . 41 Normalmente, quereremos empregar tódolos datos dipoñibles para entrenar o clasicador, co cal só temos unha mostra de entrenamento, pero non de testeo. Primeiro derivamos unha regra r0 da submostra X0 de X , e calculamos o erro de clasicar as observacións Xp que nos quedan de X empregando r0 . X0 é a mostra de entrenamento e Xp a mostra de testeo . Este proceso é a base da idea de entrenamento e testeo da Aprendizaxe Supervisada. Podemos escoller X0 de moitas maneiras, a máis simple sería usar tódalas observacións excepto unha. Denición 3.3. Consideramos os datos etiquetados X= [ X 1 X 2··· X n] e a regra r . Para cada i≤n tomamos X0,(−i)= [ X 1··· X i−1 X i+1 ··· X n] , ao que chamamos conxunto de entrenamento leave-one-out (deixar un fóra) i -ésimo, e Xp,(−i)= X i . Sexa r(−i) a regra construida igual que r pero baseada en X0,(−i) , considerando así X i unha nova observación. Un factor de costes k= [k−1k−2···k−n] asociado coas regras r(−i) é k−i(= 0 se X i foi clasicada correctamente por r(−i) >0 se a clasicou incorretamente . Co cal o erro leave-one-out , loo , baseado nas n regras r(−i) , os valores asignados r(−i)( X i) e o factor de costes k será loo =1 n n X i=1 k−i. Así, o método leave-one-out consistirá na elección dos conxuntos de entrenamento X0,(−i) , a construción das regras r(−i) e o cálculo do erro loo . Deste xeito cada observación déixase fóra exactamente unha vez. Adóitase tomar costes 0 e 1, xa que así contamos o número de erros de clasicación empregando r(−i) . Polo xeral, tódalas regras r(−i) serán distintas, e diren tamén de r . Un bo clasicador caracterizarase porque esta diferencia se faga pequena a medida medra o tamaño da mostra. Se comparamos mis e loo , mis debería ser máis pequeno pois cada punto é parte do conxunto de entrenamento da regra usada para clasicalo, pero loo introduce a maiores unha medida do erro da predicción. Moitas veces, será preferible tomar un conxunto Xp con máis dun elemento, para analizar mellor o comportamento da regra fóra da mostra. Outro modo de dividir os datos dispoñibles entre mostra de entrenamento e de testeo é a seguinte: Denición 3.4. Sexan X= [ X 1 X 2··· X n] un conxunto de datos etiquetados e r unha regra. Sexan k e m enteiros tales que km =n . Para j≤m , denimos o conxunto de 42 CAPÍTULO 3. ESTIMACIÓN DO ERRO E AVALIACIÓN DE REGRAS DISCRIMINANTES entrenamento j -ésimo X0,j e o conxunto de testeo j -ésimo Xp,j como X0,j = [ X 1··· X (j−1)k X jk+1 ··· X n] Xp,j = [ X (j−1)k+1 ··· X jk]. Sexa rj a regra derivada soamente de X0,j e usámola para calcular o erro de clasicación para cada X i do conxunto Xp,j . O erro j sobre Xp,j é j= jk X i=(j−1)k+1 ci sendo ci= 0 cando rj( X i=Yi e ci>0 noutro caso. O erro de validación cruzada m veces ( m -fold cross-validation), cv é cv =1 n m X j=1 j. Chámaselle entón validación cruzada m veces á partición dos datos en {X0,j,Xp,j} m veces xunto coas regras rj e os erros j e cv . Neste proceso, as mostras de testeo sempre teñen k vectores e son disxuntas, de xeito que cada observación X i dos datos orixinais X é contada unha única vez de cara ao cálculo de cv . Pódese ver fácilmente que o erro leave-one-out é un caso particular de cross-validation con m=n e k= 1 . Para este método, hai que dividir os datos e entrenar a regra m veces, o cal supón un alto coste computacional, que medra tamén canto máis grande é m . Se k é demasiado grande, entón o conxunto de entrenamento pode ser demasiado pequeno, o cal pode afectar á precisión da regra. Debemos achar valores para m e k de tal xeito que non cheguemos a ter estos problemas, por exemplo tomar m entre 10 e 20 mentres que entre o 10 e 20% dos datos sexan de testeo en cada partición. A validación cruzada é moi costosa computacionalmente, polo que moitas veces emprégase unha única partición dos datos. Este proceder non é o mellor cando queremos comparar distintos métodos ou regras. Na clasicación, o número de clases e o número de observacións que pertencen a cada unha delas tamén xogan un papel importante. Se temos dúas clases, unha delas con moitos máis vectores que a outra, debemos escoller coidadosamente o tamaño das mostras de entrenamento e testeo. Estas dúas últimas medidas do erro son de especial interese, pois achegan información sobre a capacidade predictora da regra. En moitas ocasións, unha regra estímase a partir dunha mostra de xeito que se minimiza o erro de clasicación destos datos coñecidos, pero ao aplicala a novos datos comete un erro moito maior. Isto coñécese como overtting e 43 supón un problema ao construir unha regra discriminante, sobre todo cando o número de variables é moi grande, como veremos no seguinte capítulo. Exemplo 3.5. Volvemos a considerar os catro conxuntos de datos simulados do Exemplo 2.9, dous deles normais, e outros dous non normais, con igualdade de matriz de covarianzas e sen ela. Para clasicalos usaremos regras tipo kernel, considerando distintas funcións kernel (normal, o de Epanechnikov e uniforme) para poder comparar os resultados obtidos. Empregarase a función classif.np do paquete fda.usc . Esta función calcula de maneira empírica o parámetro h que leva á clasicación óptima, empregando unha modicación do criterio de validación cruzada, e devolve tamén unha estimación da probabilidade de clasicación correcta para cada clase, indicada na seguinte táboa. Cadro 3.2: Probabilidades de clasicación correcta por clase, para cada función kernel considerada. N1 N2 P1 P2 C1C2C1C2C1C2C1C2 Normal 0.835 0.943 0.93 0.9867 0.825 0.92 0.895 0.903 Epanechnikov 0.83 0.94 0.9 0.9967 0.84 0.92 0.91 0.8767 Uniforme 0.815 0.95 0.875 0.9967 0.77 0.933 0.885 0.8967 Pódese apreciar que a probabilidade de clasicar correctamente observacións da segunda clase é maior en tódolos casos, isto débese a que esta clase é a maioritaria nos nosos conxuntos de datos. En canto a determinar cal dos tres kernel é máis adecuado para cada unha das situacións, os resultados non son concluíntes, pois todos teñen un poder predictivo similar e bastante elevado. Para o primeiro caso, podemos calcular o erro de Bayes partido da función de densidade da normal multivariante. Para N1 sería L∗= 0,097 , co que se clasicarían correctamente o 90,3% dos datos, que é máis ou menos o que se obten cas regras tipo kernel. L∗ é o valor ao que debería converxer a probabilidade de erro ao medrar o tamaño mostral se usásemos unha estimación da regra de Bayes. Este é o caso das regras tipo kernel, nas que a estimación se fai con técnicas non paramétricas, a diferenza da regra linear. 50 CAPÍTULO 4. CONSIDERACIÓNS SOBRE A ALTA DIMENSIÓN E BIG DATA 4.2. Análise dicriminante linear regularizado disperso: sparse rLDA O Sparse Regularized Linear Discriminant Analysis (sparse rLDA) presentado por Qiao et al. (2008) parte da regra de Fisher e introduce a selección de variables impoñendo que o autovector η sexa disperso, é dicir, que case todas as súas compoñentes sexan nulas. Con isto pretenden identicar aquelas variables que son determinantes para diferenciar as clases entre a estrutura de covarianzas que as variables do problema poidan presentar e así descartar información redundante que perxudique a clasicación. Cando a matriz de variabilidade within-class W é singular, substitúese no problema pola matriz regularizada ˜ W presentada ao comezo deste capítulo. Dise neste artigo que a elección do parámetro γ non inúe nos resultados obtidos, sempre que este sexa pequeno. O rLDA consiste en empregar esta matriz regularizada, achar a primeira dirección discriminante da mesma, η e proxectar os datos sobre ela. Estas proxeccións considéranse como novas variables discriminantes unidimensionais, ás que se lle aplican regras discriminantes como a do centroide máis cercano, SVM ou unha regra kNN, por exemplo. Chegando así a unha clasicación dos datos orixinais. Este vector η pode ter moitas compoñentes non nulas, de xeito que tódalas variables dos datos inuirán na clasicación. Ao introducir o sparse rLDA, impoñen que η teña só unhas poucas entradas non nulas e están facendo unha selección de variables e eliminando a información redundante. Para obter variables discriminantes dispersas, primeiro relacionan o vector η co vector de coecientes dunha regresión transformando o problema de autovalores do Teorema 2.2 nun problema de regresión, e despois plantéxano como un problema de mínimos cadrados no que introducen unha penalización da norma L1 do vector de coecientes na función obxectivo igual que no problema LASSO de Tibshirani (1996). 4.3. Análise discriminante en alta dimensión: HDDA Nesta sección, presentaremos o método introducido por Bouveyron et al. (2007), ao que chaman High Dimensional Discriminant Analysis (HDDA). Baséase en que cando a dimensión é moi grande ocurre o fenómeno do espacio vacío e podemos supoñer que os datos viven en subespacios de dimensión menor. O que fai é reducir a dimensión para cada clase Cν de forma independente e impón unha estrutura determinada sobre as matrices de covarianzas Σν para adaptar un contexto gaussiano á alta dimensión, reducindo o número 4.3. ANÁLISE DISCRIMINANTE EN ALTA DIMENSIÓN: HDDA 51 de parámetros a estimar. Suponse que as clases son esféricas (é dicir, a matriz de covarianzas é un múltiplo da identidade) nestes subespazos ou o que é o mesmo, que Σν teñen só dous autovalores distintos. Igual que nos métodos de análise discriminante clásicos, supoñemos a normalidade das clases Cν=N(µν,Σν) para ν∈ {1, ... κ} . Posto que Σν son simétricas, gracias ao Teorema espectral obtemos unha descomposición matricial Σν=Qν∆νQ| ν , onde Qν é unha matriz ortogonal cuxas columnas son unha base de autovectores de Σν e ∆ν é unha matriz diagonal formada polos autovalores de Σν . Supoñemos que ∆ν ten dous autovalores diferentes, aν> bν . Chamamos Eν ao subespacio afín de dimensión dν xerado polos autovectores asociados ao autovalor aν con µν∈Eν , e sexa E⊥ ν tal que Eν⊕E⊥ ν=Rd con µν∈ E⊥ ν . Consideramos Pν( x ) = ˜ Qν˜ Qν |( x −µν) + µν a proxección de x sobre Eν , onde ˜ Qν é a matriz formada polas dν primeiras columnas de Qν e o resto ceros. Análogamente, P⊥ ν( x )=(Qν−˜ Qν)(Qν−˜ Qν)|( x −µν) + µν é a proxección de x sobre E⊥ ν . Partindo da regra de Bayes e impoñendo a forma descrita da matriz de covarianzas, a regra discriminante asignará un vector aleatorio X á clase Cν que minimice a expresión hν( X ) = kµν−Pν( X )k2 aν +k X −Pν( X )k2 bν +dνlog aν+ (d−dν) log bν−2 log πν. Vexamos como se pode chegar a ela a partir da función de densidade dunha distribución normal. Asignarase X á clase Cν que maximice πνfν( X ) , ou o que é o mesmo, que minimice −2 log πνfν( X ) . −2 log πνfν( X ) = −2 log πν(2π)−d/2|Σν|−1/2exp −1 2( X −µν)|Σ−1 ν( X −µν) =−2 log πν+dlog(2π) + log |Σν|+ ( X −µν)|Σ−1 ν( X −µν) =−2 log πν+dlog(2π) + log adν νb(d−dν) ν+ ( X −µν)|Qν∆−1 νQ| ν( X −µν) =−2 log πν+dlog(2π) + dνlog aν+ (d−dν) log bν +( X −µν)|(˜ Qν+Qν−˜ Qν)( 1 aν ˜ Qν+1 bν (Qν−˜ Qν))|( X −µν) =−2 log πν+dlog(2π) + dνlog aν+ (d−dν) log bν +1 aν ( X −µν)|(˜ Qν+Qν−˜ Qν)˜ Qν |( X −µν)+ 1 bν ( X −µν)|(˜ Qν+Qν−˜ Qν)(Qν−˜ Qν)|( X −µν) =−2 log πν+dlog(2π)+ dνlog aν+ (d−dν) log bν+1 aνkµν−Pν( X )k2+1 bνk X −Pν( X )k2. Como dlog(2π) é constante para tódalas clases, elimínase e chegamos á expresión de hν . 52 CAPÍTULO 4. CONSIDERACIÓNS SOBRE A ALTA DIMENSIÓN E BIG DATA Introducimos a seguinte notación para simplicar a interpretación da regra: aν=σ2 ν αν e bν=σ2 ν 1−αν , con αν∈[0,1] e σ2 ν>0 . Así a expresión anterior pódese reescribir como hν( X ) = 1 σ2 νανkµν−Pν( X )k2+ (1 −αν)k X −Pν( X )k2 +2dlog(σν) + dνlog 1−αν αν−dlog(1 −αν)−2 log πν. Para certos valores dos parámetros αν e σ2 ν , téñense regras particulares das que xa falamos anteriorimente. Se αν= 1/2∀ν , esta regra non é máis que a regra cuadrática coa suposición adicional de que a matriz de covarianzas é un múltiplo da identidade, Σν=σ2 νI . Se a maiores o parámetro σ2 ν=σ2 é o mesmo para todas as clases, é a regra linear esférica, con Σ = σ2I . Deixando xos algúns, pero non todos os parámetros involucrados na regra HDDA, pódense obter moitos modelos distintos con claras interpretacións xeométricas. Presentamos dúas opcións: Regra discriminante isométrica (HDDAi): Fanse as seguintes suposicións: αν=α, σν=σ, dν=d∗ e πν=π∗∀ν≤κ, de xeito que a expresión que queremos minimizar será hν( X ) = αkµν−Pν( X )k2+ (1 −α)k X −Pν( X )k2. Se α= 0 , HDDAi asignará X a Cl se d( X , El)< d( X , Eν)∀ν6=l . É dicir, levará X á clase asociada ao subespazo Eν máis cercano. Se α= 1 , r( X ) = l cando d(µl, Pl( X )) < d(µν, Pν( X )) ∀ν6=l , é dicir, leva X á clase cuxa media está máis cerca da proxección de X sobre o subespazo. Cando 0< α < 1 , a regra será unha mestura destas dúas. Será necesario tamén estimar α , pero discutirémolo máis adiante. Regra discriminante homotécica (HDDAh): A diferencia deste método co anterior é que non se impón que σν sexa constante para tódalas clases. Así a expresión queda hν( X ) = 1 σ2 ναkµν−Pν( X )k2+ (1 −α)k X −Pν( X )k2+ 2dlog σν. Posto que σν está no denominador, se X está á mesma distancia de dúas clases, a regra favorecerá aquela cuxa varianza sexa maior. Obviamente, pode haber situacións en que as hipóteses anteriores sexan demasiado restrictivas, co cal optaremos por empregar a expresión xeral. 4.3. ANÁLISE DISCRIMINANTE EN ALTA DIMENSIÓN: HDDA 53 4.3.1. Estimación dos parámetros A estimación de parámetros é necesaria para calquera regra, pois dispoñemos dunha mostra e non da poboación total. Neste caso empregaremos os estimadores de máxima verosimilitude calculados a partir dunha mostra X de tamaño n . Igual que xemos para as regras lineais e cuadráticas, estimaremos as probabilidades a priori polas proporcións mostrais, as medias de cada clase polas medias mostrais e as matrices Σν polas matrices de covarianzas mostrais. Podemos supoñer nun principio que coñecemos a dimensión dos subespazos dν , e así obtéñense os estimadores ˆaν=1 dν dν X j=1 λν,j e bν=1 d−dν d X j=dν+1 λν,j, onde λν,j son os autovalores de Sν . A columna j -ésima de Qν estímase polo autovector de Sν asociado ao autovalor λν,j . Así aν e bν son estimados polas varianzas empíricas da clase Cν nos subespazos Eν e E⊥ ν respectivamente. Pódense entón deducir os seguintes estimadores: ˆαν=ˆ bν ˆaν+ˆ bν e ˆσν2=ˆaνˆ bν ˆaν+ˆ bν . Por último, queda achar os parámetros dν . A proposta de Bouveyron et al. (2007) é empregar o método empírico dos scree plots , que analiza a diferencia entre os autovalores para atopar unha ruptura na gráca. Este método baséase en que o autovalor λν,j representa a fracción da varianza total correspondente ao j -ésimo autovector de Σν . Escollemos a dimensión a partir da cal as diferencias son moi pequenas respecto á máxima das diferencias. Os scree plots empréganse tamén na ACP para decidir o número de compoñentes principais que se consideran para representar os datos, como se explica no segundo capítulo de Koch (2014). Un dos problemas co que nos atopamos ao facer isto é que a medida que d aumenta, os autovalores aportan menos información ao total da varianza, e é máis difícil identicar o índice dν a partir do cal o poder explicativo decrece signicativamente. No último capítulo ilustraremos con datos simulados como obter a estimación dos parámetros dν , ademáis de comparar os resultados obtidos co HDDA respecto aos obtidos con métodos clásicos cando a dimensión medra. 54 CAPÍTULO 4. CONSIDERACIÓNS SOBRE A ALTA DIMENSIÓN E BIG DATA Capítulo 5 Ilustración sobre datos simulados e reais Neste último capítulo farase unha simulación de datos empregando R e tamén se tratará unha base de datos real, que mostre como para unha dimensión d grande comezan a fallar os métodos clásicos e os introducidos para o contexto de Big Data funcionan mellor. Implementarase o método HDDA de Bouveyron et al. (2007) no software R. Neste artigo móstranse os resultados obtidos ao aplicalo ao recoñecemento de obxectos en imaxes reales en comparación con métodos de clasicación clásicos. Agora será aplicado sobre os dous exemplos estudados en Qiao et al. (2008) para poder comparar o comportamento das dúas propostas. Empregaremos o erro de clasicación mis como medida da calidade das regras, aplicadas sobre o conxunto de entrenamento e tamén sobre o de testeo. Deste xeito mediremos o poder clasicador das regras fóra da mostra. 5.1. Datos simulados Simularemos un conxunto de datos de entrenamento de tamaño 25 para cada unha das dúas clases e un conxunto de datos de proba de tamaño 100 para cada clase tamén. Os datos serán de dimensión d= 100 , polo que nos atopamos nun caso HDLSS. Das 100 variables que forman o vector aleatorio X , só as dúas primeiras serán diferentes entre as clases C1 e C2 . Os datos seguirán unha distribución normal cos seguintes parámetros: X ∼ N2(µν,Σ) Nd−2(0,Ip−2)!, para ν= 1,2, onde µ1= 0 0,9!,µ2= 0 −0,9! e Σ = 1 0,7 0,7 1 !. 55 56 CAPÍTULO 5. ILUSTRACIÓN SOBRE DATOS SIMULADOS E REAIS Claramente, a clasicación depende únicamente das dúas primeiras variables, e tódalas demáis aportan información redundante. Identicar aquelas variables sucientes para a clasicación axudará a evitar o overtting da regra sobre a mostra de entrenamento e levará a mellores resultados sobre datos fóra desta mostra. O enfoque de Qiao et al. (2008) céntrase especialmente no fenómeno chamado data piling . Ao empregar unha estimación da matriz de covarianzas ou da variabilidade withinclass, a cal non é nada able cando d>n , as direccións sobre as que se proxectan os datos para reducir a dimensión poden ser moi distintas das teóricas, levando a un sobreaxuste da regra aos datos de entrenamento. O que ocorre é que as proxeccións dos datos de cada clase están moi separados, pero ao proxectar o conxunto de testeo hai un solapamento moito maior ca usando as direccións teóricas. Eles introducen a dispersión no modelo para evitar que isto pase e ilustran como as direccións discriminantes dispersas estimadas se parecen máis ás teóricas, levando a unha pérdida de precisión sobre os datos mostrais pero aumentándoa fóra da mostra. Ao implementar o método HDDA, atopámonos co problema de determinar a dimensión intrínseca dν dos subespazos de autovectores de S . Non nos podemos axudar dos scree plots ao estar nun contexto HDLSS, co cal para identicar o mellor dν , estimaremos a regra para valores de dν entre 1 e 20 e escolleremos aquel que nos leve ao menor erro de clasicación. Repetirase a simulación dos datos e o cálculo da regra 50 veces para poder promediar os erros de clasicación. Tomarase dν igual para as dúas clases por simplicidade. As suposicións da regra discriminante isométrica teñen sentido neste caso, xa que ambas clases aparecen na mesma proporción na mostra e teñen a mesma matriz de covarianzas teórica. Tamén podemos esperar que a regra isométrica leve a mellores resultados que a homotécica, ao supoñer que σ1=σ2 . Implementaranse soamente estas dúas versións particulares do HDDA. Na Figura 5.1 observamos como o erro de clasicación no conxunto de entrenamento, en negro para HDDAi e vermello para HDDAh, decrece ao aumentar dν , pero fóra da mostra, en verde para HDDAi e azul para HDDAh, aumenta. Tomaremos polo tanto o menor valor posible dν= 1 , para evitar na medida do posible o overtting . As liñas horizontais debuxadas mostran o erro promedio ao aplicar a regra linear, que non depende de dν . Estos valores foron calculados coa función lda do paquete MASS de R, que emprega unha pseudoinversa para lidiar coa singularidade da matriz de covarianzas mostral. Está claro que empregando calquera das versións de HDDA obtemos mellores resultados que con esta regra clásica. Para o valor óptimo de dν , obtemos un erro do 14% 5.1. DATOS SIMULADOS 57 Figura 5.1: Erro de clasicación medio para as distintas regras (HDDAi, HDDAh e regra linear) sobre os conxuntos de entrenamento e testeo simulados, fronte aos distintos valores considerados para a dimensión intrínseca dν . 5 10 15 20 0.0 0.1 0.2 0.3 0.4 0.5 sobre a mostra e do 43% sobre o conxunto de testeo, en azul claro e rosa respectivamente. Como era de esperar, HDDAi da mellores resultados que HDDAh. Na táboa 5.1 resúmense os resultados obtidos con estes métodos e os que aparecen en Qiao et al. (2008). As tres primeiras columnas correspóndense co lda , e as dúas variantes do HDDA implementadas, e mostran os datos obtidos facendo unha simulación que imita a feita por Qiao et al. (2008). As dúas últimas columnas conteñen os valores proporcionados no seu artigo, pois non se implementou o método nin se repetiu a simulación neste traballo por cuestións de tempo. Pódese facer entón unha comparación dos métodos, pero sobre todo é unha ilustración das técnicas presentadas e un exercicio computacional. Podemos observar que o método implementado de Bouveyron et al. (2007) chega a uns resultados similares aos da regra linear regularizada rLDA, pero non tan bos coma os da sparse rLDA (calculado para 5 coecientes non nulos, é dicir, 5 variables signicativas). 58 CAPÍTULO 5. ILUSTRACIÓN SOBRE DATOS SIMULADOS E REAIS Cadro 5.1: Erro cometido por cada método considerado no conxunto de entrenamento e no de testeo (en tanto por cento). LDA HDDAi HDDAh rLDA sparse rLDA Entrenamento 13,9 2,6 3,7 0 12 Testeo 43,3 31,8 34 32 13.5 Tanto o HDDA coma o rLDA fan unha reducción da dimensión e solvéntase o problema de que as matrices de covarianzas sexan singulares regularizándoas, pero non hai unha selección de variables implícita. O método disperso busca reducir explícitamente o número de variables explicativas, e por iso obtén a mellor clasicación fóra da mostra. Esta simulación era de tal xeito que só dúas variables diferenciaban as clases, ten sentido que unha regra que dependa de moi poucas variables funcione mellor. 5.2. Datos de medidas de expresión xénica Consideramos o conxunto de datos Colon de Alon et al. (1999), que contén 42 mostras de tecido tumoral e 20 de tecido de colon normal. Para cada mostra hai 2000 medidas do nivel de expresión de xens. O obxectivo é clasicar as mostras normais e tumorais en función das medidas de expresión xénica. Estos son os datos empregados en Qiao et al. (2008) para ilustrar o sparse rLDA nun caso real. Fixeron unha análise de dúas etapas, repetida 50 veces para promediar a proporción de erros obtida. Posto que 2000 variables suporían un esforzo computacional moi grande, reducen a dimensión seleccionando as 200 máis signicativas para axilizar os cálculos. Só se dispón dunha mostra de 62 casos en total, que se debe particionar para obter conxuntos de entrenamento e de testeo. Divídese en 2/3 e 1/3 das observacións, respectivamente, de xeito que a proporción de datos tumorais e normais sexa equilibrada. Aos datos resultantes con d= 200 , aínda nunha situación HDLSS, aplícanselle os métodos rLDA e sparse rLDA usando o centroide máis cercano, SVM e o veciño máis cercano sobre as proxeccións obtidas para distintos números de coecientes non nulos da dirección discriminante. Pódese ver unha gráca cos resultados sobre o conxunto de testeo no seu artigo, onde os mellores foron obtidos usando a regra do centroide máis cercano con entre 10 e 20 xens signicativos, con erros de clasicación preto do 15%. Considerando tódalas variantes do sparse rLDA 5.2. DATOS DE MEDIDAS DE EXPRESIÓN XÉNICA 59 mencionadas, o erro mantense case sempre por debaixo do 20%. Para ilustrar o comportamento da regra de Bouveyron et al. (2007) sobre un conxunto de datos real, procederemos dun xeito parecido ao que se acaba de explicar empregando a implementación do HDDA feita para este traballo. Tamén se repetirá o proceso 50 veces para achar o erro promedio, e de novo outras 20 para ver a desviación deste. Tomaranse 50 das 2000 variables ao chou cada vez, para poder levar a cabo os cálculos nun ordenador de sobremesa, e empregaranse tanto HDDAi como HDDAh para dimensións intrínsecas entre 1 e 20. A selección de variables previa non pretende quedarse ca información signicativa e rexeitar a redundante, polo que non podemos esperar que os resultados que se obteñan sexan moi bos, será máis ben un exercicio computacional para ilustrar o método. Figura 5.2: Erro de clasicación medio para as variantes do HDDA sobre os conxuntos de entrenamento e testeo sacados dos datos Colon, fronte aos distintos valores considerados para a dimensión intrínseca dν . 5 10 15 20 0.0 0.1 0.2 0.3 0.4 0.5 66 APÉNDICE A. SCRIPTS UTILIZADOS PARA OS EXEMPLOS E IMPLEMENTACIÓN DO HDDA # Media e desviación típica dos erros nas 100 repeticións erros<-list(efN1,efN2,efP1,efP2,enN1,enN2,enP1,enP2,eqN1,eqN2,eqP1,eqP2) em<-rapply(erros,mean) edt<-rapply(erros,sd) Na ilustración da regra dos k veciños máis cercanos, a gráca da Figura 2.3 obtense co seguinte código. A función knn é da librería class e é necesario xar unha semente porque cando hai empate nas votacións entre os k veciños máis cercanos, a asignación a unha clase ou outra faina ao azar. set.seed(1997) ek<-c() erros<-list() for(j in 1:50){ rX<-c() for(i in 1:150){ rX<-c(rX,knn(iris[-i,1:4],iris[i,1:4],iris[-i,5],k=j)) } ek<-c(ek,sum(rX!=as.integer(iris[,5]))) erros[[j]]<-which(rX!=as.integer(iris[,5])) } plot(1:50,ek,xlab= '' ,ylab= '' ,type= ' o ' ) A representación das distintas funcións kernel da Figura 2.4 é resultado destas liñas. par(mfrow=c(1,3)) x<-seq(-3,3,by=0.01) plot(x,exp(-x^2),main= ' Kernel gaussiano ' ,type= ' l ' ,ylim=c(0,1),xlab= '' ,ylab= '' ) plot(x,1/(1+x^2),main= ' Kernel de Cauchy ' ,type= ' l ' ,ylim=c(0,1),xlab= '' ,ylab= '' ) plot(x, (1-x^2)*(abs(x)<=1),main= ' Kernel de Epanechnikov ' ,type= ' l ' ,ylim=c(0,1),xlab= '' ,ylab= '' ) E a táboa do exemplo dos métodos kernel obtense co seguinte código. mu1<-c(0,0,0) mu2<-c(2,0,-0.5) E1<-cbind(c(1,0,0),c(0,0.25,0),c(0,0,0.25)) A.1. EXEMPLOS 67 E2<-cbind(c(0.5,0.5,0.125),c(0.5,1,0),c(0.125,0,0.25)) set.seed(1997) N1<-rbind(mvrnorm(200,mu1,E1),mvrnorm(300,mu2,E1)) N2<-rbind(mvrnorm(200,mu1,E1),mvrnorm(300,mu2,E2)) P1<-rbind(cbind(rpois(200,10)+5,rpois(200,10)-5,rpois(200,10)+2), cbind(rpois(300,10),rpois(300,10),rpois(300,10))) P2<-rbind(cbind(rpois(200,10)+5,rpois(200,10)-5,rpois(200,10)+2), cbind(rpois(300,10),rpois(300,20)-10,rpois(300,30)-20)) y=factor(c(rep(1,200),rep(2,300))) N1.norm<-classif.np(y,N1,Ker =Ker.norm,metric=metric.dist) N1.epa<-classif.np(y,N1,Ker =Ker.epa,metric=metric.dist,h=seq(1.22,2,by=0.02)) N1.uni<-classif.np(y,N1,Ker =Ker.unif,metric=metric.dist,h=seq(1.22,2,by=0.02)) N1.norm$prob.classification N1.epa$prob.classification N1.uni$prob.classification N2.norm<-classif.np(y,N2,Ker =Ker.norm) N2.epa<-classif.np(y,N2,Ker =Ker.epa,metric=metric.dist,h=seq(1.1,2,by=0.02)) N2.uni<-classif.np(y,N2,Ker =Ker.unif,metric=metric.dist,h=seq(1.1,2,by=0.02)) N2.norm$prob.classification N2.epa$prob.classification N2.uni$prob.classification P1.norm<-classif.np(y,P1,Ker =Ker.norm) P1.epa<-classif.np(y,P1,Ker =Ker.epa,metric=metric.dist,h=seq(6,10,by=0.02)) P1.uni<-classif.np(y,P1,Ker =Ker.unif,metric=metric.dist,h=seq(7,10,by=0.02)) P1.norm$prob.classification P1.epa$prob.classification P1.uni$prob.classification P2.norm<-classif.np(y,P2,Ker =Ker.norm) P2.epa<-classif.np(y,P2,Ker =Ker.epa,metric=metric.dist,h=seq(8,10,by=0.02)) P2.uni<-classif.np(y,P2,Ker =Ker.unif,metric=metric.dist,h=seq(8,10,by=0.02)) P2.norm$prob.classification P2.epa$prob.classification 68 APÉNDICE A. SCRIPTS UTILIZADOS PARA OS EXEMPLOS E IMPLEMENTACIÓN DO HDDA P2.uni$prob.classification A Figura 4.2 é resultado do que segue. set.seed(1997) mu1<-c(0,0); mu2<-c(0,2) E<-matrix(c(2,1,1,0.75),2,2) N<-rbind(mvrnorm(100,mu1,E),mvrnorm(100,mu2,E)) PC1<-eigen(E)$vector[,1] prox<-N[,1]*PC1[1]+N[,2]*PC1[2] par(mfrow=c(1,2)) plot(N,col=c(rep(1,100),rep(2,100)),xlab= '' ,ylab= '' ) plot(prox,0.03*rnorm(200),xlab= '' ,ylab= '' ,col=c(rep(1,100),rep(2,100)),ylim=c(-1,1)) A.2. Implementación do HDDA e aplicación a datos simulados e reais As seguintes liñas denen unha función de R que calcula as clases ás que se asigna unha mostra de entrenamento e outra de testeo cos métodos HDDAi e HDDAh a partir da primeira, e tamén os erros de clasicación correspondentes. Esta foi a función empregada para obter os resultados mostrados no último capítulo. hdda<-function(train,test,ctrain,ctest,k,d,dnu){ alfa<-sigma<-c() mu<-S<-list() lambda<-Q<-list() classi<-classh<-ei<-eh<-list() classi$train<-classi$test<-classh$train<-classh$test<-1 # Cálculo dos estimadores dos parámetros para cada unha das k clases for(i in 1:k){ ni<-sum(ctrain==i) mu[[i]]<-colMeans(train[ctrain==i,]) S[[i]]<-(ni-1)/ni*var(train[ctrain==i,]) A.2. IMPLEMENTACIÓN DO HDDA E APLICACIÓN A DATOS SIMULADOS E REAIS 69 #Autovalores e autovectores de S lambda[[i]]<-eigen(S[[i]])$values Qnu<-eigen(S[[i]])$vectors # Matriz das proxeccións Q[[i]]<-cbind(Qnu[,1:dnu],matrix(0,d,d-dnu)) #Parámetros do HDDA a<-sum(lambda[[i]][1:dnu])/dnu b<-sum(lambda[[i]][(dnu+1):d])/(d-dnu) alfa[i]<-b/(a+b) sigma[i]<-a*b/(a+b) } alfa<-mean(alfa) #Clasificación da mostra de entrenamento con HDDAi e HDDAh HDDAi<-HDDAh<-c() for(j in 1:length(ctrain)){ for(i in 1:k){ prox<-Q[[i]]%*%t(Q[[i]])%*%matrix(as.matrix(train[j,]-mu[[i]]),d,1) +mu[[i]] HDDAi[i]<-alfa*sum((mu[[i]]- prox)^2)+(1-alfa)*sum((train[j,]-prox)^2) HDDAh[i]<-1/sigma[i]*HDDAi[i]+2*d*log(sigma[i]) } classi$train[j]<-which(HDDAi==min(HDDAi)) classh$train[j]<-which(HDDAh==min(HDDAh)) } ei$train<-sum(classi$train!=ctrain)/length(ctrain) eh$train<-sum(classh$train!=ctrain)/length(ctrain) #Clasificación da mostra de testeo con HDDAi e HDDAh HDDAi<-HDDAh<-c() for(j in 1:length(ctest)){ for(i in 1:k){ prox<-Q[[i]]%*%t(Q[[i]])%*%matrix(as.matrix(test[j,]-mu[[i]]),d,1)+mu[[i]] HDDAi[i]<-alfa*sum((mu[[i]]- prox)^2)+(1-alfa)*sum((test[j,]-prox)^2) HDDAh[i]<-1/sigma[i]*HDDAi[i]+2*d*log(sigma[i]) } classi$test[j]<-which(HDDAi==min(HDDAi)) 70 APÉNDICE A. SCRIPTS UTILIZADOS PARA OS EXEMPLOS E IMPLEMENTACIÓN DO HDDA classh$test[j]<-which(HDDAh==min(HDDAh)) } ei$test<-sum(classi$test!=ctest)/length(ctest) eh$test<-sum(classh$test!=ctest)/length(ctest) return(list(classi=classi,classh=classh,ei=ei,eh=eh)) } A continuación móstrase a simulación dos datos e os cálculos feitos sobre eles. Están tamén as liñas que xeran a Figura 5.1. set.seed(1997) dnu<-1:20 d<-100; train<-25; test<-100 mu1<-c(0,0.9); mu2<-c(0,-0.9); E<-matrix(c(1,0.7,0.7,1),2,2) cero<-rep(0,d-2); I<-diag(d-2) ctrain<-c(rep(1,train),rep(2,train)); ctest<-c(rep(1,test),rep(2,test)) err_test_l<-err_train_l<-err_test_h<-err_train_h<-err_test_i<-err_train_i<-matrix(0,50,20) for(veces in 1:50){ Xtrain<-rbind(cbind( mvrnorm(train,mu1,E), mvrnorm(train,cero,I)), cbind( mvrnorm(train,mu2,E), mvrnorm(train,cero,I))) Xtest<-rbind(cbind( mvrnorm(test,mu1,E), mvrnorm(test,cero,I)), cbind( mvrnorm(test,mu2,E), mvrnorm(test,cero,I))) for(j in dnu){ res<-hdda(Xtrain,Xtest,ctrain,ctest,2,d,j) #Almacenamento dos erros para promediar err_train_i[veces,j]<-res$ei$train err_train_h[veces,j]<-res$eh$train err_test_i[veces,j]<-res$ei$test err_test_h[veces,j]<-res$eh$test #Erros da regra linear l<-lda(Xtrain,ctrain) err_train_l[veces,j]<-sum(predict(l)$class!=ctrain)/(2*train) err_test_l[veces,j]<-sum(predict(l,Xtest)$class!=ctest)/(2*test) } } plot(dnu,colMeans(err_test_i),col=3,type= ' o ' ,ylim=c(0,0.5),xlab= '' ,ylab= '' ) A.2. IMPLEMENTACIÓN DO HDDA E APLICACIÓN A DATOS SIMULADOS E REAIS 71 lines(dnu,colMeans(err_train_i),col=1,type= ' o ' ) lines(dnu,colMeans(err_train_h),col=2,type= ' o ' ) lines(dnu,colMeans(err_test_h),col=4,type= ' o ' ) abline(h=colMeans(err_train_l),col=5) abline(h=colMeans(err_test_l),col=6) Por último, o código empregado para tratar os datos de expresión xénica de Alon et al. (1999), que están dispoñibles no paquete dprep , usando o método HDDA de Bouveyron et al. (2007). data(colon) dnu<-1:20 HDDAi<-HDDAh<-list() HDDAi$meantrain<-HDDAi$meantest<-HDDAi$sdtrain<-HDDAi$sdtest<-1 HDDAh$meantrain<-HDDAh$meantest<-HDDAh$sdtrain<-HDDAh$sdtest<-1 for(k in dnu){ # Clasificación tomando menos variables, repetida moitas veces eHDDAitrain<-eHDDAitest<-eHDDAhtrain<-eHDDAhtest<-matrix(0,50,20) for(j in 1:20){ for(i in 1:50){ # Subconxuntos de 50 expresións xénicas extraídos ao azar datos<-colon[,c(sample(1:2000,50),2001)] # 2/3 dos datos para entrenamento e 1/3 para testeo part<-sample(1:62,20) train<-datos[-part,1:50]; ctrain<-datos[-part,51] test<-datos[part,1:50]; ctest<-datos[part,51] # Chamada á función e resultados de aplicar HDDAi e HDDAh res<-hdda(train,test,ctrain,ctest,2,50,dnu) eHDDAitrain[i,j]<-res$ei$train eHDDAitest[i,j]<-res$ei$test eHDDAhtrain[i,j]<-res$eh$train eHDDAhtest[i,j]<-res$eh$test } } # Almacenamento dos resultados para cada dimensión intrínseca considerada eHDDAitrain<-colMeans(eHDDAitrain) 72 APÉNDICE A. SCRIPTS UTILIZADOS PARA OS EXEMPLOS E IMPLEMENTACIÓN DO HDDA HDDAi$meantrain[k]<-mean(eHDDAitrain); HDDAi$sdtrain[k]<-sd(eHDDAitrain) eHDDAitest<-colMeans(eHDDAitest) HDDAi$meantest[k]<-mean(eHDDAitest); HDDAi$sdtest[k]<-sd(eHDDAitest) eHDDAhtrain<-colMeans(eHDDAhtrain) HDDAh$meantrain[k]<-mean(eHDDAhtrain); HDDAh$sdtrain[k]<-sd(eHDDAhtrain) eHDDAhtest<-colMeans(eHDDAhtest) HDDAh$meantest[k]<-mean(eHDDAhtest); HDDAh$sdtest[k]<-sd(eHDDAhtest) } # Gráfica do erro promedio cometido sobre as mostras usando HDDA plot(dnu,HDDAi$meantrain,col=1,type= ' o ' ,ylim=c(0,0.5),xlab= '' ,ylab= '' ) lines(dnu,HDDAi$meantest,col=2,type= ' o ' ) lines(dnu,HDDAh$meantrain,col=3,type= ' o ' ) lines(dnu,HDDAh$meantest,col=4,type= ' o ' ) Bibliografía [1] Koch, I., Analysis of Multivariate and High-Dimensional Data , Cambridge Series in Statistical and Probabilistic Mathematics, Cambride University Press, 2014. [2] Devroye, L., Györ, L. e Lugosi, G., A Probabilistic Theory of Pattern Recognition , Applications of Mathematics, Stochastic Modelling and Applied Probability, 31, Springer- Verlag, New York, 1996. [3] James, G., Witten, D., Hastie, T. e Tibshirani, R., An Introduction to Statistical Learning , Springer Texts in Statistics, 103, Springer-Verlag, New York, 2013. [4] Fix, E. e Hodges J., Discrminatory analysis, nonparametric discrimination: Consistency properties , Tecnical report, Randolph Field, TX, USAF School of Aviation Medicine, 1951. [5] Hastie, T., Tibshirani, R. e Friedman, J., The Elements of Statistical Learning: Data Mining, Inference, and Prediction , Springer Series in Statistics, Springer-Verlag, New York, 2001. [6] Galeano, P. e Peña, D., Data science, big data and statistics , TEST, 28: 289, 2019. [7] Klamelä, J., Multivariate Nonparametric Regression and Visualization: with R and applications to nance , Wiley Series in Computational Statistics, Wiley, 2014. [8] Mardia, K., Kent, J. e Bibby, J., Multivariate Analysis , Probability and Mathematical Statistics, Academic Press, London, 1979. [9] Bouveyron, C., Girard, S. e Schmid, C., High Dimensional Discriminant Analysis , Communication in Statistics-Theory and Methods, 36, 2007. [10] Qiao, Z., Zhou, L. e Huang, J. Z., Eective Linear Discriminant Analysis for High Dimensional, Low Sample Size Data , Proceedings of the World Congress on Engineering 2008 Vol II WCE 2008, July 2 - 4, 2008, London, U.K. 73 74 BIBLIOGRAFÍA [11] Friedman, J.H., Regularized Discriminant Analysis , Journal of the American Statistical Association, 84, 1989. [12] Qin, Y., A review of quadratic discriminant analysis for high-dimensional data , Wiley Interdisciplinary Reviews: Computational Statistics. 10. e1434. 10.1002/wics.1434., 2018. [13] Tibshirani, R., Regression Shrinkage and Selection via the Lasso , Journal of the Royal Statistical Society. Series B (Methodological), vol. 58, no. 1, 1996 [14] Alon, U., Barkai, N. et al., Broad patterns of gene expression revealed by clustering analysis of tumor and normal colon tissues probed by oligonucleotide arrays , Proceedings of the National Academy of Sciences of the United States of America, 96 (12), 1999. [15] Cai, T. e Zhang, L, High dimensional linear discriminant analysis: optimality, adaptive algorithm and missing data , Journal of the Royal Statistical Society , Series B, Statistical Methodology, 2019. [16] Manuel Febrero-Bande, Manuel Oviedo de la Fuente, Statistical Computing in Functional Data Analysis: The R Package fda.usc. Journal of Statistical Software, 51(4), 1-28, 2012. [17] Venables, W. N. e Ripley, B. D. Modern Applied Statistics with S , Fourth Edition. Springer, New York, 2002.