#1306 closed task (fixed)
Overdispesed pseudo-uniform sampler under constraints
| Reported by: | Jorge | Owned by: | Víctor de Buen Remiro |
|---|---|---|---|
| Priority: | high | Milestone: | Mantainance |
| Component: | Math | Version: | 2.0.1 |
| Severity: | critical | Keywords: | |
| Cc: |
Description
In order to initialize multiple MCMC is it needed a functionality which provide the sampling of multiple points dispersed over the constrained domain.
Attached to this ticket there is a paper which describe a method to generate a sample with the previous characteristics.
Attachments (1)
Change History (17)
Changed 14 years ago by
| Attachment: | 10.1.1.153.7868.pdf added |
|---|
comment:1 Changed 14 years ago by
| Status: | new → accepted |
|---|
comment:2 Changed 14 years ago by
Estoy pensando que tanto el método del papel como los míos suponen matar moscas a cañonazos. Lo que se pretende es encontrar unos cuantos puntos factibles más o menos dispersos, para lo cual no es necesario generar una uniforme perfecta ni nada por el estilo. Basta con un paseo arbitrario que se mueva en direcciones aleatorias y a grandes pasos dentro del politopo. Para ello el radio de dispersión se multiplicará por cierto factor mayor que 1 cada vez que el punto generado esté dentro de la región factible y se dividirá por ese mismo factor cada vez que se salga de la misma.
En la siguiente figura se puede ver una cadena bidimensional de unos 500 puntos generados en centésimas de segundo con el casi trivial código del ejemplo FindArbitraryPointInPolytope.tol
Puede verse como el algoritmo tiende a cubrir toda la región factible más o menos por igual, aunque con cierta tendencia a acumular puntos en los vértices, lo cual no tiene la menor importancia para lo que se pretende hacer en este caso.
En medio segundo genera mil muestras con 100 variables y 150 restricciones y en 37 segundos para 1000 variables y 1500 restricciones.
comment:3 Changed 14 years ago by
Una vez tenemos una muestra dispersa como la del ejemplo se puede obtener una selección de unos pocos puntos lejanos entre sí siguiendo estos pasos:
- Se obtiene el primer punto como el centro de masas, o sea, la media aritmética de todos los puntos generados.
- Seleccionamos como segundo punto al más alejado del centro de masas.
- Seleccionamos como tercer punto el que dé máxima suma de distancias a los dos anteriores .1 Seleccionamos como k-ésimo punto el que dé máxima suma de distancias a los h anteriores, siendo h=Min(k-1,n) y n la dimensión del espacio
comment:6 Changed 14 years ago by
| Summary: | Uniform sampler under constraints → Overdispesed pseudo-uniform sampler under constraints |
|---|
Una forma más sencilla de obtener puntos sobre-dispersos sería tomar el centro de masas de la muestra inicial y generar direcciones aleatorias para formar rectas radiales en las que buscar puntos alejados del centro dentro de esa recta.
comment:7 Changed 14 years ago by
| Resolution: | → fixed |
|---|---|
| Status: | accepted → closed |
comment:11 Changed 14 years ago by
Se ha creado una función general para generar cuasi-uniformes arbitrarias
VMatrix BysSampler::Mcmc.ConstrainedCuasiUniform ( Code g, VMatrix x0, Real size, Real burnin, Real thinning);
Starting at a given feasible point, draws a Montecarlo Markov Chain from a cuasi-uniform distributed random variable under arbitrary constraining inequations g(x)<=0, where
- g is a funtion declared as VMatrix g(VMatrix x)
- x0 is the initial feasible point
- size is the final length of the MCMC
- burnin is the number of skiped initial points
- thinning id the length of skiped internal intervals
Y también para el casi de politopos
VMatrix Mcmc.ConstrainedCuasiUniform.Polytope ( VMatrix A, VMatrix a, VMatrix x0, Real size, Real burnin, Real thinning)
Starting at a given feasible point x0, draws a Montecarlo Markov Chain from a cuasi-uniform distributed random variable inside a polytope A*x<=a where
- A is the matrix of coefficients of the polytope
- a is the matrix of borders of the polytope
- x0 is the initial feasible point
- size is the final length of the MCMC
- burnin is the number of skiped initial points
- thinning id the length of skiped internal intervals
Por ejemplo, si queremos generar 5 puntos que estén dispersos podemos obviar los 100 primeros y luego tomar una de cada 50 iteraciones con esta sentencia:
VMatrix chain = BysSampler::Mcmc.ConstrainedCuasiUniform.Polytope( A, a, x0, 5, 100, 50);
comment:14 Changed 14 years ago by
La combinación de dos nuevas funciones produce muestras sobre-dispersas de forma mucho más eficiente que un simple thinning.
- Muestrea un paseo aleatorio en una región
g(x)<= 0con tendencia hacia la frontera. En cada paso se genera una dirección aleatoria y se busca en ambos sentidos un número máximomaxTryde veces hasta toparse con la frontera mediante el método de la bisectriz ponderada. El tamaño de paso se ajusta en cada intento según el valor del argumentofactor, multiplicando si no se ha alcanzado la frontera y dividiendo si se sobrepasa. Luego se genera el siguiente punto de forma aleatoria en el segmento pero con mayor probabilidad en los extremos cuanto más cercano a cero sea el argumentodispersionExponent. El casodispersionExponent=1correspondería de forma bastante aproximada a generar una uniforme en la región factible.
VMatrix chain = { Mat2VMat(BysSampler::CppTools::DispersedRandomWalk ( Code g, VMatrix x0, Real size = 500, Real burnin=0, Real thinning=1, Real factor=2.0, Real maxTry=5, Real dispersionExponent=0.02)) }; - Extrae de una colección de puntos una selección de ellos que estén lo más
alejados entre sí, empezando por el centro de masas.
VMatrix overDispersed = { BysSampler::ExtractOverDispersed(chain, Real size=5) };
El resultado se puede ver en este ejemplo de dos dimensiones de forma muy intuitiva
El método propuesto en 10.1.1.153.7868.pdf requiere resolver dos sistemas de ecuaciones en cada iteración (2.1.1 pág 5)
Hace tiempo que hice el diseño de un par de métodos que guardan cierto parecido con el propuesto para la generación de candidatos factibles en un politopo que sirven para cualquier distribución de probabilidad, lo cual incluye la uniforme, que es lo que se solicita aquí. Estos métodos son mucho más sencillos de implementar y resultan más prometedores y generales por lo que creo, que sería mejor empezar por alguno de ellos. Concretamente el Paseo aleatorio radial asimétrico sirve incluso para restricciones de desigualdad arbitrarias, no necesariamente lineales.