Arboricultura: parte 1

Leí algo que había leído hace como un millón de años: los arboles de Huffman. La pregunta que responde es simple: dado un texto formado por símbolos, ¿cual es la forma de codificar cada símbolo para minimizar el tamaño final del texto codificado? Para responder esto, se arma un árbol binario, de abajo para arriba. Pensando en los tipos de árboles que se pueden generar, me intrigó uno en particular: Desde la raíz para abajo, todos los nodos tienen: una hoja y una bifurcación. (mas abajo hay dibujitos donde queda claro). Como se genera este tipo de arboles? Si \(\left\vert{\text{simbolo}_i}\right\vert > \sum_{j > i}\left\vert{\text{simbolo}_j}\right\vert\), entonces vamos a generar ese tipo de árbol. En particular, si  \(\left\vert{\text{simbolo}_i}\right\vert = 2^{n-i}\), vamos a estar observando ese tipo de árbol. O, si vemos ese tipo de árbol, podemos hacer una transformación a un texto donde los símbolos tienen esta relación de ocurrencia: existe una ley de potencia detrás. Entonces, refinando la pregunta: sirve analizar la topología del árbol de Huffman para detectar leyes de potencia? Existe la definición de distribuciones Fat-tailed , que no es la misma que Heavy-tailed. El tema es que esta definición pide que exista una ley de potencia en la cola de la distribución. No hace falta que sea en toda la distribución. Entonces, podremos detectar este tipo de distribuciones usando los árboles? No lo sé, pero creo que sí :)

Para probarlo hice lo siguiente:

• generé muestras de tamaño \(n = 100000\) de varias distribuciones. 
• de cada distribución tomé el histograma usando nBins (\(\approx 1000\) )
• Usando este histograma armé el árbol de Huffman, pero solo tomando los bins con al menos una ocurrencia

Existen dos casos patológicos: la (predecible) uniforme, mostrando un árbol binario perfecto y la distribución de plank, que es un perfecto ejemplo del tipo de árbol que quería encontrar.

 

(para la uniforme utilicé solamente 10 bins porque de otra forma necesitaba hacer un gran zoom para notar la estructura).

Las distribuciones Heavy tailed muestran la estructura similar a la de Plank (sin ser tan extremo), pero ademas, hay distribuciones que tienen (en sus colas) una estructura similar a la de Plank. Dos casos que me llamaron la atención son: la Geometrica y la Poisson

Posiblemente esto esta relacionado con lo que dice (not-so)weak-pedia:

"Fat tail distributions have power law decay in the tail of the distribution, but do not necessarily follow a power law everywhere"

Puede ser que sean solo propiedades de muestra finita.

Ahora queda esta pregunta: ¿como detectar numéricamente esta estructura analizando el grafo? Sería agradable encontrar un estadístico (una variable aleatoria con una distribución de probabilidad) sobre la cual hacer test de hipótesis. No encontré nada de eso.
Probé varias medidas, pero ninguna me convence. Las que pongo acá son las siguientes:

• el número de Strahler: mide la complejidad de las ramas. Los resultados son interesantes. Va de \(1\) a \(\lceil\log_2 \text{nBins}\rceil\). La uniforme tiene el número mas alto
• Bifurcation Ratio (bf): En teoría hay un Teorema Central del Limite sober este estadistico y arboles aleatorios, pero los resultados no fueron para nada buenos
• \(\psi\): Es el promedio de la distancia de los leafs al root dividido por \(\lceil\log_2 \text{nBins}\rceil\) (la cantidad de bits necesaria para representar a los leafs)

resultado para nBins = 1025

*** normal ***

bf ratio: 1.30154390424

strahler: 9

distance (psi): 1.222749 upperBound: 41.100122

*** uniform ***

bf ratio: 1.82279829545

strahler: 11

distance (psi): 1.004502 upperBound: 47.727359

*** weibull ***

bf ratio: 1.33107640276

strahler: 9

distance (psi): 1.233694 upperBound: 41.600120

*** Heavyweibull ***

bf ratio: 1.25439814815

strahler: 6

distance (psi): 1.665082 upperBound: 11.500679

*** pareto ***

bf ratio: 1.2426702689

strahler: 6

distance (psi): 1.843468 upperBound: 13.875563

*** lognormal ***

bf ratio: 1.26271929806

strahler: 6

distance (psi): 1.539683 upperBound: 19.111435

*** lognormalHeavy ***

bf ratio: 1.25248954477

strahler: 5

distance (psi): 1.705980 upperBound: 6.144518

*** vonmises ***

bf ratio: 1.34570436657

strahler: 9

distance (psi): 1.279363 upperBound: 39.300127

*** expovariate ***

bf ratio: 1.32432288929

strahler: 8

distance (psi): 1.311187 upperBound: 36.700136

*** zipf ***

bf ratio: 1.16205357143

strahler: 4

distance (psi): 2.105769 upperBound: 4.336538

*** planck ***

bf ratio: 1.09090909091

strahler: 2

distance (psi): 1.750000 upperBound: 1.522727

*** burr ***

bf ratio: 1.22805961483

strahler: 6

distance (psi): 1.878981 upperBound: 9.875796

*** cauchy ***

bf ratio: 1.29671439671

strahler: 5

distance (psi): 1.650579 upperBound: 5.287645

*** levy ***

bf ratio: 1.25

strahler: 3

distance (psi): 1.250000 upperBound: 1.277778

*** geom ***

bf ratio: 1.125

strahler: 3

distance (psi): 1.965517 upperBound: 3.006897

*** t ***

bf ratio: 1.25474980384

strahler: 6

distance (psi): 1.484777 upperBound: 15.667062

*** poisson ***

bf ratio: 1.11538461538

strahler: 3

distance (psi): 1.487500 upperBound: 1.612500

*** gamma ***

bf ratio: 1.31289744848

strahler: 8

distance (psi): 1.257250 upperBound: 33.500149

Acá va el resto de los grafos

El proximo post, va a compartir con este solo las hormigas (espero)

[EX-POST: como los árboles de Huffman están relacionados con la entropía, y esta es una propiedad de las distribuciones, puede ocurrir que la estructura sea un invariante de las distribuciones. Esto sería muy bueno. Por ejemplo, se podría determinar la distribución (o familia de) mirando la estructura. Supongo que hay que hacer bastante trabajo para demostrar esto]

apéndice

