Cuantificación de incertidumbre sobre parámetros en modelos no lineales


==Cuantificación de incertidumbre sobre parámetros en modelos no lineales==

'''Lic. Rodrigo Gonzaga Sierra 

Dr. José del Carmen Jiménez Hernández 

Dr. José Andrés Christen Gracia'''

==1 Introducción==

De acuerdo con <span id='citeF-1'></span>[[#cite-1|[1]]], tradicionalmente se han modelado los problemas de ciencias, ingeniería, medioambiente y otras aplicaciones mediante modelos matemáticos deterministas que describen leyes naturales subyacentes. En la actualidad, se tiende cada vez más a incorporar algún tipo de incertidumbre para representar la falta de conocimiento sobre parámetros y datos físicos relevantes, variaciones aleatorias en las condiciones de operación o pura ignorancia sobre cómo debiera ser el modelo en realidad. 

Suponga que se cuenta con observaciones <math display="inline">\mathbf{Y}=(y_1,\dots ,y_n),</math> tomadas en los tiempos <math display="inline">\mathbf{t}=(t_1, \dots , t_n)</math>, de un fenómeno representado por medio de un sistema de ecuaciones con la siguiente estructura:

<span id="eq-1"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>y_i = \mathcal{H} (X_{\theta }(t_i)) + \varepsilon _i,   i=\overline{1,n},      </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (1)
|}

donde <math display="inline">\mathcal{H}</math> es el funcional de observaciones, es común por ejemplo tener <math display="inline">\varepsilon _i \sim _{i.i.d.} \mathcal{N}(0, \sigma ^{2})</math> y <math display="inline">X_{\theta }</math> es la solución del siguiente sistema de ecuaciones diferenciales ordinarias, es decir; el regresor o el '''modelo directo (forward)''':

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\frac{dX_{\theta }}{dt} = F(X_{\theta },t, \theta );   X_{\theta }(t_0) = X_0.     </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (2)
|}

Conociendo el verdadero valor de <math display="inline">\theta </math> y las condiciones iniciales, resolver y conocer a <math display="inline">X_{\theta }</math> se le llama un '''problema directo'''. 

Sin embargo, en este caso el objetivo es hacer inferencia sobre <math display="inline">\theta </math> a partir de las observaciones <math display="inline">Y,</math> por eso se habla de un '''problema inverso''' <span id='citeF-2'></span>[[#cite-2|[2]]]. 

Se puede considerar este problema como un mapeo:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>  \mathcal{F}_{\mathbf{t}} (\theta ) = (\mathcal{H} (X_{\theta } (t_{1})), \dots ,\mathcal{H} (X_{\theta } (t_n))), </math>
|}
|}

este es el “mapeo del modelo directo”.

El mapeo inverso está, en general; mal planteado y no tiene mucho sentido:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>  \mathcal{F}_{\mathbf{t}}^{-1} (y_1,\dots ,y_n) = \theta ,</math>
|}
|}

por eso es necesario realizar una estrategia de regularización, como la cuantificación de la incertidumbre mediante inferencia bayesiana.

Para resolver el problema de inferencia dado en la ecuación [[#eq-1|(1)]] se propone que las entradas de <math display="inline">\theta </math> sean variables aleatorias que siguen alguna distribución de probabilidad. Esto no necesariamente tiene un significado físico o intrínseco, sólo se sabe que se tiene incertidumbre acerca de los valores que pueden tomar estas variables aleatorias, y que la distribución de probabilidad cuantifica su incertidumbre. En este  sentido, en el presente trabajo se aborda el problema de cuantificar la incertidumbre desde la perspectiva Bayesiana, en donde se obtiene la distribución ''a posteriori'' de los parámetros de interés.

==2 Teoría Bayesiana==

El trabajo de Thomas Bayes, publicado de manera póstuma en 1763, ha tenido una importante consecuencia en la forma de hacer inferencia estadística, este provee una manera formal de combinar el conocimiento ''a priori'' (o inicial) que se tiene sobre un fenómeno, con el nuevo conocimiento que se adquiere a partir de nuevos datos y mediciones sobre el mismo, obteniendo así un conocimiento ''a posteriori'' (o final). Es decir, el conocimiento ''a priori'' se actualiza con la nueva información, y dicho conocimiento ''a posteriori'' se convertirá en el nuevo conocimiento ''a priori'', a la espera, otra vez, de nueva información que lo actualice. En este sección se utilizaron las siguientes bibliografías:  <span id='citeF-3'></span>[[#cite-3|[3]]], <span id='citeF-4'></span>[[#cite-4|[4]]], <span id='citeF-5'></span>[[#cite-5|[5]]], <span id='citeF-6'></span>[[#cite-6|[6]]] y <span id='citeF-7'></span>[[#cite-7|[7]]].

===2.1 Distribución ''a priori'' y ''a posteriori''===

En estadística bayesiana, el término común para referirse a la información con la que cuenta el investigador es el de '''información subjetiva''', y es importante aclarar, al menos brevemente, qué se entiende en este contexto por el adjetivo “subjetiva”, ya que puede tener una connotación distinta a la que se requiere bajo el enfoque bayesiano. 

Al hablar de información subjetiva se refiere a toda aquella información ''a priori'' que se tiene en relación al fenómeno aleatorio de interés, antes de recolectar o realizar nuevas mediciones sobre el mismo, y esto incluye: datos históricos, teorías, opiniones y conjeturas de expertos, conclusiones basadas en estudios previos.

El primer paso en la inferencia estadística bayesiana es traducir todo lo anterior en una distribución de probabilidad '''a priori (o inicial)'''. El segundo paso consiste en recolectar o realizar nuevas mediciones, y actualizar la distribución de probabilidad a priori, para obtener, mediante el teorema de Bayes, una distribución de probabilidad '''a posteriori (o final)'''.

Será esta última la mejor descripción posible de la incertidumbre, de acuerdo a toda la información disponible, y por tanto, será la herramienta fundamental a partir de la cual se realiza inferencia estadística.

Para referirse a un '''modelo probabilístico paramétrico general''' se denota como <math display="inline">p_{X|\Theta }(x | \theta )</math>, donde la función <math display="inline">p_{X|\Theta }( \cdot | \theta )</math> puede ser una función de masa de probabilidades de una variable (o vector) aleatoria discreta o bien una función de densidad de una variable aleatoria continua. El escribir dicha función condicional en el parámetro (o vector de parámetros) <math display="inline">\theta </math> se debe al hecho de que, una vez dado un valor específico de <math display="inline">\theta </math>, la función de probabilidad queda totalmente determinada. 

Para referirse a '''una muestra aleatoria (m.a.)''' se utilizará la notación <math display="inline">\mathbf{X} = (X_1, \ldots , X_n)</math> y para referirse a una '''observación muestral''' se utilizará <math display="inline">\mathbf{x} = (x_1, \ldots , x_n)</math>. Por '''espacio paramétrico''' se entenderá como el conjunto <math display="inline">\Theta </math> de todos los valores que puede tomar <math display="inline">\theta </math>, y por familia paramétrica se entenderá como un conjunto <math display="inline">P = \{ p_{X|\Theta }(x | \theta ) : \theta \in \Theta \} </math>.

Siguiendo a <span id='citeF-7'></span>[[#cite-7|[7]]], la estadística bayesiana modela la incertidumbre que se tiene sobre <math display="inline">\theta </math> probabilísticamente. Esto es, considere el valor de <math display="inline">\theta </math> como una variable (o vector) aleatoria con '''una distribución de probabilidad a priori''' ('''o inicial''') <math display="inline">p(\theta )</math>, de la misma forma se denotará solo como <math display="inline">p(\cdot )</math>, sin importar si <math display="inline">\theta </math> es una variable aleatoria discreta o continua. Se trata de una distribución basada en experiencia previa (experiencia de especialistas, datos históricos, etc.) antes de obtener datos. Luego se procede a observar los nuevos datos (obtención de la muestra) <math display="inline">\mathbf{x}=(x_1, \ldots , x_n)</math> y combina esta información con la distribución ''a priori'' mediante el teorema de Bayes y se obtiene una '''distribución de probabilidad a posteriori (o final)''':

<span id="eq-3"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>p_{\Theta |\mathbf{X}}(\theta | \mathbf{x}) = \frac{p_{\mathbf{X},\Theta }(\mathbf{x}, \theta )}{p_{\mathbf{X}}(\mathbf{x})} = \frac{p_{\mathbf{X}|\Theta }(\mathbf{x} | \theta ) \cdot p_{\Theta }(\theta )}{\int p_{\mathbf{X}|\Theta }(\mathbf{x} | \theta ) \cdot p_{\Theta }(\theta ) \, d \theta }. </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (3)
|}

Note que <math display="inline">p_{\Theta |\mathbf{X}}(\theta | \mathbf{x})</math> es también una distribución de probabilidad de <math display="inline">\theta </math>, pero a diferencia de la distribución a priori <math display="inline">p_{\Theta }(\theta )</math> toma en cuenta tanto la información contemplada en <math display="inline">p_{\Theta }(\theta )</math> como la información contenida en los datos observados <math display="inline">\mathbf{x} = (x_1, \ldots , x_n)</math>. La distribución ''a posteriori'' de la variable aleatoria <math display="inline">\Theta </math> es la base para hacer inferencias sobre <math display="inline">\theta </math>.

Tenga presente que, por un lado, la función de verosimilitud <math display="inline">p_{\mathbf{X}|\Theta }(\mathbf{x} | \theta )</math> y <math display="inline">p_{\Theta }(\theta )</math> son distribuciones de probabilidad, y por otro:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math> p_{\mathbf{X}}(\mathbf{x}) = \int p_{\mathbf{X}|\Theta }(\mathbf{x} | \theta ) \cdot p_{\Theta }(\theta ) \, d \theta , </math>
|}
|}

es la probabilidad (o densidad) conjunta de la muestra <math display="inline">\mathbf{x} = (x_1, \ldots , x_n)</math> observada a partir del vector aleatorio <math display="inline">\mathbf{X} = (X_1, \ldots , X_n)</math>. Pero hay que estar consciente de que <math display="inline">p(\mathbf{x})</math> es constante respecto a <math display="inline">\theta </math>, por lo que se puede escribir:

<span id="eq-4"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>p_{\Theta |\mathbf{X}}(\theta | \mathbf{x}) \propto p_{\mathbf{X}|\Theta }(\mathbf{x} | \theta ) \cdot p_{\Theta }(\theta ), </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (4)
|}

note que <math display="inline">p_{\mathbf{X}|\Theta }(\mathbf{x} | \theta ) = p_{X|\Theta }(x_1, x_2, \dots , x_n | \theta )</math>, es la probabilidad conjunta de la muestra condicional en <math display="inline">\theta </math>, llamada función de ''verosimilitud'', denotada también por <math display="inline">L(\theta |\mathbf{x})</math>. En el caso particular de que los componentes del vector aleatorio <math display="inline">\mathbf{X} = (X_1, \ldots , X_n)</math> resulten ser independientes, se tiene que:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>p_{\mathbf{X}|\Theta }(\mathbf{x} | \theta ) = \prod _{j=1}^{n} p_{X|\Theta }(x_j | \theta{).} </math>
|}
|}

