Curso de metrologı́a por internet

Unidad 3-1 Módulo 3

Introducción a la estimación de incertidum-

bres utilizando el Método de Monte-Carlo
1. Introducción
El objetivo de esta unidad es introducir al lector en la estimación de incertidumbres de medida uti-
lizando simulación estadı́stica (Método de Monte-Carlo) de acuerdo con lo descrito en el Suplemento
1 de la Guı́a ISO-GUM [1].

Dado que el Método de Mote-Carlo es un método numérico, se hace necesario disponer de un software
adecuado que permita realizar la simulación estadı́stica. Aunque dicha simulación se puede realizar
en lenguajes de programa de amplio uso como puedan ser BASIC, VisualBASIC, C, FORTRAN
o similares se recomienda utilizar softwares especı́ficamente diseñados para el cálculo numérico o
estadı́stico como puedan ser:
Matlab (www.mathworks.es). Software especifı́camente diseñado para trabajar con matrices y
vectores. Posee un lenguaje de programación de sintaxis sencilla y potentes herramientas de para
la simulación estadı́stica, permitiendo un facil manejo de vectores de dimensiones muy elevadas
(millones de elementos). Es un software comercial y el coste de las licencias es muy elevado. Sin
embargo, es habitual que las universidades y centros de investigación dispongan de licencias.
Octave (www.octave.org). Es un clon de Matlab distribuı́do bajo licencia GNU (software de
libre difusión).
R (www.r-project.org). Es un software libre para cálculo estadı́stico.
Los ejemplos que se van a presentar a lo largo de esta unidad han sido todos ellos resueltos utilizando
Octave (versión 3.6.2). Se recomienda que el lector no familiarizado con Octave (o Matlab) lea con
anterioridad algún manual de uso de este software como pueden ser:
Aprenda Matlab 7.0 como si estuviera en primero. J. Garcı́a de Jalón, J. Ignacio Rodrı́guez, J.
Vidal. ETS Ingenierios Industrales, Universidad Politécnica de Madrid.
http://mat21.etsii.upm.es/ayudainf/aprendainf/Matlab70/matlab70primero.pdf
GNU Octave Ed. 3. J.W. Eaton, D. Bateman, S. Hauber. Febrero 2011.
http://www.gnu.org/software/octave/octave.pdf
Matemáticas en Ingenierı́a con Matlab y Octave, GNU Octave Ed. 3. G. Borrell i Nogueras.
Marzo 2011.
http://iimyo.forja.rediris.es/iimyo2.pdf

c CEM
Curso Virtual de Metrologı́a: Unidad 3-1 v0.5 1
2. Simulacion de variables aleatorias
Simulacion de distribuciones uniformes

Casi todos los lenguajes de programación y hojas de cálculo disponen de un generador aleatorio de
números uniformemente distribuı́dos entre 0 y 1. Por ejemplo, en EXCEL esa función es aleatorio(),
en VisualBASIC es Rnd() y en Matlab y Octave es rand(). Como ejemplo, a continuación se genera
en Octave1 un vector x que contiene un millón de simulaciones de una variable aleatoria uniforme-
mente entre 0 y 1, dibujándose a continuación su histograma:

> x=rand(1e6,1);
> hist(x,400)

A continuación se muestra como se puede simular una variable aletoria uniformemente distribuı́da a
lo largo de un intervalo [a,b]:

> x=a+(b-a)*rand(1e6,1);
> hist(x,400)

Los histogramas obtenidos se muestran a continuación:

Figura 1

Simulacion de distribuciones triangulares

Una variable aletoria con función de densidad de probabilidades triangular simétrica se obtiene a
partir de la suma de dos variables aletorias uniformes iguales (con igual longitud de intervalo a lo
largo del cual se distribuyen uniformemente), tal y como puede comprarse con el siguiente ejemplo
en Octave2 :

> a=0; b=1;

> y=rand(1e6,1);
> z=rand(1e6,1);
> x=y+z;
> hist(x,400)

1 El prompt de Octave se ha modificado con la instrucción PS1(’>’) para hacerlo mas similar al que presenta Matlab
2 Para hacer visible la rejilla de lı́neas de puntos se debe ejecutar la instrucción grid on después de haberse dibujado el histograma

c CEM
Curso Virtual de Metrologı́a: Unidad 3-1 v0.5 2
 Figura 2

Si se desea generar una variable aleatoria x cuya distribución sea triangular simétrica de media m y
semiamplitud a ésta se puede obtener como suma de dos variables aleatorias y y z que se distribuyan
uniformemente entre 0 y 1 utilizando la siguiente expresión:

x = m + a · (y + z − 1) (1)
Como ejemplo se muestra la siguiente simulación realizada en Octave acompañada del histograma
resultante:

> m=15; a=5;

