math - number - Muestra uniformemente al azar de una unidad n-dimensional simplex

random mathematica

Question

El muestreo uniforme al azar de una unidad n-dimensional simple es la forma elegante de decir que desea n números aleatorios de tal manera que

• todos ellos son no negativos,
• suman a uno, y
• todos los vectores posibles de n números no negativos que suman uno son igualmente probables.

En el caso n = 2 desea muestrear uniformemente el segmento de la línea x + y = 1 (es decir, y = 1-x) que se encuentra en el cuadrante positivo. En el caso n = 3, estás muestreando desde la parte en forma de triángulo del plano x + y + z = 1 que está en el octante positivo de R3:

(Imagen de http://en.wikipedia.org/wiki/Simplex .)

Tenga en cuenta que seleccionar n números aleatorios uniformes y luego normalizarlos para que sumen uno no funciona. Terminas con un sesgo hacia números menos extremos.

De manera similar, escoger n-1 números aleatorios uniformes y luego tomar la nth como uno menos la suma de ellos también introduce un sesgo.

Wikipedia proporciona dos algoritmos para hacer esto correctamente: http://en.wikipedia.org/wiki/Simplex#Random_sampling (Aunque el segundo en la actualidad dice que solo es correcto en la práctica, no en teoría. Espero limpiar eso o aclárelo cuando lo entiendo mejor. Inicialmente me puse una "ADVERTENCIA: este tipo de papel afirma que lo siguiente es incorrecto" en esa página de Wikipedia y alguien más lo convirtió en la advertencia "funciona solo en la práctica".

Finalmente, la pregunta: ¿Cuál considera que es la mejor implementación del muestreo simplex en Mathematica (preferiblemente con la confirmación empírica de que es correcta)?

Preguntas relacionadas

• Generando una distribución de probabilidad.
• porcentajes aleatorios de java

Answer 1

Aquí hay una buena implementación concisa del segundo algoritmo de http://en.wikipedia.org/wiki/Simplex#Random_sampling :

SimplexSample[n_] := Rest@# - Most@# &[Sort@Join[{0,1}, RandomReal[{0,1}, n-1]]]

Se adaptó de aquí: http://www.mofeel.net/1164-comp-soft-sys-math-mathematica/14968.aspx (Originalmente tenía Union en lugar de Sort @ Join; esta última es un poco más rápida).

(Ver comentarios para alguna evidencia de que esto es correcto!)

Answer 2

Después de investigar un poco, encontré esta página que ofrece una buena implementación de la Distribución Dirichlet. Desde allí parece que sería bastante simple seguir el método 1 de Wikipedia. Esta parece ser la mejor manera de hacerlo.

Como prueba:

In[14]:= RandomReal[DirichletDistribution[{1,1}],WorkingPrecision->25]
Out[14]= {0.8428995243540368880268079,0.1571004756459631119731921}
In[15]:= Total[%]
Out[15]= 1.000000000000000000000000

Una parcela de 100 muestras:

texto alt http://www.public.iastate.edu/~zdavkeos/simplex-sample.png

Answer 3

Este código puede funcionar:

samples[n_] := Differences[Join[{0}, Sort[RandomReal[Range[0, 1], n - 1]], {1}]]

Básicamente, simplemente elige n - 1 lugares en el intervalo [0,1] para dividirlo, luego toma el tamaño de cada una de las piezas usando Differences .

Una ejecución rápida de Timing muestra que es un poco más rápido que la primera respuesta de Janus.

Answer 4

Estoy con zdav: la distribución de Dirichlet parece ser la forma más fácil de avanzar, y el algoritmo para muestrear la distribución de Dirichlet al que se refiere zdav también se presenta en la página de Wikipedia sobre la distribución de Dirichlet .

En cuanto a la implementación, es un poco costoso hacer primero la distribución completa de Dirichlet, ya que todo lo que realmente necesita son n muestras Gamma[1,1] aleatorias. Comparar a continuación
Implementación simple

SimplexSample[n_, opts:OptionsPattern[RandomReal]] :=
  (#/Total[#])& @ RandomReal[GammaDistribution[1,1],n,opts]

Implementación completa de Dirichlet

DirichletDistribution/:Random`DistributionVector[
 DirichletDistribution[alpha_?(VectorQ[#,Positive]&)],n_Integer,prec_?Positive]:=
    Block[{gammas}, gammas = 
        Map[RandomReal[GammaDistribution[#,1],n,WorkingPrecision->prec]&,alpha];
      Transpose[gammas]/Total[gammas]]

SimplexSample2[n_, opts:OptionsPattern[RandomReal]] := 
  (#/Total[#])& @ RandomReal[DirichletDistribution[ConstantArray[1,{n}]],opts]

Sincronización

Timing[Table[SimplexSample[10,WorkingPrecision-> 20],{10000}];]
Timing[Table[SimplexSample2[10,WorkingPrecision-> 20],{10000}];]
Out[159]= {1.30249,Null}
Out[160]= {3.52216,Null}

Así que el Dirichlet completo es un factor de 3 más lento. Si necesita m> 1 puntos de muestreo a la vez, probablemente podría ganar más haciendo (#/Total[#]&)/@RandomReal[GammaDistribution[1,1],{m,n}] .

Answer 5

He creado un algoritmo para la generación aleatoria uniforme sobre un símplex. Puede encontrar los detalles en el documento en el siguiente enlace: http://www.tandfonline.com/doi/abs/10.1080/03610918.2010.551012#.U5q7inJdVNY

En pocas palabras, puede usar las siguientes fórmulas de recursión para encontrar los puntos aleatorios sobre el simplex n-dimensional:

x ₁ = 1- R ₁ ^{1 / n-1}

x _k = (1-Σ _{i = 1} ^k x _i ) (1- R _k ^{1 / nk} ), k = 2, ..., n-1

x _n = 1-Σ _{i = 1} ^n-1 x _i

Donde R_i son números aleatorios entre 0 y 1.

Ahora estoy tratando de hacer un algoritmo para generar muestras aleatorias uniformes a partir de simplex restringido. Es decir, la intersección entre un simplex y un cuerpo convexo.