parámetros de las funciones

  normal(0, 1), uniform(0, 1) weibull(1, 1.7) weibull(1, 0.3) (heavy) pareto(2.5)  lognormvariate(0, 1) lognormvariate(0, 2.5)(heavy) vonmises(0, 4) exponential(4) zipf(3) planck(1) burr(2, 1) cauchy() levy(c=4) geom(0.3) t(4) poisson(4) gamma(2, 2)

La culpa no es de los floats

Sino, del que le da de comer.

Al final del ultimo post, comenté mi preocupación por la lenta convergencia, y culpe apresuradamente al uso de los objetos de tipo float. 
Haciendo cuentas noté que no había nada raro en la convergencia.

El estimador era:
\[\tilde{\pi} = 4 * \frac{\sum^n_{t=1}I(\text{punto i esta dentro de circulo})}{n}\]
Esto se puede reescribir de la siguiente forma:
\[\tilde{\pi} = 4 * \tilde{F}\], donde \[\tilde{F} = \frac{\sum^n_{t=1}I(\text{1 si punto i esta dentro de circulo})}{n}\] Sabemos que \(I(\text{1 si punto i esta dentro de circulo})\) es una variable aleatoria Bernoulli, entonces
\(\tilde{F} = \frac{Y}{n}\) donde Y es una variable Binomial con parámetros \(n\) y \(p = \pi/4\).
\[E[\tilde{\pi}] = 4 * \frac{E[Y]}{n} = 4 * \frac{\pi}{4} * n * \frac{1}{n} = \pi\]

Ahora, la varianza es lo que nos interesa.
\[\begin{aligned}Var[Y] &= np(1-p)\\
Var[\tilde{F}] &= \frac{np(1-p)}{n^2} = \frac{p(1-p)}{n} = \frac{\tilde{F}(1-\tilde{F})}{n}\\
Var[\tilde{\pi}] &= 4^2Var[\tilde{F}] = 8\frac{\tilde{F}(1-\tilde{F})}{n}
\end{aligned}\]

Ahora, quiero generar un intervalo de confianza del \(95\%\) (usando la normal, porque \(n\) es muy grande), y quiero que ese intervalo tenga un tamaño del orden de \(10^{-4}\). Teniendo la formula de la varianza del estimador, podemos calcular el \(n\) para esa precisión: \(n = 1035989708.4897318\). Acá queda claro que floats no tenía nada que ver con la falta de precision. Fue solo una cuestión de ansiedad.
De hecho, si usamos \(n=10^9\), luego de un tiempito, tenemos el siguiente resultado


\[\tilde\pi = 3.141564804\, (3.67207693157*10^{-5})\]
y con un \(\alpha = 5\%\) podemos construir el intervalo: \( [3.14149283129\; 3.14163677671]\)

Pi con Montecarlo

Hablando con tenuki sobre como usar montecarlo (ese pariente cercano a bootstrapping) para estimar objetos matematicos no aleatorios, me comentó que, por ejemplo, se podía estimar pi. De hecho, resultaba bastante facil hacerlo.
Imaginemos un cuadrado de lado 2, y dentro un circulo de radio 1. También pensemos en puntos que se distribuyen según una distribución uniforme dentro de ese cuadrado.
La probabilidad de que un punto esté dentro del circulo es pi/4 (la integral en el recinto definido por el circulo de f(x) = 1/area_del_cuadrado).
Una forma de estimar esta probabilidad es contando los puntos que estan dentro del circulo y dividir por n, lo que tiene como esperanza pi/4
Entonces, 4*puntos_dentro/total_puntos = pi

Es notable que la convergencia no es monotona, y puede tardar bastante en lograr buena precisión. Supongo que el uso de floats tiene algo que ver en el asuntó (en limitar la uniformidad de los puntos).

 

Generación de variables normales con una matriz de covarianzas

Este es un problema que tuve hace un tiempo, que lo solucione de otra forma (no me acuerdo como). Hace poco se me apareció esta solución, y me pareció interesante.

Nosotros tenemos una matriz de NxN de covarianzas (C). Queremos generar m samples aleatorios de N normales que respeten esa matriz de covarianzas.

Como la matriz C es real y simétrica, podemos encontrar una eigen decomposition:

C = RDR^tr

donde R es una matriz orthogonal formada por los autovectores, y D una matriz diagonal con los autovalores en la diagonal.

Por otro lado, si tenemos m realizaciones de n variables normales aleatorias independientes (con media cero) en una matriz X de mxn, la estimación de la matriz de Covarianzas sería:

C = 1/mX^trX.

Si las n variables son independientes, entonces, 1/mX^trX. va a ser una matriz diagonal, y el elemento {i,i} de la diagonal va a corresponder con la varianza de la variable i-esima.

Ahora, si generamos m realizaciones de n variables aleatorias, y cada variable aleatoria tiene media cero, y varianza igual al i-esimo autovalor, tenemos que

D = 1/mX^trX

Por lo que, si definimos a

U = XR^tr

tenemos una matriz de n variables con m realizaciones, y con matriz de covarianzas igual a

= 1/mU^trU

= 1/mRX^trXR^tr

= R( 1/mX^trX)R^tr

= RDR^tr

C = RDR^tr.

Que era lo que estábamos (bue, estaba buscando yo, no se si todos estábamos buscando eso) buscando.

Este metodo es super costoso si la matriz de covarianzas cambia siempre -por lo costoso de calcular los autovalores-, pero realmente sirve para entender por que PCA funciona. En realidad, habría que tener en cuenta los residuos, espero en algún momento del próximo año hacer un post sobre eso

Feliz año

Web Semantica

Esto lo pensamos con amigos hace varios años, y lo encontre por casualidad... Puede ser cualquiera, lo tendria que releer :)

Modelo

Primero definimos semántica: Una semántica para una palabra determinada es la forma en que se relaciona esta palabra con un conjunto de palabras.

En Internet existe en cada punto de la red una semántica propia para ese punto. Esta semántica esta dada por las conexiones que tiene una pagina hacia otras, por los contenidos que tiene esta pagina, por la forma en que están organizados internamente los contenidos, por las paginas que linkean a una pagina, etc. Esta semántica existe, y no podemos alterarla, así que la podemos definir como semántica f'áctica.

