Apunte Deteccion y Estimacion

SEÑALES Y SISTEMAS II
Apunte del Curso
Prof. Marcos Orchard
Lerko Araya H.
Francisco Casado C.
Constanza Villegas M.
13 de abril de 2020
Departamento de Ingenierı́a Eléctrica

Facultad de Ciencias Fı́sicas y Matemáticas
Universidad de Chile
Índice general
Introducción VI
¿Cómo usar este apunte? VII
1 Fundamentos de Detección 1
1.1. Introducción 1
1.2. Definiciones básicas 2
1.3. Neyman-Pearson 5
1.3.1. Test de Razón de Verosimulitud 6
1.4. Curvas ROC para determinar el mejor test 9
1.4.1. Test compuestos y múltiples hipótesis 11
1.5. Riesgo Bayesiano 12
1.5.1. Clasificación por mı́nimo riesgo 12
Problemas Resueltos 14
III
Bibliografı́a del capı́tulo 17
2 Análisis de Componentes Principales 19
2.1. Motivación 19
2.2. Análisis PCA 21
2.3. Resultados de PCA 22
2.3.1. Score Plot 22
2.3.2. Aporte de cada Componente Principal a la Variabilidad de los Datos 24
2.3.3. Loading Plot 24
2.4. Elipses de Control y Test de Hotelling 25
2.4.1. Test de Hotelling 27
Ejercicios 29
3 Estimación 33
3.2. Estimadores MSE y LS 34
3.2.1. Esperanza Condicional 35
3.2.2. Estimador Lineal de Mı́nimo Error Cuadrático Medio 37
3.2.3. Least Squares (LS) 38
3.2.4. Predictores MSE para Sistemas Lineales 38
3.3. Estimador de Máxima Verosimilitud 39
3.4. Estimador Máximo a Posteriori 41
3.5. Estadı́sticos Suficientes 42
3.6. Eficiencia de un Estimador 44
3.6.1. Score Function 44
3.6.2. Cota de Cramér-Rao 45
IV
V
4 Filtro de Kalman 55
4.2. Deducción Analı́tica 57
Agradecimientos 71
INTRODUCCIÓN
En este apunte se enfoca en los contenidos asociados a la segunda y tercera unidades del curso Señales
y Sistemas II, que imparte el Departamento de Ingenierı́a Eléctrica (DIE) de la Universidad de Chile.
Pretende ilustrar conceptos básicos y tópicos relevantes que deben tenerse en cuenta a la hora de
estudiar estos contenidos. Es menester indicar que se asume que la primera parte del curso, relativa a
modelación y análisis de sistemas, ya forma parte de la base teórica del estudiante.
El primer capı́tulo de este apunte se centra en el problema de “Detección”, partiendo del enfoque de
Neymann-Pearson y el test de razón de verosimilitud. Luego se analiza la importancia de las curvas
ROC. El problema se extiende a métodos que incorporan una caracterización explı́cita del riesgo, tales
como el enfoque Bayesiano. Posteriormente, este análisis se extiende al caso multivariable, dónde el
propósito es reducir la dimensionalidad del problema y definir la estructura de un test de hipótesis
compuesto (donde arbitrariamente se establece un caso “normal” y una hipótesis alternativa) para la
detección de anomalı́as. Ası́, se presenta el análisis de componentes principales (PCA), dando cierre
al capı́tulo.
El capı́tulo de “Estimación” trata sobre fundamentos de la minimación del error cuadrático medio,
mı́nimos cuadrados, regresión lineal, máxima verosimilitud y máximo a posteriori. Después, se pre-
senta tópicos de suficiencia y eficiencia, finalizando con la cota de Crámer-Rao.
Finalmente, en el último capı́tulo se trabaja la deducción del Filtro de Kalman, indicando nociones
básicas y su relación con el observador de Luenberger; junto a algunos ejemplos ilustrativos.
VI
¿Cómo usar este apunte?
Este apunte busca ser material de apoyo, pero en ningún caso pretende ser una guı́a completa ni
autocontenida de las materias; sólo se busca ilustrar algunos conceptos fundamentales.
Los autores quieren enfatizar que este escrito no debe considerarse como “La Biblia del curso”,
puesto que corresponde solamente a una visión particular para enseñar estos contenidos, y claramente
existen versiones alternativas de cómo tratar los tópicos.
Dicho esto, se recomienda no dejar de asistir a las cátedras, pues son ellas mismas las que han nutrido
a este apunte.
Por último, los autores invitan a dejar cualquier corrección, sugerencia, reclamo y/o felicitaciones en
el siguiente correo:
apunte-ssii@googlegroups.com.
VII
CAPÍTULO 1
FUNDAMENTOS DE DETECCIÓN
1.1. Introducción
La primera unidad del curso EL4003 apunta a entender las propiedades fundamentales de un sistema
dinámico. Ası́, mucho de la discusión se centró en cómo establecer, por ejemplo, si un sistema era
estable, observable, o controlable.
A partir de estos resultados, surge naturalmente la pregunta de cómo hacer buen uso de esta informa-
ción. Desde el punto de vista del diseño e implementación de sistemas de monitoreo y supervisión de
procesos, la detección de anomalı́as cobra especial relevancia, dado que permite identificar cuándo
el proceso se aleja del punto de operación deseado. Por lo tanto, no es extraño que el problema de
Detección tenga especial relevancia. En este problema, se deben reconocer situaciones especı́ficas,
dentro de un conjunto limitado de opciones. Una extensión natural del problema, en la que se desea
determinar el valor exacto del vector de estados del proceso, o de todos los parámetros del modelo
que inetenta caracterizarlo, corresponde al problema de Estimación.
En este capı́tulo, explicamos el procedimiento necesario para, a partir de un conjunto de observacio-
nes, aceptar o rechazar hipótesis asociadas al posible valor de un parámetro. Estas observaciones, o
mediciones, corresponden a la información adquirida a través de “sensores” (sujeta a diversas fuen-
tes de incertidumbre o ruido). La naturaleza de las mediciones está determinada por algún parámetro
θ ∈ Θ. Nuestra decisión, δ, se toma en base a la información disponible, y, y a una regla o test, φ.
1
2 FUNDAMENTOS DE DETECCIÓN
Acertar o errar en la decisión tomada tiene un costo asociado, por lo que se define una función de
mérito o de costos (L). Estos costos también están relacionados con el tiempo empleado en generar δ.
Gráficamente, el proceso completo sigue la siguiente secuencia:
Medición Toma de decisión

θ ∈ Θ −−−−−−→ y ∈ Y −−−−−−−−→ φ(y) = δ ∈ ∆; L(δ, θ)
Otros problemas que se pueden abordar desde el punto de vista de la detección son aquellos en que
se quiere clasificar o categorizar el punto de operación de un sistema. La diferencia fundamental con
los problemas de estimación, es que en estos últimos se busca aproximar el valor especı́fico de una
variable. Este enfoque se abordará en el Capı́tulo 3.
Resumiendo la notación a considerar en problemas de detección, se tiene:
Θ ⊆ Rn Conjunto de todos los posibles valores de θ,

Y Conjunto de posibles observaciones y,
∆ Espacio de todas las decisiones posibles,
φ : Y 7→ ∆ Test que permite tomar una decisión δ a partir de la observación y; i.e., δ = φ(y),
L : Θ × ∆ 7→ R Costos asociados a una decisión δ dado un valor de θ; i.e., L = L(θ, δ).
1.2. Definiciones básicas
Para entender un poco mejor los conceptos introducidos en la sección anterior, abordaremos el si-
guiente ejemplo:
EJEMPLO 1.1
Se obtienen los siguientes valores al medir sucesivamente la caı́da de voltaje ([V]) en una resis-
tencia: Y = {4,1; 4,3; 4,2; . . .}. Se desea determinar si el voltaje medido en la resistencia tiene
relación con una operación “normal”del circuito.
En este ejemplo, la decisión (δ) se relaciona con “aceptar” o “rechazar” la hipótesis de “normalidad”
asociada a las mediciones (y) que entrega el voltı́metro. Se asume que la condición de normalidad del
circuito se asocia a un valor especı́fico de un parámetro (θ). Consideremos un caso en que queremos
detectar un modo anormal de operación especı́fico del circuito. Esto quiere decir que nuestro espa-
cio ∆ consiste de dos elementos {δ0 , δ1 }. La decisión δ0 asume que el circuito opera normalmente
(hipótesis H0 ), mientras que δ1 considera que hay un problema (hipótesis H1 ).
DEFINICIONES BÁSICAS 3
En otras palabras, al elegir δ0 se acepta la hipótesis H0 (rechazando la hipótesis H1 ), y al elegir la

