Números aleatórios com soma fixa

32

Sua tarefa é escrever um programa ou uma função que produza n números aleatórios do intervalo [0,1] com soma fixa s.

Entrada

n, n≥1, número de números aleatórios a serem gerados

s, s>=0, s<=n, soma dos números a serem gerados

Saída

Um número aleatório nde números de ponto flutuante com todos os elementos do intervalo [0,1] e a soma de todos os elementos iguais a s, são emitidos de qualquer maneira conveniente e inequívoca. Todos os n-tuplos válidos devem ter a mesma probabilidade dentro das limitações dos números de ponto flutuante.

Isso é igual à amostragem uniforme da interseção dos pontos dentro do ncubo da unidade tridimensional e do n-1hiperplano tridimensional que passa (s/n, s/n, …, s/n)e é perpendicular ao vetor (1, 1, …, 1)(veja a área vermelha na Figura 1 para três exemplos).

Exemplos com n = 3 e somas 0,75, 1,75 e 2,75

Figura 1: O plano de saídas válidas com n = 3 e soma 0,75, 1,75 e 2,75

Exemplos

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]

Regras

  • Seu programa deve terminar em menos de um segundo em sua máquina pelo menos com n≤10e quaisquer s válidos.
  • Se desejar, seu programa pode ser exclusivo na extremidade superior, ou seja, s<nos números de saída do intervalo semiaberto [0,1) (quebrando o segundo exemplo)
  • Se o seu idioma não suportar números de ponto flutuante, você poderá falsificar a saída com pelo menos dez dígitos decimais após o ponto decimal.
  • As brechas padrão não são permitidas e os métodos padrão de entrada / saída são permitidos.
  • Isso é , então a entrada mais curta, medida em bytes, vence.
Angs
fonte
Relacionado.
Martin Ender
Quando você diz This is equal to uniformly sampling from the intersection- eu posso ver um programa escolhendo aleatoriamente apenas nos cantos do cruzamento. Isso seria válido?
21418 JayCe
2
@KevinCruijssen Não, isso é verdade apenas para s==0 or s==3. Para todos os outros valores de s, o plano tem área diferente de zero e você precisa escolher aleatoriamente um ponto nesse plano.
User202729 18/0118
3
Exigir que o intervalo seja fechado ou semi-fechado (ao contrário de aberto) é um requisito teoricamente não observável. Muitos geradores de números aleatórios fornecem a saída em (0,1). Como testar se o intervalo de saída é [0,1) e não (0,1)? O valor 0 "nunca" ocorre de qualquer maneira
Luis Mendo
2
Tudo bem se nosso código usa amostragem por rejeição e leva muito tempo para casos de teste como s=2.99999999999, n=3? Podemos gerar reais aleatórios em múltiplos de, digamos 1e-9,?
Xnor

Respostas:

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]&

Experimente online!

Código não golfe:

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       *)

Aqui está uma solução que funciona em 55 bytes, mas por enquanto (Mathematica versão 12) é restrita n=1,2,3porque RandomPointse recusa a extrair pontos de hiperplanos de dimensões mais altas (na versão 11.3 do TIO também falha n=1). Pode funcionar para mais alto nno futuro:

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

Experimente online!

Código não golfe:

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
fonte
7

JavaScript (Node.js) , 122 115 bytes

N=>W=S=>2*S>N?W(N-S).map(t=>1-t):(t=(Q=s=>n?[r=s-s*Math.random()**(1/--n),...r>1?[++Q]:Q(s-r)]:[])(S,n=N),Q?t:W(S))

Experimente online!

l4m2
fonte
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)

Experimente online!


  • -20 bytes, graças a Kevin Cruijssen
TFeld
fonte
@LuisMendo Deve ser corrigido agora #
TFeld
eles ainda não parecem uniformes
l4m2 18/18
@ l4m2 Corri g(4, 2.0)1.000 vezes para obter 4.000 pontos e os resultados são assim, o que parece bastante uniforme.
Engineer Toast
122 bytes?
Giuseppe
2
Ainda está errado
l4m2 23/05
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 e 196 → 237) para corrigir a velocidade dos casos de teste próximos a 1, bem como para corrigir o algoritmo em geral.

Experimente online

Explicação:

Usa a abordagem de esta resposta stackoverflow , dentro de um loop, enquanto um dos itens é ainda maior do que 1.
Além disso, se 2*s>n, sserá alterado em n-s, e um sinalizador está definido para indicar que devemos usar 1-diffem vez de diffno resultado-array (obrigado pela dica @soktinpk e @ 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
fonte
Tempo limite paratest(10, 9.99);
l4m2
@ l4m2 Sim, notei o mesmo 10, 9.0logo após a edição para corrigir o n=10, s=9.999999999999caso de teste .. Não tenho certeza se existe uma correção em Java enquanto ainda mantém sua aleatoriedade uniforme .. Terá que pensar nisso por um tempo. Por enquanto, vou editá-lo para declarar o tempo limite.
Kevin Cruijssen 18/05/19
Se n-s<1você pode ligar f(n,n-s)e virar todos os números 1/2(ou seja, substituir xpor 1-x), como l4m2 fez. Isso pode resolver o problema de números spróximos a n.
Soktinpk 19/05/19
@soktinpk Obrigado pela dica. Na verdade, é em s+s>nvez de n-s<1, mas quando olhei para as outras respostas do JavaScript, de fato fazia sentido. Tudo está consertado agora, incluindo outro bug que ainda estava presente. Os bytes subiram um pouco, mas agora todos funcionam. Trabalhará a contagem de bytes a partir daqui. :)
Kevin Cruijssen
Não conheço uma prova geral, mas acredito que esse algoritmo funcione porque um hipercubo N-dimensional pode ser cortado em hiper-pirâmides N-dimensionais.
Neil
3

