The university of los Andes

Book

Métodos Computacionales

Manuel Alejandro Segura Delgado

PhD
ii
Contents

List of figures vii

List of Tables xi

1 Shell-scripting 3
1.1 Bash . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
1.2 Git y GitHub . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4

2 Python 7
2.0.1 Números factoriales . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
2.0.2 Máximos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
2.0.3 Sucesión de Fibonnaci . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
2.0.4 Números primos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
2.0.5 Espiral de arquı́medes . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
2.0.6 Figura de Lissajous . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
2.0.7 Fitting . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
2.0.8 Choques de duración finita . . . . . . . . . . . . . . . . . . . . . . . . . . 12
2.0.9 Cinemática . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13

3 Derivación e integración numérica 1

3.1 Algoritmo babilónico para la función raı́z cuadrada . . . . . . . . . . . . . . . . 1
3.2 Taylor Series . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
3.3 Error de redondeo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
3.4 Error de truncamiento . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
3.5 Derivación . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
3.5.1 Derivada Progresiva . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4

iii
iv Contents

3.5.2 Derivada Regresiva . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5

3.5.3 Derivada Central . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
3.5.4 Error Local . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
3.6 Segunda Derivada . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
3.7 Ejercicios: Derivación . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
3.8 Método de Newton-Raphson . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
3.8.1 Criterio de parada . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
3.9 Método de bisección . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
3.10 Ejercicios: Raı́ces de polinomios . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
3.11 Interpolación de Lagrange . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
3.11.1 Error . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
3.12 Interpolación de Newton-Gregory . . . . . . . . . . . . . . . . . . . . . . . . . . 19
3.13 Ejercicios: Interpolación de Lagrange . . . . . . . . . . . . . . . . . . . . . . . . 21
3.14 Integración . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
3.14.1 Método de trapecio simple . . . . . . . . . . . . . . . . . . . . . . . . . . 24
3.14.2 Método de trapecio compuesto . . . . . . . . . . . . . . . . . . . . . . . . 25
3.14.3 Método de Simpson simple 1/3 . . . . . . . . . . . . . . . . . . . . . . . 26
3.14.4 Método de Simpson compuesto . . . . . . . . . . . . . . . . . . . . . . . 28
3.14.5 Cuadratura Gaussiana . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
3.14.6 Cuadratura de Gaus-Legendre . . . . . . . . . . . . . . . . . . . . . . . . 32
3.14.7 Doble cuadratura . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
3.14.8 Ejercicios: Integración . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36

4 Método de MonteCarlo 45
4.1 Generación pseudo-aleatoria de números . . . . . . . . . . . . . . . . . . . . . . 45
4.2 Integración MonteCarlo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46
4.3 Ley de grandes números . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
4.4 Método de la transformada inversa . . . . . . . . . . . . . . . . . . . . . . . . . 48
4.5 Método de aceptación y rechazo . . . . . . . . . . . . . . . . . . . . . . . . . . . 49
4.6 Incertidumbre del método de MonteCarlo . . . . . . . . . . . . . . . . . . . . . . 49
4.7 Ejercicios: Método de MonteCarlo . . . . . . . . . . . . . . . . . . . . . . . . . . 50

5 Álgebra lineal 55
5.1 Método de Jacobi . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55
Contents v

5.2 Norma matricial . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57

5.2.1 Radio espectral . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58
5.3 Método de Gauss-Seidel . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59
5.4 Método de relajación sucesiva . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60
5.5 Método generalizado de Newton-Raphson . . . . . . . . . . . . . . . . . . . . . . 62
5.6 Método del descenso gradiente . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63
5.7 Propagacion del error en una red neuronal . . . . . . . . . . . . . . . . . . . . . 64
5.8 Ejercicios: Álgebra lineal . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66

6 Probabilidad 77
6.1 Definición de probabilidad . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77
6.1.1 Definición frecuentista de probabilidad . . . . . . . . . . . . . . . . . . . 79
6.2 Ejercicios: Axiomas de la probabilidad . . . . . . . . . . . . . . . . . . . . . . . 80
6.2.1 Independencia de eventos . . . . . . . . . . . . . . . . . . . . . . . . . . 81
6.2.2 Probabilidad condicional . . . . . . . . . . . . . . . . . . . . . . . . . . . 81
6.2.3 Teorema de la probabilidad total . . . . . . . . . . . . . . . . . . . . . . 82
6.3 Ejercicios: Probabilidad condicional y total . . . . . . . . . . . . . . . . . . . . . 84
6.4 Teorema de Bayes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85
6.5 Ejercicios: Teorema de Bayes . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87
6.6 Técnicas de conteo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88
6.7 Ejercicios: Técnicas de conteo . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90
6.8 Ejercicios: Generales de probabilidad . . . . . . . . . . . . . . . . . . . . . . . . 91
6.9 Distribuciones discretas de probabilidad . . . . . . . . . . . . . . . . . . . . . . 94
6.9.1 Covarianza . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96
6.9.2 Correlación . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97
6.9.3 Distribución de Bernoulli . . . . . . . . . . . . . . . . . . . . . . . . . . . 97
6.9.4 Distribución Binomial . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97
6.9.5 Distribución de Poisson . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98
6.10 Ejercicios: Distribuciones discretas de probabilidad . . . . . . . . . . . . . . . . 98
6.11 Distribuciones continuas de probabilidad . . . . . . . . . . . . . . . . . . . . . . 99
6.11.1 Distribución normal . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99
6.11.2 Distribución χ2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100
6.11.3 Distribución t-Student . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100
vi Contents

6.12 Ejercicios: Distribuciones continuas de probabilidad . . . . . . . . . . . . . . . . 100

6.13 Procesos de Markov . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102
6.14 Ejercicios: Cadenas de Markov . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105
6.15 Metropolis-Hasting algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105
6.16 Ejercicios: Metropolis-Hastings . . . . . . . . . . . . . . . . . . . . . . . . . . . 107
6.17 Hidden Markov models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111
6.18 Ejercicios: Hidden Markov models . . . . . . . . . . . . . . . . . . . . . . . . . . 114

7 Estadı́stica 117
7.1 Muestreo (sampling) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 117
7.1.1 Aleatorio (probabilı́stico) . . . . . . . . . . . . . . . . . . . . . . . . . . . 118
7.1.2 Muestreo aleatorio simple . . . . . . . . . . . . . . . . . . . . . . . . . . 119
7.1.3 Muestreo aleatorio estratificado . . . . . . . . . . . . . . . . . . . . . . . 119
7.1.4 Muestreo aleatorio sistemático . . . . . . . . . . . . . . . . . . . . . . . . 120
7.1.5 Muestreo no aleatorio . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121
7.2 Distribución de frecuencias (Histogramas) . . . . . . . . . . . . . . . . . . . . . 121
7.3 Estadı́sticos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123
7.4 Valor medio muestral . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123
7.5 Desviación estándar . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 124
7.6 Teorema del lı́mite central . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 125
7.7 Ejercicios: Muestreo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 125

8 Estimación de parámetros 129

8.1 Mı́nimos cuadrados . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 129
8.1.1 Espacio Columna . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 129
8.2 Varianza de la estimación de parámetros . . . . . . . . . . . . . . . . . . . . . . 132
8.3 Función de costo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 133
8.4 Ejercicios: Mı́nimos cuadrados . . . . . . . . . . . . . . . . . . . . . . . . . . . . 133
8.5 Función de verosimilitud (Likelihood) . . . . . . . . . . . . . . . . . . . . . . . . 137
8.6 Método de máxima verosimilitud . . . . . . . . . . . . . . . . . . . . . . . . . . 138
8.6.1 Unbinned likelihood fit . . . . . . . . . . . . . . . . . . . . . . . . . . . . 138
8.6.2 Binned likelihood fit . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 139
8.6.3 Varianza de los estimadores . . . . . . . . . . . . . . . . . . . . . . . . . 140
8.7 Ejercicios: Parameter Estimation . . . . . . . . . . . . . . . . . . . . . . . . . . 140
Contents vii

9 Confidence Intervals 145

9.1 Intervalos de confianza para la media . . . . . . . . . . . . . . . . . . . . . . . . 145
9.2 Intervalo de confianza para la diferencia de medias . . . . . . . . . . . . . . . . . 146
9.3 Intervalo de confianza para una proporción . . . . . . . . . . . . . . . . . . . . . 146
9.4 Intervalo para la varianza . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 146
9.5 Ejercicios: Intervalos de confianza . . . . . . . . . . . . . . . . . . . . . . . . . . 147

10 Pruebas de Hipótesis 151

10.1 Ejercicios: Tablas de contingencia . . . . . . . . . . . . . . . . . . . . . . . . . . 152
10.2 Prueba de hipótesis para generadores de números aleatorios . . . . . . . . . . . . 154
10.3 Ejercicios: Hypotesis testing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 155

References 163
viii Contents
List of Figures

2.1 Serie de datos y sus máximos locales. . . . . . . . . . . . . . . . . . . . . . . . . 8

2.2 Primero 20 números de Fibbonacci (izquierda), Estimación de la relación aurea
y su valor exacto (derecha). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
2.3 Primeros 1000 números primos. . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
2.4 Espiral de arquı́medes. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
2.5 Figuras de Lissajous para los desfases: δ = [0, π/4, π/2] respectivamente. . . . . 12

3.1 Campo de velocidades cerca de un cilindro sólido. . . . . . . . . . . . . . . . . . 8

3.2 Electric field of an electric charge close to grounded conducting plane. Note that
the Gauss’ Law is satisfied because of the electric field is normal in the surface.
Nice! . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
3.3 Propagación de la luz en dos medios minimizando el tiempo de viaje. . . . . . . 16
3.4 Interpolación usando el método de Newton-Gregory para el conjunto de datos
descritos en el problema. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
3.5 Domo esférico . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
3.6 Accurracy vs n-points for the Gauss-Laguerre quadrature. . . . . . . . . . . . . 40
3.7 Banda prohibida de la fase superconductora de un material con Tc ≈ 12.1278.
El valor en ∆(T = 0), es independiente del material del superconductor. . . . . . 43

4.1 Correlaciones en los primeros k-vecinos para un mal generador (izquierda) y un

buen generador (drand48) de números aleatorios. . . . . . . . . . . . . . . . . . 47
4.2 Primeros k = 10 vecinos del generador Numpy como función del número de eventos
generados. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51
4.3 Correlaciones de los primeros k = 30 vecinos del generador Numpy como función
de k-esimo vecino. Note que el valor debe fluctuar alrededor del valor teórico
C(k) = 1/4. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52

ix
x List of Figures

4.4 Generación de eventos cartesiana (izquierda) y esférica para R = 1 (derecha). . . 52

5.1 Esquema de colisión de dos discos rı́gidos. . . . . . . . . . . . . . . . . . . . . . 68

5.2 Esquema para la conservación del momentum angular. . . . . . . . . . . . . . . 69
5.3 Conservación del momento lineal total antes y después del choque. . . . . . . . . 71
5.4 Conservación del momento angular total antes y después del choque. . . . . . . 71
5.5 Variables cinemáticas los dos discos después del choque. . . . . . . . . . . . . . . 72
5.6 Estrella naciente: Puntos originales (azul), los ejes principales del cuerpo (rojo)
y puntos sobre los ejes principales (verde.) . . . . . . . . . . . . . . . . . . . . . 72
5.7 Esquema del espejo del telescopio espacial. Los vértices del cuadrado registran
las temperaturas en grados kelvin. En la parte derecha, se ilustra una rotación
para minimizar la temperatura en el punto P = (0, 0.5). . . . . . . . . . . . . . 73
5.8 Mapa de temperatura al interior del espejo del telescopio. . . . . . . . . . . . . 74

6.1 Probabilidad P(n ≤ 80) de que tengan una edad diferente. . . . . . . . . . . . . 93

6.2 Superficie de probabilidad del evento A. . . . . . . . . . . . . . . . . . . . . . . 94
6.3 Probabilidades de transición de los estados para n = 10 pasos. . . . . . . . . . . 105
6.4 Relación entre la distancia cuadrática media y el número de pasos de un proceso
Markoviano. Entender este proceso fue uno de los pasos definitvos para aceptar
la existencia de las moléculas empezando el siglo XX. . . . . . . . . . . . . . . . 108
6.5 Nube electrónica encontrada con el algoritmo de Metropolis asociada al electrón
en cada núcleo. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110
6.6 Valor esperado de la energı́a potencial eléctrica de la molécula de Hidrógeno.
Note que existe un estado ligado que conocemos como enlace covalente, la energı́a
del pozo de potencial está alrededor de Umin ≈ 2.0 eV y se da a una separación
nuclear de dequilibrium ≈ 2.1a0 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110
6.7 Valor esperado de la energı́a potencial eléctrica de la molécula de Hidrógeno
incluyendo el efecto de Heitler y London. Note que existe un estado ligado que
conocemos como enlace covalente, la energı́a del pozo de potencial está alrededor
de Umin ≈ 2.05 eV y se da a una separación nuclear de dequilibrium ≈ 2.3a0 . . . . 111
6.8 Probabilidad de cada secuencia oculta para el estado ΩO = [V, A, R]. . . . . . . . 113
6.9 Probabilidad de cada secuencia oculta para el estado observado ΩO =
[S, C, C, C, S, C, S, C]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115
List of Figures xi

7.1 Distribución de muestreo (izquierda) y distribución acumuladad de probabilidad

(derecha). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 127

8.1 Distancia ||A~x∗ −~b|| en la región rectangular definida en el problema. La distancia

mı́nima debe corresponder con el resultado del método de mı́nimos cuadrados. . 134
8.2 Gráfica de curvas de nivel de la función L(µ, ) (izquierda) y superficie (derecha). 143
8.3 Medida del equipo de 3 canales. . . . . . . . . . . . . . . . . . . . . . . . . . . . 143

10.1 Trayectoria de la partı́cula (Toy model). En este experimento la partı́cula le

tomo 33 pasos salir de la caja. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 158
10.2 Ajuste a los datos usando la distribución de Weibull exponencial. . . . . . . . . 158
10.3 Comparación entre la distribución observada y esperada para el número de pasos
de escape de la partı́cula. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 159
10.4 Distribución χ299 , el valor observado es menor que el valor crı́tico de la
distribución. p − value = 0.111, se acepta H0 . . . . . . . . . . . . . . . . . . . . 160
10.5 Distribución χ2999 , el valor observado es menor que el valor crı́tico de la
distribución. p − value = 0.190, se acepta H0 . . . . . . . . . . . . . . . . . . . . 162
xii List of Figures
List of Tables

5.1 Parámetros del cálculo. Recuerde que el momento de inercia de un sólido respecto
a su centro de masas es I = kmr2 . . . . . . . . . . . . . . . . . . . . . . . . . . . 70

6.1 Tabla de contingencia de genero y el uso de gafas. . . . . . . . . . . . . . . . . . 85

6.2 Distribución de probabilidad . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95
6.3 Distribución de probabilidad del número de microchips defectuosos. . . . . . . . 98

7.1 Distribución de estudiantes según su carrera universitaria. . . . . . . . . . . . . 119

7.2 Distribución muestral según su carrera universitaria. . . . . . . . . . . . . . . . . 120
7.3 Distribución de pesos y distribución acumulada de pesos para un conjunto de
mediciones de 100 estudiantes. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122
7.4 Distribución de salarios para 65 empleados de la empresa anónima. . . . . . . . 126

10.1 Tabla de contingencia entre artı́culos vendidos y los dı́as de la semana . . . . . . 152
10.2 Tabla de contingencia asociada a la encuesta realizada por la universidad. . . . . 153

1
2 List of Tables
Chapter 1

Shell-scripting

Section 1.1
Bash

1. Escriba un programa que calcule el factorial de un número n.

2. Escriba un verificador de parámetros desde consola:

a) Escriba una funcion de ayuda (help), que le diga a la usuario que debe ingresar 3
parámetros: posición inicial, velocidad inicial y tiempo total.

b) Use una sentencia if para verificar que el número de parámetros $# es tres. Si el

número de parámetros es diferente de tres, llamar la función de ayuda y salir (exit 1;)
de lo contrario, escribir en pantalla: corriendo programa.

3. Escriba un verificador de entradas desde consola:

a) Declare una variable (pass) y asigne el valor 0.

b) Defina una función de verificación (checkvalue), cuyo primer parámetro ($1) verifique
si es cero o uno. Si esta condición es verdadera, se asigna (pass=1). En caso contrario,
se le pide al usuario que intente de nuevo.

c) Use un ciclo while para leer una variable usando (read) y pasar ese valor como primer
parámetro a la función de verificación. El ciclo while se repite si la variable pass tiene
un valor cero, en caso contrario, el ciclo termina.

3
4 Chapter 1. Shell-scripting

4. Escriba un programa que verifique si una carpeta llamada (data) existe en su pwd. Si la
carpeta no existe crear dicha carpeta en esa ubicación actual.

5. Escriba un programa que calcule las variaciones sin repetición dadas por la siguiente
expresión:

n!
Vrn = . (1.1)
(n − r)!
donde n = 20 y r = 3. Esto puede significar el número de formas que puedo acomodar
20 partı́culas en 3 estados.

6. Cuál es la complejidad computacional del algoritmo del ejercicio anterior O(np rq )?

7. Modificar la función factorial para que reciba el número (n) como un parámetro. Usar
un for loop para calcular e imprimir en pantalla los primeros 20 números factoriales.

8. Leer un archivo de datos que tiene la ruta de 10 archivos de texto y guardar la ruta un
array de bash. Imprimir en pantalla el tercer elemento.

Section 1.2
Git y GitHub

1. Crear un video en YouTube y enviar el link describiendo las siguiente etapas de desarrollo
de un proyecto:

(a) Crear una carpeta llamada: Project

(b) Entrar a la carpeta e inicializar el repositorio local: git init
(c) Agregar un script de bash que sume los dos primeros parámetros: git add
script [Link]
(d) Realizar un commit: git commit -m "first tracking"
(e) Cambiar el nombre de la rama principal: git branch -M main
(f) Crear un el repositorio llamadado Project en GitHub
(g) Agregar el repositorio remoto: git remote add origin [Link]
(h) Subir el repositorio local: git push -u origin main
——————————————————————
1.2. Git y GitHub 5

(i) Ramificar el proyecto, creando una rama llamada Product: git branch Product
(j) Cambiar de rama Product: git checkout Product
(k) En la nueva rama crear un script de bash que multiplique los dos primeros
parámetros: git add script [Link]
(l) Realizar commit: git commit -m "script 2"
(m) Subir la rama Product al repositorio remoto: git push -u origin Product
(n) Fusionar la rama Product a la rama main: Estando en la rama principal ejecutar:
git merge Product
(o) Subir cambios a la rama main. git push origin main
(p) Eliminar la rama auxiliar Product: git branch -d Product
(q) Eliminar la rama auxiliar remota Product: git push origin --delete Product
6 Chapter 1. Shell-scripting
Chapter 2

Python

2.0.1 Números factoriales

1. Escriba una función que calcule el factorial de n, con n ∈ N.

n! = n(n − 1)(n − 2)...1 (2.1)

Calcule los primeros 20 números factoriales.

2. Escriba una función que calcule las variaciones sin repetición de n elemenos tomados de
r en r:
n!
Vrn = (2.2)
(n − r)!

a) Calcule de cuantas maneras puedo ubicar 6 carros en 3 estacionamientos. Ans: 120

3. Escriba una función que calcule las combinaciones sin repetición de n elementos tomados
de m en m, con n > m.
n n!
Cm = (2.3)
m!(n − m)!

Calcule cuantos equipos de 11 jugadores puedo formar con 22 jugadores disponibles.

Suponga que:

a) Cualquiera puede ser el arquero. Ans: 705432

b) Ya sabemos quién será el arquero. Ans: 352716

7
8 Chapter 2. Python

2.0.2 Máximos
1. Descargue los datos de: [Link]
main/MetodosComputacionalesReforma/[Link] Diseñe un algoritmo
para encontrar todos los máximos locales en esta serie de datos. La Figura [2.1] muestra
la serie de datos y los máximos locales.

1.6
1.5
1.4
1.3
1.2
1.1
1.0
0.9
3.4 3.5 3.6 3.7 3.8 3.9 4.0 4.1

Figure 2.1: Serie de datos y sus máximos locales.

2.0.3 Sucesión de Fibonnaci

1. La sucesión de Fibonnaci está definida por la siguiente definición recurrente:

f0 = 0
f1 = 1
fn = fn−1 + fn−2 (2.4)

Encuentre los primeros 20 terminos de esta sucesión.

2. Graficar la sucesión de Fibonnaci.

3. El numéro áureo está dado por:

 9

√
1+ 5
ϕ= (2.5)
2

La sucesión de Fibonnaci se relaciona con este número de la siguiente forma:

fn+1
ϕ = lim , (2.6)
n→∞ fn

Usando la sucesión de números de Fibonnaci, calcular el número aureo y comparar con

el valor exacto.

7000
Serie Fibonacci 2.0 Estimación usando la serie
Numero Aureo
6000
1.8
5000

4000 1.6

3000
1.4
2000
1.2
1000

0 1.0
0 5 10 15 20 0.0 2.5 5.0 7.5 10.0 12.5 15.0 17.5

Figure 2.2: Primero 20 números de Fibbonacci (izquierda), Estimación de la relación aurea y

su valor exacto (derecha).

2.0.4 Números primos

1. Escriba un código que calcule los primeros 1000 números primos.

2. Imprima en pantalla los primos 10 números primos.

3. Graficar los números como función de su posición, es decir, 2 es el primero, 3 es el

segundo,...
10 Chapter 2. Python

8000
7000
6000
5000
4000
3000
2000
1000
0
0 200 400 600 800 1000

Figure 2.3: Primeros 1000 números primos.

2.0.5 Espiral de arquı́medes

(a) La espiral de arquı́medes está descrita por la siguiente ecuación en coordenadas

polares:

r = a + bθ (2.7)

Haga a = 0 y b = 1 y calcule las posiciones r entre θ ∈ [0., 2π].

(b) Haga el cambio de coordenadas y gráfique la espiral.

2.0.6 Figura de Lissajous

(a) Las figuras de Lissajous se pueden obtener mediante la superposición ortogonal de

dos movimientos armónicos simples:

x = Asin(ωx t)
y = Asin(ωy t + δ) (2.8)

Podemos realizar un representación gráfica usando la siguiente relación entre las

frecuencias:
 11

15

10

0
y

10

15

15 10 5 0 5 10 15 20
x

Figure 2.4: Espiral de arquı́medes.

ωx /ωy = 1/1, 1/2, 1/3, 1/4, 1/5, 2/2, 2/3, 2/4, 2/5, 3/3, 3/4, 3/5, 4/4, 4/5, 5/5 (2.9)

Ajuste A = 1, ωx = ωy = 1 y número de lados como nsides = 5, use los siguientes

valores de desfase: δ = [0, π/4, π/2]. Note que ωx t = nx θ, donde θ ∈ [0, 2π] y
nx = [1, 2, 3, 4, 5]; con nx < ny . Debe realizar un doble for loop donde se cambie
el valor de nx y ny . Adicionalmente, debe tener un iterador entero que agregue un
nuevo sub-plot:

fig = [Link](nsides,nsides,k)
k += 1

Quitar los ejes con:

[Link]("off")

Finalmente, se obtiene:

2.0.7 Fitting
1. Bajar los datos disponibles en [Link]
main/MetodosComputacionales/[Link]
12 Chapter 2. Python

Figure 2.5: Figuras de Lissajous para los desfases: δ = [0, π/4, π/2] respectivamente.

a) Graficar los datos para visualizar la función.

b) Hacer un ajuste a los datos para encontrar los parámetros A y B.

A
f (t) = . (2.10)
1 + e−Bt

2.0.8 Choques de duración finita

Utilice el template de la clase partı́cula vista en clase para simular una colisión de duración
finita. Para la interacción de las dos partı́culas use la siguiente fuerza de contacto
(Landau-Lifshitz volumen 7).

  3

 K R + R − |~r − ~r | n̂ si |~r − ~r | < R + R
1 2 1 2 1 2 1 2
f~(~r1 , ~r2 ) = (2.11)

 0 En otro caso.

donde K = 100 N/m3 , R1 , R2 son los radios de las esferas y ~n = cos(θ)î + sen(θ)ĵ es el
vector normal al plano de contacto. Note que c = R1 + R2 − |~r1 − ~r2 | es la deformación de las
esferas. Usar m = 1 kg y R1 = R2 = 2 m, con las siguientes condiciones iniciales:
 13

Bola 1 :

~r0 = (−15, 1)
~v0 = (10, 0)
~a0 = (0, 0) (2.12)

Bola 2 :

~r0 = (0, −1.5)

~v0 = (0, 0)
~a0 = (0, 0) (2.13)

1. Crear el método de clase que calcule la fuerza entre partı́culas en la clase Particle.
Recuerde que un cuerpo no se puede aplicar fuerza a si mismo. Cada vez que se avance
en el tiempo, debe reiniciar la fuerza y la aceleración a cero, y verificar las paredes (mirar
video en bloque-neon).

2. Qué significa fı́sicamente K.

2.0.9 Cinemática
1. La ecuación de movimiento para una caı́da libre está dada por:

d2 y
m = −mg, (2.14)
dt2
donde y representa la posición de la partı́cula, m su masa y g es la aceleración debida
a la gravedad (Use g = 9.8 m/s2 ). Esta ecuación está sujeta a las condiciones iniciales
(y0 = 0, v0 = +50 m/s), que representan la posición y velocidad iniciales respectivamente.
La solución analı́tica de esta ecuación está dada por:

1
y(t) = y0 + v0 t − gt2 . (2.15)
2
Dibuje la posición y velocidad como función del tiempo, iterando hasta que el sistema
llegue hasta la posición inicial.
14 Chapter 2. Python

Ahora asuma que el sistema está sometido a una fuerza resistiva proporcional al cuadrado
de la velocidad instantánea, tal que:

d2 y 1
m 2
= −mg + CρAv 2 , (2.16)
dt 2
donde C = 0.8 es un factor experimental asociado a la resistencia del aire, ρ =
1.225 kg/m3 es la densidad del aire y A = πR2 es el área eficaz de la esfera (tomar
inicialmente R = 0.05 m).

a) Demuestre que la solución analı́tica a esta ecuación diferencial está dada por:

1
 cosh(tanh−1 ( γv
√ 0 ))
g

y(t) = y0 + 2 ln √ , (2.17)
γ cosh(−γ gt + tanh−1 ( γv √ 0 ))
g
√   
g √ −1 γv0
v(t) = tanh − γ gt + tanh √
γ g
q
CρA
donde γ = 2m
es el parámetro de amortiguamiento, con m = 1 kg. Note que existe
γv0
una restricción cinemática en este sistema −1. < √
g
< +1.
b) Dibuje la posición y velocidad como función del tiempo en la misma gráfica del caso
ideal.
c) ¿Cuál es el tiempo de vuelo para ambos casos?
d) Estudie el comportamiento del tiempo de vuelo como función de C y R (variar los
parámetros hasta encontrar diferencias). ¿Cuál parámetro tiene mayor impacto en el
tiempo de llegada de este sistema.
Chapter 3

Derivación e integración numérica

Section 3.1
Algoritmo babilónico para la función raı́z cuadrada

El algoritmo para la aproximación de la raı́z cuadrada de cualquier número real está definido
por la siguiente formula recursiva:

f0 (x) = x
 
1 x
fn (x) = + fn−1 (x) (3.1)
2 fn−1 (x)
√
donde fn (x) es la n-ésima aproximación a x, De esta forma:

√
f∞ (x) = x. (3.2)

Section 3.2
Taylor Series

Let f (x) be a function of class C ∞ . The Taylor polynomial is defined as:

∞
X
pn (x) = an (x − x0 )n , (3.3)
n=0

where x0 is a specific point, where one makes the expansion. Using x0 = 0 (Maclaurin serie),
the representation of the f (x) is therefore given by:

1
2 Chapter 3. Derivación e integración numérica

∞
X
f (x) = an x n . (3.4)
n=0

One inside for obtaining the coefficients is driven by the value of the function and its
derivatives; Hence:

f (0) = a0 + a1 0 + a2 02 + ...
f 0 (0) = a1 + 2a2 0 + ...
00
f (0) = 2a2 + 3 × 2a3 0 + ... (3.5)

Note that, a0 = f (0), a1 = f 0 (0)/1!, a2 = f 00 (0)/2!, in general, we get:

f (n) (x)
an = (3.6)
n!
This expression is known as the Maclaurin serie of f (x).

∞
X f (n) (x)
f (x) = xn . (3.7)
n=0
n!

For expanding in any arbitrary point (a), the expression is:

∞
X f (n) (a)
f (x) = (x − a)n . (3.8)
n=0
n!

We are interesting in any point of the discretization (x0 ) and its neighborhood x0 + h, where
h << 1. The expansion to around x0 is given by:

∞
X f (n) (x0 )
f (x) = (x − x0 )n . (3.9)
n=0
n!

Therefore, the value of the function in x0 + h is:

∞ ∞
0
X f (n) (x0 ) 0 0 n
X f (n) (x0 )
f (x + h) = (x + h − x ) = (h)n (3.10)
n=0
n! n=0
n!
Which is an useful formulae. the expansion up to first order is given by:

f (x + h) = f (x) + hf 0 (x) + O(h2 ) (3.11)

3.3. Error de redondeo 3

Section 3.3
Error de redondeo

Es la diferencia entre el valor exacto de un número y la aproximación calculada debida al

redondeo. Por ejemplo, π = 3.1415926535... si se aproxima a 3.1416 el error es 7.3464 × 10−6 ,
entonces la pregunta natural es: ¿cuál es el número más pequeño que podemos aproximar
usando el computador? en otras palabras ¿cuál es el valor de  para que se cumpla 1 +  = 1.

Section 3.4
Error de truncamiento

El error de truncamiento aparece cuando se usan aproximaciones en lugar de las expresiones

exactas, en general, este tipo de error depende del tipo de aproximación que se realiza. Por
ejemplo, cuando expandimos una cierta función alrededor de un punto y despreciamos términos
de orden superior, se introducen error de truncamiento a nuestras estimaciones.

sin(x) ∼
= x + O(x3 )
x3
sin(x) ∼
= x− + O(x5 ) (3.12)
3!
tiene diferente error de truncamiento para la estimación de la función sin(x).
4 Chapter 3. Derivación e integración numérica

Section 3.5
Derivación

Para construir la derivada numérica se define la siguiente discretización de puntos nodales

equi-espaciados.

xj = x0 + jh, (3.13)

donde h se denomina paso, siendo una variación ”pequeña” de la función. El conjunto de

puntos nodales con sus respectivas imágenes se denomina conjunto soporte de la discretización:

Ω = {(x0 , f (x0 )), (x1 , f (x1 )), ..., (xn , f (xn ))} (3.14)

En el caso de más dimensiones, el conjunto soporte debe estar bien definido:

Ω2 = {(x0 , y0 , f (x0 , y0 )), (x1 , y1 , f (x1 , y1 )), ..., (xn , yn , f (xn , yn ))} (3.15)

3.5.1 Derivada Progresiva

Dada esta condición podemos hacer una expansión en series de Taylor de f (x) (el dominio son
los puntos nodales).

h2 00
f (x + h) = f (x) + hf 0 (x) + f (x) + ... para h << 1 (3.16)
2
despejando la primera derivada tenemos:

f (x + h) − f (x) h 00
f 0 (x) = − f (x) (3.17)
h |2 {z }
O(h)

para algún punto de la partición:

f (xj+1 ) − f (xj )
f 0 (xj ) ∼
= (3.18)
h
La última expresión se denomina la derivada progresiva del punto xj , la cuál tiene orden
O(h) en la estimación.
3.5. Derivación 5

3.5.2 Derivada Regresiva

Para obtener la derivada regresiva se expande:

h2 00
f (x − h) = f (x) − hf 0 (x) + f (x) + ... para h << 1 (3.19)
2
despejando la primera derivada tenemos:

f (x) − f (x − h) h 00
f 0 (x) = + f (x) (3.20)
h |2 {z }
O(h)

para algún punto de la partición:

f (xj ) − f (xj−1 )
f 0 (xj ) ∼
= (3.21)
h
Para la derivada regresiva se tiene un orden de aproximación de orden O(h).

3.5.3 Derivada Central

Para estimar la derivada central se compara las expresiones de ambos desarrollos de Taylor.

h2 00 h3
f (x + h) = f (x) + hf 0 (x) + f (x) + f 000 (x) + ...
2 3!
2 3
h h
f (x − h) = f (x) − hf 0 (x) + f 00 (x) − f 000 (x) + ... (3.22)
2 3!

Restamos las dos expresiones tenemos:

f (x + h) − f (x − h) h2 000
f 0 (x) = − f (x) (3.23)
2h |3 {z }
O(h2 )

para algún punto de la partición:

f (xj+1 ) − f (xj−1 )
f 0 (xj ) ∼
= . (3.24)
2h
Notar que la estimación central tiene un orden O(h2 ) en la estimación.
6 Chapter 3. Derivación e integración numérica

3.5.4 Error Local

Una medida del error local es la distancia (euclidiana) entre el valor estimado y el valor real.

σ(Dn f (xj )) = |f n (xj ) − Dexact

n
f (xj )| (3.25)

Section 3.6
Segunda Derivada