Puede proponer como estimador puntual de <math display="inline">\theta </math> alguna medida de tendencia central, por ejemplo la mediana o la esperanza:

<span id="eq-5"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\hat{\theta } := E(\theta ) = \int \theta \, p_{\Theta |X}(\theta | x) \, d\theta{.} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (5)
|}

Y aún en el caso de que no se cuente con información muestral se puede calcular <math display="inline">\hat{\theta }</math> utilizando <math display="inline">p_{\Theta }(\theta )</math> en lugar de <math display="inline">p_{\Theta |\mathbf{X}}(\theta | \mathbf{x})</math>.

En algunos  casos en los que, en vez de conocer el vector de parámetros <math display="inline">\theta </math>, lo que interesa es describir el comportamiento de observaciones futuras del fenómeno aleatorio en cuestión, esto es, hacer '''predicción'''. 

Dado un valor de <math display="inline">\theta </math>, la distribución que describe el comportamiento de la observación futura <math display="inline">X_n</math> es <math display="inline">p_{X_n|\Theta }(x | \theta )</math>. El problema es que por lo general el valor de <math display="inline">\theta </math> es desconocido. Por lo regular, la estadística frecuentista aborda este problema estimando puntualmente a <math display="inline">\theta </math> con base en la muestra observada, y dicho estimador <math display="inline">\hat{\theta }</math> es sustituido en <math display="inline">p_{X_n|\Theta }(x | \hat{\theta })</math>. Desde la perspectiva bayesiana, el modelo <math display="inline">p_{X_n|\Theta }(x | \theta )</math> junto con la distribución a priori <math display="inline">p_{\Theta }(\theta )</math> inducen una distribución conjunta para el vector aleatorio <math display="inline">(X_n, \Theta )</math> mediante el concepto de probabilidad condicional:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math> p_{X_n,\Theta }(x, \theta ) = p_{X_n|\Theta }(x | \theta )p_{\Theta }(\theta{).} </math>
|}
|}

Si se marginaliza la distribución de probabilidad conjunta anterior se obtiene:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math> p_{X_n}(x) = \int _{\Theta } p_{X_n,\Theta }(x, \theta ) \, d\theta{.} </math>
|}
|}

De los dos resultados anteriores, se tiene:

<span id="eq-6"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>p_{X_n}(x) = \int _{\Theta } p_{X_n|\Theta }(x|\theta )p_{\Theta }(\theta ) \, d\theta{.} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (6)
|}

A <math display="inline">p_{X_n}(x)</math> se le denomina '''distribución predictiva a priori (o inicial)''' y describe el conocimiento acerca de una observación futura <math display="inline">X_n</math> basado únicamente en la información contenida en <math display="inline">p_{\Theta }(\theta )</math>. Nótese que <math display="inline">p_{X_n}(x)</math> no depende ya de <math display="inline">\theta </math>.

Para hacer estimación por regiones, por ejemplo, si desea calcular la probabilidad de que el vector de parámetros <math display="inline"> \theta </math> pertenezca a una región <math display="inline"> A \subset \Theta </math>:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math> P(\theta \in A) = \int _{A} p_{\Theta |\mathbf{X}}(\theta | \mathbf{x}) d\theta , </math>
|}
|}

