Modelado y Simulación

Proyecto Final: Texturas de Madera

2012-05-29T13:32:00.000-07:00

Presentación Proyecto Final

View more presentations from Saul Gausin.

Código en Python:

Imágen generadas:

Ecuación de Lotka-Volterra

2012-05-28T08:30:00.001-07:00

La ecuación de Lotka-Volterra, también conocida como las ecuaciones de depredador-presa, que en sí son ecuaciones de primer orden, no líneales, que son usados para describir sistemas dinamicos biologicos, en donde dos especies interactuan, donde uno es depredador y otro es presa. En donde se representa la ecuación de la siguiente manera :

Donde,

x es el número de presas.
y es el número de depredadores.
"dy/dt" y "dx/dt" representan el crecimiento de las dos poblaciones contra tiempo.
t representa el tiempo.
alfa, δ, ß y γ son parámetros que representan la interacción entre dos especies.

Ahora lo haremos en python, para posterior mente gráficarlo con gnuplot.

#sigma, alfa, beta y r floats                                                                                                                                                        
# x  y ints                                                                                                                                                                          
f =open("algo.dat","a")
sigma = float(raw_input("Sigma"))
alfa = float(raw_input("Alfa"))
beta = float(raw_input("Beta"))
r = float(raw_input("R"))
x = float(raw_input("X"))
y = float(raw_input("Y"))
iter = 100
dt=0.01
t=0.01
for i in range(iter):
    dx = r*x-alfa*x*y
    dx = dx*dt
    dy = -sigma*y-beta*x*y
    dy = dy*dt
    (x,y) =(x+dx, y+dy)
    if x < 0:
        x = 0
    if y < 0:
        y = 0
    print >> f, '%g %g %g \n'%(t,x,y)
    t+=dt
f.close()

Ahora lo que despliega en gnuplot es lo siguiente :

Fuentes de información:

http://en.wikipedia.org/wiki/Lotka%E2%80%93Volterra_equation

Sistemas caóticos: Envío de mensajes secretos por medio de caos

2012-05-11T13:43:00.003-07:00

Una de las grandes aplicaciones de caos ha sido utilizado en las comunicaciones privadas, ya sea para ocultar información enviada a un particular hasta ser una completamente clasificada.Comúmente la personas usan código para poder enviar mensajes secretos, pero al paso del tiempo estos logran ser decifrados.

Kevin Cuomo y Alan Oppenheim en el año de 1992 se enfocaron en resolver este problema basandose en el descubrimiento de Pecora y Carroll sobre caos sincronizado.La estrategia era de que cuando se transmita dicho mensaje hacia el receptor, este iba a ir enmascarado con un sonido más fuerte usando caos.Por consecuencia alguien que recepte ese mensaje solamente oirá el caos, en otras palabras escuchará solamente puro sonido sin sentido como cuando enciendes la televisión y no tienes señal.

El receptor lo que hará será solamente usar un traductor que reprodusca el caos para así substraer la máscara de caos y así sustraer el mensaje oculto.

Para hacer esta demostración Kevin Cuomo desarrolló un proyecto en el cual por medio de un circuito analógico implemento las ecuaciones de lorenz para crear una máscara caótica en mensajes, para hacer la demostración el reprodució una canción, para después mostrar la canción enmascarada la cual no tenía sentido alguno y luego de ser pasada hacia el receptor que este desenmascarado la máscara de caos y así se logró escuchar la canción que se escuchó con algo de ruido sin perder el mensaje.

En la que si lo veríamos gráficamente se podría acercar a algo como esto.

Esto puede ser demostrada por medio de algún programa por medio de un mapa lógistico :

En donde Xn es un número entre 0 y 1, dicho número lo tomaremos de manera pseudoaleatoria.r es un número positivo y este representa una tasa combinada para una reproducción y inanición.

Tomaremos un valor entre 3.57 y 4 para r ya que no todos los valores producen caos y estos valores son los más sensibles a producirlo.