> y=rand(1e6,1);
> z=rand(1e6,1);
> x=m+a*(y+z-1);
> hist(x,400)

c CEM
Curso Virtual de Metrologı́a: Unidad 3-1 v0.5 3
 Figura 3

Si lo que se desea es simular una distribución triangular asimétrica, con un máximo en c y lı́mites infe-
rior y superior en a y b respectivamente, es posible realizarlo a partir de la simulación de dos variables
aleatorias y y z que se distribuyan uniformemente entre 0 y 1 utilizando la siguiente expresión:
√
x = c + [a − c + (b − a)z] · y (2)
El código Octave que a continuación se muestra realiza la simulación correspondiente a una dis-
tribución triangular asimétrica con lı́mites en 4 y 10 y con el máximo en 5 mostrando al final su
histograma:

> a=4; b=10; c=5;

> n=1e6;
> y=rand(n,1); z=rand(n,1);
> x=c+sqrt(y).*(a-c+z*(b-a));
> hist(x,400)

c CEM
Curso Virtual de Metrologı́a: Unidad 3-1 v0.5 4
 Figura 4

En el caso de una distribución triangular asimétrica como la anterior, no es inmediata la determi-

nación de la media de dicha distribución ni mucho menos su desviación tı́pica. Ahora bien, si hemos
almacenado las simulaciones en el vector x, entonces se muy simple el cálculo de la media y de la
desviación tı́pica utilizando las rutinas mean y std de Octave:

> m=mean(x)
m =
6.3339

> s=std(x)
s =
1.3123

Debe notarse que la media y la desviación tı́pica ası́ calculadas son la media muestral y desviación
tı́pica muestral de las n simulaciones generadas. Por tanto, se encontrarán cerca de la media µ y de
la desviación tı́pica σ de la población, pero no tienen por que coincidir. La diferencia entre ambas
será tanto mayor cuanto mayor sea el número n de simulaciones. Para comprobar lo anterior, basta
con volver a simular el vector x y calcular de nuevo la media y la desviación tı́pica. Se observa que
los valores obtenidos son ligeramente distintos:

> y=rand(n,1); z=rand(n,1);

> x=c+sqrt(y).*(a-c+z*(b-a));
> m=mean(x)
m =
6.3338

c CEM
Curso Virtual de Metrologı́a: Unidad 3-1 v0.5 5
 > s=std(x)
s =
1.3115

2.1. Ejercicio:
Se sabe que una determinada magnitud x se distribuye según una distribución
triangular asimétrica que alcanza su valor máximo en c= 18 y cuyos lı́mites infe-
rior y superior son a= 13 y b= 21 .

Se pide determinar su media x̄ y su desviación tı́pica s con un error inferior a un

5 %:

x̄ =
OK

s=
OK

Simulacion de distribuciones trapezoidales

Para obtener una variable aleatoria x que se distribuya trapezoidalmente de forma simétrica con media
m, de modo que la longitud de la base mayor del trapecio sea b + a y aquella de la base menor a − b
, basta con sumar dos variables aleatorias uniformemente distribuı́das entre 0 y 1 de acuerdo con la
siguiente expresión:

x = m + a · (y − 0, 5) + b · (z − 0, 5) (3)
A continuación se muestra un ejemplo al respecto en Octave:

> m=35; a=3; b=7;

> y=rand(1e6,1);
> z=rand(1e6,1);
> x=m+a*(y-0.5)+b*(z-0.5);
> hist(x,400)

c CEM
Curso Virtual de Metrologı́a: Unidad 3-1 v0.5 6
 Figura 5

2.2. Ejercicio:
Se sabe que una determinada magnitud x se distribuye según una distribución
triangular trapezoidal simétrica, de media m= 5 , longitud de la base mayor del
trapecio L= 20 y longitud de la base menor del trapecio l= 2 .

Se pide determinar su desviación tı́pica s con un error inferior al 5 %:

s=
OK

Simulacion de distribuciones normales

Una variable x que se distribuya normalmente (con media cero y varianza unidad) también puede ser
obtenida a partir de dos variables aleatorias y y z distribuı́das uniformemente entre 0 y 1 utilizando
para ello la transformación de Box-Müller [4,5]:

x = −2 ln y · cos(2πz)
p
(4)
El siguiente código de Matlab genera un vector con un millón de valores muestrales de una normal
(con media cero y varianza unidad) utilizando la transformación de Box-Müller. El resultado obtenido
(histograma) se muestra en la figura 6.

y=rand(1e6,1); % y es una uniforme en [0,1]

z=rand(1e6,1); % z es una uniforme en [0,1]

c CEM
Curso Virtual de Metrologı́a: Unidad 3-1 v0.5 7
x=sqrt(-2*log(y)).*cos(2*pi*z); % Transformacion de Box-Muller
hist(x,400); % Se dibuja el histograma