También existen agentes (supongamos), estos agentes tienen intereses. Estos intereses los podemos representar con el significado que le dan a distintas palabras. Entonces los agentes crean significados (son significantes), y sus significados son propios. Podemos decir que tienen una semántica subjetiva. Un ejemplo de esto seria que para algunas personas la palabra c# significa una nota (do sostenido), para otras personas significa un lenguaje de programación. Acá, implícitamente se esta armando una red entre distintas palabras para armar el significado de un símbolo (ver http://en.wikipedia.org/wiki/Semantic_networks ). La información va a ser interesante para el agente en tanto responda al orden que define en su semántica subjetiva (por construcción). Entonces, podríamos decir que para que el agente tenga una buena experiencia de navegación por Internet, la información tendría que estar organizada según su semántica.

Sabemos que el agente, cuando navega por la red, necesariamente esta moviéndose por un mundo donde la información esta organizada según la semántica fáctica, sin embargo, la navegación la hace con el mapa de la semántica subjetiva. El agente va a estar tratando de maximizar todo el tiempo la proximidad con su semántica, y eso se va ver, por ejemplo, cuando decide hacer clic en un determinado link, o cuando decide ignorarlo. Por ejemplo, si alguien va a clarín.com, la información esta organizada segundo lo que a clarín le parece, pero esta persona va a leer lo que a esta persona le parece que son noticias. Por ejemplo, si a esta persona le parece que los pozos que hay en Palermo no son noticias, no va a pedir mas información al respecto. Y si esta persona tiene un interés por las noticias internacionales, va a hacer clic directamente en la sección de internacionales.

Entonces, el problema se puede resumir a conseguir una función que mapeé de la semántica fáctica a la subjetiva (es decir, que pueda transformar el orden que tiene la información en la red, en un orden que le produce interés al agente).

Supongamos que esta función existe 

Como existe esa función, vamos a aproximarla usando redes neuronales.

Estas redes van a ser tan grandes, que cada agente va a ser una simple neurona. Esta neurona va a estar relacionado con otros neuronas. La fuerza que los une va a representar lo cerca que van a estar las semánticas para algunas palabras. Incluso van a existir relaciones de segundo orden. La neurona A esta conectado con B que a la vez esta conectado con C. La fuerza con la que A esta conectado con C depende de la relación que tiene A con B y B con C. Cuando B nota que la relación A(B)C es suficientemente fuerte, puede presentar C a A, entonces se creara la relación AC. De esta forma, se va a estar creando una forma de aproximar las función de traducción semántica, usando las relaciones entre los agentes.

Esta nueva red va a tener una topología particular, y distinta de la topología fáctica y las topologías subjetivas, es decir, va a tener un significado propio, emergente de la relación entre los agentes

.

Una aplicación alternativa.

Una aplicación nueva que se le puede dar a este modelo es la enseñanza.

Una de las cosas que me parece a mi que significa enseñar, es transmitir un significado para un grupo de símbolos. (Crearle significado, o resignificarlos). Entonces, me imagino a un profesor universitario diciéndole a su grupo de alumnos: “agréguenme a su lista de contactos mandatory”. En esta situación, lo que esta haciendo el profesor es diciéndole que el va a ser el significante, el que va a cargar de significado la información que busquen los alumnos. De esta manera, el profesor puede transmitir una manera de ver el mundo a los alumnos, transmitiéndole la semántica que el supuestamente tiene.

Una aplicación similar se puede dar en un equipo de research, o en algún equipo de trabajo. Cuando empieza una persona nueva a trabajar con el equipo, el tech lead puede ofrecerse como el significante, de esta manera, cuando la persona nueva busque algo en la red, la información va a salir de los lugares a los que el tech lead acude cuando necesita información, pudiéndole transmitir ese conocimiento.

El comportamiento de esto es interesante por que, si bien el significante es definido, los símbolos a los que se les quiere dar significado van a estar ingresados por el agente final, y van a ser distintos por su subjetividad a lo que podría el significante haber pensado, lo que produce cierto significado emergente.

TODO:

Pensar sobre privacidad

Como se implementa? Simplemente indexando todo?

Un regalito

K-means y gravedad

Espero que este articulo cierre con un moño azul

 

Disclaimers: Articulo computronico. Drito, no te quejes

Mirando un poco una implementación de k-means me di cuenta que uno de los problemas que tenia era que los puntos no tenían “conciencia” de su posición en el espacio ni de su entorno. Básicamente, un punto era un placeholder de coordenadas, y nada mas.

Como resolver este problema era algo que empezó a dar vueltas por mi cabeza. No llegue a ninguna solución (me adelanto, y bajo las expectativas sobre este post). A pesar de eso, me di cuenta de que el universo esta todo el tiempo haciendo k-means. Básicamente, la gravedad es una función de la distancia entre los cuerpos (y también de la masa, pero gracias al 1, podemos olvidarnos momentáneamente de eso).

Los cuerpos en el espacio “sufren” y la atracción gravitacional de otros cuerpos (no toman conciencia de esa atracción, pero tiene un impacto sobre estos cuerpos). Entonces, el modelo que planteo es pensar a los distintos puntos como distintos cuerpos (a priori, con masa igual a 1), y en determinar los clusters en base a los conjuntos donde las relaciones gravitatorias son mas fuertes (de la misma manera que existen planetas, sistemas planetarios, sistemas de estrellas, galaxias, etc.). Ahora se produce un cambio de la relación, en vez de ser la relación de cada punto contra la media teórica de los clusters, la relación es entre cada unos de los cuerpos (creo que es un poco lo que se hace en fuzzy k-means).

Hasta ahora este cambio de metáfora no aporta nada (es mas, agranda la complejidad algorítmica, lo que no es bueno). Sin embargo, creo que se pueden hacer agregados de cuerpos. Esto seria, definir un Epsilon, y en base a ese Epsilon crear esferas, cuya densidad depende de la cantidad de puntos de esta esfera y del tamaño del Epsilon (volvemos a algo parecido a k-means). Ahora, estos cuerpos pasan a ser iguales a el resto de los cuerpos (no agregados), y solamente difieren en la densidad (y acá si podemos empezar a utilizar feature propios de la gravedad).

Ahora aparece un problema con la definición de un Epsilon. Si hacemos una esfera de un radio determinado, estamos suponiendo implícitamente que los puntos van a estar distribuidos uniformemente en las n dimensiones. Este supuesto es demasiado fuerte. Quizás, algo que se podría hacer es hacer una ponderación en las dimensiones donde los puntos tienen mayor variabilidad. Bueno, básicamente esto seria proyectar todos los puntos en el espacio de los autovectores de la matriz de varianza/covarianza. Es decir, utilizar PCA (sin perdida de información).

