Probabilidade e Estatística

Projeto computacional

O relatório intitulado “Crime in the United States by Region, Geographic Division and State, 2019” fornece uma visão detalhada da ocorrêncios de vários tipos de crimes nos Estados Unidos da América em 50 dos seus estados no ano de 2019. Parte desse relatório, disponível na página oficial do FBI, deu origem ao conjunto de dados em estudo, disponibilizado no ficheiro Data_FBI_USA_Crime.csv.

Estes dados são usualmente usados para calcular a taxa de um tipo de crime, utilizando estimativas populacionais provisórias do U.S. Census Bureau. Ou seja, a taxa de um crime numa região define-se como o número de ocorrências do crime por cada 100 mil habitantes dessa região.

Com recurso ao pacote ggplot2 em R, produza um gráfico que represente, através de diagramas de caixa (boxplots) paralelos, a variação da taxa do crime Burglary em 2019 agrupados pelo tipo de área (‘Cities outside metropolitan areas’, ‘Metropolitan statistical area’, ‘Nonmetropolitan counties’).

Submeta um ficheiro em formato PDF com uma única página A4, mantendo todo o texto da figura em inglês e incluindo:
- O código em R utilizado para produzir o gráfico, incluindo todos os comandos necessários para a leitura, seleção e visualização dos dados.
- A figura produzida.
Resolução
```
# Por exemplo:
# variavel <- "Murder"

# url <- 'http://web.tecnico.ulisboa.pt/~ist13493/Projecto26/Data_FBI_USA_Crime.csv'
# dados <- read.csv(url)

library(ggplot2)

ggplot(dados) +
  geom_boxplot(aes(x = Area, y = .data[[variavel]] / Population * 100000), fill = "cyan") +
  labs(title = "Crime in the United States by Area, 2019") +
  theme_light()
```
O ficheiro WUP2018-F02-Proportion_Urban.xlsx contém valores referentes à percentagem anual da população a residir em zonas urbanas até 2015, e projeções até 2050, por regiões e países (Fonte: Nações Unidas, 2018).

Com recurso ao pacote ggplot2, produza um cronograma que represente a evolução temporal da percentagem da população residente em zonas urbanas, desde 1950 até 2050, para os países Portugal e New Zealand, comparando também com os valores globais (WORLD). No gráfico devem ser distinguidos os dois períodos, até 2015 e após 2015. Por simplicidade mantenha o texto do gráfico em inglês.

Submeta um ficheiro em formato PDF, com uma única página A4, que inclua:
1. O código em R utilizado para produzir o gráfico, incluindo todos os comandos necessários para a leitura, seleção e visualização dos dados.
2. O gráfico produzido.
Resolução
```
# Por exemplo:
# paises <- c("Japan", "Portugal")

# dados <- readxl::read_xlsx('WUP2018-F02-Proportion_Urban.xlsx')

library(ggplot2)

subcjt <- subset(dados, Region %in% c(paises, 'WORLD'))

ggplot(subcjt, aes(Year, Population, color = Region)) +
  geom_line(data = subset(subcjt, Year <= 2015), linewidth = 1.5) +
  geom_line(data = subset(subcjt, Year >= 2015), linetype = 2,
            linewidth = 1.1, alpha = 1) +
  labs(title = 'Resident population in urban areas', y = '% of total population',
       subtitle = '(dashed years are predictions)') +
  theme_linedraw()
```
O conjunto de dados QualidadeAR.txt contém medições horárias de qualidade do ar recolhidas por cinco sensores químicos de óxidos metálicos, instalados numa zona urbana italiana fortemente poluída. Inclui as concentrações dos principais contaminantes como CO, NOx, NO2, hidrocarbonetos não metânicos e benzeno; os valores em falta estão assinalados como -200. (Fonte:UCI).

Após a leitura desse ficheiro no R, selecione todos os dados referentes à hora 10:00:00 do ano 2004.

Utilizando o pacote ggplot2, produza um único gráfico que combine o histograma das frequências relativas acumuladas para a variável T com o polígono das frequências relativas acumuladas associado a essa mesma variável. Use o número classes dado pelo método de Freedman-Diaconis.

