Números aleatorios con suma fija

32

Su tarea es escribir un programa o una función que genere n números aleatorios del intervalo [0,1] con una suma fija s.

Entrada

n, n≥1, número de números aleatorios para generar

s, s>=0, s<=n, suma de números a generar

Salida

Una ntupla aleatoria de números de coma flotante con todos los elementos del intervalo [0,1] y la suma de todos los elementos iguales a s, salida de cualquier manera conveniente y sin ambigüedades. Todas las ntuplas válidas tienen que ser igualmente probables dentro de las limitaciones de los números de coma flotante.

Esto es igual al muestreo uniforme desde la intersección de los puntos dentro del ncubo de la unidad dimensional y el n-1hiperplano dimensional que atraviesa (s/n, s/n, …, s/n)y es perpendicular al vector (1, 1, …, 1)(ver el área roja en la Figura 1 para ver tres ejemplos).

Ejemplos con n = 3 y sumas 0.75, 1.75 y 2.75

Figura 1: El plano de salidas válidas con n = 3 y sumas 0.75, 1.75 y 2.75

Ejemplos

n=1, s=0.8 → [0.8]
n=3, s=3.0 → [1.0, 1.0, 1.0]
n=2, s=0.0 → [0.0, 0.0]
n=4, s=2.0 → [0.2509075946818119, 0.14887693388076845, 0.9449661625992032, 0.6552493088382167]
n=10, s=9.999999999999 → [0.9999999999999,0.9999999999999,0.9999999999999,0.9999999999999,0.9999999999999,0.9999999999999,0.9999999999999,0.9999999999999,0.9999999999999,0.9999999999999]

Reglas

  • Su programa debería terminar en menos de un segundo en su máquina al menos con n≤10y cualquier s válido.
  • Si lo desea, su programa puede ser exclusivo en el extremo superior, es decir, s<ny los números de salida del intervalo medio abierto [0,1] (rompiendo el segundo ejemplo)
  • Si su idioma no admite números de coma flotante, puede falsificar la salida con al menos diez dígitos decimales después del punto decimal.
  • Las lagunas estándar no están permitidas y se permiten métodos de entrada / salida estándar.
  • Este es el , por lo que gana la entrada más corta, medida en bytes.
Angs
fuente
Relacionado.
Martin Ender
Cuando dices This is equal to uniformly sampling from the intersection: puedo ver un programa que elige al azar solo desde las esquinas de esa intersección. ¿Sería eso válido?
JayCe
2
@KevinCruijssen No, eso solo es cierto para s==0 or s==3. Para todos los demás valores de s, el plano tiene un área distinta de cero y tiene que elegir de manera uniforme y aleatoria un punto en ese plano.
user202729
3
Requerir que el intervalo esté cerrado o medio cerrado (en oposición a abierto) es un requisito teóricamente inobservable. Muchos generadores de números aleatorios dan la salida en (0,1). ¿Cómo comprobar que el intervalo de salida es [0,1) y no (0,1)? El valor 0 "nunca" aparece de todos modos
Luis Mendo
2
¿Está bien si nuestro código usa muestreo de rechazo, y por lo tanto toma mucho tiempo para casos de prueba como s=2.99999999999, n=3? ¿Podemos generar reales aleatorios en múltiplos de, digamos 1e-9,?
xnor

Respuestas:

1

Wolfram Language (Mathematica) , 92 90 bytes

