Projeto computacional de PE

O ficheiro RenewEnergConsumption.xls contém valores referentes ao consumo de energias renováveis (% do total do consumo de energia) em vários países do mundo (base de dados: “SE4ALL Global Tracking Framework” do “Energy Sector Management Assistance Program”, Banco Mundial e Agência Internacional de Energia).

Compare o consumo de energias renováveis nos países Romania e United States, com base nos dados do período 1990–2019, construindo um diagrama de dispersão.

Submeta um ficheiro em formato PDF, com uma única página, que inclua:
1. O código em R.
  
  Nota: no código devem também constar os comandos para leitura e seleção dos dados do ficheiro.
2. O gráfico produzido.
Resolução
```
library(ggplot2)
theme_set(theme_light())

# renovaveis <- readxl::read_excel("RenewEnergConsumption.xls", skip = 3)
# pais <- c("Turkiye", "Japan") # por exemplo

renovaveis <- renovaveis[renovaveis$`Country Name` %in% pais,]

primeira <- which(names(renovaveis) == "1990")
ultima <- which(names(renovaveis) == "2019")
dados <- data.frame(t(renovaveis[, primeira:ultima]))

ggplot(dados) +
  geom_point(aes(x = X1, y = X2), color = "magenta") +
  labs(x = pais[1], y = pais[2],
       title = "Consumo de energias renováveis (% do consumo total)")
```
No ficheiro de dados hicp.xlsx encontra informação fornecida pela Eurostat relativa à HICP (Harmonised Index of Consumer Prices). Este índice mede as alterações nos preços ao consumidor nos 27 países de União Europeia. Os valores são mensais, desde agosto 2021 até agosto 2022, e estão divididos por categoria de produto alimentar.

Considere as categorias All-items HICP, Butter e Fresh low fat milk. Com recurso ao pacote ggplot produza um único gráfico que lhe permita fazer uma análise da evolução dos preços neste período.

Submeta um ficheiro em formato PDF, com uma única página, que inclua:
1. O código em R.
  
  Nota: no código devem também constar os comandos para leitura e seleção dos dados do ficheiro.
2. O gráfico que achar mais adequado para analisar a evolução dos preços.
3. Comentários sobre os resultados obtidos.
Resolução
```
library(ggplot2)
theme_set(theme_light())

# hicp <- readxl::read_excel("hicp.xlsx", sheet = "Sheet 1", range = "A2:C379")
# categorias <- c("All-items HICP", "Eggs", "Poultry") # por exemplo

names(hicp) <- c("Tempo", "Categoria", "HICP")
hicp$Tempo <- as.Date(paste0(hicp$Tempo, "-01"), format = "%Y-%m-%d")

hicp <- hicp[hicp$Categoria %in% categorias,]

ggplot(hicp) +
  geom_line(aes(x = Tempo, y = HICP, color = Categoria)) +
  geom_point(aes(x = Tempo, y = HICP, color = Categoria)) +
  ggtitle("Evolução do HICP entre agosto 2021 e agosto 2022")
```
O ficheiro TIME_USE_24092022.csv contém uma compilação de dados enviados por diversos países para a OCDE (Organização para a Cooperação e Desenvolvimento Económico) sobre o tempo médio diário (em minutos) despendido pelas pessoas entre os 15 e os 64 anos em diferentes tipos de ocupações.

Com recurso ao pacote ggplot produza um único gráfico de barras que permita comparar as ocupações médias do tempo, conjuntamente para homens e mulheres, registadas nos seguintes países: México, Dinamarca e Suécia.

Submeta um ficheiro em formato pdf, com uma única página, que inclua:
1. O código em R.
  
  Nota: no código devem também constar os comandos para leitura e seleção dos dados do ficheiro.
