Download Introduccion a la Simulacion y a la Generacion de Datos

Introduccion a la Simulacion y a la
Generacion de Datos Pseudoaleatorios
Departamento de Computación
Facultad de Ciencias Exactas y Naturales
Universidad de Buenos Aires
Taller de Informática I - 1er Cuatrimestre 2014
Andrea Manna
Definición
 Es el proceso de diseñar un modelo de un
sistema real, que servirá para dirigir
experimentos con el propósito de
entender el comportamiento del sistema
o evaluar nuevas estrategias – dentro de
los límites impuestos por un cierto
criterio o un conjunto de ellos –para el
funcionamiento del sistema (R.E.
Shannon)
Simulación en Ciencias Geológicas:
Para que?
Existen múltiples aéreas de la investigación en Ciencias de la Tierra, así
como de su aplicación, que recurren intensivamente a la simulación por
computadora y, en algunos casos, no serían posibles en absoluto sin ella.
 Tectónica de placas.
 Movimiento de fallas
 Predicciones:
 Erupción de un volcán
 Desbordamiento de ríos
 Terremotos
 Geología de combustibles fósiles:
 Petróleo: prospección y perforación
 Gas
 Exploración minera
 Recursos naturales (preservación)
Experimentación Real y Simulación
La experimentación directa sobre la realidad puede traer algunos
problemas:
 Costo muy alto
 Gran lentitud
 En ocasiones las pruebas son destructivas
 A veces puede no ser ética (sobre todo si están involucrados
seres humanos)
 Puede resultar imposible (por ejemplo para predecir sucesos
futuros)
Origen
 La construcción de modelos tiene su origen en la época del
renacimiento, aunque la palabra „Simulación‟ aparece cerca de 1940,
cuando los científicos Von Neumann y Stanislau Ulam, que trabajaban en
el proyecto Montecarlo, durante la Segunda Guerra Mundial. El
proyecto consistía en desarrollar la primera bomba atómica antes que
Alemania. Trabajaron en conjunto muchos científicos de EEUU, Reino
Unido y Canadá, entre ellos Albert Einstein.
 En el proyecto se hizo referencia a la simulación Montecarlo o método
Montecarlo. Se trata de un método usado para aproximar expresiones
matemáticas complejas y costosas de evaluar con exactitud. El método se
llamó así en referencia al Casino de Montecarlo, por ser la capital de los
juegos de azar y por ser la ruleta el método más simple de generación de
números aleatorios. Con la utilización de las computadoras en los
experimentos de simulación, surgieron muchísimas aplicaciones y
además una mayor cantidad de problemas teóricos y prácticos.
Introducción a la generación de
números pseudoaleatorios
 Se denomina variable aleatoria discreta aquella que sólo
puede tomar un número finito de valores dentro de un intervalo.
Ejemplo: el resultado de lanzar un dado o una moneda
 Una variable aleatoria discreta X se dice que es uniforme en
el intervalo [1,n] si la probabilidad de ocurrencia es la misma para
cada posible valor de la variable. Es decir,
P[X=j] = 1/n para j=1,…n.
 Casi todos los métodos de simulación se basan en la posibilidad de
crear números aleatorios con distribución U(0,1).
 Dato: Hasta la aparición de las computadoras, los números
aleatorios se obtenías de procedimientos experimentales como
lotería o ruleta y se almacenaban en tablas.
Método del Cuadrado Medio
 Es uno de los métodos más clásicos y más simples de generación
de números pseudoaleatorios.
 Fue creado por Von Neumann en 1940. La técnica parte desde una
semilla inicial, un número enteros de 2n cifras propuesto por el
usuario y los pasos a seguir son:
1.
2.
3.
4.
5.
Se toma un X0 como semilla, tal que X0 es entero y de 2n cifras
(por ej: x0=3234)
Se calcula y como x02 (Ejemplo y=10458756)
Se toman los 2n dígitos centrales. Este número servirá para formar
el primer número y para continuar calculando números
pseudoaleatorios (Ejemplo: 4587)
El primer número pseudoaleatorio es y/102n (Ejemplo: 0.4587)
La nueva semilla es 4587, o sea, se comienza el procedimiento
nuevamente para crear otro número pero con la nueva semilla
Métodos por computadora
 Los números generados por computadora se llaman números