Ahora tomamos un mensaje, el cual queremos evniar el cual llamaremos m.Luego sumamos el valor anterior generado.En esta ecuación obtendremos un valor nuevo llamado señal s.

Esta es la forma en la que nosotros enmascaramos el mensaje por medio de caos, si nosotros escuchamos dicho mensaje nosotros solamente oiremos ruido.

Ahora lo que necesitamos es reproducir el mensaje para ello supongamos que hacemos otro mapa lógistico.En esta ocación en lugar de iterar el valor actual de nuetro sistema para obtener el siguiente estado, promediaremos el estado actual por medio de la señal que ha sido transmitido.Llamaremos señal transmitida s y nuestro nuevo mapa será así.

Tomamos el mismo valor que tomamos en r, y lo que haremos será subtraer x' de la señal s para obtener un señal reconstruida .

Si recordamos s era la suma de m y x, así que sustituimos.

Para después multiplicarlo por 100 y así obtener el mensaje aproximado.

Ahora demostrare esto en un código de Python el cual encontre en una de las páginas que tengo como referencia en la parte de abajo, al código le modifique para que este nos gráficara que es lo que esta sucediendo.

El codificador.

El decodificador.

Video del mensaje, codificación y decodificación

Gráfica:

Fuentes:

http://jellymatter.com/2012/01/04/a-secret-message-from-another-dimension/

http://en.wikipedia.org/wiki/Logistic_map

Nonlinear Dynamics and Chaos por Steven H. Strogatz

Pruebas estadísticas para los números pseudoaleatorios

2012-04-20T15:42:00.000-07:00

En esta entrada hablaremos sobre los números pseudoaleatorios, generaremos un número pseudo-aleatorio por medio de un programa hecho en Ruby y de ahí nos pasaremos a hacer unas pruebas estadísticas sobre dichos datos.

Así es como se generan números aleatorios en Ruby:

archivo = File.new("rand.dat","w+")
for s in 1...100
      print s,"\n";
      archivo.print s," ",rand(),"\n";
end

Este genera un archivo llamado rand.dat en el cual se van guardando los datos, para que posteriormente sean gráficados.

Donde genera 100 números de manera pseudoaleatoria en un rango entre 0 y 1.

Para hacer dicha modificaré el código de Ruby para generar números que sean 0 ó 1:

archivo = File.new("rand.dat","w+")
for s in 1...100
      print s,"\n";
      archivo.print s," ",rand(2),"\n";
end

Ahora haremos una prueba llamada "Random Excursions Test" para comprobar que realmente son pseudoaleatorios, donde el propósito de está prueba es concentrarse en el número de ciclos que tienen exactamente K visitas en un paseo al azar la suma acumulada, para más información entra a está liga:

http://es.scribd.com/doc/11305936/133/Random-Excursions-Test

Donde se programó esta prueba en Python:

import math
import random
def tabla(x, k):
    val = 1.0 / (2.0 * math.fabs(x))
    return val * (1.0 - val)**(k - 1)

def lookup(lista):
    ciclo = list()
    c = None
    for e in lista:
        if e == 0:
            if c is not None:
                ciclo.append(c)
            c = list()
        else:
            c.append(e)
    return ciclo

file = open("rand.dat")
lista = list()
line = file.readline()
line = line.strip()
for elem in line.split():
    lista.append(int(elem)*2 -1)
n=len(lista)
s=list()

print lista
for e in range(0,n+1):
    s.append(sum(lista[:e]))

s.append(0)
print s
ciclo=lookup(s)
J = len(ciclo)

v=dict()
for x in range(-4, 5):
    if x != 0:
        for c in ciclo:
            oc = c.count(x)
            if oc > 5:
                oc = 5 # los mayores a cinco van con cinco
            if (x, oc) in v:
                v[(x, oc)] += 1
            else:
                v[(x, oc)] = 1
for x in range(-4, 5):
    if x != 0:
        ver = 0
        stat = 0.0
        for k in range(0, 6):
            if (x, k) in v:
                val = v[(x, k)] # frecuencia observada
            else:
                val = 0.0
            teor = J * tabla(x, k) # frecuencia esperada
            stat += (val - teor)**2 / teor
            ver += val
        print x, stat