Considere o seguinte:
- O texto no ficheiro de dados está em inglês; por simplicidade, mantenha todos os elementos gráficos em inglês.
- Se a sua variável apresentar valores em falta (-200), substitua esses valores pela média dos restantes valores dessa variável.
Submeta um ficheiro em formato PDF, com uma única página A4, que inclua:
1. O código em R utilizado para produzir o gráfico, incluindo todos os comandos necessários para a leitura, seleção, transformação e visualização dos dados.
2. O gráfico produzido.
Resolução
```
# Por exemplo:
# ano <- 2004
# hora <- "00:00:00"
# variavel <- "PT08.S1(CO)"

# url <- "http://web.tecnico.ulisboa.pt/~ist13493/Projecto26/QualidadeAR.txt"
# dados <- read.delim(url, dec = ",", check.names = FALSE)

library(ggplot2)

dados[dados == -200] <- NA
dados$Date <- as.Date(dados$Date,  format = "%d-%m-%Y")

subcjt <- subset(dados, Date >= as.Date(paste0(ano, "-01-01")) &
                         Date <= as.Date(paste0(ano, "-12-31")) &
                         Time == hora)

subcjt[[variavel]][is.na(subcjt[[variavel]])] <- mean(subcjt[[variavel]], na.rm = TRUE)

nbins <- nclass.FD(subcjt[[variavel]])

ggplot(subcjt, aes(.data[[variavel]])) +
  geom_histogram(aes(y = after_stat(cumsum(count) / sum(count))), 
                 color = "black", fill = "lightblue", bins = nbins) +
  geom_freqpoly(aes(y = after_stat(cumsum(count) / sum(count))), 
                color = "tomato", bins = nbins, linewidth = 1) +
  labs(title = paste("Daily distribution of", variavel, "at", hora, "in", ano),
       y = "Relative frequency") +
  theme_classic()
```
Baseando-se em experiências passadas, pode afirmar-se que a função densidade de probabilidade das vendas semanais de gasolina, $X$ (em dezenas de milhares de litros), num determinado posto de abastecimento é dada por:

\[ f_X(x) = \begin{cases} \frac{x-a}{b^2}, &a \leq x < a+b, \\ \frac{a+2\,b-x}{b^2} , &a+b \leq x \leq a+2\,b, \\ 0, &\text{caso contrário}. \end{cases} \]

onde $a = 1$ e $b = 3$ são os parâmetros da distribuição de $X$.

Visando gerar uma amostra da distribuição $f_X(x)$, considere a semente em $129$ de um processo iterativo com os seguintes passos:
1. Simule um valor $x$ de uma distribuição Normal com média $\mu = 4$ e variância $\sigma^2 = 2.25$, com função densidade de probabilidade denotada por $g(x)$.
2. Simule $u$ de uma distribuição uniforme no intervalo ($0,1$).
3. Se $u<\frac{f_X(x)}{2\,g(x)}$, aceite $x$. Caso contrário, rejeite-o.
4. Repita os passos anteriores até obter uma amostra de dimensão $n = 522$.
Indique a proporção observada de semanas em que o volume de gasolina vendida não ultrapasse $25$ mil litros, arredondada a 4 casas decimais. Selecione a resposta correta de entre as seguintes opções:
1. 0.0041
2. 0.8014
3. 0.1143
4. 0.3305
5. 0.1494
Resolução
```
# Por exemplo:
# seed = 123 
# a = 1
# b = 2
# n = 501
# media = 3
# variancia = 1
# valor = 2

set.seed(seed)

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

amostra <- vector(length = n)
aceites <- 0

while(aceites < n) {
  x <- rnorm(1, mean = media, sd = sqrt(variancia))
  u <- runif(1)
  if(u < (f(x)/(2*dnorm(x, mean = media, sd = sqrt(variancia)))) ) {
    aceites <- aceites + 1
    amostra[aceites] <- x
  }
}

solution <- length(amostra[amostra <= valor])/n
```
O valor pedido é igual a 0.1494.
Usando o R, fixe a semente em $1041$ e, com recurso ao pacote MASS, gere $n=310$ pares de observações provenientes de uma população normal bivariada, $(X,Y)$, de valores médios $\mu_X=2.1$, $\mu_Y=3$, desvios padrão $\sigma_X=4.7$, $\sigma_Y=1.8$ e coeficiente de correlação $\rho=-0.1$.

