Probabilidade e Estatística

Projeto computacional

O ficheiro diabetes.txt contém um conjunto de dados de um estudo sobre a progressão do diabetes.

Comece por classificar cada observação de acordo com a sua pertença a um dos três intervalos definidos pelos quantis de ordens 0.25 e 0.75 da variável BP.

Com recurso ao pacote ggplot2, produza um único diagrama de dispersão que represente os valores da variável S6 em função da variável RISK conjuntamente com a classificação atrás referida.

Tendo em conta que o texto no ficheiro de dados se encontra em Inglês, mantenha todo o texto do gráfico em Inglês.

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

Todo o código em R, incluindo os comandos para leitura dos dados do ficheiro.
O gráfico produzido.

Resolução

library(ggplot2)
theme_set(theme_grey())

# url <- 'https://web.tecnico.ulisboa.pt/paulo.soares/pe/projeto/diabetes.txt'
# dados <- read.delim(url)
# # Por exemplo:
# # RISK vs BMI, classificação por AGE
# quantis <- c(0.25, 0.75)
# vars <- c("BMI", "RISK", "AGE")

limites <- quantile(dados[[vars[3]]], probs = c(0, quantis, 1))
dados[[vars[3]]] <- cut(dados[[vars[3]]], breaks = limites, include.lowest = TRUE)
# Para uma variável fixa use:
# limites <- quantile(dados$AGE, probs = c(0, quantis, 1))
# dados$AGE <- cut(dados$AGE, breaks = limites, include.lowest = TRUE)

ggplot(dados) +
  geom_point(aes(.data[[vars[1]]], .data[[vars[2]]], colour = .data[[vars[3]]])) +
  # Para variáveis fixas use:
  # geom_point(aes(BMI, RISK, colour = AGE)) +
  labs(title = "Diabetes progression study")

O ficheiro bvl.csv contém um conjunto de dados sobre as transações diárias entre outubro de 2022 e agosto de 2024 de ações de 10 empresas cotadas na Bolsa de Valores de Lisboa (Fonte: EURONEXT).

Após a leitura desse ficheiro no R, selecione todos os dados referentes à empresa GALP ENERGIA-NOM do dia 2022-10-01 até ao dia 2022-11-30.

Com recurso ao pacote ggplot2, produza um único gráfico que ilustre a variação da cotação da empresa, em €, no fecho de cada sessão diária (variável Close) ao longo do período de tempo referido. Adicionalmente, pretende-se enquadrar essa cotação pelas cotações mínima e máxima registadas em cada dia (variáveis Low e High). Todo o texto do gráfico deverá estar em Português.

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

# url <- 'https://web.tecnico.ulisboa.pt/paulo.soares/pe/projeto/bvl.csv'
# dados <- read.csv(url)
# # Por exemplo:
# inicio <- as.Date("2023-09-01")
# fim <- as.Date("2023-11-30")
# empresa <- "EDP"

dados$Date <- as.Date(dados$Date)

subcjt <- subset(dados, Name == empresa & Date >= inicio & Date <= fim)

ggplot(subcjt) +
  geom_ribbon(aes(Date, ymin = Low, ymax = High), fill = "limegreen",
              alpha = 0.25) +
  geom_line(aes(Date, Close), color = "tomato", linewidth = 0.25) +
  geom_point(aes(Date, Close), color = "tomato", size = 0.5) +
  labs(title = paste("Cotação diária da", empresa, "no fecho da BVL\nentre",
                     inicio, "e", fim),
       x = "Data", y = "Cotação no fecho (mínima e máxima)  em €")