2. O gráfico produzido.
3. Comentários sobre os resultados obtidos.
Resolução
```
library(ggplot2)
theme_set(theme_light())

# url <- 'https://fenix.tecnico.ulisboa.pt/downloadFile/845043405571253/TIME_USE_24092022.csv'
# tempo <- read.table(url)
# paises <- c("México", "Portugal", "Itália") # por exemplo

dados <- tempo[tempo$Sexo == "Total" & tempo$País %in% paises,]
dados$Tempo <- as.numeric(dados$Tempo)

ggplot(dados) +
  geom_col(aes(x = País, y = Tempo, fill = Ocupação), position = "fill") +
  scale_y_continuous() +
  labs(title = "Uso do tempo diário em três países da OCDE",
       y = "Proporção média do dia")
```
Considere o seguinte jogo em que se lança um dado equilibrado com seis faces numeradas de 1 a 6. Se sair $4$ ou $6$, perde-se $41$ euros. Se sair $1$ ou $3$ nada acontece. Se sair $5$, ganha-se $41$ euros. Se sair $2$, ganha-se $9 \times 41$ euros. Gostaríamos de saber quanto dinheiro se pode esperar ganhar por jogo, em média.

Fixando a semente em $511$, simule $5200$ repetições do jogo. Calcule a média dos ganhos nessas repetições, como uma estimativa do valor esperado do ganho no jogo, e indique esse valor arredondado a 2 casas decimais.

Resolução
```
# seed <- 435 # por exemplo
# const <- 9
# valor <- 28
# reps <- 4316

set.seed(seed)
valores <- c(-valor, 0, valor, const * valor)
probs <- c(2/6, 2/6, 1/6, 1/6)
ganhos <- sample(valores, size = reps, replace = TRUE, prob = probs)

ganho.medio <- mean(ganhos)
```
O valor estimado do ganho médio no jogo é de 52.72 euros.
O objectivo desta questão é gerar distribuições arbitrárias, a partir de distribuições uniformes (por exemplo), recorrendo ao chamado método “Aceitação-Rejeição”. Neste exercício vai simular uma distribuição triangular.

A função densidade de probabilidade de uma v. a. $X$ com distribuição triangular de parâmetros $(a,b,c)$ é dada por

\[f(x)=\begin{cases} 0, & x<a\\ \frac{2(x-a)}{(b-a)(c-a)}, & a\leq x \leq c\\ \frac{2(b-x)}{(b-a)(b-c)}, & c<x\leq b\\ 0, & x>b \end{cases}\]

No caso concreto, vai gerar uma distribuição triangular de parâmetros $(a=15, b=40, c=27.5)$. Para tal deverá seguir os seguintes passos:
1. Simular um valor, $x$, proveniente de uma distribuição uniforme no intervalo $(a,b)=(15, 40)$.
2. Simular $u$, proveniente de uma v. a. com distribuição uniforme no intervalo (0, 1).
3. Se $u<\frac{f(x)}{\frac{2}{(b-a)}}=25\frac{f(x)}{2}$, aceitar $x$. Caso contrário, rejeitar.
Fixando a semente em $348$, implemente este método repetindo os passos anteriores até obter uma amostra de dimensão $n = 193$.