pseudoaleatorios, dado que son predecibles a partir del primer número
denominado semilla
 Para poder utilizar un generador automático de números
pseudoaleatorios, éste debe cumplir con ciertas propiedades:
1.
2.
3.
4.
5.
6.

Producir muestras según la distribución U(0,1)
Pasar los contrastes de aleatoriedad e independencia más habituales
Que la sucesión generada sea reproducible a partir de la semilla
Tener una longitud de ciclo tan grande como se desee
Generar valores a alta velocidad
Ocupar poca memoria
Ejemplo propuesto por Anderson en 1990, basado en el reloj de la
computadora:
X0=iyear+100*(imonth-1 +12*(iday-11+31*(ihour+24*(imin-60*isec))))
Métodos por computadora
 En Matlab, se utiliza la función rand() para generar números
pseudoaleatorios. La sintaxis es la siguiente:
R=rand(n)
% Genera una matriz de NxN de numeros pseudoaleatorios
entre 0 y1
R=rand(m,n) % Genera una matriz de mxn
R=rand
% Genera un solo número pseudoaleatorio entre 0 y 1
 Así, podemos entonces simular la tirada de una moneda utilizando esta
función de Matlab. Sabemos que hay dos posibles resultados para esto:
Cara ó Ceca, o sea, la probabilidad para cada valor es 0.5
 Dividimos el intervalo [0,1] del siguiente modo: Si al utilizar rand nos da
una valor menor a 0.5, decimos que salió Cara. De otro modo, salió Ceca
Ejemplo
function s = moneda(n)
% función que simula n tiradas de una moneda.
% el vector s retorna 1 si sale Cara y 0 si sale Ceca
s = [];
for i =1:n
m =rand(); %también puede ser m=rand;
if m < 0.5
s = [s 1];
else
s = [s 0];
end
end
end
Histograma
 Un histograma es una representación gráfica de una variable
en forma de barras.
 La superficie de cada barra es proporcional a la frecuencia de
los valores representados.
 En el eje vertical se representan las frecuencias, y en el eje
horizontal los valores de las variables, normalmente
señalando las marcas de clase, es decir, la mitad del intervalo
en el que están agrupados los datos.
 Se utiliza cuando se estudia una variable continua, como
franjas de edades o altura de la muestra, y, por comodidad,
sus valores se agrupan en clases, es decir, valores continuos.
Histograma
 Ejemplo:
Histograma
 En Matlab:
>> v=[12 19 25 26 32 18 33 27 34 28 35 21 36 29 23 37 30
33 14 36 26 34 13 37 28 34 30 32 31 33 15 44 29 48 16
31 13 27 28 25 ]
9
8
>> hist(v)
7
6
5
4
3
2
1
0
10
15
20
25
30
35
40
45
50
Histograma
>>centers= [14 21 28 35 42 49]
>>hist(v,centers)
15
10
5
0
14
21
28
35
42
49
Ejemplo
 Desde la línea de comandos de Matlab:
>> m= moneda(100); % Simulamos la tirada de 100
monedas
 ¿Cómo podemos visualizar la distribución de las tiradas? Para
eso utilizamos la función hist(m)
>> N = hist(Y)
% Toma cada elemento del vector Y y lo
ubica en alguna columna de las 10
disponibles
En nuestro caso, realizamos 100 tiradas de una moneda y
tenemos dos resultados posibles dentro del vector: 0 ó 1
Ejemplo
Realizamos hist(m) y vemos lo que sucede…
60
El intervalo [0,1] se divide en 10
partes y se ubican en la primera
barra todos los valores 0 del vector y
en la segunda todos los valores 1
50
40
30
20
10
0
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
Como vemos, no parece tan
uniforme….
Esto es porque sólo realizamos 100
tiradas. A medida que el número de
muestras sea mayor, el resultado se
acerca más a la uniformidad.
En este caso, marcamos los centros:
>>centers= [0.25 0.75]
>>hist(m,centers)
Ejemplo
Repetimos el experimento para 10000 tiradas
6000
En este ejemplo se puede
observar como de 10000 tiradas
tenemos un resultado uniforme
donde la mitad de las tiradas
corresponde al valor 1 (cara) y
la otra mitad al valor 0 (ceca)
En
este
caso
usamos
directamente el comando
hist(m) sin determinar los
centros, con lo cual el intervalo
(0,1) se dividió en 10 partes
iguales
5000
4000
3000
2000
1000
0
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
Ejemplo 2
 ¿Como hacemos para simular la tirada de un dado? En este
