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 < b < 2, no tenemos gran idea de lo que va a pasar, pero el folklore fractal nos permite sospechar que va a haber una ley de escala entre el contenido del intervalo y su longitud. Una dimension fractal, vamos. Venga, pues ejecutamos un rato el codigo y vemos que leyes de potencias salen. Fijaos sobre todo en las inversas del ajuste de minimos cuadrados, que son las dos ultimas columnas: [code lang="raroquenoesta"] 0.1 21 10.2215654543 10.2220691497 9.21554333643 | 0.098 0.098 0.1 21 10.6414413083 10.6340220642 9.66758183443 | 0.094 0.094 0.1 21 7.65785791454 7.64289916897 6.66616080022 | 0.131 0.131 0.1 21 9.03031553255 9.03397493797 8.41987390083 | 0.111 0.111 0.2 21 4.65155228853 4.65019960409 3.70803482651 | 0.215 0.215 0.2 21 4.7092066135 4.71368302913 3.66266507205 | 0.212 0.212 0.2 21 5.31074565522 5.28702270808 4.3601627827 | 0.188 0.189 0.2 21 5.45839569952 5.46993848673 4.59458429733 | 0.183 0.183 0.2 21 5.83462477135 5.8564820055 4.8274032778 | 0.171 0.171 0.4 21 1.83451060471 1.81096494583 0.717692574423 | 0.545 0.552 0.4 21 2.47018412976 2.45082565822 1.3613718474 | 0.405 0.408 0.4 21 2.49601035805 2.49387109342 1.52922470856 | 0.401 0.401 0.4 21 3.36861402177 3.436989876 2.54679249394 | 0.297 0.291 0.5 21 1.63770991123 1.61747116549 0.760356432681 | 0.611 0.618 0.5 21 1.68566575589 1.63785973424 0.670415437561 | 0.593 0.611 0.5 21 1.78332522902 1.71580018618 0.717280512823 | 0.561 0.583 0.5 21 2.40239982042 2.47533986413 1.42345780308 | 0.416 0.404 0.7 21 1.36275890665 1.31614771144 0.332481839006 | 0.734 0.760 0.7 21 1.48185928701 1.51426384241 0.596712802431 | 0.675 0.660 0.7 21 1.56630315806 1.60428744524 0.558623488121 | 0.638 0.623 0.7 21 1.69957421131 1.76821248138 0.734001583234 | 0.588 0.566 0.9 21 1.06908594435 0.933263576648 0.173359067467 | 0.935 1.072 0.9 21 1.17471966923 1.11284062068 0.165422583512 | 0.851 0.899 0.9 21 1.24097842931 1.28089286355 0.268913193027 | 0.806 0.781 0.9 21 1.24521989096 1.25095473827 0.204964077263 | 0.803 0.799 1.0 21 1.04630285531 1.00331140079 -2.15669570539e-16 | 0.956 0.997 1.0 21 1.07769987303 1.06240437302 -1.43779713692e-16 | 0.928 0.941 1.0 21 1.08120311841 0.976718941315 0.125182594747 | 0.925 1.024 1.0 21 1.09264515958 1.03451857085 0.174950463183 | 0.915 0.967 1.1 21 1.0058975402 0.750527707987 0.0961440263979 | 0.994 1.332 1.1 21 1.01660273885 0.726634095302 0.0739695911276 | 0.984 1.376 1.1 21 1.03432731085 0.98588293891 0.263265432537 | 0.967 1.014 1.1 21 1.07115363408 0.90503139179 0.0631928991542 | 0.934 1.105 1.3 21 0.991933368553 0.769953867525 -1.19816428077e-16 | 1.008 1.299 1.3 21 1.02891798659 0.816123578885 0.0208273031947 | 0.972 1.225 1.3 21 1.03174200349 0.815488117252 0.0786105411361 | 0.969 1.226 1.3 21 1.05128184958 0.808258753134 0.0294745071911 | 0.951 1.237 1.5 21 1.00515263544 0.579819941948 -5.99082140385e-17 | 0.995 1.725 1.5 21 1.01303768898 0.673988410337 0.0215068079652 | 0.987 1.484 1.5 21 1.02587547379 0.663660567458 -4.1935749827e-17 | 0.975 1.507 1.5 21 1.04824778494 0.733044916946 0.0562548668574 | 0.954 1.364 1.9 21 1.00111938437 0.461678724256 -2.39632856154e-17 | 0.999 2.166 1.9 21 1.02520517843 0.58669081982 -4.1935749827e-17 | 0.975 1.704 1.9 21 1.02522041532 0.52934586514 -5.99082140385e-18 | 0.975 1.889 1.9 21 1.03455532005 0.609178814618 0.0211572781037 | 0.967 1.642 2.0 21 0.996494113688 0.458529058988 -2.39632856154e-17 | 1.004 2.181 2.0 21 1.00973714499 0.49041126089 -8.98623210578e-18 | 0.990 2.039 2.0 21 1.01024230667 0.506934761729 -2.24655802644e-18 | 0.990 1.973 2.0 21 1.01202314635 0.48770648771 -1.49770535096e-18 | 0.988 2.050 2.1 21 1.00182622643 0.576169020827 -2.99541070193e-17 | 0.998 1.736 2.1 21 1.00284106918 0.512319911546 -5.99082140385e-17 | 0.997 1.952 2.1 21 1.00889558877 0.484739406135 -2.62098436419e-18 | 0.991 2.063 2.1 21 1.02196668605 0.511140047796 -3.74426337741e-18 | 0.979 1.956 2.2 21 0.97589734141 0.47555151319 -9.58531424616e-17 | 1.025 2.103 2.2 21 1.00735679224 0.488186293291 0.00226927797562 | 0.993 2.048 2.2 21 1.01059268424 0.44823257229 2.24655802644e-18 | 0.990 2.231 2.2 21 1.01218046491 0.456664291781 -3.74426337741e-18 | 0.988 2.190 2.5 21 1.00090235046 0.401806983522 7.48852675482e-19 | 0.999 2.489 2.5 21 1.00197124232 0.464038218364 -2.39632856154e-17 | 0.998 2.155 2.5 21 1.00912887228 0.468148804093 0.0186860940474 | 0.991 2.136 2.5 21 1.00916804822 0.370121556675 4.49311605289e-18 | 0.991 2.702 2.9 21 0.972394716299 0.255410402177 -1.79724642116e-17 | 1.028 3.915 2.9 21 0.989122681052 0.369322087019 -4.49311605289e-18 | 1.011 2.708 2.9 21 0.995254739229 0.358256016446 0.0014406745445 | 1.005 2.791 2.9 21 1.00074419185 0.367401325739 0.00687184912409 | 0.999 2.722 3.0 21 0.998601234573 0.322878557304 -1.12327901322e-18 | 1.001 3.097 3.0 21 1.00141126952 0.358920253703 1.19816428077e-17 | 0.999 2.786 3.0 21 1.00521760062 0.319382388595 0.0100505368038 | 0.995 3.131 3.0 21 1.01061046245 0.378741341224 8.42459259917e-19 | 0.990 2.640 3.1 21 0.983718843546 0.315996058698 2.62098436419e-18 | 1.017 3.165 3.1 21 0.997954016029 0.296572714283 9.36065844352e-19 | 1.002 3.372 3.1 21 1.00122589274 0.369748338667 0.00229312488497 | 0.999 2.705 3.1 21 1.00580877392 0.365142965802 0.00663213587297 | 0.994 2.739 4.0 21 1.0004116485 0.244518485921 5.99082140385e-18 | 1.000 4.090 4.0 21 1.00050900344 0.296338796516 0.00165401924686 | 0.999 3.375 4.0 21 1.00130111749 0.271827343336 2.99541070193e-17 | 0.999 3.679 4.0 21 1.00663887701 0.247636779564 0.00849941533798 | 0.993 4.038 5.0 21 0.9943337524 0.222995672828 2.39632856154e-17 | 1.006 4.484 5.0 21 0.995315611176 0.264445123538 0.00524689699634 | 1.005 3.782 5.0 21 1.00701066078 0.227780600217 6.58990354424e-17 | 0.993 4.390 5.0 21 1.01401145796 0.249559622575 0.00578083319398 | 0.986 4.007 8.0 21 0.98346612918 0.134672452357 2.99541070193e-17 | 1.017 7.425 8.0 21 0.984714702399 0.166741954129 2.99541070193e-17 | 1.016 5.997 8.0 21 0.987523204813 0.130386848118 5.99082140385e-18 | 1.013 7.669 8.0 21 1.01184575427 0.1833677755 4.79265712308e-17 | 0.988 5.454 12.0 21 0.99290709814 0.141661447191 0.0 | 1.007 7.059 12.0 21 0.998197475165 0.143316659926 1.19816428077e-17 | 1.002 6.978 12.0 21 0.999387185328 0.172183181435 9.58531424616e-17 | 1.001 5.808 12.0 21 1.00165466091 0.161164575256 3.59449284231e-17 | 0.998 6.205 20.0 21 1.00450195603 0.142938681948 0.00827642053648 | 0.996 6.996 20.0 21 1.00469595135 0.16088611642 0.00561075814881 | 0.995 6.216 20.0 21 1.00739431872 0.12650962307 8.38714996539e-17 | 0.993 7.905 20.0 21 1.01266206973 0.156867141824 0.0127964121666 | 0.987 6.375 [/code] En efecto, hay algo como una transicion de fase, el momento en que b es igual o mayor que 2 (los 1.0 de la tabla, con la parametrizacion que usa numpy para las distribuciones "de Pareto") y comienza a existir un valor medio finito para la distribucion. A partir de ahi el intervalo tendra tipicamente una longitud simplemente multiplo del numero de vuelos. Curiosamente no detectamos nada nuevo para b=3, donde de verdad se activa el teorema central del limite. Mirando con detalle, parece ser que tanto la longitud total como el mayor salto generado escalan con el numero de puntos con un exponente similar al inverso del que parametriza la distribucion, b-1. Cuando llegamos a b=2 nos encontramos con que la longitud deja de escalar; hemos "saturado" la dimension fractal, y ya no puede ser mayor que uno. En cambio el mayor salto sigue creciendo con el numero de puntos que generamos. Yo habia esperado que esto era lo que iba a anularse en b=3 (los 2.0 de la tabla), pero no ocurre asi. Da la impresion de que se frena un poco, pero no mucho mas.