Para estimar la segunda derivada numérica, se suma los dos desarrollos de Taylor en la
Ecuación (3.22).

f (x + h) − 2f (x) + f (x − h)
f 00 (x) = + O(h2 ) (3.26)
h2
para algún punto de la partición:

f (xj+1 ) − 2f (xj ) + f (xj−1 )

f 00 (xj ) ∼
= (3.27)
h2
Notar que la estimación tiene un orden O(h2 ) en la estimación.

Section 3.7
Ejercicios: Derivación

1. (Theoretical) Demuestre la formula alternativa para la estimación de la segunda derivada

discreta:

d2 f (xi ) f (xi+2 ) − 2f (xi ) + f (xi−2 )

2
= (3.28)
dx 4h2
2. Usando la definición de derivada central (con h = 0.05) estime la derivada de la función:

1
f (x) = √ , (3.29)
1 + e−x2

a) En el intervalo −10 ≤ x ≤ 10.

b) Para el intervalo anterior, estimar el error en cada punto nodal.

3.7. Ejercicios: Derivación 7

3. (Machine Learning analogy) La convolución discreta entre dos señales y[m] y x[m] a
tiempo discreto está dada por:

∞
X
y[m] ∗ x[m] = h[n]x[m − n], (3.30)
n=−∞

donde h[n] es una función de entrada y se denomina kernel de convolución. El operador

derivada central se puede definir a través de una operación de convolución con el siguiente
kernel K:

K = [1, 0, −1] (3.31)

De modo que el operador primera derivada central queda definido por:

∞
1 1 X
D f (xm ) = K[n + 1]f (xm−n ). (3.32)
2h n=−∞

Como el kernel tiene dimensión 3, la suma corre sobre los valores n = −1, 0, 1. No
obstante, en Python las listas tienen como primer valor n = 0, por esa razón el kernel se
translada a n + 1. Implemente este algoritmo para calcular la derivada de la función del
punto 2).

4. Diseñe un kernel de convolución para definir el operador segunda derivada central.

a) Dar la expresión matemática: D2 f (xm ) =?.

b) Implementar el cálculo de esta derivada a la función del punto 2).

5. (Theoretical) Show that the D4 f operator is given by:

f (xj+2 ) − 4f (xj+1 ) + 6f (xj ) − 4f (xj−1 ) + f (xj−2 )

D4 f (xj ) ∼
= (3.33)
h4

For this operator, what is the order (O(hk )) of the approximation?

6. Calcular el campo de velocidades cerca de la superficie de un cilindro de radio R = 2 cm.

Para esta tarea realizar los siguientes pasos:
8 Chapter 3. Derivación e integración numérica

a) Definir una discretización en los ejes x e y, donde la región es: A ∈ [−4, 4] con 25
puntos en cada eje.

b) Definir la función potencial del flujo dada por:

R2
 
φ(x, y) = V x 1 − (3.34)
x2 + y 2

donde V = 2 cm/s

c) Calcule y guarde adecuadamente el campo de velocidades usando la definición de

derivada parcial central como:

∂φ
vx =
∂x
∂φ
vy = − (3.35)
∂y

use h = 0.001. Note que al interior del cilindro el campo de velocidades debe ser igual
a cero.

d) Dibuje el campo de velocidades usando el método: [Link](x[i],y[j],Vx[i,j],Vy[i,j]).

Deberı́a obtener algo como 3.2:

1
y[cm]

4 3 2 1 0 1 2 3 4
x[cm]

Figure 3.1: Campo de velocidades cerca de un cilindro sólido.

3.7. Ejercicios: Derivación 9

7. A classical electrodynamics problem is the following: suppose a point charge q held a

distance d above an infinite grounded conducting plane. The electric potential is given
by:

 
q 1 1
V (x, y) = p −p (3.36)
4π0 (x − r~q [0])2 + (y − d)2 (x − r~q [0])2 + (y + d)2

For this numerical estimation the position of the charge is r~q = (0.51, 0.21) m. Note that
d = r~q [1].

q
a) Calculate the electric field in arbitrary units where 4π0
= 1.

~ = −∇V (x, y)
E (3.37)

The derivation region is R ∈ [0., 1.] × [0., 1.] m2 and the step of discretization is
h = 0.05.
b) Plot the vectorial field using: [Link](x[i],y[j],Ex[i,j],Ey[i,j])
The result looks like this 3.2:

1.0 Conducting plane

Electric Charge

0.8

0.6
y[m]

0.4

0.2

0.0
0.2 0.4 0.6 0.8 1.0
x[m]

Figure 3.2: Electric field of an electric charge close to grounded conducting plane. Note that
the Gauss’ Law is satisfied because of the electric field is normal in the surface. Nice!

8. Es posible construir una aproximación de orden O(h2 ) para la derivada progresiva. Para
tal propósito, se escribe el polinomio de interpolación de grado 2 para el conjunto soporte
10 Chapter 3. Derivación e integración numérica

Ω = {(x0 , f (x0 )), (x1 , f (x1 )), (x2 , f (x2 ))}, y posteriormente se calcula la derivada de este
polinomio.

a) Calcular analı́ticamente el polinomio que interpola el conjunto soporte.

b) Derivar el polinomio interpolador para encontrar la derivada en el punto x0 :

1
f 0 (x0 ) ≈ p0 (x0 ) =(−3f (x0 ) + 4f (x1 ) − f (x2 )). (3.38)
2h
Si la discretización es equidistante, tenemos:

1
f 0 (x) ∼
= (−3f (x) + 4f (x + h) − f (x + 2h)). (3.39)
2h
p
c) (Python) Para f (x) = tan(x) estimar la derivada progresiva de orden O(h2 )
(expresión anterior) en el intervalo [0.1, 1.1] con h = 0.01.
p
d) (Python) Para f (x) = tan(x) estimar la derivada central de orden O(h2 ) en el
intervalo [0.1, 1.1] con h = 0.01.
e) Calcule analı́ticamente la derivada de la función f (x), y grafique con la estimación
central y progresiva de orden O(h2 ).
f) Grafique el error nodal para ambas aproximaciones. ¿Tienen efectivamente el mismo
orden de precisión ambos resultados?

9. (SymPy) El operador derivada central de orden superior se puede ver como una
combinación lineal de las imágenes del conjunto soporte:

n
df (x) ∼ X
= ci+n f (x + ih). (3.40)
dx i=−n

donde n representa el número de vecinos del punto de la estimación, por ejemplo, para
n = 2 se tiene: X = {−2h, −h, 0, h, 2h}. La primera derivada es entonces:

df (x) ∼
= c0 f (x − 2h) + c1 f (x − h) + c2 f (x) + c3 f (x + h) + c4 f (x + 2h). (3.41)
dx
El procedimiento de interpolación permite calcular el valor de los coeficientes y establecer
el orden de aproximación de la derivada (O(hk )). Los coeficientes del operador están dados
por las derivadas de las funciones cardinales evaluadas en el punto central X[2] = 0:
3.7. Ejercicios: Derivación 11


dLi (x) 
ci = . (3.42)
dx x=0

Para cinco puntos nodales los coeficientes son:

 
1 2 2 1
~c = + , − , 0, + , − (3.43)
12h 3h 3h 12h

Usando elementos de álgebra lineal pueden ver: [Link]

En esta parte del curso, para encontrar los coeficientes les propongo la siguiente estrategia:

(a) Crear un sı́mbolo para evaluar la base cardinal:

x = [Link](’x’,real=True)

(b) Crear un sı́mbolo para el paso de derivación:

h = [Link](’h’,real=True)

(c) Crear una lista de sı́mbolos:

X=[-2*h,-1*h,0*h,1*h,2*h]

(d) Use la función Lagrange de la clase de interpolación para calcular la base cardinal
de cada punto. Debe cambiar en la lı́nea prod=1.0 a prod=1 para que no aproxime
los coeficientes a decimales.

(e) Cree una función que calcule los coeficientes y los guarde en otra lista. La función
debe tener como entrada el sı́mbolo x, la lista de simbolos X y el punto de evaluación
p.
def GetCoefficients(x,p,X):
Coefficients = []

(f) En la función, haga un for loop encontrando la base cardinal de cada punto de la
discretización.
for i in range(len(X)):
Li = Lagrange(x,X,i)

(g) Luego calcule la derivada:

dLi = [Link](Li,x,1)
12 Chapter 3. Derivación e integración numérica

(h) Calcule el coeficiente como la derivada evaluada en el punto p = 2.

C=[Link](x,X[p])
(i) Guarde el coeficiente.
[Link](C)
(j) Al iterar todos los puntos de la discretización retorne la lista de coeficientes.
(k) Deberı́a obtener el resultado deseado, Ecuación (3.43).
(l) Cuál es el orden de aproximación O(hk )?
Hint: Recuerde el error asociado al polinomio de interpolación de orden n.
3.8. Método de Newton-Raphson 13

Section 3.8
Método de Newton-Raphson

Es un método iterativo para encontrar las raı́ces reales polinomios usando conceptos de cálculo
diferencial. Tomemos un punto cualquiera xn y construimos la ecuación de la recta usando la
derivada de f (x) en xn .

f (xn+1 ) − f (xn )
Df (xn ) = (3.44)
xn+1 − xn
Se pretende que el siguiente punto en la iteración sea un raı́z de f (xn+1 ) = 0, por tanto:

f (xn )
xn+1 = xn − (3.45)
f 0 (xn )
Este es conocida como la forma iterativa de Newton-Raphson. Otro camino para deducir
esta formula consiste en expandir f (x) alrededor de xn .

(x − xn )2 00
f (x) = f (xn ) + f 0 (xn )(x − xn ) + f (xn ) + ... (3.46)
2!
Si se trunca la función hasta orden O(x2 ) y se evalúa en el siguiente punto xn+1 , el cuál se
considera un raı́z de f (x); se llega a la formula deseada.

3.8.1 Criterio de parada

Podemos usar el error relativo en cada iteración para detener el método.

|xk+1 − xk |
= (3.47)
|xk+1 |
el cual detiene el método cuando sea menor a una tolerancia, i.e,  < 10−6 .

Section 3.9
Método de bisección

Este método consiste en obtener una estimación de la raı́ces de una función, la cuales se suponen
existe en cierto intervalo [a,b]. La existencia de las raı́ces (el teorema no establece cuántas raı́ces
existen en dicho intervalo) se garantizan por medio del teorema de Bolzano, que establece lo
siguiente:
14 Chapter 3. Derivación e integración numérica

Sea f : [a, b] → R una función continua en [a,b] tal que f (a) < 0 < f (b). Entonces, existe
al menos un punto c ∈ [a, b] tal que f (c) = 0.

Si se cumple el teorema de Bolzano en nuestro intervalo de estudio, podemos acercarnos a

la raı́z de la solución dividiendo el intervalo en su punto medio iterativamente.

a+b
xm = (3.48)
2
De los dos intervalos elegimos el intervalo que cumpla la condición de Bolzano: f (a)f (xm ) <
0 o f (xm )f (b) < 0. Repitiendo el proceso hasta que el intervalo sea lo suficientemente pequeño
para acotar la raı́z solución dentro de una vecindad arbitrariamente pequeña. La noción de
distancia está dada por:

|xk − xk−1 | < k (3.49)

El algoritmo de bisección se puede describir como sigue:

a) En el intervalo [a,b] se debe cumplir la condición de Boltzano: f (a)f (b) < 0. De lo contrario,
no existe la raı́z en dicho intervalo.

b) De Divide el intervalo en su punto medio, donde tenemos dos intervalos: [a, xm ] y [xm , b].

c) Si f (xm ) = 0, el punto medio es la raı́z, con un cierto error y la iteración termina.

d) Si f (a)f (xm ) < 0. la raı́z está en [a, xm ]. Hacer b = xm y repetir paso b). De lo contrario,
la raı́z está en [xm , b]. Hacer a = xm y repetir el paso b).

Dado que en cada paso la longitud del intervalo se divide por dos, la longitud del intervalo
b−a
después de k pasos se reduce a 2k
. La aproximación a la raı́z será el punto medio de dicho
intervalo, de este modo, el error en la aproximación está dada por:

1b−a
|k | < . (3.50)
2 2k
|k | es asignado arbitrariamente. Para cierta precisión se tiene la siguiente expresión del
número de pasos que acotan la solución en esa vecindad (|xk − xexacta | < k ):

ln((b − a)/k )
k> −1 (3.51)
ln(2)
3.10. Ejercicios: Raı́ces de polinomios 15

Section 3.10
Ejercicios: Raı́ces de polinomios

1. ¿De qué tipo es el error asociado a la estimación de raı́ces usando el método de

Newton-Raphson?

2. ¿Cómo ajustar la precisión para estimar raı́ces con el método de Newton-Raphson?

3. Calcular todas las raı́ces reales de:

f (x) = 3x5 + 5x4 − x3

(3.52)

4. (SymPy) Calcular todas las raı́ces reales de los primeros 20 polinomios de Legendre. La
formula de rodrigues es:

1 dn 2
pn (x) = (x − 1)n . (3.53)
2n n! dxn
El intervalo donde existen las raı́ces es: −1 ≤ x ≤ 1.

5. (SymPy) Calcular todas las raı́ces reales de los primeros 20 polinomios de Laguerre. La
formula de rodrigues es:

ex dn −x n
Ln (x) = (e x ). (3.54)
n! dxn
El intervalo donde existen las raı́ces es: 0 ≤ x ≤ ∞.

6. Un transmisor de un barco requiere enviar un pulso de luz a un submarino que se encuentra

sumergido. El transmisor de encuentra en la posición T~ = (−3, 2) m, el submarino se
~ = (2, −2) m y la superficie del agua de encuentra en el plano
encuentra en la posición R
y = 0. Como ustedes saben la luz se propaga a velocidad constante en el vacı́o y cambia
de velocidad al cambiar de medio. El comportamiento de la refracción está definido por
la ley de Snell:

n0 sin(α0 ) = n1 sin(α1 ) (3.55)

16 Chapter 3. Derivación e integración numérica

donde n es el ı́ndice de refracción del medio, y es la relación entre la velocidad de la luz

en el vació y la velocidad de la luz en el medio: n = c/v. Esta caracterı́stica óptica hace
imposible mandar el pulso de luz en lı́nea recta entre el transmisor y el emisor.

En su parcial deben lograr la comunicación con alta precisión, no pueden usar la

ley de refracción dado que no se tiene el ángulo de incidencia (α0 ) del pulso de luz.
Adicionalmente, la ley de Snell no involucra directamente los puntos donde están el barco
y el submarino. No obstante, la luz es muy bondadosa minimizando el tiempo de viaje
(principio de Fermat). Les propongo la siguiente estrategia:

Air
Transmitter Water

0
y[m]

Receiver
x[m]

Figure 3.3: Propagación de la luz en dos medios minimizando el tiempo de viaje.

(a) (Theoretical) Muestre que el tiempo de viaje del pulso de luz está dado por:

p p
t(x) = n0 (x − T [0])2 + T [1]2 + n1 (x − R[0])2 + R[1]2 (3.56)

donde x es el punto donde se debe apuntar el láser para alcanzar el receptor. Use
n0 = 1 para el aire y n1 = 1.33 para el agua.

(b) Implemente el tiempo de viaje en Python y dibújelo para algunos valores de x. Antes
de nada, identifique cualitativamente el lugar donde se minimiza el tiempo de viaje
del pulso.
3.10. Ejercicios: Raı́ces de polinomios 17

(c) Encuentre el mı́nimo del tiempo usando el método de Newton-Raphson con una
precisión de  = 1 × 10−6 . Recuerde que la condición del mı́nimo es:

dt(x) 
=0 (3.57)
dx x=xmin

Hint: Deberá implementar la primera y segunda derivada centrales (h = 1 × 10−5 ).

(d) Con la información vectorial del punto de ingreso del láser (xmin , 0), estime el ángulo
de incidencia α0 y ángulo de refracción α1 respecto al eje y (Figura 3.3).
(e) Se cumple la ley de Snell?
Hint: Puede calcular:

n0 sin(α1 )
= (3.58)
n1 sin(α0 )
18 Chapter 3. Derivación e integración numérica

Section 3.11
Interpolación de Lagrange

Teorema 1. (Teorema de aproximación de Weierstrass). Sea f ∈ [a, b], dado  > 0 existe un
polinomio p(x) definido en [a, b] tal que:

|f (x) − p(x)| <  ∀x ∈ [a, b]. (3.59)

Descubierto por Edwarg Waring en 1779 y redescubierto por Leonhard Euler en 1783, fue
publicado por Lagrange en 1795. Sea el conjunto soporte Ω = {(x0 , y0 ), (x1 , y1 ), ..., (xn , yn )} de
n + 1 puntos diferentes, existe un polinomio interpolador de grado n:

n
X
pn (x) = f (xi )Li (x), (3.60)
i=0

donde Li (x) es la base de Lagrange (también conocidas como funciones cardinales).

n
Y x − xj
Li (x) = (3.61)
j=0,j6=i
xi − xj

Este polinomio cumple que p(xk ) = yk para todo k en {0, ..., n}.

Ejemplo:
Encontrar las funciones cardinales (Li (x)) y el polinomio interpolador para el siguiente conjunto:
Ω = {(5, 10), (10, 15)}.

x − 10 1
L0 (x) = = − (x − 10) (3.62)
5 − 10 5

x−5 1
L1 (x) = = (x − 5) (3.63)
10 − 5 5
Por tanto, el polinomio interpolador es:

p1 (x) = 10L0 (x) + 15L1 (x)

p1 (x) = x + 5 (3.64)
3.12. Interpolación de Newton-Gregory 19

3.11.1 Error
Sea f : I → R, {xi }ni=0 ⊆ I, xi 6= xj para i 6= j y suponemos que f es derivable n + 1 veces. El
error asociado a la interpolación está dado por:

f n+1 (ξx )
E = f (x) − p(x) = (x − x0 )(x − x1 )...(x − xn ) (3.65)
(n + 1)!
donde p(x) es el polinomio interpolador en {xi }ni=0 y ξx ∈ al intervalo que contiene los puntos.

proof:
Si x es un punto xk la identidad se satisface para cualquier ξ. De lo contrario, si x es fijo y
diferente xk se considera una función auxiliar:

f (x) − p(x)
F (t) = f (t) − p(t) − cL(t), donde c= (3.66)
L(x)
Qn
L(x) = i=0 (x−xi ). Si evaluamos la función auxiliar en los puntos xk , F (xk ) = yk −yk −0 =
0 para todo k. Por tanto, F tiene n+1 ceros. Adicionalmente F (x) = f (x) − p(x) − cL(x) = 0,
dada la definición de c. entonces la función F tiene n + 2 ceros en el intervalo I. Ahora, por
el teorema de Rolle, F 0 debe tener al menos n+1 ceros en el intervalo que tiene a los puntos
xk ; entonces la (n+1)-ésima derivada debe tener al menos un cero. Sea ξx ese cero. Entonces
derivamos (n + 1) veces y evaluamos en ξx :

F n+1 (ξx ) = f n+1 (ξx ) − c(n + 1)! = 0 (3.67)

debido a que la (n+1)-ésima derivada de p(x) es cero. Entonces:

f n+1 (ξx ) 1
c= → cL(x) = f (x) − p(x) = f n+1 (ξx )L(x)  (3.68)
(n + 1)! (n + 1)!

Section 3.12
Interpolación de Newton-Gregory

Sean los n + 1 puntos de la discretización x0 , x1 , ..., xn , tal que xi − xi−1 = h y sus respectivas
imágenes: f (x0 ), f (x1 ), ..., f (xn ). Se propone el polinomio de interpolación:

p(x) = a0 + a1 (x − x0 ) + a2 (x − x0 )(x − x1 ) + ... + an (x − x0 )(x − x1 )...(x − xn−1 ) (3.69)

20 Chapter 3. Derivación e integración numérica

Se requiere el cálculo de los coeficientes del este polinomio. Se establece la notación f (x0 ) =
f0 , evaluando el polinomio en x0 :

p(x0 ) = f0 (3.70)
Entonces a0 = f0 : evaluando en x1 :

p(x1 ) = f1
f1 = f0 + a1 (x1 − x0 )
f1 = f0 + a1 h (3.71)

Note que a1 es la derivada derecha. Evaluando en x2 :

p(x2 ) = f2
f2 = f0 + a1 (x2 − x0 ) + a2 (x2 − x0 )(x2 − x1 )
f1 − f0
f2 = f0 + (2h) + a2 (2h)(h) (3.72)
h
Entonces a2 se relaciona con la segunda derivada central:

f2 − 2f1 + f0
a2 = (3.73)
2h2
Pensemos en un término general para los coeficientes:

∆n f0
an = (3.74)
n!hn
donde ∆n+1 fi = ∆n fi+1 − ∆n fi , que se conoce como las diferencias finitas de la función.
Por ejemplo: podemos ver los primeros términos:

∆0+1 fi = ∆0 fi+1 − ∆0 fi
= fi+1 − fi (3.75)

∆1+1 fi = ∆1 fi+1 − ∆1 fi
= fi+2 − fi+1 − (fi+1 − fi )
= fi+2 − 2fi+1 + fi (3.76)
3.13. Ejercicios: Interpolación de Lagrange 21

∆2+1 fi = ∆2 fi+1 − ∆2 fi
= fi+3 − 2fi+2 + fi+1 − (fi+2 − 2fi+1 + fi )
= fi+3 − 3fi+2 + 3fi+1 − fi (3.77)

∆3+1 fi = ∆3 fi+1 − ∆3 fi
= fi+4 − 3fi+3 + 3fi+2 − fi+1 − (fi+3 − 3fi+2 + 3fi+1 − fi )
= fi+4 − 4fi+3 + 6fi+2 − 4fi+1 + fi (3.78)

El cálculo de estas diferencias resulta más evidente en el marco de matrices. Por ejemplo:
 
∆0 ∆1 ∆2 ∆3 ∆4 ∆5 ∆6
 
 39 0 0 0 0 0 0
 
−20
 
 19 0 0 0 0 0
 
−21 −40 −20 0 0 0 0
∆= (3.79)
 

−57 −36 4 24 0 0 0
 
−65 −8 28 24 0 0 0
 
 
−21 44 52 24 0 0 0
 
99 120 76 24 0 0 0
En este caso, el polinomio interpolador será de orden tres.

Section 3.13
Ejercicios: Interpolación de Lagrange

1. (Theoretical) Demuestre que el polinomio interpolador es único.

2. (Theoretical) Compruebe que las funciones cardinales son base (i.e, Li (x) = δij para
cada j ∈ {0, 1, ..., n}).
√
3. ¿Con qué grado de exactitud podemos calcular 114 mediante la interpolación de de
√
Lagrange para la función f (x) = x, si elegimos los puntos x0 = 100, x1 = 121, x2 = 144.
Rpta: |E| ' 1.8 × 10−3 .
22 Chapter 3. Derivación e integración numérica

4. En el lanzamiento de una bala, una cámara fotográfica registra las siguientes posiciones
en metros respecto al arma homicida (tome ~g = −9.8 m/s2 ĵ):
[Link]
[Link]
Estime el vector velocidad inicial, que estarı́a definido por la magnitud y dirección. Rpta:
V0 = 10 m/s y θ = 20◦ . Hint: Encuentre el termino lineal y cuadrático de la interpolación
y compare con la ecuación de trayectoria de la bala.

5. Interpolación Newton-Gregory: Para el siguiente conjunto de puntos:

[Link]
[Link]
Encuentre el polinomio interpolante de menor grado usando el método Newton-Gregory.

p(x) = a0 + a1 (x − x0 ) + a2 (x − x0 )(x − x1 ) + ... + an (x − x0 )(x − x1 )...(x − xn−1 ) (3.80)

Recuerde que las pendientes están definidas de forma recursiva:

f0 (xi ) = términos de la secuencia

f0 (x1 ) − f0 (x0 )
f1 (x0 , x1 ) =
x1 − x0
f1 (x1 , x2 ) − f1 (x0 , x1 )
f2 (x0 , x1 , x2 ) =
x2 − x0
fi−1 (x1 , ..., xi−1 , xi ) − fi−1 (x0 , x1 , ..., xi−1 )
fi (x0 , x1 , ..., xi−1 , xi ) = (3.81)
xi − x0
Usando las pendientes, el polinomio de grado n queda definido por:

p0 (x) = f0 (x0 )
p1 (x) = p0 (x) + f1 (x0 , x1 )(x − x0 )
p2 (x) = p1 (x) + f2 (x0 , x1 , x2 )(x − x0 )(x − x1 )
i−1
Y
pi (x) = pi−1 (x) + fi (x0 , x1 , ..., xi−1 , xi ) (x − xj ) (3.82)
j=0

El resultado se muestra en la Figura [3.4]:

3.13. Ejercicios: Interpolación de Lagrange 23

10 Interpolacion
Data
0

10

20

30

0 1 2 3 4 5 6

Figure 3.4: Interpolación usando el método de Newton-Gregory para el conjunto de datos

descritos en el problema.
24 Chapter 3. Derivación e integración numérica

Section 3.14
Integración

Para el calculo de integrales definidas existe una familia de métodos denominados Métodos de
Newton-Côtes, los cuales se basan en cambiar el integrando f (x) a un polinomio interpolador
que aproxima a f (x) en el intervalo de integración. El grado del polinomio interpolador queda
definido por el número de puntos a considerar, por ejemplo, en los casos más simples de
interpolación lı́neal (Regla de trapecio) e interpolación cuadrática (Regla de simpson), se tiene
expresiones sencillas que son fáciles de implementar computacionalmente.

3.14.1 Método de trapecio simple

Para la integral definida:
Z b
I= f (x)dx, (3.83)
a

el método de trapecios simple cambia el integrando por un polinomio interpolador de grado

uno. Este polinomio interpolador esta definido en el conjunto Ω = {(a, f (a)), (b, f (b))} que
finalmente conduce a la siguiente aproximación:

x−b x−a
f (x) ≈ p1 (x) = f (a) + f (b), ∀x ∈ [a, b]. (3.84)
a−b b−a
De modo que la integral tiene un valor aproximado:

b b
b−a
Z Z
I= f (x)dx ∼
= p1 (x)dx = (f (a) + f (b)) (3.85)
a a 2
El error en la estimación está asociado al procedimiento de interpolación. Suponiendo que
f (x) es continua y derivable de clase C 2 en el intervalo [a,b]:

f (x) = p1 (x) + (x), (3.86)

donde

f 00 (ξ)
(x) = (x − a)(x − b), a ≤ ξ ≤ b. (3.87)
2
Entonces el error asociado a la integración por le método del trapecio está dado por (h =
b − a):
3.14. Integración 25

b
h3 00
Z
E= (x)dx = − f (ξ) (3.88)
a 12
De esta forma el error alcanza un valor máximo para algún valor de ξ, donde la segunda
derivada de f (x) se maximice; de modo que podemos acotar el error máximo en esta estimación.

h3
|E| ≤ max |f 00 (ξ)|. (3.89)
12 |{z}
a≤ξ≤b

Notar que el error es de orden O(h3 ).

3.14.2 Método de trapecio compuesto

Para la integral definida:
Z b
I= f (x)dx, (3.90)
a

el método de trapecios compuesto genera una partición equi-espaciada tal que: xi+1 −xi = h,
∀i = [1, ..., n] sobre el conjunto Ω = {(x0 , f (x0 )), (x1 , f (x1 )), ..., (xn , f (xn ))}. Las condiciones
de borde corresponden con los lı́mites de integración (x0 = a y xn = b), definiendo el paso de
b−a
integración h = n
. La integral se puede escribir como:
Z b Z x1 Z x2 Z xn
f (x)dx = f (x)dx + f (x)dx + ... + f (x)dx, (3.91)
a x0 x1 xn−1

aplicando el método de trapecios simple se tiene:

Z b
h h h
f (x)dx ≈ (f (x0 ) + f (x1 )) + (f (x1 ) + f (x2 )) + ... + (f (xn−1 ) + f (xn )), (3.92)
a 2 2 2

tenemos la expresión para la regla de trapecio compuesta:

Z b n−1
h X
f (x)dx ≈ (f (a) + f (b)) + h f (xi ). (3.93)
a 2 i=1

La estimación del error se calcula sumando los errores en cada sub-intervalo.

n
X h2 00
E= Ei = − (f (ξ1 ) + f 00 (ξ2 ) + ... + f 00 (ξn )) (3.94)
i=1
12
26 Chapter 3. Derivación e integración numérica

El error puede ser acotado por el valor máximo promedio que tome la segunda derivada en
el intervalo [a, b].

h2 (b − a)
|E| ≤ max |f 00 (ξ)|. (3.95)
12 |{z}
a≤ξ≤b

Notar que el error es de orden O(h2 ).

3.14.3 Método de Simpson simple 1/3

Para la integral definida:
Z b
I= f (x)dx, (3.96)
a

el método de Simpson simple cambia el integrando por un polinomio interpolador

de grado dos. Este polinomio interpolador esta definido en el conjunto Ω =
a+b
{(a, f (a)), (xm , f (xm )), (b, f (b))}, donde xm = 2
es el punto medio del intervalo. La
interpolación conduce a la siguiente aproximación:

(x − b)(x − xm ) (x − a)(x − b) (x − a)(x − xm )

f (x) ≈ p2 (x) = f (a)+ f (xm )+ f (b), ∀x ∈ [a, b].
(a − b)(a − xm ) (xm − a)(xm − b) (b − a)(b − xm )
(3.97)
De modo que la integral tiene un valor aproximado:
Z b Z b
h
f (x)dx ∼
= p2 (x)dx = (f (a) + 4f (xm ) + f (b)). (3.98)
a a 3
Adicionalmente, es importante mencionar que la discretización debe ser par. El error de la
aproximación tiene una caracterı́stica interesante, dado que el error de la interpolación genera
un resultado idénticamente nulo. Suponiendo que f (x) es continua y derivable de clase C 3 en
el intervalo [a,b]:

f (x) = p2 (x) + (x), (3.99)

donde

b b
f 000 (ξ)
Z Z
E= (x)dx = (x − a)(x − b)(x − (a + b)/2)dx = 0, a ≤ ξ ≤ b. (3.100)
a a 4!
3.14. Integración 27

Esto significa que la regla de Simpson es exacta a orden O(h3 ), entonces necesitamos otro
camino para calcular el error de la estimación a orden O(h4 ). Para tal propósito, se va a
expandir en series de Taylor alrededor del punto medio xm a f (x) y al polinomio interpolador
p2 (x). Suponiendo que f (x) es continua y derivable de clase C 4 en el intervalo [a,b], a orden
O(h4 ) se tiene:

f 00 (xm ) f 000 (xm )

f (x) = f (xm ) + f 0 (xm )(x − xm ) + (x − xm )2 + (x − xm )3 + (x), (3.101)
2! 3!
donde el error es de orden O(h4 ).

f (4) (ξ)
(x) = (x − xm )4 . (3.102)
4!
Para reemplazar en la regla simple de Simpson, se debe encontrar f (a) y f(b), entonces:

0 f 00 (xm ) 2 f 000 (xm )

f (a) = f (xm − h) = f (xm ) + f (xm )(−h) + (−h) + (−h)3 + (xm − h)
2! 3!
00 000
f (x m ) f (x m)
f (b) = f (xm + h) = f (xm ) + f 0 (xm )(+h) + (+h)2 + (+h)3 + (xm + h)
2! 3!
(3.103)

Por tanto, la regla de Simpson queda expresada como (Notar como el término de orden
O(h3 ) no contribuye):

h h
(f (a) + 4f (xm ) + f (b)) = (6f (xm ) + f 00 (xm )h2 + (xm + h) + (xm − h))
3 3
h 1
= (6f (xm ) + f 00 (xm )h2 + f (4) (ξ)h4 ) (3.104)
3 12
Ahora integrando el desarrollo de Taylor dado por la Ecuación (3.101):

b
f 00 (xn ) 3 f (4) (ξ) 5
Z
f (x)dx = 2hf (xm ) + h + h. (3.105)
a 3 60
De modo que es posible calcular el error como la diferencia en las dos integrales.

Z b Z b
E = f (x)dx − p2 (x)dx
a a
(4) (4)
f (ξ) 5 f (ξ) 5 1
= h − h = − f (4) (ξ)h5 (3.106)
60 36 90
28 Chapter 3. Derivación e integración numérica

El error puede ser acotado por el valor máximo promedio que tome la cuarta derivada en el
intervalo [a, b]; notar que el error es de orden O(h5 ).

1 5
|E| ≤ h max |f (4) (ξ)|. (3.107)
90 |{z}
a≤ξ≤b

3.14.4 Método de Simpson compuesto

Para la integral definida:
Z b
I= f (x)dx, (3.108)
a

el método de Simpson compuesto genera una partición equi-espaciada tal que: xi+1 −xi = h,
∀i = [1, ..., n] sobre el conjunto Ω = {(x0 , f (x0 )), (x1 , f (x1 )), ..., (xn , f (xn ))}. Las condiciones
de borde corresponden con los lı́mites de integración (x0 = a y xn = b). Definiendo el paso de
b−a
integración como h = n
, la integral se puede escribir:
Z b Z x2 Z x4 Z b
f (x)dx = f (x)dx + f (x)dx + ... + f (x)dx (3.109)
a a x2 xn−2

Notar que el número de puntos de la partición deber ser par y que los puntos medios serı́an
el sub-conjunto {x2n+1 }. Aplicando el método simple a cada integral tenemos:

