Probabilidade e Estatística

Projeto computacional

  1. O ficheiro Paises_PIB_ICH.csv contém um conjunto de dados com as seguintes variáveis:

    Após a leitura do ficheiro no R, selecione todos os países dos seguintes continentes: Africa e Europe (Nota: os nomes dos países e dos continentes estão em Inglês no ficheiro).

    Com recurso ao pacote ggplot2, produza um único diagrama de dispersão que represente os valores do Índice de Capital Humano em função do PIB per capita e tal que:

    1. seja possível identificar o continente a que cada país pertence;

    2. o PIB esteja representado em escala logaritmica;

    3. todos os paises sejam representados por um símbolo e, além disso, os seguintes países: Rwanda, Gabon, Portugal e Serbia, apareçam representados pelo símbolo e identificados pelo nome.

    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 contenha:

    1. O código em R, incluindo os comandos para leitura dos dados do ficheiro e a construção do gráfico propriamente dito.

    2. O gráfico produzido.


    Resolução

    library(ggplot2)
theme_set(theme_classic())

# url <- 'https://web.tecnico.ulisboa.pt/~paulo.soares/pe/projeto/Paises_PIB_ICH.csv'
# dados <- read.csv(url)
# # Por exemplo:
# continentes <- c("Americas", "Africa")
# paises <- c("Antigua and Barbuda", "Grenada", "Ivory Coast", "Comoros")

subcjt <- subset(dados, Continent %in% continentes)
etiquetados <- subset(subcjt, Country %in% paises)

ggplot(subcjt, aes(x = GDP, y = HCI)) +
  geom_point(aes(color = Continent), size = 2, alpha = 0.8) +
  geom_label(data = etiquetados, aes(label = Country), size = 2.5, alpha = 0) +
  scale_x_log10() +
  labs(title = "Human Capital Index versus GDP per capita",
       x = "GDP per capita adjusted (international dollars, IMF, 2023)",
       y = "Human Capital Index (World Bank, 2020)")


  2. O ficheiro master.csv consiste de um conjunto de dados (obtidos no kaggle) sobre suicídios em 101 países, incluindo nomeadamente as variáveis:

    Após a leitura desse ficheiro no R, selecione os dados referentes ao ano de 2001 e ao grupo etário 15-24 years (em Inglês no ficheiro).

    Com recurso ao pacote ggplot2, produza um único gráfico com dois diagramas de extremos e quartis do número de suicídios por 100000 habitantes da população de cada país que permitam comparar homens e mulheres.

    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_minimal())

# url <- 'https://web.tecnico.ulisboa.pt/~paulo.soares/pe/projeto/master.csv'
# dados <- read.csv(url, check.names = FALSE)
# # Por exemplo:
# ano <- 2016
# idade <- "15-24 years"

subcjt <- subset(dados, year == ano & age == idade)

ggplot(subcjt) +
  geom_boxplot(aes(sex, `suicides/100k pop`)) + 
  labs(title = paste("Rate of suicides for age group", idade, "in", ano))


  3. O ficheiro electricity.xlsx contém dados coletados pela Agência Internacional da Energia sobre a produção mensal de eletricidade desde 2010 até 2022 em vários países. A energia produzida foi medida em gigawatt-horas (GWh) e é proveniente de diferentes fontes incluindo as produções hídrica, eólica, solar, geotérmica, nuclear e a partir de combustíveis fósseis, entre outras.

    Recorrendo ao pacote ggplot2, produza um único gráfico que permita ilustrar a evolução mensal da proporção (share) de energia elétrica produzida a partir de fontes renováveis (Renewables, no ficheiro) desde o início de 2015. O gráfico deve incluir os valores médios para todos os países incluídos (IEA Total, no ficheiro) e os países Czech Republic e Colombia (em Inglês no ficheiro). O gráfico deve ter ainda as seguintes características:

    1. a proporção de energia elétrica produzida mensalmente a partir de fontes renováveis deve ser expressa em percentagem;

    2. a escala do eixo das ordenadas deve cobrir todo o intervalo \([0, 100]\).

    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_light())

# url <- 'https://web.tecnico.ulisboa.pt/~paulo.soares/pe/projeto/electricity.xlsx'
# download.file(url, destfile = "local.xlsx")
# eletro <- readxl::read_xlsx("local.xlsx")
# # Por exemplo:
# paises <- c("Portugal", "Mexico", "IEA Total")

dados <- subset(eletro, COUNTRY %in% paises & YEAR >= 2015 & PRODUCT == "Renewables")
dados$DATE <- as.Date(paste(1, dados$MONTH, dados$YEAR), format = "%d %m %Y") 
dados$share <- as.numeric(dados$share) * 100