Una vez que realizamos esta proyección, nos queda una topología donde, en principio, tiene mas sentido utilizar una esfera para crear los agregados.

 

Creo que una de las ventajas que puede tener esta metáfora es con respecto a las proyecciones.

Pensando por ejemplo en el tiempo. Si el tiempo no es una variable endógena (no es una de las n dimensiones), y quisiéramos hacer una extrapolación. Es decir, quisiéramos hacer una proyección sobre la posición futura de los distintos cuerpos, creo que tiene sentido estadístico utilizar los mismos modelos que se utilizan en la física para establecer la posición futura de un cuerpo en movimiento (es decir, para donde caen las manzanas). Creo que los datos pasados nos podrían ayudar para calcular la velocidad de los cuerpos, y quizás nos ayudarían a detectar “agujeros negros” en nuestra muestra.

De cualquier manera, estas son especulaciones teóricas, recién tengo planeado cursar física del CBC el cuatrimestre próximo, ahí quizás tenga alguna idea mas acabada sobre este tema

(algo de ) Lo que no me enseñaron sobre la correlacion

Ante todo, buenas tarde

Hace un par de meses, Marianita me paso un paper sobre analisis de lenguaje natural, o algo por el estilo. 

Lo interesante (para mi) de ese paper es que utilizaban como medida de distancia de dos variables, el coseno del angulo que forman los vectores de sus realizaciones (normalizadas).

Esto realmente no tenia demasiado sentido para mi. 

         1.         Angulo entre vectores en espacios n-dimensionales?! Donde estaba el cateto opuesto y la hipotenusa (ok, podian ser las soluciones a ecuaciones diferenciales, pero no me simplificaba el panorama)

         2. Que era lo que media el coseno ahi?

Bueno, el punto uno me costo un tiempo entenderlo (el dos mas), pero mas o menos la idea es asi:

Suponga que queremos ver como en distintos documentos (n documentos) aparecen distintas palabras (m palabras). -si tuve que uzar un ejemplo concreto, por que me estaba haciendo un rollo con la abstraccion-

Para cada palabra, tenemos un vector (en realidad, un punto) con las apariciones en los n documentos, es decir, tenemos un vector n-dimensional por cada palabra. Esto nos da m vectores n dimensionales.

Ahora, queremos sacar el coseno entre dos palabras

Un vector, ademas de ser una n-upla (algo que tiene n elementos),  define un segmento (y una recta) que pasa por el 0 (si, claro, un vector... :) ) 

Entonces, nosotros teniamos dos vectores, es decir, que por ahi oculto, y bastante callado, habia otro punto inportante (que a priori, es el unico que tienen en comun todos los vectores), el 0.

Agregando el 0, pasamos a tener 3 puntos! (de vuelta, de dimension n). Pero aca esta la magia, de la misma manera que en un espacio euclideo, dos puntos definen una recta, 3 puntos definen un plano! (si, obvio, pero…)

Entonces, si vemos los vectores en el plano definido por los 3 puntos (las dos palabras mas el 0), vemos que, por construccion, los vectores estan en el mismo plano, y pasan en (por lo menos) un punto. Es decir, podemos calcular el coseno!!!

Bueno, esto ahora tenia un poco mas de sentido, pero todabia tenia un par de cosas haciendome ruido en la cabeza.

Por ejemplo, cual era el espacio original en el que estabamos trabajando? Era un espacio n-dimensional, haci que tenia que ser el espacio de las posibles realizaciones, donde cada uno de los puntos definia un conjunto complete de realizaciones para una variable (si estaba normalizado, no importaba para que variable). 

La existencia de este espacio me parece un poco rara, e interesane –voy a pensar mas al respecto J -

 

Ahora nos queda el Segundo punto, que es el significado de esta medida.

Si vemos la forma compacta -sin hacer explicita la formula del plano- de calcular el coseno para dos angulos (la copiaria, pero major vallan a wikipedia), nos vamos a dar cuenta que es la misma que la formula de la funcion de correlacion (para valores centrados)!!

Una de las cosas a la que le da sentido esto, es que cuando tenemos variables independientes, los vectores son ortogonales, y el angulo es un angulo de 90º!!

 

Esto esta buenisimo, por que todo lo que dije sobre la creacion del plano en el que aparecen los vectores es algo que hacemos implicitamente cuando calculamos la correlacion!

Me da lastima que despues de tantos años de facultad, eso sea algo de lo que nunca he odio, pero me parece que puede ser interesante pensar un poco sobre eso.

 

(puf, post innecesariamente largo, como este fin de semana…)

saludos

/jb

 

Accion farmacologica

"el pantopazol inhibe selectivamente la bomba de protones localizada en la membrana apical de las celulas parietales de la mucosa gastrica. Este bloqueo irreversible provoca una inhibicion mas pronunciada de la secrecion acida en comparacion con los antagonistas de los receptores histaminergicos H2.
En el interior de la celula parietal el pantoprazol alcanza su forma activa, una sulfenamida ciclica cationica, que luego se une a los residuos de cisteina de la bomba protonica, inhibiendo asi la accion de la misma"
y esto se puede conseguir en cualquier farmacia...

finance google no funciona tan bien

disclaimer: Este no es un post computronico

Y eso me pone contento :)

Le estaba mostrando a un amigo (juan pablo) finance.google.com, y el descubrio un pequeño bug.

Por curiosidad, quizimos ver como andaba tenaris ( http://finance.google.com/finance?q=NYSE%3ATS ), y al ver el Market Cap, me dijo: "che, eso no esta bien. el numero es 111.56B, pero esta mal", de hecho, en http://finance.yahoo.com/q?s=ts, el Market Cap era 22.31B, es decir, una pequeña diferencia de 89.25B, casi nada, digamos :)

El me dijo que no estaban teniendo e cuena que tenaris era un ADR, pero no pude encontrar otro titulo donde el error sea tan grosero.

Realmente me alegra que a la gente de google no le salgan siempre tan perfecto las cosas (si, es envidia :) )

Por si no se entiende, donde dice B se lee Billions, que es mil millon dollars (mucha mucha plata)

kmeans y c++ template metaprogramming

