[MODELADO Y SIMULACIÓN] Pruebas estadísticas para números pseudoaleatorios

En esta tarea se seleccionara una librería de cualquier tipo que genere números pseudoaleatorios siguiendo una distribución, y comprobaremos que tan aleatorios son estos números generados.

Librería generadora de números pseudoaleatorios

Seleccioné la distribución logística generada por la librería numpy (numpy.random.logistic)

Prueba estadística

Cuando se analizan datos medidos por una variable cuantitativa continua, las pruebas estadísticas de estimación y contraste frecuentemente empleadas se basan en suponer que se ha obtenido una muestra aleatoria de una distribución de probabilidad de tipo normal o de Gauss.

Pero en muchas ocasiones esta suposición no resulta válida, y en otras la sospecha de que no sea adecuada no resulta fácil de comprobar, por tratarse de muestras pequeñas.

En estos casos disponemos de dos posibles mecanismos:

Los datos se pueden transformar de tal manera que sigan una distribución normal.
O bien se puede acudir a pruebas estadísticas que no se basan en ninguna suposición en cuanto a la distribución de probabilidad a partir de la que fueron obtenidos los datos, y por ello se denominan pruebas no paramétricas (distribución free), mientras que las pruebas que suponen una distribución de probabilidad determinada para los datos se denominan pruebas paramétricas.

Seleccioné la prueba de Anderson-Darling de la librería scipy para comprobar que la distribución en los números generados es realmente la seleccionada. También se realizó un test de frecuencia para saber que tan aleatorios son los números generados por la librería, en este caso monobit.

Programa

En la práctica me basé en el código de Cecilia y de Juan Carlos para generar los números pseudoaleatorios, utilizando las herramientas anteriormente descritas.

import numpy
from scipy import stats

ls = []

f = open("datos.txt","w")

for i in range(1000):
    dist = numpy.random.logistic(1, 1, None)
    ls.append(dist)
    f.write("%0.3f\n" %dist)

anderson, critical, sig = stats.anderson(ls, dist="logistic")

for i in range(len(critical)):
    if anderson > critical[i]:
        print "Hipotesis rechazada con valor de significancia "+str(sig[i])
    else:
        print "Hipotesis no rechazada con valor de significancia "+str(sig[i])

El código anterior genera 1000 números pseudoaleatorios siguendo la distribución logística. y dando como resultadoestos números graficados, que por lo que se ve, parese que son muy aleatorios:

El programa hace un test de anderson dando como resultado el siguiente:

Parece ser que la distribución es la correcta, puesto que las hipotesis o se rechazan.

Ahora veremos que tan aleatorios son los números mediante monobit.

import math, numpy
from scipy.special import erf
                                                                                                                                                                                                                 
def mono(lista):
    lista2 = [-1 for i in lista if i==0]
    suma = sum(lista) + sum(lista2)
    sumaAbs = abs(suma)/math.sqrt(len(lista)) 
    p_value = erf(sumaAbs/math.sqrt(2))
    
    if(p_value > 0.01):
        print "%f: P" % p_value
    else:
        print "%f: NP" % p_value

def main():
    cantidad = 1000
    L1 = []
    
    for i in range(cantidad):
        x = numpy.random.logistic(1,1, None)
        x = abs(int(round(x)))
        if(x > 1.0):
                x = 1.0
        elif(x < 0.0): 
                x = 0.0
        L1.append(x)
    
    mono(L1)

main()

El código es basado en el del Compañero Juancarlos para realizar la prueba de frecuencia monobit.

El resultado es el siguiente:

Noté que al correr la prueba varias veces, de vez en cuando fallaba. Sin embargo la mayoría de las veces esta pasa. Se puede concluir que genera números aleatorios de forma muy eficiente.

[MODELADO Y SIMULACIÓN] Sistemas caóticos

Diagrama de Bifurcación:

Un diagrama de bifurcación de un sistema dinámico es una estratificación de su espacio de parámetros inducida por la equivalencia topológica, junto con los retratos de fase representativos de cada estrato. Las soluciones estables suelen representarse mediante líneas continuas, mientras que las soluciones inestables se representan con líneas punteadas.

Código en python

import math
max = 480
a = open("datos.txt", "w")

xa = 2.9
xb = 4.0
maxit = 250

for i in range(max):
    r = xa + (xb - xa) * float(i) / (max - 1)
    x = 0.5
    for j in range(maxit):
        x = r * x * (1 - x)
     a.write(str(x)+"\n")