ggplot(dados) +
  geom_line(aes(DATE, share, color = COUNTRY), linewidth = 1) +
  ylim(c(0, 100)) +
  labs(title="Monthly share of renewable sources in electricity production",
       y = "RENEWABLE SOURCES (% of total production)")


  4. Um sistema é formado por \(6\) circuitos eléctricos, em que cada um emite sinais codificados no conjunto \(\{1,2,\ldots,12\}\). Cada circuito emite o sinal \(i\) com probabilidade \(\frac{i}{78}\), para \(i \in \{1,2,\ldots,12\}\), independentemente dos restantes circuitos. Se pelo menos um dos circuitos emitir o sinal 2 é produzido um aviso sonoro e, caso pelo menos um dos circuitos emita o sinal 1, o sistema é desligado.

    Fixando a semente em \(2030\), simule \(150\) realizações do estado de um sistema com as características acima descritas e calcule a proporção de vezes em que é produzido um aviso sonoro num sistema que não é desligado. Essa proporção arredondada a 2 casas decimais é igual a:


    1. 0.00
    2. 0.14
    3. 0.07
    4. 0.04
    5. 0.06

    Resolução

    # # Por exemplo:
# circuitos <- 8
# sinais <- 6
# simulacoes <- 500
# seed <- 1234

set.seed(seed)

ligado <- 0
aviso <- 0

probs <- 1:sinais * 2 / (sinais * (sinais + 1))

for (i in 1:simulacoes) {
  replica <- sample(sinais, size = circuitos, replace = TRUE, prob = probs)
  minimo <- min(replica)
  if (minimo > 1) {
    ligado <- ligado + 1
    if (minimo == 2) aviso <- aviso + 1
  }
}
solution <- aviso / ligado

    A proporção pedida é igual a 0.14.


  5. Sejam \((Z_1,Z_2,\ldots,Z_{n+1})\) variáveis aleatórias independentes e com distribuição normal padrão. Prova-se que a variável aleatória, \[T=\sqrt{n}\frac{Z_1}{\sqrt{\sum_{i=2}^{n+1} Z_i^2}}\] tem uma distribuição \(t\)-Student com \(n\) graus de liberdade.

    Como confirmação empírica, pretende-se que avalie a probabilidade \(p = P(T\leq -1.4)\), procedendo da seguinte forma:

    1. considere \(n = 24\) e, fixando a semente em \(1704\), gere \(r=100\) amostras formadas por \(m=140\) valores de \(T\), à custa da geração de variáveis normais padrão;

    2. para cada amostra gerada, determine a proporção de valores menores ou iguais a \(-1.4\) e, finalmente, calcule a média de essas proporções como uma aproximação de \(p\).

    Calcule o valor absoluto da diferença entre a aproximação anterior e o valor obtido através da função do R que permite calcular diretamente a probabilidade de uma variável aleatória com distribuição \(t_{(24)}\) tomar valores menores ou iguais do que \(-1.4\). Multiplique o valor obtido por 100 e indique esse resultado final arredondado a 5 casas decimais.


    Resolução

    # # Por exemplo:
# valor <- 2
# r <- 300
# n <- 15
# m <- 100
# seed <- 1234

set.seed(seed)
prop <- vector(length = r)

for(j in 1:r) {
  prop[j] <- 0
  for(i in 1:m) {
    z <- rnorm(n + 1)
    t <- z[1] / sqrt(sum(z[-1]^2) / n)
    prop[j] <- prop[j] + (t <= valor)
  }
}
prop <- prop / m

solution <- abs(pt(valor, n) - mean(prop)) * 100

    O valor pedido é 0.15602.


  6. Considere que \(X_1, \ldots, X_{35}\) são variáveis aleatórias independentes e identicamente distribuídas de acordo com uma distribuição exponencial com valor esperado \(a = 3\), onde \(X_i\) denota a duração da componente eletrónica \(i\), \(i=1,...,35\).

    Para o cálculo da fiabilidade da duração total das 35 componentes eletrónicas em 81, ou seja, o cálculo de \(P(Y > 81)\), onde \(Y=\sum_{i=1}^{35}{X_i}\), considere as seguintes abordagens:

    1. Valor simulado: Fixando a semente em 1985, gere 1000 amostras de dimensão \(n=35\) da distribuição exponencial referida e calcule quer o valor simulado de \(Y\) para cada uma das amostras quer a proporção de valores simulados de \(Y\) que são maiores do que 81.

    2. Valor exato: Sabendo que a distribuição exata de \(Y\) é uma distribuição gama com parâmetros de forma \(35\) e taxa (ou rate) \(1/{3}\), obtenha o valor da fiabilidade da duração total das 35 componentes eletrónicas em 81.

    Obtenha a probabilidade em causa em cada uma das abordagens e indique o valor absoluto da diferença entre os resultados obtidos com base nas abordagens 1. e 2., multiplicado por 100 e com 4 casas decimais.


    Resolução

    # # Por exemplo:
# seed <- 1966
# n <- 30
# a <- 2
# valor <- 63

set.seed(seed)
y <- vector(length = 1000)

for (i in 1:1000) {
  y[i] <- sum(rexp(n, rate = 1 / a))
}
sim <- sum(y > valor) / 1000

ext <- 1 - pgamma(valor, n, rate = 1 / a)

solution <- abs(sim - ext) * 100

    O valor pedido é igual a 0.6663.


  7. Um gestor de risco de uma companhia de seguros assume que o montante de indemnização (em milhares de euros) de determinada carteira de apólices é representado pela variável aleatória \(X\) com função de densidade de probabilidade dada por \(\theta x^{-\theta-1} a^\theta\), para \(x \geq a\), onde \(a=5\) e \(\theta\) é um parâmetro com um valor positivo desconhecido relacionado com a ocorrência de indemnizações elevadas.

    A concretização de uma amostra aleatória de dimensão \(n=10\) proveniente de \(X\) conduziu ao seguinte conjunto de observações:

    5.05, 8.56, 5.65, 5.28, 5.32, 5.12, 5.16, 7.16, 5.30, 6.12

    1. Recorra à função mle do pacote base stats4 para obter a estimativa de máxima verosimilhança de \(\theta\), usando o valor 4.6 como valor inicial da pesquisa numérica e sem alterar qualquer outro argumento opcional de tal função.

    2. Utilize a estimativa obtida em 1. para determinar a estimativa de máxima verosimilhança do quantil de probabilidade \(p = 0.75\) de \(X\).

    3. Calcule o desvio absoluto entre a estimativa obtida em 2. e o verdadeiro valor do quantil, quando \(\theta=4.6\). Apresente o resultado com 4 casas decimais.


    Resolução

    # # Por exemplo:
# data <- c(4.98, 4.69, 6.01, 5.78, 7.21, 4.88, 4.28, 5.75, 4.41, 5.8, 4.26, 4.93)
# valor_inicial <- 4.6
# p <- 0.25
# theta <- 4.6
# a <- 4

library(stats4)

soma_logdata <- sum(log(data))
n <- length(data)

fun <- function(theta) {
  res <- -n * log(theta) + (theta + 1) * soma_logdata - n * theta * log(a)
  return(res)
}

result <- mle(minuslog = fun, start = list(theta = valor_inicial))
est_theta <- as.vector(result@coef)

quantil <- function(p, theta) {
  res <- a * (1 - p)^(-1 / theta)
}

estimativa <- quantil(p, est_theta)
valor_exato <- quantil(p, theta)

solution <- abs(estimativa - valor_exato)

    O desvio absoluto pretendido é 0.6360.


  8. Uma engenheira mecânica assume que a distância (em metros) de travagem, ao utilizar-se determinadas pinças de travão, é uma variável aleatória \(X\) com distribuição normal de valor esperado desconhecido \(\mu\) e variância desconhecida \(\sigma^2\).

    1. Fixando a semente em 1751, extraia ao acaso e sem reposição \(n= 10\) observações da seguinte amostra casual de distâncias de travagem em condições controladas:

      33.5, 31.5, 35.4, 36.3, 33.3, 36.9, 35.3, 31.5, 31.8, 32.9, 34.0, 36.0, 33.0, 39.4, 36.5, 32.9, 31.8, 37.6

    2. Obtenha o intervalo de confiança a um nível de confiança \(\gamma= 0.94\) para \(\sigma^2\) com base na amostra selecionada em 1. e nos quantis de probabilidade \(a=F_{\chi_{(n-1)}^2}^{-1}\left(\frac{1-\gamma}{2}\right)\) e \(b=F_{\chi_{(n-1)}^2}^{-1}\left(\frac{1+\gamma}{2}\right)\).

    3. De modo a minimizar a amplitude esperada do intervalo aleatório de confiança para \(\sigma^2\) ao nível de confiança anterior, considere o par de quantis de probabilidade \((c,d)\) satisfazendo as duas equações seguintes: \[\begin{cases} F_{\chi_{(n-1)}^2}(d)-F_{\chi_{(n-1)}^2}(c) = \gamma\\ f_{\chi_{(n+3)}^2}(d)-f_{\chi_{(n+3)}^2}(c) = 0 \end{cases}\] Para obter \((c, d)\) instale o pacote pracma e use a função fsolve considerando \((a,b)\) como o valor inicial da pesquisa numérica e sem utilizar qualquer outro argumento opcional dessa função.

      Utilize estes quantis para determinar um novo intervalo de confiança para \(\sigma^2\) com base nos mesmos dados usados em 2.

    4. Indique o valor da diferença, arredondada a três casas decimais, entre as amplitudes dos intervalos de confiança calculados em 2. e 3.


    Resolução

    # # Por exemplo:
# seed <- 1966
# amostra <- c(31.9, 34.7, 30.5, 38.4, 37.1, 34.7, 32.7, 31.7, 36.3, 33.1, 28.9, 33.5,
# 32,  29.2, 33.4, 32.2, 32.8, 33.4, 29.9)
# gama <- 0.96
# n <- 9

set.seed(seed)
dados <- sample(amostra, n)

a <- qchisq((1 - gama) / 2, n - 1)
b <- qchisq((1 + gama) / 2, n - 1)
amp_ic <- (n - 1) * var(dados) * (1 / a - 1 / b)

library(pracma)

restricoes <- function(x) {
  F1 <- pchisq(x[2], n - 1) - pchisq(x[1], n - 1) - gama
  F2 <- dchisq(x[2], n + 3) - dchisq(x[1], n + 3)
  return(c(F1, F2))
}

raiz <- fsolve(restricoes, c(a, b))$x

amp_ic_min <- (n - 1) * var(dados) * (1 / raiz[1] - 1 / raiz[2])

solution <- amp_ic - amp_ic_min

    O valor pretendido é 2.155.


  9. Seja \(\left(X_1, \ldots, X_{250}\right)\) uma amostra aleatória de uma população com distribuição de Poisson de parâmetro \(\lambda\). Para testar a hipótese \(H_0: \lambda=2.70\) contra \(H_1: \lambda = 2.90\) construiu-se um teste de hipóteses que conduz à rejeição de \(H_0\) quando \(\bar{X}>k\) com \(k=2.871\).

    Fixando a semente em \(5621\), gere \(m=10000\) pares de amostras de dimensão \(n=250\) da distribuição de Poisson, uma sob \(H_0\) e outra sob \(H_1\) em cada par, e calcule as frequências relativas de cada um dos dois tipos de decisões erradas a que o teste conduz. Com base nos resultados obtidos, calcule um valor aproximado do quociente entre a probabilidade de erro de 2ª espécie e a probabilidade de erro de 1ª espécie. Selecione a resposta correta de entre as seguintes opções:


    1. 2.22
    2. 13.20
    3. 7.22
    4. 6.22
    5. 9.72

    Resolução

    # # Por exemplo:
# seed <- 4600
# n <- 200
# m <- 70000
# lambda0 <- 2.60
# lambda1 <- 2.85
# k <- 2.764

set.seed(seed)

falsa_aceitacao <- 0
falsa_rejeicao <- 0

for(i in 1:m) {
  falsa_rejeicao <- falsa_rejeicao + (mean(rpois(n, lambda0)) > k)
  falsa_aceitacao <- falsa_aceitacao + (mean(rpois(n, lambda1)) <= k)
}