Estoy haciendo una implementacion de kmeans tratando de hacer la maxima cantidad de cosas posibles con templates metaprogramming, para transferir en el codigo el conocimiento que ya tengo en tiempo de compilacion.
Me aparecio este petit problema (gracias g++, sos lindo... por suerte ya lo resolvi):
iChing:~/.../code/kmeans jb$ make maing++ main.cpp -o maini686-apple-darwin8-g++-4.0.1: Internal error: Illegal instruction(program cc1plus)Please submit a full bug report.See <URL:http://developer.apple.com/bugreporter> for instructions.make: *** [main] Error 1

Information Entropy y k-means

Hace un par de semanas me encontre con la medida de la entropia de Shannon. Es simplemente esta formula:

I = -∑p_i*ln₂(p_i) donde p_i = N_i / N

N es la cantidad de elementos en todo el conjunto, y N_i es la cantidad de elementos en cada uno de los k i-subconjuntos

Parece una formula simple, pero en realidad, simplemente mirandola, no podia entender su significado.

En principio parece una esperanza, en particular, la -E(ln₂p_i), esto seria algo asi como la cantidad de bits necesarios para representar, en promedio, los elementos existentes en cada uno de los i grupos.

Pero seguia sin decirme demasiado, asi que gracias a la ayuda de Manuk, empezamos a destripar esta formula.

Partimos de:

I = -∑p_i*ln₂(p_i) 

= -∑p_i*ln₂(N_i/N)

= -∑p_i*[ln₂(N_i) - ln₂(N)]

= ∑p_iln₂(N) - ∑p_iln₂(N_i)

como ln₂(N) no depende de i, es una constante, y puede ser extraido de la suma.

= ln₂(N)*∑p_i - ∑p_i*ln₂(N_i)

como ∑p_i = 1, y el segundo termino es una esperanza, nos queda:

=ln₂(N) - E(ln₂(N_i))

El primer termino es constante, siemre que se mantengan constantes las cantidades de elementos del conjunto, entonces, si queremos maximizar o minimizar la entropia, tenemos que modificar el segundo termino (la distribucion de los elementos dentro de los conjuntos).

De hecho, la entropia va a ser maxima cuando cada conjunto tenga exactamente la misma cantidad de elementos. Es decir, cuando los N elementos esten distribuidos entre los k conjuntos. 

Cuando eso ocurra vamos a tener que N_i = N / k

Entonces,

I = ln₂(N) - ∑N/k*ln₂(N / k))

= ln₂(N) - ∑N/k*[ln₂(N) - ln₂ (k))]

= ln₂(N) - [ln₂(N) - ln₂ (k))]

= ln₂(N) - ln₂(N) + ln₂ (k))

=ln₂(k)

es decir, a la cantidad de bits necesarios para representar a todos los conjuntos

Todo esto con respecto a la entropia, ahora vayamos a k-means.

k-means , es un algoritmo para agrupar un conjunto de datos en k subconjuntos, minimizando alguna distancia (en general, la distancia euclidea).

la funcion a minimizar es la siguiente

V = ∑_i∑_jdist(x_ij;c_i)

donde c_i es el centro del subconjunto.

Si utilizamos el error cuadratico medio, lo que vamos a estar minimizando es esperanza de la varianza intragrupos (es decir, dentro de los grupos)

Esto de encontrar grupos que tengan algun sentido (minimizen alguna distancia), parecia que tenia bastante que ver con la teoria de la informacion, asi que empezamos a pensar un poco al respecto (aca vortex y tenuki me ayudaron un poco).

Como la funcion distancia, puede ser cualquier funcion distancia, que cumpla con las condiciones de una funcion distancia, utilize la distancia discreta multiplicada por ln(N_i)

dist(x_ij;c_i) = 0 sii x_ij = c_i

dist(x_ij;c_i) = ln₂(N_i) sii x_ij != c_i

como ninguno de los N_i elementos es c_i,  (si alguno fuese, multiplicamos y dividimos la distancia por algun valor, y el problema se resuelve), entonces,

∑dist(x_ij;c_i) = N_i*ln(N_i)

reemplazando en V, tenemos

V = ∑_iN_i*ln(N_i)

dividiendo m. a m. por N, tenemos

V / N =∑_iN_i*ln(N_i)/N

como N_i = p_i, entonces

V / N =∑_ip_i*ln(N_i)

