GET00130 - Métodos Computacionais para Estatística II

Conteúdo da aula

Estruturas de Loopings e condicionamentos;
Criando funções;
Distribuições de probabilidade.

1 - Condicionamento

Em certas ocasiões, queremos ou precisamos executar parte de um código apenas se alguma determinada condição for atendida. O R fornece algumas opções básicas para estruturar seu código: if, else e if_else. A seguir, vamos discutir cada uma das três opções.

1.1 - O if

A estrutura básica do if é a seguinte:

if(condicao){
  # conjunto de comandos que serão executados
  # caso a condição seja TRUE
}

Notem que, o início do código se dá com o comando if seguido de parênteses e chaves;
Dentro do parênteses é preciso conter uma condição lógica (ou um conjunto de condições), que deverá ter como resultado ou TRUE ou FALSE;
Dentro das chaves temos um conjunto de comandos que serão executado se a condição do parênteses for TRUE.

Veja o exemplo abaixo:

# Criando o objeto x
x = 10

# só executa se a condição lógica for atendida
if(x < 20){
  print("x é menor do que 20")
}

[1] "x é menor do que 20"

Percebemos que a mensagem só é printada pois o objeto x é menor do que 20 (a condição foi satisfeita).

# Criando o objeto h
h = 8

# só executa se a condição lógica for atendida
if(x == h){
  x = x + h
}

# verificando o valor de x
x

[1] 10

Aqui percebemos que a condição lógica não foi satisfeita e o conjunto de códigos especificado dentro do if não foi executado.

1.2 - O if em conjunto com o else

Vamos ver agora como acrescentar a estrutura de condicionamento if o opcional else.

A estrutura básica do if com o opcional else else é a apresentada a seguir:

if(condicao){
  # conjunto de comandos que serão executados
  # caso a condição seja TRUE
} else{
  # conjunto de comandos que serão executados
  # caso a condição seja FALSE
}

Percebam que o início do código se dá do mesmo modo que no caso anterior, com o comando if seguido de parênteses e chaves;
Novamente, dentro do parênteses é preciso conter uma condição lógica (conjunto de condições), que deverá ter como resultado ou TRUE ou FALSE;
Dentro das chaves do if temos um conjunto de comandos que serão executado se a condição do parênteses for TRUE;
Logo em seguida temos a presença do opcional else, seguido de chaves;
Dentro das chaves do else temos um conjunto de comandos que serão executado se a condição do parênteses for FALSE.

# Criando o objeto z
z = 2

if(z < 20){
  z = z^h
  cat("O valor de z foi alterado para", z)
} else{
  z = z * h
  cat("O valor de z foi alterado para", z)
}

O valor de z foi alterado para 256

Podemos encadear diversos if e else em sequência como mostra o exemplo a seguir:

# Criando o objeto x
x = 10

if(x < 10){
  x = x + h
  cat("O valor de x foi alterado para", x)
} else if(x >=10 & x<20){
  x = x/h
  cat("O valor de x foi alterado para", x)
} else{
    x = x
  cat("O valor de x não foi alterado:", x)
}

O valor de x foi alterado para 1.25

Se x for um vetor ao invés de um número, será que o conjunto de comandos criados anteriormente executa as ações de interesse para todos os componentes de x?

# Criando o objeto x
x = c(310,20,230,40)

if(x < 10){
  y = 1
  cat("O número y criado é igual a", y)
} else if(x >=10 & x<20){
  y = 10
  cat("O número y criado é igual a", y)
} else{
    y = 100
  cat("O número y criado é igual a", y)
}

O número y criado é igual a 100

1.3 - A função if_else

Os comandos if com ou sem o opcional else não são vetorizados, como podemos perceber no exemplo criado na seção anterior. O conjunto de comandos executado, só levou em consideração o primeiro elemento do vetor x e exibiu uma mensagem de warning. Uma alternativa para estes casos, é utilizar a função if_else.

A função if_else tem a seguinte estrutura básica:

