{% comment %}Custom stylesheet, it must be in the same directory as the html file {% endcomment %} {% comment %}Loading mathjax macro {% endcomment %} {% comment %}Load mathjax {% endcomment %} {% comment %}MathJax configuration {% endcomment %} {% comment %}End of mathjax configuration {% endcomment %}

Sorteios de bolas em baldes: com e sem restrições e degenerescência

Carregando o pacote de fazer gráficos:

In [1]:
using Plots

Serão sorteadas 20 bolas, onde cada uma delas pode estar em um de 20 baldes:

In [2]:
nbolas = 20;
In [3]:
nbaldes = 20;

Um sorteio vai consistir em sortear em que balde está cada bola. Isto é feito com o comando que gera um vetor de nbolas posições, com números que variam aleatoriamente de 1 a nbaldes:

In [4]:
balde = rand(1:nbaldes,nbolas)
Out[4]:
20-element Array{Int64,1}:
 19
 19
  1
 20
 18
 12
 18
 19
 20
  7
  2
 20
 18
 12
 10
 13
  1
  6
 12
  7

Vamos começar fazendo 100 sorteios:

In [5]:
nsorteios = 100;

O resultado final será guardado no vetor nbolas_em_cada_balde, no qual vamos acumular quantas bolas cairam em cada balde. Portanto, é um vetor com nbaldes posições, que incialmente tem zero em todas posições:

In [6]:
nbolas_em_cada_balde = zeros(Int,nbaldes);

Agora vamos especificamente fazer os sorteios. Para cada sorteio geramos o vetor aleatório balde e, em seguida, contamos quantas vezes o número associado a cada balde apareceu nesse vetor aleatório. Essa contagem é feita com o comando count(isequal(5),balde) que, neste exemplo, contaria quantas vezes o número 5 aparece no vetor balde.

In [7]:
for i in 1:nsorteios
  balde = rand(1:nbaldes,nbolas)
  for j in 1:nbaldes
    nbolas_em_cada_balde[j] = nbolas_em_cada_balde[j] + count(isequal(j),balde)
  end
end

Por fim, vamos fazer a média, já que somamos o resultado de todos os sorteios. Vamos fazer o número médio de bolas que foi observado em cada balde, nesse conjunto de sorteios:

In [8]:
nbolas_em_cada_balde = nbolas_em_cada_balde / nsorteios ;

E, finalmente, fazemos um gráfico de barras do resultado:

In [9]:
bar(nbolas_em_cada_balde)
Out[9]:
0 5 10 15 20 0.0 0.2 0.4 0.6 0.8 1.0 1.2 y1

Como temos 20 bolas e 20 baldes, esperamos que em media tenhamos 1 bola por balde no sorteio, o que é aproximadamente o que foi observado. Essa media vai ser tão mais próxima de 1 quanto maior o número de sorteios. A população de número de bolas médio por balde converge para o valor esperado.

Vamos ver isso, primeiro, aumentando o número de sorteios para 1000, e repetindo o procedimento:

In [10]:
nsorteios = 1000;
nbolas_em_cada_balde = zeros(Int,nbaldes);
for i in 1:nsorteios
  balde = rand(1:nbaldes,nbolas)
  for j in 1:nbaldes
    nbolas_em_cada_balde[j] = nbolas_em_cada_balde[j] + count(isequal(j),balde)
  end
end
nbolas_em_cada_balde = nbolas_em_cada_balde / nsorteios ;
bar(nbolas_em_cada_balde)
Out[10]:
0 5 10 15 20 0.00 0.25 0.50 0.75 1.00 y1

E, em seguida, aumetnamos o número de sorteios para 1000.000 (um milhão), e temos um resultado que é praticamente o resultado esperado, apenas com pequenas flutuações:

In [11]:
nsorteios = 1000000;
nbolas_em_cada_balde = zeros(Int,nbaldes);
for i in 1:nsorteios
  balde = rand(1:nbaldes,nbolas)
  for j in 1:nbaldes
    nbolas_em_cada_balde[j] = nbolas_em_cada_balde[j] + count(isequal(j),balde)
  end