A matriz de covariância de $(X,Y)$ é expressa por: \[ \Sigma = \begin{pmatrix} \operatorname{Var}(X) & \operatorname{Cov}(X, Y) \\ \operatorname{Cov}(Y, X) & \operatorname{Var}(Y) \end{pmatrix} \]

Calcule empiricamente as probabilidades conjuntas $P(X>\mu_X+\sigma_X, Y>\mu_Y+\sigma_Y)$ e $P(X<\mu_X+\sigma_X, Y<\mu_Y+\sigma_Y)$ isto é, determine a proporção de observações da amostra que satisfazem as respetivas condições.

Indique o valor da soma das probabilidades empíricas obtidas anteriormente, arredondada a 4 casas decimais.

Resolução
```
# seed <- 1234
# n<-100
# muX<-2; sX<-2
# muY<-4; sY<-2
# rho<-0.6

set.seed(seed)

library(MASS)

mu <- c(muX,muY)
sigma <- matrix(c(sX^2,sX*sY*rho,sX*sY*rho,sY^2),2)

bivsample <- mvrnorm(n,mu,sigma)

solution <- sum(bivsample[,1]>muX+sX & bivsample[,2]>muY+sY)/n +
            sum(bivsample[,1]<muX+sX & bivsample[,2]<muY+sY)/n
```
O valor pedido é 0.7387.
Considere variáveis aleatórias $X_1, X_2, \ldots, X_n$ independentes e identicamente distribuídas, cada uma seguindo uma distribuição exponencial com parâmetro $\lambda = 19.56$.

Considere a soma ponderada:

\[ S_n = \sum_{i=1}^{n} w_i X_i \] onde $w_i = 1 + \frac{i-1}{n}$ para $i = 1, \ldots, n$.
1. Considere $n=90$ e calcule aproximadamente, usando o Teorema do Limite Central o quantil de probabilidade $\chi_p$, para $p=0.64$. Note que, neste caso, sob certas condições de regularidade, o Teorema do Limite Central continua a ser aplicável, conservando as propriedades assintóticas usuais. Ou seja, \[ \frac{S_n - \mathbb{E}[S_n]}{\sqrt{\operatorname{Var}(S_n)}} \xrightarrow{d} N(0, 1) \]
2. Fixe a semente em $908$ e realize uma simulação com $m =400$ réplicas (amostras). Para cada réplica, gere os $n=90$ valores de $X_i$, calcule a soma ponderada $S_n$ e, no final, estime empiricamente o quantil de probabilidade $\chi_p$, para $p=0.64$, através dos resultados simulados (usando a função quantile, com a opção type=4).
Obtenha o valor absoluto da diferença entre os valores que obteve em 1. e em 2., arredondado a 4 casas decimais e selecione a resposta correta de entre as seguintes opções:
1. 0.0627
2. 0.0585
3. 0.0630
4. 0.0572
5. 0.0616
Resolução
```
# # Por exemplo:
# seed <- 666
# lb <- 5.56
# n <- 100
# m <- 250
# pq <- 0.13

w <- 1 + ((1:n)-1)/n

media_X <- 1/lb
var_X <- 1/lb^2
media_Sn <- sum(w)*media_X
var_Sn <- sum(w^2)*var_X
chi_p_TLC <- qnorm(pq, mean=media_Sn, sd=sqrt(var_Sn))

set.seed(seed)
Sn_sim <- replicate(m, {
  X <- rexp(n, rate=lb)
  sum(w*X)
})

chi_p_sim <- quantile(Sn_sim, probs=pq, type=4)

solution <- abs(chi_p_TLC - chi_p_sim)
```
A diferença pedida é igual a 0.0616.
Numa experiência com 2020 pessoas, classificadas quanto ao daltonismo e ao sexo, obtiveram‑se os seguintes resultados:

Masculino e normal Masculino e daltónico Feminino e normal Feminino e daltónico

844 116 1028 32

Admita que esta situação pode ser descrita por uma experiência aleatória composta por $n=2020$ ensaios independentes, em cada um dos quais pode ocorrer um de quatro acontecimentos mutuamente exclusivos. Denote por $p_i$ a probabilidade de ocorrência do acontecimento $i$ em cada ensaio, com $\sum_{i=1}^4 p_i = 1$.

Seja $X_i$ a variável aleatória que representa o número de ensaios em que o acontecimento $i$ ocorre ($i\!=\!1,2,3,4$), de modo que $\sum_{i=1}^4 X_i = n$.

Os quatro acontecimentos correspondem a: $i=1$: Masculino e normal; $i=2$: Masculino e daltónico; $i=3$: Feminino e normal e $i=4$: Feminino e daltónico.

O vetor aleatório $(X_1,X_2,X_3,X_4)$ segue uma distribuição Multinomial com função (massa) de probabilidade dada por \[ f_{X_1,X_2,X_3,X_4}(x_1,x_2,x_3,x_4|p_1,p_2,p_3,p_4) \ = \ g(x_1,x_2,x_3,x_4) \times \prod_{i=1}^4 p_i^{x_i} \] De acordo com o modelo genético, as probabilidades dos quatro acontecimentos são dadas por: $p_1=\frac{p}{2}$, $p_2=\frac{(1-p)}{2}$, $p_3=p-\frac{p^2}{2}$ e $p_4=\frac{(1-p)^2}{2}$, com $0<p<1$.

Com base nos dados observados, considere a função de verosimilhança associada à distribuição Multinomial e determine a estimativa de máxima verosimilhança (EMV) do parâmetro $p$, maximizando o logaritmo da função de verosimilhança.

Apresente o valor estimado de $p$, arredondado a quatro casas decimais, selecionando a opção correta de entre as respostas abaixo:
1. 0.1438
2. 0.8591
3. 0.9377
4. 0.9827
5. 0.1576
Resolução

\[ \begin{align*} \log L(p|x_1,\ldots,x_4) & = x_1 \log\left(\frac{p}{2}\right) + x_2\log\left(\frac{1-p}{2}\right) + x_3 \log\left(\frac{p-p^2}{2}\right) + \\ &+ x_4 \log\left(\frac{(1-p)^2}{2}\right) + \log g(x_1,x_2,x_3,x_4) \end{align*} \]
```
# Por exemplo:
# x <- c(884, 76, 1028, 12)

log_vero <- function(p) {
  x[1] * log(p/2) + x[2] * log((1-p)/2) + x[3] * log(p-p^2/2) + x[4] * log(((1-p)^2)/2)
}

emv <- optimize(log_vero, c(0, 1), maximum = TRUE)

solution <- emv$maximum
```
O valor pedido é igual a 0.8591.
Pretende-se construir intervalos de confiança para a diferença entre os valores médios de duas populações independentes, $X_1$ e $X_2$, utilizando o Teorema do Limite Central aplicado a cada uma das médias amostrais, assumindo que as variâncias populacionais são conhecidas.