Z b  
h
f (x)dx ≈ f (a)+4f (x1 )+f (x2 )+f (x2 )+4f (x3 )+f (x4 )+...+f (xn−2 )+4f (xn−1 )+f (b)
a 3
(3.110)
En general, podemos escribir la expresión final para el método de Simpson compuesto:

Z b  n−1 n−2 
h X X
f (x)dx ≈ f (a) + 4 f (xi ) + 2 f (xi ) + f (b) (3.111)
a 3 i=1,impares i=2,pares

Para la estimación del error total, se debe considerar el error cometido en cada una de las
integrales. De este modo:
 
1 5 (4) (4) (4)
|E| = h f (ξ1 ) + f (ξ2 ) + ... + f (ξn/2 ). (3.112)
90
Finalmente, acotando superiormente la estimación tenemos el error de la regla de Simpson
compuesta a orden O(h4 ).
3.14. Integración 29

b−a 4
|E| ≤ h max |f (4) (ξ)|. (3.113)
180 |{z}
a≤ξ≤b

Cabe resaltar que la regla de Simpson compuesta que aproxima el integrando por un
polinomio de orden tres, reduce en dos ordenes el error cometido por método de trapecio
compuesto.

Ejemplo:
Para la integral:
Z 1 p
1 + e−x2 dx, (3.114)
−1

Usando la regla de Simpson compuesta, ¿cuál debe ser el número de puntos para aproximar
la integral con un error menor a 0.001?

Usando el operador discreto D4 f (xj ) se estima el máximo en M = 3.183. Entonces:

r
4 (b − a)5
n≥ M ≈ 4.87 = 6 Intervalos (3.115)
E180
debido que el número de intervalos debe ser par!

3.14.5 Cuadratura Gaussiana

La idea es aproximar la integral a una suma ponderada de la función f (x), que use los valores
más óptimos posibles para su determinación. En general, la formula de cuadratura está dada
por:

n Z 1
∼
X (n) (n)
GM (f ) = wi f (xi ) = f (x)dx. (3.116)
i=1 −1

(n) (n)
Donde {wi } son los pesos de ponderación y los {xi } son los puntos de Gauss, ambos
cantidades dependen del número de puntos que se elijan. A diferencia de los métodos anteriores,
queremos encontrar una regla que nos permita maximizar el grado de precisión en nuestras
estimaciones. En general, para calcular los coeficientes se exige que la regla sea exacta para el
polinomio de mayor grado que se tenga. Tenemos entonces:
30 Chapter 3. Derivación e integración numérica

Z 1 n
X (n) (n)
xk dx = wi (xi )k k = 0, 1, ..., N (3.117)
−1 i=1

donde N es tan grande como sea posible. Necesitamos 2n condiciones para integrar
polinomios de grado 2n − 1, esto se denomina la regla de n puntos, para el grado polinomial
2n − 1.

(n)
Lema 1: Si N = 2n, entonces no existe el conjunto {wi } tal que se verifique la igualdad.

Qn n 2
Supongamos que existe dicho polinomio, si definimos L(x) = j=1 (x − xj ) entonces
R1 P (n) (n)
L(x) ≥ 0 y por tanto la integral −1 L(x)dx > 0. Sin embargo, wi L(xi ) = 0, entonces
hay una contradicción.

(n) (n)
Lema 2: Sea {wi } los pesos y {xi } los puntos de Gaus tal que se cumpla la
Ecuación (3.117) para k = 0, 1, 2, ..., N = 2n − 1, los pesos deben satisfacer:

Z 1
(n) (n)
wi = Li (x)dx, (3.118)
−1

con

n (n)
(n)
Y x − xk
Li (x) = (n) (n)
(3.119)
k=1,k6=i xi − xk
(n)
Las bases cardinales son de grado ≤ 2n − 1 y Li (xj ) = δij , como la regla de cuadratura
debe ser exacta para los polinomios de grado 2n − 1, entonces:

Z 1 n
X n
X
(n) (n) (n) (n) (n)
Li (x)dx = wj Li (xj ) = wj δij = wi (3.120)
−1 j=1 j=1

Qn (n)
Ahora, tenemos un polinomio de grado n, Pn (x) = i=1 (x − xi ) y supongamos que
tenemos un polinomio cualquier Q(x) de grado ≤ n − 1, de modo que la multiplicación tiene
grado ≤ 2n − 1. Por tanto la identidad de la cuadratura se debe satisfacer.

Z 1 n
X (n) (n)
Pn (x)Q(x)dx = wi Pn (xni )Q(xi ). (3.121)
−1 i=1
3.14. Integración 31

Por la definición de Pn (x) la integral es nula y los polinomios deben formar una base
ortonormal. Adicionalmente, el conjunto {xni } deben ser las raı́ces de los Pn (x).

En particular, tomemos una regla de cuadratura con 2 puntos nodales:

Z 1
pn (x)dx = w0 pn (x0 ) + w1 pn (x1 ). (3.122)
−1

Adicionalmente, se toma una base genérica b = {1, x, x2 , x3 }, imponiendo que la regla de

cuadratura sea exacta hasta el grado polinomial más alto.

Z 1
x0 dx = w0 + w1
−1
Z 1
x1 dx = w0 x0 + w1 x1
−1
Z 1
x2 dx = w0 x20 + w1 x21
−1
Z 1
x3 dx = w0 x30 + w1 x31 (3.123)
−1

Lo que conduce al siguiente sistema de ecuaciones no-lineales:

2 = w0 + w1 (3.124)
0 = w0 x0 + w1 x1 (3.125)
2/3 = w0 x20 + w1 x21 (3.126)
0 = w0 x30 + w1 x31 (3.127)

De las ecuaciones homogéneas, adecuadamente divida la Ecuación (3.127) por la

Ecuación (3.126) para encontrar que x0 = −x1 . De la Ecuación (3.125) muestre que w0 = w1 ,
con la Ecuación (3.126) muestre que w0 = w1 = 1, y finalmente, muestre que la ecuación (3.127)
√ √
conduce a x0 = −1/ 3 y x1 = 1/ 3. Con los anteriores resultados, la regla de cuadratura es:
Z 1 √ √
pn (x)dx = pn (−1/ 3) + pn (+1/ 3), n = 0, 1, 2, 3. (3.128)
−1

Ejemplo:
32 Chapter 3. Derivación e integración numérica

Obtener la regla de cuadratura para de f (x) en el intervalo [−3, 3] usando la siguiente

partición P = {−1, 0, 1}.

Calculando las funciones de Lagrange:

1
L1 (x) = x(x − 1)
2
L2 (x) = −(x + 1)(x − 1)
1
L3 (x) = x(x + 1) (3.129)
2
Integrando dichas funciones en el intervalo [−3, 3] obtenemos el conjunto de pesos:

wi3 = {9, −12, 9} (3.130)

Entonces la regla de cuadratura es:

Z 3
f (x)dx ≈ 3[3f (−1) − 4f (0) + 3f (1)]. (3.131)
−3

Note que es exacta para f (x) = x2 .

3.14.6 Cuadratura de Gaus-Legendre

De la necesidad del uso de polinomios ortonormales y el cálculo de sus raı́ces, surge la cuadratura
de Gaus-Legendre. Este calculo usa los polinomios de Legendre como conjunto ortonormal.

1 dn 2
pn (x) = (x − 1)n , (3.132)
2n n! dxn
Los primeros polinomios son:

p0 (x) = 1
p1 (x) = x
1
p2 (x) = (3x2 − 1)
2
1
p3 (x) = (5x3 − 3x) (3.133)
2
cumplen con la siguiente relación de completez:
3.14. Integración 33

(
1
0 m 6= n
Z
pn (x)pm (x)dx = 2
(3.134)
−1 2n+1
m=n
y las siguientes relaciones de recurrencia.

(n + 1)pn+1 (x) = (2n + 1)xpn (x) − npn−1 (x). (3.135)

dpn
(x2 − 1)= nxpn − npn−1 . (3.136)
dx
Adicionalmente se tiene el siguiente coeficiente director.

(2n)!
an = (3.137)
2n (n!)2
Para calcular los pesos de la cuadratura de Gauss-Legendre, se requiere escribir las funciones
cardinales en términos de los polinomios de Legendre. Si definimos las funciones base de
Lagrange como:

pn (x)
L(x) = (3.138)
x − xk
Tenemos una forma indeterminada cuando x = xk . Para conocer dicho lı́mite, si existe,
usemos la regla de L’Hôpital.

dpn (x) 
= p0n (xk )
dx

lim L(x) = 
d(x−xk ) 
(3.139)
x→xi
dx x=xk

como las funciones de Lagrange deben ser base y valor 1 justo en cada raı́z, se tiene la
siguiente expresión.

1 pn (x)
L(x) = (3.140)
p0n (xk ) x − xk
Por tanto, los pesos pueden ser evaluados por:
Z 1
1 pn (x)
wk = 0 dx (3.141)
pn (xk ) −1 x − xk
Ejemplo:

Calcular los pesos de ponderación para dos puntos (n = 2). El polinomio de Legendre serı́a
p2 (x) = 12 (3x2 − 1), cuyas raı́ces son: x0 = √1 ,
3
x1 = − √13 . Los pesos están dados por:
34 Chapter 3. Derivación e integración numérica

1 1
(3x2 − 1)
Z
1 2
w0 = dx = 1. (3.142)
3( √13 ) −1 x− √1
3

1 1
(3x2 − 1)
Z
1 2
w1 = dx = 1. (3.143)
3(− √13 ) −1 x+ √1
3

Por tanto, la regla de cuadratura para dos puntos es:

Z 1
1 1
f (x)dx = f ( √ ) + f (− √ ) (3.144)
−1 3 3
para integrales que tiene un lı́mite de integración diferente al de los polinomios de Legendre,
es posible hacer un cambio de variable de modo que -1 coincida con el lı́mite inferior y 1 con el
lı́mite superior de la integral. Planteando la siguiente proporción:

x−a t − (−1)
= , (3.145)
b−a 1 − (−1)
de manera que:

1
x =[t(b − a) + a + b]
2
1
dx = (b − a) (3.146)
2
Usando este cambio de variable, podemos escribir la regla de cuadratura para cualquier
intervalo de integración.
Z b n
1 X 1
f (x)dx = (b − a) wk f ( [tk (b − a) + a + b]) (3.147)
a 2 k=0
2
La Ecuación (3.141) aunque es correcta, resulta poco útil para calcular los pesos a alto orden
computacionalmente. Es posible encontrar una versión integrada de dicha expresión usando el
siguiente teorema:

Teorema 2. (Identidad de Christoffel-Darboux) Sea {pn(x) }∞

n=0 el sistema ortonormal de
polinomios respecto a un peso ω(x). Entonces, para todo x, y ∈ R, x 6= y y n ≥ 1 se cumple.

n−1  
X an−1 pn (x)pn−1 (y) − pn−1 (x)pn (y)
pk (x)pk (y) = (3.148)
k=0
an x−y
siendo pn (x) = an xn + ... y an > 0
3.14. Integración 35

Usando esta identidad y haciendo y = xk las raı́ces de los polinomios tenemos.

n−1
X an−1 pn (x)pn−1 (xk )
pk (x)pk (xk ) = (3.149)
k=0
an x − xk
Integrando, tenemos:
Z 1
an−1 pn (x)
1= pn−1 (xk ) dx (3.150)
an −1 x − xk
Por tanto, podemos escribir la Ecuación (3.141) como:

an
wk = (3.151)
an−1 pn−1 (xk )p0n (xk )
an
Un cálculo sencillo muestra que an−1
= 2/n y usando la Ecuación (3.136).

2
wk = , k = 1, ..., n (3.152)
(1 − x2k )[p0n (xk )]2
Esta última expresión permite calcular los pesos usando librerias como Simpy.

3.14.7 Doble cuadratura

En el caso de integración doble:
Z d Z b
I= f (x, y)dxdy (3.153)
c a
Tenemos un extensión de la regla de cuadratura denominada doble cuadratura. Recordemos
el cambio de variable para llevar la integral al intervalo [−1, 1].

d 1
b−a
Z Z
1
I= f ( (t1 (b − a) + a + b), y)dt1 dy (3.154)
c 2 −1 2
Aplicando el cambio de variable nuevamente para la variable y.

1 1
(b − a)(d − c)
Z Z
1 1
I= f ( (t1 (b − a) + a + b), (t2 (d − c) + c + d))dt1 dt2 . (3.155)
4 −1 −1 2 2
Aplicando la regla de cuadratura de n puntos, tenemos:

n n
(b − a)(d − c) X X (n) (n) 1 (n) 1 (n)
I∼
= wi wj f ( (xi (b − a) + a + b), (xj (d − c) + c + d)). (3.156)
4 i=1 j=1
2 2
36 Chapter 3. Derivación e integración numérica

En el caso del cáscaron esférico, los limites de integración son variables. En efecto, x2 =
R1
R2 − y 2 y dA = 2πxdx. El volumen del cáscaron será 0 A(y)dy. En términos de integrales
dobles tenemos:

Z 1 Z √R2 −y2
2
2πxdxdy = π, R = 1. (3.157)
0 0 3
Esto involucra una dificultad adicional para aplicar la regla de cuadrátura, dado que:
q
(n) (n)
b(xi ) = R2 − (xi )2 . (3.158)

Por inspección, usted puede comprobar que para lı́mites variables la cuadratura está dada
por:

 n X
X n
d−c 1 (n) 1 (n)
I∼
(n) (n) (n)
= wi (b(xi )−a)wj f ( (xi (b−a)+a+b), (xj (d−c)+c+d)). (3.159)
4 j=1 i=1
2 2

3.14.8 Ejercicios: Integración

1. (Theoretical) Hacer pasos intermedios para regla de trapecio simple, Ecuación (3.85).

2. (Theoretical) Encontrar el error para regla de trapecio simple, Ecuación (3.88).

3. (Theoretical) Hacer los pasos intermedios para encontrar la regla de Simpson simple,
Ecuación (3.98).

4. (Theoretical) Verificar el resultado presentado en la Ecuación (3.100).

5. Resolver la siguiente integral con el método del trapecio:

Z 1
2
e−x dx (3.160)
0

usando la condición que el error debe ser menor a 0.005 en el cálculo (si el número de
puntos (n) es decimal tomar el mayor entero).

6. Para encontrar la inductancia exacta de una bobina toroidal es necesario calcular la

integral:
3.14. Integración 37

√
Z a
a2 − x 2 √
dx = π(R − R2 − a2 ), (3.161)
−a R+x

donde R = 0.5 cm es el radio al centro del toriode y a = 0.01 cm es el radio de la

sección transversal del toriode. Estimar la integral con el método del trapecio y la regla
de Simpson 1/3, con un error menor al 0.5%. Realizar la integral es un verdadero reto
analı́tico.

7. Una forma de generalizar el método de integración del trapecio para una integral doble
de una función f (x, y) consiste en dividir el plano xy en un grilla de cuadrados iguales y
calcular el promedio del valor de la función de cada uno de los 4 vértices de cada cuadrado
pequeño de la grilla. Calcule numéricamente el volumen de una semiesfera [3.5] de radio
R = 1 como sigue:

a) Cree una grilla entre −R y R en el plano xy, donde el número de cuadrados en cada
lado de la grilla sea n. Es decir, la grilla tendrı́a n + 1 puntos en cada eje, y n2
cuadrados pequeños.

b) Para cada cuadrado pequeño calculo el promedio de la función en los cuatro vértices
y multiplique por el área del cuadrado pequeño. Si el punto queda fuera de la esfera
asuma que el valor de la función f (x, y) es cero.

1.0
0.8
0.6
0.4
0.2
0.0
1.00
0.75
0.50
0.25
1.000.75 0.00
0.500.25 0.25
0.000.25 0.50
0.500.75 0.75
1.00 1.00

Figure 3.5: Domo esférico

38 Chapter 3. Derivación e integración numérica

8. Implemente la regla de doble cuadratura dada por la Ecuación (3.159) para solucionar el
problema del cáscaron esferico.

9. (Sympy) La regla de Simpson 3/8 consiste en aproximar el integrando por un polinomio

interpolador de orden 3.

a) Encontrar las funciones cardinales de dicha interpolación e integrar para demostrar

que:
Z b  
∼ 3h 2a + b a + 2b
f (x) = f (a) + 3f ( ) + 3f ( ) + f (b) . (3.162)
a 8 3 3
El número de puntos de la discretización debe ser múltiplo de tres.
Hint: Use la siguiente discretización {0, h, 2h, 3h}.
2a+b a+2b
b) Dado que h = (b − a)/3, muestre que los puntos intermedios son 3
y 3
respectivamente.

10. (Sympy) Muestre que el error asociado a la regla de Simpson 3/8 simple está dado por:

b
f (4) (ξ)
Z
3 5 (4)
E= (x − x0 )(x − x1 )(x − x2 )(x − x3 )dx = − h f (ξ). (3.163)
4! a 80

Hint: Considere la siguiente integral:

Z 3h
I= (x)(x − h)(x − 2h)(x − 3h)dx (3.164)
0

11. (Theoretical) Muestre que para la regla de Simpson 3/8 compuesta, el error puede ser
acotado por:

(b − a) 4 (4)
|E| ≤ |{z} |f (ξ)|.
h max (3.165)
80
a≤ξ≤b

12. Using the Simpson’s 3/8 rule, calculate the following integral:

Z 1 p
I= 1 + e−x2 dx ≈ 2.6388571169 (3.166)
−1
3.14. Integración 39

13. Evaluar:

Z 2
1
dx, (3.167)
1 x2

por medio de la cuadratura de Gauss-Legendre con dos y tres puntos. Rpta: I2 =

0.497041, I3 = 0, 499874.

10
14. Escribir el polinomio p(x) = 3 + 5x + x2 en la base de Legendre. Rpta: p(x) = 3 0
p (x) +
5p1 (x) + 32 p2 (x)

15. (Sympy) Dada la aproximación de cuadratura gausiana:

Z 1 n
X
f (x) = wk f (xk ), (3.168)
−1 k=0

donde w0 , w1 , ...wn son los coeficientes ponderados ó pesos.

(a) Halle los ceros de los primeros 20 polinomios de Legendre.

(b) Halle los pesos de ponderación para los primeros 20 polinomios de Legendre.

16. Estime la siguiente integral usando el método de cuadratura de Gaus-Legendre:

Z ∞
1
≈ 1.110721 (3.169)
0 x4 +1

Hint: Dividir la integral para tener dos integrales con lı́mites [-1,1] y [0,1].

17. In the black-body radiation problem the following integral appears:

∞
x3 π4
Z
dx = . (3.170)
0 ex − 1 15

a) Compute this integral using the Gauss-Laguerre quadrature method for n=3 evaluation
points.

b) For this estimation, plot the relative error (r (n) = Iestimated (n)/Iexact ) as a function
of the evaluation points, with n = [2, 3, ..., 10] [3.6].
40 Chapter 3. Derivación e integración numérica

Hint: For the Gauss-Laguerre method, the weights are given by:

xk
wk = (3.171)
(n + 1) [Ln+1 (xk )]2
2

1.000

0.998

0.996

0.994
(n)

0.992

0.990

0.988 Laguerre quadrature acurracy

2 3 4 5 6 7 8 9 10
n

Figure 3.6: Accurracy vs n-points for the Gauss-Laguerre quadrature.

18. (Sympy) La cuadratura de Gaus-Hermite está definida para integrales de la forma:

Z ∞
2
I= f (x)e−x dx (3.172)
−∞

que tiene la siguiente representación en cuadraturas:

N
I∼
X
= wi f (xi ) (3.173)
i=1

donde los puntos xi son las raı́ces de los polinomios de Hermite dado por la formula de
Rodrigues:

2 dn −ξ2
Hn (ξ) = (−1)n eξ e (3.174)
dξ n
La formula de los pesos de Gaus-Hermite está dada por:

√
2n−1 n! π
wi = 2 (3.175)
n Hn−1 (xi )2
3.14. Integración 41

(a) Encontrar los primeros 20 ceros de los polinomios y los correspondientes pesos de la
cuadratura.
(b) El estado de un oscilador armónico en mecánica cuántica está dado por la funciones
de probabilidad:

1 mω 1/4 −ξ2 /2
ψn (ξ) = √ ( ) e Hn (ξ) (3.176)
2 n! π~
n
p mω p mω
donde ξ = ~
x. Haga ~
= 1, es decir, ξ = x para la aplicación numérica.
Estime numéricamente el valor cuadrático medio de la posición de la partı́cula en el
primer estado exitado (n = 1). El valor exacto de la integral está dado por:
Z ∞
2
< x >= |ψ1 (x)|2 x2 dx = 3/2 (3.177)
−∞

El polinomio de Hermite de orden 1 está dado por:

H1 (x) = 2x (3.178)

19. Superconductividad: BSC(Bardeen,Schrieffer,Cooper) Temperatura crı́tica

En este modelo de superconductividad a temperaturas cerca del cero absoluto, surge una
ecuación denominada la ecuación de banda prohibida.

p
~ωD
2 + ∆(T )2 /2kB T )
Z
1 1 tanh(
= pk dk (3.179)
N0 V 2 −~ωD 2k + ∆(T )2
donde las constantes definen un sistema termodinámico a nivel cuántico y ∆(T ) define
el valor de la banda prohibida (Note que la función no tiene anti-derivada y requiere
un procedimiento completamente numérico). Definamos la temperatura de Debye como
~ωD
TD = kB
= 300 K y N0 V = 0.3. Por inspección, la integral toma la forma:

1
p
tanh( x2 + ∆0 (T )2 TD /2T )
Z
1 1
= p dx (3.180)
N0 V 2 −1 x2 + ∆0 (T )2
Note que ∆0 (T ) = ~ωD ∆(T ) (banda prohibida primada). Existe una temperatura a la
cuál la banda prohibida tiende a cero ∆(T ) → 0, dicha temperatura es la temperatura
crı́tica Tc donde el material se convierte en superconductor. Para estimar Tc de este
material les propongo la siguiente estrategia:
42 Chapter 3. Derivación e integración numérica

(a) Implemente la función a integrar en Python. Esta función debe tener la temperatura
T y el tamaño de la banda prohibida ∆0 (T ) como parámetros.

(b) Cargar los puntos y los pesos de la cuadratura de Gaus-Legendre a orden n = 20.

(c) Debe variar la temperatura del material entre 1 a 20 y hacer el cálculo de la integral.
La temperatura crı́tica será justamente cuando la integral sea aproximadamente igual
a:

1
≈ Integral(T = Tc , ∆0 (Tc ) = 0). (3.181)
N0 V
(d) En ese sentido deberá usar algún criterio de parada del algoritmo. Por ejemplo:
if [Link](I-1/(N0V)) < 1e-3: return Tc. La temperatura crı́tica es: Tc ≈
12.1278 K (Pag 97) [1].

20. Superconductividad 2: BSC(Bardeen,Schrieffer,Cooper) banda prohibida

superconductora
En este problema se calculará la banda prohibida como función de la temperatura en la
fase superconductora. Usando el problema anterior, implemente la siguiente rutina:

(a) Debe variar los parámetros T y ∆0 (T ). La temperatura debe variar como:

t = [Link](0.05,20.,50) y la banda prohibida primada como:
delta = [Link](0.,1.,10000)
Para cada valor de temperatura, debe encontrar el valor de la banda prohibida
primada donde la integral sea aproximadamente igual a:

1
≈ Integral(T = t, ∆0 (T ) = delta). (3.182)
N0 V
(b) Guarde los valores de temperatura y banda como: T = t/Tc y ∆(T ) = delta(TD /Tc ).
Este cálculo se debe al cambio de variable en la integral y a la forma de graficar la
banda prohibida.

(c) Dibuje la banda prohibida ∆(T )/kb Tc como función de la temperatura crı́tica
T /Tc . Si guardo correctamente los datos del numeral anterior, deberı́a obtener algo
como (Pag 96) [1]:

(d) Estime el valor de la banda prohibida en el cero absoluto:

3.14. Integración 43

1.75
1.50
1.25

(T)/kbTc
1.00
0.75
0.50
0.25
0.00
0.2 0.4 0.6 0.8 1.0
T/Tc
Figure 3.7: Banda prohibida de la fase superconductora de un material con Tc ≈ 12.1278. El
valor en ∆(T = 0), es independiente del material del superconductor.

∆T (0)/kb Tc ≈ 1.763 (3.183)

Este valor no depende del material del que está hecho el superconductor y es uno de
los resultados más hermosos de la fı́sica del siglo XX.
44 Chapter 3. Derivación e integración numérica
Chapter 4

Método de MonteCarlo

Section 4.1
Generación pseudo-aleatoria de números

En general, no es posible generar auténticos números aleatorios por computador. Para obtener
una secuencia aparentemente aleatoria, se utilizan métodos deterministas con altos periodos de
repetición. Los métodos de congruencia lineal están definidos por:

xn+1 = (axn + c)mod m n ≥ 0, (4.1)

donde 0 ≤ a < m se denomina multiplicador, 0 ≤ c < m es el incremento y m es

el módulo del método. En la secuencia x0 se denomina la semilla del generador y debe
ser ajustada inteligentemente para obtener secuencias diferentes en cada generación. En la
expresión congruencial, el modulo representa el resto entre (axn + c) y m, Por ejemplo si
x0 = 1, a = 16, c = 0 y m = 7 tenemos la siguiente secuencia x1 = 2, x2 = 4, .... Uno de los
generadores más conocidos y usados tanto en C++ como en Python está definido por:

a = 5DEECE66D (base 16)

c = B (base 16)
m = 248 (4.2)

Para obtener un secuencia de números normalizada entre cero y uno (0 ≤ r ≤ 1), se debe
normalizar al modulo (m) la ecuación congruencial.

45
46 Chapter 4. Método de MonteCarlo

xn+1 (axn + c)mod m

r= = (4.3)
m m
Para obtener una distribución de número en otro intervalo cualquiera (A ≤ x ≤ B) se tiene
debe escalar la distribución U(0, 1).

x = A + (B − A)r. (4.4)

Por otro lado, es importante garantizar la aleatoriedad de la secuencia, para tal propósito,
se tienen algunas pruebas a la función de distribución de dicha secuencia [3], por ejemplo,
podemos usar el momento de la distribución (P (x)):

N
1 X k∼ 1 k
Z
1 1
xi = x P (x)dx ∼
= + O( √ ) (4.5)
N i=1 0 k+1 N

También se puede estudiar las correlación entre los k-vecinos cercanos, donde k es pequeño:

N
1 X
C(k) = xi xi+k , k = 1, 2, ... (4.6)
N i=1

Si la secuencia es aproximadamente uniforme podemos aproximar la suma como:

N Z 1 Z 1
1 X ∼ 1
xi xi+k = dx dyxyP (x, y) = . (4.7)
N i=1 0 0 4

De modo que, si la distribución de números pseudo-aleatorio sigue una función de

distribución uniforme, C(k) debe tener como máximo 1/4. La Figura [4.1] muestra las
correlaciones de los primeros 30 vecinos para un generador de números aleatorios malo y uno
bueno. Note que un buen generador de eventos debe mantener las correlaciones alrededor de
1/4.

Section 4.2
Integración MonteCarlo

La integración Montecarlo es un método de integración de funciones generales que se sustenta

en la ley de grandes números.
4.3. Ley de grandes números 47

0.34 0.34

0.32 0.32

0.30 0.30

0.28 0.28
C(k)

0.26 0.26

0.24 0.24

0.22 0.22

0.20 0.20
0 5 10 15 20 25 30 0 5 10 15 20 25 30
k-esimo, Mal Generador k-esimo, Generador drand48

Figure 4.1: Correlaciones en los primeros k-vecinos para un mal generador (izquierda) y un
buen generador (drand48) de números aleatorios.

Section 4.3
Ley de grandes números

Sea {xi } un conjunto de puntos muestrales (observaciones o datos simulados por computador)
y el promedio muestral de dicho conjunto.

N
1 X
x̄ = xi . (4.8)
N i=1
Al aumentar el tamaño de la muestra, el promedio muestra converge al valor esperado con
probabilidad i.

Muestral → P oblacional
x̄ → E(x)
R1
De esta manera, el operador integración 0 f (x)dx puede verse como el valor esperado de
E(f (x)) donde x ∼ U(0, 1) tiene una distribución uniforme. De este modo, tomemos una
muestra x1 , x2 , ... ∼ U(0, 1) y evaluamos en la función f (x1 ), f (x2 ), ..., entonces:

N
1 X
→ E(f (x))
f (xi ) |{z} (4.9)
N i=1
N →∞
48 Chapter 4. Método de MonteCarlo

El valor promedio de una función está dado por:

Z b
1
f¯ = f (x)dx, (4.10)
b−a a
en consecuencia:

N Z b
1 X 1
lim f (xi ) = f (x)dx. (4.11)
N →∞ N
i=1
b − a a

Otro resultado interesante es que si x es una variable aleatoria con densidad de probabilidad
f y g es una función g : R → R, entonces:

Z ∞ N
X
E(g(x)) = g(x)f (x)dx = g(xi )f (xi ). (4.12)
−∞ i=1

Section 4.4
Método de la transformada inversa

Rx
Sea x una variable aleatoria con densidad de probabilidad f , y sea F (x) = −∞
f (s)ds la función
acumulada de f . Dado que, la función acumulada de probabilidad es monotona creciente resulta
ser inyectiva. De manera que, existe un único valor en el rango tal que:

F −1 (u) = x. (4.13)

Por tanto, podemos generar una distribución uniforme y tomar la pre-imagen en la

distribución acumulada. El resultado es el siguiente:

Sea U ∼ U(0, 1) y F una distribución acumulada, entonces F es injectiva en algún intervalo

pre-imagen de [0,1]. Si X se define como F −1 (U ) = X, entonces X tiene distribución F .

proof:

Queremos ver si P (X ≤ x) = F (x), lo que significa que x tiene distribución F . Por definición
de X:

P (X ≤ x) = P (F −1 (U ) ≤ x) (4.14)
4.5. Método de aceptación y rechazo 49

Como F es monotona creciente en el intervalo [0,1] es invertible:

P (X ≤ x) = P (F (F −1 (U )) ≤ F (x))
= P (U ≤ F (x))
= F (x) (4.15)

El último paso se justifica dado que U tiene distribución uniforme U ∼ U(0, 1). El algoritmo
se resume como sigue:

1. Generar un número aleatorio que siga una distribución uniforme U ∼ U(0, 1).

2. Tomar la pre-imagen de U , X = F −1 (U ).

Section 4.5
Método de aceptación y rechazo

En la gran mayoria de los casos es imposible encontrar la función inversa para usar el método
de transformada inversa. En estos casos, se utiliza el método de aceptación y rechazo, el
cual mediante la distribución uniforme de números puede calcular el área bajo la curva de
distribuciones complicadas. Este método consiste se base en los siguientes pasos:

1. Generamos un xi que siga una distribución uniforme entre contenido en el intervalo de

integración [a,b].

2. Para xi , generamos un yi que siga una distribución uniforme entre 0 y el máximo de f(x).

3. Si yi < f (xi ) incluimos el valor xi en la lista, de otro modo, no incluimos el valor xi .

Section 4.6
Incertidumbre del método de MonteCarlo

La incertidumbre en el método de MonteCarlo se estima propagando el error cometido en cada

uno de los puntos muestrales xi . Recordemos que la media muestral está dada por:

N
1 X
x̄ = xi (4.16)
N i=1
50 Chapter 4. Método de MonteCarlo

entonces:

∂ x̄ 1
= (4.17)
∂xi N
Teniendo en cuenta que cada medición es independiente la varianza asociada a la media
muestra es:

N  2 N
X ∂ x̄ X 1 2
V ar(x̄) = δi2 = δ (4.18)
i=1
∂xi i=1
N2 i

Por tanto, la incertidumbre está dada por:

δ
V ar(x̄) ∼
p
=√ (4.19)
N
δ es despreciable con el número de eventos simulados, en general se ajusta a 1 para tener
√
una relación exacta con 1/ N .

√
En esta ecuación es importante resaltar que la precisión del Método disminuye con N,
lo cuál debe tenerse en cuenta para cada calculo realizado con esta técnica. En el cálculo del
número aureo usando la sucesión Fibonnaci aunque los puntos no provienen de un proceso
aleatorio, la precisión del método tiene una dependencia similar.

Section 4.7
Ejercicios: Método de MonteCarlo

1. Demostrar la Ecuación (4.4).

2. Programar un mal generador de números aleatorios y usar drand48 en Python para

reproducir la Figura [4.1].

3. Un test simple para probar la calidad de un generador de eventos es evaluar los momentos
de la distribución de datos, la cual difiere del valor de la distribución uniforme en un orden
de aproximación O( √1N ); esto se debe al error asociado al método de MonteCarlo.

N
1 X k∼ 1 k
Z
1 1
xi = x P (x)dx ∼
= + O( √ ) (4.20)
N i=1 0 k+1 N
4.7. Ejercicios: Método de MonteCarlo 51

a) Implemente un código que estime los primeros k = 10 momentos de la distribución

de datos generados por Numpy. Genere [Link](2,6,5) eventos y grafique el eje
x en escala logaritmica. Lo momentos de la distribución de datos de muestran en la
Figura [4.2]

0.6 k=1
k=2
k=3
k=4
k=5
0.5 k=6
k=7
k=8
k=9
k = 10

0.4
k-moment value

0.3

0.2

0.1

0.0

