Fundamentos teóricos de los métodos Monte Carlo, aplicados en el software PHOEBE para encontrar los modelos más probables para un conjunto de datos de curvas de luz.
métodos MonteCarlo Curvas de luz en estrellas binarias y efectos en los tiempos de tránsitos FaMAF - OAC Marcelo Lares, clase invitada, 2 de octubre de 2020 Curvas de luz en estrellas binarias y efectos en los tiempos de tránsitos Inferencia aproximada con métodos MonteCarlo 1 / 87
Carlo methods Momentos de una distribución Integral de Riemann–Stieltjes Estadística inferencial y estimadores 2 Model selection Cuadrados mínimos Regresión con modelos no lineales 3 Muestreo (Sampling) Importance sampling MonteCarlo Sampling Sampling algorithms Curvas de luz en estrellas binarias y efectos en los tiempos de tránsitos Inferencia aproximada con métodos MonteCarlo 2 / 87
quiere decir "Monte Carlo Markov Chains" MCMC NO ES: Un método de optimización Un método para encontrar el mínimo de ξ2 Un método para encontrar el máximo de Likelihood Un método para encontrar el mejor modelo Un método de ajuste de modelos "Monte Carlo" viene de hacer cuentas con números randoms... Curvas de luz en estrellas binarias y efectos en los tiempos de tránsitos Inferencia aproximada con métodos MonteCarlo 4 / 87
distribución Integral de Riemann–Stieltjes Estadística inferencial y estimadores Aplicaciones de los PRNs Una de los principales usos de los pseudo random numbers (PRNs) se encuentra en los métodos de Monte Carlo: Definition Los métodos Monte Carlo son un conjunto de algoritmos computacionales que se basan en el muestreo aleatorio para estimar resultados numéricos. e.g.: Integración Optimización Simulación Curvas de luz en estrellas binarias y efectos en los tiempos de tránsitos Inferencia aproximada con métodos MonteCarlo 5 / 87
distribución Integral de Riemann–Stieltjes Estadística inferencial y estimadores Aplicaciones de PRNGs: Integración MonteCarlo La estimación de integrales es una de los principales usos de los PRNs. Curvas de luz en estrellas binarias y efectos en los tiempos de tránsitos Inferencia aproximada con métodos MonteCarlo 6 / 87
distribución Integral de Riemann–Stieltjes Estadística inferencial y estimadores Aplicaciones de PRNGs: Integración MonteCarlo Estimación de π: Curvas de luz en estrellas binarias y efectos en los tiempos de tránsitos Inferencia aproximada con métodos MonteCarlo 7 / 87
distribución Integral de Riemann–Stieltjes Estadística inferencial y estimadores Aplicaciones de PRNGs: Integración MonteCarlo Las probabilidades de ocurrencia de una variable aleatoria contínua en un dado intervalo se calculan como una integral. La misma se puede hacer con métodos Monte Carlo. Ésta es la idea de porqué funciona: Curvas de luz en estrellas binarias y efectos en los tiempos de tránsitos Inferencia aproximada con métodos MonteCarlo 8 / 87
distribución Integral de Riemann–Stieltjes Estadística inferencial y estimadores Definición de variable aleatoria (VA) Curvas de luz en estrellas binarias y efectos en los tiempos de tránsitos Inferencia aproximada con métodos MonteCarlo 9 / 87
distribución Integral de Riemann–Stieltjes Estadística inferencial y estimadores Función distribución de una VA La función de distribución de una VARIABLE ALEATORIA X es la función F(x) : R → [0, 1] dada por F(x) = P(X ≤ x) A veces se la denota por FX (x) y también se la llama "función acumulada". Esta función contiene toda la información sobre X y sobre P. Curvas de luz en estrellas binarias y efectos en los tiempos de tránsitos Inferencia aproximada con métodos MonteCarlo 10 / 87
distribución Integral de Riemann–Stieltjes Estadística inferencial y estimadores Propiedades matemáticas de F Las propiedades de F se deducen de las propiedades de las medidas de probabilidad: F(x) = P(X ≤ x) limx→∞ F(x) = 1 limx→−∞ F(x) = 0 x1 ≤ x2 =⇒ F(x1 ) ≤ F(x2 ) lim x→x+ 0 F(x) = F(x0 ) Curvas de luz en estrellas binarias y efectos en los tiempos de tránsitos Inferencia aproximada con métodos MonteCarlo 11 / 87
distribución Integral de Riemann–Stieltjes Estadística inferencial y estimadores Propiedades matemáticas de F Sean: F(x+) = lim x→0+ F(x + y) F(x−) = lim x→0− F(x + y) Entonces: P(X < a) = F(a−) P(X = a) = F(a) − F(a−) P(a < X ≤ b) = F(b) − F(a) P(a ≤ X ≤ b) = F(b) − F(a−) P(a < X < b) = F(b−) − F(a) P(a ≤ X < b) = F(b−) − F(a−) Curvas de luz en estrellas binarias y efectos en los tiempos de tránsitos Inferencia aproximada con métodos MonteCarlo 12 / 87
distribución Integral de Riemann–Stieltjes Estadística inferencial y estimadores Discontinuidades de F Curvas de luz en estrellas binarias y efectos en los tiempos de tránsitos Inferencia aproximada con métodos MonteCarlo 13 / 87
distribución Integral de Riemann–Stieltjes Estadística inferencial y estimadores Propiedades matemáticas de F En particular, si F es una función contínua: F(b) − F(a) = P(a < X < b) (1) = P(a ≤ X < b) (2) = P(a < X ≤ b) (3) = P(a ≤ X ≤ b) (4) Curvas de luz en estrellas binarias y efectos en los tiempos de tránsitos Inferencia aproximada con métodos MonteCarlo 14 / 87
distribución Integral de Riemann–Stieltjes Estadística inferencial y estimadores VARIABLE ALEATORIA es absolutamente contínua Una VARIABLE ALEATORIA es absolutamente contínua si ∃ f ≥ 0 e integrable tal que ∀ x se cumple que F(x) = x −∞ f(u) du En tal caso f se llama función de densidad de probabilidad. Si X es absolutamente contínua, por el teorema fundamental del cálculo tenemos que F (x) = f(x) Curvas de luz en estrellas binarias y efectos en los tiempos de tránsitos Inferencia aproximada con métodos MonteCarlo 15 / 87
distribución Integral de Riemann–Stieltjes Estadística inferencial y estimadores VARIABLE ALEATORIA es absolutamente contínua F(x) = x −∞ f(u) du Curvas de luz en estrellas binarias y efectos en los tiempos de tránsitos Inferencia aproximada con métodos MonteCarlo 16 / 87
distribución Integral de Riemann–Stieltjes Estadística inferencial y estimadores La integral como medida de probabilidad Curvas de luz en estrellas binarias y efectos en los tiempos de tránsitos Inferencia aproximada con métodos MonteCarlo 17 / 87
distribución Integral de Riemann–Stieltjes Estadística inferencial y estimadores Propiedades de f fX (x) ≥ 0 ∞ −∞ fX (x) dx = 1 fX (x) es contínua a pedazos P(a ≤ X ≤ b) = b a fX (x) dx Curvas de luz en estrellas binarias y efectos en los tiempos de tránsitos Inferencia aproximada con métodos MonteCarlo 18 / 87
distribución Integral de Riemann–Stieltjes Estadística inferencial y estimadores Integral de Riemann–Stieltjes Es una integral de la forma: b a h(x) dF(x) donde dF es la función integradora, que debe ser contínua por derecha y no decreciente. Además debe satisfacer que: lim x→∞|F(x)−F(−x)|<M Estas integrales tienen propiedades similares a las integrales de Riemann, y dado que las funciones de distribución satisfacen las condiciones de la integral de RS, se pueden usar para calcular cantidades relacionadas con las VA. Curvas de luz en estrellas binarias y efectos en los tiempos de tránsitos Inferencia aproximada con métodos MonteCarlo 19 / 87
distribución Integral de Riemann–Stieltjes Estadística inferencial y estimadores Integral de Riemann–Stieltjes Cuando F es diferenciable, b a h(x) dF(x) = b a h(x) F (x) dx Si h es una función contínua y F es constante salvo en un número finito de puntos {xi }n i=1 de saltos positivos de tamaño p1 , p2 , . . . , pn, b a h(x) dF(x) = ∞ i=1 h(xi ) p(xi ) Curvas de luz en estrellas binarias y efectos en los tiempos de tránsitos Inferencia aproximada con métodos MonteCarlo 20 / 87
distribución Integral de Riemann–Stieltjes Estadística inferencial y estimadores Momentos de una F La integral de RS permite definir los momentos de una función de distribución. Mn = ∞ −∞ xn dF(x) Curvas de luz en estrellas binarias y efectos en los tiempos de tránsitos Inferencia aproximada con métodos MonteCarlo 21 / 87
distribución Integral de Riemann–Stieltjes Estadística inferencial y estimadores Momentos de una F El primer momento se denomina esperanza o valor de expectación de X, E(X) = ∞ −∞ x dFX (x) y es equivalente a calcular Ω X(ω)dP(ω) Para los casos contínuo y discreto, el valor de expectación es, respectivamente: E = ∞ −∞ x F (x) E = x x p(x) Curvas de luz en estrellas binarias y efectos en los tiempos de tránsitos Inferencia aproximada con métodos MonteCarlo 22 / 87
distribución Integral de Riemann–Stieltjes Estadística inferencial y estimadores Valor de expectación de una función Es posible además calcular el valor de expectación de una función: E[g(x)] = ∞ −∞ g(x) dF(x) Curvas de luz en estrellas binarias y efectos en los tiempos de tránsitos Inferencia aproximada con métodos MonteCarlo 23 / 87
distribución Integral de Riemann–Stieltjes Estadística inferencial y estimadores Inferencia estadística Curvas de luz en estrellas binarias y efectos en los tiempos de tránsitos Inferencia aproximada con métodos MonteCarlo 24 / 87
distribución Integral de Riemann–Stieltjes Estadística inferencial y estimadores Inferencia estadística Curvas de luz en estrellas binarias y efectos en los tiempos de tránsitos Inferencia aproximada con métodos MonteCarlo 25 / 87
distribución Integral de Riemann–Stieltjes Estadística inferencial y estimadores Muestras aleatorias Def./ Una muestra aleatoria (MA) es una colección de variables aleatorias que son independientes e idénticamente distribuidas. Una muestra aleatoria entonces es un conjunto de VA: {X1 , X2 , . . . , Xn } que suele denotarse también como {Xi }n i=1 Curvas de luz en estrellas binarias y efectos en los tiempos de tránsitos Inferencia aproximada con métodos MonteCarlo 26 / 87
distribución Integral de Riemann–Stieltjes Estadística inferencial y estimadores Estimadores Def./ Un estimador es una función cualquiera de una muestra aleatoria Se suele denotar genéricamente a un parámetro de una distribución como θ, y al correspondiente estimador como ˆ θ. Un estimador sirve para estimar el valor de un parámetro. e.g.: Algunos parámetros tienen una notación propia, por ejemplo, θ = µ → ˆ θ = ¯ X. Curvas de luz en estrellas binarias y efectos en los tiempos de tránsitos Inferencia aproximada con métodos MonteCarlo 27 / 87
distribución Integral de Riemann–Stieltjes Estadística inferencial y estimadores Estimadores puntuales Cuando se hace una estimación numérica de un parámetro a partir de una MA, se obtiene un estimador puntual. Propiedades deseables de un estimador puntual Insesgabilidad: se dice que el estimador es insesgado si E[ˆ θ] = θ. Consistencia: si al aumentar la muestra, el estimador se aproxima al parámetro. Eficiencia: se calcula comparando su varianza con la de otro estimador. Suficiencia: utiliza toda la información obtenida de la muestra. Curvas de luz en estrellas binarias y efectos en los tiempos de tránsitos Inferencia aproximada con métodos MonteCarlo 28 / 87
distribución Integral de Riemann–Stieltjes Estadística inferencial y estimadores e.g.: Estimador de la media muestral Dadas n observaciones: X1 , X2 , . . . , Xn, con una misma distribución, la media muestral se define por X(n) = X1 + X2 + · · · + Xn n . La media muestral se utiliza como un estimador de la media θ. Suele denotarse a la media de la población como µ y por lo tanto ¯ X = ˆ θ = ˆ µ. Curvas de luz en estrellas binarias y efectos en los tiempos de tránsitos Inferencia aproximada con métodos MonteCarlo 29 / 87
distribución Integral de Riemann–Stieltjes Estadística inferencial y estimadores e.g.: Estimador de la media muestral El estimador ¯ X es insesgado, consistente y tiende asintóticamente a la normal. • Es insesgado, es decir, µ = E[¯ X], si la media es finita. E[X(n)] = E n i=1 Xi n = n i=1 E[Xi ] n = nθ n = θ. • Es un estimador consistente, pues, por la ley de los grandes números, lim n→∞ [X(n)] = θ. • Y tiende asintóticamente a la normal, lo cual se deduce de aplicar el teorema del límite central. Z = √ n(X(n) − θ)/σ =⇒ lim n→∞ [FZ (n)(x)] = Φ(x) Curvas de luz en estrellas binarias y efectos en los tiempos de tránsitos Inferencia aproximada con métodos MonteCarlo 30 / 87
distribución Integral de Riemann–Stieltjes Estadística inferencial y estimadores Aplicaciones de PRNGs: Integración MonteCarlo Dado que ¯ X es un estimador insesgado de E(X): E(X) = b a f(x) dx ¯ X = 1 N N i=1 f(xi ), con xi ∼ f, a < xi < b Luego, I = b a f(x) dx I ≈ 1 N N i=1 f(xi ), con xi ∼ f, a < xi < b donde los números xi siguen una distribución f (son un sorte de esa distribución) y están en el intervalo [a, b]. Curvas de luz en estrellas binarias y efectos en los tiempos de tránsitos Inferencia aproximada con métodos MonteCarlo 31 / 87
distribución Integral de Riemann–Stieltjes Estadística inferencial y estimadores e.g.:Integración MonteCarlo Sup. que queremos calcular (para un N muy grande): I = N i=0 i Analíticamente, I = N(N + 1)/2. Pero si tuviéramos que calcularlo necesitaríamos sumar N términos. Esta cuenta se puede pensar como el valor de expectación de la función g(x) = 1/N de la variable X ∼ U(0, N). Como es una VA discreta, en lugar de una integral tenemos una sumatoria. Este valor de expectación se puede estimar mediante el estimador ¯ X de una M.A. de la variable X: I 1 N M i=0 xi , xi ∼ U(0, N) (ver notebook) Curvas de luz en estrellas binarias y efectos en los tiempos de tránsitos Inferencia aproximada con métodos MonteCarlo 32 / 87
con modelos no lineales The model selection problem Ahora, datos y modelos! Curvas de luz en estrellas binarias y efectos en los tiempos de tránsitos Inferencia aproximada con métodos MonteCarlo 33 / 87
con modelos no lineales The model selection problem Si se tiene un conjunto de datos y se propone un modelo paramétrico, en realidad hay toda una familia de modelos que se pueden usar para describir los datos. La PH permite determinar la valoración de la hipótesis de que un parámetro tenga cierto valor. El poder de predicción de un modelo depende de la exactitud del parámetro propuesto. Como problema general, se quiere saber cuál es la probabilidad de un modelo para un conjunto de datos. Pero eso no se puede saber, lo que sí se puede saber es la probabilidad de un conjunto de datos para un modelo dado. Curvas de luz en estrellas binarias y efectos en los tiempos de tránsitos Inferencia aproximada con métodos MonteCarlo 34 / 87
con modelos no lineales Teorema de Bayes Teorema Sea (Ω, F, P) un espacio de probabilidad, y {Ai }∞ i=1 una partición de Ω con P(Ai ) > 0 ∀i. Para cualquier evento B ∈ F, P(Ai |B) = P(B|Ai )P(Ai ) ∞ k=1 P(B|Ak ) P(Ak ) Curvas de luz en estrellas binarias y efectos en los tiempos de tránsitos Inferencia aproximada con métodos MonteCarlo 35 / 87
con modelos no lineales Teorema de Bayes Interpretación visual del teorema de Bayes: Curvas de luz en estrellas binarias y efectos en los tiempos de tránsitos Inferencia aproximada con métodos MonteCarlo 36 / 87
con modelos no lineales Teorema de Bayes Demostración Sea {Ai }∞ i=1 una partición de Ω con P(Ai ) > 0 ∀i y B ∈ F un evento cualquiera con P(B) = 0. El teorema de Bayes se puede demostrar fácilmente en dos pasos, usando la definición de probabilidad condicional y ell teorema de la probabilidad total: P(Ai |B) = P(Ai ∩ B) P(B) = P(B|Ai )P(Ai ) P(B) = P(B|Ai )P(Ai ) ∞ i=1 P(B|Ai ) P(Ai ) Curvas de luz en estrellas binarias y efectos en los tiempos de tránsitos Inferencia aproximada con métodos MonteCarlo 37 / 87
con modelos no lineales Bayesian model selection El teorema de Bayes se puede reescribir en términos de modelos (M) y datos (D): P(M|D) = P(D|M)P(M) P(D) Si se tiene una familia paramétrica de modelos, M θ P(M θ |D) = P(D|M θ )P(M θ ) P(D) Cabe preguntarse cuál es el conjunto de parámetros θ que producen la mayor probabilidad P(M θ |D) Curvas de luz en estrellas binarias y efectos en los tiempos de tránsitos Inferencia aproximada con métodos MonteCarlo 38 / 87
con modelos no lineales Bayesian model selection P(M|D) = P(D|M)P(M) P(D) Para calcular el Likelihood L = P(D|M), hay que tener en cuenta que si tenemos una MA, todas las observaciones son independientes, y por los tanto: P(D|M) = P({Xi }n i=1 |M) = P({X = X1 }|M) P({X = X2 }|M) ...P({X = Xn }|M) = πn i=1 P({X = Xi }|M) Curvas de luz en estrellas binarias y efectos en los tiempos de tránsitos Inferencia aproximada con métodos MonteCarlo 39 / 87
con modelos no lineales Bayesian model selection Esto es una función de los parámetros que se quieren maximizar: L(θ) = Πn i=1 P({X = Xi }|M) Sin pérdida de generalidad, se puede maximizar el logaritmo del Likelihood, ya que el logaritmo es una función no decreciente: log(L(θ)) = log(Πn i=1 P({X = Xi }|M)) log(L(θ)) = n i=1 log(P({X = Xi }|M)) Curvas de luz en estrellas binarias y efectos en los tiempos de tránsitos Inferencia aproximada con métodos MonteCarlo 40 / 87
con modelos no lineales e.g.: Errores gaussianos Sea X una MA, y sea f(x) = y(x) + un modelo lineal con errores gaussianos, es decir, si la relación teórica es y(x), f(x) = N[y(x; a, b), σ2] donde E( ) = 0 y E( 2) = σ2. =⇒ f(x) = A exp − (yi − y(x; a, b))2 √ 2σ donde A = 1 √ 2πσ es una constante de integración para que la distribución tenga area igual a 1. Curvas de luz en estrellas binarias y efectos en los tiempos de tránsitos Inferencia aproximada con métodos MonteCarlo 41 / 87
con modelos no lineales e.g.: Errores gaussianos Entonces, para un dato xk, P(yk ) dyk = N(yk ; σk )(xk ) = A exp − (yi − y(x; θ))2 √ 2σ y por lo tanto, log(P(D|M)) = log( n k=1 P(yk )) = log( n k=1 N(yk ; σk )(xk ) = n k=1 log(N(yk ; σk )(xk )) Curvas de luz en estrellas binarias y efectos en los tiempos de tránsitos Inferencia aproximada con métodos MonteCarlo 42 / 87
con modelos no lineales e.g.: Errores gaussianos Gracias a que el logaritmo es una función no decreciente, tenemos que: max(P(D|M)) = min(− log(P(D|M))) donde: − log(P(D|M)) = n log √ 2πσ dy + n k=1 (yi − y(x; θ))2 √ 2σ Por lo tanto, los parámetros que buscamos son los que minimizan los cuadrados de las diferencias entre los valores observados y los valores predichos. θFIT = minθ n k=1 (yi − y(x; θ))2 √ 2σ este procedimiento se denomina de cuadrados mínimos. Curvas de luz en estrellas binarias y efectos en los tiempos de tránsitos Inferencia aproximada con métodos MonteCarlo 44 / 87
con modelos no lineales Least squares El problema consiste ahora en hallar una forma práctica de encontrar el conjunto de parámetros que minimizan el χ2. χ2 = n i=1 (yi − M(xi , θ))2 2σ2 i La fórmula depende de la forma funcional de M(x, θ). Si M(x, θ) depende de manera lineal de los parámetros, se realiza un ajuste por cuadrados mínimos lineal. Esto requiere resolver un sistema de ecuaciones matricial, para lo cual se usan métodos de análisis numérico aplicado al álgebra lineal. Si M(x, θ) depende de manera no lineal de los parámetros, se realiza un ajuste por cuadrados mínimos no lineal. Estos métodos son iterativos y se realizan de manera numérica. Curvas de luz en estrellas binarias y efectos en los tiempos de tránsitos Inferencia aproximada con métodos MonteCarlo 45 / 87
con modelos no lineales Linear least squares: e.g.:ajuste lineal Para minimizar χ2 se debe plantear una forma para el modelo: χ2 = n i=1 (yi − M(xi , θ))2 2σ2 i E(Y|X) = M(x; θ) + donde la variable es el ruido, con cierta distribución. Sea un modelo lineal de la forma Y = α + βx, queremos entonces minimizar: χ2 = n i=1 (yi − M(xi ; α, β))2 2σ2 i = n i=1 (yi − (α + βx))2 2σ2 i Curvas de luz en estrellas binarias y efectos en los tiempos de tránsitos Inferencia aproximada con métodos MonteCarlo 46 / 87
con modelos no lineales Linear least squares: e.g.:ajuste lineal Curvas de luz en estrellas binarias y efectos en los tiempos de tránsitos Inferencia aproximada con métodos MonteCarlo 47 / 87
con modelos no lineales Linear least squares: e.g.:ajuste lineal o sea, buscamos α y β tales que: ∂χ2 ∂α = 0; ∂χ2 ∂β = 0 Si definimos S = n i=1 1 σ2 n , Sx = n i=1 xn σ2 n , Sy = n i=1 yn σ2 n , Sxx = n i=1 x2 n σ2 n , Sxy = n i=1 xn yn σ2 n , obtenemos que: ∂χ2 ∂α = Sy + Sα − Sx β, ∂χ2 ∂β = Sxy − Sx α − Sxx β =⇒ α = Sxx Sy − Sx Sxy S Sxx − S2 x , β = Sxy S − Sx Sy S Sxx − S2 x Curvas de luz en estrellas binarias y efectos en los tiempos de tránsitos Inferencia aproximada con métodos MonteCarlo 48 / 87
con modelos no lineales Linear least squares Una generalización de la regresión lineal es aquella que considera un modelo que es una combinación lineal de muchas funciones. Por ejemplo, si f1 (x) = 1, f2 (x) = x, f3 (x) = x2, fM (x) = xM−1, y = a1 + a2 f2 (x) + a3 f3 (x) + . . . + aM fM (x) y = a1 + a2 x + a3 x2 + . . . + aM xM−1 Es decir, M(x; θ) = M k=1 ak fk (x) Curvas de luz en estrellas binarias y efectos en los tiempos de tránsitos Inferencia aproximada con métodos MonteCarlo 49 / 87
con modelos no lineales Linear least squares Para minimizar χ2 se debe planear una forma para el modelo: χ2(θ) = n i=1 (yi − M(xi , θ))2 σi M(x; θ) = M k=1 ak fk (x) χ2(θ) = n i=1 (yi − M k=1 ak fk (s))2 σi χ2(θ) = n i=1 yi σi − 1 σi M k=1 ak fk (x) 2 Curvas de luz en estrellas binarias y efectos en los tiempos de tránsitos Inferencia aproximada con métodos MonteCarlo 50 / 87
con modelos no lineales Linear least squares χ2(θ) = n i=1 yi σi − 1 σi M k=1 ak fk (x) 2 se puede escribir como χ2 = |Ax − b|2 con A = f1 (x1 )/σ1 . . . . . . fM (x1 )/σ1 f1 (x2 )/σ2 ... . . . . . . . . . f1 (x2 )/σM . . . . . . fM (xn )/σM a = a1 a2 . . . aM b = y1 /σ1 y2 /σ2 . . . yM /σM Curvas de luz en estrellas binarias y efectos en los tiempos de tránsitos Inferencia aproximada con métodos MonteCarlo 51 / 87
con modelos no lineales Linear least squares El mínimo de χ2 ocurre cuando las derivadas parciales respecto de todos los (M) parámetros se anulan, ∂χ2 ∂ai = 0; ∀ai χ2(θ) = n i=1 yi σi − 1 σi M k=1 ak fk (x) 2 Esto produce un sistema de ecuaciones: 0 = N i=1 1 σ2 i yi − M k=1 ak fk (xi ) fk (xi ) k = 1, . . . , M Curvas de luz en estrellas binarias y efectos en los tiempos de tránsitos Inferencia aproximada con métodos MonteCarlo 52 / 87
con modelos no lineales Linear least squares Se puede definir: αkj = N i=1 fj (xi ) fk (xi ) σ2 i βk = N i=1 yi fk (xi ) σ2 i que son una matriz MxM ([α] = AT · A) y un vector de tamaño N ([β] = AT · b). Con estas definiciones, la ecuación 0 = N i=1 1 σ2 i yi − M k=1 ak fk (xi ) fk (xi ) k = 1, . . . , M se puede reescribir como: M j=1 αkj aj = βk Curvas de luz en estrellas binarias y efectos en los tiempos de tránsitos Inferencia aproximada con métodos MonteCarlo 53 / 87
con modelos no lineales Linear least squares M j=1 αkj aj = βk Existen métodos numéricos para resolver este sistema de ecuaciones: Descomposición LU Descomposición de Cholesky Eliminación de Gauss-Jordan Descomposición en valores singulares ... y en general cualquier método numérico que sirva para el caso particular. Curvas de luz en estrellas binarias y efectos en los tiempos de tránsitos Inferencia aproximada con métodos MonteCarlo 54 / 87
con modelos no lineales Non–linear least squares: e.g.: Qué pasa si el modelo no es lineal con los parámetros? χ2 = n i=1 (yi − M(xi , θ))2 2σ2 i donde M(xi , θ) no es lineal con los parámetros θi. M(x, θ) = θ1 + θ2 exp(−θ3 x2) M(x, θ) = θ1 + θ2 x 1 + θ3 x) M(x, θ) = θ1 sin(θ2 x) + θ3 cos(θ3 x) M(x, θ) = θ1 xθ2 Curvas de luz en estrellas binarias y efectos en los tiempos de tránsitos Inferencia aproximada con métodos MonteCarlo 55 / 87
con modelos no lineales Non–linear least squares Existen varios métodos para encontrar los parámetros θ que dan el mínimo. Por ejemplo: Métodos analíticos Método de Newton (o Gauss–Newton) Método del gradiente descendiente Método “steepest descent” Método “stochastic gradient descent” Método de Levenberg–Marquardt Método Simplex, o de Nelder–Mead (multidimensional) Método de Powell Método de los gradientes conjugados (Fletcher) Curvas de luz en estrellas binarias y efectos en los tiempos de tránsitos Inferencia aproximada con métodos MonteCarlo 56 / 87
con modelos no lineales Non–linear least squares Hay problemas en donde no es posible implementar esos métodos, ya sea porque no se pueden encontrar las derivadas, porque tiene muchos parámetros, etc. En esos casos se recurre a los métodos Monte Carlo. Factores que influyen en la elección del método: Función de costo: es fácil de calcular? Derivadas: es fácil/posible de calcularlas? Hay límites conocidos para los parámetros? No hay UN método perfecto Curvas de luz en estrellas binarias y efectos en los tiempos de tránsitos Inferencia aproximada con métodos MonteCarlo 57 / 87
con modelos no lineales Non–linear least squares Algunas consideraciones: Resolver analíticamente siempre que se pueda Hay métodos para una o varias dimensiones Algunos métodos numéricos requieren el cálculo de derivadas Algunos métodos numéricos son iterativos El método más conveniente depende del problema Dominar el arte de resolver este problema significa: Encontrar el mejor método Determinar si el método convergió a la mejor solución Curvas de luz en estrellas binarias y efectos en los tiempos de tránsitos Inferencia aproximada con métodos MonteCarlo 58 / 87
Sampling Sampling algorithms Sampling methods El método de expectation–maximization es un caso particular de maximum likelihood estimation, cuando hay variables latentes u ocultas (datos incompletos). Recordemos que el maximum likelihood estimation se basa en el Teorema de Bayes, que relaciona: la probabilidad posterior: P(θ|D, M), el Likelihood: L(θ) = P(D|M; θ), los priors: P(θ|M), y la evidencia: P(D|M). Curvas de luz en estrellas binarias y efectos en los tiempos de tránsitos Inferencia aproximada con métodos MonteCarlo 59 / 87
Sampling Sampling algorithms Sampling methods P(θ|D, M) = L(θ) P(θ|M) P(D|M) La evidencia es irrelevante en cualquier procedimiento de inferencia paramétrica, pero adquiere relevancia cuando se realiza una selección de modelos. Una vez que por algún procedimiento se obtienen las estimaciones de los parámetros θ, para dar la estimación de un parámetro θi, dado el resto de los parámetros θ[i] , tal que θ = (θi , θ[i] ), hay que hacer: P(θi |D, M) = L(θi , θ[i] )P(θi , θ[i] |M)dθ[i] Se desea dar una estimación de θi usando la media, la mediana, etc., y su error estándar. Cuando sea posible, lo ideal es dar su distribución. Curvas de luz en estrellas binarias y efectos en los tiempos de tránsitos Inferencia aproximada con métodos MonteCarlo 60 / 87
Sampling Sampling algorithms Sampling methods P(θi |D, M) = L(θi , θ[i] )P(θi , θ[i] |M)dθ[i] En los casos reales, lo normal es que este procedimiento no se pueda realizar de manera analítica, por lo que hay que recurrir a métodos numéricos para calcular el Likelihood y la probabilidad posterior. Dado que lo que se busca es dar una medida de la tendencia central de la distribución del parámetro, se puede usar una estimación del valor de expectación. El mismo se puede estimar usando métodos MonteCarlo. Estudiamos ahora métodos de inferencia aproximada usando métodos de muestreo numérico. Los métodos de muestreo permiten tener una estima de la probabilidad posterior, pero se usan fundamentalmente para evaluar valores de expectación. Curvas de luz en estrellas binarias y efectos en los tiempos de tránsitos Inferencia aproximada con métodos MonteCarlo 61 / 87
Sampling Sampling algorithms Estimación de valores de expectación Como planteo general, queremos evaluar el valor de expectación de una función f(z) de una VA que tiene distribución p(z). E[ f ] = f(z)dp(z) Si tengo un conjunto de muestras independientes tomadas de la distribución p(z), {zl}L l=1 , entonces ˆ f 1 L L l=1 f(z(l)) Este es un estimador MonteCarlo de E[f], y si el estimador es insesgado, E[ˆ f] = E[f] La varianza del estimador está dada por: Var[ˆ f] = 1 L E[(ˆ f − E[ˆ f])2] Curvas de luz en estrellas binarias y efectos en los tiempos de tránsitos Inferencia aproximada con métodos MonteCarlo 62 / 87
Sampling Sampling algorithms Estimación de valores de expectación Podemos usar varios métodos de muestreo (samplers) para hacer esto, por ejemplo: Método de la función inversa Método de aceptación rechazo (Box–Muller) (importance sampling) Markon Chain MonteCarlo: Metrópolis–Hastings algorithm Markon Chain MonteCarlo: Gibbs sampler Markon Chain MonteCarlo: Slice sampling Estos métodos se podrán aplicar en función del problema. Curvas de luz en estrellas binarias y efectos en los tiempos de tránsitos Inferencia aproximada con métodos MonteCarlo 63 / 87
Sampling Sampling algorithms Importance Sampling Una de las razones para sortear números randoms de una distribución complicada es calcular una valor de expectación: E[ f ] = f(z)dp(z) que se puede aproximar a partir de una muestra aleatoria z(l) tomada de la distribución p(z), por: ˆ f 1 L L l=1 f(z(l)) Curvas de luz en estrellas binarias y efectos en los tiempos de tránsitos Inferencia aproximada con métodos MonteCarlo 64 / 87
Sampling Sampling algorithms Importance Sampling En algunas situaciones es difícil, impráctico o costoso sortear randoms de p(z), pero es fácil evaluar p(z). En tal caso sería fácil aproximar la integral de f por una suma de Riemman: E[ f ] L l=1 p(z(l))f(z(l)) El problema con esta estrategia es que para muchas dimensiones el tamaño del grid se hace muy grande, y además si la función es compacta, la mayoría contribuirá muy poco a la suma. Sería ideal elegir puntos en donde el producto p(z(l))f(z(l)) es grande... Curvas de luz en estrellas binarias y efectos en los tiempos de tránsitos Inferencia aproximada con métodos MonteCarlo 65 / 87
Sampling Sampling algorithms Importance Sampling La técnica de importance sampling permite calcular una aproximación a un valor de expectación directamente, sin mediar el sorteo de la distribución p(z). Para eso se propone una función, q(z) (la proposal distribution) de la cual es fácil sortear randoms (parecido al método de aceptación/rechazo). Curvas de luz en estrellas binarias y efectos en los tiempos de tránsitos Inferencia aproximada con métodos MonteCarlo 66 / 87
Sampling Sampling algorithms Importance Sampling E[f] = f(z)p(z)dz = f(z) p(z) q(z) q(z)dz 1 L L l=1 p(z(l)) q(z(l)) q(z(l))dz Las cantidades p(z(l)) q(z(l)) son los pesos que compensan el hecho de que estamos sorteando randoms de la distribución ”equivocada”. Esta técnica funciona bien cuando la proposal distribution es parecida a la función p(z). Curvas de luz en estrellas binarias y efectos en los tiempos de tránsitos Inferencia aproximada con métodos MonteCarlo 67 / 87
Sampling Sampling algorithms Sampling methods: MonteCarlo Para evaluar una integral con el método MonteCarlo, es necesario tener una secuencia de puntos en el espacio de parámetros cuya densidad sea proporcional a la función que se quiere integrar, en nuestro caso la probabilidad posterior P(θ|D, M). Los métodos MonteCarlo funcionan bien aún en condiciones "adversas", por ejemplo: El Likelihood es difícil de calcular, o sólo puede obtenerse mediante simulación hay muchos parámetros la hipersuperficie de la probabilidad posterior tiene muchos picos la hipersuperficie de la probabilidad posterior tiene una estructura complicada Curvas de luz en estrellas binarias y efectos en los tiempos de tránsitos Inferencia aproximada con métodos MonteCarlo 68 / 87
Sampling Sampling algorithms MonteCarlo Markov Chains Los métodos MonteCarlo permiten obtener una secuencia de puntos que siguen la distribución P(θ|D, M). Estos conjuntos de puntos no son independientes, sino que presentan una dependencia explícita con el punto anterior en la secuencia, formando una especie de cadena. Esta propiedad de ”no memoria” se denomina ”de Markov” o ”Markoviana”. e.g.:El movimiento Browniano es un proceso de Markov. Estos métodos se encuentran abarcados por los procedimientos de Cadenas de Markov MonteCarlo. Curvas de luz en estrellas binarias y efectos en los tiempos de tránsitos Inferencia aproximada con métodos MonteCarlo 69 / 87
Sampling Sampling algorithms MonteCarlo Markov Chains Una ventaja de estos métodos es que una vez que se tiene la cadena, es posible estimar el valor de expectación de cualquier función de los parámetros: E(θ) = (dθ)θ P(θ|D) E(θ) 1 M M i=1 θ(i) donde los θ(i) son las muestras de la cadena. E(G(θ)) 1 M M i=1 G(θ(i)) Es importante entender que el objetivo de una MCMC no es hacer un ajuste, sino dar una estimación de los parámetros del ajuste. Curvas de luz en estrellas binarias y efectos en los tiempos de tránsitos Inferencia aproximada con métodos MonteCarlo 70 / 87
Sampling Sampling algorithms MonteCarlo Markov Chains: formulación del problema La idea de estos métodos es usar también una proposal function (como en el caso del importance sampling). Sin embargo, esta vez no vamos a sortear cada vez un número random nuevo de la proposal distribution, sino que vamos a generar uno a partir del anterior. Es decir, si en el estado actual tenemos un valor zτ , usamos una proposal function q(z|zτ ) para generar zτ+1, de tal forma que la secuencia z1, z2, . . . forma una cadena de Markov. La proposal distribution se elige de tal forma que sea fácil sortear un número directamente. Una vez que se genera un candidato z∗, se acepta de acuerdo a determinados criterios. Curvas de luz en estrellas binarias y efectos en los tiempos de tránsitos Inferencia aproximada con métodos MonteCarlo 71 / 87
Sampling Sampling algorithms MCMC: Metrópolis En el algoritmo propuesto por Metrópolis, se asume que la distribución propuesta es simétrica: q(x | x ) = q(x | x) Dado un estado x se sortea otro estado x∗, el cual se incorpora a la cadena con probabilidad A(x(τ), x∗) = min 1, P(x∗) P(x(τ)) Esto se puede hacer sorteando un número aleatorio entre 0 y 1, Si no se acepta un punto nuevo, se repite el anterior. Curvas de luz en estrellas binarias y efectos en los tiempos de tránsitos Inferencia aproximada con métodos MonteCarlo 72 / 87
Sampling Sampling algorithms MonteCarlo Markov Chains Curvas de luz en estrellas binarias y efectos en los tiempos de tránsitos Inferencia aproximada con métodos MonteCarlo 73 / 87
Sampling Sampling algorithms MCMC: Metrópolis Curvas de luz en estrellas binarias y efectos en los tiempos de tránsitos Inferencia aproximada con métodos MonteCarlo 74 / 87
Sampling Sampling algorithms MCMC: why Metrópolis? Porqué funciona el algoritmo de Metrópolis? Recordemos que el algoritmo busca hacer un muestreo de la probabilidad posterior (target density) de un modelo paramétrico, para un dado conjunto de datos: P(θ|D) Usando el teorema de Bayes: P(θ|D) = L(D|θ)P(θ) P(D) Para lo cual vamos a calcular el Likelihood, L(D|θ) y usar una proposal function (jump function), que debe ser simétrica: q(θ∗|θ) Curvas de luz en estrellas binarias y efectos en los tiempos de tránsitos Inferencia aproximada con métodos MonteCarlo 75 / 87
Sampling Sampling algorithms MCMC: why Metrópolis? En el i-ésimo paso del algoritmo de Metrópolis, se propone un punto candidato que se acepta con probabilidad: A(θ∗|θi−1 ) = min 1, P(θ∗|D) P(θi−1 |D) Si P(θ∗|D) > P(θi−1 |D), entonces A(θ∗|θi−1 ) = 1 y se acepta el valor propuesto, θi = θ∗ Si P(θ∗|D) < P(θi−1 |D), entonces se acepta el valor propuesto θi = θ∗, con probabilidad igual a P(θ∗|D) P(θi−1 |D) . Si no se acepta entonces θi = θi−1. Curvas de luz en estrellas binarias y efectos en los tiempos de tránsitos Inferencia aproximada con métodos MonteCarlo 76 / 87
Sampling Sampling algorithms MCMC: why Metrópolis? Supongamos que en el i-ésimo paso, tenemos dos puntos sorteados de la probabilidad posterior, θa y θb, y supongamos que P(θb |D) > P(θa |D). Entonces, P(estar en θa y saltar a θb en el paso i) = P(θi−1 = θa , θi = θb ) = P(estar en θa )P(proposal)P(aceptar) = P(θa |D) q(θb |θa ) y además, P(estar en θb y saltar a θa en el paso i) = P(θi−1 = θb , θi = θa ) = P(estar en θb )P(proposal)P(aceptar) = P(θb |D) q(θb |θa ) P(θa |D) P(θb |D) = P(θa |D) q(θb |θa ) Curvas de luz en estrellas binarias y efectos en los tiempos de tránsitos Inferencia aproximada con métodos MonteCarlo 77 / 87
Sampling Sampling algorithms MCMC: why Metrópolis? Pero como la función de salto es simétrica: q(θb |θa ) = q(θa |θb ) debe ser: P(θi−1 = θa , θi = θb ) = P(θi−1 = θb , θi = θa ) Entonces, las distribuciones conjuntas de dos muestras consecutivas son simétricas, θi−1 y θi tienen la misma distribución marginal y P(θ|D) es la distribución estacionaria de la cadena de Markov. Curvas de luz en estrellas binarias y efectos en los tiempos de tránsitos Inferencia aproximada con métodos MonteCarlo 78 / 87
Sampling Sampling algorithms MCMC: why Metrópolis? Algunas referencias para lectura adicional: https://arxiv.org/pdf/1710.06068.pdf Curvas de luz en estrellas binarias y efectos en los tiempos de tránsitos Inferencia aproximada con métodos MonteCarlo 79 / 87
Sampling Sampling algorithms MonteCarlo Markov Chains: Algoritmos Hemos visto las propiedades de las cadenas de Markov y cómo usarlas para estimar parámetros a partir de la probabilidad posterior. pero... cómo se hace para construir una cadena? Hay dos algoritmos principales para construir cadenas de Markov en espacios de estados contínuos: Algoritmo de Metrópolis–Hastings: El algoritmo de Metrópolis–Hastings es una generalización del algoritmo de Metrópolis, en donde la función de salto no necesariamente es simétrica. Algoritmo de Gibbs: Usa las distribuciones marginales. Curvas de luz en estrellas binarias y efectos en los tiempos de tránsitos Inferencia aproximada con métodos MonteCarlo 80 / 87
Sampling Sampling algorithms MonteCarlo Markov Chains: Gibbs sampler En el algoritmo de muestreo de Gibbs, se parte de un dado estado de la cadena y se reemplaza el valor de uno de los parámetros por otro valor sorteado de la distribución marginal en ese parámetro, condicionado a los valores de los demás parámetros en el paso actual. Es decir, si z es el vector de los parámetros, se toma uno de los elementos zi y se sortea de la distribución marginal en zi fijando todos los otros parámetros z¬i (zi , z2 , . . . , zi , . . . , zK ) → (zi , z2 , . . . , zi , . . . , zK ) Este procedimiento se repite de manera cíclica sobre todas las variables en algún order particular, o bien eligiendo de manera aleatoria una variable cada vez. Curvas de luz en estrellas binarias y efectos en los tiempos de tránsitos Inferencia aproximada con métodos MonteCarlo 81 / 87
Sampling Sampling algorithms MonteCarlo Markov Chains: Gibbs sampler Curvas de luz en estrellas binarias y efectos en los tiempos de tránsitos Inferencia aproximada con métodos MonteCarlo 83 / 87
Sampling Sampling algorithms Sampling methods P(θ|D, M) = L(θ) P(θ|M) P(D|M) Una vez que por algún procedimiento se obtienen las estimaciones de los parámetros θ, para dar la estimación de un parámetro θi, dado el resto de los parámetros θ[i] , tal que θ = (θi , θ[i] ), hay que hacer: P(θi |D, M) = L(θi , θ[i] )P(θi , θ[i] |M)dθ[i] Se desea dar una estimación de θi usando la media, la mediana, etc., y su error estándar. Cuando sea posible, lo ideal es dar su distribución. Curvas de luz en estrellas binarias y efectos en los tiempos de tránsitos Inferencia aproximada con métodos MonteCarlo 84 / 87
Sampling Sampling algorithms Sampling methods En los casos reales, lo normal es que este procedimiento no se pueda realizar de manera analítica, por lo que hay que recurrir a métodos numéricos para calcular el Likelihood y la probabilidad posterior. Dado que lo que se busca es dar una medida de la tendencia central de la distribución del parámetro, se puede usar una estimación del valor de expectación. El mismo se puede estimar usando métodos MonteCarlo. Estudiamos ahora métodos de inferencia aproximada usando métodos de muestreo numérico. Los métodos de muestreo permiten tener una estima de la probabilidad posterior, pero se usan fundamentalmente para evaluar valores de expectación. Curvas de luz en estrellas binarias y efectos en los tiempos de tránsitos Inferencia aproximada con métodos MonteCarlo 85 / 87
Sampling Sampling algorithms Estimación de valores de expectación Como planteo general, queremos evaluar el valor de expectación de una función f(z) de una VA que tiene distribución p(z). E[ f ] = f(z)dp(z) Si tengo un conjunto de muestras independientes tomadas de la distribución p(z), {zl}L l=1 , entonces ˆ f 1 L L l=1 f(z(l)) Este es un estimador MonteCarlo de E[f], y si el estimador es insesgado, E[ˆ f] = E[f] La varianza del estimador está dada por: Var[ˆ f] = 1 L E[(ˆ f − E[ˆ f])2] Curvas de luz en estrellas binarias y efectos en los tiempos de tránsitos Inferencia aproximada con métodos MonteCarlo 86 / 87
Sampling Sampling algorithms MonteCarlo Markov Chains: Bibliografía Bibliografía: Simulation, Cap. X A concise course on statistics, Cap. 24 Bayes book, Cap. 4-6 Pattern recognition, Cap. 11 Elements of statistical learning, Cap. 8 Introduction to probability model, Cap. 4 Machine Learning, a probabilistic perspective, Cap. 17 Procesos 2012, Cap. 3 Curvas de luz en estrellas binarias y efectos en los tiempos de tránsitos Inferencia aproximada con métodos MonteCarlo 87 / 87
mlares
sip3ristex
trycycle
shimacos
lw
saltcooky12
kentaitakura
konakalab
techmathproject
addyosmani
palkan
pedronauck
brettharned
smashingmag
brad_frost
sugarenia
marcduiker
vanstee
skipperchong
soudai
62gerente