end
nbolas_em_cada_balde = nbolas_em_cada_balde / nsorteios ;
bar(nbolas_em_cada_balde)
Out[11]:
0 5 10 15 20 0.00 0.25 0.50 0.75 1.00 y1

Sorteios com restrições

Agora vamos ver o que acontece quando introduzimos uma restrição sobre o sorteio. Vamos imaginar que o balde 1 tem energia 1, o balde 2 tem energia 2, o balde 3 tem energia 3, e assim sucessivamente.

Desta forma, a "energia" total é a soma das energias dos baldes, que neste caso é simplesmente a soma dos números dos baldes que apareceram no sorteio.

Por exemplo, neste sorteio:

In [12]:
balde = rand(1:nbaldes,nbolas)
Out[12]:
20-element Array{Int64,1}:
 13
  2
 15
 14
  6
 20
 16
  1
 10
  9
 19
 20
  8
 16
 20
 17
 12
 17
 15
  2

A energia total é:

In [13]:
sum(balde)
Out[13]:
252

Vamos então introduzir uma restrição, e só aceitar os sorteios em que a energia esteja entre 160 e 180. Queremos 1000 sorteios com essas características, mas cada sorteio pode ou não satisfazê-las. Portanto, vamos mudar um pouco a lógica do loop que faz os sorteios.

O loop será um while, ou seja, "enquanto", i for menor que nsorteios. i começa valendo 1, e só aumenta se o sorteio foi válido, isto é, se a energia caiu no intervalo desejado. Esta condição é testada pelo if abaixo (notando que && significa "e"):

In [22]:
nsorteios = 1000;
nbolas_em_cada_balde = zeros(Int,nbaldes);
i = 1
while i < nsorteios
  # O sorteio é feito aqui:
  balde = rand(1:nbaldes,nbolas)
  # Se a energia está no intervalo ... 
  if sum(balde) > 160 && sum(balde) < 180
    # Atualizamos o resultado como antes:
    for j in 1:nbaldes
      nbolas_em_cada_balde[j] = nbolas_em_cada_balde[j] + count(isequal(j),balde)
    end
    # E, somente neste caso, aumentamos o valor de i    
    i = i + 1 
  end
end

Finalmente, vamos calcular a media e fazer o gráfico de barras:

In [23]:
nbolas_em_cada_balde = nbolas_em_cada_balde / nsorteios ;
bar(nbolas_em_cada_balde)
Out[23]:
0 5 10 15 20 0.0 0.5 1.0 1.5 y1

Note que a restrição sobre a energia fez com o que os baldes de menor energia aparecessem com mais frequência que os baldes de maior energia. A restrição sobre a energia total, portanto, induziu uma mudança na população.

Acrescentando degenerescência:

Vamos agora supor que há 1 balde do tipo 1, 2 baldes do tipo 2, 3 baldes do tipo 3, e assim por diante. Portanto, há mais baldes quanto maior o número de baldes.

Para fazer o sorteio, agora, vamos definir um vetor que tem o valor dos baldes. São muitos baldes. Mais especificamente, são $$1 + 2 + 3 + 4 + 5 + \cdots + 20 = \frac{21\times 20}{2} = 210 $$ baldes (a soma acima é facil de fazer se percerbermos que estamos somando todos os elementos da parte triangular superior, sem a diagonal, de uma matriz de dimensão $21\times 21$).

Vamos criar um vetor de 210 posições, com o valor de cada balde em cada posição. Isto pode ser feito como segue. Primeiro, criamos um vetor de números inteiros com 210 posições:

In [24]:
nbaldes = 210;
In [25]:
tipo_de_balde = zeros(Int,nbaldes);

E, agora, vamos preenchar cada posição com o valor dos baldes, de acordo com o que decidimos:

In [26]:
ibalde = 0
for i in 1:20 # este loop é sobre os "tipos" de balde
  for j in 1:i # este loop é sobre o valor dos baldes, vai de 1 até o tipo.
    ibalde = ibalde + 1  # este contador vai de 1 a 210
    tipo_de_balde[ibalde] = i   # o valor do balde corresponde ao tipo
  end
end

Podemos conferir que o vetor tipo_de_baldes tem a informação que queremos:

In [27]:
tipo_de_balde
Out[27]:
210-element Array{Int64,1}:
  1
  2
  2
  3
  3
  3
  4
  4
  4
  4
  5
  5
  5
  ⋮
 20
 20
 20
 20
 20
 20
 20
 20
 20
 20
 20
 20

Note, como esperado, que há 1 elemento com valor 1, 2 com valor 2, etc., como queríamos.

O sorteio agora vai consistir em sortear as 20 bolas não entre 20 baldes, como antes, mas entre 210 baldes. Só que a restrição sobre a soma dos valores vai continuar a mesma, associada ao valor de cada balde, dado pelo vetor acima. São 20 tipos de baldes,

In [28]:
ntipos = 20;

Um sorteio consiste em gerar nbolas (que é 20, aqui) números inteiros aleatórios entre 1 e 210:

In [29]:
balde = rand(1:nbaldes,nbolas)
Out[29]:
20-element Array{Int64,1}:
 140
 182
   5
  43
  60
  39
 129
 142
  92
 145
 147
  47
  24
  83
 167
  94
   3
 163
  27
  86

O resultado acima é um vetor, de tal forma, por exemplo, que o balde do elemento 5 do vetor é:

In [30]:
balde[5]
Out[30]:
60

O tipo de balde deste elemento é dado pelo vetor tipo_de_balde:

In [31]:
tipo_de_balde[60]
Out[31]:
11

A notação acima pode ser generalizada, no sentido de que podemos direamente recuperar o tipo de balde do balde 5 usando:

In [32]:
tipo_de_balde[balde[5]]
Out[32]:
11

Mais ainda, podemos definir um vetor, que chamaremos tipo, que contém o tipo de todos os baldes do sorteio:

In [33]:
tipo = [ tipo_de_balde[balde[j]] for j in 1:nbolas ]
Out[33]:
20-element Array{Int64,1}:
 17
 19
  3
  9
 11
  9
 16
 17
 14
 17
 17
 10
  7
 13
 18
 14
  2
 18
  7
 13

A soma das energias pode ser calculada, agora, pela soma dos tipos:

In [34]:
sum(tipo)
Out[34]:
251

Quando fizermos muitos sorteios, vamos guardar a informação do resultado em um vetor que vai contar quantas vezes foi observado cada tipo de balde. Este vetor começa com zeros em todas as posições, e cada vez que uma bola cair em um balde de cada tipo, esta informação será adicionada:

In [35]:
nbolas_em_cada_tipo_de_balde = zeros(Int,ntipos);

Por fim, vamos fazer os sorteios. Vamos sortear em que balde está cada bola, mas a restrição continuará sobre a soma dos tipos:

In [36]:
nsorteios = 1000;
i = 1
while i < nsorteios
  balde = rand(1:nbaldes,nbolas) # gerando 20 números aleatórios entre 1 e 210
  tipos = [ tipo_de_balde[balde[j]] for j in 1:nbolas ] # Este é o vetor que contém os tipos
  # Se a energia está no intervalo, mas a soma é sobre os tipos: 
  if sum(tipos) > 160 && sum(tipos) < 180
    # Atualizamos o resultado como antes:
    for j in 1:ntipos
      nbolas_em_cada_tipo_de_balde[j] = nbolas_em_cada_tipo_de_balde[j] + count(isequal(j),tipos)
    end
    # E, somente neste caso, aumentamos o valor de i    
    i = i + 1 
  end
end

Como antes, vamos fazer a média por sorteio, e fazer o gráfico de barras:

In [38]:
nbolas_em_cada_tipo_de_balde = nbolas_em_cada_tipo_de_balde / nsorteios ;
bar(nbolas_em_cada_tipo_de_balde)
Out[38]:
0 5 10 15 20 0.0000 0.0005 0.0010 0.0015 y1

O aumento do número de baldes de cada tipo levou ao aparecimento de um máximo na distribuição. Baldes de valores muito baixo são pouco prováveis porque sua degenerescência é pequena (há poucos baldes desses tipos). Baldes de valores altos são pouco prováveis pela restrição sobre a energia total.