Probabilidade e Estatística

Projeto computacional

O ficheiro LE vs GDP.csv contém um conjunto de dados que incluem as variáveis:
- Life expectancy – esperança de vida à nascença (em anos) e
- GDP per capita – produto interno bruto per capita (em dólares dos E. U. A. de 2011),
para um grande número de países no período de 1950 a 2018 (adaptado de Kaggle).

Após a leitura desse ficheiro no R, selecione os dados referentes ao ano de 1989 e a todos os países dos seguintes continentes: Europe, Oceania e North America (em Inglês no ficheiro).

Com recurso ao pacote ggplot2, produza um único diagrama de dispersão que:
1. represente os valores da esperança de vida à nascença em função do produto interno bruto per capita;
2. represente os valores do produto interno bruto per capita numa escala logarítmica de base 10;
3. permita distinguir o continente a que cada país pertence.
Tenha em conta que o texto no ficheiro de dados se encontra em Inglês e, por simplicidade, mantenha todo o texto do gráfico nessa língua.

Submeta um ficheiro em formato PDF com uma única página A4, que inclua:
1. O código em R, que deve incluir os comandos para leitura e seleção dos dados do ficheiro.
2. O gráfico produzido.
Resolução
```
library(ggplot2)
theme_set(theme_gray())

# url <- 'https://web.tecnico.ulisboa.pt/~paulo.soares/pe/projeto/LE vs GDP.csv'
# dados <- read.csv(url, check.names = FALSE)
# Por exemplo:
# ano <- 1966
# continentes <- c("Africa", "Europe", "Asia")

dados |>
  subset(Year == ano & Continent %in% continentes) |>
  ggplot() +
  geom_point(aes(`GDP per capita`, `Life expectancy`, color = Continent)) +
  scale_x_log10() +
  labs(title = paste("Life expectancy vs GDP per capita in", ano),
       x = "GDP per capita (2011 US dollars)",
       y = "Life expectancy (years)")
```
O ficheiro countries of the world.xlsx contém um conjunto de dados com informação diversa sobre 227 países (adaptado de Kaggle).

Após a leitura desse ficheiro no R, use o pacote ggplot2 para construir um histograma de frequências relativas acumuladas (com um polígono das mesmas frequências sobreposto) dos valores da variável Service (em Inglês no ficheiro).

Na construção do gráfico:
1. fixe o número de classes ou intervalos em que os dados são agrupados recorrendo à regra de Sturges;
2. ignore o que representa a variável e a sua unidade;
3. tenha em conta que a variável não toma valores negativos;
4. dado que o texto no ficheiro de dados se encontra em Inglês, mantenha, por simplicidade, todo o texto do gráfico nessa língua.
Submeta um ficheiro em formato PDF com uma única página A4, que inclua:
1. O código em R, que deve incluir os comandos para leitura dos dados do ficheiro.
2. O gráfico produzido.
Resolução
```
library(ggplot2)
theme_set(theme_light())

# dados <- readxl::read_xlsx("countries of the world.xlsx")
# Por exemplo:
# variavel <- "Literacy (%)"

nbins <- nclass.Sturges(dados[[variavel]])
# Alternativa:
# nbins <- ceiling(log(nrow(dados), base = 2) + 1)

ggplot(dados, aes(x = .data[[variavel]],
                  y = cumsum(after_stat(count)) / nrow(dados))) + 
  geom_histogram(bins = nbins, fill = "cornflowerblue",
                 color = "cornflowerblue", boundary = 0) +
  geom_freqpoly(bins = nbins, color = "magenta", boundary = 0) +
  labs(title = paste("Distribution of", variavel, "in 227 countries"),
       y = "Cumulative relative frequency")
```

O ficheiro Turbine.csv contém um conjunto de valores simulados da potência gerada em cada hora por uma turbina elétrica sob diferentes condições meteorológicas para todos os dias de um ano (adaptado de Kaggle).

Após a leitura desse ficheiro no R, selecione os dados referentes ao mês de novembro e calcule a média e o desvio padrão diários da variável Air temperature | (ºC).

Com recurso ao pacote ggplot2, produza um único gráfico que permita ilustrar:

a evolução das médias diárias ao longo do mês;
a variabilidade diária da variável, representando, em cada dia, o intervalo média $\pm$ 2 $\times$ desvio padrão.

Tenha em conta que o texto no ficheiro de dados se encontra em Inglês e, por simplicidade, mantenha todo o texto do gráfico nessa língua.

Submeta um ficheiro em formato PDF com uma única página A4, que inclua:

O código em R, que deve incluir os comandos para leitura e seleção dos dados do ficheiro.
O gráfico produzido.

Resolução

library(ggplot2)
theme_set(theme_dark())

# url <- 'https://web.tecnico.ulisboa.pt/~paulo.soares/pe/projeto/Turbine.csv'
# aero <- read.csv(url, check.names = FALSE)
# Por exemplo:
# mes <- "August"
# variavel <- "Air temperature | (ºC)"
# k <- 1.5

dados <- subset(aero, substr(aero$Day, 1, 3) == substr(mes, 1, 3))
dados$Day <- as.numeric(substr(dados$Day, 5, 6))

ggplot(dados, aes(Day, .data[[!!variavel]])) +
  stat_summary(fun.data = mean_sdl, fun.args = list(mult = k), geom = "pointrange", color = "lightblue", linewidth = 1) +
  stat_summary(fun.y = mean, geom = "line", color = "magenta", linewidth = 1) +
  labs(y = variavel, title = paste("Daily mean values and variability in", mes))

# Alternativa

# mes.values <- subset(aero, substr(aero$Day, 1, 3) == substr(mes, 1, 3))
# mes.mean <- aggregate(mes.values[, variavel], by = list(Day = mes.values$Day), mean)
# mes.sd <- aggregate(mes.values[, variavel], by = list(Day = mes.values$Day), sd)
# dados <- merge(mes.mean, mes.sd, by = c("Day"))
# names(dados) <- c("Day", "Mean", "Sd")
# 
# dados$Day <- as.numeric(substr(dados$Day, 5, 6))
# 
# ggplot(dados, aes(Day, Mean)) +
#   geom_pointrange(aes(ymin = Mean - k * Sd, ymax = Mean + k * Sd),
#                   color = "lightblue", linewidth = 1) +
#   geom_line(color = "magenta", linewidth = 1) +
#   labs(y = variavel, title = paste("Daily mean values and variability in", mes))

Suponha que qualquer transgressão grave do código da estrada é registada pelas autoridades de trânsito com uma probabilidade $p=0.23$ e que as ocorrências de transgressões distintas são acontecimentos independentes. Qualquer condutor que atinja o registo de 3 trangressões graves é penalizado, iniciando-se um processo que pode levar à inibição temporária de conduzir.

Considere que $X$ é a variável aleatória que representa o número de transgressões graves cometidas por um condutor até ser penalizado pela primeira vez. Fixando a semente em 1622 e baseando-se apenas na geração de números aleatórios de distribuições de Bernoulli, produza uma amostra simulada de $n=1300$ valores de $X$ e indique o valor observado da mediana amostral com uma casa decimal.

Nota: não use a função sample nas resolução desta questão.

Resolução
```
# Por exemplo:
# p <- 0.25
# n <- 1500
# seed <- 1234

set.seed(seed)

x <- vector(length = n)

for(i in 1:n) {
  n_registos <- 0
  n_transg <- 0
  while(n_registos < 3) {
    n_transg <- n_transg + 1
    n_registos <- n_registos + rbinom(1, 1, p)
  }
  x[i] <- n_transg
}

# Resposta
solution <- median(x)
```
O valor observado da mediana amostral é 12.0.
A função de densidade de probabilidade conjunta da velocidade do vento ($U$, em metros por segundo) e da direção do vento ($V$, em graus a partir de uma direção de referência), durante determinado período de tempo num local onde se pretende construir um arranha-céus, é dada por \[f_{U,V} (u,v) = \frac{\pi b - b \theta v + 2 \theta \, u v}{2 \pi^2 b^2},\] para $u \in [0, b]$ e $v \in ]-\pi, \pi]$, onde $\theta = 0.95$ e $b = 10$.

Tendo em conta que $E[V]=0$ tem-se, neste caso, que $Cov[U,V]=E[UV]$. Para calcular e comparar valores aproximados da covariância, defina convenientemente no R a função integranda do integral duplo $E[UV]$ e implemente os seguintes dois procedimentos:

A. Integração numérica
1. Instale o pacote pracma.
2. Use a função integral2 (sem alterar qualquer argumento opcional) para obter um primeiro valor aproximado do integral duplo $E[UV]$.
B. Integração de Monte Carlo
1. Fixando a semente em 2028:
  1. gere $n=10000$ números pseudoaleatórios $(u)$ da distribuição uniforme no intervalo $]0, 10[$;
  2. gere $n=10000$ números pseudoaleatórios $(v)$ da distribuição uniforme no intervalo $]-\pi, \pi[$.
2. Avalie a função integranda nos $n=10000$ pontos aleatórios $(u,v)$ e calcule a média desses valores.
3. Multiplique a média anterior pela área da região de integração, obtendo assim uma aproximação de Monte Carlo do integral duplo $E[UV]$.
Indique o desvio absoluto entre os valores obtidos em A. e B., arredondado a 6 casas decimais.

Resolução
```
# Por exemplo:
# seed <- 1234
# theta <- 0.55
# b <- 7

f_integranda <- function(u, v) {
 val <- u * v * (pi * b - b * theta * v + 2 * theta * u * v) / (2 * (pi * b)^2)
 return(val)
}

# A.
library(pracma)

aproxNum <- integral2(f_integranda, 0, b, -pi, pi)$Q

# B.
set.seed(seed)

n <- 10000
u_vals <- runif(n, min = 0, max = b)
v_vals <- runif(n, min = -pi, max = pi)

f_vals <- f_integranda(u_vals, v_vals)
media <- mean(f_vals)

aproxMC <- media * (b * 2 * pi)

# Resposta
solution <- abs(aproxNum - aproxMC)
```
O desvio absoluto entre os dois valores aproximados de $Cov[U,V]$ é igual a 0.045606.
Considere que $X_1, \ldots, X_{20}$ são variáveis aleatórias independentes e identicamente distribuídas de acordo com uma distribuição uniforme no intervalo $]-1, 1[$.

Com o objetivo de calcular um valor aproximado de $P(Y > 3)$, em que $Y=\sum_{i=1}^{20}{X_i}$, pretende-se comparar os seguintes dois métodos:
1. Simulação
  
  Fixando a semente em 1983:
  1. gere 10000 amostras de dimensão $n=20$ da distribuição uniforme no intervalo $]-1, 1[$;
  2. calcule um valor simulado de $Y$ para cada uma das amostras geradas;
  3. calcule a proporção de valores simulados de $Y$ que são maiores do que 3.
2. Recurso ao Teorema do limite central
Usando apenas funções do R básico, implemente e aplique ambos os métodos. Indique o valor absoluto da diferença entre os dois resultados obtidos multiplicado por 100 com 4 casas decimais.

Resolução
```
# Por exemplo:
# seed <- 3185
# n <- 30
# l <- 3
# a <- 2

set.seed(seed)

aceites <- 0
for (i in 1:10000) {
  aceites <- aceites + (sum(runif(n, min = -a, max = a)) > l)
}
sim <- aceites / 10000

tlc <- 1 - pnorm(sqrt(3 / n) * l / a)

solution <- abs(sim - tlc) * 100
```
O valor pedido é igual a 0.1761.
Considere as médias $\bar{X}_1$ e $\bar{X}_2$, referentes a duas amostras aleatórias independentes de dimensões $n_1=116$ e $n_2=247$ recolhidas de duas populações normais com o mesmo valor esperado $\mu$ mas variâncias diferentes $v_1=4$ e $v_2=6$.

Pretende-se estimar $\mu$ através de uma combinação linear entre as médias das duas amostras aleatórias, \[\hat\mu=\alpha\bar{X}_1+(1-\alpha)\bar{X}_2, \alpha\in]0,1[.\]

Use a função optimize do R para determinar o valor de $\alpha$ que minimiza a variância de $\hat\mu$. Apresente o resultado com 4 casas decimais.

Nota: use uma tolerância (tol) igual a $10^{-6}$ e não altere qualquer outro argumento opcional da referida função.

Resolução
```
# Por exemplo:
# n1 <- 100
# n2 <- 280
# v1 <- 2
# v2 <- 6

func <- function(alpha) {
  val <- alpha^2 * v1 / n1 + (1 - alpha)^2 * v2 / n2
  return(val)
}

solution <- optimize(func, interval = c(0, 1), tol = 10^-6)$minimum
```
O valor de $\alpha$ que minimiza a variância de $\hat{\mu}$ é 0.4133.
Assuma que, em ensaios de determinado tipo de viga de cimento pré-esforçado, a deflexão (em mm) até à fissuração da viga é uma variável aleatória $X$ com distribuição normal de valor esperado desconhecido $\mu$ e desvio padrão $\sigma= 1.2$.
1. Determine a dimensão mínima da amostra, $n_{min}$, que garante que o intervalo de confiança a $\gamma\times 100\% = 95\%$ para $\mu$ possua amplitude inferior ou igual a $0.6$.
2. Fixando a semente em 1498, seleccione sem reposição $n_{min}$ observações da seguinte amostra casual:
  
  39.8, 40.5, 39.9, 38.4, 40.3, 41.1, 40.0, 43.0, 39.4, 40.7, 40.8, 39.0, 41.4, 38.9, 37.0, 40.1, 40.2, 40.4, 43.4, 41.0, 40.4, 40.9, 41.0, 39.4, 39.1, 38.3, 38.3, 40.6, 42.4, 39.2, 39.5, 39.1, 42.7, 40.8, 39.3, 39.9, 41.4, 39.1, 40.4, 40.1, 39.7, 39.4, 41.2, 42.4, 42.3, 40.5, 39.1, 43.4, 41.2, 39.7, 41.4, 38.5, 40.2, 39.4, 39.4, 38.6, 39.6, 40.0, 38.0, 37.7, 40.1, 39.5, 39.3, 40.5, 39.3, 42.2, 40.2, 41.8, 41.9, 41.3, 41.1, 43.2, 39.2, 38.1, 41.7, 41.1, 39.2, 40.2, 39.4, 39.4, 40.4, 38.3, 38.1, 43.2, 40.6, 43.0, 38.4
3. Obtenha o intervalo de confiança a $\gamma\times 100\% = 95\%$ para $\mu$ com base na amostra selecionada em 2.
4. Indique o rácio, arredondado a 4 casas decimais, entre a amplitude do intervalo de confiança calculado em 3. e a média da amostra selecionada em 2.
Resolução
```
# Por exemplo:
# dados <- c(37.9, 40.2, 40.3, 41.9, 40.6, 41.5, 40.4, 39.8, 40.1, 42.1, 39.9, 39.9,
# 39.2, 39.6, 39.8, 39.3, 40.2, 41.0, 40.2, 39.9, 40.2, 42.3, 40.3, 39.6, 41.0, 39.8, 
# 41.3, 40.9, 42.4, 40.9, 38.2, 40.9, 41.6, 38.2, 41.0, 39.8, 38.3, 42.1, 39.9, 39.7,
# 40.9, 38.3, 40.3, 41.6, 40.1, 40.4, 40.7, 42.4, 40.1, 40.1, 41.4, 39.2, 42.0, 42.1,
# 40.2, 40.5, 41.7, 40.2, 39.7, 41.6, 40.4, 40.3, 39.9, 39.3, 41.4, 42.0, 40.0, 41.4,
# 41.0, 39.4, 39.3, 38.6, 40.6, 42.2, 39.2, 37.5)
# seed <- 1234
# sigma <- 1.1
# gama <- 0.96
# amp <- 0.7

quantil <- qnorm((1 + gama) / 2)
nmin <- ceiling((2 * quantil * sigma / amp)^2)

set.seed(seed)
dadosmin <- sample(dados, nmin)

amp_ic <- 2 * quantil * sigma / sqrt(nmin)

solution <- amp_ic / mean(dadosmin)
```
A amplitude relativa do intervalo de confiança a $\gamma \times 100\% = 95\%$ para $\mu$ baseado nas $62$ observações selecionadas é 0.0149.
Seja $\left(X_1, \ldots, X_{15}\right)$ uma amostra aleatória de uma população com distribuição exponencial de parâmetro $\lambda$. Para testar a hipótese $H_0: \lambda=2$ contra $H_1: \lambda > 2$ pode usar-se a variável fulcral \[Q = 2 \lambda \sum_{i=1}^{15}X_i\] cuja distribuição é conhecida.

Considere um teste de hipóteses que conduz à rejeição de $H_0$ quando $Q_0<16.31$, em que $Q_0$ é a estatística de teste obtida da variável fulcral anterior sob a validade de $H_0$.

Fixando a semente em $3082$, gere $9000$ amostras de dimensão $15$ da distribuição exponencial de parâmetro $\lambda=2$ e calcule a frequência relativa de rejeições de $H_0$. Com base no resultado obtido selecione o único valor que completa corretamente a seguinte afirmação:

O teste de hipóteses aplicado garante que a probabilidade da hipótese $H_0$ não ser rejeitada caso $H_0$ seja verdadeira é aproximadamente igual a:
1. 0.0854
2. 0.9754
3. 0.8721
4. 0.0246
5. 0.9146
Resolução
```
# Por exemplo:
# seed <- 3185
# n <- 30
# repl <- 5000
# lambda0 <- 1
# k <- 18.24

set.seed(seed)

rejeicoes <- 0
for(i in 1:repl) {
  rejeicoes <- rejeicoes + (2 * lambda0 * sum(rexp(n, rate = lambda0)) < k)
}

solution <- rejeicoes / repl
```
O valor calculado é igual a 0.0246, ou seja, o teste de hipóteses aplicado garante que a probabilidade da hipótese $H_0$ não ser rejeitada caso seja verdadeira é aproximadamente igual a 0.9754.
Uma engenheira civil pretende averiguar se o tempo (em horas) para carregar (total ou parcialmente) um camião basculante numa obra segue uma distribuição exponencial com parâmetro $\lambda=0.8$ (hipótese $H_0$).

A recolha de uma amostra casual de dimensão $50$ conduziu ao seguinte conjunto de dados:

1.486, 1.327, 2.208, 0.481, 0.339, 1.049, 0.083, 0.034, 0.215, 1.585, 0.487, 1.531, 2.774, 1.874, 0.232, 0.954, 2.990, 0.522, 0.374, 0.395, 0.685, 3.117, 4.153, 0.405, 0.273, 0.746, 2.596, 0.008, 0.481, 0.281, 0.016, 0.976, 0.929, 1.172, 0.009, 2.472, 0.314, 0.906, 0.058, 1.730, 1.785, 1.774, 4.907, 0.220, 1.567, 2.170, 7.315, 3.235, 0.858, 3.714
1. Divida o suporte da variável aleatória $\left(\mathbb{R}^+\right)$ em $k = 8$ classes equiprováveis sob $H_0$.
2. Agrupe os dados observados nas classes definidas em 1., obtendo o conjunto de frequências absolutas observadas associadas a essas classes.
3. Recorra às frequências obtidas em 2. para testar $H_0$ e indique o valor-p do teste de ajustamento do qui-quadrado, arredondado a 4 casas decimais.
Resolução
```
# Por exemplo:
# dados <- c(0.956, 0.966, 0.79, 0.666, 0.247, 0.423, 0.641, 0.948, 0.015, 0.03,
# 0.745, 2.606, 0.061, 7.375, 1.346, 0.426, 1.61, 0.542, 1.054, 0.122, 0.076,
# 3.532, 2.76, 2.114, 3.047, 0.659, 0.153, 0.142, 0.838, 3.543, 0.615, 1.903,
# 1.166, 0.332, 0.025, 0.666, 0.299, 0.872, 2.717, 0.709, 1.472, 1.232, 0.281,
# 1.086, 0.227, 1.781, 0.422, 1.105, 0.087, 0.541, 0.064, 0.4, 1.724, 0.327, 1.42,
# 0.077, 1.948, 0.085, 0.391, 0.387)
# k <- 6
# lambda <- 0.5

probs <- 0:k / k
limites <- qexp(probs, rate = lambda)

classes <- cut(dados, breaks = limites)
freq_observadas <- table(classes)

probs <- rep(1 / k, k)
solution <- chisq.test(freq_observadas, p = probs)$p.value
```
O valor-p do teste do qui-quadrado é 0.9286.