if_else(condicao, valor_se_TRUE, valor_se_FALSE)

Notem que, o primeiro argumento trata-se de uma condição (conjunto de condições), que retornará TRUE e FALSE para cada componente do vetor;
Já o segundo argumento trata-se do valor que será retornado quando o elemento do condicao for TRUE;
E o terceiro argumento trata-se do valor que será retornado quando o elemento do condicao for FALSE.

#Carregando o pacote dplyr para termos acesso a função if_else
library(dplyr)

# Criando o objeto x
x = c(310,20,230,40)

if_else(x >130, 1, x/10)

[1] 1 2 1 4

Podemos perceber pelo resultado da função if_else que foi executada em cada elemento do vetor x a operação compatível, isto é, se a condição for satisfeita, como resultado é retornado o valor 1 e se a condição não for satisfeita, como resultado é retornado o valor de x dividido por 10.

2 - Loopings

2.1 - O for

Um looping que utiliza o for no R tem a seguinte estrututa básica:

for(j in conjunto_de_valores){
  # conjunto de comandos que 
  # serão repetidos
  # tantas vezes quanto for o tamanho do conjunto de valores
}

Notem que, o início do loop se dá, de forma bem similar ao if, com o comando for seguido de parênteses e chaves;
Dentro do parênteses, temos um indicador que será usado durante o loop (no caso escolhemos o j para representar este indicador) e um conjunto de valores que será iterado (chamado aqui de conjunto_de_valores);
Dentro das chaves, temos o conjunto de comandos que será executado durante o loop.

#Definindo x como um objeto nulo que será preenchido futuramente
x = NULL

#Definindo a como o valor 2
a = 2

#Executando um conjunto de ações 10 vezes
for(j in 1:10){
  x[j] = a^j
}

#Visualizando x
x

 [1]    2    4    8   16   32   64  128  256  512 1024

#Apresentando uma mensagem com os 5 primiros valores de x
for(j in 1:5){
  cat("O valor do", j,"elemento de x é:", x[j], "\n")
}

O valor do 1 elemento de x é: 2 
O valor do 2 elemento de x é: 4 
O valor do 3 elemento de x é: 8 
O valor do 4 elemento de x é: 16 
O valor do 5 elemento de x é: 32

O “\(\n\)” serve para forçar uma quebra de linha entre as mensagens apresentadas em cada interação.

2.2 - O while

Um looping que utiliza o while no R tem a seguinte estrututa básica:

while(condicao_for_satisfeita){
  # conjunto de comandos que 
  # serão repetidos
  # enquanto a condicao continuar retornando TRUE
}

Notem que, o início do loop se dá com o comando while seguido de parênteses e chaves;
Dentro do parênteses, temos uma condição que será verificada a cada iteração;
Dentro das chaves, temos o conjunto de comandos que será executado enquanto a condição retornar TRUE.

#Definindo o objeto x
x = 1

#Somando a x o valor 1 enquanto sua soma não exceder 5
while( x < 5){
  x = x + 1
  print(x)
}

[1] 2
[1] 3
[1] 4
[1] 5

3 - Criando funções no R

Uma das grandes vantagens de softwares estatísticos como o R é o facilidade em criar funções.

Uma função, pode ser vista, como uma sequência de passos que você deseja que o programa execute. Se você faz uso, rotineiramente, de uma sequência de passos, ao invés de ter que reescrever os passos todas as vezes, você pode simplesmente criar uma função.

Quando fazemos o uso de um pacote, na verdade estamos usando funções que outras pessoas já criaram.

Ao criar uma função é muito importante definir suas entradas (argumentos) e saídas (resultados a serem retornados).

Para criar uma função no R, precisamos utilizar a função function que tem a seguinte estrutura básica:

nome.funcao = function(argumentos){
  # conjunto de comandos que 
  # serão realizados
}

Vamos supor que queiramos criar uma função que retorne a mensagem “Viciados no R!” toda vez que for evocada:

#Criando a função e chamando de mensagem

mensagem = function(){
  print("Viciados no R!")
}

