#===============================================================================
#. Título: Revisão de Variáveis Aleatórias Contínuas - VAC
#. Curso : Introdução à estatística
#. Autor : José Cláudio Faria/UESC/DCET
# Data  : 2023-05-13 11:37:01
# Objetivos:
#===============================================================================
# - Apresentar as funções do R para as principais distribuições
#   - Densidade de probabilidade                         (d)
#   - Probabilidade acumulada                            (p)
#   - Quantil correspondente a uma dada probabilidade    (q)
#   - Amostras aleatórias da distribuição                (r)
#===============================================================================
#.. Padrão do R
#... ddist - densidade
#... pdist - probabilidade
#... qdist - quantil
#... rdist - random
#===============================================================================

#.. Exemplos
#... Normal
str(dnorm)
str(pnorm)
str(qnorm)
str(rnorm)

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

#... Qui-quadrado
str(dchisq)
str(pchisq)
str(qchisq)
str(rchisq)

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

#. Distribuição Normal
#.. Função de densidade de probabilidades - FDP
f_norm <- function(x,
                   m,
                   s)
{
  1/(s*sqrt(2*pi))*exp(-1/2*((x-m)/s)^2)
}

#... Argumentos (f_norm):
# x: quantil
# m: média
# s: desvio padrão

nv <- seq(0,
          20,
          len=1e3)

plot(f_norm(x=nv,
            m=10,
            s=2) ~ nv,
     type='l',
     col='orange',
     lwd=5)

curve(dnorm(x,    # o x aqui é necessário (mesmo que o objeto não exista)
            m=10,
            s=2),
      add=TRUE,
      n=1e3)

# Integrando f_norm
integrate(f_norm,
          -Inf,
          Inf,
          m=10,
          s=2)$value  # Considerando m=10, s=2

integrate(f_norm,
          -Inf,
          10,
          m=10,
          s=2)$value

integrate(f_norm,
          10,
          Inf,
          m=10,
          s=2)$value

# Para o R f_norm e dnorm são uma mesma função
# é posssivel integrar (também) a função dnorm
integrate(dnorm,
          -Inf,
          Inf,
          m=10,
          s=2)$value  # Considerando m=10, s=2

integrate(dnorm,
          -Inf,
          10,
          m=10,
          s=2)$value

integrate(dnorm,
          10,
          Inf,
          m=10,
          s=2)$value

# Comparando a integral com uma função de alto nível que fornece probabilidade: pnorm
integrate(f_norm,
          -Inf,
          10,
          m=10,
          s=2)$value

pnorm(10,
      m=10,
      s=2)

# Exemplos
# X ~ N(m=100, s=10), calcule as probabilidades:
# a) P[X < 95]
pnorm(95,
      100,
      10)

# via integral
integrate(dnorm,
          -Inf,
          95,
          m=100,
          sd=10)$value

# b) P[90 < X < 110]
pnorm(110,
      100,
      10) -
pnorm(90,
      100,
      10)

# via integral
integrate(dnorm,
          90,
          110,
          m=100,
          sd=10)$value

# c) P[X > 95]
1 - pnorm(95,
          100,
          10)

## TRUE:  Calcula as probabilidades da forma P(X <= x).
## FALSE: Calcula as probabilidades da forma P(X > x).
pnorm(95,
      100,
      10,
      lower=F)

# via integral
integrate(dnorm,
          95,
          Inf,
          m=100,
          sd=10)$value

# Calcular o quantil que delimita 40% dos menores valores
# O valor de n tal que P(X <= n) = 0.4
qnorm(40/100,
      m=10,
      s=2)

# Utilizando a rnorm para gerar 100 números
# de uma normal com media=10 e sd=2
rnorm(100, # aceita também 1e2, 1E2
      m=10,
      s=2)

# Simulando dados, verificando a densidade e sobrepondo a curva
# da função densidade de probabilidades
hist(rnorm(1e4,
           m=10,
           sd=2),
     prob=TRUE,
     nclass=100)

curve(f_norm(x,
             m=10,
             s=2),
      n=1e3,
      col='orange',
      lw=5,
      add=T)

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


# Várias sobrepostas
plot(NA,
     xlab='x',
     ylab='f(x)',
     xlim=c(-20, 40),
     ylim=c(0, .4),
     type='n')

# Número de curvas
nc <- seq(1,
          10,
          by=.5)

# Define as cores
cores <- rainbow(length(nc))

for (i in nc)
  curve(dnorm(x,
              m=10,
              sd=i),
        col=cores[i],
        add=TRUE,
        n=1e3)