Donde se normalizan los números pseudoaleatorios generados entre 0 y 1 donde 0 se cambia a -1, después de ello se hacen sumatorias parciales comenzando de X1 así sucesivamente hasta Xn, para luego examinar la cantidad de repeticiónes de frecuencias y compararlos apartir de la estadística Chi-square, lo que se esperaría una probabilidad de una distribución aleatoria.

Como resultado obtenemos lo siguiente:

Después comparamos con esta applet si realmente tienen una alta probabilidad de distribución aleatoria:

Con esto concluimos que a mayor cantidad sean los valores, hay mayor probabilidad de que se rechazen los valores y que a menor cantidad sean hay mayor probabilidad de que se acepten.

Se comprobó que los datos generados son uniformes y que puede mejorarse la aleatoridad. Se llegó a esta conclusión gracias a que dividimos los -1 y 1 entre los valores esperados para así mismo comprobarlo con el applet y esta misma applet nos determinó en todas las pruebas que se hicieron que nos daba una probabilidad de 0,9999.

Link de referencia:

Información de Random Excursions Test
http://es.scribd.com/doc/11305936/133/Random-Excursions-Test
Applet
http://www.danielsoper.com/statcalc3/calc.aspx?id=11

Procesos estocásticos discretos y contínuos

2012-03-13T06:32:00.000-07:00

En este tema tomaré el lanzamiento de un dado de seis lados como base de explicación de este tema, en donde se buscara la probabilidad de que obtengamos el número 6 .

Para tener un mejor entendimiento del problema se explicarán las distintas distribuciones discretas y continuas que existen, ya que con ello se entendería el por que el usó de la distribución que se eligió.

Distribuciones Discretas

En estas la variable puede tomar un número determinado de valores que pueden ser finitos o infinitos numerables; donde con infinito numerable se refiere a la cantidad que es muy grande y por serlo nadie se da lujo de enumerar, sin embargo tiene un límite y se puede cuantificar que es lo mismo a calcular.

Ejemplo:

La cantidad de galaxias que existen en el universo.
La cantidad de granos de arena en el mar.
Veces en las que cae un número par o impar en un dado.
El lanzamiento de una moneda en las que sale una cara o cruz.
Cantidad de hombres o mujeres en un salón.

En las distribuciones discretas se destacan los siguientes:

La distribución de Bernouilli es el modelo que sigue un experimento que se realiza una sola vez y que puede tener dos soluciones:

Éxito; donde toma el valor de 1

Fracaso; donde toma el valor de 0

Donde únicamente hay dos soluciones de sucesos complementarios:

A la probabilidad de éxito se le denomina "p"

A la probabilidad de fracaso se le denomina "q"

Por lo tanto:

p + q = 1

Ejemplo:

Probabilidad de acertar una quiniela:

Probabilidad de acertar p = 0.00001

Probabilidad de no acertar q = 0.99999

Por lo tanto:

p + q = 0.000001 + 0.99999

La distribución binomial es parte de la distribución de Bernouilli:

La distribución de Bernouilli se aplica cuando se realiza una sola vez un experimento que tiene

únicamente dos posibles resultados éxito o fracaso, por lo que la variable sólo puede tomar dos

valores: el 1 y el 0.

La distribución binomial se aplica cuando se realizan un número "n" de veces el experimento de

Bernouilli, siendo cada ensayo independiente del anterior, donde la variable puede tomar los

valores 0 (si todos los experimentos han sido fracaso) y n (si todos los experimentos han sido

éxitos).

Donde se sigue el siguiente modelo:

Para entenderlo veremos el siguiente ejemplo:

¿Cuál es la probabilidad de obtener cuatro veces el número 3 al lanzar un dado 8 veces?

