En primer lugar, no estoy seguro de lo que está pidiendo (aunque creo que lo sé) cuando quiere que su cadena “converja al estado de peso máximo”. ¿Quiere decir que desea que su muestra termine en el estado (un estado del espacio de estado completo) con la mayor probabilidad? Porque de hecho tu muestra lo hará. Incluso una muestra de fuerza bruta lo hará. Y ese estado es a menudo el estado menos interesante porque puedes obtenerlo de manera muy simple mediante otros experimentos numéricos. Una de mis áreas de trabajo es la termodinámica estocástica teórica, donde estos estados altamente probables corresponden a compuestos moleculares estables, y lo que nos interesa son los intermedios de reacción (que determinan las velocidades de reacción) que son mucho menos estables y, por lo tanto, menos probables.
Pero sospecho que lo que quiere decir es que desea que su cadena alcance la distribución de probabilidad entre estados con máxima probabilidad / convergencia / estabilidad de tiempo (para un sistema ergódico, esto será lo mismo). Sí, tiene razón en que hay una energía libre que esta distribución minimiza:
[matemáticas] F (\ {p_i \}) = \ sum_i p_i \ ln \ left (\ frac {p_i} {p_i ^ e} \ right) [/ math]
- ¿Qué describe un lagrangiano y por qué es una función de posición y velocidad?
- Explicaciones de Layman: ¿Qué es la amortiguación de Landau?
- Una barra medidora que pesa 0.6N está soportada por un punto de apoyo en la marca de 45 cm. Para mantener el equilibrio, se debe colocar un peso de 0,5 N a: A.45 cm o B.6 cm.
- ¿Cuál es la ley de viscosidad de Newton?
- ¿Es [matemáticas] 1 cm ^ 2 [/ matemáticas] igual a [matemáticas] 10 ^ {- 2} m ^ 2 [/ matemáticas] o [matemáticas] (10 ^ {- 2}) ^ 2 m ^ 2 [/ matemáticas] ?
donde [math] p_i ^ e [/ math] son las probabilidades de equilibrio. Esta energía libre es estrictamente no negativa y es cero siempre que [math] p_i = p_i ^ e [/ math] para todos [math] i [/ math].
Pero eso presupone que conoces tus probabilidades de equilibrio, ¡lo cual frustra el propósito! No obstante, hay una derivada dependiente del tiempo de la energía libre que no:
[matemáticas] \ frac {dF} {dt} = – \ sum_ {i \ neq j} p_i k_ {ij} \ ln (p_i k_ {ij}) = – \ sum_ {i> j} (p_i k_ {ij} – p_j k_ {ji}) \ ln \ left (\ frac {p_i k_ {ij}} {p_j k_ {ji}} \ right) [/ math]
donde [math] k_ {ij} [/ math] son las tasas de transición del estado [math] i [/ math] al estado [math] j [/ math]. Sin embargo, una vez más, si tiene las tasas de transición explícitas, es trivial recuperar la distribución estacionaria. La virtud de definir la energía libre y su tasa es, de hecho, justificar el método utilizado en el muestreo de Monte Carlo: resulta que un sistema con muestreo equilibrado detallado (que no requiere tasas de transición explícitas) tiene la misma derivada de tiempo de energía libre como lo anterior, y por lo tanto con el tiempo se establece en su energía libre mínima, que se logra mediante la distribución estacionaria.
Entonces, si está muestreando Gibbs, ya está tratando de minimizar la energía libre de la distribución muestreada. Para acelerar la convergencia, es posible que desee probar una técnica llamada muestreo general , que esencialmente introduce un sesgo para obligar al sistema a explorar su espacio de fase menos probable y luego explica el sesgo.
Introduce este sesgo en la sección de muestreo de Metrópolis de su código. Si su muestra de Gibbs hace lo que creo que hace, hay una parte en la que acepta aleatoriamente un nuevo estado candidato con probabilidad [matemática] \ alfa [/ matemática] que está determinada por su distribución de muestreo. Si [math] \ alpha \ geq 1 [/ math] siempre se acepta el estado; si [math] \ alpha <1 [/ math] la muestra obtiene un número aleatorio y acepta solo si el número aleatorio es mayor que [math] \ alpha [/ math]. Lo que hace el muestreo general es reemplazar [math] \ alpha [/ math] por [math] \ alpha / w (x_i) [/ math], haciendo que cada estado [math] x_i [/ math] sea más probable por un factor de [matemáticas] w (x_i) [/ matemáticas].
(Como un ejemplo simple: suponga que está ejecutando una muestra de Gibbs sobre una distribución conjunta de (x, y), y encuentre que los estados en los que y> 1 son aproximadamente cien veces menos comunes que los estados sin. Establecería w = 1 si y 1, y el área submuestreada ahora debe estar sesgada para que se muestree de manera equitativa).
Después de obtener una distribución ponderada [math] p_w (x_i) [/ math], debe eliminar los efectos de los pesos, lo que hace simplemente des-ponderando:
[matemáticas] p (x_i) = \ frac {p_w (x_i) / w (x_i)} {\ sum_i p_w (x_i) / w (x_i)} [/ matemáticas]
¡Pan comido!
En la práctica, esto converge mejor cuando [matemáticas] w (x_i) = 1 / p ^ e (x_i) [/ matemáticas] (su sesgo hace que la variable aleatoria sea completamente uniforme). Por supuesto, no puede obtener la distribución estacionaria exacta (o de lo contrario se haría su cálculo), pero puede hacer una carrera inicial no ponderada de fuerza bruta y tomar los pesos iniciales desde allí, o usar su propio (pre) juicio e imponer una distribución de probabilidad inicial, como un gaussiano. Luego, puede ejecutar sucesivamente iteraciones cortas para obtener mejores pesos, y después de unas pocas ejecuciones puede hacer una larga “producción” con un conjunto de pesos bastante buenos para obtener resultados precisos.
Para darle una idea de qué tan bien funciona esto, rutinariamente muestro sistemas químicos en los que los estados menos probables son aproximadamente 20-30 [matemáticas] k_B T [/ matemáticas] más enérgicos que los estados más probables, o aproximadamente [matemáticas] 10 ^ 9 [/ matemáticas] menos probable. Por lo general, me toma entre dos y tres carreras de aproximadamente [matemáticas] 10 ^ 7 [/ matemáticas] para tener una idea de la distribución inicial, seguido de una producción de aproximadamente [matemáticas] 10 ^ 8 [/ matemáticas] a [ matemática] 10 ^ 9 [/ matemática] pasos en los que los estados generalmente se submuestrean como máximo a la mitad. Las distribuciones resultantes normalmente tienen una incertidumbre en la probabilidad logarítmica de aproximadamente el 5%.
Para una convergencia aún mejor, puede limitar sus simulaciones a ventanas en las que la probabilidad de muestreo no varía demasiado, y luego combinar el muestreo general de las ventanas vecinas usando el método de análisis de histograma ponderado, o WHAM para abreviar.