#.  Distribuição t
#.. Função de densidade de probabilidades - FDP
f_student <- function(x,
                      df)
{
  (gamma((df+1)/2))/(sqrt(df*pi)*gamma(df/2))*((1+x^2/df))^-((df+1)/2)
}

#... Argumentos (f_student):
#  x: quantil
# df: graus de liberdade (n-1)

tv <- seq(-10,
           10,
           len=1e3)

plot(f_student(x=tv,
               df=10) ~ tv,
     type='l',
     col='orange',
     lwd=5)

curve(dt(x,    # o x aqui é necessário (mesmo que o objeto não exista)
         df=10),
      add=TRUE,
      n=1e3)

# Integrando
integrate(f_student,
          -Inf,
          Inf,
          df=10)$value  # Considerando df=10

integrate(f_student,
          -Inf,
          0,
          df=10)$value

integrate(f_student,
          0,
          Inf,
          df=10)$value

# Para o R f_student e dt são uma mesma função
# é posssivel integrar (também) a função dt
integrate(dt,
          -Inf,
          Inf,
          df=10)$value  # Considerando df=10

integrate(dt,
          -Inf,
          0,
          df=10)$value

integrate(dt,
          0,
          Inf,
          df=10)$value

# Comparando a integral com uma função de alto nível: pt
round(integrate(f_student,
                -Inf,
                0,
                df=10)$value,
      5) ==
round(pt(0,
         df=10),
      5)

# Calcular o quantil referente a probabilidade resultante da pt
(x = round(pt(0,
              df=10),
           5))

qt(x,
   df=10)

# Calcular o quantil que delimita 15% dos menores valores com df=10
qt(15/100,
   df=10)

# Utilizando a rt para gerar 100 números com df=10
rt(100,
   df=10)

# Simulando dados, verificando a densidade e sobrepondo a curva
# da função densidade de probabilidades
hist(rt(1e4,
        df=10),
     prob=TRUE,
     nclass=100)

curve(f_student(x,
                df=10),
      n=1e3,
      col='orange',
      lw=5,
      add=T)

curve(dt(x,
         df=10),
      n=1e3,
      add=T)

# Várias sobrepostas
plot(NA,
     xlab='t',
     ylab='f(t)',
     xlim=c(-10, 10),
     ylim=c(0, .4),
     type='n')

# Número de curvas
nc <- seq(1,
          20,
          by=1)

# Define as cores
cores <- rainbow(length(nc))

for (i in nc)
  curve(dt(x,
           df=i),
        col=cores[i],
        add=TRUE,
        n=1e3)

# Sobrepondo uma normal padrão
curve(dnorm(x),
      lw=3,
      col=gray(.2),
      add=TRUE,
      n=1e3)

#. Distribuição Qui-quadrado
#.. Função de densidade de probabilidades - FDP
f_chisq <- function(x,
                    df)
{
  (1/(2^(df/2)*gamma(df/2)))*x^(df/2-1)*exp(-x/2)
}

#... Argumentos (f_chisq):
#  x: quantil
# df: graus de liberdade (n-1)

qv <- seq(0,
          50,
          len=1e3)

plot(f_chisq(x=qv,
             df=10) ~ qv,
     type='l',
     col='orange',
     lwd=5)

curve(dchisq(x,
             df=10),
      add=TRUE,
      n=1e3)

# Integrando
integrate(f_chisq,
          0,
          Inf,
          df=10)$value  # Considerando df=10

(p1 <- integrate(f_chisq,
                 0,
                 20,
                 df=10)$value)

(p2 <- integrate(f_chisq,
                 20,
                 Inf,
                 df=10)$value)

p1 + p2
round(p1 + p2,
      5) == 1

# Para o R f_chisq e dchisq são uma mesma função!
# Assim, é posssivel integrar (também) a função dchisq
integrate(dchisq,
          0,
          Inf,
          df=10)$value  # Considerando df=10

integrate(dchisq,
          0,
          20,
          df=10)$value

integrate(dchisq,
          20,
          Inf,
          df=10)$value

# Comparando a integral com uma função de alto nível: pchisq
round(integrate(f_chisq,
                0,
                20,
                df=10)$value,
      3)

round(pchisq(20,
             df=10),
      3)

# Calcular o quantil que delimita 95% dos menores valores
qchisq(95/100,
       df=10)

# O quantil 20 (df=10) delimita qual percentual dos maiores valores?
1 - pchisq(20,
           df=10)

pchisq(20,
       df=10,
       lower=FALSE)

integrate(dchisq,
          20,
          Inf,
          df=10)$value