Donde k es el número de aciertos en este caso se toma el valor 4, n es el número de veces lanzado por lo tanto el valor de n será 8 y p es la probabilidad de que salga un 3 al tirar el dado en este caso será 1/6 que es igual a 0.1666 .

Ya desarrollado el problema nos quedaría así:

Por lo tanto obtenemos:

P(x = 4) = 0.026

Con esto queremos decir que tiene una probabilidad de 2.6% de obtener 4 veces el número 3 al tirar un dado 8 veces.

La distribución de Poisson parte de la distribución binomial:

Cuando en una distribución binomial se realiza el experimento un número n muy elevado de

veces y la probabilidad de éxito p en cada ensayo es reducida, entonces se aplica el modelo

de distribución de Poisson:

En donde tiene que cumplir con lo siguiente:

p < 0.10

p * n < 10

En el cual se sigue el siguiente modelo:

Donde:

Tenemos que e es una constante que es igual a 2.71828, lambda es igual a número de veces n que se realiza el experimento por la probabilidad p de éxito en cada ensayo, y k es el número de éxito cuya probabilidad se está calculando.

Ahora veamos el ejemplo:

La probabiliad de que un niño naza pelirrojo es de 0.012, ¿Cuál es la probabilidad de que entre 800 recien nacidos haya 5 pelirrojos?

Podemos identificar que p es menor que 0.1 y el producto de n * p es menor que 10, entonces aplicamos el modelo de distribución de Poisson.

Por lo tanto:

P(x=5) = 4.602

Entonces la probabilidad de que haya 5 pelirrojos entre 800 recien nacidos es del 4.6% .

Distribuciones Continuas

Estos pueden tomar cualquiera de los infinitos valores existentes dentro de un intervalo.

Ejemplo:

Mantenimiento de alguna computadora, en las que interviene el tiempo que existe hasta que se vuelva a descomponer.
Crecimiento relativo anual de alguna población de un país o región.
El periodo de incubación de una enfermedad, ya que una enfermedad nunca es exacta aún si está es la misma.
Nivel de líquido de un embalse.
Cantidad de decimales que se pueden encontrar entre 0 y 1 de una regla; siendo decimales los números después de la coma.

En las distribuciones continuas se destacan las siguientes:

La distribución uniforme es aquella que puede tomar cualquier valor dentro de un intervalo, todos ellos con la misma probabilidad, se dice que es una distribución continua porque puede tomar cualquier valor y no únicamente un número determinado.

La función de densidad es aquella que nos permite conocer la probabilidad que tiene cada punto del intervalo, y viene modelado de la siguiente manera:

En donde b es el extremo superior(si hablamos del 1 al 100, 100 sería nuestro valor extremo superior) y a es el extremo inferior(si hablamos del 1 al 100, 1 sería nuestro valor extremo inferior).

El valor medio de dicha distribución se calcula de la siguiente manera:

Ejemplo:

El volumen de precipitaciones estimado para el próximo año en la cuidad de Sevilla va a oscilar entre 400 y 500 litros por metro cuadrado.Calcular la función de distribución y la precipitación media esperada:

Es decir, que el volumen de precipitaciones esté entre 400 y 401 litros tiene un 1% de probabilidades; que esté entre 401 y 402 litros, otro 1%, etc.

El valor medio esperado es:

Donde la precipitación media estimada en Sevilla para el próximo año es de 450 litros.

La distribución normal o gaussiana es el modelo de distribución más utilizado en la práctica, ya que multitud de fenómenos se comportan según una distribución normal, dicha distribución se caracteriza porque los valores se distribuyen formando una campana de Gauss, en torno a un valor central que coincide con el valor medio de la distribución:

Donde hay un 50% del lado izquierdo y derecho.Está distribución viene definida por dos parámetros:

X: N(m, sˆ2)

Donde m es el valor medio de la distribución y es precisamente donde se sitúa el centro de la curva, sˆ2 es la varianza, esto indica si los valores están más o menos alejados del valor central, si la varianza es baja los valores están próximos a la media y si es alta, entonces los valores están muy alejados de ella, se representa por sˆ2 por que su raiz cuadrada, s, es la denominada desviación estandar.