solution <- falsa_aceitacao / falsa_rejeicao

    O valor calculado é igual a 7.22, ou seja, a probabilidade do teste não rejeitar erradamente \(H_0\) é aproximadamente 7.22 vezes maior do que a probabilidade de rejeitar erradamente essa hipótese.


  10. Seja \(X\) o tempo, em horas, que uma dada tarefa demora a ser executada numa empresa. O director desta empresa considera que \(X\) tem distribuição triangular simétrica com valor mínimo \(a = 4.75\) e valor máximo \(b = 14\), cuja função de distribuição é \[F_X(x)=\begin{cases} 0, & x < a\\ 2\left(\frac{x-a}{b-a}\right)^2, & a\leq x< \frac{a+b}{2}\\ 1-2\left(\frac{b-x}{b-a}\right)^2, & \frac{a+b}{2}\leq x<b\\ 1, & x \geq b \end{cases}.\]

    Para testar a sua hipótese (\(H_0\)) o director colecionou \(n = 110\) tempos de execução da tarefa, retirados ao acaso do historial da empresa, o que conduziu aos seguintes resultados:

    6.27, 10.37, 11.39, 9.04, 7.21, 8.09, 9.74, 13.28, 11.76, 6.38, 7.16, 11.20, 6.73, 10.81, 10.43, 9.97, 11.31, 9.27, 7.62, 9.68, 13.48, 9.72, 6.23, 8.18, 8.76, 7.59, 12.09, 7.02, 10.15, 7.86, 12.12, 8.03, 11.48, 9.65, 10.57, 8.07, 5.34, 12.79, 9.38, 12.57, 9.92, 13.40, 10.12, 11.77, 7.67, 8.02, 9.21, 10.81, 9.43, 9.68, 11.50, 8.93, 7.36, 5.86, 9.97, 6.91, 8.79, 10.22, 7.63, 13.29, 11.34, 9.74, 10.56, 8.02, 7.44, 8.55, 11.32, 9.83, 7.34, 9.47, 8.72, 8.93, 9.99, 11.05, 11.79, 11.43, 13.24, 7.07, 10.98, 8.68, 7.30, 8.65, 10.15, 9.24, 7.64, 6.46, 10.64, 7.82, 6.38, 7.85, 11.04, 7.58, 11.75, 10.47, 9.95, 10.08, 7.57, 9.34, 8.68, 9.39, 10.24, 8.52, 8.05, 11.06, 7.94, 10.93, 10.50, 10.69, 11.50, 10.00

    1. Divida o suporte da variável aleatória \(\left[a, b\right]\) em \(k = 6\) classes, com igual amplitude.

    2. Agrupe os valores da amostra observada 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.


    1. 0.2398
    2. 0.5885
    3. 0.7608
    4. 0.1795
    5. 0.9879

    Resolução

    # # Por exemplo:
# k <- 8
# a <- 5
# b <- 12
# amostra <- c(9.01, 9.27, 8.32, 9.80, 9.61, 5.40, 8.79, 7.86, 7.08, 8.91, 11.57,
# 10.10, 6.67, 10.24, 7.38, 8.26, 5.61, 7.31, 10.84, 6.37, 8.53, 11.18, 6.79, 9.20,
# 11.62, 9.01, 8.24, 7.76, 8.55, 8.74, 11.92, 7.09, 9.67, 11.03, 9.88, 9.93, 5.98,
# 9.91, 7.14, 7.26, 10.12, 10.13, 9.71, 8.26, 6.89, 7.61, 10.71, 6.68, 7.83, 8.87,
# 12.14, 7.87, 10.61, 10.23, 6.62, 8.98, 8.79, 8.51, 6.41, 9.33, 7.63, 7.65, 7.10,
# 7.98, 9.05, 10.66, 7.27, 11.71, 6.36, 6.08, 11.78, 9.45, 9.00, 8.99, 12.22, 6.80,
# 7.47, 8.18, 7.99, 7.48, 10.08, 11.42, 6.86, 9.78, 9.81, 9.19, 7.50, 11.80, 11.18,
# 9.19, 10.22, 11.97, 10.55, 10.55, 8.62, 10.15, 7.35, 9.57, 10.82, 10.00, 8.13,
# 5.44, 6.30, 6.83, 7.41, 8.85, 7.07, 9.93, 9.79, 6.65, 11.83, 10.78, 7.52, 7.77,
# 7.68, 7.41, 9.39, 5.79, 10.24, 7.56, 6.65, 10.10, 7.95, 8.28, 9.17, 9.91, 8.79,
# 8.90, 7.39, 10.86, 8.06, 11.43, 7.54, 6.31, 8.03, 7.40, 7.90, 11.39, 9.02, 8.32,
# 5.87, 8.57, 9.61, 9.99, 12.22, 9.69, 11.97, 9.53, 6.86, 6.92, 7.52, 9.44, 8.49,
# 6.93, 6.87, 7.51, 7.48, 8.94, 9.21, 6.56, 10.17, 10.15, 9.35, 9.92, 9.74, 10.71,
# 8.54, 10.66, 7.93, 9.88)

delta <- (b - a) / k
lim.classes <- a + delta * 0:k
classes <- cut(amostra, breaks = lim.classes)
freq.obs <- table(classes)

Fd <- function(x) {
  res <- ifelse(x < (a + b) / 2,
                2 * ((x -a) / (b - a))^2,
                1 - 2 * ((b - x) / (b - a))^2)
  return(res)
}
prob.H0 <- diff(Fd(lim.classes))

solution <- chisq.test(freq.obs, p = prob.H0)$p.value

    O valor-p do teste do qui-quadrado, arredondado a 4 casas decimais, é 0.7608.