Solução dos desafios 9 e 10
Desafio 9: Considere que cada entidade gerada no primeiro exemplo desta seção tem um peso em gramas dado por uma distribuição normal de média 10 e desvio padrão igual a 3. Crie um critério de parada para quando a média dos pesos das entidades geradas esteja no intervalo entre 9,5 e 10,5.
Este primeiro desafio envolve poucas modificações no programa original. Acrescentamos três variáveis novas: media, contador e pesoTotal; o laço while foi substituído pelo critério de parada e algumas linhas foram acrescentadas para o cálculo da média de peso até a última entidade gerada. O peso de cada entidade é sorteado pela função random.normalvariate(mu, sigma) da biblioteca random.
import simpy
import random
def geraChegada(env, numEntidades):
media, contador, pesoTotal = 0, 0, 0
while media > 10.5 or media < 9.5: # critério de parada
yield env.timeout(1)
contador += 1 # conta entidades geradas
peso = random.normalvariate(10, 3) # sorteia o peso da entidade
pesoTotal += peso # acumula o peso total até agora
media = pesoTotal/contador # calcula média dos pesos
print("%4.1f nova chegada\tPeso: %4.1f kg\tMédia atual: %4.1f"
%(env.now, peso, media))
random.seed(100)
env = simpy.Environment()
chegadas = env.process(geraChegada(env, 5)) # gere apenas 5 entidades
env.run() # executa até o fim de todos os processos
Quando executado, o modelo anterior apresenta como resultado:
1.0 nova chegada Peso: 6.7 kg Média atual: 6.7
2.0 nova chegada Peso: 14.7 kg Média atual: 10.7
3.0 nova chegada Peso: 12.1 kg Média atual: 11.2
4.0 nova chegada Peso: 13.3 kg Média atual: 11.7
5.0 nova chegada Peso: 6.0 kg Média atual: 10.6
6.0 nova chegada Peso: 13.5 kg Média atual: 11.0
7.0 nova chegada Peso: 5.8 kg Média atual: 10.3
Desafio 10: Modifique o critério anterior para que a parada ocorra quando a média for 10 kg, com um intervalo de confiança de amplitude 0,5 e nível de significância igual a 95%. Dica: utilize a biblioteca
numpypara isso (consulte o Stack Overflow.
Esta situação exige um pouco mais no processo de codificação, contudo é algo muito utilizado em modelos de simulação de eventos discretos.
Como agora queremos o Intervalo de Confiança de uma dada amostra, os valores dos pesos serão armazenados em uma lista (pesosList, no caso do desafio).
A biblioteca numpy fornece um meio fácil de se estimar a média e o desvio padrão de uma amostra de valores armazenada numa lista:
numpy.mean(pesosList): estima a média da listapesosList;numpy.std(pesosList): estima o desvio-padrão da listapesosList
Para o cálculo do intervalo de confiança, devemos lembrar que, para amostras pequenas, a sua expressão é dada por:
$$ [\bar{x}- t{1-\alpha , n-1} \frac{s}{\sqrt{n}} , \bar{x} + t{1-\alpha , n-1} \frac{s}{\sqrt{n}}]
$$
A biblioteca scipy.stats possui diversas funções estatísticas, dentre elas, a distribuição t de student, necessária para o cálculo do intervalo de confiança. Como está será uma operação rotineira nos nossos modelos de simulação, o ideal é encapsular o código em uma função específica:
def intervaloConfMedia(a, conf=0.95):
# retorna a média e a amplitude do intervalo de confiança dos valores contidos em a
media, sem = numpy.mean(a), scipy.stats.sem(a)
m = scipy.stats.t.ppf((1+conf)/2., len(a)-1)
h = m*sem
return media, h
A função anterior calcula a média e amplitude de um intervalo de confiança, a partir da lista de valores e do nível de confiança desejado.
O novo programa então ficaria:
import simpy
import random
import numpy # biblioteca numpy para computação científica http://www.numpy.org/
import scipy.stats # bilbioteca scipy.stats de funções estatísticas
def intervaloConfMedia(a, conf=0.95):
# retorna a média e a amplitude do intervalo de confiança dos valores contidos em a
media, sem = numpy.mean(a), scipy.stats.sem(a)
m = scipy.stats.t.ppf((1+conf)/2., len(a)-1)
return media, h
def geraChegada(env):
pesosList = [] # lista para armazenar os valores de pesos gerados
while True:
yield env.timeout(1)
# adiciona à lista o peso da entidade atual
pesosList.append(random.normalvariate(10, 5))
# calcula a amplitude do intervalo de confiança, com nível de significância = 95%
if len(pesosList) > 1:
media, amplitude = intervaloConfMedia(pesosList, 0.95)
print("%4.1f Média atual: %.2f kg\tAmplitude atual: %.2f kg"
%(env.now, media, amplitude))
# se a amplitude atende ao critério estabelecido, interronpe o processo
if amplitude < 0.5:
print("\n%4.1f Intervalo de confiança atingido após %s valores! [%.2f, %.2f]"
% (env.now, len(pesosList), media-amplitude, media+amplitude))
#termina o laço while
break
random.seed(100)
env = simpy.Environment()
chegadas = env.process(geraChegada(env))
env.run() # executa até o fim de todos os processos
O programa anterior leva 411 amostras para atingir o intervalo desejado:
...
410.0 Média atual: 10.17 kg Amplitude atual: 0.50 kg
411.0 Média atual: 10.18 kg Amplitude atual: 0.50 kg
411.0 Intervalo de confiança atingido após 411 valores! [9.68, 10.68]
Existem diversas maneiras de se estimar o intervalo de confiança utilizando-se as bibliotecas do Python. A maneira aqui proposta se baseia no numpy e no scipy.stats. Eventualmente tais bibliotecas não estejam instaladas na seu ambiente Python e eu antecipo: isso pode ser um baita problema para você =(
A questão aqui é que os modelos de simulação usualmente têm grande demanda por processamento estatístico de valores durante ou mesmo ao final da simulação. A biblioteca numpy facilita bastante esta tarefa, principalmente quando se considera o suporte dado pelos usuários do Stack Overflow .
Como sugestão, habitue-se a construir funções padronizadas para monitoramento e cálculos estatísticos, de modo que você pode reaproveitá-las em novos programas. Em algum momento, inclusive, você pode criar sua própria biblioteca de funções para análise de saída de modelos de simulação e compartilhar com a comunidade de software livre.