Cuando la media de la distribución es 0 y la varianza es 1, se denomina normal tipificada, y su ventaja reside en que hay tablas o rutinas de cálculo que permiten obtener esos valores, donde se recoge la probilidad acumulada para cada punto de la curva de esta distribución.

El modelo de la distribución normal es el siguiente:

Ahora ya explicado lo anterior eligiré una distribución binomial, ya que quiero saber la probabilidad de que nos salga un número impar o par y necesitamos mas de una prueba.

Problema planteado

¿Cuál es la probabilidad de obtener 5 veces el número 6 al lanzarlo 30 veces?

Primero tenemos que identificar la probabilidad que tenemos de conseguir un numero par:

El cual será 1/6, por lo tanto p vale 0.1666, n es igual a 30 y k es igual a 4 que es las veces que deseamos obtener dicho resultado.

Para ello usaremos el siguiente modelo:

Por medio de octave se hizo una función en la que se usa está formula:

function binomial(n, p, filename)
  q = 1 - p;
  output = fopen(filename, "w");
  for k = 0 : n
    fprintf(output, "%d %f\n", k, bincoeff(n, k) * p**k * q**(n - k));
  endfor
  fclose(output);
endfunction

Donde k es el número de intentos y bincoeff(n, k) es una forma simplificada de octave que nos permite resolver el coeficiente binomial, a lado de ello multiplicamos p elevada a la k por 1-p que en este caso es q elevada a la n-k.

Luego en octave lo usamos de la siguiente manera(siempre y cuando estes en la misma dirección del archivo):

binomial(30,0.1666,"lanzamientos.dat")

Donde el primer valor es n número de veces, el segundo valor es la p probabilidad de éxito y por último el nombre que deseamos del archivo.Y como resultados obtuve lo siguiente:

Ahora si vemos que el resultado de obtener 5 veces el número 6 al tirar 30 veces el dado nos da como resultado 0.1920 por lo tanto hay un 19.20% de probabilidad de obtener ese resultado.

Ahora generamos nuestra gráfica con Gnuplot:

Ahora deseamos hacerla continua, viendo las distribuciónes anteriores la única que se me hizo mas semejante al momento de generar las gráficas fue la normal, ya que me genera una campana gaussiana.

Para ello debemos calcular nuestra media y desviación estándar.

Donde la media es igual a n * p donde n sería el total de veces que quieres que aparezca y que p sería la probabilidad.

Ahora nuestra desviación estándar donde necesitamos sacar la raíz cuadrada n p q donde q es la resta de 1- p.

Luego determinamos el valor de Z que es el número de desviaciones estándar de x respecto a la media de distribución.Donde su formula es la siguiente:

Donde x es el valor de la variable aleatoria, np es la media de la distribución de esta variable aleatoria y s es la desviación estándar de esta distribución.

Ahora si lo resolvemos nos daremos cuenta que 0.00097995 nuestra Z, ahora hay que comprobar que nuestra gráfica puede ser aproximada con nuestra distribución binomial en donde por medio de una función normcdf en Octave determinamos la probabilidad de obtener 5 veces el número 6 al lanzarlo 30 veces, en donde obtenemos 0.50039 por lo tanto tenemos un 50.03% de probilidad de obtener dicho resultado, comparandolo con nuestro resultado de la distribución binomial vemos que tenemos 0.61683 por lo que tenemos un 61.68% de probabilidad.

Aquí vemos la gráfica con gnuplot de la distribución normal.

Ahora comparamos.

Para llegar a dichos resultados me base en el código de mi compañero Juan Carlos, ya que tuve ciertos problemas al momento de aproximar dichas distribuciones.

function Normal(n, p)
distBin = fopen("binomial.dat", "w");
distNor = fopen("normal.dat", "w");

#Determinamos nuestra media, desviacion estandar y q
media = n * p;
q = 1.0 - p;
desviacion_estandar = sqrt(n*p*q);
k = 0.0;