o bien, dado un valor <math display="inline"> \delta \in (0, 1) </math>, se busca un <math display="inline"> A \subset \Theta </math> tal que <math display="inline"> P(\theta \in A) = \delta </math>. Con frecuencia la solución para <math display="inline"> A </math> no es única. Cabe aclarar que si <math display="inline"> \hbox{dim}(\Theta ) = 1 </math> las regiones son subconjuntos de <math display="inline"> \mathbb{R} </math> y que un caso particular de estas regiones son los intervalos. En este sentido, la estimación por regiones en estadística bayesiana es más general que la estimación por intervalos de la estadística frecuentista. 

Una vez obtenida la muestra, el modelo <math display="inline">p_{X|\Theta }(x | \theta )</math> y la distribución ''a posteriori'' <math display="inline">p_{\Theta |\mathbf{X}}(\theta |\mathbf{x})</math> inducen una distribución conjunta para <math display="inline">(\mathbf{X}, \Theta )</math> condicional en los valores observados <math display="inline">\mathbf{x} = (x_1, \dots , x_n)</math>:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>p_{X,\Theta |\mathbf{X}}(x, \theta | \mathbf{x})  = \frac{p_{X,\Theta ,\mathbf{X}}(x, \theta , \mathbf{x})}{p_{\mathbf{X}}(\mathbf{x})} </math>
|-
| style="text-align: center;" | <math>  = \frac{p_{X|\Theta ,\mathbf{X}}(x | \theta , \mathbf{x})p_{\Theta ,\mathbf{X}}(\theta , \mathbf{x})}{p_{\mathbf{X}}(\mathbf{x})} </math>
|-
| style="text-align: center;" | <math>  = p_{X|\Theta ,\mathbf{X}}(x | \theta , \mathbf{x})p_{\Theta |\mathbf{X}}(\theta | \mathbf{x}) </math>
|-
| style="text-align: center;" | <math>  = p_{X,\Theta }(x | \theta )p_{\Theta |\mathbf{X}}(\theta | \mathbf{x}).     </math>
|}
|}

En lo anterior, <math display="inline">p_{X|\Theta ,\mathbf{X}}(x | \theta , \mathbf{x}) = p_{X,\Theta }(x | \theta )</math> se justifica por la independencia condicional de <math display="inline">X</math> y <math display="inline">\mathbf{X} = (X_1, . . . ,X_n)</math> dado <math display="inline">\theta </math>. Si se marginaliza la distribución conjunta condicional anterior:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math> p_{X|\mathbf{X}}(x | \mathbf{x}) = \int p_{X,\Theta |\mathbf{X}}(x, \theta | \mathbf{x}) \, d\theta{.} </math>
|}
|}

Combinando los dos resultados anteriores:

<span id="eq-7"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>p_{X,\mathbf{X}}(x | \mathbf{x}) = \int p_{X|\Theta }(x | \theta )p_{\Theta |\mathbf{X}}(\theta | \mathbf{x}) \, d\theta . </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (7)
|}

A <math display="inline">p_{X|\mathbf{X}}(x|\mathbf{\hbox{x}})</math> se le llama '''distribución predictiva a posteriori (o final)''', y describe el conocimiento acerca de una observación futura <math display="inline">X</math> basado tanto en la información contenida en <math display="inline">p_{\Theta }(\theta )</math> como en la información muestral <math display="inline">\mathbf{x} = (x_1, . . . , x_n)</math>. Nótese nuevamente que <math display="inline">p_{X|\mathbf{X}}(x | \mathbf{x})</math> no depende de <math display="inline">\theta </math>.

Así que para hacer predicción sobre observaciones futuras del fenómeno aleatorio que esté modelando se usa <math display="inline">p_{X}(x)</math> o bien <math display="inline">p_{X|\mathbf{X}}(x | \mathbf{x})</math>, según sea el caso. Y de manera análoga a lo mencionado sobre inferencia bayesiana, una manera simple de hacer predicción puntual, por ejemplo, de una observación futura <math display="inline">X</math> podría ser mediante alguna medida de tendencia central, como la mediana o la esperanza:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math> \hat{x} = E(X) = \int _{\hbox{R}(X)} x \, p_{X|\mathbf{X}}(x | \mathbf{x}) \, dx, </math>
|}
|}

donde <math display="inline">\hbox{R}(X)</math> es el rango de la v.a. <math display="inline">X</math>. También, una manera de calcular la probabilidad de que una observación futura se encuentre en un conjunto <math display="inline">A \subseteq \hbox{R}(X)</math> sería:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math> P(\{ X \in A\} ) = \int _A p_{X|\mathbf{X}}(x | \mathbf{x}) dx. </math>
|}
|}

Las ecuaciones [[#eq-3|(3)]], [[#eq-6|(6)]] y [[#eq-7|(7)]] constituyen ''el modelo general de la estadística bayesiana''. Cualquier problema estadístico tratado bajo el enfoque bayesiano implica la obtención y utilización de las distribuciones correspondientes.

==3 Análisis bayesiano del problema inverso==

Una vez conocidas algunas bases sobre probabilidad y estadística, estadística bayesiana, se puede realizar inferencia bayesiana para resolver problemas inversos. El caso particular que se abordará, será el del modelo de crecimiento logístico. Para el desarrollo de este capítulo se utilizaron las siguientes bibliografías:  <span id='citeF-8'></span>[[#cite-8|[8]]], <span id='citeF-7'></span>[[#cite-7|[7]]], <span id='citeF-9'></span>[[#cite-9|[9]]], <span id='citeF-10'></span>[[#cite-10|[10]]] y <span id='citeF-11'></span>[[#cite-11|[11]]].

===3.1 Aproximación bayesiana a inferencia===

Se considera un fenómeno representado por medio de un sistema de ecuaciones diferenciales ordinarias.   Suponga que se cuenta con observaciones <math display="inline">\mathbf{y}=(y_1,\dots ,y_n)</math> tomadas en los tiempos <math display="inline">t=(t_1, \dots , t_n)</math>, con la siguiente estructura:

<span id="eq-8"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>y_i = \mathcal{H} (X_{\theta }(t_i)) + \varepsilon _i,   i=\overline{1,n} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (8)
|}

Como se mencionó en la introducción, el objetivo es hacer inferencia sobre <math display="inline">\theta </math> a partir de las observaciones. Para resolver el problema de inferencia ([[#eq-8|8]]) se propone que las entradas de <math display="inline">\theta </math> sean variables aleatorias que siguen alguna distribución de probabilidad.  La incertidumbre se cuantifica con una medida de probabilidad. El agente interesado en conocer el parámetro <math display="inline">\theta </math>, establece una variable aleatoria <math display="inline">\Theta </math> con su '''densidad de probabilidad'''

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>p_{\Theta }(\cdot ),</math>
|}
|}

Los valores que <math display="inline">\Theta </math> toma son los posibles valores para los parámetros, en este caso <math display="inline">\Phi =(\Theta ,\Sigma )</math> toma valores (<math display="inline">\theta ,\sigma </math>). Esta medida de probabilidad cuantifica la incertidumbre que tiene el agente respecto a los posibles valores de los parámetros en el modelo, <math display="inline">P_{\Phi }(\theta ,\sigma )</math> es la distribución a priori .

En presencia de datos <math display="inline">\mathbf{Y} = \mathbf{y}</math>, y suponiendo un modelo para la distribución conjunta:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>p_{\mathbf{Y}|\Phi }(\mathbf{y}|\theta ,\sigma ),</math>
|}
|}

para <math display="inline">\mathbf{Y}</math>, donde <math display="inline"> E(y_i \mid \theta , \sigma ) = \mathcal{H}(X_\theta (t_i)) </math>. Al observar los datos <math display="inline"> Y </math>, interesa inferir el valor de <math display="inline"> \theta </math>. La teoría bayesiana prescribe que calcular la distribución condicional de las incógnitas de interés dados los datos, se calcula utilizando el teorema de Bayes para variables aleatorias:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>p_{\Phi |\mathbf{Y}}(\theta ,\sigma | \mathbf{y}) = \frac{p_{\mathbf{Y}|\Phi }(\mathbf{y}|\theta ,\sigma ) p_{\Phi }(\theta ,\sigma )}{p_{\mathbf{Y}}(\mathbf{y})},</math>
|}
|}