V / N = E(ln(N_i)

Que es el segundo termino de la medida de la entropia.

Por lo tanto, en el espacio metrico definido por esta funcion distancia, es equivalente minimizar la distancia, a maximizar la informacion (No hay que olvidarse que en la formula de la informacion, el segundo termino aparece restando). Y es ahi donde las dos cosas tienen algo que ver.

Pero sin embargo, este resultado no me dejo completamente contento.

Por un lado, no refleja lo que yo queria (que era mostrar la cantidad de informacion que se concentra en los c_i; de hecho, estos ni siquiera se utilizan).

Ademas, quise utilizar otra funcion de distancia, como por ejemplo, la distancia de Hamming, pero no pude obtener el resultado que queria.

Espero que a alguien que lea este post se el ocurra alguna funcion distancia mas entretenida

saludos

/jb

ps: Tambien tengo que agradecer a tomas, que desde el mas alla (mas alla del meridiano de greenwich), siempre me da una mano con estas cosas

Acoplamiento y No Convexidad

Disclaimer:

• Este articulo tenia ganas de escribirlo hace un tiempo, pero habia un par de cosas que no me cerraban del todo, y queria esperar a que tenga el concepto cerradito y monono. Eso no paso, pero lo que si paso fue que me canse.
  
• Este es un articulo computronico

Lo que queria hacer era relacionar funciones, metodos, clases,  y namespaces con conjuntos.

Supongamos que tenemos una clase que sea una estacion climatica. Esta estacion climatica tiene sensores, como por ejemplo, termometros. 

Nosotros somos consumidores exernos de informacion, y solamente nos interesa lo

 que la estacion climatica nos puede ofrecer.

La imagen anterior vendria a ser una representacion de nuestro mundo.

A travez de los metodos, se producen relaciones entre los objetos.

Por ejemplo, si el consumidor externo quiere saber la temperatura, le pregunta a la estacion climatica la temperatura, se produce una relacion entre el consumidor y la estacion climatica.

Como la clase que define el consumidor tiene como variable de instancia una estacion climatica, la relacion ya esta establecida en la propia clase, por lo tanto, la relacion que genera el metodo es una relacion que va de la clase a la clase. 

A pesar de que por claridad yo dibuje las flechitas como saliendo de las clases, cuando estos nos encontramos con estos metodos, en realidad lo que ocurre es que todas las relaciones se producen dentro del conjunto que define la clase. Lo que se parece bastante al concepto de convexidad.

Por otro lado, puede ocurrir que el consumidor externo, en vez de preguntarle a la estacion climatica la temperatura, le pida su sensor de temperatura, y a la vez a este le pida la temperatura actual.

Esto es lo que vemos en este grafico:

En este caso, la relacion que generaria el metodo para obtener la informacion de la temperatura, iria mas alla de lo que define el conjunto de la clase consumidor externo (aunque volveria al consumidor externo, "salio por un rato" para buscar la temperatura).

Si trazacemos las flechas que define las relaciones producidas por este metodo, veriamos que pasa por la clase "sensor de temperatura", el cual no es conocido por el consumidor externo.

En este caso, las relaciones son relaciones no convexas, y por otro lado, este es un buen ejemplo de alto acoplamiento entre clases.

En que sentido esto no esta piola? 

Bueno, si por ejemplo cambiamos el sensor de temperatura por uno construido en el imperio, con las temperaturas medidas en grados farengein, en vez de en grados celcios, tendriamos un gran problema, ya que el consumidor externo no sabria nada sobre eso (problema que seguro estaria resuelto por la estacion climatica).

Por otro lado, hay algo interesante. A pesar de que estas tres clase estan acopladas entre si, las relaciones que producen sus metodos se mantienen dentro del conjunto definido por esas tres clases.

Si quisieramos definir un namespace que englobe esas tres clases, podriamos hacerlo, y tendriamos las caracteristicas de convexidad para ese namespace. 

Por lo que estariamos encapsulando este micromundo en ese namespace. Y maravillosamente, el minimo conjunto convexo que contiene a una serie de puntos se llama capsula convexa (o clausura convexa), por lo que el nombre de encapsulamiento parece tener aun mas sentido :)

Bueno, espero que a alguien le guste este post, a mi realmente no me gusta, pero ya no queria pensar mas en esto, y esta es la mejor manera de sacar un tema de la cabeza :)

saludos

/jb 

Feedback Positivo y crisis financieras

Ultimamente estuve pensando bastante sobre retroalimentacion positiva y negativa. Los sistemas con retroalimentacion negativa (autoregulados), son particularmente interesantes, y por el otro lado los sistemas con retroalimentacion positiva son los que se ven generalmente en los modelos caoticos y con soluciones fuera del equilibrio. (se van se van se van...)

Con respecto a los modelos financerios, estuve pensando un poco en esto:
Muchos operadores financieros utilizan modelos que suponen eficiencia de los mercardos, para la toma de deciciones.
Supongamos que el mercado es levemente no eficiente, es decir, algunas de las condiciones de eficiencia no se esta cumpliendo (esto es mas facil que suponer que el mercado es eficiente ;) ).
Entonces, los modelos, aplicados sobre este tipo de mercados, van a ser una aproximacion (quizas muy buena, es verdad), pero solo una aproximacion a la realidad.
Como estos modelos son una aproximacion, pueden producir una solucion que no sea la optima. Esto implicaria que haya agentes que esten tomando deciciones suboptimas, y creando ineficiencias en el mercado.
Ahora, estas ineficiencias, se agregan a las que ya existian, generando un sistema aun mas ineficiente, y la solucion de los modelos seria aun mas equivocada, y el sistema se retroalimentaria, divergiendo completamente, hasta llegar a una crisis, u otra solucion de equilibrio (Hysteresis). En esta crisis, los modelos se re ajustan, y tomo comienza de cero...
Ahora, a pregunta es, por que no es que todo el tiempo no estamos llendo vertiginosamente hacia una crisis?
Bueno, creo que parte de las respuestas bienen por varios lados:
a) Puede ser que las ineficiencias del mercado produscan deciciones equivocadas en los modelos, pero la direccion del error puede ser contraria a los errores anteriores, y neutralizar el efecto divergente
b) Animal spirits: Hay agentes que "no utilizan modelos", salvo el de "obtener el mayor resultado". Al no hacer "asumpciones formals" sobre el comportamiento del mercado, aprobechan las deciciones suboptimas de los agentes mas "racionales", y neutralizan esas deciciones
c) Distintos modelos internos producen respuestas contrarias, donde las fuerzas de esas respuestas se neutralizan, y se mantiene un equilibro.

Bueno, ahora, con la existencia de a, b, y c, pasamos de un modelo de crisis permanente, a un modelo de no crisis (que tampoco es cierto). La pregunta es, por que la realidad es algo intermedia?
Bueno, creo que a, b, y c explican un poco el modelo gaussiano-browniano.
El modelo gaussiano, supone, que si hay una cantidad suficientemente grande de errores, todos independientes, y con la misma distribucion (y si, variancia finita), vamos a estar en presencia de una variable distribuida Normalmente (de aca surgen los teoremas centrales en el limite)
El modelo browniano tiene como finalidad, representar el movimiento de una particula que es afectada por la interecaccion sobre esta particula, de una cantidad suficientemente grande (si es infinita, mejor), de particulas aleatorias (tambien, independientes, identicamente distribuidas). El ejemplo tipico, es una pelota, y un monton de focas pegandole a la pelota.. :P (bueno, si, no era el tipico).
Este modelo, parece que cierra bastante con la idea del mercado, donde hay miles de agentes empujando cada uno para un lado particular, y donde cada uno de estos agentes son "independientes", pero (y aca hay un gran pero).
Como vimos, la divergencia producida por los errores de los modelos estan neutralizadas, salvo que aparesca la sincronizacion. Cuando la masa de agentes se sincronizan (por ejemplo, por que todos usan modelos que ante un estimulo en particular, responden en la misma direccion -y ademas, es una decicion suboptima-), la fuerza que pueden hacer el resto de los agentes no llega a neutralizarla, y permite que el sistema salga del equilibrio. Una vez fuera del equilibrio, mas modelos van a ser parte de la sincronizacion (ya no necesitan el estimulo comun, ya que van a responder de forma erronea por la situacion de no equilibrio del mercado). Esto produce el famoso efecto bola de nieve (que lindo ejemplo de retroalimentacion positiva), desencadenando una crisis. Lo interesante, es que por esta sincronizacion, se rompe la idea de modelo browniano, y de normalidad. Ya que ocurren eventos que tendrian una probabilidad de ocurrencia casi nula segun esos supuestos. Es por eso que se utilizan los modelos de colas pesadas (donde los eventos raros, en realidad no son tan raros como parecen)