Figura 6

Para ver la bondad de la simulación obtenida es muy util recurrir a la función histfit 3 la cual dibuja
sobre el histograma la curva de la distribución normal que mejor se ajusta al vector de simulaciones x.
Utilizando dicha función ( histfit(x,400) se ha obtenido la figura 7, donde se puede observar un
ajuste bastante bueno entre la curva roja (distribución normal exacta) y el histograma. Si se aumenta
el número de simulaciones de un millón a diez millones, las diferencias entre ambos se reducen
significativamente (figura 8).
3 La función histfit pertenece a la toolbox Statistics de Matlab y al paquete (package) Statistics de Octave. Por tanto, para poder utilizar

dicha función es necesario haber instalado previamente la toolbox o el paquete (según el caso) Statistics.

c CEM
Curso Virtual de Metrologı́a: Unidad 3-1 v0.5 8
 Figura 7. Simulación de una distribución normal con N = 106 simulaciones

Figura 8. Simulación de una distribución normal con N = 107 simulaciones

De todos modos, en lenguajes de programación del tipo Octave/Matlab existen funciones (p.e función
randn) que permiten simular una normal de manera directa. El siguiente código genera un vector con
un millón de valores muestrales de una normal con media m y desviación tı́pica s:

c CEM
Curso Virtual de Metrologı́a: Unidad 3-1 v0.5 9
 > m=20; s=3;
> x=m+s*randn(1e6,1);
> mean(x)
ans =
19.9950
> std(x)
ans =
2.9967

Se observa que existe una pequeña diferencia entre la media y desviación tı́pica de los valores mues-
trales y aquellos de la distribución (20.0000 y 3.0000 respectivamente). Esta diferencia se reduce si
se aumenta el número de simulaciones a diez millones.

> m=20; s=3;

> x=m+s*randn(10e6,1);
> mean(x)
ans =
19.9995
> std(x)
ans =
3.0002

Simulacion de distribuciones en forma de U

Otra distribución estadı́stica que se usa en metrologı́a es la distribución en forma de U. Ésta se observa
en magnitudes de influencia (como la temperatura o la humedad) controladas analógicamente. Si
el sistema de control es de gran calidad y durante el proceso de medida no existen perturbaciones
externas es habitual observar comportamientos senoidales en las magnitudes controladas. En estas
situaciones, la magnitud controlada x sigue una distribución en forma de U. La simulación de x,
supuesto que oscila entre los valores mı́nimo xmin y máximo xmax es sencilla de obtener a partir de la
simulación de una variable aleatoria y que se distribuye uniformemente a lo largo del intervalo [0,1]
utilizando la siguiente expresión:
xmáx + xmı́n xmáx − xmı́n
x= + sin(2πy) (5)
2 2
A continuación se presenta una simulación con Matlab de una variable que se distribuye en forma de
U entre 19 y 21. El resultado se muestra en la figura 9.

> xmin=19; xmax=21;

> y=rand(1e6,1);
> x=(xmax+xmin)/2+(xmax-xmin)/2*sin(2*pi*y);
> hist(x,400);

c CEM
Curso Virtual de Metrologı́a: Unidad 3-1 v0.5 10
 Figura 9

Simulacion de una χ2

Una χ2 con ν grados de libertad es igual a la suma de los cuadrados de ν distribuciones normales,
todas ellas de media nula y varianza unidad:

χ2ν = x12 + x22 + · · · + xν2 (6)

El código Matlab/Octave que a continuación se muestra, permite simular una χ2 con cuatro grados de
libertad. El resultado se muestra en la figura 10:

> N=1e6;
> x1=randn(N,1); x2=randn(N,1); x3=randn(N,1); x4=randn(N,1);
> X2=x1.ˆ2+x2.ˆ2+x3.ˆ2+x4.ˆ2;
> hist(X2,400);

Realmente lo que se muestra en la figura 10 es el resultado de utilizar la función histfit(X2,400,’gamma’)4 ,

donde la rutina histfit ha ajustado la distribución gamma que mejor se ajusta al histograma obtenido
5
. La curva roja muestra la distribución chi2 exacta mientras que azul se muestra el histograma pro-
ducto de la simulación. Se puede observar el gran ajuste entre ambos.
4 Esta opción está disponible únicamente en Matlab. En el caso de Octave histfit únicamente es capaz de ajustar distribuciones normales
5 la distribución χ2 es un caso particular de la distribución gamma

c CEM
Curso Virtual de Metrologı́a: Unidad 3-1 v0.5 11
 Figura 10

2.3. Ejercicio:
Se pide determinar, con un error inferior al 5 %, la media x̄ y la desviación tı́pica
s de una distribución χ2 con ν= 5 grados de libertad:

x̄ =
OK

s=
OK

Simulacion de una t-Student

La distribución t-Student, con ν grados de liberdad, responde a la siguiente expresión:

x
t= q . (7)
χ2ν ν
donde χ2ν es una χ2 con ν grados de libertad y x es una normal de media nula y varianza unidad. El
código Matlab/Octave que se incluye a continuación muestra como simular una t con cuatro grados
de libertad. El histograma correspondiente al millón de simulaciones se presenta en la figura 11:

> N=1e6;
> x1=randn(N,1); x2=randn(N,1); x3=randn(N,1); x4=randn(N,1);
> X2=x1.ˆ2+x2.ˆ2+x3.ˆ2+x4.ˆ2;
> y=randn(N,1);

c CEM
Curso Virtual de Metrologı́a: Unidad 3-1 v0.5 12
 > t=y./sqrt(X2/4);
> hist(t,400);

De nuevo en la figura 11 lo que realmente se muestra en la figura 10 es el resultado de utilizar la

función histfit(t,400,’tlocationscale’) donde histfit ha ajustado una t-Student al his-
tograma obtenido.

Figura 11

No obstante, si se dispone del paquete (o de la toolbox) Statistics se puede utilizar la función

random para simular directamente tanto la t-Student como la chi2 . Por ejemplo, las siguientes in-
strucciones en Matlab/Octave permite simular una t-Student y una chi2 ambas con cuatro grados de
libertad:

> N=1e6;
> t=random(’t’,4,[N 1]);
> X2=random(’chi2’,4,[N 1]);

2.4. Ejercicio:
Se pide determinar, con un error inferior al 5 %, la media la desviación tı́pica s
de una distribución t-Student con ν= 2 grados de libertad:

s=
OK

c CEM
Curso Virtual de Metrologı́a: Unidad 3-1 v0.5 13
Simulacion de una distribución cualquiera

Finalmente queda por resolver el problema de simular una variable aletoria x que posea una distribu-
ción cualquiera. En este caso, si se conoce su función de distribución de probabilidad F(x) (por tano
su función de densidad de probabilidades serı́a f (x) = F 0 (x)), x se puede simular a partir de una
variabe aleatoria y que se distribuye uniformemente en el intervalo [0,1] utilizando la función inversa
F −1 (y) :

x = F −1 (y) (8)
Para demostrarlo basta que recordar que si la relación entre dos variables aletorias es y = G(x) la
relación entre sus funciones de densidad de probabilidades fy (y) y f x (x) serı́a:

fy (y)dy = f x (x)dx (9)

Es decir, la probabilidad de que x pertenezca al intervalo (x, x + dx) deber ser idéntica a la de que y
pertenzca al intervalo (y, y + dy) ya que ambas están relacionadas por la expresión y = G(x).

Si y se distribuye uniformemente en el intervalo [0,1], entonces fy (y) = 1. Además, dy/dx = G0 (x).

Por tanto:

f x (x) = G0 (x) (10)

Si se elige G(x) = F(x), entonces G0 (x) = F 0 (x) = f (x) y consecuentemente habrı́amos conseguido
que la función de densidad de probabilidades de x es la que inicialmente habı́amos pretendido que
tuviera.

Imaginemos que se desea simular una variable aleatoria x que se distribuye cosenoidalmente en el
intervalo [-1,1] de acuerdo con la siguiente función de densidad de probabilidades:
π πx
f (x) =cos (11)
4 2
La función de distribución de probabilidades es la integral siguiente:
π πx 1  πx
Z x Z x 
F(x) = f (x)dx = cos dx = sin +1 (12)
−1 −1 4 2 2 2
La función inversa x = F −1 (y) serı́a x = π2 arcsin(2y − 1).

El siguiente código de Matlab permite simular dicha variable aleatoria a partir de la generación de
números aleatorios distribuı́dos uniformemente entre 0 y 1 (variable y). En la figura 12 se muestra el
histograma obtenido. Podrı́a comprobarse sin mucha dificultad que el histograma se ajusta bastante
bien a la función de densidad de probabilidades f (x) = π4 cos πx
2
.

N=1e6; % Número de simulaciones

y=rand(N,1); % y es una uniforme en [0,1]
x=2/pi*asin(2*y-1); % Aplicacion de la funcion inversa
hist(x,400); % Dibujo del histograma

c CEM
Curso Virtual de Metrologı́a: Unidad 3-1 v0.5 14
 Figura 12

3. Elección del número N de simulaciones

A lo largo del apartado anterior se ha observado que las estimaciones que proporciona el Método de
Monte-Carlo tanto de la media como de la desviación tı́pica de una distribución dada están afectadas
de una cierta dispersión. Para centrar algo mas el problema, considérese el siguiente problema:

Se desea estimar el valor medio y la desviación tı́pica de una variable aleatoria que se distribuye
normalmente con valor medio verdadero µ = 175, 0 y con desviación tı́pica verdadera σ = 9, 5. Para
ello se recurre al Método de Monte-Carlo utilizando N simulaciones.

A continuación se muestra el código Matlab/Octave utilizado para ello utilizando un millón de simu-
laciones:

mu=175.0; sigma=9.5;
N=1e6; % Número de simulaciones
x=mu+sigma*randn(N,1); % Simulación de la variable normal
mx=mean(x) % Estimación de su valor medio
sx=std(x) % Estimación de su desviación tı́pica

La ejecución de dicho código proporciona los siguientes resultados:

mx =
174.9919
sx =
9.4958

Pero si se repite la ejecución, se obtienen valores ligeramente diferentes:

c CEM
Curso Virtual de Metrologı́a: Unidad 3-1 v0.5 15
 mx =
175.0094
sx =
9.4994

Ello es debido a que x̄ (mx) es el el valor medio de N valores muestrales de la variable aleatoria x
(normal de media µ y desviación tı́pica σ). Por tanto, x̄ no es una constante
√ real (como lo son µ y σ),
sino que es una variable aleatoria de media µ y desviación tı́pica σ/ N. Debido a ello, cada nueva
ejecución del código conduce a un nuevo valor muestral de x̄ (mx).

Lo mismo ocurre con s(x) (sx). s(x) es un valor muestral de la siguiente variable aleatoria:
s
χ2N−1
σ· (13)
N−1
donde χ2N−1 es una χ2 con N − 1 grados
. √ de libertad. El valor medio de s(x) es σ y su desviación tipica
es muy aproximadamente igual a σ 2N.

A continuación se muestra un código Matlab/Octave que permite repetir la simulación anterior diez
veces con el fin de visualizar la variabilidad de x̄ (mx) y s(x) (sx) entre simulaciones:

c CEM
Curso Virtual de Metrologı́a: Unidad 3-1 v0.5 16
mx =[]; sx =[];
for i =1:10
mu =175.0; sigma =9.5;
N=1e6; % Número de simulaciones
x=mu+sigma* randn (N ,1); % Simulación de la variable normal
mx(i)=mean(x); % Estimación de su valor medio
sx(i)=std(x); % Estimación de su desviación tı́pica
disp( sprintf (’ mx = %8.4f sx = %7.4f’,mx(i),sx(i)));
end
disp(’ ’);
disp( sprintf (’s(mx) = %8.4f s(sx) = %7.4f’,std(mx),std(sx)));

Este código, además de repetir diez veces la simulación, almacena los resultados obtenidos y final-
mente calcula las desviaciones tı́picas de las series de diez valores de x̄ (mx) y s(x) (sx). Los resultados
que se obtienen se muestran a continuación:

mx = 174.9982 sx = 9.5013
mx = 174.9846 sx = 9.5013
mx = 174.9907 sx = 9.4993
mx = 174.9813 sx = 9.4945
mx = 174.9941 sx = 9.5064
mx = 175.0036 sx = 9.4927
mx = 174.9853 sx = 9.4905
mx = 174.9865 sx = 9.5043
mx = 174.9919 sx = 9.5035
mx = 174.9923 sx = 9.4967

s(mx) = 0.0067 s(sx) = 0.0053

A continuación se muestra un código Matlab/Octave que repite 40 veces cada simulación para un
número variable de simulaciones N comprendido entre N = 103 y N = 5 × 106 . La variabilidades
(desviaciones tı́picas) observadas en x̄ (mx) y s(x) (sx) se muestran en la figura 13. Los cı́rculos rojos
muestran las variabilidades observadas en x̄ (mx) y los cı́rculos azules las variabilidades
√ observadas
√ teórica σ[ x̄]/σ = 1/ N. La lı́nea azul
en s(x) (sx). La lı́nea roja de puntos respresenta la ecuación
de puntos representa la ecuación teórica σ[s(x)]/σ = 1/ 2N.

c CEM
Curso Virtual de Metrologı́a: Unidad 3-1 v0.5 17
Nj =1 e3*[ 1 2 5 10 20 50 100 200 500 1000 2000 5000];

% Datos de la variable aleatoria a simular

mu =175.0; sigma =9.5;

s_mx= zeros ( length (Nj) ,1);

s_sx= zeros ( length (Nj) ,1);
Nrep =40; % Número de repeticiones

for j=1: length (Nj)

N=Nj(j);
mx= zeros (Nrep ,1);
sx= zeros (Nrep ,1);
for i=1: Nrep
x=mu+sigma* randn (N ,1); % Simulación de la variable normal
mx(i)=mean(x); % Estimación de su valor medio
sx(i)=std(x); % Estimación de su desviación tı́pica
end
s_mx(j)=std(mx); % Variabilidad ( desviación tı́pica ) de mx
s_sx(j)=std(sx); % Variabilidad ( desviación tı́pica ) de sx
end

% Dibujo de los resultados obtenidos

loglog (Nj ,s_mx/sigma ,’ro’,Nj ,s_sx/sigma ,’bo’);
hold on
loglog (Nj ,1./ sqrt(Nj),’r:’,Nj ,1./ sqrt (2* Nj),’b:’);
grid on
xlabel (’Número de Simulaciones N’);
ylabel (’s(mx)/\ sigma y s(sx)/\ sigma ’);

Si se desearan obtener resultados redondeados a dos cifras significativas en la desviación tı́pica s(x)
habria que asegurar que las variabilidades observadas son un orden de magnitud inferior. En el ejem-
plo, redondear a dos cifras significativas en la desviación tı́pica s(x) es equivalente a mantener una
sola cifra decimal y, por tanto, habria que asegurar que la variabilidad observada en x̄ y en s(x) es
inferior a 0.01. En este caso, dado que sigma = 0,95, se tendrı́a:
σ( x̄) 0, 01 σ[s(x)] 0, 01
≤ = 1, 1 × 10−3 ≤ = 1, 1 × 10−3 (14)
σ 9, 5 σ 9, 5
Para conseguir ambas cosas simultáneamente, el número de simulaciones deberia ser del orden de un
millón (ver figura 13).

El el caso de que el valor de σ hubiera sido del tipo σ = 11, hubiera bastado redondear al entero
más próximo para mantener dos cifras significativas en la desviación tı́pica. Por ello, también habrı́a
bastado con limitar la variabilidad en en x̄ y en s(x) a 0,1. En ese caso, se hubiera tenido:
σ( x̄) 0, 1 σ[s(x)] 0, 1
≤ = 9, 1 × 10−3 ≤ = 9, 1 × 10−3 (15)
σ 11 σ 11
Ahora, para conseguir ambas cosas simultáneamente, el número de simulaciones deberia ser del or-
den de 12 × 103 simulaciones (ver figura 13). Ası́ pues, en el caso de distribuciones normales, el
número de simulaciones requeridas (pensando en el objetivo de mantener dos cifras significativas en
la estimación de la desviación tı́pica) variarı́a entre N = 12 × 103 y N = 106 .

c CEM
Curso Virtual de Metrologı́a: Unidad 3-1 v0.5 18
 Figura 13

En la página siguiente se muestra un estudio similar al anterior pero donde la variable aleatoria se
distribuye en forma de U en vez de normalmente. En la figura 14 puede observarse que la evolución de
la variabilidad de la estimación valor
√ medio x̄ (mx) sigue la misma pauta que la ya vista anteriormente
(lı́nea roja de puntos, σ[ x̄]/σ = 1/ N). Sin embargo, ahora la variabilidad observada en s(x) (s(x))
es claramente inferior. En el supuesto de que σ = 9, 5 y que el único objetivo fuera estimar dicha
desviación tı́pica con dos cifras significativas, bastarı́a con asegurar:
σ[s(x)] 0, 01
≤ = 1, 1 × 10−3 (16)
σ 9, 5
lo cual se lograrı́a con tan solo N = 2 × 105 simulaciones. En el caso en el cual σ = 11, bastarı́a con
asegurar:
σ[s(x)] 0, 1
≤ = 9, 1 × 10−3 (17)
σ 11
y para ello serı́an suficientes N = 1000 simulaciones.

3.1. Ejercicio:
Se pide determinar el número N mı́nimo de simulaciones para poder determinar
con dos cifras significativas, mediante simulación, la desviación tı́pica de una
distribución normal sabiendo que dicha desviación tı́pica se encuentra cercana
a 5,3.

s=
OK

c CEM
Curso Virtual de Metrologı́a: Unidad 3-1 v0.5 19
Nj =1 e3*[ 1 2 5 10 20 50 100 200 500 1000 2000 5000];

% Datos de la variable aleatoria a simular

mu =175.0; A =13.4; sigma=A/sqrt (2);

s_mx= zeros ( length (Nj) ,1);

s_sx= zeros ( length (Nj) ,1);
Nrep =40; % Número de repeticiones

for j=1: length (Nj)

N=Nj(j);
mx= zeros (Nrep ,1);
sx= zeros (Nrep ,1);
for i=1: Nrep
x=mu+A*sin (2* pi*rand(N ,1)); % Simulación de la variable en
forma de U
mx(i)=mean(x); % Estimación de su valor medio
sx(i)=std(x); % Estimación de su desviación
tı́pica
end
s_mx(j)=std(mx);
s_sx(j)=std(sx);
end

loglog (Nj ,s_mx/sigma ,’ro’,Nj ,s_sx/sigma ,’bo’,Nj ,1./ sqrt(Nj),’r:’);

grid on
xlabel (’Número de Simulaciones N’);
ylabel (’s(mx)/\ sigma y s(sx)/\ sigma ’);

c CEM
Curso Virtual de Metrologı́a: Unidad 3-1 v0.5 20
 Figura 14

3.2. Ejercicio:
Se pide determinar el número N mı́nimo de simulaciones para poder determinar
con dos cifras significativas, mediante simulación, la desviación tı́pica de una
distribución en forma de U sabiendo que dicha desviación tı́pica se encuentra
cercana a 3,3.

s=
OK

3.3. Ejercicio:
Se pide determinar el número N mı́nimo de simulaciones para poder determi-
nar con dos cifras significativas, mediante simulación, la desviación tı́pica de
una distribución uniforme. Se sabe que los valores mı́nimo y máximo de dicha
distribución son, aproximadamente, a− = 3 y a+ = 4.

s=
OK

c CEM
Curso Virtual de Metrologı́a: Unidad 3-1 v0.5 21
4. Propagación de Distribuciones mediante Simulación
Supóngase que la función modelo que caracteriza un procedimiento de medida, de acuerdo con la
guı́a ISO-GUM [1], es la siguiente:

y = f (x1 , x2 , · · · , xm ) (18)
Donde la variable de salida y representa el resultado de la medida y las variables de entrada x1 , x2 , · · · , xm
representan las lecturas directas obtenidas en los instrumentos de medida, valores certificados o cor-
recciones de calibración, otras correcciones, estimaciones de las magnitudes de influencia, etc...

Para propagar las funciones de distribución de las variables de entrada con el fin de obtener la función
de distribución de la variable de salida basta con simular un número elevado N de veces (del orden de
un millón) las variables de entrada, evaluar para cada simulación la variable de salida (se tendrı́an N
simulaciones) y dibujar el histograma de y utilizando las N simulaciones disponibles.

Dado que todo el proceso se realiza numéricamente, no es necesario que la función modelo y =
f (x1 , x2 , · · · , xm ) sea estrictamente una función analı́tica (como era necesario al aplicar la ley de propa-
gación de las incertidumbres utilizando el procedimiento descrito en la guı́a ISO-GUM). En este caso
batarı́a con que la función modelo f fuera un algoritmo susceptible de ser programable en algún
lenguaje de programación.

Como ejemplo, considerese que y es la suma de diez variables de entrada x1 , x2 , · · · , x10 que se dis-
tribuyen uniformemente entre 0 y 1. La función modelo serı́a:

y = x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8 + x9 + x10 (19)
La distribución de y podrı́a obtenerse por simulación utilizando el siguiente cógido de Matlab/Octave:

N=1 e6; % Numero de simulaciones

x=zeros (N ,10); % Se crea una matriz nx10
for i =1:10
x(:,i)=rand(N ,1); % Se simulan las xi n veces como
end % uniformes en [0 ,1]
y=0;
for i =1:10
y=y+x(:,i); % Se calcula y como suma de las xi
end
hist(y ,400); grid on; % Se dibuja el histograma
my=mean(y) % Media
sy=std(y) % Desviacion tipica

La ejecución del código anterior proporciona los siguientes resultados numéricos junto con el his-
tograma que se muestra en la figura 15:

my =
4.9998
sy =
0.9125

c CEM
Curso Virtual de Metrologı́a: Unidad 3-1 v0.5 22
 Figura 15

Puede observarse que la distribución resultante (figura 15) se parece mucho a una normal aun cuando
las distribuciones de partida sean uniformes. Se comprueba en este caso que el Teorema del Lı́mite
Central [5] casi se cumple con tan solo diez distribuciones.

La media de las N simulaciones de y coincide (no exactamente porque se trata de una simulación)
con la suma de las medias de las diez distribuciones (10 × 0, 5 = 5). Igualmente, la desviación tı́pica
de las N simulaciones de y coincide casi exactamente con la suma cuadrática de las√diez desviaciones
tı́picas de las x1 , x2 , · · · , x10 (la desviación tı́pica de una uniforme entre 0 y 1 es 1/ 12 ) :
v
t 10
X 1 !2
s = 0, 9125  √ = 0, 9129 (20)
i=1 12
Supóngase ahora que el resultado y de una medida es la diferencia entre el máximo y el mı́nimo de
diez lecturas x1 , x2 , · · · , x10 obtenidas con un instrumento. Supóngase que esas diez lecturas tienen
media nula y desviación tı́pica unidad. La función modelo podrı́a expresarse como:

y = máx {x1 , x2 , · · · , x10 } − mı́n {x1 , x2 , · · · , x10 } (21)

Ahora bien, la función anterior no es derivable y por tanto no podrı́a aplicarse el procedimiento de
estimación de incertidumbres descrito en la guı́a ISO-GUM. Sin embargo, es muy sencillo propagar
las distribuciones de x1 , x2 , · · · , x10 para obtener la distribución de y. En concreto, podrı́a utilizarse el
siguiente código de Matlab/Octave para ello:

N=1 e6; % Numero de simulaciones

x=zeros (10,N); % Se reserva espacio ( matriz de ceros)
% para una matriz de dimensión Nx10
for i =1:10
x(i ,:)= randn (1,N); % Se simulan las xi N veces como

c CEM
Curso Virtual de Metrologı́a: Unidad 3-1 v0.5 23
end % normales de media nula y s=1
y=max(x)-min(x); % y como max - min
hist(y ,400); grid on; % Se dibuja el histograma
mx=mean(y) % Media
sy=std(y) % Desviacion tipica

La ejecución de este código produce los resultados numéricos que se muestran a continuación y el
histograma representado en la figura 16. Nótese la asimetrı́a de la distribución resultante. El valor
mı́nimo de y es cero no estando acotado el valor máximo que puede alcanzar.

mx =
3.0782
sy =
0.7974

Figura 16

c CEM
Curso Virtual de Metrologı́a: Unidad 3-1 v0.5 24
 4.1. Ejercicio:
La función modelo que describe a un mensurando es la siguiente:

y = x1 · x2 (22)
Se sabe que x1 se distribuye uniformemente entre −a y +a (a= 7 ). Asimismo,
x2 se distribuye uniformemente entre −b y +b (b= 11 ).

Se pide determinar con dos cifras significativas la incertidumbre tı́pica de y.

u(y) =
OK

4.2. Ejercicio:
La función modelo que describe a un mensurando es la siguiente:

y = x12 (23)
Se sabe que x1 se distribuye en forma de U entre −a y +a (a= 1 ).

Se pide estimar con dos cifras significativas el valor de y y su incertidumbre tı́pica

u(y).

y=
OK

u(y) =
OK

4.3. Ejercicio:
La función modelo que describe a una corrección de calibración es la siguiente:

y = 1 − cos x1 (24)
Se sabe que x1 se distribuye uniformemente entre −a y +a (a= 0.01 rad).

Se pide estimar con dos cifras significativas el valor de y y su incertidumbre tı́pica

u(y).

y=
OK

u(y) =
OK

c CEM
Curso Virtual de Metrologı́a: Unidad 3-1 v0.5 25
5. Referencias
[1] JCGM/WG 1: JCGM 100:2008 Evaluation of measurement data - Guide to the expression of
uncertainty in measurement, GUM 1995 with minor corrections, First edition 2008, corrected version
2010, 14+120 p. Existe traducción al español de la edición en inglés de 2008, realizada por el CEM
y publicada como edición digital 1 en español (3a edición en español 2009), NIPO 706-10-001-0,
12+130 págs., accesible por Internet en la página web del CEM.
[2] JCGM/WG 1: JCGM 101:2008 Evaluation of measurement data - Supplement 1 to the Guide
to the expression of uncertainty in measurement - Propagation of distributions using a Monte Carlo
method, First edition 2008, 8+82 p. Existe traducción al español de la edición en inglés de 2010,
realizada por el CEM y publicada como edición digital (1a edición en español 2010), NIPO 706-10-
014-9, 7+99 págs., accesible por Internet en la página web del CEM.
[3] Cox, M.G.; Harris, P.M.: Uncertainty Evaluation, Software Support for Metrology, Best Practice
Guide No. 6. NPL, September 2006, 175 p. Disponible digitalmente a través de la página web http:
//www.npl.co.uk
[4] Wikipedia: Box-Muller Transform. http://en.wikipedia.org/wiki/Box-Muller_transform
[5] Box, G.E.P.; Muller, M. E.: A Note on the Generation of Random Normal Deviates. The Annals
of Mathematical Statistics (1958), Vol. 29, No. 2 pp. 610-611.

Califica

Curso elaborado para el CEM por el LMM-ETSII-UPM

a partir de los textos preparados por los profesores
que se relacionan al principio de cada Módulo.

c Centro Español de Metrologı́a

NIPO: 074-12-016-X
Se prohibe la reproducción total o parcial de este documento, cualquiera que sea el medio o tecnologı́a que se
utilice, sin permiso escrito del Centro Español de metrologı́a. Como excepción se autorizan:
1. La reproducción en papel para uso personal de los estudiantes registrados.
2. Las citas breves, siempre con expresión de la fuente, en publicaciones divulgativas, docentes, cientı́ficas o
profesionales.

26