vuelos subexponenciales y fractales

Una forma muy comoda de generar distribuciones de Poisson es aprovechar que la distribucion de los intervalos, de los tiempos de espera, es la exponencial. Vamos sacando aleatorios t1, t2, t3… de esa distribucion y simplemente ponemos «al vuelo» los puntos t1, t1+t2, t1+t2+t3,… El numero de puntos que generados de esa manera entren en un intervalo dado va a seguir la distribucion de Poisson.

En cierto modo la exponencial es algo forzoso. En Poisson, la decision de si un punto va a entrar en la cuenta es independiente de que hayan entrado otros puntos. Eso significa que la distribucion de intervalos no puede tener memoria: si nos dicen que han entrado los puntos x1 y x3, cualquier punto intermedio deberia tener la misma probabilidad. Pero desde el punto de vista de vuelos, la probabilidad de que hubiera un punto x2 es el producto de probabilidades del vuelo de x1 a x2 y del vuelo de x2 a x3. No queda mas narices que sea \(e^{-k(x3-x2)} e^{-k(x2-x1)}\), cualquier otra funcion va a privilegiar unos puntos x2 respecto a otros.

Dicho esto, me ha entrado la curiosidad de tipo de nubes de puntos salen cuando uno usa vuelos con distribuciones de estas que llamamos «invariantes de escala», las del tipo \(1/x^b\)
de mayor alcance que los exponenciales, Quizas tendria que llamarlas «supraexponenciales», dado que lo que ocurre es que hay un «heavy tail» y una probabilidad bastante alta de un vuelo mas largo que en el caso exponencial. Por otro lado, uno normalmente piensa de estos polinomios como una aproximacion progresiva a la caida exponencial, asi que en la cabeza las estoy llamando «sub» en vez de «supra»

Al grano. Lo que queria mirar en este post es si en algun momento la distribucion de puntos es un fractal. Una forma de verlos es comprobar como escala el contenido respecto al tamaño. Esto es, contamos por un lado el numero de puntos que hemos metido para generar un intervalo y por otro lado su longitud, la suma total. Esto lo podemos hacer con un codigo python

#!/usr/bin/python
import sys
import numpy as np
exp=float(sys.argv[1])
n=0
maximo=0
maxpromedio=0
suma=0
print
print
x=[]
y=[]
z=[]
m=[]
while True:
  n+=1
  intento=np.random.pareto(exp)
  suma+=intento
  if suma/n > maxpromedio:
     maxpromedio=suma/n
  if intento > maximo:
     maximo=intento
  if np.log2(n).is_integer() and np.log2(n)>4:
     x.append(np.log2(n))
     y.append( np.log2(suma))
     z.append(np.log2(maximo))
     m.append(np.log2(maxpromedio))
     p,V=np.polyfit(x,y,1)
     pz,V=np.polyfit(x,z,1)
     pm,V=np.polyfit(x,m,1)
     print exp, len(x), p, pz, pm,
     print  "| %1.3f %1.3f " % (1/p, 1/pz)
     #print exp, n, maximo, np.log2(n), np.log2(maximo),
     # np.log2(suma), np.log2(maxpromedio)

Y en realidad lo que estamos haciendo es para n puntos, sumar n variables independientes todas con la distribucion de Pareto correspondiente. Esperamos que entre en juego el teorema central del limite cuando las distribuciones tengan b > 3 y que ya antes, con b>2, el hecho de que exista valor medio de la distribucion haya hecho entrar en accion a la desigualdad de Markov.

Cuando 1