#===============================================================================
#. Name   : Funções de densidade de probabilidades
#. Author : José Cláudio Faria (Tinn-R Team)
# Date   : 23/05/2022 10:32:39
# Aim    : Entendimento da abordagem das funções de densidade de
#          probabilidades no R
# URL    : https://tinn-r.org/en
# Mail   : joseclaudio.faria@gmail.com
#===============================================================================
#. Distribuições abordadas:
# Funções de densidade de probabilidades:
# - Normal
# - t de Student
# - F de Snedecor
# Funções probabilidades:
# Tukey
#===============================================================================

#. Forma como o R trata o assunto:
#... d_dist*: retorna a densidade associada a um valor "x"
#... p_dist*: retorna a probabilidade associada a um quantil "q"
#... q_dist*: retorna o quantil associado a uma probabilidade "p"
#... r_dist*: retorna uma amostra da distribuição de tamanho "n"

# *: dist - uma das distribuições (Normal, t, F, chisq, etc.)

#. Observando os principais argumentos
# Observar que os argumentos das funções variam de acordo com a distribuição!

#.. Normal: norm
str(dnorm)
str(pnorm)
str(qnorm)
str(rnorm)

# mean: média
# sd: desvio padrão

#.. t de Student: t
str(dt)
str(pt)
str(qt)
str(rt)

# df: graus e liberdade

#.. F de Snedecor: f
str(df)
str(pf)
str(qf)
str(rf)

# df1: graus de liberdade do numerador
# df1: graus de liberdade do denominador

#. Plotando curvas de densidade de probabilidades
library(fdth)

#... Normal
plot(fdt(rnorm(1e6,   # Histograma (k=100) de uma amostra (n=1e4) na escala de densidades
               m=10,
               sd=2),
         k=100),
     type='d',
     ylim=c(0, .3))

curve(dnorm(x,        # Sobrepondo a função de densidade de probabilidades
            m=10,
            sd=2),    # Observar os mesmos parâmetros usados em rnorm: m e sd
      n=1e3,
      lwd=3,
      col=2,
      add=T)

curve(dnorm(x,
            m=8,
            sd=1.5),  # Observar parâmetros diferentes dos usados em rnorm: m e sd
      n=1e3,
      lwd=2,
      col=1,
      add=T)

curve(dnorm(x,
            m=12,
            sd=2),    # Observar parâmetros diferentes dos usados em rnorm: m e sd
      n=1e3,
      lwd=2,
      col=4,
      add=T)

#... t de Student
plot(fdt(rt(1e6,
            df=20),
         k=100),
     type='d',
     ylim=c(0, .5))

curve(dt(x,
         df=20),   # Observar o mesmo parâmetro usados em rt: df=20
      n=1e3,
      lwd=3,
      col=2,
      add=T)

curve(dt(x,
         df=2),    # Observar parâmetro diferente do usado em rt: df=2
      n=1e3,
      lwd=3,
      col=1,
      add=T)

curve(dt(x,
         df=50),   # Observar parâmetro diferente do usado em rt: df=50
      n=1e3,
      lwd=3,
      col=4,
      add=T)

#... F de Snedecor
plot(fdt(rf(1e6,
            df1=6,
            df2=20),
         k=100),
     type='d',
     ylim=c(0, 1))

curve(df(x,
         df1=6,
         df2=20),  # Observar os mesmos parâmetros usados em rf: df1 e df2
      n=1e3,
      lwd=3,
      col=2,
      add=T)

curve(df(x,
         df1=5,
         df2=12),  # Observar parâmetros diferentes dos usados em rf: df1 e df2
      n=1e3,
      lwd=2,
      col=1,
      add=T)

curve(df(x,
         df1=8,
         df2=52),  # Observar parâmetros diferentes dos usados em rf: df1 e df2
      n=1e3,
      lwd=2,
      col=4,
      add=T)

#. Calculado probabilidades associadas aos quantis
#.. Normal: X ~ N(Média = 10, Desvio Padrão = 2)
plot('',
     type='n',
     xlim=c(2, 18),
     ylim=c(0, .3),
     xlab='x',
     ylab='f(x)')

curve(dnorm(x,
            m=10,
            sd=2),
      n=1e3,
      lwd=2,
      col=1,
      add=T)

#... P(X < 10)
pnorm(q=10,
      m=10,
      s=2)

integrate(dnorm,
          low=10,
          upp=Inf,
          m=10,
          s=2)[1]

#... P(X < 8)
pnorm(q=8,
      m=10,
      s=2)

integrate(dnorm,
          low=-Inf,
          upp=8,
          m=10,
          s=2)[1]

#... P(X > 8)
pnorm(q=8,
      m=10,
      s=2,
      lower=F)

1- pnorm(q=8,
         m=10,
         s=2)

integrate(dnorm,
          low=8,
          upp=Inf,
          m=10,
          s=2)[1]

#... P(8 < X < 12)
pnorm(q=12,
      m=10,
      s=2) -
pnorm(q=8,
      m=10,
      s=2)

integrate(dnorm,
          low=8,
          upp=12,
          m=10,
          s=2)[1]

#.. t de Student: X ~ t (Graus liberdade = 20)
plot('',
     type='n',
     xlim=c(-5, 5),
     ylim=c(0, .5),
     xlab='t',
     ylab='f(t)')

curve(dt(x,
         df=20),
      n=1e3,
      lwd=2,
      col=1,
      add=T)

#... P(X < 0)
pt(q=0,
   df=20)