Submeta um ficheiro em formato PDF, com uma única página, que inclua:
1. O código em R.
2. O gráfico da função densidade de probabilidade teórica $f$ sobreposto a um histograma de frequências relativas unitárias dos valores simulados.
3. Comentários sobre os resultados obtidos.
Resolução
```
# seed <- 435 # por exemplo
# a <- 5
# b <- 24
# c <- 14.5
# n <- 150

set.seed(seed)

f <- function(x) {
  if(x <= c & x >= a) {
    res <-  2 / (b - a) * (x - a) / (c - a)
  } else {
    if(x <= b & x > c) {
      res <- 2 / (b - a) * (b - x) / (b - c)
    } else {
      res <- 0
    }
  }
  return(res)
}

sim <- vector(length = n)
aceites <- 0
M <- (b - a) / 2

while(aceites < n) {
  x <- runif(1, min = a, max = b)
  u <- runif(1)
  if(u < M * f(x)) {
    aceites <- aceites + 1
    sim[aceites] <- x
  }
}

limits <- c(a - (b - a) / 4, b + (b - a) / 4)
hist(sim, probability = TRUE, xlim = limits, ylim = c(0, 1.5 * f(c)), xlab = "X", 
     ylab = "Densidade", main = "Simulação de uma distribuição triangular por A/R")

g <- Vectorize(f)
curve(g, col = "red", add = TRUE)
```
Considere uma variável aleatória com distribuição de Cauchy com parâmetros de localização $x_0 = -0.9$ e de escala $\gamma = 6.8$, com a função densidade de probabilidade dada por: \[ f(x)= \frac{1}{\pi \gamma} \left[1 + \left( \frac{x-x_0}{\gamma} \right)^2 \right]^{-1},\, x\in\mathbb{R}.\] Usando o R e fixando a semente em $1133$, gere uma amostra dessa variável para cada dimensão $n\in \{3000, 3500, \ldots, 10000\}$ e calcule a média e a mediana de cada amostra gerada. Com recurso ao pacote ggplot construa um gráfico que apresente as médias e as medianas das amostras em função da dimensão da amostra.

Submeta um ficheiro em formato PDF, com uma única página, que inclua:
1. O código em R.
2. O gráfico produzido.
Resolução
```
# seed <- 1215 # por exemplo
# x0 <- 1
# gamma <- 0.95

library(ggplot2)
theme_set(theme_light())

set.seed(seed)

n <- seq(3000, 10000, 500)
res <- data.frame()

for (i in n) {
  amostra <- rcauchy(i, location = x0, scale = gamma)
  res <- rbind(res, c(i, mean(amostra)))
  res <- rbind(res, c(i, median(amostra)))
}

colnames(res) <- c("n", "Valor")
res$Medida <- rep(c("Média", "Mediana"), length.out = nrow(res))

ggplot(res, aes(x = n, y = Valor, color = Medida)) +
  geom_point() +
  geom_line() +
  labs(title = "Simulação de média e mediana de uma distribuição de Cauchy ")
```
Pretende-se simular o erro na estimação do valor médio de uma variável aleatória $X$ com distribuição Normal de valor esperado $\mu$ e variância $\sigma^2$ desconhecidos. Para isso, considere os seguintes três estimadores:
1. A média amostral
2. A mediana amostral
3. $\dfrac{\min\{X_i\}+\max\{X_i\}}{2}$
Usando o R e fixando a semente em $447$, gere $m=3935$ amostras de dimensão $n=28$ dessa variável, admitindo que $\mu = -10$ e $\sigma = 3.4$.
1. Para cada amostra, calcule o valor de cada um dos três estimadores, obtendo assim $m=3935$ estimativas associadas a cada um dos estimadores.
2. Calcule as médias das $m=3935$ estimativas associadas a cada um dos três estimadores.
3. Para cada estimador, calcule o valor absoluto da diferença entre a média respetiva calculada no ponto 2. e o valor teórico do parâmetro.
4. Calcule o menor desses três valores absolutos.
5. Apresente o último valor calculado arredondado a 6 casas decimais.
Resolução
```
# seed <- 508 # por exemplo
# m <- 7011
# n <- 9
# mu <- -11
# sigma <- 4.5

# Simulação
set.seed(seed)

res <- data.frame()
for (i in 1:m) {
  amostra <- rnorm(n, mean = mu, sd = sigma)
  res <- rbind(res, c(mean(amostra), median(amostra), mean(range(amostra))))
}

out <- min(abs(colMeans(res) - mu))
```
O menor desvio obtido é 0.005814.
Considere uma variável aleatória $X$ com distribuição Normal de valor esperado $\mu$ e variância $\sigma^2$, com valores desconhecidos. Determine o intervalo aleatório de confiança para $\mu$ com um nível de confiança $\gamma = 0.9$.

