PDF de programación - Cálculo de probabilidades con Scilab

Cálculo de probabilidades con Scilab
En este tutorial se va a asumir que ya sabemos crear funciones y trabajar con bucles en scilab.
Sin embargo, también se puede usar como práctica para aprender estos conceptos mientras se
hace algunas prácticas de cálculo de probabilidad.


Preparación

Fabricamos un dado
Una de las utilidades prácticas de la programación es simular situaciones de la vida
real para obtener datos que nos llevarían mucho tiempo de otra forma. Por ejemplo, si
queremos “comprobar” las probabilidades de jugar con dados diferentes “juegos” deberíamos
realizar esos juegos miles de veces y contar los resultados obtenidos. Mucho más rápido es
que la computadora tire los dados por nosotros y nos cuente lo que da.

Para empezar, entonces, vamos fabricar un dado digital. Osea, una función que al ejecutarla
nos devuelva un número aleatorio entre 1 y 6 con igual probabilidad. Para ello vamos a utilizar
la función del sistema rand()

rand() devuelve en cada ejecución un número aleatorio mayor e igual a cero y menor que
uno.

-->rand()
ans =

0.2937103

Si le introducimos como argumento números de filas y columnas, nos da una matriz aleatoria:

-->rand(2,3)
ans =

0.0442440 0.3614060 0.0986330
0.8900476 0.8806820 0.2604110

La computadora no puede “inventar” números por sí misma. Lo que realiza cada vez que
ejecutamos esa función es un algoritmo de generación de números pseudoaleatorios. El
algoritmo parte de un número y en cada paso hace un conjunto de cálculos que el entregan
otro número “bastante aleatorio”. Si queremos que no salgan los mismos número aleatorios
cada vez que prendemos la compu de cero y ejecutamos el programa, debemos cargarle una
semilla nueva. Eso se hace con el comando:
rand('seed',123);

para cargar la semilla 123 . En la ayuda del comando (help rand) recomienda utilizar el valor
del segundero de la computadora como semilla, así en cada ejecución de la rutina la semilla
será diferente:

rand('seed',getdate('s'));

Con esto tenemos una herramienta para fabricar una función dado. Por ejemplo, si tomamos un
número en 0 y 1 , lo multiplicamos por 6:
-->rand(1,5)*6
ans =

1.283994 3.4722662 0.1433145 4.7805427 2.6354033

y luego redondeamos para abajo con la función floor
-->floor(rand(1,5)*6)
ans =

3. 3. 0. 3. 2.

Tenemos algo muy cercano a lo deseado. Como rand() nunca nos devuelve 1, no aparecerá
nunca el 6. Nos quedan números entre 0 y 5 . Con sumar uno ya tenemos nuestro dado.
Armando la función en el editor, quedaría:

function rta=dado()
rta=floor(rand()*6)+1;
endfunction


Probamos la función en la consola ejecutando dado() muchas veces:

->dado()
ans =

5.

-->dado()
ans =

4.


-->dado()
ans =

2.

-->dado()
ans =

6.


Probamos el dado
Si nuestro dado es de buena calidad, la probabilidad de que cada número salga debería
ser parecida. Una forma de analizar esto es tirando el dado muchas veces y graficando un

histograma.

Para tirar el dado muchas veces nos valemos del bucle for y guardamos los datos en alguna
variable (por ejemplo, datos):

for i=1:10000 //tiramos el dado 10mil veces
datos(i)=dado();
end

Para graficar el histograma existe la función histplot() . La forma rápida de usarla es:
histplot(n,vector)
donde n es el número de bines del histograma y vector es el vector de datos.

Ejemplo:
histplot(6,datos)

También podemos pasar un vector con los límites donde comienza y termina cada bin. Esto
permite centrarlos mejor:

histplot(0.5:6.5,datos)

Se le está pasando como límites de los bines : [0.5 1.5 2.5 3.5 4.5 5.5 6.5]

Finalmente, ponerle nombre a los ejes para cumplir con las formalidades:

histplot(0.5:6.5,datos)
xlabel("Cara del dado")
ylabel("Probabilidad")
Se puede apreciar que la frecuencia de cada valor de dado es muy similar. Tenemos un dado
aceptable:





Algunos juegos con dados

Tirar dos dados y sumarlos
De forma similar a cómo se probó el funcionamiento del dado, podemos simular que tiramos
dos veces el dado y lo sumamos. Los números que pueden resultar irán entre 2 (1+1) y 12
(6+6), por lo que hay que acomodar los límites de nuestro histograma.



Tirar dos dados y multiplicarlos
Siguiendo el mismo procedimiento:

for i=1:10000 //tiramos el dado 10mil veces
datos(i)=dado()*dado();
end
histplot(0.5:36.5,datos)
xlabel("Producto de dados")
ylabel("Probabilidad")


Se aprecia de inmediato que hay valores que no aparecen. Esto puede ser porque son número
primos (ej: 11, 13, 17) o porque son múltiplos de números que no están en los dados (ej: 14,
21). Para tener una mejor noción de la distribución de probabilidad se puede agrandar el ancho
de los bines haciendo que los límites especificados tengan 6 unidades de separación:

histplot(0:6:36,datos)
xlabel("Producto de dados")
ylabel("Probabilidad")



Generala