<math display="inline">p_{\Phi |\mathbf{Y}}(\theta ,\sigma | \mathbf{y})</math> es la distribución ''a posteriori''; además <math display="inline">p_{\Phi }(\theta ,\sigma )</math> es la distribución ''a priori'' de <math display="inline">(\Theta ,\Sigma )</math> y

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>p_{\mathbf{Y}}(\mathbf{y}) =\int  p_{\mathbf{Y}, \Phi } (\mathbf{y} , \theta ,\sigma ) d \theta d \sigma , </math>
|}
|}

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>p_{\mathbf{Y}}(\mathbf{y}) =\int  p_{\mathbf{Y}| \Phi } (\mathbf{y} | \theta ,\sigma ) p_{\Phi }(\theta ,\sigma ) d \theta d \sigma , </math>
|}
|}

es la constante de normalización, también llamada '''verosimilitud marginal'''.

En el caso en que el error de cada observación representa un ruido aditivo gaussiano, la función de verosimilitud es:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math> p_{Y| \Phi } (\mathbf{y} | \theta , \sigma ) = \sigma ^{-n} (2\pi )^{-n/2} \exp \left(-\frac{1}{2\sigma ^2} \sum _{i=1}^{n} (y_i - \mathcal{H}(X_{\theta }(t_i)))^2 \right). </math>
|}
|}

Note que, cada vez que se evalué <math display="inline">p_{Y| \Phi } (\mathbf{y} | \theta , \sigma )</math> debe resolverse a <math display="inline">X_\theta </math>, lo cual se hace de forma aproximada por medio de un método numérico.

Como consecuencia,  el sistema de ecuaciones diferenciales ordinarias se resuelve utilizando un método numérico y la inferencia se realiza, no en el modelo exacto anterior, sino en un modelo aproximado, a saber:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math> y_i = \mathcal{H}(X_\theta ^h(t_i)) + \varepsilon _i, \quad \varepsilon _i \underset{i.i.d}{\sim } \mathcal{N}(0, \sigma ^2)  </math>
|}
|}

donde <math display="inline"> X_\theta ^h </math> denota la solución aproximada  proporcionada por el método numérico. La nueva verosimilitud derivada del modelo  es:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math> p^{h}_{\mathbf{Y}|\Phi }(\mathbf{y}|\theta , \sigma ) = \sigma ^{-n} (2\pi )^{-n/2} \exp \left(-\frac{1}{2\sigma ^2} \sum _{i=1}^{n} (y_i - \mathcal{H}(X_\theta ^h(t_i)))^2 \right). </math>
|}
|}

Para calcular <math display="inline">p_{\Phi |\mathbf{Y}} (\theta , \sigma | \mathbf{y} )</math> comúnmente se usan métodos tipo Monte Carlo vía cadenas de Markov (MCMC). Este cálculo a su vez va a estar afectado por la precisión del método numérico usado para calcular a <math display="inline">X_\theta{.}</math> 

La distribución ''a posteriori'' numérica es:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>p^{h}_{\Phi |Y} (\theta , \sigma | \mathbf{y} ) = \frac{p^{h}_{Y| \Phi } ( \mathbf{y}|\theta , \sigma ) p_{\Phi }(\theta , \sigma ) }{p^{h}_{\mathbf{Y}}(\mathbf{y})},     </math>
|}
|}

donde <math display="inline">p(\theta , \sigma )</math> es la distribución a priori en <math display="inline">(\theta , \sigma )</math> y

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>p^{h}_{\mathbf{Y}}(\mathbf{y})=\int  p^h_{\mathbf{Y}| \Phi } (\mathbf{y} | \theta , \sigma ) p_{\Phi }(\theta , \sigma ) d \theta d \sigma , </math>
|}
|}

es la constante de normalización, también llamada '''verosimilitud marginal'''. Note que dado que no hay otra alternativa que utilizar la distribución ''a posteriori'' numérica, hay una necesidad real de comprender y controlar el error incurrido al trabajar con <math display="inline">p_{\Phi |\mathbf{Y}} (\theta , \sigma | \mathbf{y} )</math> y la aproximada numéricamente <math display="inline">p^{h}_{\Phi |\mathbf{Y}} (\theta , \sigma | \mathbf{y} ).</math>

En el anexo A se menciona que un método numérico es de orden <math display="inline">p</math>, si <math display="inline">\varepsilon _h(\theta )=\mathcal{O}(h^p)</math>,  es decir, <math display="inline">\varepsilon _h(\theta ) \leq K h^p,</math> con K una constante global que no depende de <math display="inline">h</math>.  <span id='citeF-12'></span>[[#cite-12|[12]]] demuestran que bajo un tamaño de paso <math display="inline">h,</math> se garantiza que prácticamente no existe diferencia entre la distribución a posteriori teórica <math display="inline">p_{\Phi |\mathbf{Y}} (\theta , \sigma | \mathbf{y} )</math> y la aproximada numéricamente <math display="inline">p^{h}_{\Phi |\mathbf{Y}} (\theta , \sigma | \mathbf{y} ).</math> Con lo cual se le da validez al análisis de inferencia que se realice a partir de la distribución ''a posteriori'' aproximada numéricamente.

Los métodos estadísticos tradicionales se centran en la estimación puntual. El estimador de máxima ''a posteriori'' (MAP) se considera como una versión regularizada del estimador de máxima verosimilitud  o el estimador de mínimos cuadrados.

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\log (p^h_{\Phi |\mathbf{Y}}(\theta , \sigma |\mathbf{y}))= C + \log p^h_{\mathbf{Y}|\Phi }(\mathbf{y}|\theta , \sigma ) + \log p_{\Phi }(\theta , \sigma ) </math>
|-
| style="text-align: center;" | <math> = C - \frac{1}{2\sigma ^2} \sum _{i=1}^{n} (y_i - \mathcal{H}(X_\theta ^h(t_i)))^2 + \log p_{\Phi }(\theta , \sigma ), </math>
|}
|}

donde <math display="inline">C</math> es constante respecto a <math display="inline">\Phi </math>. Si se supone que hay independencia entre las variables <math display="inline">\theta </math> y <math display="inline">\sigma ,</math> se tiene que:

<span id="eq-9"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\log (p^h_{\mathbf{Y}|\Phi }(\mathbf{y}|\theta , \sigma ))=C - \frac{1}{2\sigma ^2} \sum _{i=1}^{n} (y_i - \mathcal{H}(X_\theta ^h(t_i)))^2 + \log p_{\Theta }(\theta )+\log p_{\Sigma }(\sigma ) </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (9)
|}

Esta última ecuación se llamará '''objetivo''', y se utilizará más adelante.

===3.2 El modelo logístico===