Bueno, creo que este post me termina de cerrar un par de cavos sueltos
pongo un par de links para expandir el post (si, me dio paja poner href):
http://en.wikipedia.org/wiki/Positive_feedback
http://en.wikipedia.org/wiki/Negative_feedback
http://en.wikipedia.org/wiki/Brownian_motion
http://en.wikipedia.org/wiki/Carl_Friedrich_Gauss (que grande este tipo :) )
http://en.wikipedia.org/wiki/Synchronization_of_chaos

saludos
/jb

Jerarquias (Object Based)

Disclaimer:

Este post es (particularmente) computronico, si quieren no lo lean :)

Historia previa:

El otro dia (asi empiezan todas las historias), con mat y despues con tenuki estabamos pensando sobre el concepto de classe en la programacion basada en objetos. No se si ellos estaban de acuerdo con lo que voy a escribir, pero ayudaron mucho a que estas cosas se organizen asi en mi cabeza :P.

Idea:

Basicamente, pensamos que existen dos tipos distintos de penzar el concepto de clase cuando se programa pensando en objetos.
La primer forma de penzar las clases (y la jerarquia de clases) es desde un punto de vista conceptual. Una clase representa un concepto, y sus subclases representan particularidades de este concepto (la idea de especie y subespecies). Este es un approach Top-Down, ya que primero hacemos el circulo grande donde definimos a la clase, y adentro de este circulo grande hacemos los circulos pequeños que definen a las subespecies. El consumidor de este concepto es un consumidor externo al objeto que es de la clase (o subclase).
Ejemplo: a un mamifero le podemos preguntar la cantidad de tetas que tiene. Eso lo podemos hacer, por que independientemente del mamifero que sea, el concepto de mamifero implica tener tetas.
Como el consumidor es externo, y el concepto es un concepto que se presenta al mundo, tiene sentido que los metodos (los mensajes) de las instancias de esta clase sean publicos, ya que esperamos que desde afuera se utilize este objeto a travez de la representacion conceptual que esta jerarquia esta ofreciendo. Mas aun, desde este punto de vista, tiene todo el sentido del mundo el concepto de classes abstractas (pure concept), y estas clases definen una interface para ofrecer al mundo, pero los elementos que estan en niveles mas particulares (inferiores ?) de la jerarquia tienen que definir como es el comportamiento puntual. Volviendo al ejemplo, sabemos que los mamiferos tienen tetas, pero dependiendo de cada especie la cantidad de tetas que (mas o menos) tienen.
Una de las principales ventajas que tiene esta forma de pensar las clases es la facilidad del uso, yo no necesito saber de que subespecie es tal objeto en particular, ya que puedo utilizar los conceptos que me brinda la especie. Ademas, permite programar comportamiento utilizando metaforas, que es una forma comoda para nuestro cerebro de manejar unos y ceros de forma compleja.


Por otro lado, existe otra forma de pensar el concepto de clases, que es desde un punto de vista mas matematico o funcional. Tenemos un par de objetos (que corresponden a clases distintas), y podemos pensar que puede existir alguna relacion entre estas clases cuando comparten algun tipo de comportamiento. Entonces vamos a agrupar a estas clases en una clase que va a ser "comparten tal comportamiento". Aca el approach es Bottom-Up, ya que primero tenemos los circulitos chiquitos con las clases particulares, y despues armamos un circulo mas grande en donde las agrupamos gracias a este comportamiento compartido. Los usuarios de esta jerarquia son las mismas subclasses, cuando se pretende reusar comportamiento, y agruparlo en un mismo lugar. Me parece que desde este punto de vista, tiene todo el sentido del mundo la existencia de mensajes o herencia "protected", ya que esta jerarquia no ayuda a trabajar con las clases a los consumidores externos, pero permite la reutilizacion de codigo.
Para que este tipo de jerarquia funcione bien, creo que los metodos que se comparten deberian ser metodos donde existe pure behaviour, que esta mas relacionado con el paradigma funcional que con el de Basado en Objetos.
Una de las desventajas de esta forma de ver las clases (me parece), es poder identificar una forma univoca de crear las superclases (un nombre), pero quizas si se pienza como el comportamiento compartido, y no como un concepto compartido, sea mas facil de encontrarlo (al nombre).

Hace un tiempo tenuki me paso un paper sobre la implementacion de traits en smalltalk, y encontre este link: http://www.iam.unibe.ch/~scg/Research/Traits/ Creo que esto tiene mucho que ver con la segunda forma de ver a las clases.
En particular, creo que la implementacion de traits en smalltalk, puede producir un cambio en la dimension del paradigma. Si pensamos a los paradigmas como un especio, smalltalk esta claramente sobre la recta de Object Based. El poder generar una jerarquia de clases en base a comportamiento compartido desplaza un poco el lenguaje sobre el espacio Object Based-funncional.
Bueno, espero opiniones, por que esto me interesa :) (y ademas espero que no haya sido un lio)

saludos
/jb

Calculando RSA con una regla

4 * 5
Esta historia se me puede hacer un poco larga...
Hace un tiempo, tenuki, me conto que se podia multiplicar usando rectas.
El asunto es asi (hice dibujitos):
Supongamos que quiero multiplicar 4 * 5 (señorita señorita, esa la se!). Lo primero que tenemos que hacer es trazar un segmento de tamaño 5, que empieza en el punto 1@0, y llega hasta el punto 1@5 (segmento azul). Despues trazamos una recta, que pasa por el punto 0@0, y por el punto 1@5 (recta roja). Finalmente, a baño maria, hacemos un segmento de tamaño 4, sobre el eje de las x, y que salga del punto 0@0 (segmento verde). Cuando unimos el segmento verde con la recta roja (a travez del segmento amarillo), encontramos dos puntos. El que esta sobre la recta roja (4@20) nos da en su coordenada y el valor que estabamos buscando (y si.. 20).
El metodo parece piola, y la explicacion es simple. Cuando trazamos la recta roja, en realidad, estamos trazando una recta que tiene como ecuacion y = 5 * x. Luego hacemos el segmento verde, y nos encontramos con que x = 4. Entonces y = 5 * 4 = 20. Lo mismo que queriamos, pero expresado obtenido de una manera distinta.
La pregunta seria que tiene esto que ver con RSA. Para eso necesitamos poder elevar un numero a una potencia (elevar es multiplicar un numero por si mismo n veces).
Eso tambien lo podemos hacer con regla, y para eso usamos mapas iterativos.