#evocando a função
mensagem()

[1] "Viciados no R!"

A função que criamos acima é uma que não possuem argumentos (com a função getwd).

Suponha que queiramos criar uma função que retorna a soma de dois números.

#Criando uma função que soma dois números

soma.minha = function(x,y){
  valor = x + y
  return(valor)
}

soma.minha(2,4)

[1] 6

O que acontece se ao invés de somarmos dois números, somarmos um número com uma letra?

soma.minha(2,"h")

Teremos um erro na nossa função, mas esse erro não foi especidicado. Podemos melhorar a função usando os condicionamentos estudados anteriormente.

#Melhorando a função

soma.minha = function(x,y){
  if(class(x) != "numeric" | !is.numeric(y)){
    print("Você deve inserir dois números para que a função possa ser executada!")
  }else{
  valor = x + y
  return(valor)    
  }
}

soma.minha(2,"h")

[1] "Você deve inserir dois números para que a função possa ser executada!"

soma.minha("h",2)

[1] "Você deve inserir dois números para que a função possa ser executada!"

soma.minha(2,4)

[1] 6

class(x) != “numeric” e !is.numeric(y) são duas formas distintas de perguntar se um objeto é um número no R.

Podemos pré-definir alguns argumentos já na nossas suas funções.

#Criando uma função que calcula a potência de um número (a priori o número 10)
potencia = function(x,base = 10){
 valor = base^x
 return(valor)
}

#Calculando 10^3
potencia(3)

[1] 1000

#Calculando 5^3
potencia(3,5)

[1] 125

Atividade: Crie uma função que calcula o valor de \(f(z) = 3e^{-3z}, z > 0\). A função deve se chamar f.

Vamos plotar a função \(f(z)\) no intervalo de 0 a 3.

#Carregando o ggplot2
library(ggplot2)

#Plotando o gráfico
ggplot(data = tibble(z = c(0,3)), 
       mapping = aes(x = z)) +
  stat_function(fun = f)

Para plotarmos a função criada, usamos o geom stat_function para realizar o plot de uma função específica. Veja que usamos do artifício de definir um tibble no ggplot com uma única variável que assume dois valores, o mínimo e o máximo do eixo x (os valores que definem a região na qual foi pedido para que o plot fosse realizado).

4 - Distribuições de probabilidade

O R inclui a funcionalidade para operações com distribuições de probabilidades. Para cada distribuição há 4 operações básicas indicadas pelas letras:

d - calcula a densidade de probabilidade de \(f(x)\) no ponto.
p - calcula a função de probabilidade acumulada \(F(x)\) no ponto.
q - calcula o quantil correspondente a uma dada probabilidade.
r - sorteia uma amostra aleatória da distribuição.

Para melhor entendermos as quatro operações, vamos considerar que \(X \sim N(25,36)\), \(Y \sim Poisson(8)\) e \(Z \sim Gama(1,5)\).

Suponha que desejamos obter \(P(X < 28)\). Para calcularmos essa probabilidade usaremos a função pnorm.

Principais argumentos da função pnorm :

- q - o quantil da distribuição Normal que desejamos obter a área abaixo dele;

- mean - a média da distribuição;

- sd - o desvio padrão da distribuição.

#Calculando P(X < 28) em uma N(25,36)
pnorm(q = 28,
      mean = 25, 
      sd = sqrt(36))

[1] 0.6914625

Cuidado para não confundir, é preciso informar o desvio padrão na função pnorm (bem como nas demais funções p, q e r - norm).

Suponha que desejamos obter uma amostra aleatória de tamanho 10 de uma normal padrão.

Principais argumentos da função rnorm :

- n - o número de valores que desejamos gerar da distribuição Normal;

- mean - a média da distribuição;

- sd - o desvio padrão da distribuição.

#Obtendo uma amostra aleatória de tamanho 10 de uma N(0,1)
rnorm(n = 10)

 [1] -0.1495486  0.3178268 -0.3049138 -0.1747784  0.5224511  1.0113570
 [7]  0.5385454 -0.2150279 -1.0183621  1.1244323