JavaScript (Node.js) , 153 bytes

(n,s)=>s+s>n?g(n,n-s).map(r=>1-r):g(n,s)
g=(n,s)=>{do(a=[...Array(n-1)].map(_=>Math.random()*(s<1?s:1))).map(r=>t-=r,t=s);while(t<0|t>=1);return[...a,t]}

Experimente online!

Neil
fonte
3

C ++ 11, 284 267 bytes

-17 bytes graças ao Zacharý
Usa biblioteca aleatória C ++, saída na saída padrão

#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 ligar, você só precisa fazer isso:

g<2>(0.0);

Onde o parâmetro do modelo (aqui, 2) é N e o parâmetro real (aqui, 0,0) é S

HatsuPointerKun
fonte
Eu acho que você pode remover o espaço entre <z>eu
Zachary
Eu tenho-la ainda mais: 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]<<' ';}. A nova linha não tem que ser o separador entre os itens
Zachary
11
Sugiro se livrar dcompletamente, alterando d=s/Na s/=NSugerir reformulação da segunda volta for(z c;i<N;a[++i%N]-=c)a[i]+=c=u(e);, em vez de for(;i<N;){z c=u(e);a[i]+=c;a[++i]-=c;}(note o agregado %N, a fim de tornar a calcular programa o primeiro número corretamente)
Max Yekhlakov
2

Limpo , 221 201 bytes

Números limpos, de código de golfe ou aleatórios. Escolhe dois.

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)))

Experimente online!

Função parcial literal :: (Int Real -> [Real]). Só produzirá novos resultados uma vez por segundo.
Precisas até pelo menos 10 casas decimais.

Furioso
fonte
2

R , 99 bytes (com gtoolspacote)

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}

Experimente online!

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
fonte
2

C, 132 127 125 118 110 107 bytes

-2 bytes graças 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);}

Experimente online!

OverclockedSanic
fonte
Infelizmente, esta resposta não atende à especificação de desafio. Os números aleatórios de saída não são restritos [0,1]e sua distribuição conjunta não é uniforme.
Nitrodon 16/07
@Nitrodon ei, você poderia fornecer uma entrada para a qual a saída não está restrita a [0,1]? Tentei alguns exemplos diferentes e todos pareciam estar corretos, a menos que eu não entendesse o objetivo.
OverclockedSanic
Com o estado RNG no TIO, e mantendo seus n=4, os valores s=3.23e as s=0.89saídas fora do intervalo. Mais precisamente, a distribuição de X-s/ndeve depender s, mas não.
Nitrodon 17/07
@ Nitrodon Opa, meu mal. Eu estava convertendo algumas partes da resposta C ++ acima e esqueci de adicionar uma coisa. Deve ser corrigido agora? Também jogou alguns bytes no processo.
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

Experimente online!

Às vezes, as respostas estão um pouco erradas devido, presumo, a um erro de ponto flutuante. Se for um problema, posso corrigi-lo a um custo de 1 byte. Também não tenho certeza de quão uniforme isso é (tenho certeza de que está bem, mas não sou tão bom nesse tipo de coisa), então vou descrever meu algoritmo.

A idéia básica é gerar um número x, subtraí-lo se repetir até que tenhamos nelementos e os embaralhe. Gero xcom um limite superior de 1 ou s(o que for menor) e um limite inferior de s-n+1ou 0 (o que for maior). Esse limite inferior existe para que na próxima iteração sainda seja menor ou igual a n(derivação: s-x<=n-1-> s<=n-1+x-> s-(n-1)<=x-> s-n+1<=x).

EDIT: Obrigado a @ michi7x7 por apontar uma falha na minha uniformidade. Acho que o corrigi com a reprodução aleatória, mas deixe-me saber se há outro problema

EDIT2: contagem de bytes aprimorada mais restrição de tipo fixo

user1472751
fonte
3
Encadeamento amostras uniformes nunca vai levar a uma distribuição uniforme (a última coordenada quase sempre maior do que 0,99 é no seu exemplo)
michi7x7
@ michi7x7 Entendo o seu ponto. E se eu embaralhar a ordem da lista depois de gerá-la? Eu deveria ter feito mais aulas de estatísticas
user1472751
Os números não parecem muito uniformes. Aqui , 8 resultados são> 0,99, 1 é 0,96 e o ​​último é 0,8. Isso é o que parece.
Stewie Griffin
@ user1472751 existem várias respostas boas aqui: stackoverflow.com/q/8064629/6774250
michi7x7
11
Uniformness ainda tem algum problema, veja aqui - há muitos zeros gerados (lote de valores ordenados a partir de 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

Ungolfed:

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

Experimente online!

Angs
fonte