3 ^ 5
Receta para calcular 3 ^ 5 (4 porciones):
Primero le cambiamos el nombre a los ejes, el eje X ahora se llama x(t) y el eje Y ahora se llama x(t+1). Despues trazamos dos rectas que salen desde el origen 0@0. La primer recta es una recta de 45 grados (recta verde). La segunda recta (recta roja) es una recta de pendiente 3 (si estubiesemos elevando 17 ^ n, la recta seria de pendiente 17).
Partiendo del punto 1@1 (que esta en la recta verde), trazamos una recta paralela al eje x(t+1), hasta encontrar a la recta roja (esto va a ocurrir en el punto 1@3). Maravillosamente encontramos el valor de 3 ^ 1. Luego trazamos una recta horizontal hasta volver a encontrar la recta verde, y subimos hasta la recta roja, ahora el punto va a ser el 3@9 (jua, 3 ^ 2 !!). Si hacemos esto 3 veces mas, nos vamos a encontrar con el punto 81@243 (y 243 es 3^5 !!). Esto significa que aprendimos a potenciar con regla!
Yo me pregunte por que esto funcionaba... Bueno, parte de la respuesta es que cuando trazabamos la recta de pendiente 3, podiamos escribir esta ecuacion en diferencias:
x(t+1) = 3 * x(t) que es la ecuacion de la recta roja.
La solucion a esta ecuacion en diferencias es: x(t) = C(3^t) donde C es alguna constante que depende del punto inicial (en nuestro caso, y para el bien de todos nosotros, es 1).
Asi que esto funcionaba :)
Ahora nada mas nos falta solucionar un problema: el modulo
El ejemplo anterior lo vamos a cambiar un poco, vamos a calcular 3 ^ 6 % 251 (esto se lee: 3 a la sexta modulo 251). Al grafico anterior, trazamos una recta horizontal para x(t+1) = 251.
Ahora, cuando elevamos un numero, seguimos los mismos paso que antes, pero si, cuando estamos trazando la recta vertical, encontramos la recta azul en vez de la roja, cambiamos un poco el asunto. En ese caso, trazamos una recta paralela a la recta roja, y que pase por este punto (por el punto donde nos encontramos con la recta azul). Seguimos esta nueva recta hasta que x(t+1) = 0 (es decir, hasta el eje x(t)), y volvemos a seguir los pasos anteriores (con la precaucion de que podemos volver a encontrarnos la recta azul).
Asi vemos como en el ejemplo (3 ^ 6 % 251), nos encontramos dos veces con la recta azul antes de encontrar a la recta roja en el punto que estavamos buscando (227/3@227) donde 227 es el resultado que buscabamos.
3^6 % 251
Finalmente me queda pendiente algo interesante para pensar, que es la existencia de atractores. Un atractor en este caso seria algun conjunto de dibujos que sean iguales. Lo interesante es que una vez que se toca un punto de la recta roja mas de una vez, sabemos que estamos en presencia de un atractor, es decir, vamos a estar un dibujo que ya habiamos hecho. Hay un pequeño teorema que garantiza la presencia de atractores, si se cumplen algunas situaciones, pero el numero al que tenemos que elevar es un numero muy alto. Seria piola achicar ese numero.

Mas on Random (no-tan) Walk

Siguiendo con el processo del articulo pasado, un amigo, felipe, me sugirio que trate de probar el mismo proceso, variando los parametros, pero usando siempre los mismos numeros aleatorios. La idea de la prueba esta es: Si con un parametro igual a 2, en un determinado momento t habia un movimiento hacia arriba, con parametro igual a 3, en el mismo momento, tendria que existir el mismo movimiento (en la misma direccion), pero con el cambio propio que le da el processo.
El cambio fue simple (que grande squeak), y en el momento del cambio, me di cuenta de algo interesante. Yo estaba usando variables aleatorias con distribuccion Normal, y lo unico que hacia era cambiar el desvio de las variables x(t) dependiendo de la variable y(t-1). Aca hay que recordar que si queremos cambiar el desvio podemos hacer lo siguiente: x(t) = sigma * N(0;1) (lo contrario seria normalizar la variable N(0;1) = x(t) / sigma). Bueno, entonces para hacer el cambio, lo unico que tenia que hacer era generar n variables N(0;1), y despues, para cada uno de los parametros, aplicaba el sigma. Mmmm, creo que esto no quedo claro, bueno, esto queda asi:

x(t) = b * e ^ (a * sin(y(t-1))) * N(0;1)

y(t) = x(t) + y(t-1)

Y esto me parecio bastante piola. Basicamente una cosa que se puede leer de estas ecuaciones, es que el siguiente movimiento va a ser una variable aleatoria, multiplicada por una constante asociada a cada valor y. Esto seria como si nos imaginaramos que hay un numero para cada recta paralela al eje donde graficamos el tiempo (el horizontal), entonces, por ejemplo, a la paralela correspondiente a: a = 2; y = 3 el numero asociado a esa recta es k =2.0830282438600087, esto quiere decir que si en el momento t-1, y era igual a 3, la variable aleatoria que salga sorteada, va a ser multiplicada k.

Esto me triggereo otro pensamiento. En los modelos donde se tiene algo asi como Y = aX donde X es una variable Normal con una media y una varianza, al buscar el parametro a, en realidad lo unico que se esta haciendo es buscar cuales son los parametros (media y varianza) de la variable Y, ya que Y se podria simular como una normal y = N(a*E(X);a*sqrt(Var(X))) donde sqrt(x) es la raiz cuadrada de x. Bueno, es otra forma de pensar lo mismo...
saludos
/jb

ps: dejo un par de imagenes usando el comentario de felipe, donde todas las imagenes se hicieron con las mismas 10000 N(0;1), pero con distinto parametro a

ps2: quedan pendientes los colorsitos
a=1
a=2
a=3
a=4