while(k < n)
    #Determinamos Z
  z = (k - media) / desviacion_estandar;
    # Con esta funcion creamos grafica de una distribucion normal
normal = normpdf(z, 0, 1);
    # Ajustamos nuestra grafica en base a la altura de las lineas
y = desviacion_estandar * (bincoeff(n, k) * p**k * q**(n - k));
    #Se ajusta la altura con refencia a la distribucion normal
fprintf(distBin, "%f %f\n",z , y);
    # Se imprime nuestra Distribucion normal
fprintf(distNor, "%f %f\n",z , normal);
k = k + 0.1;

  endwhile
  fclose(distBin);
  fclose(distNor);

endfunction

Links referencia:

www-users.math.umd.edu

http://nutriserver.com/Cursos/Bioestadistica/Distribuciones_Continuas.html

http://nutriserver.com/Cursos/Bioestadistica/Distribuciones_Discretas.html

http://www.gestiopolis.com/recursos/experto/catsexp/pagans/eco/44/distrinormal.htm

Caso de modelo determinístico: Elevador

2012-02-09T11:51:00.000-08:00

Para la primer entrada de este blog voy a comentar sobre ¿Qué es un modelo determinista?, y explicaré un caso de modelo determnista, elegí como ejemplo un elevado.

¿Qué es un modelo determinista?

Según la wikipedia un modelo determinístico es un modelo matemático donde las mismas entradas producirán invariablemente las mismas salidas, no contemplándose la existencia del azar ni el principio de incertidumbre. Está estrechamente relacionado con la creación de entornos simulados a través de simuladores para el estudio de situaciones hipotéticas, o para crear sistemas de gestión que permitan disminuir la incertidumbre.

Ahora según mi conocimiento un modelo determinístico es un modelo donde siempre habrá una salida esperada en la que no hay una variabilidad de valores aleatorios.

Ya que tenemos definido algo, explicaré sobre el caso que elegí que es el de un elevador el cual es un sistema de transporte que sirve para movilizarnos a ciertos niveles aunque en la mayoría de los casos lo conocemos como pisos.

Elegí el elevador por que tiene los aspectos determínisticos y según lo que definimos tiene salidas definidas y ninguna es aleatoria, entre los cuales se pueden mencionar las siguientes acciones del modelo:

Presionar botón.
Abrir puerta.
Presionar botón(elegir piso).
Cerrar puerta.
Ir a piso.
Abrir Puerta.
Si hay pisos en cola ir al siguiente, si no, quedarse en la misma.
Cerrar Puerta.

Una de las cosas que se observaron fue la secuencia de instrucciones, de las cuales tienen que seguir dicho patrón para que funcione como tal.

Código

class Elevador():
 def __init__(colaPisos,posicionActual):
  self.colaPisos=colaPisos
  self.posicionActual=posicionActual
 
 def run(self):
  while True:
   direccionElevador = colaPisos[1]
   if(direccionElevador > posicionActual):
    moverArriba(direccionElevador)
    colaPisos[1].pop
   else:
    moverAbajo(direccionElevador)
    colaPisos[1].pop
   break

 def moverAbajo(self,direccion):
  self.posicionActual=direccion
  print'Posicion Actual ',self.posicionActual
  print'Cerrando Puerta'
 
 def moverArriba(self,direccion):
  self.posicionActual=direccion
  print'Posicion Actual ',self.posicionActual
  print'Cerrando Puerta'

def main():
 print 'Presiona boton'
 m = raw_input()
 while m != '':
  print'Abre puerta'
  print 'Selecciona un piso'
  n = raw_input()
  colaPisos = []
  colaPisos.append(n)
  posicionActual = [5]
  print'Estas en el piso ',posicionActual
  elevador = Elevador(colaPisos,posicionActual)
  elevador.start()

Links de referencia:

Definición corta de determinista wikipedia

es.wikipedia.org/wiki/Modelo_determinístico

Imagen del elevador

http://blog.wegowise.com/Portals/107469/images/