If[2#2>#,1-#0[#,#-#2],While[Max[v=Differences@Sort@Join[{0,#2},RandomReal[#2,#-1]]]>1];v]&

Pruébalo en línea!

Código sin golf:

R[n_, s_] := Module[{v},
  If[s <= n/2,             (* rejection sampling for s <= n/2:                        *)
    While[
      v = Differences[Sort[
            Join[{0},RandomReal[s,n-1],{s}]]];         (* trial randoms that sum to s *)
      Max[v] > 1           (* loop until good solution found                          *)
    ];
    v,                     (* return the good solution                                *)
    1 - R[n, n - s]]]      (* for s > n/2, invert the cube and rejection-sample       *)

Aquí hay una solución que funciona en 55 bytes, pero por ahora (la versión 12 de Mathematica) está restringida a n=1,2,3porque se RandomPointniega a dibujar puntos de hiperplanos de dimensiones superiores (en la versión 11.3 de TIO también falla n=1). Sin nembargo, puede funcionar para más alto en el futuro:

RandomPoint[1&~Array~#~Hyperplane~#2,1,{0,1}&~Array~#]&

Pruébalo en línea!

Código sin golf:

R[n_, s_] :=
  RandomPoint[                           (* draw a random point from *)
    Hyperplane[ConstantArray[1, n], s],  (* the hyperplane where the sum of coordinates is s *)
    1,                                   (* draw only one point *)
    ConstantArray[{0, 1}, n]]            (* restrict each coordinate to [0,1] *)
romano
fuente
6

Python 2 , 144 128 119 bytes

from random import*
def f(n,s):
 r=min(s,1);x=uniform(max(0,r-(r-s/n)*2),r);return n<2and[s]or sample([x]+f(n-1,s-x),n)

Pruébalo en línea!


  • -20 bytes, gracias a Kevin Cruijssen
TFeld
fuente
@LuisMendo debería arreglarse ahora
TFeld
todavía no parecen uniformes
l4m2
@ l4m2 Corrí g(4, 2.0)1,000 veces para obtener 4,000 puntos y los resultados se ven así, lo que parece bastante uniforme.
Engineer Toast
122 bytes?
Giuseppe
2
Todavía mal
l4m2
4

Java 8, 194 188 196 237 236 bytes

n->s->{double r[]=new double[n+1],d[]=new double[n],t;int f=0,i=n,x=2*s>n?1:0;for(r[n]=s=x>0?n-s:s;f<1;){for(f=1;i-->1;)r[i]=Math.random()*s;for(java.util.Arrays.sort(r);i<n;d[i++]=x>0?1-t:t)f=(t=Math.abs(r[i]-r[i+1]))>1?0:f;}return d;}

+49 bytes (188 → 196 y 196 → 237) para corregir la velocidad de los casos de prueba que se acercan a 1, así como también arreglar el algoritmo en general.

Pruébalo en línea

Explicación:

Usa el enfoque de esta respuesta Stackoverflow , dentro de un bucle, siempre y cuando uno de los elementos es aún mayor que 1.
Además, si 2*s>n, sserá transformado en n-s, y se establece un indicador para indicar que debemos utilizar 1-diffen lugar de diffen el resultado de matriz (gracias por el consejo @soktinpk y @ l4m2 ).

n->s->{              // Method with integer & double parameters and Object return-type
  double r[]=new double[n+1]
                     //  Double-array of random values of size `n+1`
         d[]=new double[n],
                     //  Resulting double-array of size `n`
         t;          //  Temp double
  int f=0,           //  Integer-flag (every item below 1), starting at 0
      i=n,           //  Index-integer, starting at `n`
      x=             //  Integer-flag (average below 0.5), starting at:
        2*s>n?       //   If two times `s` is larger than `n`:
         1           //    Set this flag to 1
        :            //   Else:
         0;          //    Set this flag to 0
  for(r[n]=s=        //  Set both the last item of `r` and `s` to:
       x>0?          //   If the flag `x` is 1:
        n-s          //    Set both to `n-s`
       :             //   Else:
        s;           //    Set both to `s`
      f<1;){         //  Loop as long as flag `f` is still 0
    for(f=1;         //   Reset the flag `f` to 1
        i-->1;)      //   Inner loop `i` in range (n,1] (skipping the first item)
      r[i]=Math.random()*s;
                     //    Set the i'th item in `r` to a random value in the range [0,s)
    for(java.util.Arrays.sort(r);
                     //   Sort the array `r` from lowest to highest
        i<n;         //   Inner loop `i` in the range [1,n)
        ;d[i++]=     //     After every iteration: Set the i'th item in `d` to:
          x>0?       //      If the flag `x` is 1:
           1-t       //       Set it to `1-t`
          :          //      Else:
           t)        //       Set it to `t`
      f=(t=Math.abs( //    Set `t` to the absolute difference of:
            r[i]-r[i+1])) 
                     //     The i'th & (i+1)'th items in `r`
        >1?          //    And if `t` is larger than 1 (out of the [0,1] boundary)
         0           //     Set the flag `f` to 0
        :            //    Else:
         f;}         //     Leave the flag `f` unchanged
  return d;}         //  Return the array `d` as result
Kevin Cruijssen
fuente
Tiempo de espera paratest(10, 9.99);
l4m2
@ l4m2 Sí, noté lo mismo 10, 9.0justo después de editar para corregir el n=10, s=9.999999999999caso de prueba. No estoy seguro de si hay una solución en Java mientras aún mantiene su aleatoriedad uniforme. Tendré que pensarlo por un tiempo. Por ahora lo editaré para indicar que se agota el tiempo de espera.
Kevin Cruijssen
Si n-s<1puede llamar f(n,n-s)y voltear cada número 1/2(es decir, reemplazar xcon 1-x) como lo hizo l4m2. Esto podría resolver el problema de los números scercanos n.
soktinpk
@soktinpk Gracias por el consejo. En realidad es en s+s>nlugar de n-s<1, pero cuando miré las otras respuestas de JavaScript, realmente tenía sentido. Todo está arreglado ahora, incluido otro error que todavía estaba presente. Los bytes subieron bastante, pero todo funciona ahora. Trabajará el byte-cuenta atrás desde aquí. :)
Kevin Cruijssen
No conozco una prueba general, pero creo que este algoritmo funciona porque un hipercubo N-dimensional se puede cortar en N hiperpirámides N-dimensionales.
Neil
3

C ++ 11, 284 267 bytes

-17 bytes gracias a
la biblioteca aleatoria Zacharý Utiliza C ++, salida en la salida estándar

#include<iostream>
#include<random>
typedef float z;template<int N>void g(z s){z a[N],d=s/N;int i=N;for(;i;)a[--i]=d;std::uniform_real_distribution<z>u(.0,d<.5?d:1-d);std::default_random_engine e;for(;i<N;){z c=u(e);a[i]+=c;a[++i]-=c;}for(;i;)std::cout<<a[--i]<<' ';}

Para llamar, solo necesita hacer eso:

g<2>(0.0);

Donde el parámetro de plantilla (aquí, 2) es N, y el parámetro real (aquí, 0.0) es S

HatsuPointerKun
fuente
Creo que puedes eliminar el espacio entre <z>yu
Zacharý
Lo tengo más abajo: typedef float z;template<int N>void g(z s){z a[N],d=s/N;int i=N;for(;i;)a[--i]=d;std::uniform_real_distribution<z>u(.0,d<.5?d:1-d);std::default_random_engine e;for(;i<N;){z c=u(e);a[i]+=c;a[++i]-=c;}for(;i;)std::cout<<a[--i]<<' ';}. Una nueva línea no tiene que ser el separador entre elementos
Zacharý
1
Sugerir deshacerse por dcompleto cambiando d=s/Na s/=NSugerir volver a trabajar el segundo ciclo en for(z c;i<N;a[++i%N]-=c)a[i]+=c=u(e);lugar de for(;i<N;){z c=u(e);a[i]+=c;a[++i]-=c;}(tenga %Nen cuenta lo agregado para que el programa calcule el primer número correctamente)
Max Yekhlakov
2

Limpio , 221 201 bytes

Limpio, código de golf o números aleatorios. Elige dos.

import StdEnv,Math.Random,System._Unsafe,System.Time
g l n s#k=toReal n
|s/k>0.5=[s/k-e\\e<-g l n(k-s)]
#r=take n l
#r=[e*s/sum r\\e<-r]
|all((>)1.0)r=r=g(tl l)n s

g(genRandReal(toInt(accUnsafe time)))

Pruébalo en línea!

Función parcial literal :: (Int Real -> [Real]). Solo producirá nuevos resultados una vez por segundo.
Preciso hasta al menos 10 decimales.

Οurous
fuente
2

R , 99 bytes (con gtoolspaquete)

f=function(n,s){if(s>n/2)return(1-f(n,n-s))
while(any(T>=1)){T=gtools::rdirichlet(1,rep(1,n))*s}
T}

Pruébalo en línea!

A~={w1,,wn:i,0<wi<1;wi=s}wisA={w1,,wn:i,0<wi<1s;wi=1}

s=1Dirichlet(1,1,,1) s1<1/ss

s>n/2

Robin Ryder
fuente
2

C, 132 127 125 118 110 107 bytes

-2 bytes gracias a @ceilingcat

i;f(s,n,o,d)float*o,s,d;{for(i=n;i;o[--i]=d=s/n);for(;i<n;o[++i%n]-=s)o[i]+=s=(d<.5?d:1-d)*rand()/(1<<31);}

Pruébalo en línea!

OverclockedSanic
fuente
Desafortunadamente, esta respuesta no cumple con la especificación de desafío. Los números aleatorios de salida no están restringidos a [0,1], y su distribución conjunta no es uniforme.
Nitrodon
@Nitrodon oye, ¿podrías proporcionar una entrada para la que la salida no esté restringida a [0,1]? Probé un par de ejemplos diferentes y todos parecían ser correctos, a menos que entendiera mal el objetivo.
OverclockedSanic
Con el estado RNG en TIO, y manteniendo su n=4, los valores s=3.23y las s=0.89salidas están fuera del rango. Más concretamente, la distribución de X-s/ndebería depender de s, pero no lo hace.
Nitrodon
@Nitrodon Oops, mi mal. Estaba convirtiendo algunas partes de la respuesta C ++ anterior y olvidé agregar una cosa. Debería arreglarse ahora? También jugó unos pocos bytes en el proceso.
OverclockedSanic
1

Haskell , 122 217 208 bytes

import System.Random
r p=randomR p
(n#s)g|n<1=[]|(x,q)<-r(max 0$s-n+1,min 1 s)g=x:((n-1)#(s-x)$q)
g![]=[]
g!a|(i,q)<-r(0,length a-1)g=a!!i:q![x|(j,x)<-zip[0..]a,i/=j]
n%s=uncurry(!).(n#s<$>).split<$>newStdGen

Pruébalo en línea!

A veces las respuestas están ligeramente desviadas debido, supongo, a un error de coma flotante. Si es un problema, puedo solucionarlo a un costo de 1 byte. Tampoco estoy seguro de qué tan uniforme es esto (bastante seguro de que está bien, pero no soy tan bueno en este tipo de cosas), así que describiré mi algoritmo.

La idea básica es generar un número, xluego restarlo sy repetirlo hasta que tengamos nelementos y luego mezclarlos. Genero xcon un límite superior de 1 o s(el que sea menor) y un límite inferior de s-n+1o 0 (el mayor). Ese límite inferior está ahí, de modo que en la próxima iteración sseguirá siendo menor o igual que n(derivación: s-x<=n-1-> s<=n-1+x-> s-(n-1)<=x-> s-n+1<=x).

EDITAR: Gracias a @ michi7x7 por señalar una falla en mi uniformidad. Creo que lo solucioné con barajar pero avíseme si hay otro problema

EDIT2: recuento de bytes mejorado más restricción de tipo fijo

usuario1472751
fuente
3
El encadenamiento de muestras uniformes nunca conducirá a una distribución uniforme (la última coordenada es casi siempre mayor que 0.99 en su ejemplo)
michi7x7
@ michi7x7 Veo tu punto. ¿Qué sucede si barajé el orden de la lista después de generarla? Debería haber tomado más clases de estadísticas
user1472751
Los números no se ven muy uniformes. Aquí , 8 resultados son> 0,99, 1 es 0,96 y el último es 0,8. Así es como se ve.
Stewie Griffin
@ user1472751 hay varias buenas respuestas aquí: stackoverflow.com/q/8064629/6774250
michi7x7
1
La uniformidad todavía tiene algún problema, mira aquí : se generan demasiados ceros (gráfico de valores ordenados del 1000% 500)
Angs
1

Haskell , 188 bytes

import System.Random
import Data.List
n!s|s>n/2=map(1-)<$>n!(n-s)|1>0=(zipWith(-)=<<tail).sort.map(*s).(++[0,1::Double])<$>mapM(\a->randomIO)[2..n]>>= \a->if all(<=1)a then pure a else n!s

Sin golf:

n!s
 |s>n/2       = map (1-) <$> n!(n-s)       --If total more than half the # of numbers, mirror calculation 
 |otherwise   = (zipWith(-)=<<tail)        --Calculate interval lengths between consecutive numbers
              . sort                       --Sort
              . map(*s)                    --Scale
              . (++[0,1::Double])          --Add endpoints
              <$> mapM(\a->randomIO)[2..n] --Calculate n-1 random numbers
              >>= \a->if all(<=1)a then pure a else n!s   --Retry if a number was too large

Pruébalo en línea!

Angs
fuente