De acuerdo con <span id='citeF-9'></span><span id='citeF-13'></span>[[#cite-9|[9,13]]], una '''población''' es un grupo de organismos vivos (plantas, animales, microorganismos, etc.) que está compuesto por individuos con un comportamiento dinámico similar. Las poblaciones cambian de tamaño (crecen o disminuyen) debido al nacimiento, muerte y migración. 

La '''dinámica de poblaciones''' estudia las leyes que rigen los cambios de la población en el espacio y el tiempo. Se centra en cómo las poblaciones cambian con el tiempo. Además, una población se describe por su número de individuos.

En 1798, Thomas Robert Malthus, propuso su modelo bajo las siguientes hipótesis:

<ol>

<li>La población es homogénea (todos los individuos son idénticos).     </li>
<li>El medio es homogéneo, es decir, las características físicas, biológicas, etcétera, son las mismas en el hábitat.     </li>
<li>No hay limitaciones ni de espacio ni de alimento para el crecimiento de la población (la tasa de cambio de la población en el tiempo <math display="inline">t</math> es proporcional a la población en ese instante de tiempo).     </li>
<li>La población está aislada (no hay migración).     </li>
<li>Las tasas de natalidad y de mortalidad son constantes. </li>

</ol>

Denote por:

* <math display="inline">X(t)</math> el número de individuos en el tiempo <math display="inline">t</math>,
* <math display="inline">\beta </math> la tasa de natalidad, y
* <math display="inline">\mu </math> la tasa de mortalidad,

con <math display="inline">\beta </math> y <math display="inline">\mu </math> positivos, entonces, por las hipótesis dadas, se tiene el siguiente modelo:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math> \frac{dX(t)}{dt} = \beta X(t) - \mu X(t), </math>
|}
|}

o bien,

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math> \frac{dX(t)}{dt} = rX(t), </math>
|}
|}

donde <math display="inline">r = \beta - \mu > 0</math> y se le conoce como tasa de crecimiento instantáneo o tasa de crecimiento per cápita.

Cuarenta años más tarde, en 1838, el matemático belga Pierre François Verhulst (1804-1849) modificó el modelo de Malthus, cambiando la hipótesis 3:

* Los recursos (alimentos o tamaño del medio) son finitos.

Con esta nueva hipótesis habrá competencias entre la misma especie, así Verhulst propuso el modelo de crecimiento logístico:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math> \frac{dX(t)}{dt} = rX(t) \left(1 - \frac{X(t)}{K}\right),  </math>
|}
|}

o bien,

<span id="eq-10"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\frac{dX(t)}{dt} = L X(t) \left(K-X(t)\right),  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (10)
|}

donde <math display="inline"> L = \frac{r}{K} </math> y se llama '''capacidad de carga'''.

Dada una condición inicial <math display="inline"> X(0) = X_0 </math>, la ecuación ([[#eq-10|10]]) tiene como solución:

<span id="eq-11"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>X(t) = \frac{K X_0 e^{LKt}}{K + X_0\left(e^{LKt} - 1\right)}, </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (11)
|}

o bien,

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>X(t) = \frac{K}{1 + \left(\frac{K}{X_0} - 1\right)e^{-rt}}. </math>
|}
|}

A la ecuación ([[#eq-10|10]]) se le llama '''ecuación logística'''. A pesar de que el modelo de crecimiento logístico tiene a esta como solución explícita, se utilizará el programa ''odeint'', que se encuentra en la paquetería ''scipy.integrate'' implementado en Python; pues realizarlo de esta manera podrá ser replicable para cualquier otro modelo no lineal al que se pueda solucionar con un método numérico.

===3.3 Simulación estocástica y aplicación===

De acuerdo con <span id='citeF-8'></span>[[#cite-8|[8]]], las evaluaciones de incertidumbre también pueden adoptar la forma de intervalos o regiones de credibilidad, similares a los intervalos de confianza utilizados en la estadística clásica. En general, para realizar inferencias sobre las incógnitas en el modelo y finalmente responder preguntas relevantes de investigación, se debe ser capaces de analizar la distribución ''a posteriori''. 

Conforme a <span id='citeF-14'></span>[[#cite-14|[14]]], aunque no es prudente intentar analizar directamente las propiedades de la distribución posterior, existen métodos indirectos que puede proporcionarnos información considerable. En lugar de pensar en la distribución posterior como una función, se puede utilizar el hecho de que es una distribución de probabilidad y, por lo tanto, puede analizarse mediante métodos estadísticos, siempre que exista una manera de obtener una muestra. 

Para eso, se utilizará un método de '''Monte Carlo vía Cadenas de Markov''' (MCMC). Los métodos de MCMC representan un conjunto de algoritmos que permiten obtener muestras aleatorias de una determinada distribución de probabilidad objetivo, de la cual es difícil muestrear directamente. Estos métodos se basan en construir una cadena  de Markov cuya distribución de equilibrio, es la distribución objetivo. De esta forma, los estados de la cadena de Markov después de que esta ha alcanzado el estado estacionario representan muestras de la distribución objetivo. 

La principal ventaja de los métodos de Monte Carlo es que podemos muestrear de una medida de probabilidad solo conocida hasta una constante de normalización. La principal limitación de este enfoque es que los métodos de Monte Carlo se deterioran con el aumento de la dimensión del parámetro.

El algoritmo base para realizar MCMC, es el ''' algoritmo Metropolis-Hasting''', este tiene la siguiente estructura:


{| style="margin: 1em auto;border: 1px solid darkgray;"
|-
|
:Sea <math display="inline">f</math> la distribución de interés. Como pre-proceso se genera el valor incial <math display="inline">X_{(0)} \sim \mu _0</math> (con <math display="inline">\mu _0</math> una distribución de probabilidad con el mismo soporte que <math display="inline">f</math>). Para generar los valores <math display="inline">t=1,2, \dots </math> de la cadena de Markov, se hace: Siendo <math display="inline">X_{t}=x_t</math> el estado actual de la cadena de Markov se propone como candidato para el siguiente elemento de la cadena a <math display="inline">y_t \sim q(\cdot | x_t),</math> con <math display="inline">q</math> una distribución de probabilidad instrumental. Se hace <math display="inline"> X_{(t+1)} = \begin{cases} y_t & \hbox{ con probabilidad }\rho (x_t, y_t) \\ x_t & \hbox{ con probabilidad }1-\rho (x_t, y_t) \end{cases}</math> 

con <math display="inline">\rho (x_t, y_t) =\displaystyle \min \left\{1, \frac{f(y)q(x|y)}{f(x)q(y|x)}\right\}</math>. 


|-
| style="text-align: center; font-size: 75%;"|
<span id='algorithm-1'></span>'''Algoritmo. 1''' Algoritmo Metropolis -Hasting
|}

En general, los métodos de MCMC son muy complejos, requiriendo una calibración cuidadosa por parte de un experto, tanto para optimizar la velocidad de convergencia , como para identificar cuanto tiempo debe simularse la cadena antes de extraer la muestra. Una solución a este problema es el uso de una biblioteca de Python escrita por <span id='citeF-15'></span>[[#cite-15|[15]]], llamada el ''t-walk'', disponible en http://www.cimat.mx/&nbsp;jac/twalk/. El ''t-walk'' utiliza un tipo especial de algoritmo Metropolis-Hastings de propósito general que se ajusta automáticamente para muestrear de cualquier distribución cuando se le proporcionan las funciones de soporte y el logaritmo de la función objetivo, en este caso la ecuación ([[#eq-9|9]]). Para iniciar el MCMC se debe añadir la cantidad de muestras a obtener y dos puntos iniciales.

Simulando un conjunto de datos sintético con la ecuación ([[#eq-11|11]]), con el modelo de error <math display="inline"> y_i = X(t_i) + \varepsilon _i </math>, donde <math display="inline"> \varepsilon _i \sim N(0, \sigma ^2) </math>, y los siguientes parámetros:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math> X(0) = 100, \quad L = 1/1000, \quad K = 1000, \quad \sigma = 30. </math>
|}
|}

Se consideran 26 observaciones en los tiempos <math display="inline"> t_i </math> distribuidos regularmente entre 0 y 10.

Para realizar la cuantificación de incertidumbre para esta simulación, se tiene:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>    y_i = \mathcal{H} (X_{\theta }(t_i)) + \varepsilon _i,   i \in \{ 1,2, \dots ,n\} , </math>
|}
|}

donde <math display="inline">\mathcal{H}(x)=x</math> y <math display="inline"> \varepsilon _i \sim \mathcal{N}(0, \sigma ^2) </math>, además el modelo directo es,

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\frac{dX_{\theta }}{dt}= F(X_{\theta },t, \theta ), </math>
|}
|}

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>F(X_{\theta },t, \theta )=LX(K-X),</math>
|}
|}

donde <math display="inline">\theta = (L, K).</math> Suponiendo independencia ''a priori'' sobre los parámetros <math display="inline">K, L</math> y <math display="inline">\sigma </math>, las distribuciones ''a priori'' propuestas son:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>K \sim Gamma(2,1/0.001),</math>
|}
|}

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>L \sim Gamma(2,1/0.001),</math>
|}
|}

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\sigma \sim Gamma(2,1/0.001),</math>
|}
|}