102 103 104 105 106

N Points

Figure 4.2: Primeros k = 10 vecinos del generador Numpy como función del número de eventos
generados.

4. Otro método para probar la calidad de un generador de eventos es evaluar las correlaciones
con los k -vecinos más cercanos, donde k ∼ 30.

N
1 X
C(k) = xi xi+k , (k = 1, 2, 3...) (4.21)
N i=1

a) Implemente un código que estime los coeficientes de correlación para los primeros k=30
vecinos, con N = 104 eventos de la distribución de datos generados por Numpy. Las
correlaciones se muestran en la Figura [4.3].

5. Escriba un código que genere una secuencia uniforme de eventos, contenidas en una esfera
de radio R = 2 como se muestra en la Figura [4.4].

6. Usando la generación de puntos sobre una esfera estime los momentos de inercia de la
esfera respecto a los ejes de rotación Ixx , Iyy , Izz . Los momentos están dados por:
52 Chapter 4. Método de MonteCarlo

Generador Numpy
0.30

0.28

0.26

C(k) 0.24

0.22

0.20
0 5 10 15 20 25 30
k-esimo vecino

Figure 4.3: Correlaciones de los primeros k = 30 vecinos del generador Numpy como función de
k-esimo vecino. Note que el valor debe fluctuar alrededor del valor teórico C(k) = 1/4.

1.0 1.0

0.8
0.5

0.6
0.0
0.4

0.5
0.2

0.0 1.0
0.0 1.0
0.2 0.5
0.4
0.6 0.2 0.0 0.0 0.5 1.0
0.8 0.6 0.4 0.5 0.5 0.0
1.01.0 0.8 1.01.0

Figure 4.4: Generación de eventos cartesiana (izquierda) y esférica para R = 1 (derecha).

Z
Ixx = (y 2 + z 2 )dxdydz
ZV
Iyy = (x2 + z 2 )dxdydz
ZV
Izz = (x2 + y 2 )dxdydz (4.22)
V
4.7. Ejercicios: Método de MonteCarlo 53

En este problema estamos asumiendo que la densidad de la esfera es ρ = 1 kg/m3 con

volumen unitario. Usando el método de MonteCarlo las integrales se pueden estimar
como (por ejemplo para Ixx ):

N
1 X
Ixx ≈ (y[i]2 + z[i]2 ) (4.23)
N i=1

Compare con el resultado exacto Iij = 2/5 para i = x, y, z con i = j. Adicionalmente,

estime el producto de inercia Ixy :

Z
Ixy = − xydxdydz. (4.24)
V

Qué puede concluir sobre la simetrı́a de rotación de la esfera?

7. Usando la generación de puntos sobre una esfera estime la siguiente integral para {x2 +
y 2 + z 2 ≤ 1}:

Z Z Z √
x2 +y 2 +z 2
e dxdydz = 4π(e − 2) (4.25)

8. Haga una gráfica entre la precisión (Iestimated /Iexact ) y el número de puntos muestrales,
el cuál debe tener un rango de valores adecuado para ver el comportamiento. ¿Qué
dependecia se observa?

9. La distribución Beta está dada por:

Γ(α + β) α−1
f (x; α, β) = x (1 − x)β−1 , 0≤x≤1 (4.26)
Γ(α)Γ(β)
donde Γ(n) = (n − 1)!. Para f(x; 2,4), halle el área bajo la curva usando el método de
aceptación y rechazo con una incertidumbre del 1%.

10. La siguiente integral multidimensional:

Z 1 Z 1 8
X 2
−7
··· 2 xi dx1 dx2 ...dx8 , (4.27)
0 0 i=1

25
tiene el valor exacto 192
, usando el método de MonteCarlo estime esta integral con tres
cifras de precisión.
54 Chapter 4. Método de MonteCarlo
Chapter 5

Álgebra lineal

Section 5.1
Método de Jacobi

El método de Jacobi es usado para resolver sistemas lineales del tipo Ax = b, el cuál construye
una sucesión descomponiendo la matriz A.

A=D+R (5.1)

donde D es la matriz diagonal y R es la suma de la matriz triangular inferior y la matriz

triangular superior.

Dx = Rx + b
x = D−1 (b − Rx) (5.2)

si aii 6= 0 tenemos un método iterativo.

 
1 X
xk+1
i = bi − k
aij xj (5.3)
aii i6=j

A nivel iterativo se necesitan todos los valores del vector x en el paso k para calcular el valor
del vector en k + 1.

Ejemplo: Sea el sistema lineal:

55
56 Chapter 5. Álgebra lineal

3x − y − z = 1
−x + 3y + z = 3
2x + y + 4z = 7 (5.4)

encontrar la ecuación iterativa del método de Jacobi.

    
1/3
1 0
1/3 0
−1
    
D b= 0 1/3 0  3 =  1 
    (5.5)
0 0 1/4 7 7/4
    
1/3 0 0 0 −1 −1 0 −1/3 −1/3
D−1 R = 
    
0 1/3 0  −1 0 1  = −1/3 0 1/3  (5.6)
    
0 0 1/4 2 1 0 1/2 1/4 0
Finalmente:
   
  
xk+1 1/3 0 −1/3 −1/3 xk
      
y k+1  =  1  − −1/3 0 1/3   k (5.7)
      y 
z k+1 7/4 1/2 1/4 0 zk
Este método se puede interpretar como una sucesión de vectores que ocurre en un punto
arbitrario del espacio vectorial al vector solución, la cuál es generada por un operador de
traslación T definido por el sistema de ecuaciones.

xk+1 = c + Txk (5.8)

A nivel computacional para detener el método se calcula una métrica entre el vector solución
y el vector calculado. Se define el residuo como:
v
u n
uX
Re = ||x − y|| = t (xi − yi )2 (5.9)
i=1

para detener el método se usa la siguiente condición:

||b − Ax|| < δ (5.10)

5.2. Norma matricial 57

Convergencia débil: El método de Jacobi converge a la solución del problema lineal si la

matriz asociada al sistema de ecuaciones es diagonal dominante.

Convergencia Fuerte: La sucesión xk+1 = Txk + c, para k ≥ 0 converge a la sucesión

única x = Tx + c si y sólo si ρ(T) < 1.

Section 5.2
Norma matricial

Una norma matricial sobre el conjunto de las matrices n × n (Mn (R)) es una aplicación:

||.|| : Mn (R) −→ R

que satisface las siguientes propiedades:

a) ||A|| ≥ 0, ∀A ∈ Mn (R).

b) ||A|| = 0, si y solo si aij = 0, 1 ≤ i, i ≤ n.

c) ||αA|| = |α|||A|| ∀α ∈ R y ∀A ∈ Mn (R).

d) ||A + B|| ≤ ||A|| + ||B||, ∀A, B ∈ Mn (R).

e) ||A · B|| ≤ ||A|| · ||B||, ∀A, B ∈ Mn (R).

Tenemos un par de ejemplos útiles:

n
X
||A||1 = max
|{z} |aij | (5.11)
1≤j≤n i=1
n
X
||A||∞ = max
|{z} |aij | (5.12)
1≤i≤n j=1

Hagamos un ejemplo especı́fico:

 
−1 0 3
 
A =  4 −2 1

 (5.13)
0 6 5
Hence,
58 Chapter 5. Álgebra lineal

||A||1 = max[5, 8, 9] = 9 (5.14)

||A||∞ = max[4, 7, 11] = 11 (5.15)

5.2.1 Radio espectral

El radio espectral de una matriz A ∈ Mn (R), con notación ρ(A), se define como:

ρ(A) = max{|λ| : λ valor propio de A} (5.16)

Recordar que los valores propios de una matriz A son las raı́ces de la ecuación caracterı́stica
det(A − λI) = 0. Note que, si el valor propio es un número complejo, digamos λ = α + iβ, con
p
α, β ∈ R, entonces |λ| = α2 + β 2 .

El radio espectral no es una norma! puesto que ρ(A) = 0 con A 6= 0

Teorema: Si A ∈ Mn (R) entonces,

ρ(A) ≤ ||A||, (5.17)

para cualquier norma matricial inducida definida por:

||Ax||
||A|| = max (5.18)
|{z} ||x||
x6=0

Elijamos el auto-vector x asociado al radio espectral de A.

||Ax|| ||λx|| |λ|||x||

ρ(A) = = = =λ (5.19)
||x|| ||x|| ||x||
cuando elijamos un valor propio diferente al valor propio mas grande tenemos:

ρ(A) ≤ ||A||. (5.20)

Necesitamos un último resultado:

Teorema: Sea A ∈ Mn (R) entonces,
5.3. Método de Gauss-Seidel 59

lim Ak = 0 si y solo si ρ(A) < 1. (5.21)

k→∞

consideremos que limk→∞ Ak = 0 y tomemos algún valor de propio de A, entonces Ak x = λk x

con x 6= 0. Entonces, necesariamente limk→∞ λk = 0; lo cual se cumple solo sı́ |λ| < 1. Como
λ es un valor propio cualquiera y por la definición de radio espectral tenemos que ρ(A) < 1.
Finalmente, para entender la convergencia fuerte definamos el error en la aproximación
k-ésima:

ek = x(k) − x (5.22)

x representa el vector solución. Por otro lado, x = T x + c de manera que:

ek+1 = xk+1 − x = Txk + c − (Tx + c) = T(xk − x) = Tek (5.23)

por lo que podemos ver el error en el paso k como la sucesión matricial por el error inicial
e = Tk e0 . Debido al teorema anterior, ek = 0 solo sı́ limk→∞ Tk = 0 y eso solo es posible si
k

ρ(T ) < 1.

Section 5.3
Método de Gauss-Seidel

EL método de Gauss-Seidel se basa en el método de Jacobi con la única diferencia que cada
variable se actualiza con cada iteración interna. Usando el ejemplo anterior se evidencia como
funciona:

1+y+z
x =
3
3+x−z
y =
3
7 − 2x − y
z =
4
(5.24)

Si la iteración empieza en el vector nulo tenemos que los valores en la siguiente iteración
son:
60 Chapter 5. Álgebra lineal

1+0+0
x = = 1/3
3
3 + 1/3 + 0
y = = 10/9
3
7 − 2/3 − 10/9
z = = 47/36
4
(5.25)

Teniendo la misma condición de convergencia del método de Jacobi, el método de

Gauss-Seidel acelera el tiempo de convergencia a la solución. En general se tiene:

1 + yk + zk
xk+1 =
3
3 + xk+1 − z k
y k+1 =
3
7 − 2xk+1 − y k+1
z k+1 =
4
(5.26)

Section 5.4
Método de relajación sucesiva

Estos métodos realizan un ponderación de la solución de un sistema lineal en una iteración

especı́fica. La expresión general para un peso w dado es:

xji = wxji + (1 − w)xj−1

i (5.27)

donde el ı́ndice i se relaciona al valor incógnita y j al valor de la iteración. Dependiendo

del valor del peso, se tiene la siguiente calificación.

• Sub-relajación 0 < w < 1. No convergente (para movimiento de fluidos puede converger).

• No modificado w = 1.

• Sobre-relación 1 < w < 2. Acelera convergencia.

• Divergente w > 2.
5.4. Método de relajación sucesiva 61

Encontrar el valor adecuado del peso dada una discretización es un asunto complicado del
álgebra matricial numérica. Por ejemplo, para una discretización cuadrada de N puntos, el
peso optimo está dado por:

2
w= (5.28)
1 + Nπ
Ejemplo:

Encontrar x13 usando el método de sobre-relajación del siguiente sistema de ecuaciones

lineales. Considere w = 1.25. y ~x0 = (0, 0, 0).

4x1 − x2 = 0
−x1 + 4x2 − x3 = 6
−x2 + 4x3 = 2 (5.29)

Usando el acercamiento por Gauss-Seidel tenemos:

2 + x02 2+0
x11 = = = 0.5
4 4
6 + x11 + x03 6 + 0.5 + 0
x12 = = = 1.625
4 4
2 + x12 2 + 1.625
x13 = = = 0.906 (5.30)
4 4
Ahora usemos sobre-relajación (para mostrar como funciona).

2 + x02 2+0
x11 = = = 0.5
4 4
x11 = (1.25)(0.5) + (1 − 1.25)(0) = 0.625
6 + x11 + x03 6 + 0.625 + 0
x12 = = = 1.656
4 4
x12 = (1.25)(1.656) + (1 − 1.25)(0) = 2.070
2 + x12 2 + 2.070
x13 = = = 1.017
4 4
x13 = (1.25)(1.017) + (1 − 1.25)(0) = 1.271 (5.31)

Aunque la solución exacta para x3 = 27/28, el método de sobre-relajación parece más

distantes de la solución exacta. Sin embargo, se deja al lector la implementación de ambas
62 Chapter 5. Álgebra lineal

rutinas para verificar la potencial del método de relajación sucesiva. Adicionalmente, es posible
imponer un criterio de convergencia del método:
 j
 xi − xj−1

i

 < (5.32)
 xji 
donde  es arbitrariamente pequeño. Podemos pensar en la siguiente factorización para la
implementación de la relajación sucesiva en ecuaciones diferenciales parciales:

uli,j = ul−1 l l−1

i,j + w[ui,j − ui,j ] (5.33)

definiendo ri,j := uli,j − ul−1

i,j se tiene:

uli,j = ul−1
i,j + wri,j (5.34)

Generalmente ri,j ≤ 10−10 .

Section 5.5
Método generalizado de Newton-Raphson

Sea el sistema de ecuaciones lineales o no lineales:

f1 (x, y) = 0
f2 (x, y) = 0 (5.35)

Suponemos que las funciones son diferenciables, de modo que podemos expandir en series
de Taylor alrededor de (x0 , y0 ) a primer orden. Tenemos:

∂f1 (x0 , y0 ) ∂f1 (x0 , y0 )

f1 (x0 , y0 ) + ∆x + ∆y ≈ 0
∂x ∂y
∂f2 (x0 , y0 ) ∂f2 (x0 , y0 )
f2 (x0 , y0 ) + ∆x + ∆y ≈ 0 (5.36)
∂x ∂y
con x = x0 + ∆x y y = y0 + ∆y. Vectorizando estas ecuaciones ~x = (x, y), tenemos:

F~ (~x0 ) + J(x0 )(~x − ~x0 ) ≈ 0 (5.37)

Donde definimos el Jacobiano del sistema:

5.6. Método del descenso gradiente 63

!
∂f1 ∂f1
∂x ∂y
J= ∂f2 ∂f2
(5.38)
∂x ∂y

Si el Jacobiano del sistema es invertible, tenemos:

~x = ~x0 − J−1 (~x0 )F~ (~x0 ) (5.39)

La solución al sistema de ecuaciones significa encontrar la sucesión de Cauchy dada por:

~xn+1 = ~xn − J−1 (~xn )F~ (~xn ). (5.40)

Section 5.6
Método del descenso gradiente

~ x):
Vamos a pensar el problema de sistemas de ecuaciones como un vector de funciones. Sea G(~
 
f1 (x1 , x2 , ..., xn )
 
 f (x , x , ..., x ) 
 2 1 2 n 
 
~ x) = 
 . 
G(~  (5.41)

 . 

 

 . 

fn (x1 , x2 , ..., xn )

Este vector tendrá una norma asociada. Se define una función a minimizar usando dicha
norma.

1 ~T
F~ (~x) = G ~ x)
(~x)G(~ (5.42)
2
La solución al sistema de ecuaciones estará definida por:

||F~ (~x)|| = ~0. (5.43)

Debemos movernos en contra del gradiente de la función vectorial.

~x1 = ~x0 − γ∇F~ (~x0 ), (5.44)

64 Chapter 5. Álgebra lineal

donde γ se denomina tasa de aprendizaje. Lo que brinda la información direccional del

gradiente de una función vectorial es la matriz jacobiana.

∂fi
Jij = (5.45)
∂xj
El jacobiano es la matriz de derivadas parciales de primer orden de las funciones del vector
~ x). Suponemos F~ : Rn −→ Rm y que es diferenciable, entonces x ∈ Rn y F~ (~x) ∈ Rm .
G(~
Finalmente, el descenso del gradiente sobre una función vectorial está dada por:

~ xn ).
~xn+1 = ~xn − γJT (~xn ) · G(~ (5.46)

Note que la función gradiente es el transpuesto del Jacobiano, debido a que el gradiente es
un vector fila definido sobre el complemento ortogonal; mientras que el vector de soluciones es
un vector columna.

Section 5.7
Propagacion del error en una red neuronal

Comencemos por las definiciones. Tenemos que el vector δ l está definido como:

∂C
δil = ,
∂zil
donde z l se define más adelante. Sea al los valores de las funciones de activación en la
l-ésima capa. Tenemos que la función de costo es C = C(aL ), ya que esta solo depende de la
función de activación en la última capa. Adicionalmente, sea z l el input de la capa l. Siendo f
la función de activación, tenemos que

al = f l (z l ), z l = W l f l−1 (z l−1 ).

Utilizamos el vector de derivadas como la derivada con respecto al vector. Esto es:
 
∂C ∂C
= .
∂z l i ∂zil
Usando la regla de la cadena, tenemos que

∂C ∂C ∂aL
= .
∂ziL ∂aL ∂ziL
5.7. Propagacion del error en una red neuronal 65

En notación de vectores, por ende, tenemos que

∂C ∂aL
δL = .
∂aL ∂z L
Para las demás capas, nótese que z l−1 solo afecta a C a través de z l . Esto quiere decir que
tenemos

∂C X ∂C ∂zjl
δil−1 = l−1 = .
∂zi j
∂zjl ∂zil−1

Ahora, tenemos que los al son, en general, vectores. alk serı́a la función de activación de
l
la k-ésima neurona de la capa l. Similarmente, tenemos que Wjk serı́a el peso que la j-ésima
neurona de la capa l le asigna a la k-ésima salida de la capa l − 1. Por ende, tenemos que

X
l l−1
zjl = Wjk ak .
k

Derivando con respecto a zil−1 :

∂zjl X
l ∂al−1
k
= Wjk .
∂zil−1 k
∂zil−1
Reemplazando en nuestra expresión para δil :

X X ∂C l−1
l ∂ak
δil−1 = W .
j k
∂zjl jk ∂zil−1

Reconocemos δjl en la suma, tal que

l−1
l ∂ak
XX
δil−1 = l
δj Wjk l−1 .
j k
∂zi
Generalizando la notación de derivadas de vectores introduciendo antes, definimos la
derivada de un vector con respecto a otro como

∂al−1 ∂al−1
 
k
= .
∂zil−1 ∂z l−1 ki

Representando δ l como una matriz y reemplazando con esta definición, tenemos que

∂al−1
XX  
δil−1 = (δ l l
)j1 Wjk .
j k
∂z l−1 ki
66 Chapter 5. Álgebra lineal

Transponemos δ l para obtener un producto matricial:

∂al−1
XX  
δil−1 = l T
((δ ) l
)1j Wjk .
j k
∂z l−1 ki

Reconociendo el producto matricial, llegamos al resultado deseado:

l−1
l−1 l T l ∂a
δ = (δ ) W .
∂z l−1

Section 5.8
Ejercicios: Álgebra lineal

1. Implemente el algoritmo de Gaus-Seidel para resolver el sistema de ecuaciones visto en

clase.

2. En Python, implemente una clase para solucionar sistemas lineales. El constructor de clase
deber recibir la matriz y el vector independiente. El primer método de clase resuelve el
sistema usando el método de Jacobi, y el segundo método resuelve el sistema usando el
método de Gaus-Seidel. En la pantalla debe aparecer el número de iteraciones necesarias
para resolver el problema.

3. Implemente un algoritmo que realice la multiplicación de dos matrices. Use el algoritmo

para calcular:

  
1 0 0 4 −2 1
  
AB = 
 5 1 0 0 3 7
  (5.47)
−2 3 1 0 0 2

4. (Theoretical) Muestre con detalle que la sustitución hacia adelante se expresa como:

i−1
X
x i = bi − Aij xj (5.48)
j=0

5. (Theoretical) Muestre con detalle que la sustitución hacia atrás se expresa como:

Pn
bi − j=i+1 Aij xj
xi = (5.49)
Aii
5.8. Ejercicios: Álgebra lineal 67

donde i = n, n − 1, ..., 0. Note que la diagonal de la matriz triangular superior puede

tener cualquier valor.

6. Usando los métodos de Newton-Raphson y descenso del gradiente, encuentre la solución

de los siguientes sistemas de ecuaciones no lineales:

ln(x21 + x22 ) − sin(x1 x2 ) = ln(2) + ln(π),

ex1 −x2 + cos(x1 x2 ) = 0. (5.50)

Use x(0) = (2, 2)

6x1 − 2cos(x2 x3 ) − 1 = 0,
q
9x2 + x21 + sin(x3 ) + 1.06 + 0.9 = 0,
60x3 + 3e−x1 x2 + 10π − 3 = 0. (5.51)

Use x(0) = (0, 0, 0)

Ejercicios tomados de Numerical Analysis - Burden and Faires, Ninth Edition, Chapter
10.

7. (Machine Learning application) La 3-esfera está definida por todos los puntos
equidistantes de un punto dado en R4 , generalmente definida por:

S 3 = {(x1 , x2 , x3 , x4 ) ∈ R4 | x21 + x22 + x23 + x24 = 1} (5.52)

Dado que no podemos visualizar la 3-esfera, se usa el método de descenso del gradiente
para encontrar sistemáticamente los puntos que pertenecen a S 3 . Posteriormente, se hace
una proyección sobre R3 . Para hacer el mapeo de S 3 se sugiere la siguiente estrategia:

a) Modifique el código de la clase para que S 3 sea el único elemento del vector de
funciones.

b) Calcule la métrica de minimización como la magnitud del vector de funciones.

c) Modifique el Jacobiano a una dimensión (1 × 4).

68 Chapter 5. Álgebra lineal

d) Genere un punto aleatorio uniformemente distribuido en R4 (es decir, un vector

aleatorio de 4 componentes) en el intervalo [-1,1].

e) Corra el algoritmo de aprendizaje para que el vector aleatorio se ubique en algún punto
de S 3 . El código aprende a mover el punto minimizando la métrica.

f) Repita el proceso de aprendizaje para al menos 103 puntos.

g) Verifique que los puntos obtenidos pertenezcan a la 3-esfera.

h) Con los puntos obtenidos de S 3 , dibuje la proyección espacial (X, Y, Z), es decir, las
primeras tres columnas. ¿Qué figura geométrica resulta?

Hint: Primero haga el caso de la 2-esfera para entender el problema. Note que la
proyección de la 2-esfera en el plano es el circulo.

8. La Figura [5.1] representa la colisión de dos discos de masa m1 y m2 , y radios r1 y r2

respectivamente. El primer disco tiene una velocidad inicial ~u1 antes del choque y el
segundo disco se encuentra en reposo ~u2 = 0. El parámetro de impacto b, que es la
distancia entre la dirección de la velocidad ~u1 y el centro del segundo disco caracteriza la
interacción.

Figure 5.1: Esquema de colisión de dos discos rı́gidos.

Los siguientes son resultados teóricos que deben deducir y justificar fı́sicamente.

a) De la conservación del momento lineal demuestre que:

m1 u1 = m1 v1x + m2 v2x (5.53)

0 = m1 v1y + m2 v2y (5.54)
5.8. Ejercicios: Álgebra lineal 69

Figure 5.2: Esquema para la conservación del momentum angular.

b) De la conservación del momento angular respecto del punto de contacto demuestre

que (ver Figura [5.2]):

−m1 r1 u1 sin(θ) = I1 ω1 − m1 r1 (v1x sen(θ) + v1y cos(θ)) (5.55)

0 = I2 ω2 + m2 r2 (v2x sen(θ) + v2y cos(θ)) (5.56)

c) De la definición del coeficiente de restitución demuestre que:

(v1x cos(θ) − v1y sen(θ)) − (v2x cos(θ) − v2y sin(θ))

e=− (5.57)
u1 cos(θ)

d) Si no se considera deslizamiento de un disco respecto al otro en el punto de contacto P .

Las velocidades de los discos en el punto P son iguales. Muestre que esta restricción
conduce a:

r1 w1 + v1x sen(θ) + v1y cos(θ) = −r2 w2 + v2x sen(θ) + v2y cos(θ) (5.58)

El sistema lineal queda descrito como:

     
sinθ cosθ −senθ −cosθ 1 1 v1x 0
     
 0 M 0 1 0 0  v1y   0 
 
   
     
 M 0 1 0 0 0   v2x   M u1 
−cosθ senθ cosθ −senθ 0 0 · 
   = 
   v2y  eu1 cosθ
 
     
 senθ cosθ
 0 0 −k 0  
  r1 ω1   u1 senθ 
  
0 0 senθ cosθ 0 k r2 ω2 0
70 Chapter 5. Álgebra lineal

Parámetros
m1 3 kg
m2 1 kg
M 3
r1 0.1 m
r2 0.2 m
k 1/2
e 0.8
~u1 2î m/s

Table 5.1: Parámetros del cálculo. Recuerde que el momento de inercia de un sólido respecto
a su centro de masas es I = kmr2 .

Para la aplicación numérica definamos los parámetros del problema:

Variar el ángulo de incidencia (0 < θ < π/2) para describir algunos estados finales
después del choque. Debe verificar que los principios de conservación se cumplan en este
fenómeno, esto significa que el momento lineal y angular se mantengan constantes. La
energı́a no se conserva porque hay restitución (e) en la colisión, si todo va bien, graficar
las seis variables cinemáticas de los discos para entender como se mueven los cuerpos
después del proceso de ”scattering”. Serı́a muy interesante realizar la animación de este
fenómeno, sin embargo, las dificultades para animar la traslación + rotación de los discos
son enormes.

Para que tengan una guı́a de los resultados, las graficas de este problema son:

9. En un instante de tiempo una estrella naciente está compuesta por un conjunto

de masas unitarias (m = 1 kg). Las posiciones de cada una de las
partı́culas están escritas en: [Link]
MetodosComputacionales/[Link]

Para un sistema de N partı́culas con masa mk y posiciones rk , el tensor de inercia se

define como:
5.8. Ejercicios: Álgebra lineal 71

1.00
6.3 Momento X inicial Momento Y inicial
Momento X final 0.75 Momento Y final
6.2
0.50
6.1 0.25
px[kgm/s]

py[kgm/s]
6.0 0.00

5.9 0.25
0.50
5.8
0.75
5.7
1.00
0 20 40 60 80 0 20 40 60 80

Figure 5.3: Conservación del momento lineal total antes y después del choque.

0.0 Momento Angular Inicial

Momento Angular Final
0.1

0.2
Lz[kgm2/s]

0.3

0.4

0.5

0.6
0 20 40 60 80

Figure 5.4: Conservación del momento angular total antes y después del choque.

N
X
I= mk ((rk · rk )E − rk ⊗ rk ) (5.59)
k=1

donde E es el tensor unitario e1 ⊗ e1 + e2 ⊗ e2 + e3 ⊗ e3

a) Muestre que el tensor de inercia para esta estrella está dado por:
 
4 0 −1
1.1163 × 10 2.0524 × 10 7.4286 × 10
 
0 3
I=
 2.0524 × 10 8.9520 × 10 −2.3483 × 103 

−1 3 3
7.4286 × 10 −2.3483 × 10 4.2341 × 10

Note que el tensor es simétrico Iij = Iji como es esperado.

b) Encuentre los auto-valores y auto-vectores. Si dos valores propios son cercanos o
iguales ¿Qué interpretación tienen los auto-valores respecto a las simetrı́as de este
sistema?
72 Chapter 5. Álgebra lineal

0 1
2

2
[rad/s]

5
0 20 40 60 80

v1x v2x
1.75 v1y 2.5 v2y
1.50 2.0
1.25 1.5
v1[m/s]

v2[m/s]
1.00 1.0
0.75 0.5

0.50 0.0

0.25 0.5

0.00 1.0
0 20 40 60 80 0 20 40 60 80

Figure 5.5: Variables cinemáticas los dos discos después del choque.

c) Haga una gráfica en 3D de las estrella, los auto-vectores de la estrella y las estrella
sobre sus ejes principales. (ver Figura [5.6])

3
2
1
Z 0

1
3
2 2
1
0
3 1 X
3 2 2
1 0 1 2 3 3
Y

Figure 5.6: Estrella naciente: Puntos originales (azul), los ejes principales del cuerpo (rojo) y
puntos sobre los ejes principales (verde.)
5.8. Ejercicios: Álgebra lineal 73

10. Imagine que uno de los espejos del telescopio espacial James Webb tiene una forma
cuadrada como ilustra la Figura [5.7]. Los vértices están ubicados en las siguientes
posiciones con sus respectivas temperaturas:

P1 = (+1, +1); T (P1 ) = 1 K

P2 = (−1, +1); T (P2 ) = 2 K
P3 = (−1, −1); T (P3 ) = 0.5 K
P4 = (+1, −1); T (P4 ) = 0.3 K (5.60)

Lastimosamente en el punto P = (0, 0.5) se daño el sensor de la temperatura, lo que

genera preocupación en los ingenieros de la NASA.

T= 1 K

(-1,1) (1,1)
T= 2 K T= 1 K

T= 2 K
P=(0,0.5) P=(0,0.5)

T = 0.3 K

T = 0.5 K T = 0.3 K
(-1,-1) (1,-1)
T = 0.5 K

Figure 5.7: Esquema del espejo del telescopio espacial. Los vértices del cuadrado registran las
temperaturas en grados kelvin. En la parte derecha, se ilustra una rotación para minimizar la
temperatura en el punto P = (0, 0.5).

La NASA quiere hacer una estimación de la temperatura en el punto donde se daño el

sensor, usando la restricción de la temperatura en los vértices. Adicionalmente, se quiere
rotar el espejo hasta que la temperatura tome el valor más bajo posible y garantizar la vida
útil del espejo. Los estudiantes del curso de Métodos computacionales proponen calcular
una interpolación espacial de la temperatura en el espejo para estimar la temperatura
74 Chapter 5. Álgebra lineal

en el punto P , y luego aplicar un operador de rotación para evaluar como cambia la

temperatura en P a medida que se rota el espejo. Les propongo la siguiente estrategia:

a) Proponga que la temperatura en la región del espejo está dada por una aproximación
bilineal:

1 X
1
T (x, y) ∼
X
= aij xi y j (5.61)
i=0 j=0

Programar esta función que define la interpolación de la temperatura en el espejo. Usar

Sympy para verificar que la función es de la forma T (x, y) = a00 + a10 x + a01 y + a11 xy.
La función debe tener un parámetro p para llenar la matriz de coeficientes, tal que:
a00 = p0 , a10 = p1 , a01 = p2 , y a11 = p3 .

b) Defina el vector de posición de los vértices como: position = [Link]((4,2)) .

Llene los vectores en el orden que aparecen en (5.60)

c) Para calcular la matriz de coeficientes de la interpolación, debe resolver el problema

matricial que aparece de exigir que la temperatura en los vértices tenga los valores
presentados en (5.60). Por ejemplo, se debe cumplir 1 = a00 +a10 (1)+a01 (1)+a11 (1)(1).

d) Verifique que la temperatura en los vértices cumpla con las condiciones de frontera.

e) Grafique el mapa de temperatura en la región donde está el espejo. Deberı́a obtener

algo como: Figura [5.8].

1.00

0.75

0.50 2.00
1.75
0.25 1.50
T(x,y)

1.25
1.00
0.00 0.75
0.50
0.25 0.25
1.0
0.50 0.5
1.0 0.0
0.75 0.5
0.0 0.5 y
x 0.5
1.0 1.0
1.00
1.0 0.5 0.0 0.5 1.0

Figure 5.8: Mapa de temperatura al interior del espejo del telescopio.

f) Estime la temperatura en el punto donde se daño el sensor: T (0, 0.5) ≈ 1.225 K. (Si
llegas hasta este punto tienes 3.0)
5.8. Ejercicios: Álgebra lineal 75

g) Para el proceso de minimización, defina una función que rote el vector de los vértices
del espejo.

! ! !
Pxi (θ) cos(θ) −sin(θ) Pxi
= (5.62)
Pyi (θ) sin(θ) cos(θ) Pyi

donde i = 0, 1, 2, 3 es cada vértice del espejo. Dentro de la función debe realizar la

rotación de los cuatro vértices.

h) Una vez rote los cuatro vértices, realice nuevamente la interpolación usando los mismos
valores de temperatura de los vértices y calcule la temperatura en el punto donde se
daño el sensor.

i) Defina una sucesión de ángulos, por ejemplo:

theta = [Link](0,2*[Link],200),
¿Para cuál ángulo se minimiza la temperatura del espejo en el punto P . (Si encuentras
el ángulo donde la temperatura en P es mı́nima tienes 5.0)

11. Descargue los datos del generador parabólico ubicados en: [Link]
asegura4488/DataBase/blob/main/MetodosComputacionales/[Link].

a) En cualquier fuente confiable de aprendizaje automático puede encontrar el algoritmo

de k-vecinos más cercanos. Implemente una clase en Python que realice la clasificación
de los puntos que caen fuera del cı́rculo (0) y los puntos internos (1). Esta clase
debe tener como mı́nimo un constructor, un método de entrenamiento y un predictor.
Podrı́a tener un método para calcular el score del método; pero es opcional.

b) Mida el score del clasificador durante el entrenamiento (Training sample), que significa
medir cuál es el porcentaje de eventos que se clasifican bien del total (My result is
around 98% using k = 3, which is ”quite good”).

