Referência: S. Ross. Simulation, 2nd edition. New York: Academic Press, 1997.

1 Introdução

Neste material estamos interessados em estudar métodos para resolver o seguinte problema:

Problema. Seja \(X\) uma variável aleatória (v.a.) contínua cuja função de distribuição acumulada é dada por \[ F_{X}(x)=\Pr(X\leq x)=\int_{-\infty}^{x}f_{X}(t)dt \] onde \(f_{X}\) é a função densidade de probabilidade de \(X\), que é dada por \[ f_{X}(x)=\frac{dF_{X}(x)}{dx}. \] Desejamos então simular valores de \(X\). Ou seja, desejamos instruir um computador para gerar automaticamente valores de \(X\) de acordo com a “lei” \(F_{X}\) (ou \(f_{X}\)).

A notação descrita acima é geral e será utilizada ao longo de toda exposição. Três métodos serão estudados:

2 Simulação de números aleatórios

Tal como discutido no caso discreto, o bloco construtivo básico para simulação de v.a.’s contínuas é a função runif() que gera números aleatórios no intervalo \((0,1)\). Tais números são realizações de v.a.’s uniformemente distribuídas no intervalo \((0,1)\).

3 Método da transformação inversa

O método da transformação inversa já foi discutido para o problema de simular valores de uma v.a. discreta. O método é geral é pode ser aplicado sem modificações para o problema de simular valores de uma v.a. contínua. O algoritmo é como segue: 

Algoritmo: método da transformação inversa
Entrada: invGF() (inversa generalizada de F)
  1 Gere um número aleatório u 
  2 Calcule x <- invGF(u)
Saída: x

Resultado. Seja \(U\sim U(0,1)\). Para qualquer função de distribuição \(F\), a v.a. \(X\) definida por \(X=F^{-1}(U)\) possui distribuição \(F\), onde \(F^{-1}\) é a inversa generalizada de \(F\).

É simples verificar o resultado acima. Observe que \[ \Pr(X\leq x)=\Pr(F^{-1}(U)\leq x)=\Pr(U\leq F(x))=F(x). \] Logo \(X\sim F\).

3.1 Simulação da distribuição exponencial