```
O ficheiro Cereals Production, Area & Yield_Usable.xlsx contém dados coletados pela Comissão Europeia sobre a produção anual de cereais desde 1993 até 2021 nos vários países da EU-27. Os dados incluem a área cultivada (area em $10^3$ hectares), a produção útil (Usable production em $10^3$ toneladas) e a produtividade (Yield em toneladas por hectare) para as culturas de vários cereais.

Recorrendo ao pacote ggplot2, produza um único gráfico de barras empilhadas que permita comparar a distribuição da variável Yield entre as diferentes culturas de cereais nos anos de 1995, 2007 e 2020 em Spain (em Inglês no ficheiro).

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

# url <- 'https://web.tecnico.ulisboa.pt/paulo.soares/pe/projeto/Cereals.xlsx'
# download.file(url, destfile = "local.xlsx")
# cereals <- readxl::read_xlsx("local.xlsx")
# # Por exemplo:
#pais <- "Portugal"
#anos <- c(1993, 2007, 2020)
variavel <- "area"

subcjt <- subset(cereals, `Member State` == pais & Type == variavel &
                 Product != "Total cereals")

final <- data.frame()
for(ano in anos) {
  temp <- subcjt[, c("Product", as.character(ano))]
  names(temp) <- c("Product", "Area")
  temp$Year <- ano
  final <- rbind(final, temp)
}
final$Year <- as.factor(final$Year)

ggplot(final) +
  geom_col(aes(x = Year, y = Area, fill = Product)) +
  labs(title = paste("Cereals production in", pais),
       y = "Area (1000 hectares)")
```
Considere o seguinte jogo que envolve duas moedas distintas e em que, em cada lançamento de qualquer moeda, o jogador aumenta uma unidade ao seu saldo se obtiver uma “cara” e diminui uma unidade no caso contrário. As moedas têm as seguintes características:

Moeda 1: P(“cara”) = $p_1 = 0.091$

Moeda 2: P(“cara”) = $p_2 = 0.741$

Se o saldo do jogador antes de uma jogada for um múltiplo de 3, ele lança a moeda 1. Caso contrário, lança a moeda 2. Considere que o jogador começa com um saldo nulo e que um saldo negativo não impede o jogo de prosseguir.

Fixando a semente em $1516$, simule $n=1600$ jogadas e indique o saldo final do jogador.

Resolução
```
# # Por exemplo:
# seed <- 1234
# p1 <- 0.095
# p2 <- 0.745
# n <- 2000

jogo <- function(saldo){
  prob.ganhar <- ifelse(saldo %% 3 == 0, p1, p2)
  tmp <- ifelse(rbinom(1, size = 1, prob = prob.ganhar) == 1, 1, -1)
  return(tmp)
}

saldo <- function(k){
  saldo <- 0
  for(i in 1:k)
    saldo <- saldo + jogo(saldo)
  return(saldo)
}

set.seed(seed)

solution <- saldo(n)
```
O valor pedido é -34.
Admita que as variáveis aleatórias não negativas $X$ e $Y$ representam os níveis de cheia de certo rio e da intensidade das chuvas. Admita também que $~ P(X > x, Y > y) = e^{- x - y - \theta xy}, ~$ para $x, y \geq 0,~$ onde $~\theta~$ é uma constante no intervalo $[0,1]$.

O coeficiente de correlação entre $X$ e $Y$ é dado por $~\rho = - 1 + \int_0^{\infty } h(z) \, dz,~$ onde $~h(z) = \frac{e^{-z}}{\theta z + 1}$.

Para obter valores aproximados de $~\rho$ quando $~\theta= 0.4$, implemente os dois procedimentos seguintes.

A. Coeficiente de correlação amostral

A concretização de uma amostra aleatória de dimensão $n=15$ proveniente de $(X, Y)$ com $~\theta= 0.4~$ conduziu ao seguinte conjunto de observações.

x: 0.163, 1.845, 1.012, 1.246, 0.163, 0.026, 5.467, 3.808, 1.607, 2.006, 0.859, 2.156, 0.244, 0.462, 0.033

y: 0.6, 2.216, 0.439, 1.705, 3.524, 0.112, 0.759, 0.95, 0.066, 0.011, 1.866, 1.041, 5.351, 0.28, 0.064

Tire partido destas observações e obtenha o coeficiente de correlação amostral.

B. Integração de Monte Carlo
1. Fixando a semente em 1796, gere $m = 49500$ números pseudoaleatórios $(u)$ da distribuição exponencial com parâmetro unitário.
2. Avalie a função $g(u) = \frac{1}{\theta u + 1}$, onde $\theta = 0.4$, nos $m = 49500$ pontos aleatórios $(u)$ e calcule a média desses valores. Ao resultado subtraia uma unidade.
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.5
# data <- data.frame(x = c(3.257, 0.439, 0.376, 0.019, 0.456, 0.126, 1.219, 1.121,
#                          0.547, 0.13, 0.848, 0.442),
#                    y = c(1.418, 0.124, 1.24, 2.613, 3.341, 0.086, 0.004, 2.45,
#                          1.524, 1.092, 0.327, 2.85))