c) Haga una optimización hiper-paramétrica del número de vecinos. Lo que significa

variar el valor k hasta que el score sea lo más grande posible.

d) Usando el código del generador parabólico genere una muestra nueva de 1000 eventos
(Testing sample) y mida el score de la muestra de prueba. Esta sobre-entrenado
(overfitting) su algoritmo? Justifique.
Este video podrı́a ayudar: [Link]
76 Chapter 5. Álgebra lineal

12. Modifique el código sobre redes neuronales para entrenar una red neuronal sencilla. La
red neuronal debe comportarse como una compuerta XNOR. Dicha compuerta tiene la
siguiente tabla de verdad:

A B Output
0 0 1
1 0 0
0 1 0
1 1 1

La red debe tener la topologı́a más simple posible (un par de neuronas en un capa o dos
neuronas en dos capas) y entrenada usando solo números aleatorios.
Chapter 6

Probabilidad

Section 6.1
Definición de probabilidad

Para definir intuitivamente el concepto de probabilidad debemos definir el concepto de evento

aleatorio. Un evento aleatorio es un suceso que puede ocurrir de diversas maneras, por ejemplo,
¿va a llover?, ¿cuánto va a tardar el próximo transporte?, ¿qué cara del lado va a salir?, etc.
Para los antiguos griegos la ocurrencia de eventos se atribuye a la voluntad divina y por tanto
no era posible determinar lo que iba a suceder. Por otro lado, los romanos consideraban que
los eventos ocurren por cuestiones de suerte o azar. Esta naturaleza intrı́nseca de los eventos
aleatorios llevan asociada un incertidumbre en el resultado. La probabilidad busca cuantificar
el grado de plausibilidad de cada evento que puede suceder, en ese sentido, esta disciplina nos
permite conocer aplicando algunas reglas, los eventos que pueden suceder con mayor frecuencia
en un experimento.

Un evento en probabilidad está definido por el concepto de conjunto, el conjunto que contiene
a todos los conjuntos se denomina espacio muestral Ω. Tomemos un conjunto genérico A
del espacio muestral, este conjunto debe pertenecer a una estructura denominada σ−álgebra;
para ser considerado un evento válido. El concepto de σ−álgebra aparece en una rama de la
matemática denominada teorı́a de la medida.

Sea Ω un espacio muestral y definamos a F una colección de sub-conjuntos de Ω. Esta

colección se denomina σ−álgebra, si cumple las siguientes condiciones.

77
78 Chapter 6. Probabilidad

• ∅ ∈ F.

• A ⊂ F, entonces Ac ⊂ F.

• {An }n∈N ∈ F, entonces ∪∞

i=1 Ai ∈ F.

De esta manera, (Ω, F) se denomina un espacio medible. Por ejemplo, tomemos el siguiente
espacio medible:

Ω = {0, 1, 2}
F = {∅, Ω, {0, 1}, {2}} (6.1)

El lector se puede convencer que se cumplen todas las condiciones de espacio medible. Dada
la definición de espacio medible, podemos definir lo que es una medida de probabilidad. Sea P
una aplicación que mapea desde la σ−álgebra a los números reales R.

P : F −→ R. (6.2)

Esta aplicación asigna un número real a un conjunto de la σ−álgebra.

A −→ P(A) (6.3)

En 1933, Kolmogórov formula el sistema de axiomas para que los espacios medibles sean
espacios de probabilidad. Los axiomas son los siguientes:

1. P(Ω) = 1

2. ∀A ∈ F, P(A) ≥ 0
Pn
3. P(∪∞
i Ai ) = i=1 P(Ai ), si Ai ∩ Aj = ∅ ∀i 6= j.

La estructura (Ω, F, P) se denomina espacio de probabilidad. A nivel práctico, estamos

interesados en contestar las siguientes preguntas:

1. ¿Cuál es la probabilidad de que ocurra A1 y A2 ?, P(A1 ∩ A2 ).

2. ¿Cuál es la probabilidad de que ocurra A1 o A2 ?, P(A1 ∪ A2 ).

6.1. Definición de probabilidad 79

3. ¿Cuál es la probabilidad de que ocurra A1 pero no ocurra A2 , P(A1 − A2 ).

4. ¿Cuál es la probabilidad de que no ocurra A, P(Ac ).

Para contestar estas preguntas, se requiere establecer la regla de cálculo en esos espacios
de probabilidad. Esencialmente, la cardinalidad de un conjunto es una medida del número de
elementos de un conjunto y permite calcular la probabilidad de un evento A en un espacio
muestral Ω. En particular, cuando los resultados de un evento aleatorio son equiprobables se
tiene la regla de Laplace. La regla de Laplace establece que la probabillidad del evento A, es
la cardinalidad del evento dividido entre el cardinal del espacio muestral completo:

|A|
P(A) = . (6.4)
|Ω|
La regla de Laplace también se cumplira en la composición de experimentos aleatorios
independientes con resultados equiprobables. Por ejemplo, el lanzamiento de una moneda
equilibrada, seguido del lanzamiento de un dado equilibrado.

Ejemplo 1: En una ciudad el 30% de sus ciudadanos leen sobre ciencia, el 60% leen sobre
arte, y el 80% uno de los dos. ¿Cuál es la probabilidad de que una persona elegida al azar lea
ambos?

Nombremos A = Leer sobre ciencia B = Leer sobre arte. Entonces P(A) = 0.3 y P(B) = 0.6.
Leer uno o el otro es P(A ∪ B) = 0.8. Usamos:

P(A ∩ B) = P(A) + P(B) − P(A ∪ B) = 0.3 + 0.6 − 0.8 = 0.1 (6.5)

6.1.1 Definición frecuentista de probabilidad

En términos prácticos, la probabilidad se mide en términos de la frecuencia de aparición de
un evento especı́fico A. Sea N el número de veces que se repite una experiencia aleatoria.
Denominamos la frecuencia absoluta del evento A, al número de veces que ocurre A y se
denomina como f (A). De la misma forma, definimos la frecuencia relativa como la proporción
de eventos que aparece A, la cual está dada por:

f (A)
fr (A) = (6.6)
N
80 Chapter 6. Probabilidad

Cuando se realiza la experiencia aleatoria un número grande de veces y siguiendo la ley

de los grandes números, la frecuencia relativa del evento A se aproxima asintóticamente a la
probabilidad de ocurrencia.

lim fr (A) = P(A). (6.7)

N →+∞

Esto liga las experiencias aleatorias numéricas como: MonteCarlo, Bootstrapping,

estimación de parámetros, etc; con las predicciones teóricas que provienen de la probabilidad.

Section 6.2
Ejercicios: Axiomas de la probabilidad

1. Sean P1 y P2 dos medidas de probabilidad. Definamos P = a1 P1 + a2 P2 , donde a1 + a2 = 1

y a1 , a2 ∈ R+ . ¿Es P una medida de probabilidad?
Hint: Verifique los axiomas de Kolmogorov para P.

2. Sea Ω = {1, 2}, F = σ(Ω) y P una aplicación definida sobre F dada por:



 0 si A = {∅}


 1/3 si A = {1}
P(A) = (6.8)


 2/3 si A = {2}

1 si A = {1, 2}


1
Muestre detalladamente que P es una medida de probabilidad.

3. Sea (Ω, F, P) un espacio de probabilidad. Demuestre las siguientes propiedades básicas

de esta medida usando los axiomas de Kolmogorov y diagramas de Venn:

a) P(∅) = 0.

b) P(Ac ) = 1 − P(A).

c) Si dos eventos A y B son tales que A ⊂ B, entonces P(B) = P(A) + P(B − A).

d) Dado un evento A, P(A) ≤ 1, ∀A ∈ F. Use el teorema anterior.

1
Tomado de Probabilidad, Liliana Blanco Castañeda UNAL.
6.2. Ejercicios: Axiomas de la probabilidad 81

e) A ⊆ B, entonces P(A) ≤ P(B).

f) P(A ∪ B) = P(A) + P(B) − P(A ∩ B).

g) P(A ∪ B ∪ C) = P(A) + P(B) + P(C) − P(A ∩ B) − P(A ∩ C) − P(B ∩ C) + P(A ∩ B ∩ C).

h) Probabilidad de la diferencia: Use A − B = A ∩ B c , para mostrar: P(A − B) =

P(A) − P(A ∩ B).

i) Probabilidad de la diferencia simétrica: Use (A − B) ∪ (B − A), para mostrar: P((A ∩

B c ) ∪ (B ∩ Ac )) = P(A) + P(B) − 2P(A ∩ B). Esto significa la probabilidad que ocurra
A o B pero no ambos eventos al tiempo.

6.2.1 Independencia de eventos

Uno de los conceptos más importantes en probabilidad es la independencia de eventos. Se
establece que dos eventos A y B son independientes, si la probabilidad de la intersección está
dada por:

P(A ∩ B) := P(A) · P(B). (6.9)

Esto significa que el hecho de que ocurra A no influye en B y recı́procamente.

6.2.2 Probabilidad condicional

En otras ocasiones el resultado del evento B puede influir en el resultado del evento A. En este
caso, se define la probabilidad de A condicionado al evento B como:

P(A ∩ B)
P(A/B) := (6.10)
P(B)
Esto puede interpretarse siguiendo la noción de cardinalidad. Dado que ya ocurrió B, el
espacio completo de redujo al conjunto B y el numerador indica el cardinalidad del conjunto
donde ocurre A y B a la vez. Note que es posible, dada la probabilidad condicional, calcular
la probabilidad de la intersección.

P(A ∩ B) = P(A/B)P(B). (6.11)

Esta formula es válida para eventos no independientes.

82 Chapter 6. Probabilidad

6.2.3 Teorema de la probabilidad total

Supongamos que tenemos una partición del espacio muestral Ω tal que:

A1 , A2 , ..., An ⊆ Ω, (6.12)

que cumple las siguientes condiciones:

a) Ω = A1 ∪ A2 ∪ ... ∪ An

b) Ai ∩ Aj = ∅, i 6= j

Si tenemos un suceso B de una partición de un espacio muestral. El teorema de la

probabilidad total establece que la probabilidad P(B) es:

n
X
P(B) = P(B/Ai )P(Ai ). (6.13)
i=1

donde P(B/Ai ) es la probabilidad del suceso B restringido a cada uno de los elementos de
la partición y P (Ai ) es la probabilidad de que ocurra ese elemento de la partición.

Ejemplo 1: En una sala de aeropuerto, 15 hombre y 20 mujeres. De los 15 hombres, 5

llevan gafas y de las 20 mujeres, 8 llevan gafas.

a) Si se elige una persona al azar, ¿cuál es probabilidad de que tenga gafas?

Identificamos las cantidades P(H) = 15/35, P(M ) = 20/35. Note que la probabilidad
condicional de que lleve gafas dado que es hombre es: P(G/H) = 5/15 y la probabilidad
condicional de que lleve gafas dado que es mujer es: P(G/M ) = 8/20. La probabilidad total
está dada por:

P(G) = P(G/H)P(H) + P(G/M )P(M )

= (5/15)(15/35) + (8/20)(20/35) = 13/35 ≈ 0.37 (6.14)

b) Si se elige una persona al azar, ¿cuál es probabilidad de que no tenga gafas?

6.2. Ejercicios: Axiomas de la probabilidad 83

Identificamos las cantidades P(H) = 15/35, P(M ) = 20/35. Note que la probabilidad
condicional de que no lleve gafas dado que es hombre es: P(Ḡ/H) = 10/15 y la probabilidad
condicional de que no lleve gafas dado que es mujer es: P(Ḡ/M ) = 12/20. La probabilidad
total está dada por:

P(Ḡ) = P(Ḡ/H)P(H) + P(Ḡ/M )P(M )

= (10/15)(15/35) + (12/20)(20/35) = 22/35 ≈ 0.63 (6.15)

Ejemplo 2: Se lanzan tres monedas al aire. Calcular la probabilidad de:

1. Obtener 3 caras.
Identificamos las probabilidades de obtener cara o cruz, P(c) = P(x) = 1/2
respectivamente. En el caso, P(ccc) = P(c) · P(c) · P(c) = (1/2)(1/2)(1/2) = 1/8, dado
que el lanzamiento de cada moneda es independiente.

2. Obtener 2 caras y una cruz.

Identificamos las probabilidades de obtener cara o cruz, P(c) = P(x) = 1/2
respectivamente, y la probabilidad de obtener dos caras y una cruz con P(A). En este
caso, hay tres ramas diferentes del diagrama de árbol que cumplen con la condición:
P(ccx), P(cxc) y P(xcc). La probabilidad total será:

P(A) = P(ccx) + P(cxc) + P(xcc)

= P(c) · P(c) · P(x) + P(c) · P(x) · P(c) + P(x) · P(c) · P(c)
= 1/8 + 1/8 + 1/8 = 3/8 (6.16)

Este problema también se puede solucionar usando combinaciones sin repetición dado que
no importa el orden.

  2  2
3 1 1 3! 1 1
P(A) = = = 3/8. (6.17)
2 2 2 (3 − 2)!2! 2 2
Formula que resulta más útil en los problemas donde el conteo no es sencillo de calcular
usando el diagrama de árbol.
84 Chapter 6. Probabilidad

Ejemplo 3: Se tienen 2 urnas A y B, en la urna A se tienen 3 bolas negras y 4 blancas y

en la urna B se tienen 5 bolas negras y 3 blancas. Se saca una bola de la urna A, y sin verla,
se introduce en la urna B. Si se saca una segunda bola de la urna B, calcule la probabilidad
de los siguientes sucesos:

a) La segunda bola sea negra: P(N2 ).

b) La segunda bola sea blanca: P(B2 ).

Calculemos las probabilidades de la primera bola: la probabilidad de que la primera bola sea
negra es: P(N1 ) = 3/7 y la probabilidad de que la primera sea blanca es: P(B1 ) = 4/7. Notar
que si la primera bola es negra, en la urna B quedan 6 bolas negras y 3 blancas; de este modo
la probabilidad condicional de obtener una bola negra y una bola blanca es P(N2 /N1 ) = 6/9 y
P(B2 /N1 ) = 3/9. Si por el contrario, la primera bola es blanca, en la urna B quedan 5 bolas
negras y 4 blancas; de este modo la probabilidad condicional de obtener una bola negra y una
bola blanca es P(N2 /B1 ) = 5/9 y P(B2 /B1 ) = 4/9. Lo que conduce a la regla de probabilidad
total:

P(N2 ) = P(N2 /N1 )P(N1 ) + P(N2 /B1 )P(B1 )

= (6/9)(3/7) + (5/9)(4/7) = 38/63 (6.18)

La probabilidad de que la segunda bola sea blanca es:

P(B2 ) = P(B2 /N1 )P(N1 ) + P(B2 /B1 )P(B1 )

= (3/9)(3/7) + (4/9)(4/7) = 25/63 (6.19)

Section 6.3
Ejercicios: Probabilidad condicional y total

1. En una conferencia hay 1000 participantes distribuidos del siguiente modo: 185 hombres
usan gafas, 415 hombres no usan gafas, 115 mujeres que usan gafas. Determine la
probabilidad de:

a) Ser hombre P(H) = 3/5.

6.4. Teorema de Bayes 85

b) Ser Mujer P(M ) = 2/5.

c) Usar gafas P(G) = 3/10.

P(G∩M )
d) Que lleve gafas si sabemos que es mujer P(G/M ) = P(M )
= 23/80.

Hombres Mujeres Totales

Usa gafas 185 115 300
No usa gafas 415 285 700
Totales 600 400 1000

Table 6.1: Tabla de contingencia de genero y el uso de gafas.

2. Lanzamos un dado de 6 caras. Si sale 1 o 2 extraemos una bola de la urna 1. Si sale

3,4,5,6 extraemos una bola de la urna 2. En la urna 1 tenemos 3 bolas rojas, 1 negra y 6
verdes, y en la urna 2 tenemos 6 bolas rojas, 2 negras y 2 verdes. Calcular la probabilidad
de que la bola obtenida:

a) Sea roja: P(R) = 1/2.

b) sea negra: P(N ) = 1/6.

c) Sea de la urna 1 si se ha obtenido una bola negra: P(1/N ) = 1/5.

d) Sea de la urna 2 si se ha obtenido una bola negra: P(2/N ) = 4/5.

3. Tengo una bolsa con dos dulces de limón y tres de fresa. Si me como un dulce y después
otro, ¿cuál es la probabilidad de que ambos hayan sido de fresa. P(F ∩ F ) = 3/10

Section 6.4
Teorema de Bayes

El teorema de Bayes responde la pregunta inversa de la probabilidad total. Dado que ya ocurrió
el evento B, este teorema establece la regla para calcular la probabilidad de que ocurra el i-ésimo
elemento de la partición.
86 Chapter 6. Probabilidad

P(B/Ai )P(Ai )
P(Ai /B) =
P(B)
P(B/Ai )P(Ai )
= Pn (6.20)
i=1 P(B/Ai )P(Ai )

Ejemplo 1: En el curso de Métodos computacionales, el 20% de la población estudia

Matemáticas, 30% estudia Geociencias y el resto estudia Fı́sica. El 30% de los estudiantes de
matemáticas son mujeres, el 50% en Geociencias y el 45% en Fı́sica.

a) Si un estudiante es elegido al azar y es mujer, ¿qué probabilidad habrá de que sea estudiante
de matemáticas?

Identificamos las cantidades: P(M ate) = 0.2, P(Geo) = 0.3, P(F isi) = 0.5. Adicionalmente,
tenemos las probabilidades condicionales de ser mujer. P(M/M ate) = 0.3, P(M/Geo) = 0.5,
P(M/F isi) = 0.45.

Lo que se pide calcular es:

P(M/M ate)P(M ate)

P(M ate/M ) = (6.21)
P(M )

Usemos el teorema de probabilidad total para calcular que sea mujer:

P(M ) = P(M/M ate)P(M ate) + P(M/Geo)P(Geo) + P(M/F isi)P(F isi)

= (0.3)(0.2) + (0.5)(0.3) + (0.45)(0.5) = 0.435 (6.22)

Entonces:
P(M/M ate)P(M ate) (0.3)(0.2)
P(M ate/M ) = = = 4/29 ≈ 0.14 (6.23)
P(M ) 0.435

b) Si un estudiante es elegido al azar y es hombre, ¿qué probabilidad habrá de que sea estudiante
de Matemáticas?
6.5. Ejercicios: Teorema de Bayes 87

Identificamos las cantidades: P(M ate) = 0.2, P(Geo) = 0.3, P(F isi) = 0.5. Adicionalmente,
tenemos las probabilidades condicionales de ser hombre. P(H/M ate) = 0.7, P(H/Geo) =
0.5, P(H/F isi) = 0.55.

Lo que se pide calcular es:

P(H/M ate)P(M ate)

P(M ate/H) = (6.24)
P(H)
Usemos el teorema de probabilidad total para calcular que sea hombre:

P(H) = P(H/M ate)P(M ate) + P(H/Geo)P(Geo) + P(H/F isi)P(F isi)

= (0.7)(0.2) + (0.5)(0.3) + (0.55)(0.5) = 0.565 (6.25)

Entonces:
P(H/M ate)P(M ate) (0.7)(0.2)
P(M ate/H) = = = 28/113 ≈ 0.25 (6.26)
P(H) 0.565

Ejemplo 2: En el ejemplo 3 de la probabilidad total. Calcule la probabilidad de que la

bola extraı́da en la urna A sea negra, dado que el resultado de la urna B es blanca P(N1 /B2 ).
Usando los valores calculados anteriormente:

P(B2 /N1 )P(N1 )

P(N1 /B2 ) =
P(B2 /N1 )P(N1 ) + P(B2 /B1 )P(B1 )
(4/9)(4/7)
= = 16/25 (6.27)
(25/63)

Section 6.5
Ejercicios: Teorema de Bayes

1. Un médico ha observado que el 40% de sus pacientes fuma. De este 40% el 75% son
hombres. Entre los pacientes que no fuman, el 60% son mujeres. Si se selecciona un
paciente al azar:

a) ¿Cuál es la probabilidad de que sea mujer? P(M ) = 23/50.

88 Chapter 6. Probabilidad

b) ¿Cuál es la probabilidad de que sea hombre y fume? P(F ∩ H) = 3/10.

c) Si el paciente es mujer, ¿cuál es la probabilidad de que fume? P(F/M ) = 5/23.

2. Suponga que se quiere estimar el número medio de estudiantes que se aproximan a la

secretaria de fı́sica siguiendo una distribución de Poisson. El problema podrı́a tener
información inicial que puede ser utilizada como una función a priori de probabilidad. Esto
asigna probabilidades a los posibles valores del parámetro λ. Supongamos que el número
medio a priori de estudiantes es: 1,2,3,4 con las probabilidades Π(λ) = [0.4, 0.3, 0.2, 0.1]
según las secretarias 1,2,3 y 4 (información subjetiva). El número de estudiantes
observados durante cierto intervalo de tiempo es x = 4.

a) Haga la estimación inicial del parámetro λ̂.

b) Calcule la función de verosimilitud de la observación para cada modelo P(x/λi ) =
L(x/λi ) donde i = 1, 2, 3, 4.
c) Aplicando el teorema de Bayes, calcule la función de distribución a posteriori para
cada modelo:

L(x/λi )Π(λi )
P(λi /x) = P i, j = 1, 2, 3, 4 (6.28)
∀j L(x/λj )Π(λj )

d) Verifique que la distribución posterior está normalizada.

e) Indique el modelo que más probablemente explique los datos.
f) Usando la información subjetiva de cada secretaria y la verosimilitud de la observación
(i,e. la distribución a posteriori). Calcule el mejor parámetro λ̂.

Section 6.6
Técnicas de conteo

Dado que es necesario contar el número de elementos asociados a conjuntos finitos, necesitamos
los axiomas fundamentales de conteo y algunas herramientas de combinatoria. Los axiomas
fundamentales de conteo son los siguientes:

1. Principio multiplicativo: Si un evento A puede ocurrir de n maneras diferentes y otro

evento B puede ocurrir de m diferentes, entonces ambos eventos pueden ocurrir de n · m
maneras distintas.
6.6. Técnicas de conteo 89

2. Principio sumativo: Si un evento A puede ocurrir de n maneras diferentes o m maneras

diferentes, siendo incompatibles las unas con las otras, entonces existen n + m maneras
de ocurrir A.

Adicionalmente, las posibles configuraciones depende si importa o no el orden en los

sub-conjuntos. Si importa el orden se calculan las variaciones y permutaciones, y si no importa
el orden, calculamos la combinaciones. Además, podrı́an o no repetirse elementos en cada
sub-conjunto.

Variaciones sin repetición:

n!
Vrn = (6.29)
(n − r)!
Variaciones con repeticion:

Vrn = nr (6.30)

Permutaciones sin repetición:

P n = n! (6.31)

Permutaciones con repetición:

n!
Pnn1 ,n2 ,...,n2 = (6.32)
n1 !n2 !...nn !
Por otro lado, cuando no importa el orden de las configuraciones:

Combinaciones sin repetición:

 
n n!
Crn = = (6.33)
r r!(n − r)!
Combinaciones con repetición:

 
n+r−1 (n + r − 1)!
Crn = = (6.34)
r r!(n − 1)!
90 Chapter 6. Probabilidad

Section 6.7
Ejercicios: Técnicas de conteo

Realizar los siguientes cálculos estableciendo si es variación, permutación o combinación; con o

sin repetición.

1. Carlos, Manuel, Sandra correrán los 100 metros planos. ¿De cuántas formas puede quedar
el podio de primer y segundo lugar? Solo competirán ellos tres. R:6

2. ¿De cuántas formas se puede preparar una ensalada de frutas con solo 2 ingredientes, si
se cuenta con plátano, manzana y uva? R:3

3. ¿De cuantas formas pueden hacer cola 5 amigos para entrar al cine? R:120

4. ¿De cuántas formas puede un juez otorgar el primero, segundo y tercer premio en un
concurso que tiene 8 participantes R:336

5. El capitán de un barco solicita 2 marineros para realizar un trabajo, sin embargo, se

presentan 10. ¿De cuántas formas podrá seleccionar a los 2 marineros? R:45

6. Eduardo tiene 7 Libros, ¿De cuántas maneras podrá acomodar cinco de ellos de un
estante? R:2520

7. En un salón de 10 alumnos, ¿de cuántas maneras se puede formar un comité formado por
2 de ellos? R:45

8. ¿Cuántas palabras diferentes se puede formar con las letras de la palabra REMEMBER?
R:1680

9. Un club de basketball tiene 12 jugadoras, una de ellas es la capitana Marı́a. ¿Cuántos

equipos diferentes de 6 jugadoras se pueden formar, sabiendo que en todos ellos siempre
debe estar la capitana Marı́a. R:462

10. Con 4 frutas diferentes, ¿cuántos jugos surtidos se pueden preparar?. Un jugo surtido se
debe preparar con al menos 2 frutas. R:11

11. En un curso de 10 estudiantes se desea escoger presidente, vicepresidente y secretario.

¿De cuántas formas de pueden seleccionar los 3 estudiantes? R:720
6.8. Ejercicios: Generales de probabilidad 91

12. En un campeonato compiten 8 equipos ¿de cuántas maneras diferentes se podrı́an ganar
los premios de campeón y sub-campeón? R:56

13. ¿Cuántos números de 3 cifras distintas se pueden formar con los dı́gitos del 1 al 7? R:210

14. ¿Cuántos números de 3 cifras se pueden formar con los dı́gitos del 1 al 7? R:343

15. De un grupo de 10 estudiantes se quiere seleccionar un comité al azar de 3 estudiantes.

¿De cuántas maneras diferentes se puede seleccionar el comité? R:120

16. ¿Cuántas placas diferentes se pueden hacer con 3 letras y 3 dı́gitos? R:17576000

17. n Personar van a jugar cartas alrededor de una mesa, ¿de cuántas maneras diferentes se
pueden sentar? R:(n − 1)!

18. En una heladerı́a ofrecen 7 diferentes sabores, ¿cuántas combinaciones de helado de 3

sabores se pueden hacer? R:84

19. En un almacén venden 6 diferentes sabores de gaseosas, ¿de cuántas formas se pueden
seleccionar 3 gaseosas?. ¿De cuántas formas diferentes se pueden seleccionar 3 gaseosas?
R:56, R:20

20. Demostrar la formula de combinaciones con repetición.

 
n n+r−1
Cr = (6.35)
r

21. Calcular la probabilidad de ganar el Baloto comprando un solo billete. R:7.1511 × 10−8

22. Cuántas sumas de 3 enteros no negativos dan 10. R:66

23. Se tienen 9 llaves: 3 rojas, 3 azules y 3 verdes. Si elegimos 4, ¿de cuántas formas se
pueden distribuir los colores? R:12

Section 6.8
Ejercicios: Generales de probabilidad

Algunos tomados del libro de probabilidad y estadı́stica de Walpole [5].

1. Se lanza un dado equi-probable dos veces. Dados los siguientes eventos:

92 Chapter 6. Probabilidad

A) La suma de los resultados es menor o igual a 3.

B) El resultado del primer lanzamiento es impar.

Calcular P(A) = 1/12, P(B) = 1/2, P(A ∪ B) = 19/36 y P(Ac ) = 11/12

2. Se prueban 5 celulares de un lote de 50 equipos donde existen 2 defectuosos. Si el

muestreo se realiza sin reposición. Hallar la probabilidad que al menos un celular sea
defectuoso. P(A) = 47/245. Piense en Ac y en el principio multiplicativo o en una
estrategia combinatoria.

3. En una cierta ciudad el 60% de los propietarios están suscritos al diario y el 80% al cable.
Adicionalmente, el 50% están suscritos a ambos. Si un propietario es elegido al azar:

(a) ¿Cuál es la probabilidad que esté suscrito a uno de los dos servicios? P(A∪B) = 0.9.
(b) ¿Cuál es la probabilidad que esté suscrito al diario o al cable, pero no a ambos
servicios? P((A ∩ B c ) ∪ (B ∩ Ac )) = 0.4.

4. Calcular la probabilidad que n personas (n ≤ 365) tengan fechas diferentes de cumpleaños,

i.e, escribir la formula general de cálculo. Grafique la probabilidad P(n ≤ 80) como
función de n. Los números son demasiado grandes, pero Python puede manejar dichas
cantidades.
[Link]

5. Se lanzan dos dados equi-probables y se observan los siguientes eventos:

a) La suma es 8: P(A) = 5/36.

b) El segundo dado es impar: P(B) = 1/2.

Calcule P(A ∩ B) y P(A) · P(B). ¿Qué podrı́a concluir?

6. Se lanza simultáneamente 3 dados de 6 caras. ¿Cuál es la probabilidad de obtener 1 par?

P(A) = 5/12. Realice el cálculo de esta probabilidad usando un experimento virtual con
N = 105 eventos.

7. Se lanza simultáneamente 5 dados de 6 caras. ¿Cuál es la probabilidad de obtener?

a) 1 par: P(A) = 25/54.

6.8. Ejercicios: Generales de probabilidad 93

1.0

0.8

Pn 0.6

0.4

0.2

0.0
0 10 20 30 40 50 60 70 80
Personas
Figure 6.1: Probabilidad P(n ≤ 80) de que tengan una edad diferente.

b) 2 pares distintos: P(B) = 25/108.

c) 4 de la misma cara: P(C) = 25/1296.

Realice el cálculo de la probabilidad de los numerales a) y b) usando un experimento

virtual con N = 105 eventos.

8. Se lanzan simultáneamente 4 monedas. Determine la probabilidad de obtener dos caras

y dos sellos P(A) = 3/8. Realice el cálculo de esta probabilidad usando un experimento
virtual con N = 105 eventos, y etiquetando los resultados con +1 y -1 para cara y sello
respectivamente.

9. En el ejercicio anterior, imagine que las monedas están truncadas de tal manera que la
probabilidad de que la moneda 1 sea cara es p1 y que sea sello es 1−p1 . Usando el árbol de
probabilidad, ¿cuál es la expresión de la probabilidad de obtener dos caras y dos sellos de
este evento? Si el truncamiento de las monedas 1 y 2 puede variar como: 0.1 < p1 < 0.9 y
0.1 < p2 < 0.5, use el árbol de probabilidad para graficar la superficie de probabilidad del
evento A. ¿En qué punto la probabilidad es mı́nima y máxima, y cuáles son sus valores?

10. En un juego de Poker que consta de 5 cartas, encuentre la probabilidad de tener:

94 Chapter 6. Probabilidad

0.44
0.42
0.40
0.38
0.36
0.34
0.32
0.30
0.10
0.15
0.20
0.25
0.30
0.35
p1
0.40
0.45
0.50
0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1
p2

Figure 6.2: Superficie de probabilidad del evento A.

a) 3 ases: P(A) = 94/54145

b) 4 cartas de corazones y 1 de bastos: P(B) = 143/39984

11. Si se seleccionan al azar 3 libros de un estante que contiene 5 novelas, 3 libros de poemas
y un diccionario, ¿cuál es la probabilidad de que:

a) se tome el diccionario: P(A) = 1/3

b) se escojan 2 novelas y un libro de poemas: P(B) = 5/14

Section 6.9
Distribuciones discretas de probabilidad

Necesitamos una expresión que nos permita representar todas las probabilidades de una variable
aleatoria X, esta función se representa por f (x) y tiene la siguiente connotación:

f (x) = P(X = x) (6.36)

el conjunto de pares ordenados (x, f (x)) se denomina función de densidad de probabilidad.

Definición 1. El conjunto de pares ordenados (x, f (x)) se denomina función de densidad de

probabilidad de la variable aleatoria X si, para cada resultado posible x:
6.9. Distribuciones discretas de probabilidad 95

f (x) ≥ 0
X
f (x) = 1
x
P(X = x) = f (x) (6.37)

x 0 1 3
P(X = x) 1/3 1/2 1/6

Table 6.2: Distribución de probabilidad

Definición 2. Sea X una variable aleatoria con distribución de probabilidad f (x). La media o
valor esperado de X es:

X
µ̂ = E(X) = xf (x) (6.38)
x

Del ejemplo anterior el valor esperado será:

2
X
µ̂ = xi f (xi ) = x0 f (x0 ) + x1 f (x1 ) + x2 f (x2 ) = 0(1/3) + 1(1/2) + 3(1/6) = 1 (6.39)
i=0

Se tienen algunos resultados útiles para la esperanza de una variable aleatoria.

E(c) = c
E(cX) = cE(X)
E(c + X) = c + E(X)
E(X + Y ) = E(X) + E(Y ) (6.40)

donde c es constante.

Definición 3. Sea X una variable aleatoria con distribución de probabilidad f (x). La varianza
de X es:

X
σˆ2 = V ar(X) = (x − E(x))2 f (x). (6.41)
x
96 Chapter 6. Probabilidad

Esto es una medida de la dispersión en los datos. Se tienen algunos resultados útiles para
la varianza de una variable aleatoria.

V ar(c) = 0
V ar(cX) = c2 V ar(X)
V ar(X + c) = V ar(X)
V ar(X + Y ) = V ar(X) + V ar(Y ), X, Y independientes (6.42)

donde c es constante. En el caso donde las variables X y Y no sean independientes, tenemos

la definición de covarianza.

6.9.1 Covarianza

Definición 4. Sea X e Y dos variables aleatoria con distribución de probabilidad conjunta