Não informamos os parâmteros da normal em rnorm neste caso, pois são os valores pré-definidos na função.

E se desejarmos obter uma amostra aleatória de \(X\) de tamanho 10?

#Obtendo uma amostra aleatória de tamanho 10 de uma N(25,36)
rnorm(n = 10,
      mean = 25, 
      sd = 6)

 [1] 38.79426 20.27322 16.74726 31.61598 16.96716 15.56635 34.42529 26.55081
 [9] 20.35786 24.99401

Qual o valor do 3o decil de \(X\)?

Principais argumentos da função qnorm :

- p - a probabilidade abaixo do quantil da distribuição Normal que desejamos obter;

- mean - a média da distribuição;

- sd - o desvio padrão da distribuição.

#Obtendo o 3o decil de uma N(25,36)
qnorm(p = 0.3,
      mean = 25, 
      sd = 6)

[1] 21.8536

Qual o valor de \(f(10)\)?

Principais argumentos da função dnorm :

- x - o valor no qual deseja avaliar o valor da função densidade de probabilidade da distribuição Normal;

- mean - a média da distribuição;

- sd - o desvio padrão da distribuição.

#Calculando f(10) em uma N(25,36)
dnorm(x = 10,
      mean = 25, 
      sd = sqrt(36))

[1] 0.002921383

Qual a probabilidade de \(Y\) ser igual a 10?

Principais argumentos da função dpois :

- x - o valor no qual deseja avaliar o valor da função densidade de probabilidade da distribuição Poisson;

- lambda - a média da distribuição de Poisson.

#Calculando P(Y = 4) se Y ~ Poisson(8). 
#P(Y = 4) = f(4) em distribuições discretas
dpois(x = 10, 
      lambda = 8)

[1] 0.09926153

Qual a probabilidade de \(Y\) ser menor do que 4?

Principais argumentos da função ppois :

- q - o quantil da distribuição Poisson que desejamos obter a área abaixo dele;

- lambda - a média da distribuição Poisson.

#Calculando P(Y < 4) se Y ~ Poisson(8). 
ppois(q = 3,
      lambda = 8)

[1] 0.04238011

que é equivalente a…

#Somando P(X = 0) + P(X = 1) + P(X = 2) + P(X = 3) 
dpois(x = 0,
      lambda = 8) + 
  dpois(x = 1,
        lambda = 8) + 
  dpois(x = 2,
        lambda = 8) + 
  dpois(x = 3,
        lambda = 8)

[1] 0.04238011

Qual a probabilidade de \(Z\) ser menor do que 0,5?

Principais argumentos da função pgamma :

- q - o quantil da distribuição Gama que desejamos obter a área abaixo dele;

- shape - parâmetro de forma da gama;

- rate - parâmtero de escala da gama.

Neste caso, seja \(X \sim Gama(a,r)\), consideramos a forma da distribuição gama dada por \[f_X(x)= \frac{1}{s^a \Gamma(a)} x^{a-1} e^{-rx}, x>0\qquad \qquad a,b>0.\]

#Calculando P(Z > 0.5) se Z ~ Gama(1,5).
pgamma(q = 0.5,
       shape = 1, 
       rate = 5,
       lower.tail = FALSE)

[1] 0.082085

A modificação do argumento lower.tail para FALSE faz com que seja calculado o complementar da distribuição acumulada, isto é, é calculado \[P(Z > 0,5) = 1 - P(Z <= 0,5) = 1 - P(Z < 0,5) = 1 - F_Z(0,5).\] A segunda igualdade é justificada pois \(Z\) é uma variável contínua.

Atividade: Calcule o que se pede:

\(P(23 < X < 30)\);
\(P(3 < Y < 30)\);
Gere uma amostra de uma distribuição \(F\) com 2 e 4 graus de liberdade;
Qual o valor de \(c\), tal que \(P(Z > c) = 0,7\);

