2012-06-12 25 views
9

Espero crear 3 números cuasialeatorios (no negativos) que sumen a uno, y repita una y otra vez.Generar 3 números aleatorios que sumen a 1 en R

Básicamente estoy intentando dividir algo en tres partes aleatorias en muchas pruebas.

Aunque soy consciente de

a = runif (3,0,1)

Estaba pensando que podría utilizar 1-A como máximo en la siguiente ejecución si, pero parece desordenado .

Pero estos por supuesto no suman a uno. ¿Algún pensamiento, oh sabio stackoverflowers?

+2

¿Es una opción para renormalizar los números aleatorios después de la generación? –

+0

¿Qué hay de generar 2 números aleatorios a y b? Entonces a + b + c = 1 => c = 1 - (a + b) –

+0

y si ayb suman a mayor que 1? – mmann1123

Respuesta

9

simplemente al azar 2 dígitos de (0, 1) y si asuma su a y b entonces usted tiene:

rand1 = min(a, b) 
rand2 = abs(a - b) 
rand3 = 1 - max(a, b) 
+0

Además, debe repetir la generación del segundo número si a == b ... (debe ser MUY raro) – ddzialak

+0

@user so a = 0.85 , b = 0.99 luego tienes números: 0.85, 0.14, 0.01 (en cuanto a mí estos son muy buenos 3 números aleatorios de 0..1) – ddzialak

+3

La distribución resultante parece no ser exactamente trivial: http: //www.jstor. org/discover/10.2307/2983572? uid = 2129 & uid = 2 & uid = 70 & uid = 4 & sid = 21100849643501 y un documento posterior al que se puede acceder libremente http://doc.utwente.nl/70657/1/Sleutel67random.pdf – Christian

4

supongo que depende de lo que la distribución que desee en los números, pero aquí es una manera de :

diff(c(0, sort(runif(2)), 1)) 

uso replicate a conseguir tantos conjuntos como desee:

> x <- replicate(5, diff(c(0, sort(runif(2)), 1))) 
> x 
      [,1]  [,2]  [,3]  [,4]  [,5] 
[1,] 0.66855903 0.01338052 0.3722026 0.4299087 0.67537181 
[2,] 0.32130979 0.69666871 0.2670380 0.3359640 0.25860581 
[3,] 0.01013117 0.28995078 0.3607594 0.2341273 0.06602238 
> colSums(x) 
[1] 1 1 1 1 1 
11

Esta pregunta involucra cuestiones más sutiles de lo que podría parecer al principio. Después de ver el siguiente, es posible que desee pensar cuidadosamente sobre el proceso que está utilizando estos números para representar:

## My initial idea (and commenter Anders Gustafsson's): 
## Sample 3 random numbers from [0,1], sum them, and normalize 
jobFun <- function(n) { 
    m <- matrix(runif(3*n,0,1), ncol=3) 
    m<- sweep(m, 1, rowSums(m), FUN="/") 
    m 
} 

## Andrie's solution. Sample 1 number from [0,1], then break upper 
## interval in two. (aka "Broken stick" distribution). 
andFun <- function(n){ 
    x1 <- runif(n) 
    x2 <- runif(n)*(1-x1) 
    matrix(c(x1, x2, 1-(x1+x2)), ncol=3) 
} 

## ddzialak's solution (vectorized by me) 
ddzFun <- function(n) { 
    a <- runif(n, 0, 1) 
    b <- runif(n, 0, 1) 
    rand1 = pmin(a, b) 
    rand2 = abs(a - b) 
    rand3 = 1 - pmax(a, b) 
    cbind(rand1, rand2, rand3) 
} 

## Simulate 10k triplets using each of the functions above 
JOB <- jobFun(10000) 
AND <- andFun(10000) 
DDZ <- ddzFun(10000) 

## Plot the distributions of values 
par(mfcol=c(2,2)) 
hist(JOB, main="JOB") 
hist(AND, main="AND") 
hist(DDZ, main="DDZ") 

enter image description here

+0

Bien, estaba pensando en trazar los resultados, pero ya lo hiciste. Es interesante ver que aparentemente ninguna de las soluciones realmente hace lo que a uno le hubiera gustado intuitivamente. También es interesante que en estas tramas no se puede ver realmente que DDZ hace lo correcto según los medios mientras AND ni siquiera eso. – Christian

6

Cuando se desea generar al azar los números que se suman a 1 (o algún otro valor) entonces debería mirar el Dirichlet Distribution.

hay una función rdirichlet en el paquete gtools y funcionando RSiteSearch('Dirichlet') nos lleva a un buen número de éxitos que fácilmente podrían conducir a herramientas para hacer esto (y no es difícil de código a mano, ya sea para las distribuciones de Dirichlet simples).

2

Este problema y las diferentes soluciones propuestas me intrigaron. Hice una pequeña prueba de los tres algoritmos básicos sugeridos y qué valores promedio arrojarían para los números generados.

choose_one_and_divide_rest 
means:    [ 0.49999212 0.24982403 0.25018384] 
standard deviations: [ 0.28849948 0.22032758 0.22049302] 
time needed to fill array of size 1000000 was 26.874945879 seconds 

choose_two_points_and_use_intervals 
means:    [ 0.33301421 0.33392816 0.33305763] 
standard deviations: [ 0.23565652 0.23579615 0.23554689] 
time needed to fill array of size 1000000 was 28.8600130081 seconds 

choose_three_and_normalize 
means:    [ 0.33334531 0.33336692 0.33328777] 
standard deviations: [ 0.17964206 0.17974085 0.17968462] 
time needed to fill array of size 1000000 was 27.4301018715 seconds 

El tiempo de las mediciones se deben tomar con un grano de sal, ya que podrían ser más influenciados por la gestión de memoria Python que por el propio algoritmo. Soy demasiado flojo para hacerlo correctamente con timeit. Hice esto en un átomo de 1 GHz, así que eso explica por qué tardó tanto.

De todos modos, choose_one_and_divide_rest es el algoritmo sugerido por Andrie y el póster de la pregunta (AND): elige un valor a en [0,1], luego uno en [a, 1] y luego mira lo que te queda Se suma a uno, pero eso es todo, la primera división es dos veces más grande que las otras dos. Uno podría haber adivinado tanto ...

choose_two_points_and_use_intervals es la respuesta aceptada por ddzialak (DDZ). Toma dos puntos en el intervalo [0,1] y usa el tamaño de los tres subintervalos creados por estos puntos como los tres números. Funciona como un encanto y los medios son todos 1/3.

choose_three_and_normalize es la solución de Anders Gustafsson y Josh O'Brien (JOB). Simplemente genera tres números en [0,1] y los normaliza de nuevo a una suma de 1. Funciona igual de bien y sorprendentemente un poco más rápido en mi implementación de Python. La varianza es un poco menor que para la segunda solución.

Ahí lo tienes. No tengo idea de qué distribución beta corresponden estas soluciones o qué conjunto de parámetros en el documento correspondiente al que me refería en un comentario, pero tal vez alguien más pueda descifrarlo.

Cuestiones relacionadas