f (x, y). La covarianza de X y Y está dada por:

X
σ̂xy = Cov(X, Y ) = (x − E(x))(y − E(y))f (x, y). (6.43)
x,y

Proof: Definamos µX = E(X) y µY = E(Y ) y usando la propiedad de la esperanza

E(X + Y ) = E(X) + E(Y ) = µX + µY . Entonces:

V ar(X + Y ) = E((X + Y − (µX + µY ))2 )

= E(((X − µX ) + (Y − µY ))2 )
= E((X − µX )2 ) + E((Y − µY )2 ) + 2E((X − µX )(Y − µY )) (6.44)

En la última lı́nea se uso la linealidad de la esperanza, donde el último término representa

la covarianza entre variables. Por tanto:

V ar(X + Y ) = V ar(X) + V ar(Y ) + 2Cov(X, Y ). (6.45)

La covarianza satisface algunas propiedades:

6.9. Distribuciones discretas de probabilidad 97

Cov(aX, bY ) = abCov(X, Y )
Cov(X, X) = V ar(X)
Cov(X, Y ) = E(XY ) − E(X)E(Y ) (6.46)

6.9.2 Correlación
La covarianza nos indicar el nivel de asociación linea que existe entre las variables, sin embargo,
los valores dependen de las unidades en que se midan X y Y . En ese sentido, podemos definir
un estimador independiente de unidades denominado coeficiente de correlación. El coeficiente
de correlación de Pearson está dado por:

Cov(X, Y )
ρ= . (6.47)
σX σY

6.9.3 Distribución de Bernoulli

Sea X una variable aleatoria discreta X ∼ Bernoulli(p), descrita por un único experimento
con dos posibles resultados denominados éxito y fracaso, y con parámetro 0 < p < 1, donde p
es la probabilidad de éxito. Esta función de probabilidad está dada por:

P (X = x) = px (1 − p)1−x , x = 0, 1 (6.48)
Se puede demostrar que el valor medio está dado por E(X) = p.

6.9.4 Distribución Binomial

Sea X una variable aleatoria discreta X ∼ Bin(n, p) con parámetros n ∈ N y 0 < p < 1. Esta
función de probabilidad está dada por:
 
n x
P (X = x) = p (1 − p)1−x , x = 0, 1, ..., n (6.49)
x
Esta distribución está conectada con a distribución de Bernoulli. Si X1 , X2 , ..., Xn son n
variables aleatorias con distribución de Bernoulli Xi ∼ Bernoulli(p), entonces:

n
X
Xi ∼ Bin(n, p). (6.50)
i
98 Chapter 6. Probabilidad

Se puede demostrar que el valor medio está dado por E(X) = np.

6.9.5 Distribución de Poisson

Sea X una variable aleatoria discreta X ∼ P ois(λ) con parámetro λ ∈ R+ . Esta función de
probabilidad está dada por:

e−λ λk
P (X = k) = , k = 0, 1, ... (6.51)
k!
donde k es el número de ocurrencias del suceso.
Se puede demostrar que el valor medio está dado por E(X) = λ.

Section 6.10
Ejercicios: Distribuciones discretas de probabilidad

1. ¿Cuál es valor esperado del experimento de lanzar un dado? Ans: µ̂ = 3.5

2. ¿Cuál es valor esperado del experimento de lanzar dos dados y sumar sus resultados?
Ans: µ̂ = 7.0

3. Un embarque de 10 microchips similares que se envı́a a distribución tiene 3 aparatos

defectuosos. Si una empresa realiza un compra aleatoria de 2 de estos microchips.

a) Muestre que la distribución de probabilidad del número de microchips defectuosos.

7 3



2−x x
f (x) = 10

(6.52)
2

x 0 1 2
P(X = x) 7/15 7/15 1/15

Table 6.3: Distribución de probabilidad del número de microchips defectuosos.

b) ¿Cuál es valor esperado de microchips defectuosos? Ans: µ̂ = 3/5 = 0.6

4. Una caja cuántica tiene 3 electrones, 2 protones y 3 neutrones. Se selecciona una muestra
aleatoria de 4 partı́culas. Si x es el número de electrones e y es el número de protones.
6.11. Distribuciones continuas de probabilidad 99

a) Muestre que la distribución de probabilidad conjunta f (x, y) es:

3

2
3


x y 4−x−y
f (x, y) = 8

(6.53)
4

b) Hallar las distribuciones marginales g(x) y h(y).

c) Halle el valor esperado de electrones: E(x) = 105/70.
d) Halle el valor esperado de protones: E(y) = 1.
e) Calcular la covarianza usando: σxy = E(xy) − E(x)E(y) = −3/14
f) Calcular la covarianza usando: σxy = E((x − µ̂x )(y − µ̂y )) = −3/14
g) Son las variables x e y independientes?

Section 6.11
Distribuciones continuas de probabilidad

Definición 5. La función f (x) es una función de densidad de probabilidad para la variable

aleatoria continua X, definida en el conjunto de los números reales, si:

f (x) ≥ 0 para toda x ∈ R

Z ∞
f (x)dx = 1
−∞
Z b
P (a < X < b) = f (x)dx (6.54)
a

Definición 6. Sea X una variable aleatoria con distribución de probabilidad f (x). La media o
valor esperado de X es:
Z ∞
µ̂ = E(X) = xf (x)dx (6.55)
−∞

6.11.1 Distribución normal

La función de densidad de probabilidad normal está dada por:

1 1 x−µ 2
f (x, µ, σ) = √ e− 2 ( σ
)
, (6.56)
2πσ 2
100 Chapter 6. Probabilidad

donde µ ∈ R y σ 2 ∈ R+ son los parámetros que caracterizan las distribución. Se puede

demostrar que el valor medio está dado por E(X) = µ.

6.11.2 Distribución χ2
Si tenemos Z1 , Z2 , ..., Zk variables aleatorias con Zi ∼ N (0, 1) para i = 1, 2, ..., k, entonces la
variables aleatoria:

k
X
Zi2 ∼ χ2 (k), (6.57)
i=1

donde k > 0 es el número de grados de libertad de la distribución. La función de densidad

de probabilidad está dada por:

( 21 )k k −1 −x/2
χ2 (x, k) = x2 e . (6.58)
Γ( k2 )
Con x > 0 y Γ es la función gamma.

6.11.3 Distribución t-Student

Si tenemos una variable aleatoria X1 ∼ N (0, 1) y otra variable aleatoria X2 ∼ χ2 (n). Se define
la variable aleatoria:

X1
Xn = q (6.59)
X2
n

tiene una función de densidad t-Student dada por:

Γ((ν + 1)/2)
f (x, ν) = √ (1 + x2 /ν)−(ν+1)/2 , (6.60)
πνΓ(ν/2)
donde ν ∈ R+ es el número de grados de libertad y Γ es la función gamma.

Section 6.12
Ejercicios: Distribuciones continuas de probabilidad

1. Dada la función de probabilidad conjunta:

6.12. Ejercicios: Distribuciones continuas de probabilidad 101

(
2
3
(x + 2y) 0 ≤ x ≤ 1, 0 ≤ y ≤ 1.
f (x, y) = (6.61)
0 otro caso

encuentre analı́ticamente y a través del paquete SymPy los siguientes valores:

a) Verifique que sea una función de densidad conjunta válida.

b) Hallar las distribuciones marginales g(x) y h(y).

10
c) E(x) = 18

11
d) E(y) = 18

e) Calcular la covarianza usando: σxy = E(xy) − E(x)E(y) = −0.00617

f) Calcular la covarianza usando: σxy = E((x − µ̂x )(y − µ̂y )) = −0.00617

g) Son las variables x e y independientes?

2. Suponga que el error en la temperatura en un experimento controlado de laboratorio

es una variable aleatoria continua X, que tiene la siguiente función de densidad de
probabilidad:

(
x2
3
si −1 ≤ x ≤ 2
f (x) = (6.62)
0 en otro caso

Encuentre las siguientes probabilidades:

(a) P(0 < X ≤ 1) = 1/9

(b) P(1 < X ≤ 2) = 7/9

3. Una variable aleatoria continua X que puede asumir valores entre x = 2 y x = 3 tiene
una función de densidad f (x) = 1/2.

a) Demuestre que el área bajo la curva es igual a 1.

b) Encuentre P(2 < X < 2.5) = 1/4

c) Encuentre P(X ≤ 2.6) = 0.3

102 Chapter 6. Probabilidad

4. Una variable aleatoria continua X tiene una función densidad:

(
e−x si x>0
f (x) = (6.63)
0 en otro caso

Encuentre el valor esperado de g(X) = e2X/3 = 3. Recuerde que el valor esperado de la

variable aleatoria g(X) está dada por:

Z ∞
µ̂g(X) = g(x)f (x)dx (6.64)
−∞

Section 6.13
Procesos de Markov

Definición 7. Un proceso de Markov es un proceso estocástico discreto {xn : n = 0, 1, ...}, con

espacio de estados discreto S ⊆ {0, 1, ...} que satisfacen la propiedad de Markov.

P(Xn+1 = xn+1 /X0 = x0 , ..., Xn = xn ) = P(Xn+1 = xn+1 /Xn = xn ) (6.65)

Esto significa que la probabilidad de que la variable aleatoria tome el valor xn+1 , dado que
el proceso ha pasado por los pasos anteriores x0 , x1 , ..., xn ; solo depende del estado anterior xn .
Equivalentemente, se puede escribir la probabilidad de que la variable aleatoria tome el valor
xn+1 como:

P(X0 = x0 , X1 = x1 , ..., Xn+1 = xn+1 ) = P(X0 = x0 )P(X1 = x1 /X0 = x0 )P(Xn+1 = xn+1 /Xn = xn )

Note que el proceso completo, requiere un conjunto de probabilidades para el estado inicial.

Definición 8. La distribución de la variable aleatoria X0 se denomina distribución inicial, es

decir, {P(X0 = 0), P(X0 = 1), ...}. Se suponen que la distribución inicial es conocida.

Por otro lado, el conjunto de probabilidades condicionales que describen al proceso se

denominan probabilidades de transición.
6.13. Procesos de Markov 103

Definición 9. La probabilidad P(Xn+1 = i/Xn = j) = Pij (n, n + 1) se denomina probabilidad

de transición del estado j en el tiempo n al estado i en el tiempo n + 1.

Se supone que los valores de la probabilidades de transición son independientes del valor
especı́fico de n, en ese caso, la cadena es estacionaria.

Pij (0, 1) = Pij (1, 2) = ... = Pij (n, n + 1). (6.66)

Ejemplo 1:
Imaginemos que una partı́cula puede estar en cuatro diferentes estados 1,2,3 y 4. Esta
partı́cula puede mantener o cambiar su estado con ciertas probabilidades. Medimos el estado
de la partı́cula en intervalos de tiempo iguales t0 , t1 , ..., tn . Se tienen las siguientes reglas de
transición:

P(1/1) = 0.30, P(1/2) = 0.25, P(1/3) = 0.15, P(1/4) = 0.10

P(2/1) = 0.30, P(2/2) = 0.25, P(2/3) = 0.55, P(2/4) = 0.20
P(3/1) = 0.30, P(3/2) = 0.50, P(3/3) = 0.10, P(3/4) = 0.30
P(4/1) = 0.10, P(4/2) = 0.00, P(4/3) = 0.20, P(4/4) = 0.40 (6.67)

Por ejemplo: P(4/1) es la probabilidad de que pase al estado 4, dado que está en el estado
1.

a) Si en t0 el estado de la partı́cula es 1, ¿cuál es la probabilidad de que esté en el estado 4 en

t2 ?
Para pasar del estado 1 al estado 4 en dos pasos, la regla de probabilidad total define
cuatro caminos independientes: E1 → E1 → E4 , E1 → E2 → E4 , E1 → E3 → E4 , y
E1 → E4 → E4 . Entonces, la probabilidad es:

P(1 → 4) = P(1/1)P(4/1) + P(2/1)P(4/2) + P(3/1)P(4/3) + P(4/1)P(4/4)

= (0.3)(0.1) + (0.3)(0.0) + (0.3)(0.2) + (0.1)(0.4) = 0.13 (6.68)

b) Si en t0 el estado de la partı́cula es 1, ¿En cuál estado la partı́cula está con mayor

probabilidad?
104 Chapter 6. Probabilidad

De las reglas de transición, calculamos las probabilidades restantes:

P(1 → 1) = 0.22
P(1 → 2) = 0.35
P(1 → 3) = 0.30 (6.69)

Lo que significa que la transición más probable es P(1 → 2).

Note que en el ejemplo anterior hay una estructura intrı́nseca matricial. Por ejemplo, si
pensamos en P(1 → 4) como la multiplicación de elementos de una matriz tenemos:

X
P(1 → 4) = p41 p11 + p42 p21 + p43 p31 + p44 p41 = p4k pk1 = P41 . (6.70)
k=4

Que corresponde con la multiplicación del cuarto renglón por la primera columna de la
matriz de probabilidades de transición.

Definición 10. Para un proceso de Markov, la matriz de probabilidad de transición

(estacionaria) de un paso es:
 
P00 P01 · · · P0n
 
 P10 P11 · · · P1n 
P= . (6.71)
 
 .. .. . . .. 
 . . . 

Pn0 Pn1 · · · Pnn
donde la información de la probabilidad de transición Pij debe entenderse como la
probabilidad de ir de la columna j al renglón i.

Esta matriz debe cumplir la propiedad básica de que la suma de probabilidades de cada
columna debe ser igual a uno:

X
Pij = 1, ∀j≤n (6.72)
i

Del ejemplo 1 de una cadena de dos pasos, la probabilidad de transición del estado 1 a todos
~
los estados, puede calcularse como el producto matricial P2 E:
6.14. Ejercicios: Cadenas de Markov 105

 2    
0.3 0.25 0.15 0.10 1 0.22
     
2~
0.3 0.25 0.55 0.20 0 0.35
P(1 → 1, 2, 3, 4) = P E = 
0.3 0.50 0.10 0.30 0 = 0.30
     (6.73)
     
0.1 0.00 0.20 0.40 0 0.13

1.0 Pn(0)
Pn(1)
0.8 Pn(2)
Pn(3)
0.6
Pn

0.4

0.2

0.0
0 2 4 6 8 10
n
Figure 6.3: Probabilidades de transición de los estados para n = 10 pasos.

Section 6.14
Ejercicios: Cadenas de Markov

1. Calcule las probabilidades de transición estacionarias del ejemplo 1.

Section 6.15
Metropolis-Hasting algorithm

Suppose you wish to sample the random variable x ∼ f (x), but you can not use the direct
simulation, the inverse CDF or accept-reject methods. But you can evaluate f (x) at least up
to a proportionality constant; then you can use the Metropolis-Hastings algorithm.

Let f (x) be the (possible unnormalized) target density, xj be a current value and q(x/xj )
be a proposal distribution, which might depend on the current value xj . The algorithm is
106 Chapter 6. Probabilidad

summarized as follows:

a) Sample x∗ ∼ q(x/xj ).

b) Calculate the acceptance probability:

f (x∗ )q(xj /x∗ )

 
j ∗
ρ(x , x ) = min 1, , (6.74)
f (xj )q(x∗ /xj )

if ρ(xj , x∗ ) > r, x∗ is accepted. Where r ∼ U(0, 1).

c) Set xj+1 = x∗ with probability ρ(xj , x∗ ), otherwise set xj+1 = xj .

In general, this sequence is not independent, moreover, we concentrate on symmetry chains

(random walks) where q(x/xj ) = q(xj /x). The acceptance probability is, therefore:

f (x∗ )
 
j ∗
ρ(x , x ) = min 1, (6.75)
f (xj )
How this algorithm works?

Let Dn (x) be the points density xn around x and Dn (y) the points density yn around y. x
will be accepted or rejected. We can quantify the flux of points on the neighborhood x as:

δDn (x) = δDn+1 (x) − δDn (x) (6.76)

X X
δDn (x) = Dn (y)P (y → x) − Dn (x)P (x → y) (6.77)
y y
| {z } | {z }
Gain Lost
 
X P (y → x) Dn (x)
Dn (y)P (x → y) − (6.78)
y
P (x → y) Dn (y)
We have two properties:

1. The equilibrium is achieved asymptotically.

2. There is flux toward y, which tends Dn (x) to the equilibrium.

P (y → x) Dn (x)
< (6.79)
P (x → y) Dn (y)
6.16. Ejercicios: Metropolis-Hastings 107

As a reminder ρ(x, y) > r is the acceptance condition. We can calculate P (x → y) as pxy Axy ,
where pxy is the probability of generation of y given x, and Axy is the acceptance probability.
Since the Metropolis estrategy is symmetric (pxy = pyx ), we get:

P (y → x) Ayx f (x)
= = ρ(x, y) = (6.80)
P (x → y) Axy f (y)
Therefore, the equilibrium distribution is given by:

Deq = D(x) = f (x) (6.81)

Section 6.16
Ejercicios: Metropolis-Hastings

1. Se lanza una moneda n = 10 veces y se encuentra que r = 7 veces cae cara. Usando el
algoritmo de Metrópolis:

a) Encuentre el parámetro asociado a la probabilidad de éxito p̂ de la distribución

binomial, es decir, encuentre el máximo de la distribución posterior.

b) Usando la varianza binomial y los cuantiles de la distribución posterior, encuentre los

+
errores asociados al parámetro p̂ a un nivel de confianza del 68% (i.e, σ− ).

c) ¿Podemos decir que la moneda está truncada?

Utilice la siguiente distribución a priori para el espacio de parámetros:

(
1 si 0 < p < 1
π(p) = (6.82)
0 si Otro caso

2. Resuelva el problema anterior usando el paquete Emcee.

[Link]

3. Usando el algoritmo de Metrópolis, realice el muestreo de N = 1000 eventos de una

distribución normal: A ∼ N (x; µ = 2, σ = 0.5).

4. Usando el algoritmo de Metrópolis, realice el muestreo de N = 1000 eventos de una

distribución estándar de Cauchy: A ∼ f (x; 0, 1).
108 Chapter 6. Probabilidad

1
f (x; 0, 1) = (6.83)
π(1 + x2 )
5. Genere un ensamble de al menos 30 caminatas aleatorias como las generadas en clase.
Calcule la distancia cuadrática desde el inicio del movimiento hasta el final, es decir
R2 = x2 + y 2 dado que se inicia siempre en [0,0] [6.4]. Calcule la distancia cuadrática
media para cada uno de los siguientes pasos:

N = {10, 100, 500, 1000, 5000, 10000, 20000, 100000} (6.84)

Deberı́a obtener algo como:

Figure 6.4: Relación entre la distancia cuadrática media y el número de pasos de un proceso
Markoviano. Entender este proceso fue uno de los pasos definitvos para aceptar la existencia
de las moléculas empezando el siglo XX.

6. Sea el volumen d−dimensional (Vd ):

Z ∞
Vd = cos(||r||)exp(−||r||2 )drd , (6.85)
−∞
qP
d
donde ||r|| = i=1 x2i y drd = dx1 · · · dxd , con el factor de normalización:

d
fn = (2πσ 2 ) 2 (6.86)
6.16. Ejercicios: Metropolis-Hastings 109

a) Usando el método de Metrópolis-Hasting para N = 105 eventos y σ = √1 , muestre

2
que los volúmenes d = 2, 3 son respectivamente V2 ≈ 1.817671646 y V3 ≈ 2.167233695.

7. La molécula de Hidrógeno está compuesta por dos núcleos separados una distancia L. En
torno a cada núcleo existe un electrón que supondremos que está en el nivel más bajo de
energı́a, en un órbital s. La función de onda que describe a cada electrón está dada por:

~ = p1 e−|~r−R|/a
ψ(~r, R)
~ 0
(6.87)
3
πa0

~ es la posición del núcleo respecto al sistema

donde a0 = 0.529 Å es el radio de Borh y R
de referencia. Adicionalmente, el operador de energı́a electrostática de la molécula (en
electronvoltios), que incluye la interacción entre electrones y núcleos está dado por:

~ 1, R
~ 2) = 1 1 1 1 1 1
U (~r1 , ~r2 , R + − − − −
~1 − R
|~r1 − ~r2 | |R ~ 2 | |~r1 − R
~ 1 | |~r2 − R
~ 1 | |~r1 − R
~ 2 | |~r2 − R~ 2|
(6.88)

Usar un sistema de unidades donde a0 = 1 y poner los núcleos en las posiciones: R ~1 =

~ 2 = [0, 0, −L/2]. La idea es calcular el valor esperado de la energı́a potencial
[0, 0, L/2] y R
como función de la separación de los núcleos de la molécula. Este promedio está dado
(sin considerar el espı́n electrónico) por:

Z
< U >= ~ 1 )|2 |ψ(~r2 , R
d3 r1 d3 r2 |ψ(~r1 , R ~ 2 )|2 U (~r1 , ~r2 , R
~ 1, R
~ 2) (6.89)

a) Note que esta integral se puede ver como el valor promedio de la energı́a potencial
evaluada en el muestreo de la densidad de probabilidad de cada electrón. Usando el
algoritmo de Metropolis en coordenadas cartesianas, hacer el muestreo del cuadrado
de la función de onda fijando L = 2 y 105 pasos en la cadena. La distribución de
puntos de la nube electrónica asociada a cada electrón está dada por [6.5]:

b) Calcule el valor esperado de la energı́a potencial para el siguientes separaciones

nucleares.
l = [Link](1.0,3.0,10). Deberı́a obtener algo similar a [6.7]:

c) Incluya el efecto de Heitler y London (buscar el sentido fı́sico):

110 Chapter 6. Probabilidad

Figure 6.5: Nube electrónica encontrada con el algoritmo de Metropolis asociada al electrón en
cada núcleo.

1.90

1.92

1.94
<U(r)>

1.96

1.98

2.00

2.02

2.04
1.00 1.25 1.50 1.75 2.00 2.25 2.50 2.75 3.00
r[A]

Figure 6.6: Valor esperado de la energı́a potencial eléctrica de la molécula de Hidrógeno.

Note que existe un estado ligado que conocemos como enlace covalente, la energı́a del pozo
de potencial está alrededor de Umin ≈ 2.0 eV y se da a una separación nuclear de dequilibrium ≈
2.1a0 .

Z
< U >= ~ 1 )ψ(~r2 , R
d3 r1 d3 r2 |ψ(~r1 , R ~ 2 ) + ψ(~r2 , R
~ 1 )ψ(~r1 , R
~ 2 )|2 U (~r1 , ~r2 , R
~ 1, R
~ 2 ) (6.90)

d) Note que los valores están de acuerdo a los observados en la naturaleza: La distancia
de separación de los protones es dequlibrium = 2.4a0 y la energı́a de ligadura es 2.8 eV .
Nice isn’t it!
6.17. Hidden Markov models 111

1.86

1.88

1.90

<U(r)>
1.92

1.94

1.96

1.00 1.25 1.50 1.75 2.00 2.25 2.50 2.75 3.00

r[A]

Figure 6.7: Valor esperado de la energı́a potencial eléctrica de la molécula de Hidrógeno

incluyendo el efecto de Heitler y London. Note que existe un estado ligado que conocemos
como enlace covalente, la energı́a del pozo de potencial está alrededor de Umin ≈ 2.05 eV y se
da a una separación nuclear de dequilibrium ≈ 2.3a0 .

Section 6.17
Hidden Markov models

Un proceso de Markov oculto es un proceso estocástico definido por dos variables aleatorias,
la variable oculta x(t) toma valores en el paso t y la variable observada y(t) toma valores en
el paso t. La variable x(t) toma valores siguiendo la matriz de transición T y la variable y(t)
toma valores siguiendo la matriz de emisión E. Por ejemplo, definimos las siguientes matrices
de transición y emisión:
 
F T
 
T=
 F 0.7 0.5
 (6.91)
T 0.3 0.5
 
F T
 
R 0.8 0.2
E=
  (6.92)
V 0.1 0.3


A 0.1 0.5
Las columnas de estas matrices representan las probabilidades de los estados en el tiempo
t − 1 y las filas las probabilidades asociadas al tiempo t. La suma de todos los valores en una
columna debe ser igual a 1. Definamos los estados de transición como el estados de animo de una
persona; Feliz y Triste como: S = [F, T ] con las siguientes probabilidades a priori π = [0.4, 0.6].
112 Chapter 6. Probabilidad

De esta forma, la probabilidad de estar feliz al dı́a (t − 1) y luego estar triste al siguiente dı́a (t)
es: T10 = 0.3. Los estados de transición que siguen este proceso aleatorio no son observables
por lo que se denominan ocultos. Sin embargo, esta probabilidad de transición se propaga
influyendo en la manera de vestir de la persona. Definamos los estados de emisión como el
color de la ropa que usa la persona al tiempo t − 1; Rojo, Verde, Azul como E = [R, V, A].
En este contexto, la probabilidad de que la persona use ropa azul dado que está feliz es E20 = 0.1.

La secuencia de n eventos observados es:

ΩO = {y(t1 ), y(t2 ), ..., y(tn )} (6.93)

Usando la propiedad de los procesos de Markov, se busca encontrar la secuencia de eventos

ocultos más probable.

Ωhidden = {x(t1 ), x(t2 ), ..., x(tn )} (6.94)

Suponga la siguiente secuencia de ropa usada por la persona: ΩO = [V, A, R].

a) Encuentre la probabilidad de cada secuencia de estados (ocultos) de ánimo.

b) Encuentre el estado de animo más probable que tuvo la persona durante esos tres dı́as.

c) Calcule las probabilidades de cada estado observable (i) como la suma de las probabilidades
de todos los estados ocultos. Verifique que sumando todos los estados observables la medida
es 1.
X
Pi = 1 (6.95)

Remarks:
nh es el número de estados ocultos y no es el número de estados observados. El número de
secuencias posibles (ns ) está dada por las variaciones con repetición de estas cantidades:

ns = nnho = 8 (6.96)

De las ocho secuencias ocultas posibles, supongamos que la secuencia oculta es S1 =

[F, F, T ]. Entonces:
!
F F T
S1 = (6.97)
V A R
6.17. Hidden Markov models 113

0.0175
0.0150
0.0125
0.0100 Probabilidad por secuencia
MaxP = 0.018
0.0075
0.0050
0.0025
0.0000
1 2 3 4 5 6 7 8

Figure 6.8: Probabilidad de cada secuencia oculta para el estado ΩO = [V, A, R].

Recordemos la probabilidad condicional de dos sucesos:

P(A ∩ B) = P(A/B)P(B) (6.98)

Recordemos las matrices involucradas:

 
F T
 
T = F 0.7 0.5

 (6.99)
T 0.3 0.5
 
F T
 
R 0.8 0.2
E=   (6.100)
 V 0.1 0.3 

A 0.1 0.5
Por tanto, la probabilidad de la secuencia es:

P(V, A, R ∩ F, F, T ) = P(V /F )P(A/F )P(R/T ) × P(F ) P(F/F )P(T /F )

| {z } | {z } | {z }
Emission Matrix Prior Transition Matrix
= E10 E20 E01 P0 T00 T10
= 0.1 · 0.1 · 0.2 × 0.4 · 0.7 · 0.3 = 1.68 × 10−4 (6.101)

En general:

n
Y n
Y
P(Oi , qi ) = P(Oi /qi ) P(qi /qi−1 ) (6.102)
i=1 i=1
114 Chapter 6. Probabilidad

Section 6.18
Ejercicios: Hidden Markov models

1. Un casino tiene dos tipos monedas una justa (J) y una sesgada (B), las probabilidades
de obtener cara y sello son: PJ = [0.5, 0.5] y PB = [0.9, 0.1] respectivamente. El tipo de
moneda se escoge siguiendo esta ley de transición:

 
J B
 
T =  J 0.8 0.2

 (6.103)
B 0.2 0.8

Por otro lado, la matriz de emisión está dada por:

 
J B
 
E=
 C 0.5 0.9
 (6.104)
S 0.5 0.1

Se realiza un experimento de 8 lanzamientos y se obtiene la siguiente secuencia:

ΩO = [S, C, C, C, S, C, S, C] (6.105)

a) Use la siguiente distribución de probabilidad a-priori π = [0.2, 0.8] para la moneda

justa y sesgada.

b) Encuentre la secuencia oculta más probable del tipo de moneda que se eligió en cada
lanzamiento y su respectiva probabilidad Pi .

c) Calcule las probabilidades de cada estado observable (o) como la suma de las
P
probabilidades de todos los estados ocultos, Po = i Pi .

d) Verifique que la suma de todos los estados observables es 1.

X
Po = 1 (6.106)

e) ¿Depende el resultado de la probabilidad a-priori?

6.18. Ejercicios: Hidden Markov models 115

0.000200
Probabilidad por secuencia
0.000175 MaxP = 0.0002
0.000150
0.000125
0.000100
0.000075
0.000050
0.000025
0.000000
0 50 100 150 200 250

Figure 6.9: Probabilidad de cada secuencia oculta para el estado observado ΩO =

[S, C, C, C, S, C, S, C].

2. Las bases nitrogenadas fundamentales que componen el ADN son: Adenina (A), Citosina
(C), Guanina (G) y Timina (T). Un gen se puede representar a través de una secuencia
ordenada de dichas bases. Suponga la siguiente matriz de transición entre bases:

 
A C G T
 
A 0.4 0.2 0.2 0.2 
 
T = C 0.25 0.25 0.25 0.25 (6.107)
 
 
G 0.3 0.3 0.1 0.3 
 
T 0.1 0.1 0.1 0.7

La probabilidad a-priori está dada por:

π = [0.25, 0, 0.5, 0.25] (6.108)

Encuentre la probabilidad de obtener el siguiente gen:

g = [T, G, C, T, C, A, A, A] (6.109)

Pg = 7.5 × 10−6
116 Chapter 6. Probabilidad
Chapter 7

Estadı́stica

La estadı́stica estudia la organización y análisis de datos para dar conclusiones y tomar

decisiones válidas. Es imposible estudiar el grupo completo de datos denominado población, en
su lugar, se toma un sub-conjunto llamado muestra (que debe ser representativo). En general,
la población puede ser finita o infinita. La descripción y análisis de unicamente la muestra,
se denomina estadı́stica descriptiva. Si por el contrario se desea inferir conclusiones sobre la
población, con base a la muestra, se denomina estadı́stica inferencial.

Section 7.1
Muestreo (sampling)

Debido a que en general es imposible tener la población completa, se requiere una metodologı́a
sistemática y que sea imparcial para elegir elementos de una población. La técnica para elegir
una muestra de una población de denomina muestreo. Suponiendo un tamaño de población N
y de muestra (n) determinado, el muestreo debe permitir hacer el estudio y estimaciones de las
caracterı́sticas de la población, por ejemplo, promedios, desviaciones tı́picas, etc. La muestra
debe ser representativa para generalizar los resultados a la población.

Ventajas:

1. Permite tener un sub-conjunto de toda la población.

2. Si la muestra es representativa, el muestreo permite inferir los resultados de la muestra a

la población.

117
118 Chapter 7. Estadı́stica

3. Cuesta menos dinero.

4. Se pueden obtener resultados más rápido.

Desventajas:

1. Existe posibilidad de sesgo, i.e, preguntar sobre la preferencia de un candidato presidencial

a un solo tipo de población.

2. Margen de error.

3. Dificultad para obtener los datos.

Se tienen las siguientes técnicas de muestreo:

7.1.1 Aleatorio (probabilı́stico)

Es un método adecuado para seleccionar una muestra lo suficientemente representativa de una
población, esta información se puede inferir al conjunto completo con cierto nivel de confianza.
Este tipo de muestreo tiene las siguientes caracterı́sticas:

a) Se basa en que todos los elementos de la población tiene la misma probabilidad de ser
elegidos.

b) Este tipo de muestreo garantiza la representatividad de la muestra. (procesos de inferencia

válidos).

c) Permite la estimación de caracterı́stica poblacionales (promedios, total, proporciones, etc).

d) Las estimaciones se realizan a través de funciones matemáticas denominadas estimadores

estadı́sticos.

e) Los errores se cuantifican mediante varianzas y desviaciones de los estimadores estadı́sticos.

En particular, se tienen las siguientes técnicas aleatorias:

7.1. Muestreo (sampling) 119

7.1.2 Muestreo aleatorio simple

Es método consiste en seleccionar una muestra de tamaño n, cuyos elementos tienen igual
probabilidad dentro de una población N . Esta muestra se obtiene elemento por elemento de
forma aleatoria, adicionalmente, la construcción de la muestra no dependen del orden en que
se obtiene la secuencia. Los elementos pueden ser elegidos con y sin reposición, cuando hace
el muestreo con reposición, se puede obtener el mismo elemento más de una vez, en cuyo caso
todas las muestras son equiprobables. Cuando realiza el muestreo sin reposición, cada elemento
puede aparecer una sola vez en la muestra y cada muestra es equiprobable. La probabilidad de
que un elemento de la población pertenezca a una muestra es:

Número de muestras que contienen el elemento i

Pi = (7.1)
Número de muestras totales

N −1

(N −1)!
n−1 (n−1)!(N −n)! n
Pi = N

= N!
= (7.2)
n n!(N −n)!
N

7.1.3 Muestreo aleatorio estratificado