# Utilizando a rchisq para gerar 100 números com df=10
rchisq(100,
       df=10)

# Simulando dados, verificando a densidade e sobrepondo a curva
# da função densidade de probabilidades
hist(rchisq(1e4,
            df=10),
     prob=TRUE,
     nclass=100)

curve(f_chisq(x,
              df=10),
      n=1e3,
      col='orange',
      lwd=5,
      add=T)
curve(dchisq(x,
             df=10),
      n=1e3,
      add=T)

# Várias sobrepostas
plot(NA,
     xlab='X2',
     ylab='f(X2)',
     xlim=c(0, 80),
     ylim=c(0, .5),
     type='n')

# Número de curvas
nc <- 40

# Define as cores
cores <- rainbow(nc)

for (i in 1:nc)
  curve(dchisq(x,
               df=i),
        col=cores[i],
        add=TRUE,
        n=1e3)

#. Distribuição F
#.. Função de densidade de probabilidades - FDP
f_snedecor <- function(x,
                       df1,
                       df2)
{
  ((gamma((df1+df2)/2))/(gamma(df1/2)*gamma(df2/2)))*((df1/df2)^(df1/2))*((x^((df1-2)/2))/((1+(df1/df2)*x)^((df1+df2)/2)))
}

#... Argumentos (f_chisq):
#   x: quantil
# df1: graus de liberdade do numerador (n-1)
# df2: graus de liberdade do numerador (n-1)

fv <- seq(0,
           10,
           len=1e3)

plot(f_snedecor(x=fv,
                df1=5,
                df2=15) ~ fv,
     type='l',
     col='orange',
     lwd=5)

curve(df(x,
         df1=5,
         df2=15),
      add=TRUE,
      n=1e3)

# Integrando
integrate(f_snedecor,
          0,
          Inf,
          df1=5,
          df2=15)$value  # Considerando df1=5, df2=15

(p1 <- integrate(f_snedecor,
                 0,
                 3,
                 df1=5,
                 df2=15)$value)

(p2 <- integrate(f_snedecor,
                 3,
                 Inf,
                 df1=5,
                 df2=15)$value)

p1 + p2

# Para o R f_snedecor e df são uma mesma função
# é posssivel integrar (também) a função df
integrate(df,
          0,
          Inf,
          df1=5,
          df2=15)$value  # Considerando df1=5, df2=15

integrate(df,
          0,
          3,
          df1=5,
          df2=15)$value

integrate(df,
          3,
          Inf,
          df1=5,
          df2=15)$value

# Comparando a integral com uma função de alto nível: pf
round(integrate(f_snedecor,
                0,
                3,
                df1=5,
                df2=15)$value,
      5) ==
round(pf(3,
         df1=5,
         df2=15),
      5)

# Calcular o quantil referente a probabilidade resultante da pf
(x = round(pf(3,
              df1=5,
              df2=15),
           5))

qf(x,
   df1=5,
   df2=15)

# Calcular o quantil que delimita 60% dos menores valores
# O valor de F tal que P(X <= F) = 0.6
qf(60/100,
   df1=5,
   df2=15)

# Utilizando a rf para gerar 100 números com o df1=5 e df2=15
rf(100,
   df1=5,
   df2=15)

# Simulando dados, verificando a densidade e sobrepondo a curva
# da função densidade de probabilidades
hist(rf(1e4,
        df1=5,
        df2=15),
     prob=TRUE,
     nclass=100)

curve(f_snedecor(x,
                 df1=5,
                 df2=15),
      n=1e3,
      col='orange',
      lw=5,
      add=T)

curve(df(x,
         df1=5,
         df2=15),
      n=1e3,
      add=T)

# Várias sobrepostas
plot(NA,
     xlab='F',
     ylab='f(F)',
     xlim=c(0, 10),
     ylim=c(0, 1.2),
     type='n')

# Número de curvas
nc <- 30

# Define as cores
cores <- rainbow(nc)

#.. Variando apenas df1
for (i in 1:nc)
  curve(df(x,
           df1=i,    # <- df1
           df2=15),
        col=cores[i],
        add=TRUE,
        n=1e3)

##... Variando df1 e df2
#plot(NA,
#     xlab='F',
#     ylab='f(F)',
#     xlim=c(0, 10),
#     ylim=c(0, 1.2),
#     type='n')
#
#for (i in 1:10)
#  for (j in 1:nc)
#    curve(df(x,
#             df1=i,    # <- df1
#             df2=j),   # <- df2
#          col=cores[j],
#          add=TRUE,
#          n=1e3)