decisión δ1 se acepta la hipótesis H1 (rechazando la hipótesis H0 ). Usualmente a H0 se le llama
hipótesis nula y a H1 , hipótesis alternativa.
Definición 1.1. Un test para θ ∈ Θi con i = 0, .., k −1 se dice simple si cada Θi consiste exactamente
en un solo elemento. Si algún Θi tiene más de un elemento, el test se dice compuesto.
Dado que en este caso el espacio de decisión consiste en dos puntos, es decir, ∆ = {δ0 , δ1 } (lo que
implica aceptar H0 o H1 , respectivamente), el ejemplo trata un test simple de decisión binaria.
(
H0 : Voltaje = 4,1 [V]
Test de hipótesis simple:
(
Test de hipótesis compuesto:
H1 : Voltaje>4,1 [V]
De lo anterior se desprenden cuatro casos posibles:
Caso 1: elegir δ0 cuando H0 es verdadera (decisión correcta)

Caso 2: elegir δ1 cuando H1 es verdadera (detección correcta)
Caso 3: elegir δ0 cuando H1 es verdadera (fallo en detección)
Caso 4: elegir δ1 cuando H0 es verdadera (falsa alarma)
Considere el caso en que cada muestra y(t), asociada a un instante especı́fico de tiempo t, está con-
taminada por ruido n(t); provocando diferencias entre el valor real del voltaje (constante), θ, y la
medición del instrumento. Matemáticamente, este caso puede caracterizarse de la siguiente forma:
y(t) = θ + n(t),
lo que quiere decir que cada valor de y(t) tiene una distribución de probabilidad que depende directa-
mente de la distribución del ruido de medición.
Si asumimos que n(t) es variable aleatoria independiente, idénticamente distribuida (i.i.d.), tal que
n(t) ; N (0, σ 2 ), entonces, la esperanza y varianza de y(t) están dadas por:
E(y) = θ, V(y) = σ 2 ,
Ahora bien, la distribución de probabilidad de cada medición depende fuertemente de cuál hipóte-
sis es realmente correcta; es decir, depende de la naturaleza. Para entender el análisis posterior es
fundamental recordar que θ es un parámetro fijo, cuyo valor desconocemos.
Gráficamente, el caso del Ejemplo 1.1, las distribuciones asociadas a una muestra y(t), condicionales
al valor real del parámetro, se ven como en la Figura 1.1.
fy|θ=4.1 fy|θ=4.5
(H0 ) (H1 )
4.1 4.5
Figura 1.1: Distribuciones de probabilidad asociadas a cada muestra, condicional al valor del paráme-
tro (Ejemplo 1.1.)
Si ahora se escoge un umbral para decidir cuál hipótesis es la correcta, en función de una medición
especı́fica y(t), se genera inmediatamente la probabilidad de incurrir en dos tipos de errores.
Definición 1.2. La probabilidad de rechazar la hipotesis nula H0 cuando realmente es verdadera, se

conoce como error de tipo I y se denota α. A esta probabilidad también se le conoce como el tamaño
del test, o probabilidad de falsa alarma. Corresponde al Caso 4 del Ejemplo 1.1.
Definición 1.3. La probabilidad de aceptar la hipótesis H1 cuando es verdadera, se conoce como la

probabilidad de detección y se denota β. Corresponde al Caso 2 del Ejemplo 1.1.
Definición 1.4. La probabilidad de aceptar la hipotesis nula H0 cuando es H1 es verdadera, se conoce

como error de tipo II, equivalente a 1 − β. Corresponde al Caso 3 del Ejemplo 1.1.
Si agregamos estas definiciones a la Figura 1.1, se tiene:
umbral
fy|θ=4.1 fy|θ=4.5
(H0 ) (H1 )
1−β α
4.1 4.5
Figura 1.2: Errores de tipo I y II para el test simple del ejemplo 1.1.
En la práctica, el umbral se define en base a un α dado. Generalmente se usan valores del orden del
5 %.
No hay que olvidar que una decisión depende fuertemente de cómo procese los datos que recibo.
Por ejemplo, podrı́a considerarse para el análisis tan sólo el primer valor de las muestras adquiridas.
También podrı́a elegirse tomar una decisión en base a la última muestra. o cualquier otra escogida al
azar. Sin embargo, sabiendo que el ruido distribuye Gaussiano, el promedio de las muestras resulta
NEYMAN-PEARSON 5
ser una de las mejores formas de procesar datos. En efecto, considérese

n
1X
y= y[i],
n i=1
la distribución de probabilidad de y está dada por:
σ2

y ; N 0, .
n
Por lo tanto, si quisiera tomar una decisión en base a y, encontrarı́a que las distribuciones de la Figura
1.1 estarı́an más concentradas alrededor de θ para un elevado número de datos, tal como lo muestra
la Figura 1.3.
umbral
fy|θ=4.1 fy|θ=4.5
(H0 ) (H1 )
1−β α
Figura 1.3: Distribuciones de y según las hipótesis del ejemplo 1.1
Este fenómeno hace que mi decisión tenga un grado de confianza más elevado. De hecho, α es mucho
más pequeño y β se acerca a 1, lo que es bastante conveniente.
1.3. Neyman-Pearson
El enfoque de Neyman-Pearson busca generar un marco de toma de decisiones donde se logre maxi-
mizar la probabilidad de detección, simultáneamente acotando la probabilidad de falsa alarma.
Definición 1.5. Un test φ se dice ser el mejor, de tamaño α, para H0 vs H1 si:
Eθ0 [φ(y)] = α,
∀φ0 (y) tal que Eθ0 [φ0 (y)] ≤ α se tiene que
β = Eθ1 [φ(y)] ≥ Eθ1 [φ0 (y)] = β 0

Teorema 1.1. Lema de Neyman-Pearson
Sea Θ = {θ0 , θ1 } y X v.a. con distribución de probabilidad fX (x|θ). Sea ν > 0 un umbral.
i Sea 0 ≤ γ ≤ 1, todo test de la forma:


1 Si fX (x|θ1 ) > νfX (x|θ0 )
φ(x) = γ Si fX (x|θ1 ) = νfX (x|θ0 )
0 Si fX (x|θ1 ) < νfX (x|θ0 )

es el mejor en su tamaño para probar H0 : θ = θ0 vs. H1 : θ = θ1 .

ii Para 0 ≤ α ≤ 1 existe un test φ(x) de la forma anterior, con γ constante para el cual
Eθ0 [φ(x)] = α.
iii Si φ0 (x) es un mejor test de tamaño α para H0 vs. H1 , entonces tiene la forma anterior,
excepto tal vez para algún set X con probabilidad nula para H0 y H1 .
1.3.1. Test de Razón de Verosimulitud
Definición 1.6. Función de razón de verosimilitud

Se define la función de razón de verosimilitud `(x), como sigue:
fX (x|θ1 )
`(x) :=
fX (x|θ0 )
Es importante observar que esta función se inspira en la propuesta del Teorema de Neyman-Pearson.
Teorema 1.2. Test de Razón de Verosimilitud (Likelihood Ratio)
Se decide la hipótesis H1 por sobre H0 si:
fX (x|θ1 )
`(x) = >ν
fX (x|θ0 )
Dado que la mayorı́a de las distribuciones utilizadas en la práctica pertenecen a la familia de expo-
nenciales, resulta útil tomar el logaritmo de `(x). Notemos que no se afecta la desigualdad, ya que el
logaritmo es una función continua y estrictamente creciente.
Definición 1.7. Función de Log-Verosimilitud
La función de log-verosimilitud se define como el logaritmo de la función de verosimiliud, es decir:

fX (x|θ1 )
Λ(x) := log
fX (x|θ0 )
De la definición anterior, y el Teorema 1.2, se obtiene inmediatamente:

NEYMAN-PEARSON 7
Corolario 1.1. Test de Log-Verosimilitud

Se decide H1 por sobre H0 si:
Λ(x) = log (fX (x|θ1 )) − log (fX (x|θ0 )) > log(ν) = ν ∗
Para ilustrar la utilidad del test de razón de verosimilitud, revisemos el caso en que se quiere decidir
respecto del valor de un parámetro constante θ en un sistema que mide una señal x(t) = θ + n(t). Si
tratamos de tomar una decisión basados solamente en la primera medición, se tiene como única fuente
de información a x1 = x(1) = θ + n(1). Además, suponiendo que el ruido n(t) is i.i.d. y Gaussiano,
se puede definir nuestro test de hipótesis como sigue:
−(x1 −θ0 )2

 H : θ = θ , f (x |θ ) =
0 0 x1 1 0 √1 e 2σ 2
2πσ
−(x1 −θ1 )2
.
 H : θ = θ , f (x |θ ) =
1 1 x1 1 1 √1 e 2σ 2
2πσ
Aplicando el concepto del test de log-verosimilitud:
Λ(x) > ν ∗
1 1
− 2
(x1 − θ1 )2 + 2 (x1 − θ0 )2 > ν ∗
2σ 2σ
2x1 θ1 − 2x1 θ0 + θ02 − θ12 > ε(ν ∗ )
x1 > δ(ν ∗ ),
En los cálculos anteriores, se omitieron pasos algebraicos para mostrar que la importancia del test se
reduce a comparar el valor de la medición con un cierto umbral δ adecuado.
Lo interesante de este resultado es que si se tiene un conjunto de mediciones, suponiendo que el ruido
n(t) es i.i.d.1 , el resultado converge al promedio de las mediciones. Partiendo por el caso simple en
que se agrega una segunda muestra x2 = x(2) = θ + n(2), las densidades de probabilidad relevantes
pasan a ser las conjuntas dadas por:
1 2 −12 [(x1 −θ1 )2 −(x2 −θ1 )2 ]
fx1 ,x2 (x1 , x2 |θ1 ) = ( √ ) e 2σ ,
2πσ
1 2 −12 [(x1 −θ0 )2 −(x2 −θ0 )2 ]
fx1 ,x2 (x1 , x2 |θ0 ) = ( √ ) e 2σ ,
2πσ
1 Independiente e idénticamente distribuido.
Calculando la razón de log-verosimilitud entre ambas distribuciones, la decisión depende de

2
X 2
X
2
− (xi −θ1 ) + (xi − θ0 )2 > ε,
i=1 i=1
2
X
xi > δ(θ0 , θ1 , σ),
i=1
lo que para n muestras se extrapola a
n n
1 X 1 X
Λ(x) = − 2 2
(xi − θ1 ) + 2 (xi − θ0 )2 > ν ∗ (1.1)
2σ i=1 2σ i=1
n
1X
xi > δ
n i=1
De la Ecuación 1.1 y la hipótesis i.i.d. para n(t), resulta interesante notar que Λ(x) = γ ni=1 x(i) +
P
ε, lo que hace evidente que Λ(x) es una V.A. con distribución Gaussiana. Su media, sin embargo,
depende de cuál hipotesis es correcta, según lo decida la naturaleza:
(
Si H0 es correcta ⇒ Λ(x) ; N (θ0 γn + ε, γ 2 σ 2 n)
Si H1 es correcta ⇒ Λ(x) ; N (θ1 γn + ε, γ 2 σ 2 n)
Con lo anterior, y denominando E[Λ(x)] = µ, se concluye que nuestra decisión se basa en un umbral
η tal que:
H0 √
Λ(x) ≶ η = µ + γ nσz1−α , (1.2)
H1
donde z1−α representa el valor que acumula una probabilidad de 1 − α en la distribución N (0, 1)
(estos valores se pueden encontrar en tablas). Podemos, entonces, fijar dependencias entre el valor de
η y el error de tipo I:
Z ∞
α= fΛ(x) (Λ(x)|H0 )dΛ(x). (1.3)
η(α)
La Figura 1.4 ilustra esta probabilidad α para una distribución cualquiera (no necesariamente Gaus-
siana):
CURVAS ROC PARA DETERMINAR EL MEJOR TEST 9
fΛ(x)|Hi umbral de decisión
Figura 1.4: Distribución de probabilidades para Λ(x).
El valor de η se puede calcular según la Ecuación (1.2) en el caso Gaussiano. Una vez encontrado η,
se puede calcular la probabilidad de detección:
Z ∞
β= fΛ(x) (Λ(x)|H1 )dΛ(x). (1.4)
η
1.4. Curvas ROC para determinar el mejor test
Ahora la pregunta es cómo decidir cuál test de hipótesis es mejor. Para esto se estudia la relación
entre α y β mediante las denominadas Curvas ROC (Receiver Operating Characteristic). Estas curvas
suelen lucir de la siguiente forma:
curva ideal
1
β
equiprobabilidad
0
0 1
α
Figura 1.5: Curvas ROC.
Conceptualmente, un test “perfecto”serı́a aquel que, para una tasa de falsas alarmas arbitraria, siempre
logre detectar situaciones en las que la hipótesis nula es incorrecta, es decir, ∀α, β = 1. Este caso
corresponde a la linea horizontal que pasa por β = 1. Por otra parte, la recta que pasa por el origen y
por el punto (1;1) corresponde al peor escenario aceptable ya que muestra que una falsa alarma tiene
la misma probabilidad que una detección, por lo que se puede cambiar este test por uno más sencillo
como tomar la decisión en base al resultado de lanzar una moneda. Luego, cualquier curva que se
sitúe por debajo de esta recta tendrá mayor tasa de menor detección que de falsas alarmas.
Tomando en cuenta esto, se puede ver que las curvas ROC siempre estarán contenidas en el triángulo
ubicado por sobre la diagonal mencionada. Es por esto que tı́picamente las curvas ROC se ven como
ilustra la Figura 1.6
La elección del mejor test se realiza considerando que para una tasa de falsas alarmas α dada, es mejor
aquel que entrega una probabilidad de detección β mayor. En el caso de la Figura 1.6 es evidente que
que la Curva 1 siempre es mejor que las otras.
(1)
(3)
(2)
β
0
0 α0 1
α
Figura 1.6: Comparación de curvas ROC.
Sin embargo, es interesante revisar casos como los ilustrados por las Curvas (2) y (3) en las cercanı́as
de α = α0 . En estos casos, resulta provechoso sacrificar un poco la tasa de falsas alarmas, es decir,
aumentar α, para lograr un aumento de la probabilidad de detección. Para la Figura 1.6, esto significa
pasar de la Curva (2) a la Curva (3).
EJEMPLO 1.2
Un ejemplo en que se realiza este trade-off es el caso en que se quiere estimar el estado de salud
(SoH, por sus siglas en inglés) de baterı́as de Ion-Litio. Estas baterı́as suelen experimentar un
fenómeno denominado regeneración de capacidad, en que la baterı́a pareciera recuperar su salud
por algunos ciclos de carga-descarga (ver Fig. 1.7). Este fenómeno hace que muchas veces se
estime una capacidad superior a la real, y por tanto, una vida útil mayor a la efectiva. En casos
como este, las falsas alarmas no son crı́ticas, dado que por lo general se tiene otro algoritmo más
sofisticado que, en una etapa posterior del proceso, permite filtrar los falsos negativos generados.
CURVAS ROC PARA DETERMINAR EL MEJOR TEST 11
Degradación de una Batería

15
Regeneración
Capacidad [A-h]
10
-5
0 1 2 3 4 5
Tiempo [Ciclos]
Figura 1.7: Regeneración de capacidad en una baterı́a de Ion-Litio.
1.4.1. Test compuestos y múltiples hipótesis
Generalmente, los test no son en base a hipótesis simples. La complicación en estos casos es que no
se puede determinar la probabilidad de detección, pues no existe una única hipótesis alternativa a la
nula, ¡sino infinitas!
Es por esto que suelen construirse las curvas ROC de manera empı́rica para datos con una distribución
“conocida”. Por ejemplo, para el siguiente test:
(
H0 : θ = θ0
H1 : θ =
6 θ0
lo que se hace es “simular” datos y generar empı́ricamente la distribución asociada a la hipótesis

nula. Luego, aproximando dicha distribución a situaciones similares a la ilustrada en la Figura 1.8, se
obtiene la frontera entre H0 y H1 para un α dado:
H0
H1 H1
α α
2 2
θ0
Figura 1.8: Distribución supuesta para los datos y obtención de la frontera entre las hipótesis.
Esto mismo, llevado al caso de múltiples hipótesis puede llevar a situaciones como la ilustrada en
la Fig. 1.9. En este caso, algoritmos tales como “Random Forest” o estructuras tipo caja negra tales
como “Redes Neuronales” pueden ser de utilidad.
H3
H1
H0
H2
H4
Figura 1.9: Caso multivariable, multiples hipótesis.
1.5. Riesgo Bayesiano
En el enfoque de Riesgo Bayesiano, cada decisión está asociada a un costo, y se interesa minimizar
el valor esperado del costo total. Este enfoque se diferencia de los métodos anteriormente estudiados
por cuanto supone la existencia de un conocimiento previo respecto del valor que el parámetro θi ∀i
podrı́a tener. Dicho conocimiento se representa a través de la probabilidad “a priori” del vector de
parámetros; supuesto que implica que θi es ahora una variable aleatoria.
Se define x ∈ Ci , con i = 1, . . . , n las múltiples hipótesis del test y Lij el costo asociado a elegir la
clase j cuando realmente x ∈ ωi . 2
Dicho esto, nos interesa estudiar la probabilidad P(ωi |x), es decir, la probabilidad de estar en la
clase ωi habiendo medido un valor x. Será útil, entonces, el Teorema de Bayes, que se enuncia a
continuación:
P(B|A) · P(A)
P(A|B) = . (1.5)
P(B)
1.5.1. Clasificación por mı́nimo riesgo
Definición 1.8. Sean Lij los costos ya definidos, el riesgo rj asociado a elegir la clase j está dado por
C
X
rj = Lij P(ωi |x). (1.6)
i=1
2 Notar que Lii muestra el costo de acertar, que no necesariamente es igual a cero.
RIESGO BAYESIANO 13
Generalmente, calcular P(ωi |x) (elegir la clase ωi habiendo medido x) es difı́cil, sin embargo, el
Teorema de Bayes resulta ser útil para salir de este entuerto, por cuanto dicha probabilidad puede
expresarse como:
P(x|ωi ) · P(ωi )
P(ωi |x) = . (1.7)
P(x)
Ası́, el riesgo total asociado a decidir “j”queda expresado de la siguiente forma
C
X P(x|ωi ) · P(ωi )
rj = Lij , (1.8)
i=1
P(x)
donde es notorio que el término P(x) es común a todos los términos de la suma. Por otra parte, las
decisiones son de la forma: si r1 < r2 se elige 1; de lo contrario, se elige 2. Para el caso de de dos
variables esto se denota:
1
r1 ≶ r2 .
2
Luego, como P(x) ≥ 0, el valor de esta probabilidad no afecta la decisión. Esto permite comparar los
riesgos sin necesidad de calcular la probabilidad de recibir la información “X”.
EJEMPLO 1.3
Sean r̂1 , r̂2 los riesgos dados por:
r̂1 =L11 P(x|ω1 ) · P(ω1 ) + L21 P(x|ω2 ) · P(ω2 )

r̂2 =L12 P(x|ω1 ) · P(ω1 ) + L22 P(x|ω2 ) · P(ω2 )
Si se hace un supuesto fuerte en que no equivocarse tiene costo cero (L11 = L22 = 0), se tiene
que el criterio de decisión es
1
L21 P(x|ω2 )P(ω2 ) ≶ L12 P(x|ω1 )P(ω1 )
2
P(x|ω2 ) 1 L12 P(ω1 )

≶ (1.9)
P(x|ω1 ) 2 L21 P(ω2 )
Este resultado es importante, pues exhibe el mismo criterio de decisión de la razón de verosimilitud.
Ası́, el enfoque Bayesiano resulta ser más general y eficiente para tomar decisiones.
Problemas Resueltos
EJEMPLO 1.4
Considere un problema de detección donde se quiere establecer cuál es el mejor modelo para
describir una variable aleatoria Z en función de dos variables explicativas X e Y cuya distribución
conjunta es 1/(XY )2 ; X ≥ 1, Y ≥ 1. Para este efecto, se consideran dos posibles hipótesis:
H0 : Z = Y vs. H1 : Z = XY
Implemente el test de razón de verosimilitud que le permitirı́a rechazar Ho, indicando claramente
la definición de la región de aceptación/rechazo del test. Considere el cambio de variables U =
XY y V = X/Y para efectos de simplificar el cálculo de las marginales. Finalmente, encuentre
expresiones para el nivel de significancia deseado, PF A , y la probabilidad de detección PD , en
función del umbral del test, η.
Para obtener el test de verosimilitud se tiene que calcular:
fZ (z|H1 )
`(x) :=
fZ (z|H0 )
Es fácil observar que:
∞
1 1 1
Z
fZ (z|H0 ) = dx = = 2
1 x2 y 2 y 2 z
X
Para calcular la distribución fZ (z|H1 ) se deberá realizar el cambio de variable U = XY ,V = Y
.
De este modo:
1
fU,V u, v = 2 Ju,v
u
Ası́ se puede calcular el jacobiano por definición:
√ √
v u
√
2 u
√
2 √v
1
Ju,v = =

√1
2 uv
− 2√uv3 2v
Luego es importante observar que de los lı́mites antiguos de X, Y se concluye que 1 ≤ u ≤ ∞ y

que u1 ≤ v ≤ u. De esta manera, la distrución buscada es:
PROBLEMAS RESUELTOS 15
fZ (z|H1 ) = fu
Z u
1
= dv
1
u
2u2 v
log(u)
=
u2
log(z)
=
z2
Ası́ se puede calcular la razón de verosimilitud:

1
2
`(x) =  z  = log(z) (1.10)
 log(z) 
 1 
z2
Ası́ la implementación del test nos queda:
log(z) < η Decido H0

log(z) = η Indiferente
log(z) > η Decido H1
De este modo, para calcular el error de tipo I se debe calcular la función de distribución f`(x) (z|H0 ),
para esto considera el siguiente cambio de variables z = log(y), w = log(x).
1
fZ,W (z, w) = Jz,w
e2(z+w)
donde:
1
0 1
Jz,w = 1 ez
= ez+w ⇒ fZ,W (z, w) = z+w

ew 0 e
Finalmente, se puede calcular f`(x) (z|H0 ) integrando en w. Con este resultado se calcula PF A .
Z ∞
PF A = α = e−z dz = e−η
η
Del mismo modo, se procede para calcular la probabilidad de detección. Cambio de variable
z = log(xy), w = log(x/y).
1
fZ,W (z, w) = 2z Jz,w
e
z+w z−w

1e 2 1
e 2 1 1
Jz,w = 12 z+w 2 1 z−w = ez ⇒ fZ,W (z, w) = z

2e 2 −2e 2 2 2e
Z ∞
z −2η 1 η
PD = β = dz = e +
η e2z 4 2
EJEMPLO 1.5
Considere un problema de detección donde una señal cambia su valor esperado en algún momento
no (conocido, para simplificar el problema) y se requiere detectar la existencia de dicho cambio.
Para ello, se quiere comparar las siguientes hipótesis:
H0 : Xi ∼ N (m0 , σ 2 ), i = 1, 2, ..., n
Xi ∼ N (m0 , σ 2 ), i = 1, 2, ..., no − 1

H1 :
Xi ∼ N (m1 , σ 2 ), i = no , no + 1, ..., n
Demuestre que el mejor test, para un nivel de significancia dado, puede escribirse:
n
1 X
T (x) = (xi − mo ) > η
n − no + 1 i=n
o
para algún umbral η. Encuentre una expresión que relacione dicho umbral con el nivel de signifi-
cancia deseado, PF A .
Aplicando el test de razón de log-verosimilitud a estas hipótesis, se tiene:
Λ(x) > η ∗

fX (x|θ1 )
log > η∗
fX (x|θ0 )
0 −1
nX N N
(xi − m0 )2 X (xi − m1 )2 X (xi − m0 )2
2
+ 2
− 2
> η∗
i=1
2σ i=n
2σ i=1
2σ
0
N N
X (xi − m1 )2 X (xi − m0 )2
− > η∗
i=n0
2σ 2 i=n
2σ 2
0
N
X N
X
(x2i − 2xi m1 + m21 ) − (x2i − 2xi m0 + m20 ) > η ∗∗
i=n0 i=n0
N
X
2xi (m0 − m1 ) > η ∗∗∗
i=n0
N
X
xi > η ∗∗∗∗
i=n0
N
1 X
(xi − m0 ) > η
N − n0 + 1 i=n
0

Dada la estructura del test, se puede utilizar el lema de Neyman-Pearson (Teorema 1.1 en el
apunte) visto en clases, se puede concluir que este test es el mejor para su nivel de significancia.
De la demostración de existencia, se observa que E[T (x)] = 0 si H0 es cierto. Luego, asumiendo
independencia, el nivel de significancia deseado viene dado por la siguiente expresión:
∞
1
Z
− x
2σ 2
η
α= e T = 1 − erf ( √ )
η 2 2πσT
2
1
Donde σT2 = N −n0 +1
(N + n0 )σ 2 .
Bibliografı́a del capı́tulo
1. T. K. Moon, W.C. Stirling, “Introduction to Detection and Estimation and Mathematical

Notation”, in Mathematical Methods and Algorithms for Signal Processing. Prentice-Hall:
New Jersey, 2000.
2. T. K. Moon, W.C. Stirling, “Detection Theory”, in Mathematical Methods and Algorithms

for Signal Processing. Prentice-Hall: New Jersey, 2000.
Correcciones, sugerencias, reclamos y/o felicitaciones del capı́tulo, por favor enviarlas a: morchard@ing.uchile.cl.
CAPÍTULO 2
ANÁLISIS DE COMPONENTES PRINCIPALES
2.1. Motivación
Generalmente en la industria se dispone de registros históricos con gran cantidad de datos y una
igualmente extensa lista de variables de proceso. Esta cantidad de datos puede resultar difı́cil de
analizar, especialmente considerando que muchas veces sólo podemos acceder a ellos a través de
Microsoft Excel o software equivalente.
t x1 x2 x3 x4 x5 ...
1 0.128 0.062 0.360 1.801 2.085 ... 10

0
-10
20
30
2 0.221 1.223 1.693 0.798 2.236 20
3 0.485 1.323 0.756 0.014 1.476 ...

10
4 0.896 0.195 2.487 2.605 3.432
5 0.356 0.860 0.045 0.229 4.405 ... 0
6 0.541 1.563 2.174 0.198 0.223 0

.. .. .. .. .. .. .. 10
. . . . . . . 20
Tabla 2.1: Datos en una matriz Figura 2.1: Datos disponibles.
En este contexto, el Análisis de Componentes Principales o PCA (por su sigla en inglés Principal
Component Analysis) muchas veces pasa a ser la primera herramienta de análisis que puede utilizarse
para intentar extraer información útil de los datos almacenados. PCA apunta a visualizar rápidamente
19
la estructura de correlación de las variables de proceso y reducir la dimensionalidad del problema de
supervisión asociado.
Para desarrollar la intuición detrás de este análisis, primero consideremos datos de dos variables x1 , x2
en un rango tı́pico de operación, tal como muestra la Figura 2.2
x2
200
No respeta la correlación, pero está

dentro de los rangos
150
100
Fuera de los
rango de x 1
50
Fuera de los
rango de x 2
x1
0 50 100 150 200
Figura 2.2: Rangos de operación para dos variables.
Cualquier dato que se aleje en deması́a en una dirección distinta a la que define la máxima variabilidad
conjunta entre x1 y x2 es un dato que no respeta la correlación entre las variables. Como puede
apreciarse, no basta con estar dentro de los rangos de operación de cada variable. Si bien siempre es
posible realizar una análisis similar al aquı́ presentado, obviamente la complejidad computacional se
incrementa drásticamente si debemos estudiar lo que ocurre con 120 variables de proceso o más. Es
aquı́ donde nos apoyamos en conceptos de álgebra lineal, encontrando la transformación lineal que
define cada uno de los ejes de máxima variabilidad, proyectando los datos en ellos, de modo de rotar
los ejes de observación hasta hacerlos coincidir con dichas direcciones de máxima variabilidad.
Ahora bien, PCA va un poco más allá, puesto que reduce la cantidad de ejes a sólo aquellos que
presentan variabilidades relevantes en los datos, dejando fuera del análisis a aquellos que no. Por
ejemplo, si los datos están contenidos en una región del espacio de observación similar a una galleta,
como en la Figura 2.3, donde existen dos ejes de alta variabilidad y uno de baja variabilidad, se pueden
reducir la dimensión del problema a una proyección en el plano, despreciando el último eje.
Figura 2.3: Datos contenidos en un volumen similar a una galleta.
Es decir, x3 puede despreciarse en el análisis por poseer una baja variabilidad relativa a las direcciones
x1 y x2 . Esto es en esencia lo que hace PCA , con la salvedad que las direcciones x1 , x2 y x3 se ordenan
20
ANÁLISIS PCA 21
en función del porcentaje de variabilidad acumulada en cada una de ellas. Esto permite visualizar los
datos de manera más eficiente.
2.2. Análisis PCA
Formalmente, consideremos para efectos de implementar el análisis PCA que los datos de la Tabla 2.1
(normalizados) se agrupan en la matriz X̂ dada por:
 
x̂11 . . . x̂1m
X̂ =  ... . . . ..  ,

. 
x̂n1 . . . x̂nm
donde cada columna representa una variable de proceso distinta y cada fila un instante de tiempo u
observación. De este modo se tienen m variables, y n datos por variable. Nótese que la normalización
de los datos corresponde a realizar la operación
xij − xj
x̂ij = . (2.1)
σj
Con los datos normalizados, se procede a calcular la matriz de varianza-covarianza empı́rica como
1
S= X̂ T X̂ = V ΛV T . (2.2)
n−1
Cabe notar que la matriz de varianza-covarianza empı́rica tiene dimensión m × m y se aproxima con
un estimador insesgado, lo que explica el denominador n − 1. Por otro lado, el último término en
la igualdad corresponde a su diagonalización, donde V es la matriz de vectores propios ortogonales
ordenados de la siguiente forma:
|λ1 | > |λ2 | > |λ3 | > . . . > |λm |

V = [v1 |v2 |v3 | . . . |vm ]
La matriz de varianza-covarianza tiene vectores propios que son ortogonales entre sı́. Es más, dichos
vectores propios están alineados con las direcciones de máxima variabilidad de los datos. De hecho,
la variabilidad de las proyecciones de los datos en cada vector se relaciona directamente con el valor
del valor propio correspondiente. Al estar ordenados los vectores, en función de valores propios de-
crecientes, se está directamente ordenando la contribución de cada proyección a la variabilidad total
de los datos. Debido a esto, si se quiere reducir la dimensionalidad del problema, siempre se establece
una base ortogonal de proyección con las primeras p columnas de la matriz V .
Las coordenadas de cada vector propio indican la importancia de cada variable de proceso dentro de
la componente principla respectiva, información que será de utilidad al analizar los Loading Maps.
22 ANÁLISIS DE COMPONENTES PRINCIPALES
2.3. Resultados de PCA
En esta sección ilustraremos el tipo de resultados que puede obtenerse al implementar PCA , usando
el ejemplo “cities”que está incorporado en el software Matlab1 , donde se estudian los indicadores
socioeconómicos (variables de proceso) de la Figura 2.4.
economics
recreation
arts
education
transportation
crime
health
housing
climate
0 1 2 3 4 5
×10 4
Figura 2.4: Indicadores entregados por cities.
Lo primero que se debe notar es que estos datos no están normalizados, por lo que las varianzas entre
dos o más indicadores no está en los mismos órdenes de magnitud. Por ejemplo, los valores medidos
para climate no se pueden comparar directamente con los de housing. Es por esto que, como se dijo
anteriormente, es necesario normalizar los datos. Una vez normalizados, la desviación estándar de
cada variable será unitaria. Si bien el análisis PCA muestra su mayor valor en el caso de variables que
tienen distribución conjunta Gaussiana, esto no es una restricción para la implementación del análisis
empı́rico. De hecho, los datos en la Figura 2.4 distan mucho de tener dicha distribución.
2.3.1. Score Plot
El primer tipo de resultado que entrega PCA es el denominado Score Plot, un gráfico donde se mues-
tra cada dato proyectado sobre una selección (generalmente pares) de componentes principales. Por
ejemplo, en la Figura 2.5 se muestran los datos proyectados sobre las dos primeras componentes
principales.
Cada punto en el gráfico corresponde a uno de los n datos (filas) de la matriz X̂, proyectado en las
componentes ya mencionadas.
Para calcular los “scores”(proyecciones de cada punto en el score plot), se considera la matriz de
proyección P = [v1 |v2 | . . . |va ], con a < m y la siguiente transformación lineal:
1 El script completo se encuentra disponible en U-Cursos, con el nombre “EL4003.m”

4
2nd Principal Component

2
-2 Outliers
Puntos similares
-4
-5 0 5 10 15
1st Principal Component RESULTADOS DE PCA 23
Figura 2.5: Score Plot
Y = X̂P
= [y1 |y2 | . . . |ya ] ∈ Rn×a . (2.3)
Considerando que los datos en la matriz X̂ ya están normalizados, por tanto centrados en el orı́gen, la
operación anterior puede entenderse como una rotación del eje coordenado para visualizar los datos
usando las primeras “a” componentes principales.
Una interesante propiedad de las proyecciones en componentes principales es que las coordenadas
yti del vector ~yt = [yt1 |yt2 | . . . |yta ] no están correlacionadas entre sı́. Visto de otra forma, la ma-
triz de varianza-covarianza de los“scores” es diagonal. Si además ~y distribuye normal, entonces las
componentes principales son independientes.
En el “score plot”de la Figura 2.5 se han destacado algunos grupos de puntos. Hay un conjunto de
puntos que muestran un aglutinamiento importante alrededor del promedio de los datos (son similares
entre sı́), mientras que otros muestran diferencia significativas, considerándose Outliers. Esta sepa-
ración entre puntos normales y outliers se debe al hecho que si los datos tuviesen una distribución
Gaussiana, uno esperarı́a una gran acumulación alrededor de la media, y una menor probabilidad de
encontrar puntos conforme me alejo de dicho estadı́stico. Esta forma de imaginar el problema se tra-
duce generalmente en mapas donde se establecen lı́mites (umbrales) para una condición de operación
normal que tienen forma de elipses (ver Fig. 2.6).
3
2
1
v2
0
-1
3 2 -2
1 0 2 3
-1 -2 -1 0 1 -3
v2 -3 -3 -2 -3 -2 -1 0 1 2 3
v1
v1
Figura 2.6: Curvas de nivel de una distribución normal.

2.3.2. Aporte de cada Componente Principal a la Variabilidad de los Datos
El segundo resultado importante proviene de la comparación entre los aportes de cada componente
principal a la variabilidad total de los datos. Esta variabilidad puede medirse en función de la traza de
la matriz de varianza-covarianza empı́rica:
m
X
2
σ = tr(Λ) = λi , (2.4)
i=1
donde el aporte individual de la componente principal vj se puede escribir como:
λ
Pm j . (2.5)
i=1 λi
Este concepto puede visualizarse en un gráfico que muestre la acumulación de variabilidad incorpo-
rada secuencialmente al análisis por parte de cada componente principal; ver Figura 2.7:
100 100%
Variance Explained (%)
80 80%
60 60%
40 40%
Últimas componentes
20 suelen 20%
ser ruido
0 0%
1 2 3 4 5 6 7
Figura 2.7: Aporte a la varianza de cada componente principal
El gráfico anterior permite rápidamente establecer cuántas componentes principales estarı́a yo dis-
puesto a incorporar en mi análisis, equilibrando apropiadamente le complejidad computacional y la
capacidad de caracterizar la variabilidad de los datos. Aportes más pequeños a la variabilidad pueden
asumirse como parte del ruido de proceso.
2.3.3. Loading Plot
El tercer resultado importante de un análisis de componentes principales es el denominado Loading

Plot o Loading Map, que muestra la relevancia de cada variable de proceso en cada una de las compo-
nentes principales. Este gráfico puede superponerse al Score Plot, pues tiene los mismos ejes (aunque
su significado es distinto). Este gráfico es útil para estudiar la estructura de correlación entre las va-
riables de proceso.
ELIPSES DE CONTROL Y TEST DE HOTELLING 25
0.5 education
Componente de mayor
dependencia
health
arts
transportation
Component 2
0
climate
housing
crime
recreation
Componentes muy
economics
-0.5 correlacionadas
-0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5 0.6
Component 1
Figura 2.8: Loading Plot
Ahora bien, también pueden darse casos en que PCA no sea la herramienta más apropiada para efcuar
un análisis preliminar de datos. Considere por ejemplo los casos ilustrados en la Figura 2.92 , en
donde PCA no ayuda a distinguir condición normal versus una anomalı́a, a pesar de que existen datos
que están claramente segregados (puntos de operación dintingubles), pese a tener media nula y tener
direcciones de máxima variabilidad evidentes.
Cluster in cluster Corners Clusters

5 10
20
5
0 0 0
-5
-10 -20
-5
-5 0 5 -10 0 10 -20 0 20
Figura 2.9: Ejemplos en que PCA no es tan útil.
2.4. Elipses de Control y Test de Hotelling
Recordemos que las curvas de nivel de una distribución conjunta Gaussiana serı́an elipsoides centra-
das en la media. Esto motiva la utilización de curvas elı́pticas para establecer lı́mites a una operación
normal de proceso, ver Figura 2.10a y Figura 2.10b (esta última incorporando la rotación asociada a
la matriz de proyección P ). Si además, se utiliza una métrica ponderada para la distancia de un punto
2 Base de datos por: http://www.junuxx.net/

a la media (que incorpore explı́citamente la variabilidad que existe en uno u otro eje), se puede trans-
formar las elipses en circunferencias; ver Figura 2.10c. La distancia de Mahalanobis
√ pondera cada
coordenada de yti por el inverso de la raı́z de su varianza respectiva (σi = λi ) para estos efectos.
v2 v1
~xt
~xt ~xt
v2 /σ2
x̂2
x̂1 v2 v1 v1 /σ1
(a) (b) (c)
Figura 2.10: Agrupación de los datos normalizados (a), scores (b) y scores normalizados (c).
Luego, la distancia del dato xt al origen del score plot normalizado se mide de la siguiente forma:
x2t1 x2t2 x2ta

kxt k2 = + + · · · +
σ12 σ22 σa2
Lo que matricialmente se escribe como una forma cuadrática :
kxt k2 = xt P Λ−1 T T
a P xt (2.6)
La pregunta que uno debe hacerse es ¿para qué se quiere calcular la distancia de los datos al origen?
Y la respuesta es bastante simple: basta recordar que las distribuciones de los datos se asumen Gaus-
sianas y como se utilizan los datos normalizados, su media es cero, por lo que el origen del score plot
coincide con el promedio de los datos. Luego, saber la distancia de cada dato al origen es equivalente
a saber cuán alejado está dicho dato del promedio.
Ahora bien, ¿cuál es la distancia distancia (ponderada) η 0 de la media que asegura que existe una
probabilidad (1-α) de que encontrar un dato más alejado? Para responder esta consulta, recordemos
que la forma cuadrática en Ecuación (2.6) distribuye χ2 (a) (ver Fig. 2.11). Con lo que la probabilidad
de que un dato esté a menos de un cierta distancia del centro estará determinada por dicha distribución.
χ2
Figura 2.11: Densidad de probabilidad de la distribución χ2 .

ELIPSES DE CONTROL Y TEST DE HOTELLING 27
Este resultado, en que se pasa de un problema multivariable a uno escalar, facilita enormemente el
análisis. Sin embargo, no hay que confiarse. Si recordamos que en realidad se tiene un set de datos
con los que se estima la matriz de varianza-covarianza, Ec. (2.2), entonces esta simplificación deja de
ser válida, pues la distribución de Tt2 deja de ser χ2 (a).
¡Pero no hay que entrar en pánico! Pues, si se supone que la distribución conjunta de xt es Gaus-
siana, se puede realizar un test de hipótesis, con su respectivo umbral, basado en el hecho de que
Tt2 ; Fn,n−a , la llamada “distribución de Fisher”. Es ası́ como se da origen al test de Hotelling.
2.4.1. Test de Hotelling
Sea η un umbral. Si xt P Λ−1 T T

a P xt ≤ η, entonces se acepta H0 , para el umbral dado por:

a(n2 − 1)
Fa,n−a;α , si ~xt ∈
/ Training Set;



 n(n − a)


η= (2.7)
(n − 1)2 [a/(n − a − 1)]Fa,n−a−1;α


, si ~xt ∈ Training Set.


a

 n 1 + n−a−1 Fa,n−a−1;α
Es llamativo que existen dos casos para el cálculo del umbral. El primero representa una situación en
que se quiere comparar un dato con registros históricos y revisar si está dentro de la elipse de control.
El segundo, en cambio, sirve para purgar los datos con los que se realizará el PCA.
Este test suele ser efectivo (y de los más usados) para detectar anomalı́as (NO FALLAS) en la opera-
ción y además es común usar dos umbrales para dos valores de α distintos, definiendo ası́ dos niveles
de alerta en la elipse de control. Por ejemplo, para el ejemplo cities, considerando sólo las prime-
ras dos componentes como principales (a = 2) y haciendo el test de Hotelling para cada ciudad, se
obtienen las elipses siguientes elipses correspondientes al 95 % y roja 99 %.
4
-2
-4
-10 -5 0 5 10 15
1st Principal Component
Figura 2.12: Elipses para alerta amarilla (95 %) y roja (99 %).
Entre 95 % y 99 % Fuera del 99 %

(amarillo) (rojo)
Lafayette, LA Atlantic City, NJ

Newark, NJ Miami-Hialeah, FL
San Jose, CA Las Vegas, NV
Honolulu, HI Boston, MA
San Diego, CA Philadelphia, PA-NJ
Anaheim-Santa Ana, CA Washington, DC-MD-VA
Santa Barbara-Santa Maria-Lompoc, CA Los Angeles, Long Beach, CA
Baltimore, MD San Francisco, CA
West Palm Beach-Boca Raton-Delray Beach, FL Chicago, IL
Seattle, WA New York, NY
Pittsburgh, PA
Salinas-Seaside-Monterey, CA
Midland, TX
Tabla 2.2: Ciudades en los dos niveles de alerta
El análisis, sin embargo, no está completo, pues se debe tener en cuenta que la proyección reduce
la dimensión del espacio de variables a sólo las más significativas, por lo que existen pérdidas de
información asociadas a la proyección. Es por esto que se define el siguiente indicador:
E = (I − P P T )~xt , (2.8)
que da cuenta de las discrepancias entre ~xt y ~x. Este indicador es útil para detectar casos como el que
muestra la Figura 2.13, donde hay puntos que se proyectan dentro de la elipse, pero en realidad están
bastante lejos del punto de operación (y que pueden ir alejándose en el tiempo).
~x∗t
v2
v1
Figura 2.13: Caso útil para el indicador E.
Por ejemplo, haciendo el análisis de cities, considerando las tres primeras componentes como prin-
cipales, se obtienen puntos que no pasan el test de Hotelling, sin embargo se proyectan dentro de la
elipse:
EJERCICIOS 29

2
-2
-4
-10 -5 0 5 10 15
1st Principal Component
Figura 2.14: Caso útil para el indicador E en cities.
Este hecho se explica porque el cálculo de las distancias se hace considerando todas las componentes
principales, pero el gráfico bidimensional sólo entrega información respecto de las dos primeras,
por lo que ocurre que hay datos que tienen sus primeras dos componentes dentro de la elipse, pero la
tercera (saliendo del plano) está suficientemente lejos como para hacer que no pase el test de Hotelling.
Es por esto que para a > 2, las elipses son en realidad elipsoides en a dimensiones. Para el caso de la
Figura 2.14, lo más correcto es examinar un gráfico tridimensional.
EJERCICIOS
2.1 Realice el análisis PCA para cities, considerando las tres primeras componentes como principa-
les y muestre las elipsoides de control de cada alerta. ¿Qué factores tienen en común las ciudades que
están en ambos niveles de alerta?
2.2 Estudie la evolución del indicador E a medida que se aumenta el número de componentes que
se consideran principales para cities, es decir, cómo varı́a E para distintos valores de a. ¿Cómo se
interpreta el resultado cuando a = 9?
Problemas Resueltos
EJEMPLO 2.1
Considere el siguiente conjunto de datos bidimensionales:
C : (0, 0); (2, 1); (1, 1); (1, 2)

Encuentre los ejes de mayor varianza utilizando Analisis de Componentes Principales. Encuentre
los vectores bases y muestre gráficamente la ubicación de los ejes en el espacio original.
En efecto, para abordar este problema, se normalizarán los datos y luego se calculará la matriz de
covarianzas.
" #
0 2 1 1
S=
0 1 1 2
El promedio en cada eje está dado por:
" #
1
µ=
1
Luego la matriz de covarianzas:
Cov = (S − µ)(S − µ)T

 
" # −1 −1
−1 1 0 0  1 0 
Cov =
 
−1 0 0 1 0 0 
 

0 1
" #
2 1
Cov =
1 2
Calculamos los valores propios:
||Cov − λI|| = 0

2−λ 1
=0

1 2−λ

(2 − λ)2 − 1 = 0
(1 − λ)(3 − λ) = 0
Para λ = 1 se tiene:
" #" #
1 1 v1
= 0 ⇒ v1 = −v2
1 1 v2
De este modo el vector √12 (1, −1) es el vector propio para λ = 1 ⇒ √12 (1, 1) es el otro vector
propio asociado a λ = 3, dado que ambos vectores deben generar y ser ortonormales. Luego el
gráfico nos queda:
-1
-2
-3
-3 -2 -1 0 1 2 3
1. L. H. Chiang, E. L. Russell, R. D. Braatz, “Principal Component Analysis”, en Fault De-

tection and Diagnosis in Industrial Systems. Springer-Verlag London: London, 2001.
CAPÍTULO 3
ESTIMACIÓN
3.1. Introducción
El problema de la estimación de parámetros se centra en tratar de aproximar el valor de una determi-

nada incógnita. A diferencia del caso de detección, donde existe un número finito de alternativas, en
el problema de estimación existen infinitos posibles valores para el parámetro en cuestión.
En algunos casos, estaremos interesados en estimar el valor de un vector de parámetros que se supone
constante. En otros, se supone que existen variaciones del parámetro en el tiempo, pero que dichas
variaciones son lentas con respecto al periodo de muesteo de la señal.
Es importante destacar que existen dos corrientes de pensamiento que abordan estos problemas:
1. Frecuentista: La naturaleza fija el valor de un parámetro, pero hay incertidumbre asociada

a las mediciones realizadas.
2. Bayesiano: El parámetro es una variable aleatoria.
Esta discusión, puramente filosófica, está aún abierta. Dependiendo del caso de estudio, uno u otro
enfoque presentará ventajas y desventajas (como todo en la vida).
Una última observación es que existe una gran variedad de criterios para decidir que resultado del
proceso de estimación es “mejor”(mérito). Por ejemplo, la Figura 3.1 ilustra dos criterios de error
33
34 ESTIMACIÓN
Criterio 1 Criterio 2
3 3
2.5 2.5
2 2
v2
v2
1.5 1.5
1 1
0.5 0.5
0 0
0 1 2 3 0 1 2 3
v1 v1
(a) (b)
Figura 3.1: Dos criterios diferentes para estimar el comportamiento de un set de datos.
distintos que pueden usarse para definir la calidad de una aproximación lineal a la relación entre
dos variables de proceso. En un caso, el error se mide en la vertical, mientras que en el segundo se
considera el error de proyección en la recta. Ambos son válidos. ¿Cuál usar entonces? Pues bien,
depende de nosotros.
En este capı́tulo se abordarán cuatro enfoques: error cuadrático medio, mı́nimos cuadrados, máxi-
ma verosmilitud y máximo a posteriori. Una vez expuestos estos enfoques, se entregarán nociones
de suficiencia y eficiencia.
3.2. Estimadores MSE (Mean Square Error) y LS (Least Squares)
El enfoque MSE busca encontrar la proyección de la variable a estimar (y) sobre el espacio definido
por la información disponible (variables medidas, x); ver Figura 3.2. Dicha proyección, geométrica-
mente, asegura el mı́nimo error cuadrático medio al tratar de aproximar y.
~x2 ~y
~x1
P~x {~y }
Figura 3.2: Proyección de y sobre x para encontrar el mı́nimo error.
Para realizar una proyección es necesario primero definir un producto interno entre vectores. El pro-
blema aquı́ es que tanto x1 como x2 son variables aleatorias, por lo que dejan de ser simples vectores
en Rn para convertirse en funciones con una densidad de probabilidad asociada. La buena noticia es
que el producto interno entre v.a. no es complejo de definir. En efecto, se tiene que si X e Y son v.a.,
entonces una extrapolación del concepto de producto interno entre X e Y puede definirse en base
ESTIMADORES MSE Y LS 35
al caso tı́pico de producto punto en Rn , pero considerando un vector de infinitas componentes. Cada
componente representa una realización de las v.a., por lo que si los eventos son equiprobables se tiene:
hX, Y i = [x1 (ω1 ), x2 (ω2 ), . . . ] · [y1 (ω1 ), y2 (ω2 ), . . . ]T
Si existe una distribución de probabilidad conjunta para X e Y , entonces cada realización tiene un
ponderador
R asociado: la probabilidad del evento. Ası́, se puede pensar en una matriz de ponderadores
P = fXY , diagonal, con lo que se obtiene:
hX, Y i =~x(ω) · P · ~y (ω)T (3.1)

=x1 (ω1 )y1 (ω1 ) · P (ω1 ) + x2 (ω2 )y2 (ω2 ) · P (ω2 ) + . . .
=E[X · Y ] (3.2)
En el caso de una v.a. X con media nula, entonces se puede definir el “tamaño”de X en función de la
definición de distancia para un Espacio de Hilbert. Ası́, la varianza de la v.a. X resulta ser una métrica
y permite medir el tamaño de dicha función:
kxk2 = hx, xi = E[~x · ~xT ] = V[x] (3.3)
El estimador MMSE (Minimum Mean Squared Error) busca minimizar el tamaño de la v.a. error de
estimación usando esta métrica.
3.2.1. Esperanza Condicional
Definición 3.1. (Esperanza Condicional) La esperanza condicional, h, de una v.a. y , dado un pro-
ceso estocástico {X(τ ), a ≤ τ ≤ b} es la v.a. que:
1. Es un funcional de X(·), digamos h({X(τ ), a ≤ τ ≤ b}).

2. Satisface la condición de ortogonalidad dada por 1 :
E {[y − h(X(τ ))] · g(X(τ ))} = 0
Es bueno recalcar que, como lo dicta la definición, la esperanza condicional es una variable aleatoria,
pues es una función determinı́stica de un proceso estocástico.
Gráficamente, la situación del lado izquierdo de la condición de ortogonalidad se ilustra en la Figu-
ra 3.3, donde el vector y − h que cumple la definición es único, y coincide con la proyección de y
sobre el plano definido por la información disponible (x1 ):
1 Resulta conveniente ver la definición de esperanza como un producto interno.

36 ESTIMACIÓN
~y
~y − ~h
~h
P~x1 {~y}
Figura 3.3: Representación geométrica de la condición de ortogonalidad.
Por lo anteriormente descrito, la esperanza condicional es el estimador que minimiza el error cuadráti-
co medio, por lo que se le suele llamar MMSEE (por su sigla en inglés Minimum Mean Squared Error
Estimator). Si las distribuciones de probabilidad asociadas están bien definidas, este estimador se
calcula en base a la definición clásica de esperanza:
Z +∞
h({X(τ ), a ≤ τ ≤ b}) = Ey|x = y · fy|x (y|x)dy (3.4)
−∞
Ahora bien, esta expresión no siempre genera una solución analı́tica al problema de estimación. Más
aún, la densidad fy|x no es fácil de obtener en la práctica (con suerte, simplemente tenemos una
planilla de datos). Esos problemas motivan la búsqueda de estimadores sub-óptimos (en el sentido
MSE), pero que sean simples de implementar.
Una observación pertinente tiene relación con el hecho de que la esperanza se denomina condicional
porque depende de la realización meidad para x,. Tal como se observa en la Figura 3.4, y mostrará una
distribución condicional distinta para cada valor de x, incluso en casos donde la distribución conjunta
sea uniforme:
Distribuciones diferentes
y
fxy = cte
x1 x2 x3 x4
x
Figura 3.4: Región con distribución uniforme que muestra la condicionalidad de la esperanza.
Tal como se mencionó anteriormente, si se quiere forzar un estimador lineal, esto es, ŷ = Ax + b, en
vez de la estructura propuesta por el estimador esperanza condicional, entonces se sale del contexto
ESTIMADORES MSE Y LS 37
de minimizar el error cuadrático medio. En este caso se habla del estimador lineal de mı́nimo error
cuadrático medio. Aquı́, el concepto de mejor tiene, nuevamente, que ver con la minimización del
error cuadrático medio.
3.2.2. Estimador Lineal de Mı́nimo Error Cuadrático Medio
En este caso se propone encontrar un estimador lineal (estructuras afines pueden incluirse si se asume
una entrada constante) que cumpla con ser insesgado y con el criterio de ortogonalidad. Considerando
una estructura del tipo ŷ = Ax + b, se deben cumplir entonces, las siguientes condiciones:
1. El estimador debe ser insesgado, es decir, E[e] = 0. Esto tiene como consecuencia que, en
promedio, el estimador acierta y entrega el valor correcto del parámetro.
2. El error de estimación debe ser ortogonal al espacio de variables explicativas, es decir,
cov(e, x) = 0. Esto también puede interpretarse de la siguiente forma: no hay información
que no esté contemplada (o descrita) por el vector x que sea útil para la estimación.
De la primera condición, se debe cumplir que:
E[y − ŷ] =E[y − Ax + b] = 0
=⇒ b =E[y] − AE[x] (3.5)
Mientras que de la segunda condición, se tiene que
cov(e, x) =cov(y − Ax + b, x)
=cov(y, x) − cov(Ax, x) = 0
=⇒ A =cov(x, y) · cov(x, x)−1 (3.6)
Con esto, se llega a la conocida fórmula para el estimador LMMSE (Linear Minimum Mean Squared
Error):
ŷ = E[y] + cov(y, x)cov(x, x)−1 (x − E[x]) (3.7)

Se debe tener presente que éste es estimador es menos preciso que el MMSEE, puesto que cov(e, e)
será siempre superior que en el caso del estimador MSE. Este es el precio que se debe pagar por forzar
una estructura más simple en el estimador.
Por otra parte, cuando cov(x, x) no es invertible o existen variables altamente correlacionadas, ocurre
que para un estimador de la forma y = α1 x1 + α2 x2 , éste se puede aproximar por un sólo término
38 ESTIMACIÓN
y ≈ (α1 + α2 )x1 , lo que tiene como consecuencia que se deben encontrar dos incógnitas (α1 y α2 ),
pero se tiene, prácticamente, una sola ecuación. Este problema suele presentarse cuando las varianzas
de α1 y α2 son demasiado grandes, por lo que no se tiene certeza de los valores que se deben usar para
estimar.
3.2.3. Least Squares (LS)
Tanto en el caso MSE como en el LMSE se asume que los momentos de segundo orden asociadas
a las variables aleatorias son conocidas. En la práctica, esto no es el caso y tı́picamente dichos mo-
mentos deben estimarse de los datos (medias y matrices de varianza-covarianza empı́ricas). En esta
situación, no puede asegurarse que los estimadores desarrollados minimicen el error cuadrático me-
dio. Sólo puede asegurarse que se minimice el error cuadrático asociado al ajuste para los datos que
se adquirieron. Es ası́ como aparece el concepto de mı́nimos cuadrados (LS, Least Squares en inglés).
En caso de que y = α1 x1 + α2 x2 + ... + αn xn + tenemos que Ê[y~x] = Ê[xxT ]~

α
3.2.4. Predictores MSE para Sistemas Lineales
Supongamos que existe un modelo de la forma:
y = α1 x1 + α2 x2 + ... + αn xn + ω (3.8)
Donde, αi son conocidos y ω es ruido blanco.

Si multiplicamos 3.8 por x1 , x2 , ..., xn , tenemos un set de ecuaciones:
 2
 yx1 = α1 x1 + · · · + αn x1 xn + ωx1

..
 .
yxn = α1 x1 xn + · · · + αn x2n + ωxn

Si aplicamos el operador esperanza:
 2   
 x1 · · · xi xj
 
 α1
E[xy] = E .. ... .. .  ...  + E[x]E[ω]
 
 . . 
xj xi · · · x2n αn
 
Como E[ω] = 0, finalmente:

ESTIMADOR DE MÁXIMA VEROSIMILITUD 39
 2   
 x1 · · · xi xj
 
 α1
E[xy] = E .. .. .. .  ... 
 
(3.9)
. . .
 
xj xi · · · x2n αn
 
Una observación importante es que se debe calcular la matriz de esperanzas de la Ecuación 3.2.4,
resulta muy común que esta matriz no se puede calcular de forma explı́cita, por que se suele estimar.
En este caso se habla de un estimador LS (Least Squares).
3.3. Estimador de Máxima Verosimilitud
Conceptualmente, este estimador pretende encontrar un θ̂ tal que, realizado el experimento, sea el
más verosı́mil, es decir, que el estimador para el parámetro que convierta al experimento realizado en
el caso que se observarı́a con mayor frecuencia. Este razonamiento lleva al concepto del estimador
MLE (por sus siglas en inglés Maximun Likelihood Estimator).
Definición 3.2. Se define la función de verosimilitud, en el caso discreto, a:
`(θ, x) = `(θ, x1 , ..., xm ) = P[x1 , ..., xm |θ] (3.10)

Definición 3.3. Se define la función de verosimilitud en el caso continuo a:
`(θ, x) = `(θ, x1 , ..., xm ) = fx1 ,...,xm |θ [x1 , ..., xm |θ] (3.11)
En caso de que las variables x1 , ..., xm sean i.i.d. se tiene que la función de verosimulitud se puede
escribir como:
m
Y
`(θ, x) = fxi |θ (xi |θ) (3.12)
i=1
Como generalmente se trabaja con distribuciones de la familia de exponenciales, es usual considerar

la función ∆(θ, x) = log `(θ, x). Dada la propiedades del logaritmo (función continua, monótona y
estrictamente creciente) el estimador que maximiza ∆(θ, x), también maximiza `(θ, x).
El estimador MLE tiene las siguientes propiedades:
Es consistente.
Es asintóticamente insesgado.
Distribuye asintóticamente como una Gaussiana y además es asintóticamente eficiente: la
varianza del estimador alcanza la Cota de Cràmer-Rao.
Es invariante frente a funciones no lineales. fˆ(θ)M L = f (θ̂M L )
40 ESTIMACIÓN
EJEMPLO 3.1
Sean
y(t) = θ + n(t)
n
Y 1 −(xi −θ)2
f (y|θ) = √ e 2σ2
i=1
2πσ
Se quiere asignar θ tal que la secuencia medida {y(1), y(2), ..., y(n)} represente el resultado más
frecuente.
EJEMPLO 3.2
Considere un set de m muestras, las cuales pueden ser modeladas como xi = θ + ni , donde el el
ruido distribuye ni ∼ N (0, σ 2 ) y se quiere encontrar una expresión para el estimador de máxima
verosimilitud para el parámetro θ y para la varianza σ 2 .
Respuesta: Dado que las muestras son i.i.d., la distribución conjunta se calcula de la forma
n
Y1 (x(i)−θ)2
`(x(1), x(2), ..., x(m), θ) = e− 2σ2
√
i=1
2πσ
m P
1 m (x(i)−θ)
2
= √ e− 1 2σ2
2πσ
√ m
X (x(i) − θ)2
∆(x, θ) = log(`(x, θ)) = −m · log ( 2πσ −
1
2σ 2
Para obtener el valor del estimador se debe derivar ∆(x, θ) con respecto al (o los) parámetro(s),
ası́ entonces:
m
X 2(x(i) − θ̂M L )
∂∆(x, θ)
= 0 =⇒ − =0
∂θ 1
2σ 2
m
X m
X
⇐⇒ x(i) = θ̂M L
1 1
m
X
⇐⇒ x(i) = (m − 1 + 1)θ̂M L
1
m
X xi
⇐⇒ = θ̂M L
i=1
m
ESTIMADOR MÁXIMO A POSTERIORI 41
Procediendo de forma análoga para la estimar la varianza se maximiza ∆(x, θ) con respecto a
v = σ 2 , de esta manera se tiene:
m
√ X (x(i) − θ)2
∆(x, θ) = −m · log( 2πv) −
1
2v
X (x(i) − θ)2 m
∂∆(x, θ) 1 1 1
= 0 =⇒ −m √ √ 2π − − 2
=0
∂v 2πv 2 2πv 1
2v
m 2
m X (x(i) − θ)
− + =0
2v 1
2v 2
m m
X (x(i) − θ)2 X (x(i) − θ)2
= m =⇒ = v∗
1
v 1
m
Donde v ∗ = σ̂M
2
L
3.4. Estimador Máximo a Posteriori
A diferencia del estimador de máxima verosimilitud, donde el valor θ es un parámetro constante y

de valor único que se quiere determinar, el estimador máximo a posteriori considera que θ es una
variable aleatoria (por lo tanto en cada experimento puede tener un valor distinto).
Este estimador pretende maximizar directamente la distribución a posteriori del parámetro: fθ|x (θ|x) =
fx|θ (x|θ)fθ (θ)/f (x). El problema que surge es que f (x) no es trivial de calcular, ya que implica carac-
terizar completamente la probabilidad de recibir la información adquirida. Este problema se resuelve
rápidamente si se aplica logaritmo:
log(fθ|x (θ|x)) = ∆(θ, x) + log(fθ (θ)) − log(f (x)) (3.13)
Notemos que log(fθ|x (θ|x)) caracteriza el conocimiento a posteriori y log(fθ (θ)) el conocimiento a
priori.
Al maximizar la distribución a posteriori, log(f (x)) desaparece (al derivar con respecto a θ), por lo
que este término en la práctica no se considera en el análisis.
EJEMPLO 3.3
m
Y 1 (xi − θ)2
Sea la distribución conjunta condicional fx1 ,...,xm |θ = √ e−
i=1
2πσ 2σ 2
42 ESTIMACIÓN
2
− θ2
Si se asume que θ ∼ N (0, σθ2 ), tenemos que fθ (θ) = √ 1 e σθ , con lo que el estimador máximo
2πσθ
a posteriori es:
m
σ2 1 X
θ̂M AP = 2 θ2 xi
σθ + σ /m m i=1
Notemos que σθ es conocimiento a priori. Si hay gran incertidumbre respecto del valor del paráme-
tro, los datos son decidores a la hora de calcular el estimador. Si hay poca incertidumbre, el conoci-
miento a priori decide. Sin embargo, el estimador sólo puede mejorar si se tiene más información.
3.5. Estadı́sticos Suficientes
Un estádistico es una función de v.a.. Los estadı́sticos pretenden resumir la información disponible en
la distribución de probabilidad. El estadı́stico se declara suficiente si permite escribir: fx|T (x|T, θ) =
fx|T (x|T ). De lo anterior se desprende que “no se necesita conocer θ”si conozco el estadı́stico (pues
θ puede calcularse a partir del estadı́stico).
Definición 3.4. Sea t una función de v.a., se dice que t es estadı́stico suficiente si se puede escribir:
fx|θ (x|θ) = b(t(x), θ) · a(x)
.
Resulta interesante notar que la función a(x) puede ser la unidad; es decir, si se puede escribir
fx|θ (x|θ) = b(t(x), θ), entonces t es estadı́stico suficiente.
EJEMPLO 3.4
Consideremos la función de densidad de probabilidad fx y los estadı́sticos t1 y t2 dados por:
Pn
i=1 (xi − µ)2
1 −
fx (x|µ, σ) = √ n e 2σ 2
2πσ
 n
1X
t (x) = xi

 1

n i=1

n
1X
(xi − t1 )2

 t2 (x) =


n i=1
ESTADÍSTICOS SUFICIENTES 43
Notemos que fx (x|µ, σ) puede escribirse como:
n(t1 − µ)2 )

1 −nt2
fx (x|µ, σ) = √ n exp exp −
2πσ 2σ 2 2σ 2
Este resultado no depende de xi . Luego, t1 y t2 son estadı́sticos suficientes.
EJEMPLO 3.5
Sean x1 , x2 , ..., xn i.i.d. Bernoulli β(p) donde p es desconocido. La densidad de probabilidad

conjunta está dada por:
n
Y
fx|p (x|p) = pxi (1 − p)1−xi , xi = 0, 1
i=1
Se pide demostrar que un estadı́stico suficiente es t = ni=1 xi . Demostremos que es un estadı́stico

P
suficiente y encontremos en estimador de p donde p̂M L = p̂M L (t)
Trabajamos fx|p (x|p):
Pn Pn
fx|p (x|p) = pPi=1 xi (1 − p) i=1 1−xi
n Pn
= p i=1 xi (1 − p)n− i=1 xi
= pt (1 − p)n−t
Para encontrar el estimador de máxima verosimilitud derivamos la función de log-verosimilitud y

buscamos el estimador p̂M L que la maximice:
∆(p, t) = t log(p) + (n − t) log(1 − p)

∂∆ t (n − t)
(p̂M L , t) = − =0
∂p p̂M L 1 − p̂M L
t
=⇒ p̂M L =
n
Es interesante notar que:
44 ESTIMACIÓN
P
X P (x1 = x1 , .., xn = xn , T = xi = t/p)
f (x1 , x2 , ..., xn | xi , p) = P
P (T = xi = t|p)
pt (1 − p)n−t
= n t
t
p (1 − p)n−t
−1
n
=
t
Finalmente, no depende de p.
Ası́, se puede demostrar de las dos formas que es un estadı́stico suficiente.
3.6. Eficiencia de un Estimador
3.6.1. Score Function
La idea es entregar un puntaje a la forma de la función verosimilitud.
Definición 3.5. La función de puntaje (score function) se define como sigue:
∂log(`(x, θ))
s(x, θ) = (3.14)
∂θ
Con la definción anterior, resulta interesante recordar que s(x, θ̂M L ) = 0. Es decir, la función toma
valores cercanos a cero para valores cercanos al MLE.
Una propiedad importante de esta función es la que resulta al multiplicarla con un estadı́stico t.
Teorema 3.1. Sea s la función de puntaje y t un estadı́stico, se tiene entonces que:
∂E[tT ] ∂tT
E[s · tT ] = − E[ ] (3.15)
∂θ ∂θ
EFICIENCIA DE UN ESTIMADOR 45
Dem: Dado que t es función de x, su esperanza se calcula como:

Z
E[t ] = `(x, θ)tT dx
T
R
∂E[tT ] ∂ `(x, θ) · tT dx
=
∂θ ∂θ
∂tT
Z Z
∂`(x, θ) T
= `(x, θ) · dx + t dx
∂θ ∂θ
∂tT
= E[ ] + E[s · tT ]
∂θ

3.6.2. Cota de Cramér-Rao
La Cota de Cràmer-Rao representa la cota para la varianza de un estimador insesgado.
Nos hacemos la siguiente pregunta: si se quiere encontrar un estimador θ̂ para θ tal que E[θ̂] = θ
(insesgado) ¿Cuál es la mı́nima varianza V[θ̂] ?
Para esto, definimos lo siguiente:
Definición 3.6. La Matriz de Información de Fisher corresponde a la matriz de covarianza de la
función de puntaje, y se denota por J(θ). Dado que la función de puntaje es de media nula, se tiene
que:
" #
∂(log(`(θ, x))) ∂(log(`(θ, x))) T
J(θ) = E[s(θ, X)sT (θ, X)] = E (3.16)
∂θ ∂θ
Luego, la información de Fisher se puede escribir como:
∂ ∂(log(`(θ, x))) T
J(θ) = −E[ ] (3.17)
∂θ ∂θ
Teorema 3.2. El Teorema de Crámer-Rao indica que, si θ̂ es un estimador insesgado de θ basado
en una función de verosimilitud diferenciable, entonces:
E[θ̂ − θ][θ̂ − θ]T ≥ J −1 (θ) (3.18)
Considerando lo anterior, podemos definir un estimador eficiente.

Definición 3.7. Sea θ̂ un estimador de θ. Se dice que θ̂ es eficiente si, y solo si, es insesgado y su
error de estimación alcanza la cota de Crámer-Rao, es decir, si cumple la igualdad:
46 ESTIMACIÓN
E[θ̂ − θ][θ̂ − θ]T = J −1 (θ) (3.19)
Teorema 3.3. Un estimador insesgado θ̂ es eficiente sı́ y sólo si:
J(θ)(θ̂ − θ) = s(θ, X) (3.20)
Más aún, cualquier estimador insesgado eficiente es un estimador de máxima verosimiltud.

Problemas Resueltos
EJEMPLO 3.6
Sean X e Y v.a. que distribuyen uniformemente sobre la región triangular del plano x-y definida
por los vértices (0,0), (0,0.8), y (1,2). Encuentre el estimador LMMSE de X dado Y, y el estimador
MMSE de X dado Y. Calcule el MSE (mean squared error) para cada estimador. ¿Qué porcentaje
de reducción ofrece el estimador MMSE respecto del LMMSE?
Area de las variables aleatorias

2
1.8
1.6
1.4
1.2
Y
0.8
0.6
0.4
0.2
0
0 0.2 0.4 0.6 0.8 1
X
Figura 3.5: Región definida por las variables aleatorias
La Figura 3.5 muestra la región comprendida entro los puntos (0, 0) (0, 45 ) (1, 2), cuya área es
base · altura 0,8 · 1
= . Dado que X,Y distribuyen uniforme sobre esta región, la probabilidad
2 2
conjunta está dada por:
1
fX,Y (x, y) =
0,4
48 ESTIMACIÓN
MMSE Finalmente:
Z y −1
2 1 5
Con lo anterior, se puede calcular el MM- E[X|Y ] = y− 4 x − (y − 2) dx
SE por definición como E[X|Y ]. Nótese que 6
5 0,4 6
5
este MMSE está definido por partes, cuan- 1
do 0 ≤ Y ≤ 0,8 y 0,8 ≤ Y ≤ 2. Caso = (2y − 1)
0 ≤ Y ≤ 0,8 3
Z y Y su MSE es E[(x − x̂)2 ] = 0,00667.
2
E[X|Y ] = xfX|Y (x|y)dx
0 LMMSE
Luego se debe calcular la distribución condi-
cional, lo cual se puede calcular mediante la Para calcular el LMMSE se utilizará la formu-
siguiente identidad: la directamente:
σXY
fX,Y (x, y) = fX|Y (x|y)fy (y) x̂ = (y − ȳ) + x̄
σY2
El lado derecho de la identidad se conoce, se Para esto se debe calcular las covarianzas y las
debe calcular entonces fy (y) el cual se puede respectivas esperanzas.
calcular por definición:
Z y Cov(X, Y ) = E[XY ] − E[X] ∗ E[Y ]
2
fy (y) = fX,Y (x, y)dx V[Y ] = E[Y 2 ] − E2 [Y ]
0
y
1 5y Se procede a calcular cada termino:
Z
2
= dx =
0 0,4 4 1 6
x+ 45
1 2
Z Z
5
De esta manera, se calcula el MMSE: E[XY ] = xy dydx =

0 2x 0,4 5
6
1 x+ 45
1 1
Z y −1 Z Z
5
2 1 5y y E[X] = x dydx =
E[X|Y ] = x dx = 0,4 3
0 0,4 4 4 0 2x
6
1 x+ 54
1 14
Z Z
5
E[Y ] = y dydx =
Caso 0,8 ≤ Y ≤ 2 Para este caso, se pro- 0 2x 0,4 15
6
cede de la misma manera: 1 x+ 45
1 26
Z Z
5
E[Y 2 ] = y2 dydx =
Z y
2 0 2x 0,4 25
E[X|Y ] = y− 4 xfX|Y (x|y)dx
6
5 De este modo, se tiene:
5
4 38
Donde: Cov(X, Y ) = ; V[Y ] =
45 225
Z y
2
fy (y) = fX,Y (x, y)dx Finalmente el LMMSE está dado por:
0
Z y
1 5 10 14 1
=
2
dx = − (y − 2) x̂ = (y − ) +
y− 4
5 0,4 6 19 15 3
6
5
Y el MSE E[(x − x̂)2 ] = 0,453216
Si ahora la distribución de Y cambia a una normal N (0, 1) y que X = |Y |, los estimadores

cambian un poco:
MMSE LMMSE
Dada las condiciones del problema, nuestras Por definición E[Y ] = 0. Ahora calculemos
distribuciones queda como sigue: E[XY ]
y2 |y|2 y2
e− 2 ∞ ∞
e− 2 e− 2
Z Z
fY (y) = √ E[XY ] = √ √ dxdy
2π −∞ 0 2π 2π
|y|2
e − 2 =0
fX|Y (x|y) = √ ⇒Cov(X, Y ) = 0
2π
( q 2
2 − y2 Finalmente el estimador buscado es:
e y>0
⇒ fX (y) = π
0 si no
r
2
x̂ = E[X] =
π
4
De esta manera, MMSE se calcula como si- Donde su MSE es E[(x − x̂)2 ] = +1
π
gue:
|y|2
∞
e−
Z
2
E[X|Y ] = x√ dy = x = |y|
0 2π
Luego el MSE: E[(x − x̂)2 ] = 0
Análisis
Resulta interesante notar que para ambos casos, que el MSE del estimador MMSE es menor que
el MSE del estimador LMMSE, esto dado que el LMMSE tiene como hipóstesis algo muy fuerte,
lo cual no minimiza necesariamente el Error Cuadrático Medio.
Solo en a) tiene sentido calcular un porcentaje de reducción, dado que en b) el MSE del MMSE
es 0.
0,453216 − 0,0066666
%= · 100 = 98,52 %
0,453216
50 ESTIMACIÓN
EJEMPLO 3.7
Considere el problema donde se mide x[n] = Arn + w[n] para n = 0, 1, . . . N − 1, donde w[n]
es ruido Gaussiano blanco con varianza σ 2 , r > 0 y A quiere ser estimado. Encuentre la cota
de Cramér-Rao para A. Muestre que un estimador eficiente existe y encuentre su varianza. ¿Qué
sucede si N → ∞, considerando distintos valores de r?
Se tiene que la cota de Crammer-Rao está dada por:
1 1
C −R= = ∂ 2
I(θ) −E[ ∂A2 log(`(x))|θ]
Luego:
N −1
1 1 X ∂
log(`(x)) = (N − 1) log √ − 2 (xi − Ari )2 /
2πσ 2σ i=0 ∂A
N −1
1 X ∂
S(~x|A) = 2 2(xi − Ari ) · ri /
2σ i=0 ∂A
N −1
1 X 2i 1
=− 2 r /−
σ i=0 E[·|A]
1 1 − r2
= σ2
I(A) 1 − r2N
El estimador eficiente que buscamos es el mismo estimador de log-verosimilitud. En Efecto:

N −1 N −1
1 X i i 1 − r2 X
(xi − Âr ) · r = 0 ⇒ Â = xi r i
σ 2 i=0 1 − r2N i=0
Notemos que este estimador es insesgado, dado que:

N −1
1 − r2 X
E[Â] = E[ xi r i ]
1 − r2N i=0
N −1
1 − r2 X
= E[xi ]ri
1 − r2N i=0
N −1
1 − r2 X
= 2N
E[Ari + wi ]ri
1−r i=0
N −1
1 − r2 X 2i
= A r
1 − r2N i=0
=A
Luego solo basta observar que:

N −1
1 X
S(~x|A) = 2 (xi ri − Ar2i )
σ i=0
N −1
!
1 1 − r2 1 − r2N X
= 2 (xi ri ) − A
σ 1 − r2N 1 − r2 i=0
1 1 − r2
= Â − A
σ 2 1 − r2N

Es decir, es eficiente y su varianza está dada por la Cota de Cramér-Rao:
1 − r2 2
V[Â] = σ
1 − r2N
Resulta interesante notar que:
r < 1 ⇒ V[Â] → (1 − r2 )σ 2
σ2
r = 1 Se debe recalcular la fórmula. Luego V[Â] = I(A) = N
.
r > 1 ⇒ V[Â] → 0
52 ESTIMACIÓN
EJEMPLO 3.8
Sea X = [x1, x2, x3]T un vector Gaussiano con esperanza y matriz de covarianza dados por
[0, 0, 0]T y [1 0,7 0,5; 0,7 4 0,2; 0,5 0,2 3], respectivamente. Determine los coeficientes óptimos
para el predictor x̂3 = c1 x1 + c2 x2 y el mı́nimo error medio cuadrático asociado. ¿Cómo mo-
dificarı́a la estructura del predictor y cómo cambiarı́an sus resultados si la esperanza de X fuera
[123]T ?
~ sobre (x1 , x2 ). De este modo, se desea calcular:
Se calculará la proyección K
~ = P~
RK
Donde:
" # " #
E[X1 , X1 ] E[X1 , X2 ] E[X3 , X1 ]
R= P =
E[X2 , X1 ] E[X2 , X2 ] E[X3 , X2 ]
Para calcular E[Xi , Xj] se recuerda que:
Cov(Xi , Xj ) = E[Xi , Xj ] − E[Xi ]E[Xj ]

E[Xi , Xj ] = Cov(Xi , Xj ) + E[Xi ]E[Xj ]
Como E[Xi ] = 0∀i, se desprende que los vectores y matrices anteriores son:
" # " #
1 0,7 0,5
R= P =
0,7 4 0,2
Ası́: " #
62
~ =
K 117
5
− 117
Luego el MSE:
107
E[(x3 − xˆ3 )2 ] = S33 + µ23 − P~ T R−1 P~ =
39
Ahora bien, si las esperanzas cambian y ya no son cero, el predictor debe considerar en su estruc-
tura estas, de esta forma:
" # " # " #

2 2,7 3,5 ~ = 1,2928
R= P = ⇒K , E[(x3 − xˆ3 )2 ] = 5,3754
2,7 8 6,2 0,3387
EJEMPLO 3.9
Considere X1 , X2 , ...Xn variables aleatorias i.i.d. con f.d.p.:
f (x|θ) = θxθ−1 Para 0 ≤ x ≤ 1, 0 ≤ θ ≤ ∞
Encuentre el estimador de máxima verosimilitud para θ

Se calculará el estimador de log-verosimilitud.
N
Y
f~x (X|θ) = θxθ−1
i / log(·)
i=1
N
X d(·)
= N log(θ) + (θ − 1) log(xi ) /
i=1
dθ
N
1 X
0=N + log(xi )
θ̂ i=1
1
θ̂ = − 1
PN
N i=1 log(xi )
1. T. K. Moon, W.C. Stirling, “Estimation Theory”, in Mathematical Methods and Algorithms

2. J. V. Candy, “Bayesion Estimation”, in Bayesian Signal Processing Classical, Modern, and
Particle Filtering Methods. John Wiley & Sons, Inc., Hoboken: New Jersey, 2009.
CAPÍTULO 4
FILTRO DE KALMAN
4.1. Introducción
Este capı́tulo presenta una de las posibles formas en las que pueden deducirse las ecuaciones del Filtro
de Kalman. Iniciamos esta exposición recordando algunos conceptos básicos, relativos a sistemas
lineales y el observador de Luenberger.
Durante la primera parte del curso, se estudió que para sistemas dinámicos lineales de tiempo discreto
se puede usar la siguiente representación en variables de estado:
x(n + 1) =A(n) · x(n) + B(n) · u(n)

y(n) =C(n) · x(n) + D(n) · u(n)
Para el caso de un sistema invariante en el tiempo (i.e., A, B, C y D constantes), se demostró que el

observador de Luenberger tiene la siguiente estructura:
x̂(k) = A · x̂(k − 1) + B · u(k − 1) + F (y − ŷ) (4.1)
En este contexto, si se conocen completamente las matrices A, B, C y D, las únicas fuentes de error
para el observador se asociaban al desconocimiento de las condiciones iniciales del proceso. Sin
embargo, en la realidad ocurre que las mediciones siempre están contaminadas por ruido. Es más,
55
56 FILTRO DE KALMAN
pocas veces podemos darnos el lujo de enunciar que conocemos completamente la estructura del
modelo del proceso y las matrices ya mencionadas. De este modo, lo primero que se debe considerar
es que la incertidumbre asociada al proceso de estimación del vector de estado (objetivo último del
observador) se puede caracterizar en función de principalmente dos elementos: el ruido de sensor, v,
y la incertidumbre realativa a la transición de estados (ruido de proceso), ω. Por lo tanto, un modelo
dinámico para sistemas lineales con incertidumbre de medición puede representarse:
x(n + 1) =A(n) · x(n) + B(n) · u(n) + ω(n) (4.2)

y(n) =C(n) · x(n) + D(n) · u(n) + v(n) (4.3)
De lo anterior, es fácil concluir que, bajo estas condiciones, el filtro de Luenberguer no es suficiente
para asegurar convergencia de su salida (estimador del vector de estado).
Es por este motivo que Kalman introdujo el concepto de una ganancia de aprendizaje variante en el
tiempo, es decir F = F (k). Conceptualmente, se desea ir corrigiendo el valor la ganancia F con-
siderando la información ya recopilada a través de mediciones, hasta lograr convergencia. En
particular, se busca que el estimador de estados (denominado Filtro de Kalman) cumpla las siguientes
propiedades:
1. Minimice el error cuadrático medio de estimación. 1

2. Insesgado.
3. Lineal, especificamente de la forma x̂(n|n) = K 0 x̂(n|n − 1) + Ky(n).
4. Recursivo.
Definición 4.1. Se denota, x̂(n|i), como la mejor estimación de x(n) dada la información hasta el
instante i.
Se dice entonces que:
x(n|i) es suavizador si n < i.

x(n|i) es filtro si n = i.
x(n|i) es predictor si n > i.
El caso particular x̂(n|n−1) se denomina predicción a un paso. Con todo lo anteriomente mencionado
se definen entonces:
Definición 4.2. Error de predicción: e(n|i) = x(n) − x̂(n|i)
Definición 4.3. Varianza del error de predicción: P (n|i) = E[e(n|i)eT (n|i)]
Los tres estadios asociados a una iteración del filtro de Kalman, se pueden observar en la Figura 4.1.
1 Con esta propiedad, se suele decir que el filtro de Kalman es óptimo, pero es importante recordar que sólo es óptimo en el sentido MSE.
DEDUCCIÓN ANALÍTICA 57
cción a un paso
Predi
Modelo Corrección
Figura 4.1: Etapas del filtro de Kalman
Ası́, a través del modelo del sistema, y el conocimiento acumulado hasta n − 1, se predice el valor
del vector de estado para el instante n,. Luego, habiéndose recibido la medición y(n), se realiza una
corrección a esta predicción a un paso. Ası́ se obtiene la estimación de estados válida para el instante
n.
4.2. Deducción Analı́tica
El filtro de Kalman pretende ser una herramienta de estimación que funcione en condiciones lo más
generales posibles. Este método se puede utilizar cuando el sistema es lineal y está afecto a ruido
aditivo. En esta sección, deduciremos las expresiones que permiten implementar el Filtro de Kalman,
bajo ciertos supuestos que sólo tienen un fin pedagógico. Primero demostraremos que el estimador
buscado tiene una forma similar a la del estimador de Luenberguer. Para esto se debe considerar la
condición de estimador insesgado, es decir:
E[e(n|n)] = E[x(n) − [K 0 (n)x̂(n|n − 1) + K(n)y(n)] = 0 (4.4)
Observemos que la expresion x̂(n|n − 1) representa la mejor predicción hasta n − 1, y que:
x̂(n|n − 1) = x(n) − e(n|n − 1) (4.5)
Imponiendo esta condición y observando que2 E[e(n|n − 1)] = 0 obtenemos que:
E[(I − K 0 (n) − K(n)C(n))x(n) + K 0 (n)e(n|n − 1) − K(n)v(n) + K(n)D(n)u(n)] = 0
La expresión anterior se reduce notablemente dado que v(n) es ruido blanco y asumiendo, sin pérdida
de generalidad, que D(n) = 0.
2 Esta demostración queda para el lector.

58 FILTRO DE KALMAN
Recordando que:
E[e(n|n − 1)] = 0 ⇒ E[K 0 (n)e(n|n − 1)] = 0 (4.6)
E imponiendo que la condición de insesgo se cumpla para toda muestra x(n), se tiene:
∀x(n), (I − K 0 (n) − K(n)C(n)) = 0
Con esto finalmente se tiene que:
x̂(n|n) = (I − K(n)C(n))x̂(n|n − 1) + K(n)y(n) (4.7)

= x̂(n|n − 1) + K(n)[y(n) − C(n)x̂(n|n − 1)] (4.8)
Que es exactamente la forma que tiene el observador de Luenberguer.

Dado que ya conocemos la forma del predictor, consideremos el siguiente sistema donde se ha im-
puesto B = 0 (por motivos pedagógicos) y D = 0 (sin pérdida de generalidad). Además, asumiremos
que w(n) y v(n) son i.i.d. Es preciso explicitar que todos los supuestos aquı́ enunciados no son nece-
sariamente un requisito para que el Filtro de Kalman funcione.
x(n + 1) =A(n) · x(n) + ω(n)

y(n) =C(n) · x(n) + v(n)
Para este sistema, se puede calcular entonces:
x̂(n|n − 1) = Ax̂(n − 1|n − 1) + 0

P (n|n − 1) = E[e(n|n − 1)eT (n|n − 1)]
= E[Ae(n − 1|n − 1)eT (n − 1|n − 1)AT + Qw ]
= AP (n − 1|n − 1)AT + Qw
La expresión para la matriz de Kalman K(n) es el resultado de minimizar la varianza P .
P (n|n) = E[e(n|n)eT (n|n)]

= (I − K(n)C(n))P (n|n − 1)(I − K(n)C(n))T + K(n)Rvv K T (n)
Con este resultado y utilizando la norma de una matriz como la traza de la misma, se quiere:
Min T r[P (n|n)]

K(n)
Lo cual requiere que:
−2(I − K(n)C(n))C T (n)P (n|n − 1) + 2K(n)Rvv = 0
⇒ K(n) = P (n|n − 1)C T (n)[C(n)P (n|n − 1)C T (n) + Rvv ]−1

(4.9)
Problemas Resueltos
EJEMPLO 4.1
Se tiene el siguiente sistema ruido de sensor v(n), que se modela como ruido blanco Gaussiano:
x(n + 1) = x(n)
y(n) = x(n) + v(n)
El objetivo es iterar el algorimo del filtro de kalman, para observar el comportamiento del mismo.
Notemos que A = 1, B = 0, C = 1 y D = 0.
i)
x̂(n|n − 1) = x̂(n − 1) = x̂(n − 1|n − 1) Da igual escribir esta parte

p(n|n − 1) = p(n − 1) = p(n − 1|n − 1) Da igual escribir esta parte
ii)
p(n − 1)
K(n) =
p(n − 1) + σv2
iii)
p(n − 1)
x̂(n) = x̂(n − 1) + [y(n) − x̂(n − 1)]
p(n − 1) + σv2
σv2
p(n) = p(n − 1)
p(n − 1) + σv2
60 FILTRO DE KALMAN
para n = 1 se puede calcular, considerando x(0) y p(0) datos.
x̂(1|0) = x̂(0)
p̂(1|0) = p(0)
Luego:
p(0)
K(1) =
p(0) + σv2
p(0)
x̂(1) = x̂(0) + [y(1) − x̂(0)]
p(0) + σv2
Obs: La ganancia K es menor que p(0), es decir, es menor a la incertidumbre inicial.

Esto es debido al ruido de sensores.
Agrupando términos se tiene que:
σv2 p(0)
x̂(1) = 2
x̂(0) + y(1)
p(0) + σv p(0) + σv2
σv2
p(1) = p(0) < p(0)
p(0) + σv2
Notemos que estas ecuaciones son un promedio ponderado, es decir, se puede variar
cuanta importancia se le da a x̂(0) e y(1). Además, la incertidumbre disminuye.
En la siguiente iteración se obtendrá que, n = 2:
σv2
p(2) = p(0)
2p(0) + σv2
Se puede demostrar que en general:
σv2
p(n) = p(0)
n · p(0) + σv2
EJEMPLO 4.2
Este ejemplo es similar al anterior, pero considera ruido de proceso:
x(n + 1) = x(n) + ω(n)

y(n) = x(n) + v(n)
i)
x̂(n|n − 1) = x̂(n − 1)
p(n|n − 1) = p(n − 1) + σω2
Luego se tiene que:
p(n − 1) + σω2
K(n) =
p(n − 1) + σω2 + σv2
p(n − 1) + σω2
p(n) = 2 2
σv2
p(n − 1) + σω + σv
Obs: Notemos que en este caso, K(n) es mayor que en el caso sin ruido de proceso, es
decir no se le cree tanto al modelo, por lo que se da más peso al sensor.
PROPUESTO: Encontrar una expresión con n explı́cito para la ganancia de Kalman.
EJEMPLO 4.3
( P3.b) Examen Otoño 2015) Considere un sistema con dos sensores que toman medidas de una
constante desconocida x. Cada observación es ruidosa y puede ser modelada de la forma:
y(1) = x + v(1), y(2) = x + v(2)
donde v(1) y v(2) son variables aleatorias definidas por ruido no-correlacionado, media cero y
varianzas σ12 y σ22 . , respectivamente. En auscencia de información previa acerca del valor de x, se
busca el mejor estimador lineal de x de la forma:
x̂ = k1 y(1) + k2 y(2)
A partir del marco que proporciona el filtro de Kalman, encuentre los valores para k1 y k2 que
definen un estimador insesgado para x que minimiza el error cuadrático medio, E{(x − x̂)2 },
considere que los sensores toman sus medidas y(1) e y(2) alternada y secuencialmente.
Respuesta: Dadas las caracterı́sticas del problema, es posible trabajar el filtro de kalman con va-
lores escalares sin pérdida de generalidad, esto permite el reemplazo de matrices por coeficientes.
Ahora bien, para ajustar el problema al modelo propuesto en las expresiones 4,2 y 4,3 se debe
recordar que x es una constante y que mantendrá su valor en cada paso del sistema, con esto se
tiene que A = 1, B = 0, ω(n) = 0.
x(n + 1) = x(n)
62 FILTRO DE KALMAN
Por otro lado, se quiere una expresión para y(n) que represente la alternancia de los sensores , una
buena forma de hacerlo es mediante la expresión:
" #
v1
y(n) = x + D Donde D = [1 0] si n = 1, D = [0 1] si n = 2
v2
Para obtener el estimador que minimiza el ECM es necesario utilizar toda información disponible,
es decir, iterar en el filtro de kalman dos veces (Ver definición 4.1). Sean las ecuaciones del filtro:
1. x̂(n|n − 1) = Ax̂(n − 1) + Bu(n − 1)

2. P (n|n − 1) = AP (n − 1)AT + Qw
−1
3. K(n) = P (n|n − 1)C T C P (n|n − 1)C T + Qv

4. P (n|n) = I − K(n) C T P (n|n − 1))

5. x̂(n|n) = x̂(n|n − 1) + K(n) y(n) − C x̂(n|n − 1)
Los autores recomiendan fuertemente escribirlas en un papel, guardarlo en la billetera y revisar-

las de vez en cuando.
Recordando que A = 1, B = 0, Qw = 0 y C = 1 se opera el algoritmo con n=1.
1. x̂(1|0) = x̂(0)
2. P (1|0) = P (0)
3.
−1
K(1) = P (1|0) P (1|0) + σ12
P (0)
= → 1∗
P (0) + σ12
4.

P (1) = 1 − K(1) P (0)
σ12 P (0)
= → σ12
P (0) + σ12
Cuando no se tiene información previa sobre el parámetro, es decir x̂(0) es desconocido, es

recomendable tomar el lı́mite P (0) → ∞, de esta manera se asegura que x̂(0) puede tomar
cualquier valor y que no tendrá impacto en el estimador final, esto es debido a que el filtro
de kalman busca converger al parámetro sin importar sus condiciones iniciales. Luego por
simplicidad se utiliza x̂(0) = 0 y tomando el lı́mite mencionado se tiene que K(1) = 1 y
P (1) = σ12 .
5.

x̂(1|1) = x̂(1|0) + K(1) y(1) − x̂(1|0)
= y(1)
Siguiendo con el algoritmo para n = 2 se tiene:
1. x̂(2|1) = x̂(1) = y(1)
2. P (2|1) = P (1)
3.
P (2|1)
K(2) =
P (2|1) + σ22
P (1)
=
P (1) + σ22
σ2
= 2 1 2
σ1 + σ2
4. Notar en esta parte que el valor de P (2|2) no se utiliza en esta iteración, y dado que es la
última, se puede prescindir del cálculo.
5.

x̂(2|2) = x̂(2|1) + K(2) y(2) − x̂(2|1)
σ2
= y(1) + 2 1 2 y(2) − y(1)

σ1 + σ2
2
σ σ2
x̂ = 2 2 2 y(1) + 2 1 2 y(2)
σ1 + σ2 σ1 + σ2
Donde se ha encontrado el estimador de la forma propuesta en el enunciado. Es importante agregar

que la manera en que las varianzas de los sensores están ponderando a las salidas hace que el
estimador tenga una mayor credibilidad en aquel que tenga la menor varianza. (Esta conclusión
surge al analizar qué ocurre con el estimador si una de las varianzas tiende al infinito)
EJEMPLO 4.4
En este problema, usaremos el filtro de Kalman para identificación de sistemas. Asuma que te-
nemos un sistema lineal de tiempo discreto que se describe a través de la siguiente ecuación de
diferencias (u(n) es la entrada e y(n) es la salida):
64 FILTRO DE KALMAN
y(n) = −0,4y(n − 1) + 0,5y(n − 2) + u(n)

a) Con u(n) = 1 para n ≥ 0, genere un conjunto de salidas para y(n), para n = 0, 1, ..., 200
asumiendo condiciones iniciales nulas.
Se generarán usando las condiciones del enunciado, mediante un ciclo for dada la estructura re-
cursiva del sistema.
for i=3:n
y(i)=-0.4*y(i-1)+0.5*y(i-2)+u(i)+sqrt(sigma_v);
end
b) Con los datos generados en (a), use el filtro de Kalman para identificar los coeficientes de un
modelo de segundo orden para y(n) de la forma:
b0 + b1 z −1
H(z) =
1 + a1 z −1 + a2 z −2
Al correr el algoritmo un par de veces se nota que los resultados son relativamente similares. De
esta manera se puede observar la efectividad del filtro de kalmal para estimación de parámetros
de la función de transferencia.
Ejemplos >> p2s2 >> p2s2
ans = ans =
0.8644 0.8644
0.8644 0.8644
-0.4057 -0.4057
0.4942 0.4942
>> p2s2
ans =
0.8644
0.8644
-0.4057
0.4942
c) Grafique la salida del modelo encontrado en la parte (b) y, en el mismo gráfico, muestre la
salida del sistema real. Compare los resultados.
Sistema Real vs Sistema Estimado

2.5
1.5
y(n)
1
Sistema Real
0.5 Sistema Estimado
0
0 20 40 60 80 100 120 140 160 180 200
n
Figura 4.2: Sistema Real vs Estimado, entrada u(n) = 1
d) Repita las partes (a)-(c) cuando la entrada u(n) es una secuencia de muestras de ruido blanco
Gaussiano i.i.d. de varianza unitaria.
En esta parte se observa que los parámetros son mejor estimados que cuando la entrada es un
escalón. A continuación se muestran ejemplos y evolución del sistema.
Ejemplos >> p2s2 >> p2s2
ans = ans =
0.9931 1.0046
-0.0021 0.0016
-0.3989 -0.3973
0.4928 0.5039
>> p2s2
ans =
1.0110
0.0051
-0.3984
0.5213
Sistema Real vs Sistema Estimado
400
200
y(n) 0
-200
Sistema Real
-400 Sistema Estimado
-600
0 20 40 60 80 100 120 140 160 180 200
66 FILTRO DE KALMAN n
Figura 4.3: Sistema Real vs Estimado, entrada u(n) ∼ N (0, 100).
e) Grafique p(n) (varianza del error predicha) y s(n) (varianza real del error) versus n para n =
0, 1, ..., 150 y discuta lo que observa.
Caso u(n) constante En este caso, se observa que la varianza predicha aumenta casi de ma-
nera lineal, en cambio la varianza real se va a cero muy rápidamente.
×10 5 EMS de X1 (n) ×10 5 EMS de X2 (n)

8 8
P(n) P(n)
6 S(n) 6 S(n)
EMS
EMS
4 4
2 2
0 0
0 50 100 150 200 0 50 100 150 200
n n
×10 5 EMS de X3 (n) ×10 5 EMS de X4 (n)
6 6
P(n) P(n)
S(n) S(n)
4 4
EMS
EMS
2 2
0 0
0 50 100 150 200 0 50 100 150 200
n n
Figura 4.4: EMS para cada variable

×10 4 EMS de X1 (n) ×10 4 EMS de X2 (n)

3 3
P(n) P(n)
2 S(n) 2 S(n)
EMS
EMS
1 1
0 0
0 50 100 150 200 0 50 100 150 200
n n
×10 4 EMS de X3 (n) ×10 4 EMS de X4 (n)
3 2
P(n) P(n)
1.5
2 S(n) S(n)
EMS
EMS
1
1
0.5
0 0
0 50 100 150 200 0 50 100 150 200
n n
Figura 4.5: EMS para cada variable
EJEMPLO 4.5
Consideremos el circuito RC que muestra la Figura 4.6. Este representa la medición de la tensión
através del condensardor. Notemos que existe una fuente de incerteza asociada al instrumento de
medición.
R +
+
Ein C Eout V
-
Figura 4.6: Diagrama del circuito RC
a) Modele el sistema utilizando un modelo de ecuaciones diferenciales.

En efecto, utilizando ley de Kirchoff para las corrientes se tiene que:
68 FILTRO DE KALMAN
e(t) de(t)
Iin (t) − −C =0
R dt
Donde e(t = 0) = e0 la tensión inicial. La medida de la salida está dada por Ke ganancia asociada
al instrumento de medición.
eout (t) = Ke e(t)
b) Discretize el sistema encontrado en la parte anterior.

Para discretizar el sistema se utilizará la derivada hacia atrás. De este modo:

e(t) − e(t − 1) e(t − 1) ∆T ∆T
C =− + Iin (t − 1) ⇔ e(t) = 1 − e(t − 1) + Iin (t − 1)
∆T R RC C
c) Para el modelo discreto encontrado en la parte anterior, utilice C = 1000µF , R = 3,3kΩ,

∆T = 100ms, e0 = 2,5V , Ke = 2 e incerteza del sistema Rww = 0,0001. La presición del
voltimetro es de ±4V . Transforme el sistema a variables de estados con x(t) = e(t) y u(t) =
Iin (t). Finalmente, utilice un Filtro de Kalman para estimar la tensión en el condensador.
En efecto, el sistema modelado en variables de estados es el siguiente:
x(t) = 0,97x(t − 1) + 100u(t) + w(t − 1)

y(t) = 2x(t) + v(t)
Notemos que el filtro de kalman estará dado por:
x̂(t|t − 1) = 0,97x̂(t − 1|t − 1) + 100u(t − 1) [Estado Predicho]

P (t|t − 1) = 0,94P (t − 1|t − 1) + 0,0001 [Covarianza Predicha]
e(t) = y(t) − 2x̂(t|t − 1)
P (t|t − 1)
K(t) = 2 [Ganancia de Kalman]
4P (t|t − 1) + 4
x̂(t|t) = x̂(t|t − 1) + K(t)e(t) [Corrección de Estado]
P (t|t − 1)
P (t|t) = [Corrección de Covarianza]
P (t|t − 1) + 1
Finalmente, se observan los resultados de nuestro filtro de Kalman en la Figura 4.7.

Problema del circuito RC Error circuito RC

3 0.1
Medición
Voltaje [V]
2 Estimación
ex
0
1
0 -0.1
0 5 10 15 20 0 5 10 15 20
Tiempo [s] Tiempo [s]
Estimación de Salida Error de Salida
10 10
Medición
5 Estimación
ey
0
ŷ
-5 -10
0 5 10 15 20 0 5 10 15 20
×10 -3 Ganancia ×10 -3 Update
1 2
K
0.5 1
0 0
0 5 10 15 20 0 5 10 15 20
Figura 4.7: Filtro de Kalman del problema del circuito RC

70 FILTRO DE KALMAN
1. T. K. Moon, W.C. Stirling, “The Kalman Filter”, in Mathematical Methods and Algorithms
2. J. V. Candy, “Classical Bayesian State-Space Processors”, in Bayesian Signal Processing
Classical, Modern, and Particle Filtering Methods. John Wiley & Sons, Inc., Hoboken:
New Jersey, 2009.
3. R. F. Stengel “Optimal State Estimation”, in Optimal Control and Estimation. Dover Pu-
blications, Inc.: New York, 1994.
AGRADECIMIENTOS
Este apunte no habrı́a sido posible sin el apoyo espontáneo del profesor Marcos Orchard, ni tampoco
sin el aporte de los siguientes colaboradores:
Alejandro Veragua
Sebastián Sepúlveda
No está demás recordar que cualquier corrección, sugerencia, reclamo y/o felicitaciones se reciben en
este correo:
apunte-ssii@googlegroups.com.
71

Apunte Deteccion y Estimacion

Cargado por

Copyright:

Formatos disponibles

Apunte Deteccion y Estimacion

Cargado por

Información del documento

Descripción original:

Título original

Derechos de autor

Formatos disponibles

Compartir este documento

Compartir o incrustar documentos

Opciones para compartir

¿Le pareció útil este documento?

¿Este contenido es inapropiado?

Copyright:

Formatos disponibles

Apunte Deteccion y Estimacion

Cargado por

Copyright:

Formatos disponibles

SEÑALES Y SISTEMAS II

Apunte del Curso

Prof. Marcos Orchard

Departamento de Ingenierı́a Eléctrica

¿Cómo usar este apunte? VII

2 Análisis de Componentes Principales 19

Medición Toma de decisión

Θ ⊆ Rn Conjunto de todos los posibles valores de θ,

1.2. Definiciones básicas

En otras palabras, al elegir δ0 se acepta la hipótesis H0 (rechazando la hipótesis H1 ), y al elegir la

De lo anterior se desprenden cuatro casos posibles:

Caso 1: elegir δ0 cuando H0 es verdadera (decisión correcta)

Definición 1.2. La probabilidad de rechazar la hipotesis nula H0 cuando realmente es verdadera, se

Definición 1.3. La probabilidad de aceptar la hipótesis H1 cuando es verdadera, se conoce como la

Definición 1.4. La probabilidad de aceptar la hipotesis nula H0 cuando es H1 es verdadera, se conoce

Si agregamos estas definiciones a la Figura 1.1, se tiene:

ser una de las mejores formas de procesar datos. En efecto, considérese

la distribución de probabilidad de y está dada por:

Figura 1.3: Distribuciones de y según las hipótesis del ejemplo 1.1

β = Eθ1 [φ(y)] ≥ Eθ1 [φ0 (y)] = β 0

i Sea 0 ≤ γ ≤ 1, todo test de la forma:

es el mejor en su tamaño para probar H0 : θ = θ0 vs. H1 : θ = θ1 .

1.3.1. Test de Razón de Verosimulitud

Definición 1.6. Función de razón de verosimilitud

De la definición anterior, y el Teorema 1.2, se obtiene inmediatamente:

Corolario 1.1. Test de Log-Verosimilitud

Λ(x) = log (fX (x|θ1 )) − log (fX (x|θ0 )) > log(ν) = ν ∗

Aplicando el concepto del test de log-verosimilitud:

Calculando la razón de log-verosimilitud entre ambas distribuciones, la decisión depende de

lo que para n muestras se extrapola a

fΛ(x)|Hi umbral de decisión

Figura 1.4: Distribución de probabilidades para Λ(x).

1.4. Curvas ROC para determinar el mejor test

Figura 1.5: Curvas ROC.

Figura 1.6: Comparación de curvas ROC.

Degradación de una Batería

Figura 1.7: Regeneración de capacidad en una baterı́a de Ion-Litio.

1.4.1. Test compuestos y múltiples hipótesis

lo que se hace es “simular” datos y generar empı́ricamente la distribución asociada a la hipótesis

Figura 1.9: Caso multivariable, multiples hipótesis.

1.5. Riesgo Bayesiano

1.5.1. Clasificación por mı́nimo riesgo

Sean r̂1 , r̂2 los riesgos dados por:

r̂1 =L11 P(x|ω1 ) · P(ω1 ) + L21 P(x|ω2 ) · P(ω2 )

P(x|ω2 ) 1 L12 P(ω1 )

Luego es importante observar que de los lı́mites antiguos de X, Y se concluye que 1 ≤ u ≤ ∞ y

Ası́ se puede calcular la razón de verosimilitud:

Ası́ la implementación del test nos queda:

log(z) < η Decido H0

Bibliografı́a del capı́tulo

1. T. K. Moon, W.C. Stirling, “Introduction to Detection and Estimation and Mathematical

2. T. K. Moon, W.C. Stirling, “Detection Theory”, in Mathematical Methods and Algorithms

ANÁLISIS DE COMPONENTES PRINCIPALES

1 0.128 0.062 0.360 1.801 2.085 ... 10

3 0.485 1.323 0.756 0.014 1.476 ...

6 0.541 1.563 2.174 0.198 0.223 0

Tabla 2.1: Datos en una matriz Figura 2.1: Datos disponibles.

No respeta la correlación, pero está

En caso de que y = α1 x1 + α2 x2 + ... + αn xn + tenemos que Ê[y~x] = Ê[xxT ]~