[Esta es una pregunta de socio para calcular una probabilidad exactamente ]
Esta tarea consiste en escribir código para calcular una probabilidad de manera exacta y rápida . El resultado debe ser una probabilidad precisa escrita como una fracción en su forma más reducida. Es decir, nunca debería salir, 4/8
sino más bien 1/2
.
Para algún número entero positivo n
, considere una cadena uniformemente aleatoria de 1s y -1s de longitud n
y llámela A. Ahora concatene a A
su primer valor. Es decir, A[1] = A[n+1]
si la indexación desde 1. A
ahora tiene longitud n+1
. Ahora también considere una segunda cadena aleatoria de longitud n
cuyos primeros n
valores son -1, 0 o 1 con probabilidad 1 / 4,1 / 2, 1/4 cada uno y llámelo B.
Consideremos ahora el producto interno de A[1,...,n]
y B
y el producto interno de A[2,...,n+1]
y B
.
Por ejemplo, considere n=3
. Los posibles valores de A
y B
podrían ser A = [-1,1,1,-1]
y B=[0,1,-1]
. En este caso, los dos productos internos son 0
y 2
.
Su código debe generar la probabilidad de que ambos productos internos sean cero.
Copiando la tabla producida por Martin Büttner tenemos los siguientes resultados de muestra.
n P(n)
1 1/2
2 3/8
3 7/32
4 89/512
5 269/2048
6 903/8192
7 3035/32768
8 169801/2097152
Idiomas y bibliotecas
Puede utilizar cualquier idioma y bibliotecas disponibles gratuitamente que desee. Debo poder ejecutar su código, así que incluya una explicación completa sobre cómo ejecutar / compilar su código en Linux si es posible.
La tarea
Su código debe comenzar con n=1
y dar la salida correcta para cada n creciente en una línea separada. Debería detenerse después de 10 segundos.
El marcador
El puntaje es simplemente el más alto n
alcanzado antes de que su código se detenga después de 10 segundos cuando se ejecuta en mi computadora. Si hay un empate, el ganador será el que obtenga el puntaje más alto más rápido.
Tabla de entradas
n = 64
en Python . Versión 1 por Mitch Schwartzn = 106
en Python . Versión del 11 de junio de 2015 por Mitch Schwartzn = 151
en C ++ . La respuesta del puerto de Mitch Schwartz por kirbyfan64sosn = 165
en Python . Versión del 11 de junio de 2015, la versión de "poda" de Mitch Schwartz conN_MAX = 165
.n = 945
en Python por Min_25 usando una fórmula exacta. ¡Asombroso!n = 1228
en Python por Mitch Schwartz usando otra fórmula exacta (basada en la respuesta anterior de Min_25).n = 2761
en Python por Mitch Schwartz usando una implementación más rápida de la misma fórmula exacta.n = 3250
en Python usando Pypy por Mitch Schwartz usando la misma implementación. Esta puntuación debepypy MitchSchwartz-faster.py |tail
evitar la sobrecarga de desplazamiento de la consola.
fuente
Respuestas:
Pitón
Una fórmula de forma cerrada de
p(n)
esUna función generadora exponencial de
p(n)
esdonde
I_0(x)
está la función Bessel modificada del primer tipo.Editar en 2015-06-11:
- actualizado el código de Python.
Editar en 2015-06-13:
- se agregó una prueba de la fórmula anterior.
- arreglado el
time_limit
.- Se agregó un código PARI / GP.
Pitón
PARI / GP
Prueba:
este problema es similar a un problema de caminata aleatoria bidimensional (restringida).
Si
A[i] = A[i+1]
, podemos pasar de(x, y)
a(x+1, y+1)
[1 manera],(x, y)
[2] o formas(x-1, y-1)
[1 manera].Si
A[i] != A[i+1]
, podemos pasar de(x, y)
a(x-1, y+1)
[1 manera],(x, y)
[2] o formas(x+1, y-1)
[1 manera].Deje
a(n, m) = [x^m]((x+1)^n + (x-1)^n)
,b(n) = [x^n](1+x)^{2n}
yc(n)
sea la cantidad de formas de moverse de(0, 0)
a(0, 0)
conn
pasos.Entonces,
c(n) = \sum_{i=0}^n a(n, i) * b(i) * b(n-i).
Desde entonces
p(n) = c(n) / 8^n
, podemos obtener la fórmula de forma cerrada anterior.fuente
Pitón
Nota: ¡ Felicitaciones a Min_25 por encontrar una solución de forma cerrada!
Gracias por el interesante problema! Se puede resolver con DP, aunque actualmente no me siento muy motivado para optimizar la velocidad para obtener una puntuación más alta. Podría ser bueno para el golf.
El código llegó
N=39
en 10 segundos en esta vieja computadora portátil con Python 2.7.5.Para tupla
(a,b,s,t)
:a
es el primer elemento deA
,b
es el último elemento deB
,s
es el producto interno deA[:-1]
yB
, yt
es el producto interno deA[1:-1]
yB[:-1]
, utilizando la notación de corte de Python. Mi código no almacena las matricesA
ni enB
ninguna parte, por lo que uso esas letras para referirme a los siguientes elementos que se agregaránA
yB
, respectivamente. Esta elección de nombres variables hace que la explicación sea un poco incómoda, pero permite una apariencia agradableA*b+a*B
en el código mismo. Tenga en cuenta que el elemento que se agregaA
es el penúltimo, ya que el último elemento siempre es el mismo que el primero. He usado el truco de Martin Büttner de incluir0
dos veces enB
candidatos para obtener la distribución de probabilidad adecuada. El diccionarioX
(que lleva el nombreY
deN+1
) realiza un seguimiento del recuento de todos los arreglos posibles de acuerdo con el valor de la tupla. Las variablesn
yd
representan numerador y denominador, es por eso que cambié el nombren
de la declaración del problema comoN
.La parte clave de la lógica es que se puede actualizar a partir
N
deN+1
utilizar sólo los valores de la tupla. Los dos productos internos especificados en la pregunta están dados pors+A*B
yt+A*b+a*B
. Esto está claro si examina un poco las definiciones; tenga en cuenta que[A,a]
y[b,B]
son los dos últimos elementos de las matricesA
yB
respectivamente.Tenga en cuenta que
s
yt
son pequeños y acotados de acuerdo conN
, y para una implementación rápida en un lenguaje rápido, podríamos evitar los diccionarios a favor de las matrices.Es posible aprovechar la simetría considerando valores que difieren solo en el signo; No he investigado eso.
Observación 1 : El tamaño del diccionario crece cuadráticamente
N
, donde tamaño significa número de pares clave-valor.Observación 2 : Si establecemos un límite superior
N
, podemos podar tuplas para las cualesN_MAX - N <= |s|
y de manera similar parat
. Esto podría hacerse especificando un estado absorbente, o implícitamente con una variable para mantener el recuento de estados podados (que tendrían que multiplicarse por 8 en cada iteración).Actualización : esta versión es más rápida:
Optimizaciones implementadas:
main()
- el acceso variable local es más rápido que el globalN=1
explícitamente para evitar la verificación(1,-1) if N else [a]
(lo que obliga a que el primer elemento de la tupla sea consistente al agregar elementos alA
comenzar desde la lista vacía)c
para agregar un0
a enB
lugar de hacer esas operaciones dos veces8^N
por lo que no necesitamos hacer un seguimientoA
as1
y dividir el denominador entre2
, ya que los pares válidos(A,B)
conA[1]=1
y aquellos conA[1]=-1
pueden ser puestos en correspondencia uno a uno mediante la negaciónA
. Del mismo modo, podemos arreglar el primer elemento deB
como no negativo.N_MAX
para ver qué puntaje puede obtener en su máquina. Podría reescribirse para encontrar un apropiadoN_MAX
automáticamente mediante búsqueda binaria, pero parece innecesario? Nota: No es necesario que verifiquemos la condición de poda hasta alcanzar el límiteN_MAX / 2
, por lo que podríamos acelerar un poco iterando en dos fases, pero decidí no hacerlo por simplicidad y limpieza del código.fuente
N=57
la primera versión yN=75
la segunda.Pitón
Usando la idea de Min_25 de caminata aleatoria, pude llegar a una fórmula diferente:
Aquí hay una implementación de Python basada en Min_25:
Explicación / prueba:
Primero resolvemos un problema de conteo relacionado, donde permitimos
A[n+1] = -A[1]
; es decir, el elemento adicional concatenadoA
puede ser1
o-1
independientemente del primer elemento. Por lo tanto, no necesitamos hacer un seguimiento de cuántas vecesA[i] = A[i+1]
ocurre. Tenemos la siguiente caminata aleatoria:Desde
(x,y)
podemos pasar a(x+1,y+1)
[1 vía],(x+1,y-1)
[1 vía],(x-1,y+1)
[1 vía],(x-1,y-1)
[1 vía],(x,y)
[4 vías]donde
x
yy
soporte para los dos productos escalares, y estamos contando el número de maneras de pasar de(0,0)
que(0,0)
enn
los pasos. Ese recuento se multiplicaría por2
para tener en cuenta el hecho de queA
puede comenzar con1
o-1
.Nos referimos a permanecer en
(x,y)
un movimiento cero .Repetimos el número de movimientos
i
que no son cero , lo que tiene que ser uniforme para volver a hacerlo(0,0)
. Los movimientos horizontales y verticales forman dos caminatas aleatorias unidimensionales independientes, que pueden contarseC(i,i/2)^2
, dondeC(n,k)
es el coeficiente binomial. (Para una caminata conk
pasos hacia la izquierda y hacia lak
derecha, hayC(2k,k)
formas de elegir el orden de los pasos). Además, hayC(n,i)
formas de colocar los movimientos y4^(n-i)
formas de elegir los movimientos cero. Entonces obtenemos:Ahora, tenemos que volver al problema original. Defina un par permitido
(A,B)
para que sea convertible siB
contiene un cero. Defina un par(A,B)
para que sea casi permisible siA[n+1] = -A[1]
los dos productos de punto son cero.Lema: Para un hecho
n
, los pares casi permitidos están en correspondencia uno a uno con los pares convertibles.Podemos (reversiblemente) convertir un par convertible
(A,B)
en un par casi permisible(A',B')
al negarA[m+1:]
yB[m+1:]
dóndem
está el índice del último ceroB
. La comprobación de esto es sencilla: si el último elemento deB
es cero, no necesitamos hacer nada. De lo contrario, cuando negamos el último elemento deA
, podemos negar el último elemento deB
para preservar el último término del producto de puntos desplazados. Pero esto niega el último valor del producto de punto no desplazado, por lo que solucionamos esto al negar el penúltimo elemento deA
. Pero luego esto arroja el penúltimo valor del producto desplazado, por lo que negamos el penúltimo elemento deB
. Y así sucesivamente, hasta llegar a un elemento cero enB
.Ahora, solo tenemos que mostrar que no hay pares casi permitidos para los que
B
no contenga un cero. Para un producto de punto a igual a cero, tenemos que tener un número igual de1
y-1
términos a cancelarse. Cada-1
término se compone de(1,-1)
o(-1,1)
. Por lo tanto, la paridad del número de-1
que se produce se fija de acuerdo conn
. Si el primer y el último elemento deA
tienen signos diferentes, cambiamos la paridad, por lo que esto es imposible.Entonces obtenemos
que da la fórmula anterior (reindexando con
i' = i/2
).Actualización: Aquí hay una versión más rápida que usa la misma fórmula:
Optimizaciones implementadas:
p(n)
C(n,k)
conk <= n/2
fuente
p(n)
no es necesario que sea una función por partes. En general, sif(n) == {g(n) : n is odd; h(n) : n is even}
puede escribirf(n) == (n-2*floor(n/2))*g(n) + ((n+1)-2*(floor((n+1)/2)))*h(n)
o usar enn mod 2
lugar de(n-2*floor(n/2))
. Ver aquíExplicación de la fórmula de Min_25
Min_25 publicó una gran prueba, pero tardó un tiempo en seguirla. Esta es una pequeña explicación para completar entre las líneas.
a (n, m) representa el número de formas de elegir A de manera que A [i] = A [i + 1] m veces. La fórmula para a (n, m) es equivalente a a (n, m) = {2 * (n elija m) para nm par; 0 para nm impar.} Solo se permite una paridad porque A [i]! = A [i + 1] debe ocurrir un número par de veces para que A [0] = A [n]. El factor 2 se debe a la elección inicial A [0] = 1 o A [0] = -1.
Una vez que el número de (A [i]! = A [i + 1]) se fija en q (denominado i en la fórmula c (n)), se separa en dos recorridos aleatorios 1D de longitud q y nq. b (m) es la cantidad de formas de realizar una caminata aleatoria unidimensional de m pasos que termina en el mismo lugar donde comenzó, y tiene un 25% de posibilidades de moverse hacia la izquierda, un 50% de posibilidades de quedarse quieto y un 25% de posibilidades de moviéndose a la derecha Una forma más obvia de establecer la función generadora es [x ^ m] (1 + 2x + x ^ 2) ^ n, donde 1, 2x y x ^ 2 representan izquierda, sin movimiento y derecha, respectivamente. Pero entonces 1 + 2x + x ^ 2 = (x + 1) ^ 2.
fuente
C ++
Solo un puerto de la (excelente) respuesta de Python de Mitch Schwartz. La principal diferencia es que solía
2
representar-1
para laa
variable e hice algo similarb
, lo que me permitió usar una matriz. Usando Intel C ++ con-O3
, lo tengoN=141
! Mi primera versión tieneN=140
.Esto usa Boost. Probé una versión paralela pero me encontré con algunos problemas.
fuente
g++ -O3 kirbyfan64sos.cpp -o kirbyfan64sos -lboost_system -lboost_timer -lboost_chrono -lrt -lgmp
compilarse. (Gracias a aditsu.)