Es posible construir la muestra a través de estratos o categorı́as (L) que pueden ser etiquetadas
como: n1 , n2 , ..., ni , ..., nL Este esquema de muestreo está influenciado principalmente por tres
factores: el número total de elementos en cada estrato, la variabilidad de las observaciones
y en casos reales, por el costo de adquirir la muestra. El muestreo estratificado es óptimo
y homogéneo, cuando se minimiza la desviación estándar de cada estrato, sin embargo, esta
minimización resulta difı́cil de alcanzar en los casos reales. Para el muestreo estratificado,
tenemos asignación igual o proporcional. Tomemos el siguiente ejemplo: Para cierta universidad
tenemos el siguiente número de estudiantes por carrera:

Carrera Estrato Población

Medicina 1 1000
Ing. Industrial 2 500
Mercadeo 3 450
Fı́sica 4 50

Table 7.1: Distribución de estudiantes según su carrera universitaria.

120 Chapter 7. Estadı́stica

Se requiere formar una muestra aleatoria estratificada de n = 400 usando asignación igual
y proporcional. Para asignación igual, se tiene el siguiente tamaño de estrato:

400
ni = = 100 (7.3)
4
Para asignación proporcional tenemos los siguientes tamaños de estrato.

Ni
Carrera Estrato Población (Ni ) Muestra ni = n × N
Medicina 1 1000 200
Ing. Industrial 2 500 100
Mercadeo 3 450 90
Fı́sica 4 50 10
N = 2000 n = 400

Table 7.2: Distribución muestral según su carrera universitaria.

7.1.4 Muestreo aleatorio sistemático

Este método consiste en obtener los elementos de la muestra de una manera ordenada, donde
los elementos se seleccionan en intervalos regulares. En particular, este método arroja buenos
resultados cuando la variable de interés está ordenada de mayor a menor. Si la disposición
de los elementos es aleatoria, el muestreo sistemático equivale al método aleatorio simple. El
intervalo de selección de los elementos está dado por:

N
k= (7.4)
n
El primer valor (S0 ) se elige aleatoriamente de la lista entre el primer valor y k. Para
seleccionar los siguientes k − 1 cantidades se realiza una traslación sistemática de k unidades
desde el primer valor:

i0 = [Link](k)
in+1 = in + nk, con n = 1, 2, ..., k − 1 (7.5)

Tomemos la siguiente población:

7.2. Distribución de frecuencias (Histogramas) 121

Ω = [0, −1, −3, −5, 5, 7, 4, 10, 12, 15, 20, 19, 50, 55, 60, 45] (7.6)

Tomando k = 16/4 = 4 y empezando con la lista ordenada en Ωsorted

1 = 55, tenemos la
muestra de tamaño k = 4:

S = [55, 19, 7, −1] (7.7)

7.1.5 Muestreo no aleatorio

Son técnicas de muestreo donde la muestra se obtienen en procesos donde cada individuo no
tienen la misma probabilidades de ser seleccionado.

1. Investigaciones exploratorias (por cuotas, intensional, bola de nieve).

2. La forma de selección, hace que la muestra no sea representativa. (procesos de inferencia

no aplicables).

Section 7.2
Distribución de frecuencias (Histogramas)

Una de las operaciones más sencillas de realizar en un conjunto de datos es la ordenación en

términos de la frecuencia de repetición de los datos. Para tal fin, vamos a definir algunos
parámetros de la muestra que nos permite ordenar la información.

1. Rango de la población: Es la diferencia entre el mayor y el menor valor encontrado en la

muestra. Por ejemplo, se mide el peso de 300 estudiantes donde el valor mayor es 100 kg
y el menor es 65.4 kg. El rango de la muestra es: 34.6 kg.

2. Frecuencia de clase: Es apropiado definir categorı́as o clases (Ci ) en los datos para separar
la muestra en n clases. El número de datos pertenecientes a la clase Ci se denomina
frecuencia de clase. Las clases tienen las siguientes propiedades:

(a) Ci 6= ∅

(b) Ci ∩ Cj = ∅

(c) ni=1 Ci = Ω
S
122 Chapter 7. Estadı́stica

3. Lı́mite de clase: En el ejemplo del peso de los estudiantes, si tomáramos la clase C1

como todos los estudiantes que pesan entre 65 y 75 kg, entonces los lı́mites de clase serı́a
precisamente dichos valores.

4. Tamaño de clase: Es la diferencia entre los lı́mites de clase, por ejemplo: 75 - 65 kg = 10

kg es el tamaño de clase de C1 .

5. Marca de clase: Es el punto medio de cada clase, el cuál se obtiene promediando el valor
de los lı́mites de clase.

Con estas definiciones podemos definir el concepto de histograma. Un histograma es el

conjunto de clases, que pueden ser representadas gráficamente como los rectángulos, cuya base
esta centrada en las marcas de clase y su altura es la frecuencia de clase. Note que, el área
total histograma es proporcional al tamaño total de la muestra. En ciertos casos resulta muy
útil hacer que el área del histograma sea igual a 1 para comparar la forma de la distribución
con otras distribuciones.

Ası́ mismo, podemos definir una distribución de frecuencias acumulada, que nos da
información de como se acumula la información en los datos para valores menores a una marca
de clase especı́fica. La grafica de la distribución acumulada de frecuencia se denomina ojiva.
Imaginemos la siguiente situación:

Peso [kg] Frecuencia Acumulada

60 - 70 6 Menor a 60 0
70 - 80 35 Menor a 70 6
80 - 90 40 Menor a 80 41
90 - 100 19 Menor a 90 81
Menor a 100 100

Table 7.3: Distribución de pesos y distribución acumulada de pesos para un conjunto de

mediciones de 100 estudiantes.

Una descripción sencilla de la distribución acumulada es que 81 estudiantes tienen un peso

menor a 90 kg.
7.3. Estadı́sticos 123

Section 7.3
Estadı́sticos

Sea X una variable aleatoria y sean Ω = {X1 , X2 , ..., Xn } variables aleatorias con la misma
distribución que X. Diremos que el conjunto Ω es una muestra aleatoria de tamaño n de X.
Dada una muestra de una población, un estadı́stico es una función real de la muestra:

T = f (X1 , X2 , ..., Xn ) (7.8)

como el estadı́stico es una variable aleatoria, entonces será en si mismo una variable aleatorio.
Esto significa que tendrá una distribución, una media, etc. En general, la distribución que siga
el estadı́stica se denomina distribución muestral.

Section 7.4
Valor medio muestral

En el caso de la media aritmética si los números x1 , x2 , ..., xn ocurren f1 , f2 , ..., fn veces,

respectivamente. La media muestral está dada por:
P
f1 x1 + f2 x2 + ... + fn xn f i xi
X̄ = = P . (7.9)
f1 + f2 + ... + fn fi
Del mismo modo, es posible definir pesos estadı́sticos que ponderen cada punto muestral.
{w1 , w2 , ..., wn } tienen las siguiente caracterı́sticas:

n
X
wi = 1
i=1
wi ≥ 0 ∀i , i = 1, 2, ..., n (7.10)

La media ponderada está dada por:

n
w1 x1 + w2 x2 + ... + wn xn X
X̄ = = w i xi (7.11)
w1 + w2 + ... + wn i=1

Esta definición resulta muy útil en el caso de la media móvil, dado que, para todo punto i
de la secuencia de observaciones y n puntos anteriores, la media móvil del siguiente punto i + 1
queda definida por:
124 Chapter 7. Estadı́stica

n−1
X
m
X̄i+1 (n) = wj xj+i−n+1 (7.12)
j=0

donde i ≥ n. Esto significa que no es posible calcular la media móvil X̄2m (3), dado que en el
punto muestral x2 solo se cuentan con dos observaciones anteriores {x1 , x2 }. Entonces, para los
puntos i < n se mantienen los valores observados. Adicionalmente, es posible optimizar el valor
de los pesos estadı́sticos para tener la mejor estimación posible. Para realizar la optimización
se define un función de costo que tenga en cuenta la distancia entre las observaciones y las
predicciones de la media móvil en cada punto muestral. La función de costo está dada por:

N
1 X m
C(n) = |X̄ (n) − xi | (7.13)
N i≥n i

La optimización se realiza bajo las restricciones dadas por la ecuación (7.10); usando algún
paquete de minimización (i.e, [Link]).

La media armónica es útil en algunos casos particulares:

1 1 X1
= (7.14)
H N x

Section 7.5
Desviación estándar

1. La desviación estándar nos da una medida de la dispersión de los datos alrededor de la

media. Para el caso poblacional está dada por:

v
u
u1 X N
p
σ= t (xi − µ̂)2 = E(x2 ) − (E(x))2 (7.15)
N i=1

En el caso muestral, la desviación estándar está dada por:

v
u N
u 1 X
S=t (xi − x̄)2 (7.16)
N − 1 i=1
7.6. Teorema del lı́mite central 125

Section 7.6
Teorema del lı́mite central

Sean X1 , X2 , ..., Xn variables aleatorias independientes e idénticamente distribuidas con E(X) =

µ y Var(X) = σ 2 < ∞. Se define:

1
Pn
n
Xn − µ
i=1
Zn := √ , (7.17)
σ/ n
La función de distribución de Zn converge hacia la distribución normal estándar cuando
n → ∞. Entonces:
Z z
1 −x2
lim P(Zn < z) = Φ(z) = √ e− 2 dx. (7.18)
n→∞ −∞ 2π
El teorema exige la existencia de la media y la varianza, sin embargo, no especı́fica la
distribución de la variable X; es un resultado general [4].

Section 7.7
Ejercicios: Muestreo

1. Use los datos almacenados en: [Link]

Herramientas/blob/main/Data/Sesion2/[Link] para realizar un
muestreo estratificado con n = 200 y n = 300. Verifique que la suma total de personas
encuestadas en cada salón sea n.

2. Guarde la información final en el formato .xls para su distribución.

3. Use los datos almacenados en: [Link]

Herramientas/blob/main/Data/Sesion2/Datos_Sistematico.csv para realizar un
muestreo sistemático con n = 200.

4. Guarde la información final en el formato .xls para su distribución.

5. La calificación de matemáticas de 80 estudiantes están en el siguiente archivo: https://

[Link]/asegura4488/DataBase/blob/main/MetodosComputacionalesReforma/Matematicas.
txt. Hallar la siguiente información:

a) La calificación más alta.

126 Chapter 7. Estadı́stica

b) La más baja.
c) El rango.
d) La cinco más bajas.
e) Las cinco más altas.
f) La décima de mayor y menor.
g) El número de estudiantes con calificación de 75 o más.
h) El porcentaje de estudiantes con calificaciones mayores que 65 pero no superiores a
85.
i) Las calificaciones de 0 a 100 que no aparecen en las notas.

6. La Tabla [7.4] muestra la distribución de frecuencia de los salarios de los empleados de

una empresa anónima. Calcule la distribución acumulada de porcentajes y su respectiva
gráfica.

Salarios Número de empleados

250 - 259 8
260 - 269 10
270 - 279 16
280 - 299 14
290 - 300 10
300 - 309 5
310 - 319 2

Table 7.4: Distribución de salarios para 65 empleados de la empresa anónima.

7. Usando la distribución (empı́rica) de notas de matemáticas vista en clase, hacer las

siguientes estimaciones:

a) Calcular la media de la nota con [Link](data).

b) Calcular la mediana de la nota con [Link](data,50)
c) Calcular la distribución de frecuencia acumulada y la distribución acumulada de
probabilidad P(X <= x). En este caso, la mediana es el valor que toma la variable
aleatoria que divide en 2 partes iguales la función de masa de probabilidad p(X = x).
7.7. Ejercicios: Muestreo 127

d) Usando la técnica de re-muestreo (Bootstrapping), obtenga la función de distribución

de la media con [Link], donde el tamaño de muestra sea igual al poblacional.
En este caso, se analiza la variabilidad de la población respecto a la media.

e) Hacer el fit gaussiano para obtener el valor de la media y su desviación estándar

(Figura [7.1]).

f) Calcular la mediana de la distribución bootstrapped usando la función acumulada de

probabilidad (Figura [7.1]).

g) Comparar todas las estimaciones.

Fit results: = 75.25, = 1.15

0.35 Gaussian fit 1.0 cdf boostrapped
Bootstrapping dist
0.30
0.8
0.25
0.6
0.20

0.15 0.4

0.10
0.2
0.05
0.0
0.00
70 72 74 76 78 80 70 72 74 76 78 80
x

Figure 7.1: Distribución de muestreo (izquierda) y distribución acumuladad de probabilidad

(derecha).

8. Encuentre los valores óptimos de los pesos asociados a la predicción de las ventas. La
pregunta puede ser muy complicada de responder!

9. Una persona viaja de A a B con una velocidad media de 30 millas por hora y regresa de
B a A a una velocidad de media de 60 millas por hora. Hallar la velocidad media en el
viaje completo. Ans: 40 mi/h
128 Chapter 7. Estadı́stica

10. Un avión vuela d1 , d2 , d3 millas a velocidades v1 , v2 , v3 mi/h, respectivamente. Demuestre

que su velocidad media es V, dada por:

d1 + d2 + d3 d1 d2 d3
= + + (7.19)
V v1 v2 v3
es una media armónica ponderada.

11. Considere la variable aleatoria definida por una combinación lineal de otras variables
aleatorias: X = a1 X1 + a2 X2 + ...an Xn , donde las variables aleatorias (Xi ) no son
necesariamente idénticamente distribuidas, y las componentes [a1 , a2 , ..., an ] ∈ R. La
~ con la definición usual de producto
variable aleatoria X puede ser escrita como: X = aT X
interno entre vectores. El primer y segundo momento de la variable X queda entonces
expresados por:

~ = aT E(X)
E(X) = E(aT X) ~
~ = aT Cov(X)a
V ar(X) = V ar(aT X) ~ (7.20)

~ es la matriz de covarianza de X.
Donde Cov(X) ~ Sea X1 ∼ Γ(2, 3), X2 ∼ N (5, 2) y
X3 ∼ U (0, 10), Genere N = 104 eventos (que estabilice el valor de los momentos) para
obtener la distribución de X = X1 + 2X2 − X3 . Calcule el primer y segundo momento de
X a través de dos estrategias:

a) Usando directamente el array de X.

b) Usando las definiciones generales dadas en la Ecuación (7.20).
c) Calcule el coeficiente de correlación de Pearson para las tres variables.
d) Demuestre que la formula de la varianza de la media se puede escribir escalarmente
como:

N N N N
1 X 1 X 2 X X
V ar( Xi ) = 2 V ar(Xi ) + 2 Cov(Xi , Xj ) (7.21)
N i=1 N i=1 N i=1 j=i+1

Hint: Encuentre el resultados para N = 2 y use la intuición adquirida para el caso

general.
Chapter 8

Estimación de parámetros

Section 8.1
Mı́nimos cuadrados

El método de mı́nimos cuadrados surge como necesidad de plantear una solución aproximada
al problema lineal A~x = ~b, donde la solución exacta no existe. Para entender el planteamiento
desde el álgebra lineal se requiere el concepto de espacio columna de A.

8.1.1 Espacio Columna

Sea A un matriz genérica de dimensiones m × n, que se puede escribir como un arreglo de
vectores columna ~vm×1 :

 
Am×n = ~v1 , ~v2 , ... ~vn (8.1)

Cada uno se esos vectores tiene m componentes, por tanto ~vi ∈ Rm , donde i = 1, 2, ...n.
Entonces se puede definir al espacio columna de la matriz C(A), como todas las combinaciones
lineales de los vectores columnas de A. Tomemos un vector genérico de ese subespacio:

~a = c1~v1 + c2~v2 + ... + cn~vn (8.2)

Notemos que el vector cero pertenece al espacio columna, dado que es una combinación
lineal válida con ci = 0, con i = 1, 2, ...n. Adicionalmente, vemos que es cerrado bajo la
multiplicación por escalar.

129
130 Chapter 8. Estimación de parámetros

k~a = kc1~v1 + kc2~v2 + ... + kcn~vn (8.3)

~b = d1~v1 + d2~v2 + ... + dn~vn (8.4)

Entonces ~b ∈ C(A). Finalmente, se puede demostrar que el espacio columna es cerrado

frente a la suma. Dado que, ~a ∈ C(A) y ~b ∈ C(A) entonces ~a + ~b ∈ C(A).

Por otro lado, para que el producto A~x esté bien definido, ~x ∈ Rn . Usando la representación
de la matriz A se puede escribir:

A~x = x1~v1 + x2~v2 + ... + xn~vn , (8.5)

dado que estamos tomando todas las combinaciones lineales posibles por que ~x es un vector
genérico, la expresión anterior debe ser el espacio vectorial generado por los vectores columnas
de la matriz A, es decir, C(A). En el contexto de los mı́nimos cuadrados, el sistema lineal
A~x = ~b no tiene solución por que ~b ∈
/ C(A).

Una vez demostrado que el vector ~b ∈

/ C(A), se quiere encontrar un método que permita
encontrar una solución aproximada lo más exacta posible. Entonces definimos A~x∗ como la
proyección del vector ~b en el espacio columna:

P royC(A) (~b) := A~x∗ . (8.6)

Por construcción la distancia entre la proyección de ~b en el espacio columna y el vector ~b,

está dada por:

d~ = A~x∗ − ~b (8.7)

El método de mı́nimos cuadrados encuentra el vector ~x∗ que minimiza esta distancia.

min(||A~x∗ − ~b||) (8.8)

Note que explı́citamente se tiene:

8.1. Mı́nimos cuadrados 131

 2
P − b1
 1 
P − b 
 2 2
 
 . 
~ =
|d|2  (8.9)
 . 
 
 
 . 
 
P n − bn

Adicionalmente, note que el vector d~ es ortogonal a C(A) y por tanto, d~ ∈ C ⊥ (A) que es el
complemento ortogonal de A. Se tiene:

C ⊥ = {~x ∈ Rn | ~x · ~v = 0 ∀~v ∈ C} (8.10)

El complemento ortogonal es el espacio nulo izquierdo de la matriz transpuesta de A:

C ⊥ (A) = N (AT ) (8.11)

Recordemos que el espacio nulo de una matriz está dado por:

N (A) = {~x ∈ Rn | A~x = ~0}. (8.12)

proof: Formemos una matriz cuyas filas son los n vectores de la base de C(A) (es decir,
T
A ):
 
−~v −
 1 
 −~v2 − 
AT =  .  (8.13)
 
 .. 
 
−~vn −
de la definición de C ⊥ (A), si ~x ∈ C ⊥ (A) entonces debe ser perpendicular a todos los
elementos de la base de C(A):

~x · ~v1 = ~x · ~v2 = · · · = ~x · ~vn = ~0. (8.14)

esta expresión es la multiplicación matricial AT ~x = 0, que es la definición de espacio nulo

de AT :

N (AT ) = {~x ∈ Rn | AT ~x = ~0}. (8.15)

132 Chapter 8. Estimación de parámetros

Por tanto, C ⊥ (A) es el espacio nulo de la matriz cuyas filas son una base de C(A); es decir,
N (AT ). 

Usando el último resultado se tiene:

AT (A~x∗ − ~b) = 0 (8.16)

AT A~x∗ = AT~b (8.17)

Section 8.2
Varianza de la estimación de parámetros

En el caso de mı́nimos cuadrados, es posible encontrar expresiones cerradas y útiles

computacionalmente para estimar la varianza de los parámetros encontrados usando este
método. Se parte de la definición general de varianza:

n
2 1 X
σ̂ = (yi − ŷi )2
n − p i=1
n
1 X
= (yi − Ai β̂)2
n − p i=1
1
= (Y − Aβ̂)T (Y − Aβ̂) (8.18)
n−p

donde p es el número de parámetros que se ajustan en el modelo y β̂ son los parámetros

encontrados. De esta forma, la desviación tı́pica asociada a los parámetros está dada por:

√
σ̂ = σ̂ 2 . (8.19)

Más generalmente, la matriz de covarianza de los parámetros está dada por:

V ar(β̂) = σ̂ 2 (AT A)−1 . (8.20)

8.3. Función de costo 133

Section 8.3
Función de costo

En el contexto de la estimación de parámetros es necesario introducir funciones de costo. En un

marco conceptual Bayesiano, las funciones de costo miden una distancia entre los datos (que se
consideran fijos) y el mejor modelo ajustado, por tal razón, la función de costo es una función
de los parámetros a estimar. La función más sencilla y directamente relacionada con el método
de mı́nimos cuadrados es χ2 (θ), la cual está definida por:

N  2
2
X yi − M (xi , θ)
χ (θ) = . (8.21)
i=1
σi

donde M (xi , θ) es el modelo que se propone para ajustar los datos. Los parámetros que
mı́nimizan la función de costo, son los valores óptimos que describen el mejor modelo que se
ajusta a los datos. En general, están dados por:

N  2
2
X yi − M (xi , θ)
θ̂ = arg min χ (θ) = arg min . (8.22)
θ θ
i=1
σi

Section 8.4
Ejercicios: Mı́nimos cuadrados

1. Se tienen tres lı́neas en R2 :

2x − y = 2
x + 2y = 1
x+y = 4 (8.23)

(a) Con el método de mı́nimos cuadrados encuentre el punto común a las tres lı́neas.
Grafique las tres lı́neas y el punto solución, ¿qué interpretación puede dar?

(b) Realice una búsqueda iterativa entre −5 ≤ x ≤ 5 y −5 ≤ y ≤ 5 con un paso de

h = 0.01 para encontrar la menor distancia del problema. Grafique la distancia y
compare con el resultado obtenido con mı́nimos cuadrados.
134 Chapter 8. Estimación de parámetros

min(||A~x∗ − ~b||) (8.24)

20.0
17.5
15.0
d(x* )

12.5
10.0
7.5
5.0

4
2 4
0 2
y 2 0
2 x
4 4

Figure 8.1: Distancia ||A~x∗ − ~b|| en la región rectangular definida en el problema. La distancia
mı́nima debe corresponder con el resultado del método de mı́nimos cuadrados.

2. Descargue los datos: [Link]

[Link]. Realice el ajuste usando el método de mı́nimos cuadrados para
encontrar los parámetros de:

f (x) = a0 + a1 x (8.25)

Grafique los datos y el ajuste mostrando el valor de las constantes en la etiqueta de la

grafica.

3. Descargue los datos: [Link]

[Link]. Realice el ajuste usando el método de mı́nimos cuadrados para
encontrar los parámetros de:

f (x) = a0 + a1 x + a1 x2 (8.26)

Grafique los datos y el ajuste mostrando el valor de las constantes en la etiqueta de la

grafica.
8.4. Ejercicios: Mı́nimos cuadrados 135

4. Utilice el método curve fit de Python para obtener los dos ajustes. Compare con los
resultados anteriores.

5. Encontrar el intervalo de confianza a 1σ = 68% de los parámetros del ejercicio de mı́nimos

cuadrados de Bootstrapping.

a) Haga un remuestreo con reemplazo con un tamaño de muestra igual al tamaño de los
datos. resample = [Link](sample, size=len(sample),replace=True)
b) Realice el ajuste del remuestreo usando los pasos del problema anterior.
c) Guarde el valor de cada parámetro en una lista diferente. Realice la experiencia
aleatoria N = 105 veces.
d) Dibuje el histograma asociado a la distribución de muestreo de cada parámetro.
e) De la distribución de muestreo estime el intervalo de confianza al 68% de cada
parámetro del ajuste. ¿Cuál es más preciso?

6. En el caso de ajustes, es posible definir funciones de costo que midan la distancia entre
los puntos muestrales y el modelo de ajuste. En el caso de mı́nimos cuadrados la función
es χ2 . Para n puntos y un modelo lineal la función de costo es:

n
X
2
χ (a0 , a1 ) = (yi − (a0 + a1 xi ))2 (8.27)
i=1

Si hablamos en terminos Bayesianos, Ω = {a0 , a1 } define el conjunto de modelos

lineales que pueden explicar los n puntos muestrales. Al minimizar χ2 (a0 , a1 ) muestre
(analı́ticamente) que los parámetros están dados por:

a0 = ȳ − a1 x̄
P P
xy − xn y
P
a1 = P P 2 (8.28)
x2 − ( nx)

donde x̄ y ȳ son los valores medios de puntos y sus imágenes. Para n puntos y un modelo
cuadrático la función de costo es:

n
X
2
χ (a0 , a1 , a2 ) = (yi − (a0 + a1 xi + a2 x2i ))2 (8.29)
i=1
136 Chapter 8. Estimación de parámetros

Minimize χ2 (a0 , a1 , a2 ) para encontrar el siguiente sistema de ecuaciones. Nota alguna

regularidad?

n
X
[a0 + xi a1 + x2i a2 = yi ]
i=1
n
X
[a0 xi + a1 x2i + a2 x3i = xi yi ]
i=1
Xn
[a0 x2i + a1 x3i + a2 x4i = x2i yi ] (8.30)
i=1

Opcionalmente: Encuentre los parámetros usando estas expresiones para los problemas
2) y 3)

7. (Machine Learning: Logistic Regresion) Descargue los datos de: [Link]

[Link]/asegura4488/Database/main/MetodosComputacionalesReforma/
[Link]

a) Defina el modelo de ajuste como:

~ = θ0
M (x; θ) (8.31)
θ1 + e−θ2 x

donde el θ~ es el vector de parámetros del ajuste.

b) Defina la métrica (función de costo) a minimizar como:

N  ~ 2
~ =
2
X yi − M (xi , θ)
χ (θ) . (8.32)
i=1
σi
donde σi = 1 ∀i, es decir, no se consideran los errores de yi .
c) Muestre que las derivadas parciales de la métrica están dadas por:

~ N ~
∂χ2 (θ) X
~ ∂M (xi , θ)
= −2 (yi − M (xi , θ)) (8.33)
∂θi i=1
∂θi

d) Muestre que vectorialmente, el descenso del gradiente queda definido por:

 N
X 
θ~j+1 = θ~j − γ −2 ~ ~
(yi − M (xi , θj ))∇θ M (xi , θj ) (8.34)
i=1
8.5. Función de verosimilitud (Likelihood) 137

donde el gradiente es respecto a los parámetros:

h i
~ = ~
∂M (xi ,θ) ~
∂M (xi ,θ) ~
∂M (xi ,θ)
∇θ M (xi , θ) ∂θ0
, ∂θ1
, ∂θ2
(8.35)

e) Use una taza de aprendizaje γ = 1 × 10−3 o γ = 5 × 10−4 , θ~0 = [1, 1, 1], un error de
parada  = 0.01 y un máximo de iteraciones de 1 × 104 .
f) Grafique los datos y the best fit model con sus parámetros.

8. Calcule la proyección ortogonal del vector ~b = (−3, −3, 8, 9) sobre el sub-espacio W

generado por los vectores:

~u1 = (3, 1, 0, 1)
~u2 = (1, 2, 1, 1)
~u3 = (−1, 0, 2, −1) (8.36)

a) Usando mı́nimos cuadrados matriciales. La proyección ortogonal es pW (b) = Ax, donde

las columnas de A son los vectores base y x es la solución de mı́nimos cuadrados.
b) Con el proceso de Grand-Schmidt obtener una base ortonormal (v1 , v2 , v3 ) y luego
calcular la proyección sobre dicha base: pW (b) = c1 v1 + c2 v2 + c3 v3 , donde ci =< [Link] >
para i = 1, 2, 3. Respuesta: pW (b) = (−2, 3, 4, 0)

Recuerde que el procedimiento de Grand-Schmidt es:

k−1
X < vk , uj >
uk = vk − uj . (8.37)
j=1
< uj , uj >

Section 8.5
Función de verosimilitud (Likelihood)

Sea ~x(x1 , x2 , ..., xn ) un vector aleatorio, donde cada componente está definida como un canal
de observación. Cada canal está asociado a una función de densidad de probabilidad que
depende de m parámetros θ1 , θ2 , ..., θm dado un conjunto de datos, la función de verosimilitud
está definida como:
138 Chapter 8. Estimación de parámetros

L(~x; θ1 , θ2 , ..., θm ) = fθ1 ,θ2 ,...,θm (~x). (8.38)

la cual es la función de densidad de probabilidad evaluada en el vector aleatorio ~x. Si la

muestra de datos es independiente e idénticamente distribuida (iid ), la función de verosimilitud
puede ser escrita como:

L(~x; θ1 , θ2 , ..., θm ) = f (x1 ; θ1 , θ2 , ..., θm )f (x2 ; θ1 , θ2 , ..., θm )...f (xn ; θ1 , θ2 , ..., θm ),

Yn
= f (xi ; θ1 , θ2 , ..., θm ) (8.39)
i=1

Section 8.6
Método de máxima verosimilitud

Este método de estimación de parámetros consiste en maximizar la función de verosimilitud

L(~x; θ1 , θ2 , ..., θm ). El valor del parámetro que maximiza la función se denota por θ̂, y se
conoce como estimador máximo verosı́mil. En la práctica, el calculo de la función conjunta
de verosimilitud conduce a cantidades muy grandes debido a la productoria, en ese sentido, se
opta por calcular el logaritmo de la función de verosimilitud que transforma las productorias
en sumatorios; y realizar un proceso de minimización cambiando de signo a la función.
Adicionalmente, se tienen dos casos de estimación: la estimación de parámetros para datos
no agrupados (unbinned likelihood fit) y estimación para datos agrupados (binned likelihood
fit).

8.6.1 Unbinned likelihood fit

En este caso tenemos que minimizar:

X
θ̂ = arg min − Ln(L(xi ; θ)) (8.40)
θ
xi

donde el ı́ndice i de la suma correo sobre todo el conjunto de datos.

Ejemplo 1:
8.6. Método de máxima verosimilitud 139

Sea A(x1 , x2 , ..., xn ), donde A ∼ P ois(λ) tiene distribución de Poisson con parámetro λ.

n Pn
Y e−λ λxi e−nλ λ i=1 xi
L(~x; λ) = = Qn (8.41)
i=1
xi ! i=1 xi !
Tomando logaritmo natural tenemos:

n
X n
Y
Ln(L(~x; λ)) = −nλ + xi Ln(λ) − Ln( xi !) (8.42)
i=1 i=1

El parámetro máximo verosı́mil es la media muestral.

n
1X
λ̂ = xi = X̄. (8.43)
n i=1
Ejemplo 2:

Sea A(x1 , x2 , ..., xn ), donde A ∼ Exp(α) tiene distribución exponencial con parámetro α.

n  n P
Y 1 −xi /α 1 n
L(~x; α) = e = e− i=1 xi /α (8.44)
i=1
α α
Tomando el logaritmo natural de la función de verosimilitud.

n
X xi
Ln(L(~x; α)) = nLn(1/α) − . (8.45)
i=1
α

Al maximizar Ln(L(~x; α)), el parámetro máximo verosı́mil es la media muestral.

n
1X
α̂ = xi = X̄ (8.46)
n i=1

8.6.2 Binned likelihood fit

En el caso de datos agrupados la función de verosimilitud debe ser evaluada en las marcas de
clase de la distribución de frecuencia. Para tal propósito, se debe construir el histograma que
agrupa a los datos y se usará el ı́ndice j para aclarar que la suma no es sobre el conjunto de
datos (~x) si no sobre cada clase (x̄) del histograma. Si tenemos m clases:

x̄ = {x̄0 , x̄1 , ..., x̄m }. (8.47)

140 Chapter 8. Estimación de parámetros

El histograma por definición está dado por:

~h = {h0 , h1 , ..., hm }, (8.48)

donde hj es el número de eventos contenidos en la j-ésima clase. El método de estimación

de parámetros está dado por:

X
θ̂ = arg min − hj Ln(L(x̄j ; θ)). (8.49)
θ
j

8.6.3 Varianza de los estimadores

Vamos a expandir la función de verosimilitud alrededor del estimador, lo que significa que
estamos en el punto máximo.

2
 
∂Ln(L(θ)) 1 ∂ Ln(L(θ))
Ln(L(θ)) ∼
 
= Ln(L(θ̂)) +  (θ − θ̂) +  (θ − θ̂)2 (8.50)
∂θ  2 ∂θ 2 
θ=θ̂ θ=θ̂

Dado que estamos en el máximo:

(θ − θ̂)2
Ln(L(θ)) ∼
= Ln(L(θ̂)) − (8.51)
2σ̂ 2
donde hemos definido:
 2 −1
2 ∂ Ln(L(θ))
σ̂ = − (8.52)
∂θ2
Note que si evaluamos la función de verosimilitud en θ = θ̂ ± σ̂ obtenemos

Ln(L(θ̂ ± σ̂)) = Ln(L(θ̂)) − 1/2 (8.53)

Esto significa que a 1σ de desviación estándar el logaritmo de la función de verosimilitud

evaluada desde el máximo debe cambiar en 1/2.

Section 8.7
Ejercicios: Parameter Estimation

1. (Theoretical) Sea A(x1 , x2 , ..., xn ), donde A ∼ N (µ, σ) con parámetros µ y σ. Muestre

que los estimadores máximo verosı́miles son:
8.7. Ejercicios: Parameter Estimation 141

n
1X
µ̂ = xi
n i=1
n
1X
σ̂ 2 = (xi − x̄)2 (8.54)
n i=1

2. Para el ejercicio anterior, muestre que la matriz Hessiana evaluada en los estimadores es:

!
−n(σ̂ 2 )−1 0
H(µ̂, σ̂ 2 ) =
0 − n2 (σ̂ 2 )2

¿Es esta matriz definida positiva o negativa?

3. Implemente el algoritmo de Metrópolis Hastings para hacer el ajuste de un histograma

cuyas frecuencias relativas siguen una distribución normal.

