Carregando os dados

DAP = c(10.2,28,29.5,15.8,30.3,28.3,18.5,26.8,34.2,22.3,28,38.5,23.5,17.2,17.8,38.5,18.9,35.5,17.8,28.9,34.2,18.7,27.9,31.8,16.7,22.5,32.5,29.1,32.9,41.8)

A) Calcule: média, mediana, moda, 1o quartil, 3o quartil, valor máximo, valor mínimo, amplitude total, variância, desvio padrão, erro padrão da média e coeficente de variação, utilizando duas maneiras diferentes.

Primeira maneira usando as funções do R

require(plotrix)

## Carregando pacotes exigidos: plotrix

require(goeveg)

## Carregando pacotes exigidos: goeveg

## This is GoeVeg 0.7.5 - build: 2024-05-17

paste("Média: ",mean(DAP) )

## [1] "Média:  26.5533333333333"

paste("Mediana: ", median(DAP))

## [1] "Mediana:  28"

paste("Moda: ", moda(DAP))

## [1] "Moda:  17.8" "Moda:  28"   "Moda:  34.2" "Moda:  38.5"

c(paste("1° Quantil: ",quantile(DAP, prob = 0)),
  paste("3° Quantil: ",quantile(DAP, prob = 0.75)))

## [1] "1° Quantil:  10.2"   "3° Quantil:  32.325"

paste("Valor Máximo: " ,max(DAP))

## [1] "Valor Máximo:  41.8"

paste("Valor mínimo: ",min(DAP))

## [1] "Valor mínimo:  10.2"

paste("Amplitude total: ", diff(range(DAP)))

## [1] "Amplitude total:  31.6"

paste("Variância :" ,var(DAP))

## [1] "Variância : 62.4729195402299"

paste("Desvio Padrão: " , sd(DAP))

## [1] "Desvio Padrão:  7.90398124619675"

paste("Erro padrão da média:" ,std.error(DAP))

## [1] "Erro padrão da média: 1.44306294087992"

paste("Coeficiente de variação: ", goeveg::cv(DAP, na.rm = FALSE))

## [1] "Coeficiente de variação:  0.297664370306179"

Pela segunda maneira construindo nossas formulas

paste("Média", sum(DAP)/length(DAP))

## [1] "Média 26.5533333333333"

# Criando a função para mediana 

mediana_GP = function(x){
  n = length(x)
  xord = sort(x)
  ifelse(n%%2 == 1,xord[(n+1)/2], mean(xord[n/2+0:1]) )
}
paste("Mediana: ", mediana_GP(DAP))

## [1] "Mediana:  28"

# Criando a função para moda que corrige o problema anterior 

moda_new <- function(x) {
  tab <- table(x)
  modas <- as.numeric(names(tab)[tab == max(tab)])
  return(modas)
}

paste("Valores modais:", paste(moda_new(DAP), collapse = ", "))

## [1] "Valores modais: 17.8, 28, 34.2, 38.5"

# Calculando a 1° médiana 
quant_1 = function(x){
  ordenados = sort(x)
  posicao = ceiling(length(x)*0.25)
  return(ordenados[posicao])
}
paste("1° Quantil: ", quant_1(DAP))

## [1] "1° Quantil:  18.7"

# Calculando a 3° Mediana
quant_3 = function(x){
  ordenados = sort(x)
  posicao = ceiling(length(x)*0.75)
  return(ordenados[posicao])
}
paste("3° Quantil: ", quant_3(DAP))

## [1] "3° Quantil:  32.5"

# Valores Minimos e máximos

posicao = function(x){
  xord = sort(x)
  primeirovalor = xord[1]
  ultimovalor = xord[length(xord)]
  return(c(primeirovalor,ultimovalor))
}

paste("Valores mínimo e máximo: ", posicao(DAP))

## [1] "Valores mínimo e máximo:  10.2" "Valores mínimo e máximo:  41.8"

# Amplitude

Amplitude = function(x){
    xord = sort(x)
    ultimovalor = xord[length(xord)]
    primeirovalor = xord[1]
    return(ultimovalor-primeirovalor)
    }