# A. Coeficiente de correlação amostral

corr <- cor(data)
samplecorr <- corr[1,2]

# B. Integração de Monte Carlo

set.seed(seed)

g <- function(z) {
  1 / (theta * z + 1)
}

u_vals <- rexp(m, rate = 1)
f_vals <- g(u_vals)
aproxMCcorr <- mean(f_vals) - 1

solution <- abs(samplecorr - aproxMCcorr)
```
O desvio absoluto pretendido é 0.048667.
Sejam:
1. $\{ X_1, X_2, \dots\}$ uma sequência de variáveis aleatórias independentes e identicamente distribuídas a $X$, com valor esperado comum nulo, variância positiva finita $\sigma^2$ e terceiro momento absoluto finito $\nu_3 = E(|X|^3)$;
2. $F_n(x)$ a função de distribuição (f. d.) exacta de $\left(\sum_{i=1}^n X_i\right)/\sqrt{n \sigma^2}$;
3. $\Phi(x)$ a f. d. da distribuição normal padrão, que aproxima $F_n(x)$ de acordo com o teorema do limite central (TLC).
O teorema de Berry-Esseen estabelece que existe uma constante $C > 0$ tal que \[\left| F_n(x) - \Phi(x) \right| \leq C \, \frac{\nu_3 / \sigma^3}{\sqrt{n}} = \xi, \quad \forall \, x \in \mathbb{R}, \, n \in \mathbb{N},\] isto é, o desvio absoluto entre as f. d. $F_n(x)$ e $\Phi(x)$ nunca excederá $\xi$, sendo esta a velocidade de convergência do TLC. (A literatura recente estipula que $C < 0.4748$)

Considere que $X$ possui distribuição uniforme no intervalo $(-a, a)$ com $a = 1.63$, $C=0.4748$ e $n=12$. Obtenha o valor de $\xi$ arredondado a 4 casas decimais.

Resolução
```
# # Por exemplo:
# a <- 1.5
# n <- 18

C <- 0.4748

f <- function(x) {
  return(abs(x)^3 / (2 * a))
}

nu3 <- integrate(f, -a, a)$value
sigma <- sqrt(a^2 / 3)

solution <- C * nu3 / sigma^3 / sqrt(n)
```
O valor pretendido é 0.1781.
Seja $(X_1, \dots, X_n)$ uma amostra aleatória de dimensão $n$ proveniente de uma população normal com valor esperado nulo e desvio padrão desconhecido $\sigma$ $(\sigma > 0)$.

O estimador de $\sigma$, $S = \sqrt{\frac{1}{n-1} \sum_{i=1}^n (X_i - \bar{X})^2}$, possui valor esperado igual a $~E(S) = \sqrt{\frac{2}{n-1}} \, \frac{\Gamma(n/2)}{\Gamma [ (n-1)/2 ]} \times \sigma$, onde $\Gamma (x), x>0$ representa a função gamma definida no R.

Importa notar que $E(S) \neq \sigma, \, \forall \sigma > 0$, pelo que $S$ não é um estimador centrado de $\sigma$, apesar de $S^2$ ser um estimador centrado de $\sigma^2$.
1. Obtenha o valor exacto de $E(S)$, quando $\sigma = 1.4$ e $n =14$, e constate que é inferior a $\sigma$.
2. Estime o valor de $E(S)$, fixando a semente em 1082 e procedendo do seguinte modo:
  1. gere $m = 35$ amostras de dimensão $n =14$ da distribuição normal com valor esperado nulo e desvio padrão $\sigma = 1.4$;
  2. calcule o valor do desvio padrão amostral de cada uma das amostras geradas;
  3. obtenha a média dos $m = 35$ valores obtidos em b.
Indique o desvio absoluto (arredondado a 6 casas decimais) entre o valor exacto obtido em 1. e a estimativa determinada em 2c.

Resolução
```
# # Por exemplo:
# seed <- 1234
# n <- 13
# sigma <- 1.6
# m <- 50