a) Descargue los datos de: Gaussian Data

b) Utilice la siguiente distribución a priori uniforme:
(
1 si 3 ≤ µ ≤ 5, 0.5 ≤ σ ≤ 3.5
π(µ, σ) = (8.55)
0 si Otro caso

c) Escriba la función de Likelihood Gaussiana.

N
Y 1 2 /2σ 2
L(x/µ, σ) = √ e−(µ−xi ) (8.56)
i=1 2πσ 2

d) Calcule el logaritmo de la distribución posterior.

ln(P(µ, σ/x)) ∼ ln(L(x/µ, σ)π(µ, σ)) (8.57)

Dado que la distribución encontrada no está normalizada, la relación que existe es de

proporcionalidad.
e) Use el algoritmo de Metrópolis para realizar el muestreo de ln(P(µ, σ/x)) con N =
2 × 104 eventos.
f) Estime el mejor valor de los parámetros del modelo (µ̂, σ̂).
142 Chapter 8. Estimación de parámetros

+
g) Encuentre los errores σ− de los parámetros en un intervalo de confianza del CL = 68%.

4. Resuelva el problema anterior usando el paquete Emcee.

[Link]

5. En general la varianza de estimadores es no calculable:

V (θ̂) = E(θ̂2 ) − E(θ̂)2 (8.58)

En el caso de la distribución exponencial tenemos un valor analı́tico dado por:

Z ∞ Z ∞
 X n 2
1 1 −x1 /θ 1 −xn /θ
V (θ̂) = ... xi e ... e dx1 ...dxn
0 0 n i=1 θ θ
Z ∞ Z ∞ X n  2
1 1 −x1 /θ 1 −xn /θ
− ... xi e ... e dx1 ...dxn
0 0 n i=1 θ θ
θ2
= (8.59)
n
a) Intente encontrar este resultado analı́ticamente.
b) Con el método de MonteCarlo compruebe este resultado para un conjunto de n=20
variables aleatorias xi ...xn ∼ Exp(θ = 2). Generar varias muestras de distribuciones
exponenciales para tener un buen promedio en el ensamble, por ejemplo: N = 106 (Se
obtiene algo como V ar(θ̂) = 0.199).

6. En fı́sica de altas energı́as el siguiente Toy Model es relevante. Se tiene un número

observado de eventos totales (n = 10), un número de eventos de fı́sica conocida b = 9 y
un número de eventos esperados de nueva fı́sica s = 4. La distribución está caracterizada
por la siguiente función de Likelihood:

L(µ, ) = P oisson(n; µs + b)Gauss(; 1., 0.1) (8.60)

1 −(µs+b) 1 (−1)2
L(µ, ) = e (µs + b)n √ e− 2σ2 (8.61)
n! 2πσ 2
Donde  es la eficiencia de reconstrucción de los eventos de fı́sica conocida tiene valor
medio ¯ = 1. y desviación estándar σ = 0.1. La contribución Gaussiana modela el error
8.7. Ejercicios: Parameter Estimation 143

sistemático asociado a la estimación de b. µ se conoce como signal strength, el cuál es

muy importante para descartar o no, nuevas teorı́as es fı́sica. Con esta información:

(a) Reproduzca las gráficas de la Figura [8.2]. Los parámetros pueden variar como
0. < µ < 2. y 0. <  < 2.

(b) Usando el algoritmo de Metropolis-Hastings encuentre µ̂ y ˆ y los errores asociados

a cada parámetro.

(c) Elimine el error sistemático dado por la parte gaussiana. ¿Cuál serı́a el valor de µ̂?.
A qué conclusiones llega?

2.00 0.475
1.75 0.425
1.50 0.375
0.325
1.25 0.4
0.275 0.3
1.00
(,)

0.225 0.2
0.75
0.175 0.1
0.50 0.125 0.0 0.00
0.25
0.00
0.25 0.50
0.25 0.075 0.50
0.751.00 0.75
1.00
1.251.50 1.25
1.50
0.00 0.025 1.752.00 1.75
2.00
0.00 0.25 0.50 0.75 1.00 1.25 1.50 1.75 2.00

Figure 8.2: Gráfica de curvas de nivel de la función L(µ, ) (izquierda) y superficie (derecha).

7. Un equipo de medición de 3 canales detecta los pulsos de una posible señal extraterrestre.
En un momento determinado del tiempo se tiene la estadı́stica descrita por la siguiente
Figura [8.3]:

20.0 Data
Background
17.5 Signal
15.0
12.5
10.0
7.5
5.0
2.5
0.0
0.0 0.5 1.0 1.5 2.0 2.5 3.0

Figure 8.3: Medida del equipo de 3 canales.

144 Chapter 8. Estimación de parámetros

El conjunto de datos es: Data=[10,20,5], el background o la señal esperada sin existencia

extraterrestre bkg=[6,15,10] y un modelo de señal extraterrestre Signal=[2,10,3]. La
función de Likelihood que describe esta observación para un canal está dada por:

Li (µ, ) = P oisson(n; µs + b)Gauss(; 1., 0.1), (8.62)

1 −(µs+b) 1 (−1)2
Li (µ, ) = e (µs + b)n √ e− 2σ2 . (8.63)
n! 2πσ 2
Donde  es la eficiencia de reconstrucción de los eventos sin existencia extraterreste y
σ = 0.1 su desviación estándar. µ es la fortaleza de la señal y es una medida para
comprobar o descartar, en este caso, la señal extraterrestre. La función de Likelihood
(the joint probability function) de todos los canales está dada por :

N Channels
Y
L(µ, ) = Li (µ, ). (8.64)
i=1

a) Use el algoritmo de Metropolis-Hastings para encontrar los mejores estimadores que

describen esta observación µ̂ y ˆ con sus respectivos errores (asuma que todos los
canales tiene la misma resolución experimental (σ)). Los parámetros pueden variar
para 0. ≤ µ ≤ 3. y para 0.5 ≤  ≤ 1.5. Tenga cuidado cuando la caminata aleatoria
tome valores negativos para µ /.
b) Si µ es cercano a cero indica que los datos están bien descritos por la hipótesis de
no existencia extraterrestre. Por el contrario, si µ se acerca a uno podrı́a indicar que
existe vida extraterrestre. ¿Qué conclusiones reportarı́a usted sobre la observación
realizada?. Nota: En realidad se debe hacer una prueba de hipótesis pero se pierde la
alineación entre lo visto en clase y la evaluación.
Chapter 9

Confidence Intervals

Section 9.1
Intervalos de confianza para la media

Sean x1 , x2 , ..., xn variables aleatorias independiente con E(xi ) = µ y V ar(xi ) < ∞, entonces:

x̄ − µ
zn := √ (9.1)
σ/ n
la función de distribución de zn , converge a una función de distribución normal estándar
z ∼ N (0, 1).
El intervalo de confianza para la media con varianza conocida está dada por:

 
σ σ
IC1−α = x̄ − zα/2 √ < µ < x̄ + z1−α/2 √ (9.2)
n n
Este intervalo es válido para población aproximada normal, u otra distribución poblacional
con n >= 30.

Cuando no se conoce la varianza poblacional y la muestra es pequeña, se necesita hacer

una estimación del intervalo usando la distribución T-Student. El intervalo de confianza para
la media está dado por:

 
S S
IC1−α = x̄ − tα/2,n−1 √ < µ < x̄ + t1−α/2,n−1 √ (9.3)
n n
donde S es la varianza muestral.

145
146 Chapter 9. Confidence Intervals

Section 9.2
Intervalo de confianza para la diferencia de medias

Para encontrar un intervalo de confianza para la diferencia de medias (n ≤ 30), se tiene la

siguiente expresión:

s s
s21 s22 s21 s2
 
IC1−α = (x̄1 − x̄2 ) − tα/2,ν + < µ1 − µ2 < (x̄1 − x̄2 ) + t1−α/2,ν + 2 (9.4)
n1 n2 n1 n2
Los grados de libertad se obtienen usando:

s2 s22 2
( n11 + n2
)
ν= s2 s2
(9.5)
( n1 )2 ( n2 )2
1
n1 −1
+ 2
n2 −1

Para encontrar un intervalo de confianza para la diferencia de medias (n > 30), se tiene la
siguiente expresión:

s s
σ12 σ22 σ12 σ22
 
IC1−α = (x̄1 − x̄2 ) − zα/2 + < µ1 − µ2 < (x̄1 − x̄2 ) + zα/2 + (9.6)
n1 n2 n1 n2

Section 9.3
Intervalo de confianza para una proporción

Para estimar el intervalo de confianza de un proporción, se tiene:

 r r 
p̂(1 − p̂) p̂(1 − p̂)
IC1−α = p̂ − zα/2 < P < p̂ + z1−α/2 (9.7)
n n
donde se considera la aproximación binomial a la normal si np̂ ≥ 10 y n(1 − p̂) ≥ 10.

Section 9.4
Intervalo para la varianza

El estimador de la varianza muestral está dada por:

n
2 1 X
s = (xi − x̄)2 (9.8)
n − 1 i=1
9.5. Ejercicios: Intervalos de confianza 147

Este estadı́stico es un estimador insesgado de la varianza poblacional puesto que:

E(s2 ) = σ 2 (9.9)

Con base a los criterios de la ley de los grandes números y suponiendo que las muestras
siguen una distribución normal, X ∼ N (µ, σ 2 ). La varianza muestral como variable aleatoria
sigue una función de distribución χ2n−1 (s) de n − 1 grados de libertad.

s2
(n − 1) ∼ χ2n−1 (9.10)
σ2
En problemas donde la varianza poblacional no es conocida, la distribución permite
encontrar intervalos de confianza para la varianza poblacional; dada la varianza muestral.

Section 9.5
Ejercicios: Intervalos de confianza

1. El contenido de botellas de un litro tienen una distribución aproximadamente normal con

varianza de 0.15 litros2 . Se toma una muestra aleatoria de cajas y se mide el contenido,
obteniendo el siguiente resultado.

S = [0.974, 0.950, 0.932, 1.104, 1.038, 0.920, 0.935, 0.907, 0.810, 0.915] (9.11)

a) Construya el intervalo de confianza del 95% para la media población.

 
IC95% = 0.708 < µ < 1.188 (9.12)

b) Construya el intervalo de confianza del 95% usando el método de Bootstrapping (N =

50000 toy experiments).

 
IC95% = 0.903 < µ < 1.00 (9.13)

Notar que los IC no son consistentes, esto puede indicar que la población que describe
a los datos no es normal. Usted puede calcular el IC en el caso de varianza desconocida
y comprobar que se obtiene un mejor resultado.
148 Chapter 9. Confidence Intervals

2. Construir un IC de 90% de nivel de confianza para la estatura media de un grupo de

estudiantes. Se sabe que la estatura (en cm) es aproximadamente normal y está dada
por:

S = [175, 176, 180, 164, 170, 170, 181, 169, 163, 190, 170, 171] (9.14)

 
IC90% = 169.28 < µ < 177.21 (9.15)

3. Se tienen dos muestras de la medida de la resistencia eléctrica de un lote de dispositivo

electrónicos, las cuales están dadas por:

S1 = [6.0, 5.0, 6.5, 5.0, 4.0, 5.0, 5.0, 5.0, 7.0, 5.5, 4.5]
S2 = [7.0, 8.0, 8.5, 7.4, 8.9, 6.7, 9.0, 8.4, 7.8, 5.3, 8.1] (9.16)

Determine el intervalo al 95% de nivel de confianza. ¿Pertenecen las muestras a la misma

población?

 
IC95% = − 3.16 < µ1 − µ2 < −1.60 (9.17)

Las muestras no pertenecen a la misma población, presenta una desviación a ∼ 4σ.

4. Una moneda equilibrada es lanzada 500 veces. Sea X la variable aleatoria asociada al
número de soles que aparecen en los n lanzamientos. Usando la distribución Binomial y
la aproximación a la normal, estime:

a) P(X ≤ 260.) = 0.8261 Binomial, 0.8144 Normal.

b) P(230 ≤ X ≤ 260) = 0.926 Binomial, 0.926 Normal.

5. Se estudia la intención de voto de un candidato A a la presidencia de un paı́s. Para ello,

se toma una muestra aleatoria de 800 habitantes y 325 de ellos responden a la intención
de votar por A. Encontrar el intervalo de confianza a 98%, para la proporción de votantes
con intensión de voto por el candidato A.

 
IC98% = 0.365 < P < 0.446 (9.18)
9.5. Ejercicios: Intervalos de confianza 149

6. Una empresa de equipo médico quiere saber la proporción de clientes que prefieren su
marca. El gerente afirma que el 30% de los clientes del mercado prefieren sus productos.
Una muestra de 100 clientes indicó que 24 usan dicha marca. Calcule el intervalo de
confianza al 90% e indicar si el gerente tiene la razón.

 
IC90% = 0.169 < P < 0.310 (9.19)

El gerente está en lo correcto dentro de la confianza de la hipótesis nula.

7. Demuestre que E(s2 ) = σ 2 . Hint: Recuerde que V ar(x) = σ 2 .

8. Genere una muestra de números x ∼ N (0, 9) de tamaño 100.

a) Escriba las funciones para calcular la varianza muestral y χ2n−1 .

b) Con el método de re-muestreo (Bootstrapping) genere (N = 50000 toy experiments)
muestras con reemplazo para encontrar la distribución de la varianza.
c) Dibuje el histograma y la función de distribución de probabilidad χ2n−1 , usando
[Link](x_,df=n-1).
d) Calcule el valor medio de la distribución usando el histograma y usando
[Link](df=df).
e) Calcule el valor medio de la varianza poblacional (deberı́a ser cercano a 9).
f) Calcule el intervalo de confianza de este estimador a un nivel de confianza del 95%
(significancia α = 0.05).

(n − 1)S 2 (n − 1)S 2
 
2
IC95% = < σ < (9.20)
χ21−α/2 χ2α/2

9. A un grupo de individuos en un hospital se sometió a un medicamento experimental, al

final se midió el nivel de glóbulos blancos en sangre. En una escala especı́fica los resultados
fueron los siguientes:

S = [6.0, 6.4, 7.0, 5.8, 6.0, 5.8, 5.9, 6.7, 6.1, 6.5, 6.3, 5.8] (9.21)

Suponiendo que la población es normalmente distribuida, determine un intervalo de

confianza al 95% para la varianza poblacional del nivel de colesterol.
150 Chapter 9. Confidence Intervals

 
2
IC95% = 0.0770 < σ < 0.4428 (9.22)

10. Un fabricante de parte de automóvil garantiza que sus baterı́as duran en promedio 3 años
con una desviación estándar de 1 año, si 5 baterı́as presentan una varianza de 0.815.
Realice su estudio al 99% de nivel de confianza.

a) ¿Está la observación de acuerdo con la afirmación del fabricante de que la desviación

estándar es de 1 año?
b) ¿Cuál es el p-value de la observación? p = 0.620.
c) ¿Cuál es valor lı́mite (crı́tico) que puede tener la varianza de las baterı́as para mantener
2
la afirmación del fabricante? Sup = 3.715.

11. Encuentre la probabilidad de que una muestra aleatoria de 30 observaciones de una

población normal con varianza σ 2 = 10. Tenga una varianza muestral:

(a) P(S 2 ≥ 9.) = 0.379

(b) P(3.5 ≤ S 2 ≤ 11.) = 0.675
(c) Calcule el p-value de observar S 2 = 15, p = 0.04
Chapter 10

Pruebas de Hipótesis

El procedimiento estadı́stico para determinar si los datos de una muestra son compatibles con
las caracterı́sticas de una población. Para contextualizar el problema se realizan las siguientes
definiciones:

• Una hipótesis estadı́stica es una proposición sobre los parámetros de una población.

• Se denomina hipótesis nula a la hipótesis que se desea contrastar, la denominamos H0 .

La hipótesis alternativa la denominamos con H1 .

• Un estadı́stico de prueba es un estadı́stico utilizado para determinar si se rechaza o no la

hipótesis nula H0 .

• Región de no rechazo es el conjunto de valores del estadı́stico de prueba, para los cuales
no rechazamos la hipótesis nula. Su región complementaria se llama la región de rechazo
de H0 en favor de H1 .

• Dado que hay una probabilidad de que el estadı́stico genere un valor en la región de
rechazo, se tienen dos tipos de errores: Error tipo-1 es rechazar la hipótesis nula cuando
es verdadera y el Error tipo-2 es no rechazar la hipótesis nula cuando es falsa.

• El nivel de significación del test determina la probabilidad de que el estadı́stico genere

un valor en la región de rechazo de H0 . Se denomina α y está restringido al intervalo
0 ≤ α ≤ 1.

151
152 Chapter 10. Pruebas de Hipótesis

Section 10.1
Ejercicios: Tablas de contingencia

1. En cierta empresa, se calcula que la cantidad promedio de artı́culos vendidos por semana
es 2000. Se quiere averiguar si hay una relación entre el dı́a de la semana y la cantidad
de artı́culos que se venden. La Tabla [10.1] muestra las observaciones.

< 1800 1800-2200 > 2200

Lunes-Miércoles 38 55 52
Jueves-Viernes 39 65 55
Sábado-Domingo 43 60 63

Table 10.1: Tabla de contingencia entre artı́culos vendidos y los dı́as de la semana

a) Calcule el número total de observaciones.

b) Calcule los vectores de frecuencias marginales.
c) Calcule los vectores de probabilidad marginales.
d) Calcule la matriz de probabilidad.
e) Calcule la matriz de frecuencias esperadas.
f) Calcule la función de costo: χ2 dada por:

n
2
X (fe − fo )2
χ = (10.1)
i=1
fe

donde fe es la frecuencia esperada y fo es la frecuencia observada de cada punto de la

tabla.
g) Calcule el número de grados de libertad (ν).
h) Dibuje la distribución χ2ν con el número correcto de grados de libertad.
i) Describa la hipótesis nula y la alternativa.
j) Calcule el p-value observado.
k) ¿Cuál es valor crı́tico de la distribución para un nivel de confianza del 95%?
l) Indique si se rechaza o no, la hipótesis nula.
10.1. Ejercicios: Tablas de contingencia 153

2. La universidad de los Andes realiza una encuesta a N = 100 estudiantes sobre la intención
de voto de tres candidatos a la presidencia, los candidatos son: A, B y C. Adicionalmente,
se tuvo en cuenta el género de cada estudiante. Los resultados se pueden resumir en la
siguiente Tabla:

A B C
Mujeres 11 15 34
Hombre 23 6 11

Table 10.2: Tabla de contingencia asociada a la encuesta realizada por la universidad.

Los resultados mostrados en la Tabla [10.2] representan la frecuencias observadas (fo )

de la intención de voto de cada candidato. La universidad desea saber si las variables
género y candidato son o no independientes con un nivel de confianza del 95%. Para tal
propósito, planteamos las siguientes hipótesis:

a) H0 : Las variables género y candidato son independientes. (Hipótesis nula)

b) H1 : Las variables género y candidato no son independientes. (Hipótesis alternativa)

Note que la hipótesis alternativa es la negación de la hipótesis nula H1 = ¬H0 .

Suponemos que la hipótesis nula es verdadera para calcular las frecuencias
esperadas (fe ). Entonces:

a) Si las variables son independientes entonces P(A ∩ B) = P(A) · P(B). Entonces, deberá
calcular las probabilidades marginales y luego calcular el valor de la probabilidad de
cada entrada en la tabla. Usando el producto de Kronecker serı́a muy eficiente y claro.

b) La frecuencia esperada en cada posición es: fei = P(A ∩ B)i · N , i.e, la probabilidad de
cada evento multiplicado por el tamaño de muestra.

La matriz de probabilidad P(x, y) = G(x) ⊗ H(y), puede ser escrita como el producto de
Kronecker de los vectores de probabilidad marginal (serı́a interesante que implementaran
esta operación en su código). Finalmente, para medir la discrepancia entre las frecuencias
esperadas y observadas se tiene el estadı́stico de prueba (muy conocido) χ2 dado por:
154 Chapter 10. Pruebas de Hipótesis

k
2
X (f i − f i )2
o e
χ = (10.2)
i=1
fei

Este valor se compará con la distribución χ2α,(r−1)(c−1) . En esta definición, α = 0.05 se

refiere al nivel de significación de la prueba, r es el número de filas y c el número de
columnas. (r − 1)(c − 1) son los grados de libertad (conceptos por definir) asociados a
nuestra prueba estadı́stica; entonces para nuestro ejemplo, χ20.05,2 ≈ 5.991 nos definirá el
valor crı́tico del estadı́stico a un nivel de confianza del 95%.

En Python:

from [Link] import chi2

[Link](0.95,df=2)

ppf devuelve el percentil asociado a un (p-value) p = 0.05, para una distribución χ2 de

dos grados de libertad. Calcule el χ2 ≈ 16.50, que es un valor más extremo que el valor
crı́tico, por lo cuál, debemos rechazar la hipótesis nula con un nivel de confianza del 95%.
Las variables son dependientes estadı́sticamente.

Remark:

Efectivamente, como muestra la tabla, las mujeres prefieren el candidato C y los hombres
el candidato A.

Section 10.2
Prueba de hipótesis para generadores de números aleatorios

¿Sabı́as que el generador de números de Numpy (Mersenne - Twister 1998), tiene uniformidad en
d=623 dimensiones al 95% de nivel de confianza? Para probar la uniformidad de los números
aleatorios generados por pseudo-generadores es posible realizar una prueba de hipótesis. Por
lo cuál, se tienen las siguientes hipótesis:

a) H0 : La secuencia de números sigue una distribución uniforme en k-dimensiones.

10.3. Ejercicios: Hypotesis testing 155

b) H1 : La secuencia de números no sigue una distribución uniforme en k-dimensiones.

Section 10.3
Ejercicios: Hypotesis testing

1. (Binomial) Se realiza un experimento aleatorio del lanzamiento de 5 monedas N =

1000, cuyo espacio muestral está definido por Ω = {C, S}. Se mide el número
x de caras que se observan durante el experimento. El resultado del experimento
se resume en: [Link]
MetodosComputacionalesReforma/[Link]. Note que es un archivo para
cargar en Pandas. Se desea saber si el número de caras tiene una distribución binomial
x ∼ Bin(n = 5, p). Se tienen las siguientes hipótesis:

a) H0 : El número de caras tiene distribución binomial x ∼ Bin(n = 5, p = p̂).

b) H1 : El número de caras tiene distribución binomial.

Suponemos que la hipótesis nula es verdadera, de este modo, necesitamos realizar los
siguientes pasos para realizar la bondad del ajuste:

a) Estimar el parámetro de la distribución p̂ = 0.4946

b) Suponiendo que H0 es verdadera, estimar la probabilidad de que se presenten el número
de caras encontrados en el experimento. Ejemplo: la probabilidad de obtener 0 caras
en un experimento es:
 
5
P(X = 0|p̂ = 0.4946) = 0.49460 (1 − 0.4946)5−0 = 0.0329 (10.3)
0
c) Estimar las frecuencias esperadas asociadas a el número de errores fei (x) = Pi (x) × N .
d) Dibuje los histogramas correspondientes a las distribución esperada y observada.
e) Medir la discrepacia que existe entre las frecuencias esperadas y observadas a través
del estadı́stico de prueba χ2 .

k
2
X (f i − f i )2
o e
χ = (10.4)
i=1
fei

f) Cuál es el p-value de esta observación.

156 Chapter 10. Pruebas de Hipótesis

g) Calcular el valor crı́tico de la distribución χ2α,k−t−1 , donde la significancia es α = 0.05.

k = 10 es el número de categorı́as, t es el número de parámetros que se estimaron en
el proceso, en este caso, se estimo λ̂. Note que los grados de libertad son gl = 4.
h) Si el estadı́stico de prueba es menor que el valor crı́tico no rechazamos la hipótesis
nula.
i) Dar la conclusión.

2. (Poisson) Una empresa de motos realiza la construcción de N = 440 motos. Durante las
pruebas mecánicas y eléctricas, la empresa reportó el número de fallas que presenta cada
artı́culo. El archivo que resume esta información es: [Link]
com/asegura4488/Database/main/MetodosComputacionalesReforma/[Link].
Note que es un archivo para cargar en Pandas. La empresa quiere saber si la distribución
de errores es una distribución de Possion y cuál es el estimador máximo verosimil que
describe dicha distribución. De esta manera, se tienen las siguientes hipótesis.

a) H0 : El número de errores tiene distribución Poisson n ∼ P oiss(λ = λ̂).

b) H1 : El número de errores no tiene distribución Poisson.

Suponemos que la hipótesis nula es verdadera, de este modo, necesitamos realizar los
siguientes pasos para realizar la bondad del ajuste:

a) Estimar el parámetro de la distribución λ̂ = 3.047

b) Suponiendo que H0 es verdadera, estimar la probabilidad de que se presenten el número
de errores encontrados en las motos. Ejemplo: la probabilidad de obtener 0 errores en
una moto es:

e−3.047 3.0470
P(x = 0|λ̂ = 3.047) = = 0.04746 (10.5)
0!
c) Estimar las frecuencias esperadas asociadas a el número de errores fei (x) = Pi (x) × N .
d) Dibuje los histogramas correspondientes a las distribución esperada y observada.
e) Medir la discrepacia que existe entre las frecuencias esperadas y observadas a través
del estadı́stico de prueba χ2 .

k
2
X (f i − f i )2
o e
χ = (10.6)
i=1
fei
10.3. Ejercicios: Hypotesis testing 157

f) Cuál es el p-value de esta observación.

g) Calcular el valor crı́tico de la distribución χ2α,k−t−1 , donde la significancia es α = 0.05.
k = 10 es el número de categorı́as, t es el número de parámetros que se estimaron en
el proceso, en este caso, se estimo λ̂. Note que los grados de libertad son gl = 8.
h) Si el estadı́stico de prueba es menor que el valor crı́tico no rechazamos la hipótesis
nula.
i) Dar la conclusión.

3. Una partı́cula se mueve aleatoriamente desde el centro de una caja cuadradada de longitud
l = 5. El ángulo para realizar el paso se genera uniformemente entre 0 y 2π y luego se
actualiza la posición actual:

x += lstep * [Link](theta)
y += lstep * [Link](theta)

lstep es el tamaño del paso, use lstep = 0.4 para este experimento. Se quiere medir el
número de pasos promedio que le toma a una partı́cula salir de la caja. Genere N = 200
experimentos guardando el número de pasos que le toma a la partı́cula salir de la caja.
Una posible trayectoria se describe en la Figura [10.1]

a) Eliga una distribución de clases que tenga su valor mı́nimo y máximo en los valores
mı́nimo y máximo de la lista y un ancho de paso h = 15 steps:
bins = [Link](min_,max_+h,h)

b) Calcule el histograma de la distribución.

c) Normalice el histograma a la unidad.
H1Norm = H1 / [Link](H1*w)

d) Realice el ajuste de la distribución usando la densidad de Weibull exponencial:

from [Link] import exponweib

a,c,d,e = [Link](lista)
158 Chapter 10. Pruebas de Hipótesis

Y 3

0
0 1 2 3 4 5
X
Figure 10.1: Trayectoria de la partı́cula (Toy model). En este experimento la partı́cula le tomo
33 pasos salir de la caja.

lbox=5
Exponential Weibull a= 1.47, b=1.10, n=51.57
0.016 Observed

0.014

0.012

0.010

0.008

0.006

0.004

0.002

0.000
25 50 75 100 125 150 175 200
Nsteps

Figure 10.2: Ajuste a los datos usando la distribución de Weibull exponencial.

deberı́a obtener un resultado similar a Figura [10.2]:

e) Calcule el valor medio de pasos usando los métodos de scipy: mean = [Link](a,c,d,e)

f) Ahora se quiere realizar la bondad en el ajuste. Genere la distribución esperada usando

el método de generación de Scipy, con mismo número de eventos observados N = 200.
[Link]([Link](a,c,d,e))
10.3. Ejercicios: Hypotesis testing 159

Deberı́a obtener distribuciones similares como se muestra en la Figura [10.3].

lbox=5
Expected
Observed
50

40

30

20

10

0
25 50 75 100 125 150 175 200
Nsteps

Figure 10.3: Comparación entre la distribución observada y esperada para el número de pasos
de escape de la partı́cula.

g) Calcule el estadı́stico de prueba χ2 : (Ignore las clases que tengan menos de 5 entradas,
es decir las colas de la distribución).

k
X (f i − f i )2
o e
χ2 = (10.7)
i=1
fei

h) Calcular el valor crı́tico de la distribución χ2α,k−t−1 , donde la significancia es α = 0.05. k

es el número de categorı́as, t es el número de parámetros que se estimaron en el proceso,
en este caso, se estimo n̂. Note que los grados de libertad son gl = #Clases − 1 − 1.

i) Si el estadı́stico de prueba es menor que el valor crı́tico no rechazamos la hipótesis

nula.

j) Dar la conclusión.

4. En el caso 2D, se divide el espacio entre [0, 1] en k d celdas de igual área, luego se genera
una muestra aleatoria de tamaño 2n (note que el número de eventos es par), podemos
construir n pares no solapados (x0 , x1), (x2 , x3 ), ..., (x2n−1 , x2n ) y estimar la frecuencia
observada de puntos en cada celda. Si suponemos que la hipótesis nula es verdadera, la
frecuencia esperada de puntos en cada celda es:
160 Chapter 10. Pruebas de Hipótesis

n
fei = (10.8)
kd

El estadı́stico de prueba que mide la discrepancia entre la frecuencia esperada y observada

está dado por el estadı́stico χ2 :

k
2
X (f i − f i )2o e
χ = (10.9)
i=1
fei

Calcule el valor crı́tico de la distribución χ2α,kd −1 , donde la significancia es α = 0.05. k es

el número de celdas en las que se divide el dominio. Note que los grados de libertad son
gl = k d − 1.

(a) Genere una muestra de 2 ∗ (1000) números aleatorios con el generador de Numpy.

(b) Use una partición de k = 10 celdas.

(c) Calcule las frecuencias observadas y estadı́stico de prueba. Guardar los conteos en
objetos M = [Link]((k,k)). En este caso, es posible usar histogramas 2D para
contar los puntos en cada celda.

(d) Calcule el p-value de la observación.

(e) ¿Valida o rechaza la hipótesis al 95% de nivel de confianza de que los números tienen
una distribución uniforme en d = 2 dimensiones?

0.030
pmf 99
2 (x)
Upper Limit: 123.23
0.025 Observed: 116.40

0.020

0.015

0.010

0.005

0.000
60 80 100 120 140 160
2

Figure 10.4: Distribución χ299 , el valor observado es menor que el valor crı́tico de la distribución.
p − value = 0.111, se acepta H0 .
10.3. Ejercicios: Hypotesis testing 161

5. En el caso 3D, se divide el espacio entre [0, 1] en k d celdas de igual volumen, luego se
genera una muestra aleatoria de tamaño 3n (note que el número de eventos es multiplo de
3), podemos construir n pares no solapados (x0 , x1 , x2 ), (x3 , x4 , x5 ), ..., (x3n−2 , x3n−1 , x3n )
y estimar la frecuencia observada de puntos en cada celda. Si suponemos que la hipótesis
nula es verdadera, la frecuencia esperada de puntos en cada celda es:

n
fei = (10.10)
kd

El estadı́stico de prueba que mide la discrepancia entre la frecuencia esperada y observada

está dado por el estadı́stico χ2 :

k
2
X (f i − f i )2
o e
χ = (10.11)
i=1
fei

Calcule el valor crı́tico de la distribución χ2α,kd −1 , donde la significancia es α = 0.05. k es

el número de celdas en las que se divide el dominio. Note que los grados de libertad son
gl = k d − 1.

(a) Genere una muestra de 3 ∗ (1000) números aleatorios con el generador de Numpy.
(b) Use una partición de k = 10 celdas.
(c) Calcule las frecuencias observadas y estadı́stico de prueba. Guardar los conteos en
objetos M = [Link]((k,k,k)).
(d) Calcule el p-value de la observación.
(e) ¿Valida o rechaza la hipótesis al 95% de nivel de confianza de que los números tienen
una distribución uniforme en d = 3 dimensiones?
162 Chapter 10. Pruebas de Hipótesis

pmf 999
2 (x)

0.008 Upper Limit: 1073.64

Observed: 1038.00

0.006

0.004

0.002

0.000
850 900 950 1000 1050 1100 1150
2

Figure 10.5: Distribución χ2999 , el valor observado es menor que el valor crı́tico de la distribución.
p − value = 0.190, se acepta H0 .
References

[1] Roland Combescot. Superconductivity: An Introduction. Cambridge university press, 2022.

[2] William Ogilvy Kermack and Anderson G McKendrick. A contribution to the mathematical
theory of epidemics. Proceedings of the royal society of london. Series A, Containing
papers of a mathematical and physical character, 115(772):700–721, 1927, https://
[Link]/doi/pdf/10.1098/rspa.1927.0118.

[3] Rubin H Landau, José Páez, Manuel José Páez Mejı́a, and Cristian C Bordeianu. A survey
of computational physics: introductory computational science. Princeton University Press,
2008, [Link]

[4] Sheldon M Ross. Introduction to probability and statistics for engineers and scientists.
Academic press, 2020.

[5] Ronald E Walpole, Raymond H Myers, Sharon L Myers, and Keying Ye. Probabilidad y
estadı́stica para ingenierı́a y ciencias. 4 Edición, 2007.

163