a.close

Código gnuplot

set term postscript eps color 30
set key outside Right
set size 2, 1.5
set output "simugeom.eps"
set pointsize 1
plot "datos.txt" with points lt 2

Billar dinámico:

(Este ejemplo lo pongo por que me pareció divertido)

Un billar dinámico es un sistema dinámico en el cual una partícula alterna entre movimiento rectilíneo y reflexiones especulares en un contorno o frontera.

Los billares capturan toda la complejidad de los sistemas hamiltonianos, desde integrabilidad a movimiento caótico, sin las dificultades de integrar las ecuaciones de movimiento para determinar su mapa de Poincaré.

Código en python

import math
import random
from PIL import Image, ImageDraw
imgx = 480
imgy = 320
image = Image.new("RGB", (imgx, imgy))
draw = ImageDraw.Draw(image)

# Only 1 ball is used!
maxSteps = 20000 # of steps of ball motion (in constant speed)

n = random.randint(1, 7) # of circular obstacles
crMax = int(min(imgx - 1, imgy - 1) / 4) # max circle radius
crMin = 10 # min circle radius

# create circular obstacle(s)
cxList = []
cyList = []
crList = []
for i in range(n):
    while(True): # circle(s) must not overlap
        cr = random.randint(crMin, crMax) # circle radius
        cx = random.randint(cr, imgx - 1 - cr) # circle center x
        cy = random.randint(cr, imgy - 1 - cr) # circle center y
        flag = True
        if i > 0:
            for j in range(i):
                if math.hypot(cx - cxList[j], cy - cyList[j]) < cr + crList[j]:
                    flag = False
                    break
        if flag == True:
            break
    draw.ellipse((cx - cr, cy - cr, cx + cr, cy + cr))
    cxList.append(cx)
    cyList.append(cy)
    crList.append(cr)

# initial location of the ball must be outside of the circle(s)
while(True):
    x = float(random.randint(0, imgx - 1))
    y = float(random.randint(0, imgy - 1))
    flag = False
    for i in range(n):
        if math.hypot(x - cxList[i], y - cyList[i]) <= crList[i]:
            flag = True
            break
    if flag == False:
        break

# initial direction of the ball
a = 2.0 * math.pi * random.random()
s = math.sin(a)
c = math.cos(a)

counter = 0
j=0

for i in range(maxSteps):
    counter = counter +1
    if counter == 15:
        j=j+1
        image.save(str(j)+".png", "PNG")
        counter = 0
    image.putpixel((int(x), int(y)), (255, 255, 255))
    xnew = x + c
    ynew = y + s
    
    # reflection from the walls
    if xnew < 0 or xnew > imgx - 1:
        c = -c
        xnew = x
    if ynew < 0 or ynew > imgy - 1:
        s = -s
        ynew = y
    
    # reflection from the circle(s)
    for i in range(n):
        if math.hypot(xnew - cxList[i], ynew - cyList[i]) <= crList[i]:
            # angle of the circle point
            ca = math.atan2(ynew - cyList[i], xnew - cxList[i])
            # reversed collision angle of the ball
            rca = math.atan2(-s, -c)
            # reflection angle of the ball
            rab = rca + (ca - rca) * 2
            s = math.sin(rab)
            c = math.cos(rab)
            xnew = x
            ynew = y
    
    x = xnew
    y = ynew

Cuando encontre el código originalmente generaba una imagen al finalizar el programa, pero lo modifique para que generara una imagen cada 15 iteraciones y de esa manera observar como se desarrolla.

Para generar imagenes se utiliza la libreria PIL y para generar un video de una secuencia de imagenes se utiliza el comando "ffmpeg -f image2 -i imagen%d.png video.mpeg" donde %d es la numeración de imagenes, ejemplo: imagen1.png, imagen2.png, etc.

Este es el resultado de 1333 png's generados:

[MODELADO Y SIMULACIÓN] Procesos estocásticos discretos y contínuos

Para mi ejemplo decidí hablar sobre la probabilidad de pase un camión de cierta ruta entre las 6:00 y 6:30.

Estableciendo el experimento, hay una persona esperando el camion de la ruta 213 cuyo porcentaje de probabilidad de que pase es de 30%. ¿Cuál es la probabilidad de que el primer camión que pase sea de la ruta 213?

Para saber la distribución a la que pertenece este problema, analizaré las distribuciones disponibles.

Distribuciones Discretas: 