integrate(dt,
          low=-Inf,
          up=0,
          df=20)[1]

#... P(X > 0)
pt(q=0,
   df=20,
   lower=F)

integrate(dt,
          low=0,
          up=Inf,
          df=20)[1]

#... P(X < 2)
pt(q=2,
   df=20)

integrate(dt,
          low=-Inf,
          up=2,
          df=20)[1]

#... P(X > 2)
1 - pt(q=2,
       df=20)

pt(q=2,
   df=20,
   lower=F)

integrate(dt,
          low=2,
          up=Inf,
          df=20)[1]

#.. F de Snedecor: X ~ F (Graus liberdade 1 = 3, Graus liberdade 2 = 20)
plot('',
     type='n',
     xlim=c(0, 18),
     ylim=c(0, 1),
     xlab='x',
     ylab='f(x)')

curve(df(x,
         df1=3,
         df2=20),
      n=1e3,
      lwd=2,
      col=1,
      add=T)

#... P(X < 1)
pf(q=1,
   df1=3,
   df2=20)

integrate(df,
          df1=3,
          df2=20,
          low=0,
          upp=1)[1]

#... P(X > 1)
1 - pf(q=1,
       df1=3,
       df2=20)

pf(q=1,
   df1=3,
   df2=20,
   lower=F)

integrate(df,
          df1=3,
          df2=20,
          low=0,
          upp=1)[1]

#. Achando os quantis associados às probabilidades
#.. Normal: X ~ N (10, 2)
plot('',
     type='n',
     xlim=c(2, 18),
     ylim=c(0, .3),
     xlab='x',
     ylab='f(x)')

curve(dnorm(x,
            m=10,
            sd=2),
      n=1e3,
      lwd=2,
      col=1,
      add=T)

#... Quantil deixa 50% dos menores valores abaixo
qnorm(p=50/100,
      m=10,
      s=2)

#... Quantil deixa 30% dos maiores valores acima
qnorm(p=(100-30)/100,
      m=10,
      s=2)

qnorm(p=30/100,
      m=10,
      s=2,
      low=F)


#.. t: X ~ t (df=20)
plot('',
     type='n',
     xlim=c(-5, 5),
     ylim=c(0, .5),
     xlab='t',
     ylab='f(t)')

curve(dt(x,
         df=20),
      n=1e3,
      lwd=2,
      col=1,
      add=T)

#... Quantil deixa 50% dos menores valores abaixo
qt(p=50/100,
   df=20)

#... Quantil deixa 30% dos maiores valores acima
qt(p=(100-30)/100,
   df=20)

qt(p=30/100,
   df=20,
   low=F)

#.. F de Snedecor: X ~ F (df1=3, df1=20)
plot('',
     type='n',
     xlim=c(0, 18),
     ylim=c(0, 1),
     xlab='x',
     ylab='f(x)')

curve(df(x,
         df1=3,
         df2=20),
      n=1e3,
      lwd=2,
      col=1,
      add=T)

#... Qual quantil deixa 95% dos menores valores abaixo
qf(p=95/100,
   df1=3,
   df2=20)

qf(p=5/100,
   df1=3,
   df2=20,
   low=F)

#... Qual quantil deixa 5% dos maiores valores acima
qf(p=(100-5)/100,
   df1=3,
   df2=20)

qf(p=5/100,
   df1=3,
   df2=20,
   low=F)


#. Retirando uma amostra de dados pseudo-aleatórios das distribuições
#.. Normal: X ~ N (10, 2)
amo_n <- rnorm(n=1e6,
               m=10,
               sd=2)

plot(fdt(amo_n,
         k=150),
     type='d')

curve(dnorm(x,
            m=10,
            sd=2),   # Observar os mesmos parâmetros usados em rnorm: m e sd
      n=1e3,
      lwd=3,
      col=2,
      add=T)

#.. F de Snedecor: X ~ F (df1=3, df1=20)
amo_f <- rf(n=1e6,
            df1=3,
            df2=20)

plot(fdt(amo_f,
         k=150),
     type='d')

curve(df(x,
         df1=3,
         df2=20),   # Observar os mesmos parâmetros usados em rnorm: m e sd
      n=1e3,
      lwd=3,
      col=2,
      add=T)

#. The Studentized Range Distribution - Tukey
# É uma funções de probabilidade, não de densidade de probabilidades.
# Segue a mesma sintaxe mas não existem as funções:
# d_dist
# r_dist.

#... A forma da distribuição
curve(ptukey(x,
             nm=4,
             df=16),
      from=0,
      to=8,
      n=1e2,
      lwd=2,
      col=2,
      xlab='quantil (x)')

# Exemplo baseado no capítulo sobre Teste de comparação de médias - TCMM
# da apostla de CET076, pag. 81: q0.05(4, 16)
# (https://lec.pro.br/download/faria/apostilas/CET018_10ed_2pf.pdf)

# Obtendo o quantil associado a um probabilidade
qtukey(p=5/100,
       nmeans=4,
       df=16,
       lower.tail=FALSE)  # 4.046093

abline(v=4.046093,    # Adicionando linhas de referência
       h=1 - (5/100),
       lty=2)

# Obtendo a probabilidade associada a um quantil
ptukey(q=4.046093,
       nmeans=4,
       df=16,
       lower.tail=FALSE)

#. Duncan, Dunnet, etc
# As distribuições de Duncan, Dunnet e outras usadas em testes de comparação de médias,
# devem ser vistas nos pacotes que implementas esses testes.
# Muitas vezes são aproximações de alguma outra distribuição.