5.1 - Plotando funções das distribuições de probabilidade

Plotando a função densidade de probabilidade de uma distribuição normal com média 10 e variância 25. Usaremos do ggplot2 o stat_function.

Principais argumentos da função stat_function :

- fun - a função a ser plotada;

- args - os valores dos argmentos que a função definifa em fun depende (tem que ser definidos em uma lista).

#Plotando o gráfico da densidade de uma N(10,25)
ggplot(data = tibble(z = c(-10,30)), 
       mapping = aes(x = z)) +
  stat_function(fun = dnorm,
                args = list(mean = 10,
                            sd = 5)) +
  labs(y = "Densidade")

Plotando a função de distribuição acumulada de uma distribuição normal com média 10 e variância 25.

#Plotando o gráfico da função de distribuição acumulada de uma N(10,25)
ggplot(data = tibble(z = c(-10,30)), 
       mapping = aes(x = z)) +
  stat_function(fun = pnorm,
                args = list(mean = 10,
                            sd = 5)) +
  labs(y = "Densidade")

Gerando uma amostra aleatória de uma \(N(10,25)\). Em seguida, suponha que não conhecemos s distribuição e gostaríamos de checar se a amostra tem um formato semelhante com a distribuição \(N(10,25)\).

#Gerando uma amostra de tamanho 300 de uma N(10,25)
base = tibble(amostra = rnorm(300,10,5))

#Plotando o histograma da amostra com base na frequênica
#Controlando as amplitudes dos intervalos de classe
hist1 = base |> 
  ggplot(mapping = aes(x = amostra)) +
  geom_histogram(binwidth = 3,
                 fill = "Blue") +
    labs(y = "Frequência")

#Plotando o histograma da amostra com base na frequênica
#Controlando o número de intervalos de classe
hist2 = base |> 
  ggplot(mapping = aes(x = amostra)) +
  geom_histogram(bins = 7,
                 fill = "Blue") +
    labs(y = "Frequência")

#Plotando os gráficos em uma mesma figura
library(gridExtra)
grid.arrange(hist1,hist2,ncol = 2)

O argumento binwidth defini o tamanho que você gostaria que os intervalos do histograma possuíssem. Você pode ao invés de definir o tamanho do intervalo, definir o número de intervalos por meio do argumento bins. Podemos perceber as diferenças entre o número de intervalos de classes que cada histograma possui acima, o primeiro possui 11 e o segundo possui 7.

Para compararmos o histograma com a densidade de uma distribuição, precisamos que o mesmo não esteja definido com base nas frequências. Para realizar essa mudança, iremos acrescentar o y no aes no geom_histogram.

#Plotando o histograma da amostra com base na frequênica
base |> 
  ggplot(mapping = aes(x = amostra)) +
  geom_histogram(aes(y = ..density..),
                 binwidth = 3) +
  stat_function(fun = dnorm,
                args = list(mean = 10,
                            sd = 5),
                colour = "Red") +
  labs(y = "Densidade")

Se quisermos plotar uma função de probabilidade ao invés de uma função densidade de probabilidade. Basta fazermos as adequações no stat_function.

#Plotando o gráfico da função probabilidade de uma Poisson(1)
ggplot(data = tibble(x = 0:10), 
       mapping = aes(x = x)) +
    stat_function(geom = "point", 
                  n = 11, 
                  fun = dpois, 
                  args=list(lambda = 1), 
                  colour = "Red") +
  labs(x = "z", y = "P(Z = z)")

Os argumentos geom e n permitem que seja feito um gráfico de pontos na quantidade de pontos especificado (11 pois de 0 a 10 existem 11 números inteiros).

Outras distribuições:

dexp - densidade da exponencial;
rchisq - sorteia uma amostra da distribuição qui-quadrado;
qt - calcula o quantil de uma t-student;
pf - calcula \(F(x)\) de uma \(F\);
rbeta - sorteia uma amostra da distribuição beta;
dunif - densidade de uma uniforme (contínua);
E a uniforme discreta???