Son aquellas cuya función de probabilidad toma valores positivos en un conjunto finito o infinito para x:

Lo anterior representa la suma de todas las probabilidades desde -∞ hasta el valor de x.

Las distintas distribuciones mas importantes de la anterior son:

Distribución binomial

Distribución binomial negativa

Distribución Poisson

Distribución geométrica

Distribución hipergeométrica

Distribución de Bernoulli

entre otras derivadas de las anteriores.

Distribuciones Continuas:

Son las que su variable toma cualquier valor dentro de un intervalo existente. Está dada por la integral de la función densidad:

Las más conocidas son:

Distribución ji cuadrado

Distribución exponencial

Distribución t de Student

Distribución normal

Experimento:

Se puede decir del planteamiento que la distribución que sigue es discreta.

La probabilidad de éxito es de 0.3, y la de fracaso de 0.7 y puesto que solo tenemos de 6:00 a 6:30 para  tomar el camión, necesitamos que el primero sea el 213.

La formula para esto es:

donde:


P (X=x) = función de densidad, de la variable aleatoria con distribución geométrica.

X = Numero de experimentos hasta que aparece el 1er éxito.

p = probabilidad de éxito

q = probabilidad de fracaso (1 - p)

Para el porblema tenemos que;

X = 1

p = 0.3

q = 0.7

P(X=1) = [(0.3^(1-1)] (0.7) = 1(.0.7) = 0.7

Utilizando el código de clase surge la siguiente gráfica:

El problema puede seguir una distribución continua poniendo pequeños intervalos, en este caso utilizando la exponencial.

Modificando el código de clase a este:

function exp(q, filename)
  p = 1 - q;
  output = fopen(filename, "w");
  prob = 1;
  k = 1;
  lam = 1;
  while (prob > 0.0)
    prob = lam * e ** (-lam * k);
    #prob = q**(k-1) * p;
    k++;
    fprintf(output, "%d %f\n", k, prob);
    lam = lam ;
  endwhile
  fclose(output);
endfunction

[MODELADO Y SIMULACIÓN] Sistemas deterministas

Un sistema determinista es aquel que si se conoce su comportamiento y que cambios o acciones conllevan a un estado, se sabrá cuales serán sus futuros estados.

Ésto es que si despues de un determinado tiempo o mediante las acciones del entorno un objeto cambia de estado, este es determinista, siempre y cuando sus tareas se hayan cumplido.

Desde mi punto de vista cualquier cosa que esté hecha para seguir una serie de reglas es un sistema determinista, así sea un foco o bien un sistema más complejo como un circuito que haga determinada tarea.

Un ejemplo sería un cruce peatonal con semáforos para peatones.

Algunos cruces peatonales tienen estos semáforos con un botón para pasar. Éstos están programados para cambiar cada cierto tiempo, o bien, cambian dependiendo de cuantas veces se pulsó el botón de cruce.

Es determinista porque si dado el tiempo determinado nadie oprime el botón, este cambiara a verde para los peatones y rojo para el tráfico. Además e que si el botón es oprimido un número dado de veces, este cambiara a verde en un tiempo determinado, o bien si el tiempo dado se termina.


Si bien, el número de personas es una variable aleatoria, el semaforo está listo para cambiar ante un número definido de pulsaciones del botón. Quiere decir que su tarea no se ve afectada por esta variable, su estado cambiara ante el tiempo determinado sin fijarse en que momento llegan las personas ni el orden en que llegaron. Es solo una condicion que le da preferencia a un numero dado de personas o mayor a éste.

Como se ve en las imágenes anteriores, iniciando con el rojo, se dará un tiempo determinado dependiendo de las veces que se pulsó el botón, o bien este sera el tiempo total por defecto, y cambiara su estado a verde.

El siguiente pseudocódigo dará una idea de la forma que funcionan:

estado = verde
tiempo = 1 minuto

def cronometro():
    while x = true:
        tiempo -= 1 segundo
        if tiempo == 1 minuto:
            x = false
            return false

def semaforo():
    if cronometro == false:
        estado = rojo
    elif boton() == false:
        if tiempo <= 45 segudos:
            estado = rojo
        else:
            tiempo -= 45 segundos

def boton():
    while x = true:
        if se presiona boton:
            contador += 1
        if contador == 10:
            x = false
            return false

while true:
    semaforo()

Éste es el diagrama de estados ligeramente distinto al pseudocódigo para verse más simple: