4/10/13  

1  

       

Introdução  à  Simulação  Estocás5ca  usando  R  INF2035  –  PUC-­‐Rio,  2013.1  

Departamento  de  InformáAca  -­‐  PUC-­‐Rio    

Hélio  Lopes  Departamento  de  InformáAca  –  PUC-­‐Rio    

   

A  plataforma  R

•  R  é  uma   linguagem  de  programação  de  alto  nível   que  oferece  recursos  para   a  produção  de   gráficos  de  qualidade  bem  como  para   a   realização   de   cálculos   estaSsAcos   com   boa   precisão  numérica.    

•  Dentre  as  virtudes  de  R  não  podemos  deixar  de  mencionar  que  ela  obedece  os  padrões  de  soVware  livre  e  que  está  disponível  para  uma  boa  variedade  de  plataformas  de  hardware  e  sistemas  operacionais.    

•  Há   uma   grande   comunidade   de   desenvolvedores   contribuindo  permanentemente   para   a   criação   de   novas   funcionalidades,  bem   como   para   a   correção   dos   eventuais   problemas   que   R  apresenta.  

4/10/13  

2  

   

A  plataforma  R

•  A  comunidade  R_STAT,  disponível  em                                              h\p://br.groups.yahoo.com/group/R_STAT,              é  um  exemplo  disso  para  usuários  de  R  no  Brasil.    •  O  número  de  bibliotecas  de  qualidade  para  as  mais  diversas  aplicações  

cresce  constantemente  aumentando,  com   isso,  a  aplicabilidade  deste  sistema  em  diversas  áreas.    

•  R  pode  ser  obAdo  tanto  compilado  quanto  para  compilar  pelo  síAo                                                                        h\p://www.r-­‐  project.org.  •  Para   uma   excelente   introdução   ao   R   veja  “An   IntroducAon   to   R”   no  

próprio  síAo  onde  se  obtem  o  soVware.    

   

A  plataforma  R

•  Como   em   todas   as   outras   linguagens,   é   possível   declarar  variáveis  em  R  que  são  referenciadas  pelos  seus  nomes.  

•  Os  nomes  das  variáveis  podem  conter  qualquer  letra,  dígito  ou  ponto  mas  não  pode  começar  com  um  dígito.  

•  Os   nomes   das   variáveis   são     “case   sensi(ves”,   assim   dado,  Dado,  e  Dado  referem  a  três  variáveis  diferentes.  

•  O  operator  de  atribuição  em  R  é  “<-­‐”.  

4/10/13  

3  

   

A  plataforma  R

•  Exemplos:    >      Dado  <-­‐  4    >      Soma  <-­‐  2  +  2  >      um.nome.de.variavel.grande  <-­‐  sqrt(16)  

 Cada  linha  de  comando  acima  atribui  a  variável  que  está  do  lado  

esquerdo   do   operador   de   atribuição   o   valor   do   lado   direito.  Para  todos  esses  casos  é  atribuído  o  valor  4  às  variáveis  Dado,  Soma  e  um.nome.de.variavel.grande  

   

A  plataforma  R

•  Para   ver   o   valor   de   uma   variável   em   R,   você   deve   executar   o  comando  print,  ou  simplesmente  o  nome  da  própria  variável:  

>  print(Dado)  >  Dado  

4/10/13  

4  

   

A  plataforma  R

•  Nem  sempre  usaremos  variáveis  reais  no  R.  De  fato  na  maioria  das   aplicações   usaremos   variáveis   que   são   vetores,   que   são  uma  lista  de  objetos  do  mesmo  Apo.  

•  Existem   várias  maneiras   diferentes   de   se   criar   um   vetor   no   R.  Muitas   funções   e   operadores   retornam   uma   instância   de   um  vetor.  Aqui  veremos  algumas  formas  de  como  se  criar  um  vetor  e   algumas   funções   que   operam   esses   componente   a  componente.    

   

A  plataforma  R

•  A  função  c  “combina” ou  “coleciona” valores  em  um  vetor:  >  v  <-­‐  c(1,2,-­‐0.5,0.3)  >  v  

•  A  função  seq  cria  uma  sequencia  de  valores  em  um  vetor:  >  w  <-­‐  seq(1,4)  >  w  

•  Para  acessar  a  primeira  componente  do  vetor,  faça:  >  w[1]  

4/10/13  

5  

   

Distribuições  de  variáveis  aleatórias  

•  Neste  ponto  é  conveniente  lembrar  que  conhecer  a  distribuição  de   uma   variável   aleatória   consiste   em   ser   capaz   de   calcular   a  probabilidade  de  qualquer  evento  dentro  da  classe  de  eventos  possíveis,  suas  uniões  e  interseções.    

•  A  função  de  distribuição  acumulada,  definida  por  F(t)=  Pr(X<  t)    para   todo   t   real,   é   capaz   de   caracterizar   a   distribuição   de  qualquer  variável  aleatória  X.  

   

Distribuições  de  variáveis  aleatórias

•  As   propriedades   das   funções   de   distribuição   acumulada  são  as  seguintes:  

       

•  Reciprocamente,   toda   função   saAsfazendo   as   três  propriedades   acima   é   função   de   distribuição   acumulada,  isto   é,   ela   caracteriza   a   distribuição   de   alguma   variável  aleatória  real.  

Se xn ⇥ x entao F(xn) � F (x).Se x ⇥ �⇤ entao F (x)⇥ 0, e se x⇥⇤ entao F (x)⇥ 1Se x<y entaoF(x)� F (y).

4/10/13  

6  

   

Funções  densidade  de  probabilidade  

•  Caso  exista  uma  função                                                              para  a  qual                                                                                                                                                                                                                                                                então  dizemos  que  F  caracteriza  a  distribuição  de  uma  variável  aleatória  conSnua.  

•  Nesse  caso,  f  recebe  o  nome  de  densidade  de  probabilidade,  ou  simplesmente   de   densidade,   e   também   caracteriza   a  distribuição  da  viariável  aleatória.  

F(t) =� t�⇥ f(x) dx

f : � +

   

Funções  densidade  de  probabilidade

•  Toda  densidade  de  probabilidade  possui  duas  propriedades:  •  Ela  é  não  negaAva,  e  o  conjunto        recebe   o   nome   de   suporte   da   densidade   ou   da  distribuição,  e  

•  a  sua  integral  é  unitária,  isto  é:      

•  A   função   de   distribuição   acumulada   de   variáveis   aleatórias  conSnuas   é   estritamente   crescente   no   suporte   e,   portanto,  inversível  nesse  conjunto.  

{ x � : f(x) > 0}

�f(x) dx = 1.

4/10/13  

7  

   

Variáveis  aleatórias  discretas  

•  Caso   o   espaço   amostral   da   variável   aleatória   seja   finito   ou  enumerável,   dizemos   se   tratar   de   uma   variável   aleatória  discreta.  

•  Nesse   caso,   basta   conhecer   a   probabilidade   de   cada   evento  atômico  para  caracterizar  a  sua  distribuição,  isto  é,  a  função        para  todo  evento  elementar                .  

•  A  função                                                                                        recebe  o  nome  de  função  de  probabilidade,   e   é   possível   ver   que   a   função   de   distribuição  acumulada,  neste  caso,  é  constante  por  partes:  

(xi, pi = Pr(X = pi))xi

(xi, pi = Pr(X = pi))

F(t)=�

xi�t pi

   

Momentos  

•  Uma   informação   muito   importante   que   decorre   do  conhecimento  da  distribuição  de  uma  variável  aleatória  é  a  dos  seus  momentos.  

•  O   momento   de   ordem   k   da   variável   aleatória   X   com  distribuição   caracterizada   pela   densidade   f     é   dado   pela  integral      se  ela  exisAr.  

•  Caso  se  trate  de  uma  variável  aleatória  discreta,  a  integral  é  subsAtuída  por  uma  soma.  

E(Xk) =�

xkf(x) dx

4/10/13  

8  

   

Momentos  

•  O primeiro momento recebe o nome de esperança ou média, e a diferença entre o segundo momento e o quadrado da esperança é chamado variância, isto é:

Var(X) = E(X2)� (E(X))2.

   

Distribuição  uniforme  

A   distribuição   uniforme   no   intervalo   unitário   ocupa   um   lugar   de  destaque  no  universo  da  simulação  estocásAca.  Ela  é  definida  pela  densidade      onde                        é  a  função  indicadora  do  conjunto  A,  isto  é:        

f(x) = 1(0,1)(x)

1A

1A(x) = 0 se x ⇥� A

1A(x) = 1 se x � A

4/10/13  

9  

   

Distribuição  uniforme

A   probabilidade   de   qualquer   evento   dentro   do   intervalo   unitário  depende  apenas  do  comprimento  do  intervalo,  isto  é:      para  todos  0  <  a  <  b  <  1.    EXERCÍCIO:  Qual  é  a  função  distribuição  de  uma  variavel  uniforme  no  intervalo  (0,1)?    

Pr(a ⇥ X ⇥ b) = 1/(b� a)

   

Distribuição  exponencial

•  A   distribuição   da   variável   aleatória   X   chama-­‐se   exponencial  padrão  se  a  sua  densidade  é  

•  Com  isso,  a  sua  função  de  distribuição  acumulada  é    

•  A  distribuição   da   variável   aleatória  Y   =   λX,   com  X   exponencial  padrão   e   λ   >   0   chama-­‐se   exponencial   de   média   λ,   e   a   sua  densidade  é    

f(t)=exp{�t} +(t)

F(t)=(1-exp{�t}) +(t)

fY (t) =1�

exp{�t/�} +(t).

4/10/13  

10  

   

Distribuição  normal  

•  A  distribuição  da  variável  aleatória  X  é  chamada  normal  padrão  se  a  sua  densidade  é  

•  Não   dispomos   de   uma   forma   explícita   para   a   função   de  distribuição  acumulada  da  distribuição  normal  padrão,  mas  há  uma   relação   entre   ela   e   uma   importante   função   especial  chamada  “função  de  erro”:    

•  É  frequente  encontrar  a  notação                                                    ,  e  vale  que                                                                                                                                          ,  com  F  a  função  de  distribuição  acumulada  da  lei  normal    padrão.  

     

f(t) = 1�2�

exp�� 1

2 t2⇥

, para todo t ⇥ .

�(t) = 2��

� t0 exp{�x2} dx.

�(t) = erf(t)�(t) = 2F (

⇥2t)� 1

   

Distribuição  normal  

•  Se  X   segue  uma  distribuição  normal  padrão,  então  Y  =  σX  +  µ,  com  µ  um  número  real  qualquer  e  σ  >  0,  segue  uma  distribuição  normal   de  média   µ   e   desvio   padrão   σ   (ou,   equivalentemente,  variância  σ2  ).  Neste  caso,  a  densidade  de  Y  é    

fY (t) = 1�2�⇥

exp�� 1

2⇥2 (t� µ)2⇥

.

4/10/13  

11  

   

Distribuição  log-­‐normal  

•  Se  X  segue  uma  lei  normal  de  média  µ  e  variância  σ2  ,  então  Y  =  eX  segue  uma  distribuição  log-­‐normal  de  parâmetros  µ  e  σ2   .  A  densidade  que  caracteriza  esta  distribuição  é  

 •  Note  que  µ  e  σ2  são  apenas  os  parâmetros  que  indexam  esta  lei;  

não   são   a  média   e   a   variância   da   distribuição.   A  média   desta  distribuição  é  E(Y  )  =  exp  {  µ  +  σ2  /  2  }   ,  enquanto  a  variância  é  dada  por  

fY (y) = 1y⇥

2�⇥2 exp⇤� 1

2

�log(y)�µ

⇥

⇥2⌅.

Var(Y ) = exp�2µ + �2

⇥⇤exp

��2

⇥� 1

⌅.

   

Distribuição  Bernoulli  

•  Uma   das   distribuições   discretas   mais   simples   é   a   lei   de  Bernoulli.  

•  Uma   variável   aleatória   discreta   segue   a   lei   de   Bernoulli   com  parâmetro   p   que   pertence   ao   intervalo   [0,   1]   se   ela   adota   o  valor   1   com   probabilidade   p   ou   o   valor   0   com   probabilidade              1  −  p.  

•  Desta  forma,  a  função  de  probabilidade  que  caracteriza  esta  lei  é  {(0,  1  −  p),  (1,  p)}.  

4/10/13  

12  

   

Distribuição  Binomial  

•  A   soma   de   n   variáveis   aleatórias   independentes   e  idenAcamente   distribuídas,   cada   uma   obedecendo   uma  distribuição   de   Bernoulli   com   parâmetro   p   segue   uma   lei  Binomial  com  parâmetro  (n,  p).    

•  A  função  de  probabilidade  que  caracteriza  esta  distribuição  é  �k,

�nk

⇥pk(1� p)n�k

⇥para todo 0 ⇥ k ⇥ n.

   

Distribuição  de  Poisson  

•  A  distribuição  de  Poisson  caracteriza  a  contagem  do  número  de  eventos   que,   isoladamente,   são   raros   mas   para   os   quais   são  feitas  muitas  observações  independentes.    

•  Imaginemos,   por   exemplo,   a   situação   de   observar   em   um  pequeno   intervalo  de   tempo  a  chegada  de  um  freguês  em  um  banco,   ou   a   ocorrência   de   um   sinistro   para   uma   empresa   de  seguros   devido   a   um   roubo   de   carro;   a   probabilidade   disso  ocorrer   é   pequena,   mas   quando   a   observação   é   realizada   ao  longo  de  um  intervalo  maior  ela  se  torna  razoável.    

4/10/13  

13  

   

Distribuição  de  Poisson  

•  Nessas  condições,  e  com  algunas  restrições  adicionais,  podemos  modelar  o  número  de  eventos  com  a  variável  aleatória  X  cuja  lei  é  dada  por    onde   λ   >   0   é   a   média   e   o   parâmetro   que   caracteriza   a  distribuição.    

Pr(X = k) = �k

k! e��, para todo k � 0

   

Inversão  

•  Veremos   agora   como   transformar   variáveis   aleatórias  uniformes   em   variáveis   aleatórias   que   obedecem   outra  distribuição.    

•  O  resultado  que  será  apresentado  a  seguir  é  um  dos  pilares  da  simulação  estocásAca.    

4/10/13  

14  

   

Inversão  

Teorema:  Sejam  F  :  IR  →  [0,  1]  a  função  de  distribuição  acumulada  de  uma  variável  aleatória  e  F−  a  sua  inversa  generalizada,  dada  por  F−(t)  =  inf  {x  :  t  =  F(x)}.  Se  U  é  uma  variável  aleatória  que  segue  uma  lei  uniforme  no   intervalo   (0,  1)  então  F  é  a   função  de  distribuição  acumulada  da  variável  aleatória  resultante  da  transformação        V  =  F−(U).    Lema:  Quando  a  variável  aleatória  de  interesse  é  con[nua  vale  que                          F−(t)  =  F−1  (t)  para  todo  t  em  IR.  

   

Inversão  

•  A   aplicabilidade   deste   resultado   geral   restringe-­‐se   somente   à  disponibilidade  de  boas  implementações  de  F−1  ou  de  F−.  Casos  importantes  para  os  quais  não  há  formas  explícitas  simples  são  as  distribuições  normal  e  log-­‐normal,  por  exemplo.  

•  No  caso  da  distribuição  exponencial  com  parâmetro  λ,  a  função  de  distribuição  acumulada  é  dada  por   F(t)   =   (1   −   exp{−t/λ})1R+  (t),   logo  a  inversa  é  dada  por  F−1  (t)  =  −λ  log(1  −  t).  Com  isso,  a  variável   aleatória   Y   =   −λ   log(1   −   X)   segue   uma   distribuição  exponencial  com  média  λ  quando  X  segue  uma  lei  uniforme  no  intervalo  (0,  1).    

4/10/13  

15  

   

Inversão  

•  Para gerar ocorrências da distribuição de Bernoulli podemos usar o Teorema e fazer Y = 1(0,p) (U). Com isso, Y irá obedecer a lei desejada com parâmetro p.

•  As ocorrências de variáveis aleatórias binomiais podem ser construídas como a soma de eventos de variáveis aleatórias Bernoulli independentes.

•  Mas existe uma forma mais eficiente de gerar uma binomial…qual será?

   

Inversão  

•  A  geração  de  ocorrências   de   variáveis   aleatórias   Poisson  pode  parecer  inviável  com  este  método  pois,  em  princípio,  requereria  o  cálculo  de  infinitos  valores  de  probabilidade.    

•  Esse   problema   pode   ser   contornado   fazendo   uso   de   um  algoritmo  iteraAvo  que,  dada  a  ocorrência  da  variável  aleatória  uniforme      X(ω)  =  x,  calcule  o  valor  inicial  p0  =  Pr(Y  =  0)  =  e−λ  e  faça  a  comparação  com  apenas  o  primeiro  “degrau”  da  inversa  generalizada   da   função   de   distribução   acumulada.   Se   p0   <   x  então   retorne  0,   caso   contrário   calcule  p1  =  Pr(Y  =  1)   =   λe−λ   e  faça  a  comparação  com  p0  +  p1  ,  e  assim  por  diante…    

4/10/13  

16  

   

Método  polar  para  a  geração  de  uma  Normal  padrão  

•  De  volta  à  distribuição  normal,  podemos,  em  princípio,  uAlizar  o  seguinte  resultado:  

•  Sejam   U1   ,   U2   variáveis   aleatórias   independentes   e  idenAcamente   distribuídas   segundo   uma   lei   uniforme   no  intervalo  unitário.  Definamos  R  =  sqrt(−2  log  U1)  e  Θ  =  π(2  U2  −  1).  As  variáveis  aleatórias  X,  Y  definidas  por  X  =  R  cos  Θ  e  Y  =  R  sen   Θ   são   independentes   com   distribuição   comum   normal  padrão.    

   

Gerando  variáveis  aleatórias  no  R  

•  A   plataforma   R   oferece   diversas   roAnas   para   a   geração   de  ocorrências   de   variáveis   aleatórias   com   estas   e  muitas   outras  distribuições.    

•  Dentre   elas   podemos   mencionar   runif,   para   ocorrências   de  uniformes,  rbinom  para  binomiais   (e,  como  caso  parAcular,  de  Bernoulli),  rnorm  para  normais  e  rpois  para  Poisson.