# Valor exacto
ES <- sqrt(2 / (n - 1)) * gamma(n / 2) / gamma((n - 1) / 2) * sigma
cat(ES)
```
```
## 1.37336
```
```
# Estimativa
set.seed(seed)

sampleES <- 0
for (i in 1:m) {
  data <- rnorm(n, mean = 0, sd = sigma)
  sampleES <- sampleES + sd(data)
}
sampleES <- sampleES / m

solution <- abs(ES - sampleES)
```
O desvio absoluto pretendido é 0.054490.
Na gestão de um projecto, a variável aleatória $X$ representa o tempo de execução de uma tarefa restrito ao intervalo $(0,1)$. Admita que $X$ possui função de densidade de probabilidade dada por \[f_X(x) = \alpha \, x^{\alpha -1 }, \quad x \in (0,1),\] onde $\alpha$ é um parâmetro positivo desconhecido.

A gestora do projecto planeia recolher uma amostra casual de $n$ tempos de execução e pretende deduzir dois intervalos de confiança $\gamma$ para $\alpha$, recorrendo à variável aleatória fulcral para $\alpha$ \[Z = - 2 \, \alpha \sum_{i=1}^n \ln(X_i) \sim \chi_{(2n)}^2.\]
1. Deduza um intervalo aleatório de confiança $\gamma$ para $\alpha$ baseado nos quantis de probabilidade \[\begin{cases} a=F_{\chi_{(2n)}^2}^{-1}\left(\frac{1-\gamma}{2}\right) \\ b=F_{\chi_{(2n)}^2}^{-1}\left(\frac{1+\gamma}{2}\right). \end{cases}\]
2. Com vista à minimização da amplitude esperada do intervalo aleatório de confiança $\gamma$ para $\alpha$, um estatístico sugeriu à gestora que substituísse $(a,b)$ pelo par de quantis de probabilidade \[(c, d) : \begin{cases} f_{\chi_{(2n)}^2}(d) - f_{\chi_{(2n)}^2}(c) = 0\\ F_{\chi_{(2n)}^2}(d)-F_{\chi_{(2n)}^2}(c) = \gamma. \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.
3. Indique o rácio, arredondado às milésimas, das amplitudes esperadas dos intervalos aleatórios de confiança deduzidos em 1. e 2., para $n = 12$ e $\gamma = 0.95$.
Resolução
```
# # Por exemplo:
# n <- 6
# gama <- 0.95

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

library(pracma)

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

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

solution <- (b - a) / (raiz[2] - raiz[1])
```
O valor pretendido é 1.019.
Considere que a temperatura máxima diária em determinado local é representada pela variável aleatória $X$ com distribuição de Gumbel com parâmetro de localização (“location”) $\lambda$ desconhecido e parâmetro de escala (“scale”) $\delta$ unitário, i. e., \[F(x) = P(X \leq x) = e^{-e^{-\frac{x-\lambda}{\delta}}}, \quad x \in \mathbb{R},\] onde $\delta=1$.

Seja $(X_1, \ldots, X_n)$ uma amostra aleatória proveniente de $X$. Para testar a hipótese $H_0: \lambda=23$ contra $H_1: \lambda = 23.2$ recorre-se à estatística de teste \[T = \frac{- \ln \left( \frac{1}{n} \sum_{i=1}^n e^{-X_i} \right) - \lambda_0}{\sqrt{\frac{1.10866}{n}}} \stackrel{a}{\sim}_{H_0} \text{normal}(0,1)\] e rejeita-se $H_0$ caso $T > k = \Phi^{-1} (1-\alpha)$, com $\alpha = 0.05$.
1. Considere $n= 65$ e obtenha um valor aproximado da probabilidade de rejeição de $H_0$ sabendo que $\lambda = 23.2$.
2. Obtenha um valor simulado de tal probabilidade, procedendo do seguinte modo:
  1. Fixe a semente em $3034$ e recorra ao pacote extraDistr para gerar $m=900$ amostras de dimensão $n = 65$ da distribuição de Gumbel com $\lambda = 23.2$ e $\delta=1$.
  2. Calcule a frequência relativa de rejeições de $H_0$ a que o teste conduz sabendo que $\lambda = 23.2$.
Calcule o valor aproximado do quociente dos resultados obtidos em 2. e 1. e selecione a resposta correta de entre as seguintes opções:
1. 0.719
2. 1.086
3. 1.171
4. 1.021
5. 0.806
Resolução
```
# # Por exemplo:
# seed <- 1234
# alpha <- 0.05
# lambda0 <- 29
# lambda1 <- 29.25
# n <- 50
# m <- 700

