Probabilidade e Estatística

Projeto computacional

  1. O ficheiro hcvdat0.csv (45.1KB) contém um conjunto de dados baseado no estudo clínico descrito em:

    Lichtinghagen, Ralf, Frank Klawonn, and Georg Hoffmann. 2020. HCV data. UCI Machine Learning Repository. (https://doi.org/10.24432/C5D612)

    Neste estudo, analisaram-se 615 observações de dadores de sangue e pacientes com patologias hepáticas para compreender como biomarcadores bioquímicos se relacionam com o diagnóstico clínico.

    A variável Category representa o diagnóstico clínico dos indivíduos, permitindo comparar dadores de sangue saudáveis (Blood Donors) com pacientes em diferentes estádios de progressão da Hepatite C, Fibrose ou Cirrose. As restantes variáveis referem-se a concentrações de enzimas e proteínas medidas no soro sanguíneo, em diferentes unidades.

    Com recurso ao pacote ggplot2 produza um gráfico que represente, através de diagramas de caixa (boxplots) paralelos, como é que a variável CHOL se relaciona com a variável Category.

    Submeta um ficheiro em formato PDF com uma única página A4, mantendo todo o texto da figura em Inglês e incluindo:


    Resolução

    library(ggplot2)
    theme_set(theme_light())
    
    # # Por exemplo:
    # url <- 'https://web.tecnico.ulisboa.pt/paulo.soares/pe/projeto/hcvdat0.csv'
    # dados <- read.csv(url)
    # variavel <- "PROT"
    
    ggplot(dados) +
      geom_boxplot(aes(x = Category, y = .data[[variavel]]), fill = "lightblue1") +
      labs(
        title = paste("Clinical distribution of", variavel, "by diagnosis category"),
        y = paste(variavel, "concentration")
      ) +
      theme(axis.text.x = element_text(angle = 25, hjust = 1))


  2. O ficheiro Petroleum&OtherLiquidFuels.txt (173.9 KB) contém valores diários, observados entre 2020-01-02 e 2026-03-02, referentes à variação relativa anual (YoY variation) dos preços do petróleo bruto e de alguns dos seus derivados. Para cada um dos produtos (desde a coluna “CrudeWTI” até à coluna “MBTXPropane”), o valor registado em cada dia corresponde à taxa de variação relativa anual do preço face ao mesmo dia do ano anterior. Além disso, o ficheiro contém uma coluna extra, “MEAN”, que apresenta, para cada dia, a média das variações relativas de todos os produtos listados (Fonte: U. S. Energy Information Administration, 2026).

    Com recurso ao pacote ggplot2, produza um cronograma que represente a evolução temporal da variação relativa dos preços para LADiesel e NYHeatingOil, comparando também com a série dos valores médios (MEAN). 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

    library(ggplot2)
    theme_set(theme_minimal())
    
    # # Por exemplo:
    # server <- 'https://web.tecnico.ulisboa.pt/paulo.soares/pe/projeto/'
    # file <- 'Petroleum&OtherLiquidFuels.txt'
    # url <- paste0(server, file)
    # dados <-  read.table(url, dec=",", sep="\t", header=T)
    # variaveis <- c("NYGasoline", "CrudeBrent")
    
    dados$Date <- as.Date(dados$Date) 
    
    ggplot(dados, aes(x = Date)) +
      geom_line(aes(y = MEAN, color = "Mean (All)"), alpha = 0.75) +
      geom_line(aes(y = .data[[variaveis[1]]], color = variaveis[1]), alpha = 0.75) +
      geom_line(aes(y = .data[[variaveis[2]]], color = variaveis[2]), alpha = 0.75) +
      labs(
        title = 'Prices of Petroleum & Other Liquid Fuels',
        y = "YoY variation",
        color = "Product"
      )


  3. O ficheiro Jobs_and_Growth.csv (1.4 MB) contém dados coletados pela Comissão Europeia sobre o emprego no setor primário da economia desde 2008 até 2022 nos vários países da UE-27.

    Recorrendo ao pacote ggplot2, produza um único gráfico de barras que permita comparar a evolução ao longo do tempo da variável Employment in the food industry entre os países Italy e Lithuania (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_bw())
    
    # url <- 'https://web.tecnico.ulisboa.pt/paulo.soares/pe/projeto/Jobs_and_Growth.csv'
    # dados <- read.csv(url)
    # # Por exemplo:
    # paises <- c("Portugal", "Spain")
    # variavel <- "Employment in agriculture"
    
    subcjt <- subset(dados, Member.State.Name %in% paises &
                            Sub.indicator.Name == variavel)
    
    ggplot(subcjt) +
      geom_col(aes(Time.Period, Data, fill = Member.State.Name), position = "dodge") +
      labs(x = "Year", y = "% of total employment",
           title = paste(variavel, "(2008-2022)"), fill = "") +
      theme(legend.position = "bottom")


  4. A distribuição de Gumbel é frequentemente utilizada para modelar eventos extremos em engenharia hidráulica como, por exemplo, os picos de caudal de rios. Uma variável aleatória contínua \(X\) segue uma distribuição de Gumbel se a sua função densidade de probabilidade for definida por:

    \[f_X(x) = \frac{1}{\sigma}\exp\left[-\left(\frac{x-\mu}{\sigma}+\exp\left[-\frac{x-\mu}{\sigma}\right]\right)\right],\, x\in \mathbb{R},\] onde \(\mu\in \mathbb{R}\) é o parâmetro de localização e \(\sigma \in \mathbb{R}^+\) é o parâmetro de escala.

    Para a gestão de risco de cheias, estabelecem-se habitualmente dois patamares críticos:

    Considere o caso específico do rio Grijalva (México), onde se admite que o caudal máximo diário num dado troço do rio, \(X\) (em \(m^3/s\)), segue a distribuição de Gumbel com parâmetros \(\mu= 2000\) e \(\sigma= 600\). Recorrendo ao pacote extraDistr determine:

    1. O valor exato de \(P(X>L_E∣X>L_A)\).

    2. Para uma amostra de dimensão \(n=1000\) de \(X\), gerada fixando a semente em 4789, a proporção de observações superiores a \(L_E\) no subconjunto dos valores que excederam \(L_A\).

    Indique o valor do quociente entre o valor obtido em 2. pelo valor obtido 1., arredondado a 4 casas decimais.


    Resolução

    library(extraDistr)
    
    # # Por exemplo:
    # mu <- 2000
    # sigma <- 500
    # n <- 1000
    # seed <- 1940
    
    # 1. Valor exato
    
    LA <- qgumbel(0.90, mu = mu, sigma = sigma)
    LE <- LA + 1000
    
    pLE <- 1 - pgumbel(LE, mu = mu, sigma = sigma)
    valor_exato <- pLE / 0.10
    
    # 2. Valor Simulado
    
    set.seed(seed)
    
    amostra <- rgumbel(n, mu = mu, sigma = sigma)
    subconjunto_alerta <- amostra[amostra > LA]
    valor_simulado <- mean(subconjunto_alerta > LE)
    
    solution <-  valor_simulado / valor_exato

    O valor do quociente pretendido é 1.2406.


  5. Considere o par aleatório \((X,Y)\) com suporte no triângulo \(\mathcal{S}=\{(x,y) \in \mathbb{R}^2:0<x<y<1\}\) e função densidade de probabilidade conjunta: \[f_{X,Y}(x,y)=\frac{2x}{1-x}\,, (x,y) \in \mathcal{S}.\] Pretende-se validar a qualidade do gerador em R de números pseudoaleatórios Uniforme\((0,1)\), runif(n, 0, 1), confirmando se o algoritmo de geração reproduz fielmente a estrutura de dependência teórica entre estas variáveis. Para tal, siga os passos abaixo:

    Passo 1

    Obtenha o valor exato da \(P(Y >1.7X).\)

    Passo 2

    Implemente um algoritmo em R, utilizando exclusivamente o gerador runif(n, 0, 1), que garanta a estrutura de dependência do par \((X,Y)\) para simular \(P(Y > 1.7X)\). Fixe a semente em \(2254\) e realize \(K=100\) iterações. Em cada iteração do algoritmo, gere uma amostra do par \((X,Y)\) de dimensão \(n=1200\) seguindo as etapas:

    1. Geração de \(x\): sabendo que a função de distribuição marginal de \(X\) é \(F_X(x)=x^2 \,,0<x<1,\) gere \(u\) da Uniforme\((0,1)\) e determine \(x\) tal que \(x=F_X^{-1}(u)\);

    2. Geração de \(y\): para cada valor de \(x\) obtido, gere \(y\) a partir da distribuição condicional \((Y∣X=x) \sim\) Uniforme\((x,1).\) Para tal, gere \(z\) da Uniforme\((0,1)\) e obtenha \(y=F_{Y|X=x}^{-1}(z)\);

    3. Proporção de sucessos: contabilize a proporção de pares que satisfazem a condição \((y>1.7x).\)

    O valor simulado da \(P(Y >1.7X)\) corresponde à média das proporções obtidas nas \(K\) iterações.

    Indique o valor do quociente entre o valor simulado no Passo 2 e o valor exato obtido no Passo 1:


    1. 1.0702
    2. 1.0098
    3. 1.1606
    4. 1.1846
    5. 1.0864

    Resolução

    Cálculo da expressão da probabilidade exata:

    \[ \begin{split} P(Y > mX) & = \int_{0}^{1} \int_{0}^{y/m} \frac{2x}{1-x} \, dx \, dy =\\ & = \int_{0}^{1} \left[ -2x - 2\ln(1-x) \right]_{0}^{y/m} \, dy =\\ & = \int_{0}^{1} \left( -\frac{2y}{m} - 2\ln\left(1 - \frac{y}{m}\right) \right) \, dy =\\ & = \int_{0}^{1} -\frac{2y}{m} \, dy - 2 \int_{0}^{1} \ln(m-y) \, dy + 2 \int_{0}^{1} \ln(m) \, dy =\\ & = \left[ -\frac{y^2}{m} \right]_{0}^{1} - 2 \left[ -(m-y)\ln(m-y) + (m-y) \right]_{0}^{1} + 2\ln(m) =\\ & = 2 - \frac{1}{m} + 2(m-1) \ln\left( \frac{m-1}{m} \right) \end{split} \]

    # # Por exemplo:
    # seed <- 1234
    # n <- 800
    # m <- 1.2
    
    # Passo 1. Valor exato
    p_exata <- 2 - 1 / m + 2 * (m - 1) * log((m - 1) / m)
    
    # Passo 2.  Valor simulado
    set.seed(seed)
    K <- 100
    
    proporcoes <- replicate(K, {
      u <- runif(n)
      x <- sqrt(u) 
      z <- runif(n)
      y <- x + z * (1 - x)
      mean(y > m * x)
    })
    
    p_simulada <- mean(proporcoes)
    
    solution <- p_simulada / p_exata

    O valor pedido é 1.0098.


  6. Considere que o número de chamadas recebidas por hora num “call-center” é modelado por uma variável aleatória \(X\) com distribuição de Poisson de média \(\lambda=3.8\). Seja \((X_1, X_2,\dots,X_n)\) uma amostra aleatória de dimensão \(n=50\) desta população. Pretende-se obter uma aproximação para a probabilidade de a média amostral, \(\bar{X}\), não se desviar da média populacional mais do que \(\epsilon=0.18\): \[P\left(|\bar{X}-3.8|\leq 0.18\right).\]

    Compare os seguintes dois métodos para obter essa aproximação:

    1. Interpretação frequencista

      Fixando a semente em \(1645\), gere \(M=178\) amostras independentes, cada uma com dimensão \(n=50\). Determine a frequência relativa das médias amostrais que satisfazem a condição indicada.

    2. Aplicação do TLC

      Recorra ao Teorema do Limite Central (TLC) para obter a aproximação da probabilidade indicada.

    Indique o valor do quociente entre os resultados obtidos pelo método a. e o método b.:


    1. 0.9286
    2. 0.9231
    3. 0.7092
    4. 0.9937
    5. 0.7710

    Resolução

    # Por exemplo:
    # epsilon <- 0.05
    # n <- 150
    # lambda <- 1
    # m <- 100
    # seed <- 1234
    
    set.seed(seed)
    
    count <- 0
    for(i in 1:m) {
        sample <- rpois(n, lambda)
        count <- (abs(mean(sample) - lambda) <= epsilon) + count
    }
    pemp <- count / m
    
    ptlc <- 2 * pnorm(epsilon / sqrt(lambda / n)) - 1
    
    solution <- pemp / ptlc
    opcoes
    ## [1] "0.9937" "0.7710" "0.9286" "0.9231" "0.7092"

    O valor do quociente é 0.9937.


  7. Considere que a variável aleatória \(X\), representando a duração de uma componente eletrónica (em anos), pode ser modelada pela função densidade de probabilidade \[ f_X(x) = \frac{\alpha}{\lambda} \left(\frac{x}{\lambda}\right)^{\alpha-1} \exp\left[-\left(\frac{x}{\lambda}\right)^{\alpha}\right],\quad x>0, \] onde \(\alpha>0\) e \(\lambda>0\) são os parâmetros desconhecidos do modelo.

    Suponha que \(\textbf{X} =(X_1,\ldots,X_n)\) é uma amostra aleatória de \(X\) e que a observação de \(n = 20\) durações da componente eletrónica resultou em

    \(\textbf{x} =\) (2.69, 2.16, 2.45, 1.43, 1.87, 1.49, 2.06, 2.8, 2.78, 1.7, 2.53, 2.71, 2.15, 1.86, 2.68, 2.17, 2.69, 1.74, 2.38, 1.83)

    Sejam \(\hat{\alpha}\) e \(\hat{\lambda}\) as estimativas de máxima verosimilhança de \(\alpha\) e \(\lambda\), respetivamente, obtidas a partir destes dados. No sistema de equações de verosimilhança, \(\hat{\lambda}\) pode ser escrita como função explícita de \(\hat{\alpha}\), enquanto \(\hat{\alpha}\) não pode ser isolada analiticamente, podendo ser obtida resolvendo uma equação implícita. Para obter \(\hat{\alpha}\) numericamente, recorra à função uniroot do R no intervalo \([4, 6]\), sem utilizar qualquer outro argumento opcional dessa função.

    Determine a estimativa de máxima verosimilhança da duração mediana da componente eletrónica. Indique o resultado arredondado a 2 casas decimais.


    Resolução

    Sistema de equações de verosimilhança: \[ \begin{array}{l} \frac{\partial \ln L(\alpha,\lambda|\textbf{x})}{\partial \lambda} = 0 \quad \Rightarrow \quad \hat{\lambda} = \left(\frac{1}{n} \sum_{i=1}^n {x_i}^{\hat{\alpha}}\right)^{1/\hat{\alpha}} \\ \frac{\partial \ln L(\alpha,\lambda|\textbf{x})}{\partial \alpha} = 0 \quad \Rightarrow \quad \frac{1}{\hat{\alpha}} + \frac{1}{n} \sum_{i=1}^n \ln x_i - \frac{\sum_{i=1}^n {x_i}^\hat{\alpha} \ln x_i}{\sum_{i=1}^n {x_i}^\hat{\alpha}} = 0 \end{array} \]

    Tempo de vida mediano (\(x_d\)): \[ F_X(x_d)=1-\exp\left[-(x_d/\lambda)^\alpha\right] = 0.5 \ \implies\ x_d=\lambda (\ln 2)^{1/\alpha} \]

    # # Por exemplo:
    # x <- c(6.42, 2.8, 5.44, 2.03, 1.42, 10.11, 4.59, 1.94, 4.44, 5.09, 1.21, 5.11,
    # 3.59, 4.29, 8.67, 1.87, 6.81, 10.24, 6.07, 1.24)
    # liminf <- 1.11
    # limsup <- 10.34
    
    funcao <- function(k){
      1 / k + mean(log(x)) - sum((x^k) * log(x)) / sum(x^k)
    }
    
    alpha_est <- uniroot(funcao, lower = liminf, upper = limsup)$root
    
    lambda_est <- mean(x^alpha_est)^(1 / alpha_est)
    
    solution <- lambda_est *(log(2)^(1 / alpha_est))

    A estimativa pedida é igual a 2.25 anos.


  8. Considere uma variável aleatória \(X\) a representar a vida útil de partículas radioativas (em dias), com distribuição Exponencial de valor esperado \(\lambda>0\). Com base numa amostra aleatória de \(X\) de dimensão \(n\), denotada por \((X_1,\ldots,X_n)\), pode-se determinar intervalos aleatórios de confiança para \(\lambda\) com um grau de confiança \(\gamma\), utilizando as seguintes variáveis fulcrais:

    1. Variável fulcral 1 (distribuição exata): \(\dfrac{2\,n\,\bar{X}}{\lambda} \sim \chi^2_{(2\,n)}\),

    2. Variável fulcral 2 (distribuição aproximada): \(\dfrac{\bar{X}-\lambda}{\sqrt{\lambda^2/n}} \overset{a}{\sim} N(0,1)\),

    onde \(\bar{X}=\frac{1}{n} \sum_{i=1}^n X_i\) (média amostral).

    Proceda à geração de \(m =\) 1200 amostras de dimensão \(n =\) 61 de \(X\) no software R, fixando a semente em 4943 e admitindo que \(\lambda=\) 5.

    Para cada amostra gerada, calcule o intervalo de confiança para \(\lambda\) com base quer na variável fulcral 1 quer na variável fulcral 2, com um grau de confiança de \(\gamma =\) 0.97, bem como a proporção de intervalos de confiança gerados que contém o valor do parâmetro \(\lambda=\) 5, denotada por \(p_1\) (variável fulcral 1) e \(p_2\) (variável fulcral 2).

    Indique o quociente entre as proporções \(p_1\) e \(p_2\):


    1. 1.0964
    2. 1.0772
    3. 1.0083
    4. 0.9991
    5. 0.9118

    Resolução

    Intervalo aleatório de confiança a \(100\,\gamma\,\%\) para \(\lambda\) (exato): \[ \left]\frac{2\,n\,\bar{X}}{F^{-1}_{\chi^2_{(2\,n)}} ((1+\gamma)/2)},\ \frac{2\,n\,\bar{X}}{F^{-1}_{\chi^2_{(2\,n)}} ((1-\gamma)/2)} \right[ \]

    Intervalo aleatório de confiança a \(100\,\gamma\,\%\) para \(\lambda\) (aproximado): \[ \left] \frac{\sqrt{n}\bar{X}}{\sqrt{n} + F^{-1}_{N(0,1)} ((1+\gamma)/2)}\,,\ \frac{\sqrt{n}\bar{X}}{\sqrt{n} - F^{-1}_{N(0,1)} ((1+\gamma)/2)} \right[ \]

    # # Por exemplo:
    # seed <- 3194
    # m <- 1200
    # n <- 21
    # lambda <- 12
    # gama <- 0.9
    
    set.seed(seed)
    
    q1a <- qchisq((1 - gama) / 2, 2 * n)
    q1b <- qchisq((1 + gama) / 2, 2 * n)
    q2 <- qnorm((1 + gama) / 2)
    
    p1 <- 0
    p2 <- 0
    for (i in 1:m) {
      amostra <- rexp(n, 1/lambda)
      media <- mean(amostra)
      liminf1 <- 2 * n * media / q1b
      limsup1 <- 2 * n * media / q1a
      p1 <- p1 + ((lambda - liminf1) * (lambda - limsup1) < 0)
      liminf2 <- (sqrt(n) / (sqrt(n) + q2)) * media
      limsup2 <- (sqrt(n) / (sqrt(n) - q2)) * media
      p2 <- p2 + ((lambda - liminf2) * (lambda - limsup2) < 0)
    }
    
    solution <- p1 / p2

    A proporção pedida é igual a 0.9991.


  9. Seja \((X_1,\dots,X_n)\) uma amostra aleatória de uma população \(X\) com distribuição Beta\((a,b)\). Com base no Teorema do Limite Central, para uma amostra de dimensão \(n\) suficientemente grande, define-se a variável fulcral: \[ Z = \frac{\overline{X} - \mu}{\sqrt{\sigma^2 /n}}= \frac{\sqrt{n}\left(\overline{X} - \frac{a}{a+b}\right)} {\sqrt{\frac{a b}{(a+b)^2(a+b+1)}}}\overset{a}{\sim} N(0,1) \] Fixando o parâmetro \(b=4\), pretende-se testar a hipótese \(H_0:μ=μ_0=\frac{1.5}{1.5+4 }\) contra \(H_1:μ=μ_1=\frac{1.82}{1.82+4}\). Para tal, recorre-se à estatística de teste \(Z_0\), obtida a partir de \(Z\) sob a validade de \(H_0\). O teste é realizado para uma amostra de dimensão \(n=54\) e ao nível de significância \(\alpha = 0.1\).

    1. Calcule o valor aproximado da probabilidade do erro de 2.ª espécie, \(\beta\).

    2. Obtenha uma estimativa da probabilidade do erro de 2.ª espécie, \(\hat{\beta}\). Para tal, fixe a semente em 981 e gere \(m=1750\) amostras de dimensão \(n=54\) da distribuição Beta\((1.82, 4)\) em R. Aplique o teste de hipóteses anterior a cada uma das amostras e calcule a frequência relativa de não rejeições de \(H_0\) nas \(m=1750\) simulações.

    Indique o valor do quociente \(\frac{\hat{\beta}}{\beta}\):


    1. 1.0564
    2. 1.0894
    3. 1.1186
    4. 1.0806
    5. 1.0614

    Resolução

    # # Por exemplo:
    # seed <- 200
    # n <- 45
    # b <-2
    # a0 <- 1.5
    # a1 <- 2.0
    # alpha <- 0.1
    # m <- 2000
    
    # Parâmetros sob H0 e H1
    mu0 <- a0 / (a0 + b)
    var0 <- (a0 * b) / ((a0 + b)^2 * (a0 + b + 1))
    mu1 <- a1 / (a1 + b)
    var1 <- (a1 * b) / ((a1 + b)^2 * (a1 + b + 1))
    
    # 1. VALOR TEÓRICO (BETA)
    z_crit <- qnorm(1 - alpha)
    # Ponto crítico:
    x_crit <- z_crit*sqrt(var0 / n) + mu0
    # Beta = P(X_bar <= x_crit | H1)
    beta_teorico <- pnorm((x_crit - mu1) / sqrt(var1 / n))
    
    # 2. SIMULAÇÃO
    set.seed(seed)
    m_erros <- replicate(m, {
      amostra <- rbeta(n, shape1 = a1, shape2 = b)
      # Estatística Z0 sob H0
      z0_sim <- (mean(amostra) - mu0) / sqrt(var0 / n)
      # Contar se NÃO rejeitou (Z0 <= z_crit)
      z0_sim <= z_crit
    })
    beta_simulado <- mean(m_erros)
    
    # RESULTADO
    solution <- beta_simulado / beta_teorico

    O valor pedido é igual a 1.0564.


  10. Versão corrigida

    O ficheiro Modular_P800.txt (14.5KB) contém um conjunto de dados obtido na avaliação de uma máquina de análises clínicas. Em 10 etapas, foram registados os erros na medição das concentrações de 200 soluções preparadas obtendo-se assim uma amostra com 2000 observações.

    Sendo \(X\) o erro de uma medição feita pela máquina, pretende-se testar a hipótese \(H_0: X \sim \text{Normal}(0, \sigma^2 = 1.44)\) contra \(H_1: X \not \sim \text{Normal}(0, \sigma^2 = 1.44)\) recorrendo ao teste de ajustamento do qui-quadrado. Este procedimento requer o agrupamento dos dados em intervalos ou classes e, para avaliar o efeito que o número de intervalos poderá ter no resultado do teste de hipóteses, proceda da seguinte forma:

    1. Fixe a semente em \(2947\) e obtenha uma subamostra formada por 80% das observações retiradas ao acaso do conjunto de dados inicial ordenado de forma crescente;

    2. Agrupe a subamostra em \(k\) intervalos equiprováveis sob \(H_0\), para dois valores de \(k\) dados pelas regras de Sturges e de Freedman-Diaconis (ambas disponíveis no R);

    3. Para cada um dos conjuntos de dados agrupados calcule o valor-p aproximado do teste referido;

    4. Indique o valor do quociente entre o maior e o menor dos valores-p obtidos, arrendado a 4 casas decimais.


    Resolução

    # # Por exemplo:
    # seed <- 5248
    # sigma <- 1.2
    # percent <- 70
    
    set.seed(seed)
    
    dados <- scan("Modular_P800.txt") |> sort()
    dados <- sample(dados, percent / 100 * length(dados))
    
    ks <- c(nclass.Sturges(dados), nclass.FD(dados))
    valores.p <- c()
    
    for(k in ks) {
      probs <- 0:k / k
      limites <- qnorm(probs, sd = sigma)
      freq_esp <- table(cut(dados, limites))
      valor.p <- chisq.test(freq_esp)[["p.value"]]
      valores.p <- c(valores.p, valor.p)
    }
    
    solution <- max(valores.p) / min(valores.p)

    O valor do quociente pedido é igual a 9.1196.