Vamos a simular un juego de generala. Se trata de tirar 5 dados por turno y ver las
probabilidades de tener diferentes coincidencias. Como vamos a trabajar con muchos números,
es conveniente saber un poco de eficiencia de algoritmos. Osea, cómo escribir algo para que
funcione lo más rápido posible en la computadora.

Una forma de recopilar los datos de 10 mil tiradas de generala sería la siguiente:

for i=1:10000 //tiramos el dado 10mil veces
datos(i,:)=[dado(),dado(),dado(),dado(),dado()];
end

Existen dos comandos que permiten averiguar cuánto tarda el scilab en realizar dicho
proceso:tic y toc. Cuando ejecutamos tic, scilab empieza a contar los segundos que pasan.
Cuando ejecutamos toc, nos informa en pantalla cuanto tiempo pasó. Para medir este proceso
en particular, ejecutamos:

tic
for i=1:10000 //tiramos el dado 10mil veces
datos(i,:)=[dado(),dado(),dado(),dado(),dado()];
end
toc

ans =

1.768


En este caso se tardó 1.7 segundos

Otra forma de realizar el mismo proceso es la siguiente:

tic
for i=1:10000 //tiramos el dado 10mil veces
datos(i,:)=floor(rand(1,5)*6)+1;
end
toc

Acá, en lugar de ejecutar 5 veces la función dado(), escribimos lo que la función dado hace
pero para una matriz aleatoria generada con 5 elementos. El procedimiento matemático es el
mismo, pero el procesamiento se saltéa pasos internos que el scilab debe hacer al llamar una
función. En este caso el tiempo resultó:
ans =

0.101

Es casi 10 veces más rápido. Por último, una versión aún más rápida, que no requiere iterar, es
crear directamente una matriz aleatoria de 5 x 1000 :
tic
datos=floor(rand(10000,5)*6)+1;
toc


ans =

0.018

Se puede ver que con esta última versión el código es 100 veces más rápido. Ahora, si
queremos generar un conjunto de 100mil datos, el primer código tardará algunos minutos y el
último apenas unos segundos.

Hechas las consideraciones, podemos generar 100 mil tiradas de dados de generala:

-->datos=floor(rand(100000,5)*6)+1;

Como ejemplo, podemos ver las primeras 10 tiradas:

-->datos(1:10,:)
ans =

3. 5. 3. 3. 3.
1. 2. 5. 2. 1.
6. 4. 3. 6. 6.
3. 6. 5. 4. 3.
3. 2. 2. 5. 6.
2. 3. 1. 2. 4.
4. 2. 3. 3. 4.
6. 1. 5. 2. 2.
2. 1. 2. 3. 5.
6. 3. 5. 4. 5.


Análisis de los datos
Para saber si una tirada es generala, full o alguna otra figura debemos analizar cada fila de los
datos. Es conveniente pensar en una función para ello. Armaremos una función que nos cuente
cuantas veces aparece repetido cada número. Por ejemplo, en la tirada:
a=[ 1 3 5 3 6 ];
queremos que la función nos devuelva los números de cuántas veces se repiten cada uno.
Podemos empezar con la función find(), que nos encuentra las posiciones de un dado
elemento:
-->find(a==3)
ans =

2. 4.

Lo que nos interesa nosotros no es la posición, sino la cantidad, por lo que podemos contar
cuantas posiciones nos devuelve find:
-->length(find(a==3))
ans =

2.

Ya casi estamos. Solo debemos escribir una función que nos devuelva un vector con la
cantidad de veces que aparece cada cara del dado:




Para nuestro ejemplo a, nos dice que un número aparece 2 veces y el resto 1 o 0. La función
gsort ordena de mayor a menor el vector de la respuesta. Así es más fácil comparar, ya que no
nos interesa “qué número” se repite, sino las combinaciones de repeticiones.

Por último, solo debemos comparar las respuestas de nuestra función “tipo_de_tirada” con lo
valores típicos de las figuras de generala:

//Tipos de resultados que podemos tener:
nada= [1,1,1,1,1,0];
par= [2,1,1,1,0,0];
pardoble=[2,2,1,0,0,0];
triple= [3,1,1,0,0,0];
fulls= [3,2,0,0,0,0];
pocker= [4,1,0,0,0,0];
generala=[5,0,0,0,0,0];

Por último, armamos un bucle for que pase por todos los datos fabricados y nos guarde el tipo
de dato en un vector. Por ejemplo, podemos hacer que cuando encuentre un par nos guarde en
un vector el valor 2, y cuando encuentre un pardoble, nos guarde el valor 3, etc , siguiente una
convención que resulte conveniente:
//guardamos el análisis en el vector "tipo"
for i=1:100000
tmp=tipo_de_tirada(datos(i,:));
if isequal(tmp,nada) then
tipo(i)=1;
end
if isequal(tmp,par) then

tipo(i)=2;
end
if isequal(tmp,pardoble) then
tipo(i)=3;
end
if isequal(tmp,triple) then
tipo(i)=4;
end
if isequal(tmp,fulls) then
tipo(i)=5;
end
if isequal(tmp,pocker) then
tipo(i)=6;
end
if isequal(tmp,generala) then
tipo(i)=7;
end
end
Cabe señalar que para comparar cada fila con alguno de los tipos de tirada definidos arriba se
usa la función isequal(). isequal(a,b) nos devuelve “verdadero” si y sólo si todos los
elementos de a coinciden con los de b índice a índice.

Finalmente podemos graficar los resultados del