k <- qnorm(1 - alpha) 
paprox <- 1 - pnorm(k + (lambda0 - lambda1) / sqrt(1.10866 / n))

library(extraDistr)
set.seed(seed)

freq_rejeicao <- 0
for(i in 1:m) {
  data <- rgumbel(n, mu = lambda1, sigma = 1)
  t <- (-log(mean(exp(-data))) - lambda0) / sqrt(1.10866 / n)
  freq_rejeicao <- freq_rejeicao + (t > k)
}
prop_rejeicao <- freq_rejeicao / m

solution <- prop_rejeicao / paprox
```
O quociente pedido é igual a 1.021.
Um engenheiro civil obteve o seguinte conjunto de $n= 120$ observações ordenadas do recobrimento de armaduras ($X$, em mm) de determinado elemento estrutural de uma ponte.

21.6, 31.9, 41, 54.4, 55.8, 30.3, 10.8, 28.5, 45.3, 41.6, 10.5, 16.5, 54.8, 57.2, 40.1, 38.5, 54.4, 48.4, 51.7, 53.7, 41.8, 52.4, 65.5, 25.9, 76, 49.3, 57.6, 34.2, 48, 57.7, 42.3, 35.6, 21.9, 56.2, 35.5, 37.4, 45.3, 66.6, 42.7, 67.7, 32.1, 53, 52.5, 45.3, 31.4, 43.7, 17.4, 38, 21.6, 50.2, 43.2, 30, 48.2, 38.2, 20.7, 19, 27, 40.9, 56, 75.6, 57.1, 72.2, 24.6, 55.4, 53.2, 49.4, 40.6, 33.6, 56.2, 44.3, 34.3, 52.5, 47.8, 46, 39.4, 50.4, 52.6, 41.2, 35.5, 23, 26.4, 50.8, 37.5, 51, 37.5, 41.4, 56.8, 29.2, 34.6, 44.4, 52.8, 32.9, 32.7, 53, 38, 32.7, 26.2, 34.3, 22.6, 39.8, 35.4, 25, 43.1, 37.9, 55.6, 49.3, 57.7, 49.2, 39.8, 25.3, 46, 14.2, 5.5, 49.6, 54.4, 36, 36.6, 24.1, 37.9, 46.7

Ele conjectura que $X$ possua distribuição pertencente ao modelo $\text{Weibull}(\theta, \beta)$ com função de distribuição \[F_0(x) = 1 - \exp \left[ -\left( \frac{x}{\beta} \right)^\theta \right], \quad x>0,\] onde o parâmetro de forma (“shape”) é igual a $\theta= 1.85$ e o parâmetro de escala (“scale”) $\beta$ é uma constante positiva desconhecida.

Tendo em vista auxiliar o engenheiro a testar a hipótese $H_0: X \sim \text{Weibull}(\theta, \beta)$ contra $H_1: X \not \sim \text{Weibull}(\theta, \beta)$, proceda do seguinte modo.
1. Obtenha a estimativa de máxima verosimilhança de $\beta$, $\hat{\beta} = \left ( \frac{1}{n} \sum_{i=1}^n x_i^\theta \right)^{\frac{1}{\theta}}$, com base nas observações acima.
2. Divida o suporte da variável aleatória $X$, $\mathbb{R}^+$, em $k= 5$ classes equiprováveis sob a hipótese de as observações serem provenientes da distribuição $\text{Weibull}(\theta= 1.85, \beta= \hat{\beta})$.
3. Agrupe as abservações nas classes definidas em 2., obtendo o conjunto de frequências absolutas observadas, $o_i$ $(i=1, \dots, k)$, associadas a essas classes.
4. Recorra às frequências obtidas em 3. e calcule $t = \frac{k}{n} \left( \sum_{i=1}^k o_i^2 \right) - n$, o valor observado da estatística do teste de ajustamento do qui-quadrado para classes equiprováveis.
5. Tenha presente que . . .
  1. Não deverá rejeitar $H_0: X \sim \text{Weibull}(\theta, \beta)$ ao nível de significância (n. s.) $\alpha_0$, se $t < F_{\chi_{(k-2)}^2}(1-\alpha_0)$.
  2. Deverá rejeitar $H_0$ ao n. s. $\alpha_0$, se $t > F_{\chi_{(k-1)}^2}(1-\alpha_0)$.
  3. Deverá considerar o teste de ajustamento inconclusivo ao n. s. $\alpha_0$, se $F_{\chi_{(k-2)}^2}(1-\alpha_0) \leq t \leq F_{\chi_{(k-1)}^2}(1-\alpha_0)$.
Com base neste procedimento, indique qual das cinco decisões abaixo deverá tomar o engenheiro civil.
1. Rejeitar $H_0$ aos n. s. de 1% e 5%.
2. Rejeitar $H_0$ ao n. s. de 5% e não rejeitar $H_0$ ao n. s. de 1%.
3. Teste é inconclusivo ao n. s. de 1%.
4. Não rejeitar $H_0$ aos n. s. de 5% e 1%.
5. Teste inconclusivo ao n. s. de 5%.
Resolução
```
# # Por exemplo:
# dados <- c(38.8, 45.5, 32.4, 37.3, 32.3, 46.6, 67.8, 21.5, 27.7, 36.6, 84.8, 15.4,
#            49.9, 53.4, 30.8, 71.7, 24.4, 34.8, 54.6, 68.7, 53.2, 40.1, 57.7, 33,
#            71.2, 41.8, 52.2, 60.2, 14.8, 18.8, 25.9, 26.6, 30.9, 27.4, 51, 39.9,
#            28.9, 25.4, 80.4, 47.1, 20.1, 9.2, 59.9, 28.1, 57.8, 22.7, 36.5, 45.8,
#            44.5, 63, 28.4, 20.4, 29.4, 49.6, 62.4, 58.5, 58.4, 22.1, 33.8, 28.7,
#            17, 11.9, 28.2, 32.7, 37.6, 45.3, 20.9, 93.6, 37.8, 11.4, 50.5, 74.3,
#            52.9, 3.3, 68.7, 36.2, 58.5, 26.9, 56, 53.3, 59.3, 31.2, 38.1, 15.8,
#            37.8, 17.6, 55.1, 65.7, 52.3, 47.9, 49.6, 30.7, 54.2, 51.2, 31.5, 40.8,
#            45.9, 37.2, 10.3, 39.1, 26.8, 38, 48.7, 55.2, 34.4, 69.7, 72.4, 29.3,
#            39, 47.8, 38.6, 24.5, 37.1, 45.9, 43.2, 52.3, 67.6, 45.9, 36.9, 35.1,
#            44.1, 40.5, 72.6, 18.5, 35.5, 30.3, 14.7, 45, 27.9, 26.3, 61.2, 36.6,
#            42.8, 42.6, 37.8, 25.7, 40.6, 39.2, 51.8, 59.2, 40.5, 24.3, 15.3, 28.3,
#            38.4, 18.8, 93.7, 27, 37.9, 56.8, 45.1, 38.5, 59.7, 62.8, 48.3, 24.9,
#            26.8, 49.4, 61.2, 16.5, 38.9, 30.1, 38.1, 37.2, 30.6, 62.1, 51.4, 36.9,
#            28.4, 27.9, 24.8, 40.3, 36.8, 52.1, 14.4, 26.4, 24.4, 61.4, 7.2, 11.5)
# n <- 100
# theta <- 2
# k <- 5

betahat <- (1 / n * sum(dados^theta))^{1 / theta}

probs <- 0:k / k
limites <- qweibull(probs, shape = theta, scale = betahat)

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

t <- k / n * sum(freq_observadas^2) - n

a <- qchisq(0.95, df = k - 2)
b <- qchisq(0.95, df = k - 1)
c <- qchisq(0.99, df = k - 2)
d <- qchisq(0.99, df = k - 1)
```
```
cat("t =",t, "\n Para α=0.01: ", c, d, "\n Para α=0.05: ", a, b)
```
```
## t = 44.41667 
##  Para α=0.01:  11.34487 13.2767 
##  Para α=0.05:  7.814728 9.487729
```
A decisão a tomar é: Rejeitar $H_0$ aos n. s. de 1% e 5%.