Para tal deve realizar uma simulação onde, fixando inicialmente a semente do gerador de números aleatórios em $2107$, se repetem $N = 420$ vezes os seguintes passos:
- Gere uma amostra de dimensão $143$ proveniente de uma distribuição exponencial com parâmetro $1.8$.
- Gere uma amostra de dimensão $277$ proveniente de uma distribuição exponencial com parâmetro $3$.
- Construa um intervalo de confiança aproximado para a diferença entre os valores médios, $(\mathbb{E}[X_1] - \mathbb{E}[X_2])$, com nível de confiança $98\%$.
Após as $N=420$ repetições, calcule a proporção de intervalos obtidos que não contêm a diferença entre os valores médios. Apresente esse valor, arredondado a 4 casas decimais, selecionando a opção correta de entre as respostas abaixo:
1. 0.9086
2. 0.0015
3. 0.0310
4. 0.0020
5. 0.0136
Resolução
```
# seed <- 1234
# ns1<-200
# lambda1<-1
# ns2<-200
# lambda2<-2
# alpha<-0.05
# N<-300

set.seed(seed)

l <- 0
for (i in 1:N) {
  sample1 <- rexp(ns1, lambda1)
  m1 <- mean(sample1)
  s1 <- 1/(lambda1^2)
  sample2 <- rexp(ns2, lambda2)
  m2 <- mean(sample2)
  s2 <- 1/(lambda2^2)
  diff_hat <- m1 - m2
  semi_ampl <- qnorm(1-alpha/2) * sqrt(s1/ns1 + s2/ns2)
  ic_lwr <- diff_hat - semi_ampl
  ic_upr <- diff_hat + semi_ampl
  true_diff <- 1/lambda1 - 1/lambda2
  l <- l + !(ic_lwr <= true_diff & true_diff <= ic_upr)
}

solution <- l/N
```
O valor pedido é igual a 0.0310.
Para a construção de testes de hipóteses para o parâmetro $\lambda$ de uma distribuição de Poisson podemos recorrer à variável:

\[ T = \frac{\overline{X} - \lambda}{\sqrt{\lambda/n}} \overset{a}{\sim} N(0,1) \]

obtida pela aplicação do Teorema do Limite Central a uma amostra aleatória de tamanho $n$ suficientemente grande da referida distribuição.

Para testar a hipótese $H_0: \lambda=7.8$ contra $H_1: \lambda = 7.7$, recorre-se à estatística de teste $T_0$, obtida de $T$ admitindo que a hipótese $H_0$ é verdadeira.
1. Considere $n= 60$ e obtenha um valor aproximado da probabilidade de rejeição de $H_0$ sabendo que $\lambda = 7.7$, para um nível de significância $\alpha = 0.025$.
2. Obtenha um valor simulado de tal probabilidade, procedendo do seguinte modo:
  1. Fixe a semente em $402$, gere $m=430$ amostras de dimensão $n = 60$ da distribuição de Poisson com $\lambda = 7.7$.
  2. Calcule a frequência relativa de rejeições de $H_0$ a que o teste conduz sabendo que $\lambda = 7.7$.
Calcule o valor absoluto da diferença entre os valores que obteve em 2. e em 1., arredondado a 4 casas decimais e selecione a resposta correta de entre as seguintes opções:
1. 0.0146
2. 0.0149
3. 0.0145
4. 0.0159
5. 0.0165
Resolução
```
# # Por exemplo:
# seed <- 1234
# alpha <- 0.05
# lambda0 <- 5.4
# lambda1 <- 5.9
# n <- 50
# m <- 700

k <- qnorm(alpha)
vcritico <- k * sqrt(lambda0 / n) + lambda0
Pteorica <- pnorm((vcritico - lambda1) / sqrt(lambda1 / n))

set.seed(seed)
rejeita <- replicate(m, {
  x <- rpois(n, lambda1)
  T0 <- (mean(x) - lambda0) / sqrt(lambda0/n)
  T0 < k
})

Psimulada <- mean(rejeita)

solution <-abs(Psimulada - Pteorica)
```
O valor pedido é igual a 0.0149.
Um programa de computador deve gerar números pseudo-aleatórios correspondentes a realizações de uma variável aleatória $X$ cuja função densidade de probabilidade é dada por:

\[ f_X(x) = \begin{cases} 1 - \frac{x}{2}, &\ 0 \leq x \leq 2, \\ 0, &\text{caso contrário.} \end{cases} \]

Foram gerados e agrupados $217$ números pseudo-aleatórios, obtendo-se a seguinte tabela.

classe $[0.0,0.4[$ $[0.4,0.8[$ $[0.8,1.2[$ $[1.2,1.6[$ $[1.6,2.0]$

frequência 80 58 41 23 15

Indique o valor-p, arredondado a 4 casas decimais, correspondente ao teste estatístico adequado para verificar se o programa de computador cumpre a sua finalidade.

O valor pedido é igual a 0.2617.

Probabilidade e Estatística

Projeto computacional

Resolução

Resolução

Resolução

Resolução

Resolução

Resolução

Resolução

Resolução

Resolução