paste("Amplitude total:", paste(Amplitude(DAP), collapse = ","))

## [1] "Amplitude total: 31.6"

# Variancia

paste("Variância: " , sum((DAP-mean(DAP))^2/(length(DAP)-1)))

## [1] "Variância:  62.4729195402299"

# Desvio Padrão

paste("Desvio Padrão:", sqrt(sum((DAP-mean(DAP))^2/(length(DAP)-1))))

## [1] "Desvio Padrão: 7.90398124619675"

# Erro padrão da Média 

EPM = function(x){
  desvio = sd(x)
  n = length(x)
  return(desvio/sqrt(n))
  }

paste("Erro padrão da Média: ", EPM(DAP))

## [1] "Erro padrão da Média:  1.44306294087992"

# Coefiiciente de variação

CV = function(x){
  desvio = sd(x)
  xbar = mean(x)
  return((desvio/xbar)*100)
}


paste("Coeficiente de variação (%):",CV(DAP))

## [1] "Coeficiente de variação (%): 29.7664370306179"

B) Obtenha o histograma e o box-plot.

require(ggplot2)

## Carregando pacotes exigidos: ggplot2

DAP_data = as.data.frame(DAP)
ggplot(DAP_data, aes(x = DAP)) +
  geom_histogram(binwidth = 4, fill = "green", color = "black") +
  labs(title = "Histograma do DAP",subtitle = "Diâmetro das árvores para altura do peito", x = "Diâmetro para altura do peito (cm)", y = "Frequência (n)") +
  theme_light()

ggplot(DAP_data, aes(x = DAP, y = NULL)) + 
  geom_boxplot(fill = "green", color = "black") +
  labs(title = "Histograma do DAP",subtitle = "Diâmetro das árvores para altura do peito", x = "Diâmetro para altura do peito (cm)") +
  theme_light()

icm= function(x, confidence = 0.95){
  xbar = mean(x)
  n=length(x)
  sd = sd(x)
  qt = c(qt((1-confidence)/2,n-1),qt(1-(1-confidence)/2,n-1))
  return(xbar+qt*(sd/sqrt(n)))
  }

paste("Intervalo de confiança para média: ", paste(icm(DAP), collapse = ", "))

## [1] "Intervalo de confiança para média:  23.6019382311825, 29.5047284354841"

icv <- function(x, confidence = 0.95){
  n <- length(x)
  s2 <- var(x)
  alpha <- 1 - confidence
  chi <- qchisq(c(1 - alpha/2, alpha/2), df = n - 1)
  return((n - 1) * s2 / chi)
}
paste("Intervalo de confiança para Variância: ", paste(icv(DAP), collapse = ", "))

## [1] "Intervalo de confiança para Variância:  39.6243239987195, 112.900016966334"

icdp <- function(x, confidence = 0.95){
  n <- length(x)
  s2 <- var(x)
  alpha <- 1 - confidence
  chi <- qchisq(c(1 - alpha/2, alpha/2), df = n - 1)
  return(sqrt((n - 1) * s2 / chi))
}
paste("Intervalo de confiança para o desvio padrão: ", paste(icdp(DAP), collapse = ", "))

## [1] "Intervalo de confiança para o desvio padrão:  6.29478546089694, 10.6254419656941"

iccv <- function(x, confidence = 0.95){
  n <- length(x)
  xbar <- mean(x)
  s <- sd(x)
  cv <- s / xbar
  alpha <- 1 - confidence
  z <- qnorm(1 - alpha/2)
  se_cv <- cv * sqrt((1 + (cv^2)/2) / n)
  return(cv + c(-1, 1) * z * se_cv)
}

paste("Intervalo de confiança para o coeficiente de variação:", paste(iccv(DAP), collapse = ",") )

## [1] "Intervalo de confiança para o coeficiente de variação: 0.188814634232699,0.406514106379659"

Para o aproximado