Usando o R e fixando a semente em $1728$, gere $m=250$ amostras de dimensão $n=14$ dessa variável, admitindo que $\mu = -1$ e $\sigma = 2.9$. Para cada amostra gerada, calcule o respetivo intervalo de confiança para $\mu$. Indique a proporção de intervalos de confiança gerados que contém o valor do parâmetro $\mu = -1$, arredondando a 3 casas decimais.

Resolução
```
# seed <- 1435 # por exemplo
# n <- 15
# m <- 200
# mu <- 0
# sigma <- 2
# gamma <- 0.95

set.seed(seed)

k <- qt((1 + gamma) / 2, n - 1)
res <- data.frame()
for (i in 1:m) {
  amostra <- rnorm(n, mean = mu, sd = sigma)
  media <- mean(amostra)
  ampl <- sd(amostra) / sqrt(n) * k
  res <- rbind(res, c(media - ampl, media + ampl))
}
names(res) <- c("inferior", "superior")
res$Sucesso <- ((res$inferior - mu) * (res$superior - mu) <= 0)
suc <- sum(res$Sucesso) / m
```
A proporção de intervalos de confiança gerados que contém o valor do parâmetro $\mu$ é 0.892.
Considere uma variável aleatória $X$ com distribuição Normal de valor esperado $\mu$ e variância $\sigma^2 = 2.6$. Construa o teste para a hipótese $H_0:\mu = 1.3$ contra $H_1:\mu > 1.3$ ao nível de significância $\alpha = 0.07$.

Usando o R e fixando a semente em $1139$, gere $m=100$ amostras de dimensão $n=20$ dessa variável, admitindo que $\mu = 1.9$. Aplique o teste de hipóteses que construiu a cada amostra gerada e use o conjunto dos resultados para calcular uma estimativa da probabilidade do teste conduzir à não rejeição de $H_0$ quando $\mu = 1.9$. Indique o resultado com 3 casas decimais.

Resolução
```
# seed <- 1435 # por exemplo
# n <- 15
# m <- 200
# mu0 <- 0
# mu1 <- 1
# sigma2 <- 2
# alpha <- 0.05

set.seed(seed)

k <- mu0 + qnorm(1 - alpha) * sqrt(sigma2 / n)
res <- 0
for (i in 1:m) {
  amostra <- rnorm(n, mean = mu1, sd = sqrt(sigma2))
  res <- res + (mean(amostra) <= k)
}
res <- res / m
```
A estimativa da probabilidade do teste conduzir à não rejeição de $H_0:\mu = 1.3$ quando $\mu = 1.9$ é 0.430.
Em 2008 a empresa M&M anunciou que a proporção de cores dos diferentes chocolates produzidos por esta empresa era: 24% azul, 20% laranja, 16% verde, 14% amarelo, 13% vermelho, 13% castanho. Uns anos mais tarde a empresa retirou esta informação do seu site e, desde então, nunca mais a atualizou.

Considere os seguintes dados relativos a uma amostra recolhida em 2022, de dimensão $465$:

Azul Laranja Verde Amarelo Vermelho Castanho

120 83 79 55 69 59

Indique o valor-p, arredondado a 4 casas decimais, correspondente ao teste para verificar se as proporções de cores de chocolates observadas em 2022 estão de acordo com as anunciadas pela M&M em 2008.

Resolução
```
# cores <- c(121, 84, 79, 56, 70, 57) # por exemplo
probs <- c(0.24, 0.2, 0.16, 0.14, 0.13, 0.13)

res <- chisq.test(cores, p = probs)$p.value
```
O valor-p do teste do qui-quadrado é 0.4404.

Projeto computacional de PE

Resolução

Resolução

Resolução

Resolução

Resolução

Resolução

Resolução

Resolução

Resolução

Resolução