Primeira vez aqui? Seja bem vindo e cheque o FAQ!
x

Após \( N \) variáveis aleatórias, independentes e uniformemente distribuídas entre 0 e 1, forem classificadas em ordem não-decrescente, qual a probabilidade do r-ésimo menor valor ser \( \leq x \)?

0 votos
96 visitas
perguntada Out 18 em Estatística por VITOR B BORGES (1 ponto)  
editado Out 18 por VITOR B BORGES

A pergunta corresponde ao exercício 7 localizado na página 7 do Capítulo 5 - Sorting and Searching do livro "The Art of Computer Programming - Volume 3" de Donald Ervin Knuth.

Compartilhe

2 Respostas

0 votos
respondida Out 18 por VITOR B BORGES (1 ponto)  
editado Out 22 por VITOR B BORGES

Temos N variáveis aleatórias \(Z\) que seguem a distribuição \( Z \sim U[0, 1] \). Em notação de sequência: \( \{Z_{i}\}_{i = 1}^{N} \).
Vamos definir uma nova variável aleatória \( Y \), que depende de \( Z \) e segue a seguinte distribuição de Bernoulli:

\[ Y(Z) \sim Bernoulli \left\{ \begin{matrix} 1,&se& Z\leq x\\ 0,&se& Z > x \end{matrix} \right\} \]

Note que \( P(Y = 1) = P(Z \leq x) = x \) e \( P(Y = 0) = P(Z > x) = 1 - x \), devido à distribuição de \(Z\).
Considere a seguinte sequência \( \{Y_{i}\}_{i = 1} ^{N}\) de zeros e uns. Se ordenarmos a sequência \( \{Z_{i}\}_{i = 1}^{N} \) de maneira não decrescente, a sequência \( \{Y_{i}\}_{i = 1} ^{N}\) será o número 1 aparecendo algumas vezes, seguido de zeros. Esta sequência começa a ser 0 assim que o k-ésimo valor \( \{Z_{i}\}_{i = 1}^{N} \) fica maior que x.

O r-ésimo termo de \( \{Z_{i}\}_{i = 1}^{N} \) ordenada será menor que x se a sequência \( \{Y_{i}\}_{i = 1} ^{N}\) tiver mais uns que o valor de r. A quantidade de uns em \( \{Y_{i}\}_{i = 1}^{N} \) é \( \sum\{Y_{i}\}_{i = 1} ^{N}\).
Como todas essas Bernoullis são independentes e identicamente distribuídas, podemos usar o fato de que um somatório de variáveis de Bernoulli segue a distribuição Binomial:
\[ \sum\{Y_{i}\}_{i = 1} ^{N} = S \sim Binomial(N, x) \]
Estamos interessados em computar \( P(S > r )\) que é igual à
\[ P(S > r ) = \sum_{j = r + 1}^{N}P(S = j) =\\ \sum_{j = r + 1}^{N} \binom{N}{j} x^{j}(1 - x)^{N - j} \]

0 votos
respondida Out 18 por VITOR B BORGES (1 ponto)  
editado Nov 5 por VITOR B BORGES

Com finalidade de confirmar este resultado, foram implementadas duas funções na linguagem Python. A primeira repete o experimento 10 mil vezes e computa a proporção de experimentos em que o r-ésimo termo foi menor que x. A segunda utiliza o formato funcional encontrado.

Primeira:

import random
import scipy

random.seed(13)

def montecarlo(N, r, x):

experiment = []
#vai guardar 1 se o r do for menor que x e 0 se r for maior.

for i in range(1000):

    values = []
    #Esta lista vai guardar os valores de Z

    for i in range(N):

        values.append(random.uniform(0, 1))
        # distribuição Uniforme

    values = sorted(values)
    #ordenar os resultados

    if values[r] <= x:

        experiment.append(1)

    else:

        experiment.append(0)

return sum(experiment) / len(experiment) #proporção

A Segunda:

def functional(N, r, x):

p = 0

for i in range(r + 1, N + 1):

    p += scipy.special.binom(N, i)*(x**i)*((1 - x)**(N - i))

return p

Para avaliarmos o resultado, calcularemos o erro médio entre a implementação de força bruta e o formato funcional encontrado.

def result_efficiency():

results_1= [] #resultados de força bruta
results_2 = [] #resultados do formato funcional

for i in range(1000):

    randints = [random.randint(0, 1000), random.randint(0, 1000)]
    #sortear dois números

    N = max(randints) #o maior número é N
    r = min(randints) # o menor número é r
    x = random.uniform(0, 1)

    results_1.append(montecarlo(N, r, x))
    results_2.append(functional(N, r, x))
    #guardar os resultados

abserror_sum = [abs(z) for z in (np.array(results_1) - np.array(results_2))]
#erro absoluto
erro_medio = sum(abserror_sum)/len(abserror_sum)
#erro medio

return round(erro_medio, 5)

O erro médio absoluto encontrado foi de 0.00205.

comentou Nov 2 por Alan Antunes Rosendo (41 pontos)  
Ótima solução a idéia de comparar a resposta utilizando monte carlo foi uma boa idéia.
Segue algumas sugestões para aprimoramento do código:
É simples, mas você esqueceu de publicar os códigos importando os módulos utilizados;
Sempre que utilizo código com comandos aleatórios utilizo o comando random.seed(num) para que seja possível a reprodução do resultado exato por um terceiro.
comentou Nov 5 por VITOR B BORGES (1 ponto)  
Muito obrigado pelo comentário Alan, ja corrigi aqui as suas sugestões!
...