dado que lo único que se “conoce” sobre estos parámetros es que son positivos, por esto se tomarán valores de la distribución gamma, la cual toma solo valores positivos.

<div id='img-1a'></div>
<div id='img-1b'></div>
<div id='img-1c'></div>
<div id='img-1'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Tinoco Guerrero_500938831-K_simulacion.png|600px|Distribución ''a posteriori'' del parámetro K.]]
|[[Image:Draft_Tinoco Guerrero_500938831-L_simulacion.png|600px|Distribución ''a posteriori'' del parámetro L.]]
|- style="text-align: center; font-size: 75%;"
| (a) Distribución ''a posteriori'' del parámetro K.
| (b) Distribución ''a posteriori'' del parámetro L.
|-
| colspan="2"|[[Image:Draft_Tinoco Guerrero_500938831-sigma_simulacion.png|600px|Distribución ''a posteriori'' del parámetro σ.]]
|- style="text-align: center; font-size: 75%;"
|  colspan="2" | (c) Distribución ''a posteriori'' del parámetro <math>\sigma </math>.
|- style="text-align: center; font-size: 75%;"
| colspan="2" | '''Figura 1:''' Distribuciones ''a posteriori'' de los parámetros del modelo.
|}

En las Figuras [[#img-1a|1a]], [[#img-1b|1b]] y [[#img-1c|1c]], se tienen los histogramas de <math display="inline">100\,000</math> valores de la distribución ''a posteriori'' de <math display="inline">K, L</math> y <math display="inline">\sigma{.}</math> Se puede observar que el valor verdadero, con el que se realizó la simulación está dentro de la distribución ''a posteriori'', además de estar cercano al punto con mayor credibilidad. Además, note que la distribución ''a priori'' queda muy por debajo de la distribución ''a posteriori''.

<div id='img-2'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Tinoco Guerrero_500938831-incertidumbreTesis.png|540px|Incertidumbre del ejemplo de simulación.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figura 2:''' Incertidumbre del ejemplo de simulación.
|}

En la Figura [[#img-2|2]] se muestran los datos simulados, la solución del modelo logístico con los parámetros reales, la curva que mejor se ajusta la cual se forma utilizando los valores de mayor credibilidad en cada distribución a posteriori. Note que la curva con los parámetros reales y el mejor ajuste está una encima de la otra, salvo un pequeño error que es provocado por el error del método numérico. Por último, la parte sombreada son curvas solución a partir de valores de las distribución ''a posteriori'' para cada parámetro.

===3.4 Aplicación: Saccharomyces cerevisiae===

Las levaduras son hongos que forman sobre los medios de cultivo colonias pastosas, constituidas en su mayor parte por células aisladas que suelen ser esféricas, ovoideas, elipsoideas o alargadas. Los genetistas <span id='citeF-16'></span>[[#cite-16|[16]]]  realizaron un estudio para comparar la capacidad de sobrevivencia de cepas haploides, diploides y tetraploides de la levadura ''Saccharomyces cerevisiae'' al ser expuestas a mutaciones inducidas por metanosulfonato de etilo, evaluando las ventajas del enmascaramiento en niveles altos de ploidía, así como la eliminación de mutaciones en células haploides. 

<span id='citeF-9'></span>[[#cite-9|[9]]] recopilaron los datos del crecimiento de esta población de levaduras en diferentes horas, los datos se muestran en la Tabla [[#table-1|1]].


{|  class="floating_tableSCP wikitable" style="text-align: center; margin: 1em auto;min-width:50%;"
|+ style="font-size: 75%;" |<span id='table-1'></span>Tabla. 1 Crecimiento de la levadura Saccharomyces cerevisiae.
|- style="border-top: 2px solid;"
|          '''Tiempo''' 
| style="border-right: 2px solid;" | '''Tamaño''' 
| style="border-left: 2px solid;" | '''Tiempo''' 
| '''Tamaño''' 
|-
| '''(hrs)''' 
| style="border-right: 2px solid;" | '''(<math>\times 10^6</math>/mL)''' 
| style="border-left: 2px solid;" | '''(hrs)''' 
| '''(<math>\times 10^6</math>/mL)''' 
|- style="border-top: 2px solid;"
|          0  
| style="border-right: 2px solid;" | 0.200 
| style="border-left: 2px solid;" | 19 
| 209   
|-
| 1  
| style="border-right: 2px solid;" | 0.330 
| style="border-left: 2px solid;" | 20 
| 190   
|-
| 2  
| style="border-right: 2px solid;" | 0.500 
| style="border-left: 2px solid;" | 21 
| 210   
|-
| 3  
| style="border-right: 2px solid;" | 1.10  
| style="border-left: 2px solid;" | 22 
| 200   
|-
| 4  
| style="border-right: 2px solid;" | 1.40  
| style="border-left: 2px solid;" | 23 
| 215   
|-
| 5  
| style="border-right: 2px solid;" | 3.10  
| style="border-left: 2px solid;" | 24 
| 220   
|-
| 6  
| style="border-right: 2px solid;" | 3.50  
| style="border-left: 2px solid;" | 25 
| 200   
|-
| 7  
| style="border-right: 2px solid;" | 9.00  
| style="border-left: 2px solid;" | 26 
| 180   
|-
| 8  
| style="border-right: 2px solid;" | 10.0  
| style="border-left: 2px solid;" | 27 
| 213   
|-
| 9  
| style="border-right: 2px solid;" | 25.4  
| style="border-left: 2px solid;" | 28 
| 210   
|-
| 10 
| style="border-right: 2px solid;" | 27.0  
| style="border-left: 2px solid;" | 29 
| 210   
|-
| 11 
| style="border-right: 2px solid;" | 55.0  
| style="border-left: 2px solid;" | 30 
| 220   
|-
| 12 
| style="border-right: 2px solid;" | 76.0  
| style="border-left: 2px solid;" | 31 
| 213   
|-
| 13 
| style="border-right: 2px solid;" | 115   
| style="border-left: 2px solid;" | 32 
| 200   
|-
| 14 
| style="border-right: 2px solid;" | 160   
| style="border-left: 2px solid;" | 33 
| 211   
|-
| 15 
| style="border-right: 2px solid;" | 162   
| style="border-left: 2px solid;" | 34 
| 200   
|-
| 16 
| style="border-right: 2px solid;" | 190   
| style="border-left: 2px solid;" | 35 
| 208   
|-
| 17 
| style="border-right: 2px solid;" | 193   
| style="border-left: 2px solid;" | 36 
| 230   
|- style="border-bottom: 2px solid;"
| 18 
| style="border-right: 2px solid;" | 190   
| style="border-left: 2px solid;" | 
| 

|}

En la figura [[#img-3a|3a]] se puede observar el crecimiento de la levadura durante las primeras 36 horas, donde es claro observar que siguen el comportamiento de un modelo de crecimiento logístico donde su capacidad de carga llega aproximadamente a 250.

<div id='img-3a'></div>
<div id='img-3b'></div>
<div id='img-3'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Tinoco Guerrero_500938831-datosCerveza.png|600px|Crecimiento de la levadura Saccharomyces cerevisiae.]]
|[[Image:Draft_Tinoco Guerrero_500938831-vazquez minimo.png|600px|Crecimiento de Saccharomyces cerevisiae, con el ajuste por mínimos cuadrados.]]
|- style="text-align: center; font-size: 75%;"
| (a) Crecimiento de la levadura Saccharomyces cerevisiae.
| (b) Crecimiento de Saccharomyces cerevisiae, con el ajuste por mínimos cuadrados.
|- style="text-align: center; font-size: 75%;"
| colspan="2" | '''Figura 3:''' Análisis del crecimiento de Saccharomyces cerevisiae.
|}

<span id='citeF-9'></span>[[#cite-9|[9]]] sugieren el modelo logístico y muestran que las estimaciones de los parámetros son:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>K=211.538461, \quad L=0.0026,</math>
|}
|}

como los mejores para el ajuste de los datos al modelo. Cabe mencionar que no especifican el método utilizado para dichas estimaciones. 

<span id='citeF-17'></span>[[#cite-17|[17]]] presentó las siguientes estimaciones de los mismos parámetros calculados por mínimos cuadrados:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>K=208.855224, \quad L=0.00263224.</math>
|}
|}

La figura [[#img-3b|3b]] muestra la gráfica del crecimiento y la curva con los parámetros estimados en <span id='citeF-17'></span>[[#cite-17|[17]]].

Usando estos datos, se realiza el procedimiento de cuantificación de la incertidumbre suponiendo el modelo de crecimiento logístico como el regresor. Este proceso se realiza de manera similar al ejemplo de simulación, es decir, se toman las mismas distribuciones ''a priori.'' <div id='img-4'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Tinoco Guerrero_500938831-K_cerveza_corregido.png|480px|Distribución ''a posteriori'' del parámetro K, ''Saccharomyces cerevisiae''.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figura 4:''' Distribución ''a posteriori'' del parámetro <math>K</math>, ''Saccharomyces cerevisiae''.
|}

<div id='img-5'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Tinoco Guerrero_500938831-L_cerveza_corregido.png|480px|Distribución ''a posteriori'' del parámetro L, ''Saccharomyces cerevisiae''.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figura 5:''' Distribución ''a posteriori'' del parámetro <math>L</math>, ''Saccharomyces cerevisiae''.
|}

<div id='img-6'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Tinoco Guerrero_500938831-Sigma_cervezaBien.png|480px|Distribución ''a posteriori'' del parámetro σ, ''Saccharomyces cerevisiae''.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figura 6:''' Distribución ''a posteriori'' del parámetro <math>\sigma </math>, ''Saccharomyces cerevisiae''.
|}
En las Figuras [[#img-4|4]], [[#img-5|5]] y [[#img-6|6]], se tienen los histogramas de <math display="inline">100\,000</math> valores de la distribución ''a posteriori'' de <math display="inline">K, L</math> y <math display="inline">\sigma{.}</math> Se puede observar que los valores estimados por <span id='citeF-9'></span>[[#cite-9|[9]]], <span id='citeF-17'></span>[[#cite-17|[17]]] están dentro de la distribución ''a posteriori'', además de estar cercanos al punto con mayor credibilidad, puede parecer que los valores estimados son los mismos pero se debe a la escala, estos si difieren un poco.  Además, note que la distribución ''a priori'' queda muy por debajo de la distribución ''a posteriori''.  Por último, se puede ver la distribución ''a priori'', donde hay un gran cambio entre la distribución ''a priori'' y la distribución ''a posteriori''.

<div id='img-7'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Tinoco Guerrero_500938831-IncertidumbreCerveza.png|600px|Incertidumbre de ''Saccharomyces cerevisiae''.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figura 7:''' Incertidumbre de ''Saccharomyces cerevisiae''.
|}
En la Figura [[#img-7|7]] se muestran los datos recopilados de la Tabla [[#table-1|1]], la solución del modelo logístico con los parámetros de <span id='citeF-17'></span>[[#cite-17|[17]]], <span id='citeF-9'></span>[[#cite-9|[9]]] y la curva que mejor se ajusta, esta última se realizó utilizando los valores de mayor credibilidad en cada distribución ''a posteriori''. Note que la curva con los parámetros dados por Vázquez se aleja un poco de las otras dos curvas, además la curva de mejor ajuste con la de Stewart y Day coinciden, por ese motivo se colocó de manera punteada, para distinguir entre cada una. 

La Figura [[#img-7|7]] es un ejemplo de una aplicación de la distribución ''a posteriori.'' Además, como se tiene varios valores de esta distribución se pueden encontrar valores de interés como las medidas de localización, medidas de variabilidad, coeficiente de asimetría y coeficiente de curtosis o algún otro cálculo dependiendo del interés del investigador.

El área sombreada se formó tomando curvas solución a partir de valores de las distribución ''a posteriori'' para cada parámetro. El uso de un área sombreada alrededor de los parámetros que mejor se ajustan en las gráficas de modelos es una técnica común en el análisis de datos. Esta área sombreada suele representar la incertidumbre o la variabilidad de los parámetros del modelo, en este caso, obtenida de la distribución ''a posteriori''. 

Esta área sombreada proporciona una visualización clara de la incertidumbre en las estimaciones de los parámetros del modelo. Muestra cómo varían las predicciones del modelo debido a la variabilidad en los datos. Esto es crucial en el enfoque bayesiano, donde la incertidumbre en los parámetros es una parte integral del análisis, por ejemplo, permite comunicar de manera efectiva la precisión y la confianza en las predicciones del modelo. El investigador puede ver,  '''no sólo una línea de mejor ajuste''', sino también cuánto pueden variar las predicciones. Esto es importante en contextos de toma de decisiones, donde entender la variabilidad puede influir en las decisiones basadas en los resultados del modelo. Además, facilita la comparación entre diferentes modelos o ajustes. Al superponer áreas sombreadas de diferentes modelos, se puede ver rápidamente cuál modelo proporciona predicciones más precisas o con menos incertidumbre. 

Las áreas sombreadas pueden ayudar a identificar puntos de datos que se encuentran fuera de las predicciones esperadas, señalando posibles anomalías o la necesidad de ajustar el modelo. Por todo esto, es mejor optar en el futuro por resolver problemas inversos utilizando estadística bayesiana, no sólo usar métodos estadísticos y numéricos clasicos.

==4 Conclusiones==

El análisis de cuantificación de incertidumbre bayesiana de problemas inversos continua siendo un tema de investigación desafiante. El presente artículo es una introducción a realizar inferencia en el aspecto del análisis bayesiano de sistemas de ecuaciones diferenciales ordinarias en el contexto de problemas inversos, en particular el caso del modelo logístico.

Los problemas inversos surgen en una variedad de aplicaciones científicas y de las ingenierías, donde los parámetros del modelo deben ser estimados a partir de datos observacionales. Estos problemas se caracterizan por errores observacionales, errores de modelo y problemas de mal planteamiento que generan incertidumbres en los parámetros del modelo. Los enfoques estadísticos bayesianos permiten realizar simulaciones y predicciones con incertidumbres cuantificadas.

Se ha destacado la dificultad añadida que presentan estos problemas debido a la incapacidad de tratar analíticamente la función del regresor, lo que nos obliga a recurrir a aproximaciones numéricas. A pesar de que comúnmente se pasa por alto la sustitución de la solución teórica por una aproximación numérica. Se ha señalado el trabajo de investigación que <span id='citeF-12'></span>[[#cite-12|[12]]] ha llevado sobre realizar inferencia sobre la distribución ''a posteriori'' númerica, en lugar de la téorica.

Como se mencionó, la distribución ''a posteriori'' puede tomar la forma de un área sombreada alrededor de los ajustes del modelo, esta área sombreada muestra la incertidumbre en las estimaciones de los parámetros del modelo. Muestra cómo varían las predicciones del modelo debido a la variabilidad en los datos. El investigador al resolver el problema inverso utilizando este método no solo obtiene los mejores parámetros que se ajustan, sino toda una distribución para ellos, lo que ayuda a ver cuánto pueden variar las predicciones.

Uno de los logros destacables de este trabajo es presentar de manera detallada un código realizado en Python, para la obtención de la distribución ''a posteriori'' para los parámetros de un modelo no lineal. Además, que el mismo genera histogramas de esta distribución, lo que ayuda a comprender a la misma y analogizar entre distribuciones más conocidas. Este código no lo muestran en trabajos del mismo tema, con este nivel de detalle. Los códigos se encuentran disponibles en https://github.com/RodGon22/RepositarioTesisRodrigo o escaneando el código QR de la Figura [[#img-8|8]].

Por último, la aplicación como la que se mostró en la sección 3 ilustra cómo un enfoque probabilístico para codificar errores en el proceso de modelado puede conducir a simulaciones predictivas con medidas de incertidumbre confiables en problemas del mundo real. Sin embargo, aplicaciones más rigurosas conllevan más desafíos relacionados con la computación de alto rendimiento, alta dimensionalidad en datos y parámetros, predicción y otros que no se abordaron en este trabajo que se pueden considerar problemas a futuro.

<div id='img-8'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Tinoco Guerrero_500938831-QRrepositario.png|300px|QR del repositario, generado con Google.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figura 8:''' QR del repositario, generado con Google.
|}

===BIBLIOGRAFÍA===

<div id="cite-1"></div>
'''[[#citeF-1|[1]]]''' Carpio, Ana. (2023) "Incertidumbre y problemas inversos"

<div id="cite-2"></div>
'''[[#citeF-2|[2]]]''' Molina, Muñoz J. D. (2023) "Cuantificación de la Incertidumbre desde el enfoque bayesiano en el contexto de EDO"

<div id="cite-3"></div>
'''[[#citeF-3|[3]]]''' Lee, Peter M. (2009) "Bayesian statistics: an introduction. 3rd". Wiley, New York

<div id="cite-4"></div>
'''[[#citeF-4|[4]]]''' Bernardo, José M and Smith, Adrian FM. (2009) "Bayesian theory", Volume 405. John Wiley & Sons

<div id="cite-5"></div>
'''[[#citeF-5|[5]]]''' Mendoza, Manuel and Regueiro, Pedro. (2011) "Estadística Bayesiana". Instituto Teconológico de México

<div id="cite-6"></div>
'''[[#citeF-6|[6]]]''' Correa Morales, Juan Carlos and Barrera Causil, Carlos Javier. (2018) "Introducción a la estadística Bayesiana". Editorial ITM. Textos Académicos

<div id="cite-7"></div>
'''[[#citeF-7|[7]]]''' Erdely, Arturo and Gutiérrez-Peña, Eduardo. (2023) "Monografía de Estadística Bayesiana". arXiv preprint arXiv:2309.06601

<div id="cite-8"></div>
'''[[#citeF-8|[8]]]''' Christen, J Andrés. (2020) "Cuantificación de Incertidumbre bayesiana (Bayesian UQ)"

<div id="cite-9"></div>
'''[[#citeF-9|[9]]]''' Stewart, James and Day, Troy. (2015) "Biocalculus: Calculus for Life Sciences". Cengage Learning

<div id="cite-10"></div>
'''[[#citeF-10|[10]]]''' Kuschinski Kathmann, Nicolás Elio. (2014) "Improvements in Chronology Building from 14C Measurements using Bayesian Inference on Autoregressive Gamma Processes". CIMAT

<div id="cite-11"></div>
'''[[#citeF-11|[11]]]''' Torres, Maria Luisa Daza. (2018) "Numerical Solution of the Inverse Scattering Problem using High Level Representations"

<div id="cite-12"></div>
'''[[#citeF-12|[12]]]''' Capistran, Marcos and Christen, J. and Daza, Maria and Flores Arguedas, Hugo and Montesinos-López, J. (2021) "Error Control of the Numerical Posterior with Bayes Factors in Bayesian Uncertainty Quantification", Volume -1. Bayesian Analysis

<div id="cite-13"></div>
'''[[#citeF-13|[13]]]''' Iannelli, Mimmo and Pugliese, Andrea. (2015) "An introduction to mathematical population dynamics: along the trail of volterra and lotka", Volume 79. Springer

<div id="cite-14"></div>
'''[[#citeF-14|[14]]]''' Kuschinski Kathmann Nicolás Elio. (2014) "Improvements in Chronology Building from 14C Measurements using Bayesian Inference on Autoregressive Gamma Processes". CIMAT

<div id="cite-15"></div>
'''[[#citeF-15|[15]]]''' Christen, J Andrés and Fox, Colin. (2010) "A general purpose sampling algorithm for continuous distributions (the t-walk)"

<div id="cite-16"></div>
'''[[#citeF-16|[16]]]''' Mable, Barbara K and Otto, Sarah P. (2001) "Masking and purging mutations following EMS treatment in haploid, diploid and tetraploid yeast (Saccharomyces cerevisiae)", Volume 77. Cambridge University Press. Genetics Research 1 9&#8211;26

<div id="cite-17"></div>
'''[[#citeF-17|[17]]]''' Virgilio Vázquez. (2023) "Notas de modelación matemática". Universidad Tecnológica de la Mixteca