MCM
Propagación de incertidumbres aplicando el método de Monte carlo
Aplicable a un modelo con cualquier número de magnitudes de entrada y una sola magnitud de salida, parte de un muestreo aleatorio de distribuciones de probabilidad. es una alternativa práctica al enfoque GUM, se aplica cuando:
El modelo es no lineal
La FDP para la magnitud de salida se aparta de una normal o de una t, debido a una marcada asimetría.
Es difícil proporcionar las derivadas parciales del modelo, tal como requiere la ley de propagación de la incertidumbre
Los modelos son arbitrariamente complicados
“El método MCM es aplicado para comprobar el marco de referencia GUM”
Comparación
Recorderis Demostración:
la incertidumbre tipo A proviene de la incertidumbre del promedio
\[\bar x=\frac{x_1+x_2+x_3+...+x_n}{n}\]
vamos a obtener su incertidumbre:
\[u(\bar x)=\sqrt{\left (\frac{1}{n}*u(x_1)\right)^2+\left (\frac{1}{n}*u(x_2)\right)^2+...+\left (\frac{1}{n}*u(x_n)\right)^2}\]
\[u(\bar x)=\sqrt{n \left (\frac{1}{n}*s\right)^2}\]
\[u(\bar x)=\sqrt{ \frac{n}{n^2}*s^2}\]
\[u(\bar x)=\sqrt{ \frac{n}{n^2}*s^2}\]
\[u(\bar x)= \frac{s}{\sqrt{n}}\]
¿Cuántas simulaciones son necesarias?
Según GUM S1 (2008) Un valor de \(M = 10^6\) suele proporcionar un intervalo de cobertura del 95 % para la magnitud de salida, de forma que la amplitud del intervalo es correcta con una o dos cifras decimales significativas.
\[ \begin{cases} N>>\frac{1}{1-p}, & \\ N>10^4\frac{1}{1-p}, \end{cases} \]
\(N>219781\) donde p es el área de cobertura usualmente p=95.45%
Simulación de las distribuciones de probabilidad
Distribución rectangular ó uniforme
En una distribución rectangular cada valor en un intervalo dado tiene la misma probabilidad, o sea la función de densidad de probabilidad es constante en este intervalo.
Sea \(X∼U(a,b)\), es decir, una variable aleatoria con distribución uniforme en el intervalo (a,b), con \(a,b \quad \varepsilon \quad \mathbb{R}\):
| Descripción | b y a con b>a | a y -a [centrada en cero] |
|---|---|---|
| fdp f(x) | \(\frac{1}{b-a}\) | \(\frac{1}{2a}\) |
| F(x) | \(\frac{x-a}{b-a}\) | \(\frac{x-a}{2a}\) |
| prom E(X) | \(\frac{a+b}{2}\) | 0 |
| \(E(x^2)\) | \(\frac{b^2+ba+a^2}{3}\) | \(\frac{a^2}{3}\) |
| v(x) | \(\frac{(b-a)^2}{12}\) | \(\frac{a^2}{3}\) |
| sd=u(x) | \(\frac{(b-a)}{\sqrt{12}}\) | \(\frac{a}{\sqrt3}\) |
El rango de una distribución de probabilidad se refiere al conjunto de todos los posibles valores que puede tomar una variable aleatoria en esa distribución. Es la diferencia entre el valor máximo y el valor mínimo de esta variable.
Simulación de la distribución uniforme
simulemos los resultados de una incertidumbre por resolución
\[u=\frac{a}{\sqrt 3}=\frac{1}{\sqrt 3}=0.577350 \] a es la semiamplitud=1, criterio necesario para definir la simulación
## [1] 0.5776757
Simulando la incertidumbre tipo A
como la desviación estándar proviene de una distribución uniforme hacemos uso de esta misma para la estimación de la incertidumbre.
\[u=\frac{s}{\sqrt {10}}=\frac{0.5773}{10}=0.182574\]
## [1] 0.1824945Distribución triangular
Generar una variable aleatoria continua con distribución triangular con parámetros a=-2 y b=2
Distribución normal
Distribución t de student
Simulación Distribución tipo A
sean los valores
## [1] 0.0003955015## [1] 0.0003071509Propagación de incertidumbres GUM
Propagación de distribuciones MCM
Ejemplo
considere la siguiente ecuación que describe al mensurando
\[y=f(x_1,x_2,x_3)=5x_1-3x_2+x_3 \]
\[x_1= uniforme\quad u_{x_1}=\frac{1}{\sqrt 3}=0.5773\] \[x_2=triangular= \quad u_{x_2}=\frac{1}{\sqrt 6}=0.4082\]
\[x_3=normal\quad u_{x_3}=1\]
Por la GUM
\[u_C^2(y)=\sum_{i=1}^N c_i^2u^2(x_i)\]
\[u(y)=\sqrt{(5*u(x_1))^2+(3*u(x_2))^2+(u(x_3))^2}=3.29\]
Por MCM
## [1] -0.0001099584## [1] 3.289513
## [1] -10.27849## [1] -10.27849## [1] 6.107311Estimación incertidumbre con la función propagate de la librería propagate
Algunos de los argumentos que pide la función propagate:
Matriz de covarianzas
Se utiliza cuando hay correlación entre variables
\[\begin{bmatrix} u_1^2(xi) & u(x_1)u(x_2)r(x_1,x_2)\\ u(x_2)u(x_1)r(x_2,x_1) & u_2^2(x_2)\\ \end{bmatrix}\]
Matriz de transformación R
Para incluir la correlación y desviación estándar de cada magnitud de entrada descomposiciòn de cholesky para a matriz de covarianza \[Cov(X)=R^TR\], esta matriz debe ser no definida positiva
Salidas de la función
Gradiente \[\nabla (f)=\left [\frac{\partial f}{\partial x_1 },...,\frac{\partial f}{\partial x_n} \right]\]
Hessiana \[\begin{bmatrix} \frac{\partial f}{\partial x_1^2} & \frac{\partial ^2f}{\partial x_1x_2} & \dots & \frac{\partial ^2f}{\partial x_1x_n}\\ \frac{\partial f}{\partial x_2 x_1} & \frac{\partial ^2f}{\partial x_2^2} & \dots & \frac{\partial ^2f}{\partial x_2x_n}\\ \vdots & \vdots & \dots & \vdots\\ \frac{\partial f}{\partial x_n x_1} & \frac{\partial ^2f}{\partial x_nx_2} & \dots & \frac{\partial ^2f}{\partial x_n^2}\\ \end{bmatrix}\]
Ejemplo
## Results from error propagation:## Mean.1 Mean.2 sd.1 sd.2 2.5% 97.5%
## 3.000000 3.000000 3.291129 3.291129 -3.450498 9.450498## Results from Monte Carlo simulation:## Mean sd Median MAD 2.5% 97.5%
## 2.999510 3.293733 2.999416 3.297018 -3.453876 9.460244## Welch-Satterthwaite degrees of freedom:## [1] 1614466## Coverage factor (k):## [1] 1.959965## Expanded uncertainty:## [1] 6.450498## Symbolic gradient matrix:## [1] "5" "-3" "1"## Evaluated gradient matrix (sensitivity):## [1] 5 -3 1## Symbolic hessian matrix:## [1] "0" "0" "0"
## [1] "0" "0" "0"
## [1] "0" "0" "0"## Evaluated hessian matrix:## [1] 0 0 0
## [1] 0 0 0
## [1] 0 0 0## Covariance matrix:## x1 x2 x3
## x1 0.3332753 0.0000000 0
## x2 0.0000000 0.1666272 0
## x3 0.0000000 0.0000000 1## Relative contribution:## x1 x2 x3
## x1 0.7692251 0.0000000 0.00000000
## x2 0.0000000 0.1384519 0.00000000
## x3 0.0000000 0.0000000 0.09232308## Skewness / Excess Kurtosis of MC evaluations:## 0.001934064 / -0.001485569## Shapiro-Wilk test for normality:## 0.7968585 => normal## Kolmogorov-Smirnov test for normality:## 0.3474236 => normal¿Cuál distribución de probabilidad se ajusta mejor a los datos?
AIC \[AIC = -2 * log(L) + 2 * k\] BIC \[BIC = -2 * log(L) + k * log(n)\]
L representa la probabilidad maximizada del modelo, que mide qué tan bien se ajusta el modelo a los datos.
k representa el número de parámetros en el modelo
log(n) representa el logaritmo del tamaño de la muestra (n)
## 1 of 32: Fitting Normal distribution...## .........## 2 of 32: Fitting Skewed-normal distribution...## .........10.........20.......## 3 of 32: Fitting Generalized normal distribution...## .........10.........20.......## 4 of 32: Fitting Log-normal distribution...## .........## 5 of 32: Fitting Scaled/shifted t- distribution...## .........10.........20.......## 6 of 32: Fitting Logistic distribution...## .........## 7 of 32: Fitting Uniform distribution...## .........## 8 of 32: Fitting Triangular distribution...## .........10.........20.......## 9 of 32: Fitting Trapezoidal distribution...## .........10.........20.........30.........40.........50
## .........60.........70.........80.## 10 of 32: Fitting Curvilinear Trapezoidal distribution...## .........10.........20.......## 11 of 32: Fitting Gamma distribution...## .........## 12 of 32: Fitting Inverse Gamma distribution...## .........## 13 of 32: Fitting Cauchy distribution...## .........## 14 of 32: Fitting Laplace distribution...## .........## 15 of 32: Fitting Gumbel distribution...## .........## 16 of 32: Fitting Johnson SU distribution...## .........10.........20.........30.........40.........50
## .........60.........70.........80.## 17 of 32: Fitting Johnson SB distribution...## .........10.........20.........30.........40.........50
## .........60.........70.........80.## 18 of 32: Fitting 3P Weibull distribution...## .........10.........20.......## 19 of 32: Fitting 2P Beta distribution...## .........## 20 of 32: Fitting 4P Beta distribution...## .........10.........20.........30.........40.........50
## .........60.........70.........80.## 21 of 32: Fitting Arcsine distribution...## .........## 22 of 32: Fitting von Mises distribution...## .........## 23 of 32: Fitting Inverse Gaussian distribution...## .........## 24 of 32: Fitting Generalized Extreme Value distribution...## .........10.........20.......## 25 of 32: Fitting Rayleigh distribution...## .........## 26 of 32: Fitting Chi-Square distribution...## ...## 27 of 32: Fitting Exponential distribution...## ...## 28 of 32: Fitting F- distribution...## .........## 29 of 32: Fitting Burr distribution...## ...## 30 of 32: Fitting Chi distribution...## ...## 31 of 32: Fitting Inverse Chi-Square distribution...## ...## 32 of 32: Fitting Cosine distribution...## .........
## Best fit is Normal Distribution.## Parameters:## mean sd
## 2.998587 3.297084## Standard errors:## mean sd
## 0.002884458 0.002355150## Goodness of fit:## BIC = -866.5849El ejercicio del parcial hecho en R
## Results from error propagation:## Mean.1 Mean.2 sd.1 sd.2 2.5% 97.5%
## 89.587580 89.587579 1.475181 1.475181 86.696274 92.478883## Results from Monte Carlo simulation:## Mean sd Median MAD 2.5% 97.5%
## 89.585321 1.474963 89.586887 1.476517 86.700010 92.472330## Welch-Satterthwaite degrees of freedom:## [1] 1000010## Coverage factor (k):## [1] 1.959966## Expanded uncertainty:## [1] 2.891304## Symbolic gradient matrix:## [1] "-(1/y * 100)" "x/y^2 * 100"## Evaluated gradient matrix (sensitivity):## [1] -0.062907325 0.006550175## Symbolic hessian matrix:## [1] "0" "1/y^2 * 100"
## [1] "1/y^2 * 100" "-(x * (2 * y)/(y^2)^2 * 100)"## Evaluated hessian matrix:## [1] 0.000000e+00 3.957332e-05
## [1] 3.957332e-05 -8.241080e-06## Covariance matrix:## x y
## x 549.9025 0.0000
## y 0.0000 0.2601## Relative contribution:## x y
## x 0.9999949 0.000000e+00
## y 0.0000000 5.128092e-06## Skewness / Excess Kurtosis of MC evaluations:## 0.0009975888 / -0.003091823## Shapiro-Wilk test for normality:## 0.6319615 => normal## Kolmogorov-Smirnov test for normality:## 0.6470225 => normal## 1 of 32: Fitting Normal distribution...## .........## 2 of 32: Fitting Skewed-normal distribution...## .........10.........20.......## 3 of 32: Fitting Generalized normal distribution...## .........10.........20.......## 4 of 32: Fitting Log-normal distribution...## .........## 5 of 32: Fitting Scaled/shifted t- distribution...## .........10.........20.......## 6 of 32: Fitting Logistic distribution...## .........## 7 of 32: Fitting Uniform distribution...## .........## 8 of 32: Fitting Triangular distribution...## .........10.........20.......## 9 of 32: Fitting Trapezoidal distribution...## .........10.........20.........30.........40.........50
## .........60.........70.........80.## 10 of 32: Fitting Curvilinear Trapezoidal distribution...## .........10.........20.......## 11 of 32: Fitting Gamma distribution...## .........## 12 of 32: Fitting Inverse Gamma distribution...## .........## 13 of 32: Fitting Cauchy distribution...## .........## 14 of 32: Fitting Laplace distribution...## .........## 15 of 32: Fitting Gumbel distribution...## .........## 16 of 32: Fitting Johnson SU distribution...## .........10.........20.........30.........40.........50
## .........60.........70.........80.## 17 of 32: Fitting Johnson SB distribution...## .........10.........20.........30.........40.........50
## .........60.........70.........80.## 18 of 32: Fitting 3P Weibull distribution...## .........10.........20.......## 19 of 32: Fitting 2P Beta distribution...## .........## 20 of 32: Fitting 4P Beta distribution...## .........10.........20.........30.........40.........50
## .........60.........70.........80.## 21 of 32: Fitting Arcsine distribution...## .........## 22 of 32: Fitting von Mises distribution...## .........## 23 of 32: Fitting Inverse Gaussian distribution...## .........## 24 of 32: Fitting Generalized Extreme Value distribution...## .........10.........20.......## 25 of 32: Fitting Rayleigh distribution...## .........## 26 of 32: Fitting Chi-Square distribution...## ...## 27 of 32: Fitting Exponential distribution...## ...## 28 of 32: Fitting F- distribution...## .........## 29 of 32: Fitting Burr distribution...## ...## 30 of 32: Fitting Chi distribution...## ...## 31 of 32: Fitting Inverse Chi-Square distribution...## ...## 32 of 32: Fitting Cosine distribution...## .........
## Best fit is Normal Distribution.## Parameters:## mean sd
## 89.584915 1.477004## Standard errors:## mean sd
## 0.001450904 0.001184658## Goodness of fit:## BIC = -835.3804Otro ejemplo
Calibración de un multímetro DMM
Punto de medición: 40V
Especificaciones del calibrador: \(\pm\) (18 ppm set+150\(\mu\) V)
| Lecturas DMM | 40.000 | 39.999 | 39.998 | 40.000 | 39.999 | 39.999 |
|---|---|---|---|---|---|---|
| Promedio | 39.999 | sd | 0.00075277 | n | 6 |
Resultados calibración del patrón (tensión continua)
| Intervalo V | 329.9999 |
|---|---|
| Valor indicado V | 40.0000 |
| Valor patrón V | 40.0020 |
| Error V | -0.0020 |
| Factor de cobertura k(95.45%) | 2 |
| Incertidumbre expandida V | 0.0002 |
| Incertidumbre típica | 0.0001 |
Se desea estimar la incertidumbre del error dada por:
\[Error= V_{medido} - V_{Patrón} \] \[E:(V_x+\delta V_{xres})-(V_{certif}+\delta V_{Esp})\]
Método GUM
| Componente | Valor | Tipo | Estima Incer | Valor incer | Distribución |
|---|---|---|---|---|---|
| \(V_x\) | Medido | A | \(u_A=\frac{0.000075277}{\sqrt 6}\) | 0.00030732 V | t student (gl=5) |
| \(\delta V_{res}\) | Resolución | B | \(u_{res}=\frac{0.001}{\sqrt {12}}\) | 0.00028868 V | Uniforme |
| \(V_s\) | Patrón | B | \(u_{pat}=\frac{0.0002}{2}\) | 0.0001 V | Normal |
| \(V_{esp}\) | Especificaciones | B | \(u_{esp}=\frac{0.00087}{\sqrt{3}}\) | 0.0005023 V | Uniforme |
\[u_c=\sqrt{0.00030732^2+0.00028868^2+0.0001^2+0.0005023^2}=0.0006634=6.6*10^{-4}\] ¿cumple con el criterio de la Incertidumbre dominante?
\[\frac{\sqrt{0.00030732^2+0.00028868^2+0.0001^2}}{0.0005023}= 0.86>0.3\] !!!no cumple!!!
\[v_{ef}=\frac{0.0006634^4}{\frac{0.00030732^4}{5}}=108\]
con un nc=0.9545 =95.45%
1-0.9545=0.455/2=0.02275
\(k=t(\alpha= 0.02275,gl=108)=2.023\)
U=2.023*0.0006634=0.001342
## [1] 2.023416## $ws.df
## [1] 108.4775
##
## $k
## [1] 2.273081
##
## $u.exp
## [1] 0.001507941Método Monte carlo
## [1] 0.0003962538## [1] 0.0002886565## [1] 9.990505e-05## [1] 0.0005021773## [1] 0.0007081529
## [1] 0.001371322Librería propagate
## Results from error propagation:## Mean.1 Mean.2 sd.1 sd.2 1.25%
## 2.0000000000 2.0000000000 0.0006633905 0.0006633905 1.9985130737
## 98.75%
## 2.0014869263## Results from Monte Carlo simulation:## Mean sd Median MAD 1.25%
## 1.9999995658 0.0006629561 1.9999950493 0.0006640144 1.9985239703
## 98.75%
## 2.0015044842## Welch-Satterthwaite degrees of freedom:## [1] 2432424## Coverage factor (k):## [1] 2.241404## Expanded uncertainty:## [1] 0.001486926## Symbolic gradient matrix:## [1] "1" "1" "-1" "1"## Evaluated gradient matrix (sensitivity):## [1] 1 1 -1 1## Symbolic hessian matrix:## [1] "0" "0" "0" "0"
## [1] "0" "0" "0" "0"
## [1] "0" "0" "0" "0"
## [1] "0" "0" "0" "0"## Evaluated hessian matrix:## [1] 0 0 0 0
## [1] 0 0 0 0
## [1] 0 0 0 0
## [1] 0 0 0 0## Covariance matrix:## Vm Vres Vcer VEsp
## Vm 9.444558e-08 0.000000e+00 0e+00 0.000000e+00
## Vres 0.000000e+00 8.333614e-08 0e+00 0.000000e+00
## Vcer 0.000000e+00 0.000000e+00 1e-08 0.000000e+00
## VEsp 0.000000e+00 0.000000e+00 0e+00 2.523053e-07## Relative contribution:## Vm Vres Vcer VEsp
## Vm 0.2146066 0.0000000 0.00000000 0.0000000
## Vres 0.0000000 0.1893629 0.00000000 0.0000000
## Vcer 0.0000000 0.0000000 0.02272278 0.0000000
## VEsp 0.0000000 0.0000000 0.00000000 0.5733077## Skewness / Excess Kurtosis of MC evaluations:## 0.02193534 / 0.0111011## Shapiro-Wilk test for normality:## 0.4835289 => normal## Kolmogorov-Smirnov test for normality:## 0.3438292 => normalComparación
| factor | u GUM | u MCM |
|---|---|---|
| u | 0.0006634 | 0.00071 |
| k | 2.023 | — |
| U | 0.001342 | 0.00137 |
Otro ejemplo teniendo en cuenta la correlación entre variables
\[y=f(x_1,x_2,x_3)=5x_1-3x_2+x_3\]
\[r(x_1,x_2)=-0.5\]
| variable | Distribución | promedio | \(u(x_i)\) | \(c_i\) |
|---|---|---|---|---|
| \(x_1\) | Normal | 1.5 | 1.2 | 5 |
| \(x_2\) | Normal | 3 | 0.8 | 3 |
| \(x_3\) | Uniforme | \(\bar x=0\) a=1 | \(u=1/\sqrt 3=0.5773\) | 1 |
Por el método GUM
\[u_c(y)=\sqrt{(c_1u_{x_1})^2+(c_2u_{x_2})^2+(c_3u_{x_3})^2+2c_1c_2 u_{x_1} u_{x_2} r(x_1,x_2))}\] \[u_c(y)=\sqrt{(5*1.2)^2+(3*0.8)^2+(1*0.5773)^2+(2*1.2*0.8*-0.5)}=7.51\]
sin correlación=6.48
Simulación de \(x_1\)
## [1] 0.9986677## [1] 0.9986677encontrar matriz de var-cov
\[cov=\begin{bmatrix} u_1^2(xi) & u(x_1)u(x_2)r(x_1,x_2)\\ u(x_2)u(x_1)r(x_2,x_1) & u_2^2(x_2)\\ \end{bmatrix}\]
\[cov=\begin{bmatrix} 1.2^2 & 1.2*0.8*-0.5\\ 1.2*0.8*-0.5 & 0.8^2\\ \end{bmatrix}\]
\[cov=\begin{bmatrix} 1.44 & -0.48\\ -0.48 & 0.64\\ \end{bmatrix}\]
\[x=\bar x+x_s*R\]
\[x_{i_{corr}}= \begin{bmatrix} \bar x_1 & \bar x_2 \\ \bar x_1 & \bar x_2\\ \bar x_1 & \bar x_2\\ \bar x_1 & \bar x_2\\ \bar x_1 & \bar x_2\\ \vdots & \vdots \\ \bar x_1 & \bar x_2 \end{bmatrix} +\begin{bmatrix} x_{1s} & x_{2s} \\ x_{1s} & x_{2s}\\ x_{1s} & x_{2s}\\ x_{1s} & x_{2s}\\ x_{1s} & x_{2s}\\ \vdots & \vdots \\ x_{1s} & x_{2s} \end{bmatrix}* \begin{bmatrix} 1.44 & -0.48\\ -0.48 & 0.64\\ \end{bmatrix}\]
## 2 x 2 Matrix of class "Cholesky"
## [,1] [,2]
## [1,] 1.2000000 -0.4000000
## [2,] . 0.6928203## Formal class 'Cholesky' [package "Matrix"] with 5 slots
## ..@ x : num [1:4] 1.2 0 -0.4 0.693
## ..@ Dim : int [1:2] 2 2
## ..@ Dimnames:List of 2
## .. ..$ : NULL
## .. ..$ : NULL
## ..@ uplo : chr "U"
## ..@ diag : chr "N"
## [1] -0.4991404## [1] 7.489189## [1] 7.489189## [1] 14.93451Bibliografía
GUM S1. (2008). JCGM 101:2008 Evaluación de datos de medición-Suplemento 1 de la “Guía para la expresión de la incertidumbre de medida”-Propagación de distribuciones aplicando el método de Monte Carlo. www.cem.es,