caso, tenemos 6 valores posibles que deben tener la misma
probabilidad de salir, o sea 1/6
 En general, para generar un número pseudoaleatorio entre a
y b la fórmula a aplicar es la siguiente:
r = a + (b-a).*rand;
En el caso del dado, el intervalo es [1,6], por lo que la
fórmula se reduce a:
r = 1 + 5.*rand
Con esta fórmula, generamos un numero al azar entre 1 y
6 pero en decimal!!!
Ejemplo 2
 En Matlab hay varias funciones para transformar un decimal a
un entero:
 ceil(x)
% Redondea x al entero siguiente. Ejemplo:
>>ceil(3.5)
ans
4
>>ceil(3.1)
ans
4
Ejemplo 2
 floor(x)
>>floor(3.5)
ans
3
>>floor(3.1)
ans
3
>>floor(3.75)
ans
3
% Redondea x al entero anterior. Ejemplo:
Ejemplo 2
 Con estos datos, volvamos a la función para generar un número
entero entre 1 y 6 para lanzar dados. Teníamos:
r = 1 + 5.*rand
Que generaba un número entre 1 y 6 pero decimal. Ahora, vamos
a usar funciones de truncamiento. Usamos floor:
r = 1 + floor(5.*rand)
Si lo dejamos así, dado que floor redondea siempre “hacia abajo”,
se generarán números entre 1 y 5. Hacemos un pequeño cambio:
r = 1 + floor(6.*rand)
Si queremos generar 100 números al azar entre 1 y 6 hacemos:
r = 1 + floor(6.*rand(100,1))
Ejemplo 2
 Vamos a probar que con esa función obtenemos una lista de
números con distribución uniforme. Para eso “tiramos” el
dado unas 100000 veces:
>> r = 1 + floor(6.*rand(100000,1));
>>hist(r)
18000
16000
14000
12000
10000
8000
6000
4000
2000
0
1
1.5
2
2.5
3
3.5
4
4.5
5
5.5
6
Nota al margen
 El comando rand genera números entre 0 y 1. Pero no están incluidos ni el cero ni el 1.
 Por lo tanto con la fórmula anterior para simular la tirada de un dado, no corremos
riesgo de generar el valor 7.
 De todas maneras, si tuviéramos que excluir valores de un vector, podemos usar el
comando find:
>>I = find(X)
%Si X es un vector, retorna en I un vector con los índices de los elementos
de X. Ejemplo:
>>X=[13 5 7 89 1]
>>I=find(X)
I=[1 2 3 4 5]

El parámetro de find puede ser una expresion lógica:
>>I=find(X<10)
%Queremos buscar todos los índices del vector X que correspondan a los
elementos menores a 10
I=[2 3 5]
>> Y=X(I)
Y=[5 7 1]
%Estamos construyendo un nuevo vector Y con los valores de menores a 10
Números pseudoaleatorios
 Ventajas:
 La sucesión generada es reproducible, usando la misma semilla.
 Longitud arbitraria de la secuencia.
 En general son algoritmos rápidos.
 Consume poca memoria.
 Desventajas:
 La distribución no es estrictamente uniforme.
 Claramente, a partir de cierto valor, los números comenzaran a
repetirse.
Números pseudoaleatorios
Probamos diferentes distribuciones:
function [u1, u2] = pruebaDist(N)
u1 = rand(N,1); %probar con N=10000 o más
subplot(1,2,1)
hist(u1)
title('Histograma Distribución Uniforme');
end
u2=randn(N,1);
subplot(1,2,2)
hist(u2)
title('Histograma Distribución Normal');
Números pseudoaleatorios
12
4
x 10Histograma Distribución Uniforme
4
5
x 10 Histograma Distribución Normal
3.5
10
3
8
2.5
6
2
1.5
4
1
2
0.5
0
0
0.2
0.4
0.6
0.8
1
0
-10
-5
0
5
Movimiento Browniano
 El movimiento browniano es el movimiento aleatorio que se