icp_AP = function(p,n, confidence = 0.95){
  prop = p/n
  z = c(qnorm((1-confidence)/2),qnorm(1-(1-confidence)/2))
  return(prop+z*(sqrt(prop*(1-prop)/n)))
}

paste("Intervalo de confiança para proporção aproximado: ", paste(icp_AP(10,50,0.95),collapse = " , "))

## [1] "Intervalo de confiança para proporção aproximado:  0.0891276940520258 , 0.310872305947974"

Para o exato

icp_EX <- function(p, n, confidence = 0.95) {
  y = p         # número de sucessos
  alpha = 1 - confidence
  
  # Cálculo dos valores críticos da F
  fcalc1 = qf(alpha/2, df1 = 2 * (n - y + 1), df2 = 2 * y, lower.tail = F)
  fcalc2 = qf(1 - alpha/2, df1 = 2 * (y + 1), df2 = 2 * (n - y))
  
  # Limites inferior e superior do IC exato 
  lower = 1 / (1 + ((n - y + 1) * fcalc1)/y)
  
  upper = (y + 1) * fcalc2 / ((n - y) + (y + 1) * fcalc2)
  
  return(c(lower, upper))
}
paste("Intervalo de confiança para proporção exato: ", paste(icp_EX(10,50,0.95),collapse = " , "))

## [1] "Intervalo de confiança para proporção exato:  0.100302237472571 , 0.337183108383488"

Para o Aproximado sem correção e o exato iremos utilizar as formulas anteriores.

icp_AP(10,50,0.95)

## [1] 0.08912769 0.31087231

icp_EX(10,50,0.95)

## [1] 0.1003022 0.3371831

Para a Aproximação normal com correção de continuidade:

icp_AP_corr = function(p,n,confidence = 0.95){
  y = p
  alpha = 1- confidence
  z_crit = qnorm(alpha/2, lower.tail = F)
  
  # Criando os limites 
  LIC = ((y-0.5)/n) - (z_crit/sqrt(n)) * sqrt(((y-0.5)/n) * (1 - ((y - 0.5)/n)))
  LAC = ((y+0.5)/n) + (z_crit/sqrt(n)) * sqrt(((y+0.5)/n) * (1 - ((y + 0.5)/n)))
  return(c(LIC,LAC))
}

icp_AP_corr(10,50,0.95)

## [1] 0.08126174 0.32289801

Para a Aproximação devido a pratt:

icp_AP_pratt = function(p,n,confidence = 0.95){
  y = p
  n = n
  alpha = 1- confidence
  z_crit = qnorm(alpha/2, lower.tail = F)
  
  # Calculando os limites 
  
  LIC = (1 + (y / (n - y + 1))^2 * ((  (81 * y) * (n - y + 1) - (9*n) - 8 + (3*z_crit) *sqrt((9*y) * (n - y + 1) * ((9*n)+5-(z_crit)^2) + n + 1)  ) / ((81 * (y)^2) - (9* y) *(2 + (z_crit)^2) + 1))^3)^-1
  
  parte_1 = ((y+1)/(n-y))^2
  parte_2 = (81*(y+1)*(n-y)) - (9*n) - 8 - (3*z_crit)*sqrt(9*(y+1)*(n-y)*((9*n) + 5 - (z_crit)^2) + n + 1 )
  parte_3 = (81*(y+1)^2) - 9*(y+1)*(2 + (z_crit)^2) + 1
  
  LAC = (1 + (((y+1)/(n-y))^2) * ((81*(y+1)*(n-y) - 9*n - 8 - 3*z_crit * sqrt(9*(y+1)*(n - y)*(9*n+5-(z_crit)^2)+n+1))/(81*((y + 1)^2) - 9 *(y+1)*(2+(z_crit)^2)+1))^3)^(-1)
    
    
  return(c(LIC,LAC))
  
}
paste("Intervalo de confiança para proporção por Pratt", paste(icp_AP_pratt(10,50,0.95), collapse = ","))

## [1] "Intervalo de confiança para proporção por Pratt 0.100173830366947,0.337188123149737"