Lista 1 Processos Estocásticos

Prévia do material em texto

UFRN - ENGENHARIA DE COMPUTAÇÃO
TÓPICOS ESPECIAIS EM ENGENHARIA DE COMPUTAÇÃO
PROFESSOR: ADRIÃO DUARTE DORIA NETO
ALUNOS: VITOR RAMOS GOMES DA SILVA E ANDRESSA STÉFANY SILVA DE OLIVEIRA
LISTA 1
1) Considere o seguinte problema envolvendo uma sequência de experimentos. O problema consiste de duasurnas contendo bolas rotulados com os números 0 ou 1. A urna 1 contém 3 bolas sendo duas rotulados com 0 euma rotulada com 1, a urna 2 contém seis bolas sendo cinco rotulados com 1 e uma rotulada com 0. Oexperimento consiste em retirar uma bola de uma das urnas, verificar se é zero ou um e repor a bola e repetir oexperimento. A escolha da urna se faz por um experimento tipo cara ou coroa.
a) Gere um diagrama em treliça que reproduza a sequência de eventos
b) Calcule a probabilidade de se gerar a sequência 001110
Observe que as sequências geradas obedecem a seguinte propriedade ,
isto é, a probabilidade do evento presente de só do último evento passado. Esta propriedade caracteriza uma cadeia uma
Markov.
a) Diagrama em treliça: 0
1
0
1
...0,5
0,5
2/3
1/6
5/6
1/3
0,5*2/3
0,5*1/3
0,5*1/6
0,5*5/6
b)
2) Um canal de comunicações transmite símbolos ternários 0,1,2, gerados por uma determinada forte digital comas seguintes probabilidades P[0]=1/2, P[1]=1/4 e P[2]=1/4. A probabilidade de receber o símbolo transmitido é0,99 e probabilidade receber um símbolo diferente é 0,01.
a) Calcule a probabilidade de ter um erro na recepção
b) Supondo que foi recebido o símbolo 0, calcule a probabilidade de ter sido transmitido um 0, ou 1 ou 2.
P[ | , , . . . , , ] = P[ | ]sn sn−1 sn−2 s1 s0 sn sn−1
P(001110) = P(0) ∗ P(0|0) ∗ P(1|00) ∗ P(1|001). . .
= P(0) ∗ P(0) ∗ P(1) ∗ P(1) ∗ P(1) ∗ P(0)
= P(0 ∗ P(1)3 )3
= (0.5 ∗ 2/3 + 0.5 ∗ 1/6 ∗ (0.5 ∗ 1/3 + 0.5 ∗ 5/6 = 0.014358750750171465)3 )3
Lista1TopicosEspeciais http://localhost:8888/nbconvert/html/Lista1Topi...
1 de 15 29/08/2017 22:05
Árvora da fonte:
0,5
0,25
0,25
0,5
0,5
0,99
0,01
0,010
,99
0,99
0,01
0,5
0,5
0,5
0,5
00
01
02
10
12
11
22
20
21
a)
P(erro)= P(Transmitir X e receber Y)
b) P(Receber | transmitido)
3) Considere o problema de decisão caracterizado por uma sequência de eventos que podem ser apresentadospor um gráfico conhecido como árvore de decisão. O problema em questão consiste das escolhas e dasdecisões por parte de uma empresa de petróleo. Uma determinada empresa petrolífera obteve a concessão paraexplorar uma certa região. Os estudos anteriores (testes preliminares) estimam a probabilidade de existirpetróleo nessa região em 0,20. A companhia pode optar por um novo teste, que custa USS 500.000,00, sendoque, se realmente existe petróleo, esse teste dirá com uma probabilidade de 0,80 que existe, e se realmente nãoexiste, dirá com probabilidade 0,70 que não existe. Considerando que o custo de perfuração será de USS3.000.000,00 e que, se for encontrado petróleo, a companhia receberá USS 150.000.000,00. Considere, portantoos seguintes eventos e os seus complementos: (i) Evento T (a companhia faz o teste); (ii) Evento F (o teste éfavorável à existência de petróleo; (iii) Evento P (a companhia perfura o poço); (iv) Evento E (existepetróleo).
a) Construa a árvore indicando os nós de decisões e os nós ao acaso (probabilidades). Simbolize os dois de forma
diferente.
b) Determine o lucro =receita - despesas em cada percurso da árvore.
c) Usando o critério do melhor valor esperado, determine o valor esperado em cada nó de decisão.
d) Qual o valor esperado do lucro da companhia se forem tomadas as melhores decisões?
= 0.5 ∗ 0.01 + 0.25 ∗ 0.01 + 0.25 ∗ 0.01 = 0.01
P(0|0) = = 0.9949748743718594(0.5∗0.99)0.5∗0.99+0.25∗0.01∗1/2+0.25∗0.01∗1/2
P(0|1) = = 0.0025125628140703522(0.25∗0.01∗1/2)0.5∗0.99+0.25∗0.01∗1/2+0.25∗0.01∗1/2
P(0|2) = = 0.0025125628140703522(0.25∗0.01∗1/2)0.5∗0.99+0.25∗0.01∗1/2+0.25∗0.01∗1/2
Lista1TopicosEspeciais http://localhost:8888/nbconvert/html/Lista1Topi...
2 de 15 29/08/2017 22:05
a) e b) -500k-3M+150M
-500k
-500k-3M+150M
-500k
-500k-3M
-500k
-500k-3M
-500k
Exis
te
!Existe
Testa
r
Perfura
!Testar
Favor
ável
!Favorável
!Perfura
!Perfura
Perfura
Testar
!Testar
Perfura
!Perfura
!Perfura
Perfura
Perfura!Perfura
!Favorável
Favoráv
el
!Perfura
Perfura
0,8
0,2
0,3
0,7
0,8
0,2-500k
-500k
0
-3M
-3M+150M
0
Lista1TopicosEspeciais http://localhost:8888/nbconvert/html/Lista1Topi...
3 de 15 29/08/2017 22:05
In [1]: pTeste= 500*10**3 # depesas teste
pPerfura= 3*10**6 # despesas perfura
Lucro= 150*10**6 # lucro se tem
Pe= 0.2 # existe petroleo
Tf_e= 0.8 # se existe teste da favoravel 
Td_ne= 0.7 # se nao existe teste dar falso
pTeste= 500*10**3 # depesas teste
pPerfura= 3*10**6 # despesas perfura
Lucro= 150*10**6 # lucro se tem
ETFP= Pe*Tf_e*(-pTeste-pPerfura+Lucro)
ETF_P= Pe*Tf_e*(-pTeste)
ET_FP= Pe*(1-Tf_e)*(-pTeste-pPerfura+Lucro)
ET_F_P= Pe*(1-Tf_e)*(-pTeste)
E_TP= Pe*(-pPerfura+Lucro)
E_T_P=0
_ETFP= (1-Pe)*(Td_ne)*(-pTeste-pPerfura)
_ETF_P= (1-Pe)*(Td_ne)*(-pTeste)
_ET_FP= (1-Pe)*(1-Td_ne)*(-pTeste-pPerfura)
_ET_F_P= (1-Pe)*(1-Td_ne)*(-pTeste)
_E_TP= (1-Pe)*(-pPerfura)
_E_T_P= 0
print("c)")
print("Lucro esperado se não fizer o teste e perfurar ", E_TP+_E_T
P) #!TP
print("Lucro esperado se fizer o teste e perfurar ", ETFP+ET_F
P+_ETFP+_ET_FP) #TP
print("Lucro esperado se fizer o teste perfurando ou não ", ETFP+ETF_
P+ET_FP+ET_F_P+_ETFP+_ETF_P+_ET_FP+_ET_F_P) #T
print("Lucro esperado se teste favorável e perfurar ", ETFP+_ETF
P) #TFP
print("Lucro esperado se teste não favorável e perfurar ", ET_FP+_ET
_FP) #T_FP
print("Lucro esperado se teste não favorável e não perfurar ", ET_F_P+_E
T_F_P) #T_F_P
print("Lucro esperado se teste favorável e não perfurar ", ETF_P+_ET
F_P) #TF_P
print("Lucro esperado se não perfurar fazendo ou não o teste", ETF_P+ET_
F_P+E_T_P+_ETF_P+_ET_F_P+_E_T_P) #!P
print("Lucro esperado se fizer o teste e não perfurar ", ETF_P+ET_
F_P+_ETF_P+_ET_F_P) #T!P
print("Lucro esperado se não fizer o teste e não perfurar ", E_T_P+_E_
T_P) #!T!P
d) A melhor decisão é perfurar não fazendo o teste, sendo o lucro de 27 milhões de reais.
c)
Lucro esperado se não fizer o teste e perfurar 27000000.0
Lucro esperado se fizer o teste e perfurar 26500000.000000004
Lucro esperado se fizer o teste perfurando ou não 26000000.000000004
Lucro esperado se teste favorável e perfurar 21480000.000000004
Lucro esperado se teste não favorável e perfurar 5019999.999999999
Lucro esperado se teste não favorável e não perfurar -140000.00000000003
Lucro esperado se teste favorável e não perfurar -360000.0
Lucro esperado se não perfurar fazendo ou não o teste -500000.0
Lucro esperado se fizer o teste e não perfurar -500000.0
Lucro esperado se não fizer o teste e não perfurar 0
Lista1TopicosEspeciais http://localhost:8888/nbconvert/html/Lista1Topi...
4 de 15 29/08/2017 22:05
4) Considere nesta questão os conceitos básicos de confiabilidade, tempo médio de falhas, taxa de falhas ecálculo da confiabilidade de sistemas. Seja T o tempo de vida de um componente, define-se a confiabilidade
como sendo a função taxa média de falha é definida como sendo 
Considere também que um determinado sistema seja composto por subsistemas em serie ou em paralelo. Sobestas condições solucione as seguintes questões:
a) Demonstre que o tempo médio de vida (MTTF) definido por MTTF= sendo R(t)
função densidade de probabilidade da V.A. T.
b) Considerando que a função taxa de falhas seja constante e igual e que R(0)=1, mostre que e t , isto é,
T possui uma distribuição exponencial.
c) Para um sistema composto por n subsistemas independentes em série e que cada um tem taxa de falhas , calcule a
confiabilidade deste sistema.
d) Calcule a confiabilidade na questãoc para o caso em que os sistemas estejam em paralelo.
e) Para n=2 subsistemas e , calcule a razão entre a confiabilidade dos subsistemas em paralelo e
confiabilidade do subsistemas em série.
R(t) = P[T > t] = 1 − (t)FT r(t) = − (t)R′R(t)
E[T] = t (t)dt = R(t)dt∫ ∞0 ft ∫ ∞0
R(t) = 1 − (t)Ft
(t) = − (t)R′ ft
E[T] = t (t)dt = −t (t)dt∫ ∞0 ft ∫ ∞0 R′
u = −t, dv = (t)R′
du = −dt, v = R(t)
−t (t)dt = −tR(t) + R(t)dt = −(0 ∗ R(0) − tR(t)) + R(t)dt = R(t)dt∫ ∞0 R′ |∞0 ∫ ∞0 limt→∞ ∫ ∞0 ∫ ∞0
(t) = λfT e−λt
λ = − (t)R′R(t)
λR(t) + (t) = 0R′
λ (t) + (t) = 0R′ R″
− λf (t) − (t) = 0f ′
x + λ = 0
x = −λ
f (t) = ce−λt
c dt = 1 = − ( − 1)∫ ∞0 e−λt cλ e−λ∞
c = λ
f (t) = λe−λt
λi
(t) = ∏ (t) = ∏Rs Ri e− tλi
(t) = 1 − ∏ (1 − R(t)) = 1 − ∏ (1 − )Rp e− tλi
= = λλ1 λ2
= = = 2 − 1(t)Rp(t)Rs
1−(1−e−λt)2
(e−λt)2
1−1+2 −e−λt e−2λt
e−2λt eλt
Lista1TopicosEspeciais http://localhost:8888/nbconvert/html/Lista1Topi...
5 de 15 29/08/2017 22:05
5) Um certo tipo de chip de uma determinada marca tem uma duração de vida que segue uma distribuiçãoexponencial com vida média de 100 horas. Cada chip tem um custo de 10 reais e, se durar menos de 200 horas,existe um custo adicional de 8 reais.
a) Qual a probabilidade de um chip durar mais de 150 horas?
Então, a probabilidade de um chip durar mais de 150 horas é de 22,3%.
b) Foi proposta a compra de uma outra marca que tem tempo de vida média de 200 horas e um custo médio de 15 reais
e existe também a incidência do custo adicional de 8 reais se durar menos de 200 horas. Considerando o custo médio
decida se deve ser feita a troca de marcas.
Fabricante atual 
Novo Fabricante 
Logo, não trocar de fabricante é a melhor opção.
6) O número de chamadas telefônicas que chegam a uma central de comutação durante um intervalo de tempo é
uma variável aleatória com distribuição de Poisson. Considere um intervalo de dez minutos com .
Determine a probabilidade de mais de três chamadas chegam durante o intervalo de 10 minutos.
Assim, há a probabilidade de 14,29% de mais de três chamadas chegarem durante o intervalo de 10 minutos.
7) Em um determinado projeto de construção, são identificadas pela experiência dos construtores 16 atividadescríticas independentes. Os tempos de execução das atividades são variáveis aleatórias com valores médios evariâncias conhecidas mostrados na tabela abaixo. Considerando que o tempo do projeto possa ser dado pelasoma dos tempos de cada uma das atividades individuais independentes, calcule:
a)O tempo esperado para concluir o projeto
b)A probabilidade do tempo de conclusão do projeto de construção não exceda 50 semanas
c)O tempo do projeto que corresponda a uma probabilidade de 0,90 de concluir o projeto
Observação: Faça uso nas considerações para solução do problema o teorema do limite central.
f (t) = e−
t
100
100
F(t) = 1 − P(t > 150) = 1 − P(t <= 150) = 1 − F(150) = = 0.223e− t100 e− 150100
λ = 100
P(t <= 200) = F(200)
E(t) = F(200) ∗ (10 + 8) + (1 − F(200)) ∗ 10 = 16.9
λ = 200
P(t <= 200) = F(200)
E(t) = F(200) ∗ (15 + 8) + (1 − F(200)) ∗ 15 = 20.05
λ = 2
P(x) = e−λλxx!
P(x > 3) = 1 − P(x <= 3) = 1 − P(x = 0) − P(x = 1) − P(x = 2) − P(x = 3) = 1 − − −e−2200! e
−221
1!
e−
Lista1TopicosEspeciais http://localhost:8888/nbconvert/html/Lista1Topi...
6 de 15 29/08/2017 22:05
In [2]: import numpy as np
import scipy.stats
med= np.array([[2.7,3.2,4.6,2.1,3.6,5.2,7.1,1.5,3.1,4.2,3.6,0.5,2.1,1.5,
1.2,2.8],
[1.0,1.3,1.0,1.2,0.8,2.1,1.9,0.5,1.2,0.8,1.6,0.2,0.6,0.7,
0.4,0.7]])
mu=sum(med[0])
sx=np.math.sqrt(sum(med[1]))
print("a) Tempo esperado:", mu)
#print("Desvio padrao ", sx)
print("b) Probabilidade de terminar antes de 50 semanas:", scipy.stats.n
orm(mu,sx).cdf(50))
p90= scipy.stats.norm(mu,sx).ppf(0.9)
print("c) 90% de certeza terminar antes de", p90, "semanas")
8) Uma distribuição bidimensional contínua (X,Y) tem a função densidade de probabilidade
determine:
a) As distribuições marginais;
Distribuição em X: 
Distribuição em Y: 
f (x, y) = x ,e−x(y+1)0,
x > 0, y > 0
casocontrário
x dy = x dy∫ ∞0 e−x(y+1) ∫ ∞0 e−xe−xy
x = [0 − 1] =e−x 1−x e−xy|∞0 e−x e−x
x dx∫ ∞0 e−x(y+1)
u = x− > du = dx
dv = dx− > v =e−x(y+1) e−x(y+1)−(y+1)
uv − ∫ vdu = x −e−x(y+1)−(y+1) |∞0 ∫ ∞0 e
−x(y+1)
y+1
x −e−x(y+1)−(y+1) 1(y+1) 1x(y+1) e−x(y+1)
x − =e−x(y+1)−y+1 1(y+1)2 e−x(y+1)|∞0 1(y+1)2
a) Tempo esperado: 49.0
b) Probabilidade de terminar antes de 50 semanas: 0.598706325683
c) 90% de certeza terminar antes de 54.1262062622 semanas
Lista1TopicosEspeciais http://localhost:8888/nbconvert/html/Lista1Topi...
7 de 15 29/08/2017 22:05
b) O coeficiente de correlação ;
, onde 
c) A média condicional . 
=ρx,y COV(X,Y)σxσy
E(XY) = xy dxdy∫y ∫x e−x(y+1)
u = − > du = 2xdxx2
dv = e−tx t = y + 1− > v = −e−txt
uv − ∫ vdu = + ∫ 2xdxx2 e−txt e−txt
u = 2x− > du = 2x;
dv = − > v =txe−t −e
−tx
t2
− + ( + ∫ 2dx) = − − −x2 e−txt −2xe
−tx
t2
e−tx
t2 x2 e
−tx
t
2xe−tx
t2
2e−tx
t3 |∞0
dx =∫ ∞0 x2e−tx 2(y+1)3
xy dxdy = dy∫ ∞0 ∫ ∞0 e−x(y+1) ∫ ∞0 2y(y+1)3
= + +2y(y+1)3 Ay+1 B(y+1)2 C(y+1)3
C = −2
= =A(y+1 +B(y−1)−2)2(y+1)3
2y
(y+1)3
= A( + 2y + 1) + By + 1 − 2 = 2y− > A = 0y2
2A + B = 2− > B = 2
dy = − dy∫ ∞0 2y(y+1)3 ∫
∞
0
2
(y+1)2
2
(y+1)3
= − = − = 12y+1 1−1 2(y+1)2 1−2 |∞0
−2
y+1
1
(y+1)2 |
∞0
E(X) = x dx = −x + dx = −x − = 1∫ ∞0 e−x e−x ∫ ∞0 e−x e−x e−x|∞0
E(Y) = dy = − dy = ln(y + 1) + = indef inido∫ ∞0 y(y+1)2 ∫
∞
0
1
y+1
1
(y+1)2
1
y+1
= + − > A(y + 1) − 1 = Y; A = 1, B = −1.y(y+1)2 Ay+1 B(y+1)2
= = 1σx dx − E(X∫ ∞0 x2e−x )2‾ ‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾√
= = indef inidoσy dy − E(Y∫ ∞0 y
2
Y+1 )2‾ ‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾√
COV( (X, Y)) = = indef inidofx,y 1−E(Y)E( )−E(YY 2 )2√
E[X|Y] xf ( )dx∫ ∞−∞ xy
f ( ) = = = x(y + 1xy f (x,y)f (y) xe
−x(y+1)
1
(y+1)2
)2e−x(y+1)
E(X|Y) = (y + 1 dx =∫ ∞0 x2 )2e−x(y+1) 2y+1
Lista1TopicosEspeciais http://localhost:8888/nbconvert/html/Lista1Topi...
8 de 15 29/08/2017 22:05
x|y
x|y
Trabalho: Experimentos Computacionais
1) Utilizando os métodos de geração de números aleatórios, desenvolva um software que possibilite gerar asseguintes distribuições de probabilidades.
a) Exponencial
b) Normal (Gaussiana)
c) t-Student
d) Rayleight
e) Qui-quadrado
Apresente os resultados deste software sob forma de gráficos, tabelas e compare com os resultados analíticos.
In [19]: import numpy as np
import matplotlib.pyplot as plt
ft= lambda t,l: l*np.exp(-l*t)
Ft= lambda t,l: 1-np.exp(-l*t)
Fti= lambda t,l: -np.log(1-t)/l
# Utilizando o metodo da inversa
def expDist(l, Nsamples):
data= np.array(Fti(np.random.random(Nsamples),l))
return data
l=1/100
Nsamples= 99999
data= expDist(l, Nsamples)
plt.figure(figsize=(20,10))
plt.hist(data, bins=100,normed=True)
x1= np.linspace(data.min(),data.max(),1000)
y1= ft(x1,l)
plt.plot(x1,y1,'r')
plt.title('Expodencial')
plt.legend(['Teorica','Gerada'])
plt.show()
Lista1TopicosEspeciais http://localhost:8888/nbconvert/html/Lista1Topi...
9 de 15 29/08/2017 22:05
In [29]: import numpy as np
import matplotlib.pyplot as plt
fx= lambda x,mu,std: np.exp(-(x-mu)**2/(2*std**2))/np.sqrt(2*np.pi*std**
2)
#utilizando o metodo da rejeicao
def normalDist(mu, std, Nsamples):
data= np.array([])
while len(data) < Nsamples:
u1= np.random.random(Nsamples)*20-10
u2= np.random.random(Nsamples)*1
index= np.where(u2<=fx(u1,mu,std))[0]
data= np.append(data, u1[index])
#if(u2<=fx(u1,mu,std)):
# data= np.append(data, u1)
data.resize(Nsamples)
return data
mu= 0
std= 1
Nsamples= 99999
data= normalDist(mu, std, Nsamples)
plt.figure(figsize=(20,10))
plt.hist(data,normed=1,bins=100)
x= np.linspace(data.min(),data.max(),1000)
y= fx(x,mu,std)
plt.plot(x,y)
plt.title('Gaussiana')
plt.legend(['Teorica','Gerada'])
plt.show()
Lista1TopicosEspeciais http://localhost:8888/nbconvert/html/Lista1Topi...10 de 15 29/08/2017 22:05
In [30]: import numpy as np
import matplotlib.pyplot as plt
fe= lambda t,l: l*np.exp(-l*t)
flang= lambda x,k,l: ((l**k)*(x**(k-1))*np.exp(-l*x))/np.math.factorial(
k-1)
# somando varias distribuicoes expodenciais com o metodo da inversa
def erlangDist(k, l, Nsamples):
data= np.array([])
x= np.random.random((Nsamples,k))
data= np.array(-1/l*np.log(np.prod(x,axis=1)))
#for i in range(0,Nsamples):
# x=np.random.random(k)
# data= np.append(data, -1/l*np.log(np.prod(x)))
return data
Nsamples= 99999
K= 25
l= 1/15
data= erlangDist(K,l,Nsamples)
plt.figure(figsize=(20,10))
plt.hist(data, bins=100,normed=True)
x1= np.linspace(data.min(),data.max(),1000)
y1= flang(x1,K,l)
plt.plot(x1,y1,'r')
plt.title('Erlang')
plt.legend(['Teorica','Gerada'])
plt.show()
Lista1TopicosEspeciais http://localhost:8888/nbconvert/html/Lista1Topi...
11 de 15 29/08/2017 22:05
In [32]: import numpy as np
import matplotlib.pyplot as plt
import scipy.stats
from scipy.stats import chi2
Nsamples= 99999
df= 21
# ulizando a erlang com parametro lambda=1/2 e com graus de liberdade di
vidido por 2
# e somando a normal padrao ao quadrado no caso do numero de graus de li
berdade impa
def chi2Dist(df, Nsamples):
data= erlangDist(int(df/2),1/2,Nsamples)
if(df%2!=0):
data= data+normalDist(0,1,len(data))**2
#data=data+np.random.normal(0,1,len(data))**2
return data
data= chi2Dist(df, Nsamples)
 
plt.figure(figsize=(20,10))
plt.hist(data, normed=1, bins= 100)
x= np.linspace(data.max(), data.min(), 200)
plt.plot(x,scipy.stats.chi2.pdf(x, df))
plt.title('Chi2')
plt.legend(['Teorica', 'Gerada'])
plt.show()
Lista1TopicosEspeciais http://localhost:8888/nbconvert/html/Lista1Topi...
12 de 15 29/08/2017 22:05
In [35]: import numpy as np
import matplotlib.pyplot as plt
import scipy.stats
from scipy.stats import chi2
Nsamples= 99999
gl= 30
# Normal(0,1)/sqrt(Chi2(gl)/gl)
data= normalDist(0,1,Nsamples)/np.sqrt(chi2Dist(gl,Nsamples)/(gl))
#data= normalDist(0,1,Nsamples)/np.sqrt(np.random.chisquare(gl,Nsamples)
/gl)
plt.figure(figsize=(20,10))
plt.hist(data, normed=1, bins= 100)
x= np.linspace(data.max(), data.min(), 200)
plt.plot(x,scipy.stats.t.pdf(x,gl))
plt.title('T-student')
plt.legend(['Teorica', 'Gerada'])
plt.show()
Lista1TopicosEspeciais http://localhost:8888/nbconvert/html/Lista1Topi...
13 de 15 29/08/2017 22:05
In [10]: import numpy as np
import matplotlib.pyplot as plt
import scipy.stats
from scipy.stats import rayleigh
Nsamples= 99999
#Metodo da inversa
data= np.sqrt(-2*np.log(np.random.random(Nsamples)))
plt.figure(figsize=(20,10))
plt.hist(data, normed=1, bins= 100)
x= np.linspace(data.max(), data.min(), 200)
plt.plot(x,scipy.stats.rayleigh.pdf(x))
plt.title('Rayleigh')
plt.legend(['Teorica', 'Gerada'])
plt.show()
2) Utilizando método de Monte Carlo calcule as seguintes integrais definidas. Compare os resultados comvalores tabelados e ou soluções analíticas se existirem.
a) 
b) 
I = si (x) dx12π ∫
π
3
− π3 n
2 e−(x−1)2
I = xy dydx∫ 20 ∫ 20 e−( + )x2 y2
Lista1TopicosEspeciais http://localhost:8888/nbconvert/html/Lista1Topi...
14 de 15 29/08/2017 22:05
In [117]: import numpy as np
fx1= lambda x : np.sin(x)**2*np.exp(-(x-1)**2)
fx2= lambda x,y : x*y*np.exp(-(x**2+y**2))
def monteCarlo(Nsamples, a, b, f):
xi= (b-a)*np.random.random(Nsamples)+a
#xi= np.linspace(a,b,Nsamples)
I1= (b-a)*sum(fx1(xi))/len(xi)
return I1
def monteCarlo2(Nsamples, a1, b1, a2, b2, f):
xi1= (b1-a1)*np.random.random(Nsamples)+a1
xi2= (b2-a2)*np.random.random(Nsamples)+a2
#xi= np.linspace(a,b,Nsamples)
soma= 0
for i in range(0,len(xi1)):
soma= soma+f(xi1[i],xi2[i])
I2= (b1-a1)*(b2-a2)*soma/len(xi1)
return I2
I1= monteCarlo(99999, -np.pi/3, np.pi/3, fx1)
print("a)",I1)
I2= monteCarlo2(99999, 0, 2, 0, 2, fx2)
print("b)",I2)
a) 0.296733721198
b) 0.241494099042
Lista1TopicosEspeciais http://localhost:8888/nbconvert/html/Lista1Topi...
15 de 15 29/08/2017 22:05