observa en algunas partículas microscópicas que se hallan en un medio
fluido (por ejemplo, polen en una gota de agua).
 Recibe su nombre en honor al escocés Robert Brown, biólogo y
botánico que descubrió este fenómeno en 1827 y observó que
pequeñas partículas de polen se desplazaban en movimientos aleatorios
sin razón aparente.
 Por un momento pensó que se trataba de la “vida” que existía dentro del
polen, sin embargo, repitió el experimento con diferentes partículas de
polvo obteniendo resultados similares
 De sus observaciones y las de otros científicos se pudieron obtener un
par de conclusiones: que las partículas presentaban mayor movimiento
entre más pequeñas fueran y que éste aumentaba también al
incrementar la temperatura del líquido. A este tipo de movimiento
azaroso se le dio el nombre de movimiento browniano en su honor.
Ejercicio: Trayectoria de una partícula
con movimiento browniano
 Definicion: El movimiento browniano se define como un
proceso estocástico B(t): t>0 satisfaciendo:
 B(0)=0;
 B(t) es continuo o sea, las trayectorias del proceso son funciones
continuas
 Fijados n instantes 0<= t1<=….<=tn los incrementos Btn – Btn-1
- ….- Bt1 son variables aleatorias independientes y cada
incremento Bt – Bs posee una distribución ~N(0, 1), para
cualquier s>=0 y t>0.
Ejercicio: Trayectoria de una partícula
con movimiento browniano
 Realizar una función que simule el movimiento browniano
de una partícula en un cierto fluido y graficar el resultado
1
0.5
xn
0
-0.5
-1
-1.5
-2
-0.2
0
0.2
0.4
0.6
0.8
n
1
1.2
1.4
1.6
1.8
Codigo Matlab
function [X,Y]= browniano(N)
deltax=1/N;
X(1)=0;
Y(1)=0;
for i=2:N
dW =sqrt(deltax)*randn;
X(i)=X(i-1)+ dW;
dW =sqrt(deltax)*randn;
Y(i)=Y(i-1)+ dW;
end
Desde línea de comandos:
>> [x, y]=browniano(1000);
>>plot(x,y)
>>grid on
>>xlabel(‘x');
>>ylabel(‘y');
Ejercicio: Trayectoria de una partícula
con movimiento browniano
 Modificar la función anterior para que simule el
movimiento browniano de n partículas en un cierto
fluido y graficar el resultado
1.5
4
3
1
2
0.5
xn
xn
1
0
0
-0.5
-1
-1
-1.5
-0.5
-2
0
0.5
1
n
3 partículas
1.5
2
2.5
-3
-2.5
-2
-1.5
-1
-0.5
0
n
0.5
100 partículas
1
1.5
2
2.5
Codigo Matlab
function [X,Y]= Pbrowniano(N,P)
deltax=1/N;
for part=1:P
X(1,part)=0;
Y(1,part)=0;
for i=2:N
dW =sqrt(deltax)*randn;
X(i,part)=X(i-1,part)+ dW;
dW =sqrt(deltax)*randn;
Y(i,part)=Y(i-1,part)+ dW;
end
end
Desde línea de comandos:
>> [x, y]= Pbrowniano(1000,3);
>>plot(x,y)
>>grid on
>>xlabel(‘x');
>>ylabel(‘y');
Ejercicios de repaso
 Generar N números al azar entre 0 y 1 y hacer un histograma con los resultados






(rand)
¿Que representa la altura del histograma?
En el primer ejercicio, ¿Por que si la distribución es uniforme el diagrama se ve
irregular? ¿Qué pasa si aumentan el N?
¿Qué intervalo representa el ancho de cada barra? Si tiene dudas, busque en la ayuda
de Matlab la función hist.
Realice el mismo ejercicio con randn¿ Qué diferencias observa entre ambos tipos de
distribución?
¿Cómo haría para generar números al azar entre 0 y 15? ¿y entre -3 y 2? Hacerlo y
construir nuevos histogramas. Ayuda: Utilice la misma función rand y expanda el
resultado apropiadamente. Si es necesario, vea uno de los ejemplos de la Ayuda de
Matlab para la función rand
Modificar la función browniano(n) para que dibuje la simulacion en 3D