Uma v.a. contínua \(X\) possui distribuição exponencial com taxa \(\lambda>0\) se sua função densidade de probabilidade é dada por \[ f_{X}(x)=\left\{ \begin{array}{ll} 0 & x < 0 \\ \lambda e^{-\lambda x} & x \geq 0. \end{array} \right. \] É simples verificar que sua função de distribuição acumulada é dada por \[ F_{X}(x)=\Pr(X\leq x)=\left\{ \begin{array}{ll} 0 & x < 0 \\ 1-e^{-\lambda x} & x \geq 0. \end{array} \right. \] Além disso, segue que \[ E(X)=\frac{1}{\lambda}\hspace{1cm}V(X)=\frac{1}{\lambda^{2}}. \] Usamos a notação \(X\sim\mbox{Exp}(\lambda)\) para informar que \(X\) é uma v.a. com distribuição exponencial com taxa \(\lambda>0\).

Exemplo. Simule um valor da distribuição \(\mbox{Exp}(1)\). Neste caso, a função de distribuição acumulada é dada por \[ F(x)=\left\{ \begin{array}{ll} 0 & x < 0 \\ 1-e^{-x} & x \geq 0 \end{array} \right. \] e a inversa generalizada é dada por \[ F^{-1}(u)=-\ln(1-u)\hspace{1cm}0<u<1. \] Assim, para simular um valor da distribuição \(\mbox{Exp}(1)\), basta fazer \[ X=-\ln(1-U) \] onde \(U\sim U(0,1)\) é um número aleatório. Esta equação pode ser simplificada sabendo que \(U\stackrel{D}{=}1-U\), onde \(\stackrel{D}{=}\) significa igual em distribuição. Sendo assim, temos que para simular um valor da distribuição \(\mbox{Exp}(1)\), basta fazer \[ X=-\ln(U). \] No R, para simular um valor da distribuição \(\mbox{Exp}(1)\) basta fazer 

u <- runif(1)
x <- -log(u)
x
## [1] 0.8133468

Exemplo. Utilize um loop se desejar simular \(n\) valores da distribuição \(\mbox{Exp}(1)\). 

n <- 100
x <- rep(NA,n)
for(i in 1:n){
  x[i] <- -log(runif(1))
}
hist(x,freq=FALSE,xlab="x",ylab="densidade",main="")
points(seq(0,5,0.05),dexp(seq(0,5,0.05),rate=1),type="l",col="blue")

Alternativamente, basta fazer 

n <- 100
x <- -log(runif(n))
hist(x,freq=FALSE,xlab="x",ylab="densidade",main="")
points(seq(0,5,0.05),dexp(seq(0,5,0.05),rate=1),type="l",col="blue")

Exemplo. Simule um valor da distribuição \(\mbox{Exp}(\lambda)\), onde \(\lambda>0\). Sabemos que se \(X\sim\mbox{Exp}(1)\), então \[ Y=\frac{1}{\lambda}X\sim\mbox{Exp}(\lambda). \] Assim, para simular um valor da distribuição \(\mbox{Exp}(\lambda)\), basta fazer \[ Y=-\frac{1}{\lambda}\ln(U) \] onde \(U\sim U(0,1)\) é um número aleatório.

3.2 Simulação da distribuição de Erlang

Uma v.a. contínua \(X\) possui distribuição de Erlang com parâmetros \(k\) (parâmetro de forma) e \(\lambda\) (parâmetro de taxa) se sua função densidade de probabilidade é dada por \[ f_{X}(x)=\left\{ \begin{array}{ll} 0 & x < 0 \\ \frac{\lambda^{k}}{k!}x^{k-1}e^{-\lambda x} & x \geq 0 \end{array} \right. \] onde \(k\) é um inteiro positivo e \(\lambda>0\). Além disso, temos que \[ E(X)=\frac{k}{\lambda}\hspace{1cm}V(X)=\frac{k}{\lambda^{2}}. \] Usamos a notação \(X\sim\mbox{Erl}(\lambda)\) para informar que \(X\) é uma v.a. com distribuição de Erlang com parâmetros \(k\) e \(\lambda\). Note que a distribuição de Erlang generaliza a distribuição exponencial e a distribuição gama generaliza a distribuição de Erlang, visto que: \[ \mbox{Exp}(\lambda)\stackrel{D}{=}\mbox{Erl}(k=1,\lambda) \] e \[ \mbox{Erl}(k,\lambda)\stackrel{D}{=}\mbox{Ga}(\alpha=k,\beta=\lambda). \] Por fim, note também que \[ Y=\sum_{i=1}^{k}X_{i}\sim\mbox{Erl}(k,\lambda) \] se \(X_{1},\ldots,X_{k}\stackrel{idd}{\sim}\mbox{Exp}(\lambda)\).

Exemplo. Simule um valor da distribuição \(\mbox{Erl}(k,\lambda)\). Sabendo que \(\mbox{Erl}(k,\lambda)\) é igual em distribuição à soma de \(k\) exponenciais independentes, todas com taxa \(\lambda>0\), então basta fazer \[ Y=\sum_{i=1}^{k}-\frac{1}{\lambda}\ln(U_{i})=-\frac{1}{\lambda}\ln\left(\prod_{i=1}^{k}U_{i}\right) \] onde \(U_{1},\ldots,U_{k}\stackrel{idd}{\sim}U(0,1)\) são \(k\) números aleatórios independentes. Vejamos o seguinte exemplo considerando a distribuição \(\mbox{Ga}(3,2.3)\): 

n <- 300;k <- 3;lamb <- 2.3
y <- rep(NA,n)
for(i in 1:n){
  y[i] <- -(1/lamb)*log(prod(runif(k)))
}
hist(y,freq=FALSE,xlim=c(0,6),ylim=c(0,0.8),xlab="x",ylab="densidade",main="")
points(seq(0,6,0.05),dgamma(seq(0,6,0.05),shape=3,rate=2.3),type="l",col="blue")

4 Método da aceitação e rejeição

O método da aceitação e rejeição já foi discutido para o problema de simular valores de uma v.a. discreta. O método é geral é pode ser aplicado sem modificações para o problema de simular valores de uma v.a. contínua. O algoritmo é como segue: 

Algoritmo: método da aceitação e rejeição
Entrada: f, g (distribuições) e c (constante apropriada)
  1 Gere um número aleatório y da distribuição g
  2 Gere um número aleatório u 
  3 Se u < f(y)/cg(y), então faça x = y e pare. Caso contrário, volte para o passo 1.
Saída: x ~ f

Resultado. O método da aceitação e rejeição simula uma v.a. \(X\) tal que \(X\sim f\). Além disso, o número de iterações do algorítmo que são necessárias para gerar um valor simulado de \(X\) é uma v.a. com distribuição geométrica com média \(c\).

Seja \(A\) o evento “o valor \(y\) gerado no passo \(1\) é aceito”. Segue da regra de Bayes que: \[ f_{X}(y)=f_{Y|A}(y|A)=\frac{g(y)\Pr(A|Y=y)}{\Pr(A)}. \] Como \[ \Pr(A|Y=y)=\Pr\left(U\leq\frac{f(y)}{cg(y)}\right)=\frac{f(y)}{cg(y)} \] segue do teorema da probabilidade total que \(\Pr(A)\) é dada por \[ \Pr(A)=\int_{\mathbb{R}}\Pr(A|Y=y)g(y)dy=\int_{\mathbb{R}}\frac{f(y)}{cg(y)}g(y)dy=\frac{1}{c}. \] Logo, \[ f_{X}(y)=f_{Y|A}(y|A)=\frac{g(y)[f(y)/cg(y)]}{1/c}=f(y). \] e, portanto, \(X\sim f\). Além disso, a cada iteração do algoritmo, o valor gerado no passo \(1\) pode ser aceito ou não (ou seja, o evento \(A\) ocorre ou não a cada iteração). As iterações são independentes uma da outra, visto que \(Y\) (no passo \(1\)) é independente \(U\) (no passo \(2\)). Se \(N\) denota o número de iterações do algorítmo que são necessárias para gerar um valor simulado da variável \(X\), segue que \(N\sim\mbox{Geo}(p)\), com \(p=\Pr(A)=1/c\), ou seja, \[ \Pr(N=n)=\left(1-\frac{1}{c}\right)^{n-1}\frac{1}{c},\hspace{1cm}n=1,2,3,\ldots \] e \(E(N)=1/\Pr(A)=c\). Concluímos que o número de iterações que são necessárias para gerar um valor simulado da distribuição de \(X\) é uma v.a. com distribuição geométrica com média \(c\). Assim, quanto menor for o valor de \(c\) para que \(c\cdot g\) envelope \(f\), mais eficiente será a simulação.

4.1 Simulação da distribuição gama

Como já citado, a distribuição gama generaliza a distribuição de Erlang. Uma v.a. contínua \(X\) possui distribuição gama com parâmetros \(\alpha\) (parâmetro de forma) e \(\beta\) (parâmetro de taxa) se sua função densidade de probabilidade é dada por \[ f_{X}(x)=\left\{ \begin{array}{ll} 0 & x < 0 \\ \frac{\beta^{\alpha}}{\Gamma(\alpha)}x^{\alpha-1}e^{-\beta x} & x \geq 0 \end{array} \right. \] onde \(\alpha,\beta>0\). Além disso, temos que \[ E(X)=\frac{\alpha}{\beta}\hspace{1cm}V(X)=\frac{\alpha}{\beta^{2}}. \] Usamos a notação \(X\sim\mbox{Ga}(\alpha,\beta)\) para informar que \(X\) é uma v.a. com distribuição gama com parâmetros \(\alpha\) e \(\beta\). Note também que \(\mbox{Exp}(1)\stackrel{D}{=}\mbox{Ga}(1,1)\). Por fim, o processo de simulação da distribuição gama será discuto em três etapas:

  • Simulação da \(\mbox{Ga}(\alpha,1)\), quando \(\alpha>1\), usando o método da aceitação e rejeição.
  • Simulação da \(\mbox{Ga}(\alpha,1)\), quando \(0<\alpha<1\), usando o fato que \(X=YU^{1/\alpha}\sim\mbox{Ga}(\alpha,1)\) se \(Y\sim\mbox{Ga}(\alpha+1,1)\) e \(U\sim U(0,1)\) são independentes.
  • Simulação da \(\mbox{Ga}(\alpha,\beta)\), quando \(\alpha,\beta>0\), usando o fato que \(X=\beta^{-1}Y\sim\mbox{Ga}(\alpha,\beta)\) se \(Y\sim\mbox{Ga}(\alpha,1)\).

Resta claro que basta discutir a primeira etapa, pois as outras duas etapas decorrem imediatamente da primeira. Sendo assim, para simular valores da \(\mbox{Ga}(\alpha,1)\), quando \(\alpha>1\), basta usar o método da aceitação e rejeição, considerando como distribuição auxiliar \(g\) a distribuição exponencial com média \(\alpha\) (ou seja, no caso a \(\mbox{Exp}(\lambda=\alpha^{-1})\)), visto que \(\alpha\) é também a média da \(\mbox{Ga}(\alpha,1)\). É possível mostra que a exponencial com média \(\alpha\) é, dentro da classe de modelos exponenciais, a melhor escolha para a distribuição auxiliar, visto que esta distribuição minimiza o valor da constante \(c\) que permite que \(c\cdot\mbox{Exp}(\alpha^{-1})\) envelope a distribuição \(\mbox{Ga}(\alpha,1)\). Neste caso, também é possível mostrar que esta constante \(c\) é dada por \[ c=\frac{\alpha^{\alpha}}{\Gamma(\alpha)e^{\alpha-1}}. \]

Portanto, um algoritmo para simular da \(\mbox{Ga}(\alpha,1)\), quando \(\alpha>1\), é como segue: 

Algoritmo: método da aceitação e rejeição para simular da Ga(alpha,1)
  1 Gere dois números aleatórios, u1 e u2
  2 Gere y <- -alpha*log(u1) # y é gerado da exponencial com média alpha
  3 Calcule c em função do valor alpha
  4 Se u2 < Ga(y;alpha,1)/c*Exp(y;alpha^(-1)), então x <- y. Caso contrário, volte para 1
Saída: x

Vejamos agora como implementar no R este método para simular \(n>1\) valores da distribuição \(\mbox{Ga}(\alpha,1)\) quando, por exemplo, \(\alpha=3/2\). 

d.exp <- function(x,par.lambda){
  if(x>=0){out <- par.lambda*exp(-par.lambda*x)}else{out <- 0}
  return(out)
}
d.ga <- function(x,par.alpha,par.beta){
  constante <- ((par.beta^(par.alpha))/(gamma(par.alpha)))
  if(x>=0){out <- constante*x^(par.alpha-1)*exp(-par.beta*x)}else{out <- 0}
  return(out)
}
na.ga <- function(par.alpha){
  flag <- 0
  repeat{
   u1 <- runif(1); u2 <- runif(1)
   y <- -par.alpha*log(u1)
   c <- (par.alpha^(par.alpha))/(gamma(par.alpha)*exp(par.alpha-1))
   par.lambda <- 1/par.alpha
   aux <- (d.ga(y,par.alpha,1))/(c*d.exp(y,par.lambda))
   if(u2 < aux){out <- y;flag <- 1}
   if(flag==1){break}
  }
  return(out)
}
na.ga(par.alpha=1.5) # simula um único valor da Ga(1.5,1)
## [1] 0.6355191
#
# Para simular n>1 valores, use um loop
#
n <- 10000
amostra <- rep(NA,n)
for(k in 1:n){amostra[k] <- na.ga(par.alpha=1.5)}
hist(amostra,freq=FALSE,xlim=c(0,10),ylim=c(0,0.6),xlab="x",ylab="densidade",main="")
points(seq(0,10,0.05),dgamma(seq(0,10,0.05),shape=1.5,rate=1.0),type="l",col="blue")

4.2 Simulação da distribuição beta

Uma v.a. contínua \(X\) possui distribuição beta com parâmetros \(\alpha,\beta>0\) se sua função densidade de probabilidade é dada por \[ f_{X}(x)=\frac{1}{B(\alpha,\beta)}x^{\alpha-1}(1-x)^{\beta-1}\hspace{1cm}0\leq x\leq 1. \] Além disso \[ E(X)=\frac{\alpha}{\alpha+\beta}\hspace{1cm}V(X)=\frac{\alpha\beta}{(\alpha+\beta)^{2}(\alpha+\beta+1)}. \]

Fica como exercício, mostra que se \(X\sim\mbox{Ga}(\alpha,1)\) e \(Y\sim\mbox{Ga}(\beta,1)\), então \[ \frac{X}{X+Y}\sim\mbox{Beta}(\alpha,\beta). \] Este resultado pode ser usado para simular valores da distribuição beta.

5 Simulação da distribuição normal univariada

Aqui usaremos o método polar para simular valores da distribuição normal padrão. Este método é baseado na transformação Box-Muller. Antes de entrar em detalhes a respeito deste método, vamos lembrar alguns fatos importantes.

Uma v.a. contínua \(X\) possui distribuição normal padrão se sua função densidade de probabilidade é dada por \[ f_{X}(x)=\frac{1}{\sqrt{2\pi}}e^{x^{2}/2}\hspace{1cm}-\infty<x<\infty. \] É simples verificar que \[ E(X)=0\hspace{1cm}V(X)=1. \] Além disso, se \(-\infty<\mu<\infty\), \(\sigma>0\), e \(Y=\sigma X + \mu\), então \(Y\) possui distribuição normal com média \(\mu=E(Y)\) e variância \(\sigma^{2}=V(Y)\), e sua função densidade de probabilidade é dada por \[ f_{Y}(y)=f_{X}\left(x=\frac{y-\mu}{\sigma}\right)\left|\frac{dx}{dy}\right|= \frac{1}{\sqrt{2\pi}\sigma} e^{(y-\mu)^{2}/2\sigma^{2}}\hspace{1cm}-\infty<y<\infty. \] Usamos a notação \(X\sim N(\mu,\sigma^{2})\) para informar que \(X\) é uma v.a. com distribuição normal com média \(\mu\) e variância \(\sigma^{2}\). A normal padrão é denotada por \(N(0,1)\).

Outro ponto importante que precisamos lembrar é o método do Jacobiano para encontrar a densidade de um vetor aleatório \(\mathbf{Y}=(Y_{1},\ldots,Y_{d})\) contínuo que é função de um outro vetor aleatório contínuo \(\mathbf{X}=(X_{1},\ldots,X_{d})\). Considere que \(\mathbf{X}\sim f_{\mathbf{X}}\) e seja \[ \begin{array}{ll} X_{1} & = & X_{1}(Y_{1},\ldots,Y_{d}) \\ \vdots & & \vdots \\ X_{d} & = & X_{d}(Y_{1},\ldots,Y_{d}) \end{array} \] a representação componente a componente da transformação inversa da transformação \(\mathbf{Y}=\mathbf{Y}(\mathbf{X})\), que por sua vez é representada componente a componente tal como segue: \[ \begin{array}{ccc} Y_{1} & = & Y_{1}(X_{1},\ldots,X_{d}) \\ \vdots & & \vdots \\ Y_{d} & = & Y_{d}(X_{1},\ldots,X_{d}). \end{array} \] Então, a densidade de \(\mathbf{Y}\) é definida por \[ f_{\mathbf{Y}}(y_{1},\dots,y_{d})=f_{\mathbf{X}}(x_{1}(y_{1},\dots,y_{d}),\ldots,x_{d}(y_{1},\dots,y_{d})) \cdot\mbox{det}[J(y_{1},\dots,y_{d})] \] para todo \(\mathbf{y}=(y_{1},\dots,y_{d})^{\prime}\) no alcance da transformação \(\mathbf{Y}=\mathbf{Y}(\mathbf{X})\), onde \(J\) é a matrix Jacobiana da transformação inversa \(\mathbf{X}=\mathbf{X}(\mathbf{Y})\). A matrix \(J\) é definida por: \[ J(y_{1},\dots,y_{d})=\left(\frac{\partial x_{i}}{\partial y_{j}}\right)_{d\times d}= \left( \begin{array}{ccc} \frac{\partial x_{1}}{\partial y_{1}} & \cdots & \frac{\partial x_{1}}{\partial y_{d}} \\ \vdots & \ddots & \vdots \\ \frac{\partial x_{d}}{\partial y_{1}} & \cdots & \frac{\partial x_{d}}{\partial y_{d}} \end{array} \right). \] Agora estamos em posição de discutir o método polar para simular valores da normal padrão.

Como já informado, o método polar é baseado na transformação Box-Muller, que é definida por: \[ \begin{array}{rcl} X & = & X(W,\Theta)=\sqrt{W}\cos(\Theta) \\ Y & = & Y(W,\Theta)=\sqrt{W}\sin(\Theta) \end{array} \] onde \(W\sim\mbox{Exp}(1/2)\) e \(\Theta\sim U(0,2\pi)\) são v.a.’s independentes. Consequentemente, a função densidade de probabilidade conjunta de \(W\) e \(\Theta\) é dada por \[ f_{W,\Theta}(w,\theta)=\frac{1}{2\pi}\cdot\frac{1}{2}e^{-w/2},\hspace{1cm} 0<w<\infty\hspace{1cm}0<\theta<2\pi. \] Lembrando que a transformação inversa da transformação \((W,\Theta)\mapsto (X,Y)\) é dada por \[ \begin{array}{rcl} W & = & W(X,Y)=X^{2}+Y^{2} \\ \Theta & = & \Theta(X,Y)=\arctan(Y/X) \end{array} \] e usando o método do Jacobiano, segue que \[ f_{X,Y}(x,y)=f_{W,\Theta}(w=x^{2}+y^{2},\theta=\arctan(y/x))\cdot\mbox{det}[J(x,y)] \] para todo \(-\infty<x<+\infty\) e \(-\infty<y<+\infty\), onde \[ J(x,y)= \left( \begin{array}{cc} \frac{\partial w}{\partial x} & \frac{\partial w}{\partial y} \\ \frac{\partial \theta}{\partial x} & \frac{\partial \theta}{\partial y} \end{array} \right)= J(x,y)= \left( \begin{array}{cc} 2x & 2y \\ -\frac{y}{x^{2}+y^{2}} & \frac{x}{x^{2}+y^{2}} \end{array} \right). \] Assim, como \(\mbox{det}[J(x,y)]=2\), segue que \[ f_{X,Y}(x,y)=\frac{1}{2\pi}\cdot\frac{1}{2}e^{-(x^{2}+y^{2})/2}\cdot 2=\frac{1}{2\pi}e^{-(x^{2}+y^{2})/2} \] de onde se conclui imediatamente que \(X\sim N(0,1)\) e \(Y\sim N(0,1)\) são v.a.’s independentes.

Portanto, para simular um par de valores da distribuição normal padrão, basta aplicar a transformação Box-Muller \[ \begin{array}{rcl} X & = & \sqrt{W}\cos(\Theta)=\sqrt{-2\ln(U_{1})}\cos(2\pi\cdot U_{2}) \\ Y & = & \sqrt{W}\sin(\Theta)=\sqrt{-2\ln(U_{1})}\sin(2\pi\cdot U_{2}) \end{array} \] onde \(U_{1}\) e \(U_{2}\) são números aleatórios no intervalo \((0,1)\). Lembre-se que \(W\sim\mbox{Exp}(1/2)\) e \(\Theta\sim U(0,2\pi)\) são simulados como segue: \[ W=-2\ln(U_{1})\hspace{1cm}\Theta=2\pi\cdot U_{2}. \]

O algoritmo descrevendo o método polar é dado por: 

Algoritmo: método polar para simular um par de valores da N(0,1) 
  1 Gere números aleatórios, u1 e u2
  2 Calcule w <- -2*log(u1)
  3 Calcule theta <- 2*pi*u2
  4 Calcule x <- sqrt(w)*cos(theta) 
  5 Calcule y <- sqrt(w)*sin(theta)
Saída: x,y

Vejamos como implementar no R o método polar para simular \(n>1\) valores da distribuição \(N(0,1)\). 

transf.bm <- function(u){
  w <- -2*log(u[1])
  theta <- 2*pi*u[2]
  out <- sqrt(w)*c(cos(theta),sin(theta))
  return(out)
}
n <- 10000; m <- n/2
amostra <- rep(NA,n)
for(i in 1:m){
  j <- 2*i
  amostra[(j-1):j] <- transf.bm(runif(2))
}
1/sqrt(2*pi)
## [1] 0.3989423
hist(amostra,freq=FALSE,xlim=c(-5,5),ylim=c(0,0.45),xlab="x",ylab="densidade",main="")
points(seq(-5,5,0.05),dnorm(seq(-5,5,0.05),mean=0,sd=1),type="l",col="blue")

Finalmente, se desejar simular valores de uma v.a. \(Y\sim N(\mu,\sigma^{2})\), basta usar a transformação \(Y=\sigma X + \mu\), onde \(X\sim N(0,1)\). Vejamos agora como implementar no R o método polar para simular \(n>1\) valores da distribuição \(N(\mu,\sigma^{2})\) quando, por exemplo, \(\mu=10\) e \(\sigma=\sqrt{2}\). 

transf.bm <- function(u){
  w <- -2*log(u[1])
  theta <- 2*pi*u[2]
  out <- sqrt(w)*c(cos(theta),sin(theta))
  return(out)
}
na.normalp <- function(n,media,desvio){
  aux <- n%%2
  if(aux==1){nn <- n + 1}else{nn <- n}
  m <- nn/2; out <- rep(NA,nn)
  for(i in 1:m){j <- 2*i; out[(j-1):j] <- transf.bm(runif(2))}
  out <- desvio*out[1:n] + media
  return(out)
}
n <- 9999
m <- 10; dp <- sqrt(2); v <- dp^(2)
amostra <- na.normalp(n,media=m,desvio=dp)
n==length(amostra)
## [1] TRUE
aux1 <- m - 5*dp;aux2 <- m + 5*dp
aux1;aux2;1/sqrt(2*pi*v)
## [1] 2.928932
## [1] 17.07107
## [1] 0.2820948
hist(amostra,freq=FALSE,xlim=c(aux1,aux2),ylim=c(0,0.35),xlab="x",ylab="densidade",main="")
points(seq(aux1,aux2,0.05),dnorm(seq(aux1,aux2,0.05),mean=m,sd=dp),type="l",col="blue")

Uma versão mais eficiente do método polar, que evita cálculos de funções trigonométricas, é como segue: 

Algoritmo: método polar alternativo para simular um par de valores da N(0,1) 
  1 Gere dois números aleatórios, u1 e u2
  2 Calcule v1 <- 2*u1 - 1; v2 <- 2*u2 - 1 e s <- v1^(2) + v^(2) 
  3 Se S <= 1, então calcule x <- sqrt(-2*log(s)/s)*v1 e y <- sqrt(-2*log(s)/s)*v2. Senão, volte para 1
Saída: x,y