##===============================================================================
## Título: Variáveis aleatórias contínuas - VAC
## Curso : Introdução à estatística
## Autor : José Cláudio Faria/UESC/DCET
## Data  : 21/12/2017 17:38:06
## Versão: v5
## Objetivos:
##===============================================================================
## a) Compreensão das funções de densidade de probabilidades
## b) Cálculo de probabilidades via integrais dessas funções
##===============================================================================

#----------------------------------------------------------------------#
#                         Exemplos - FDP                               #
#----------------------------------------------------------------------#
#--------------------
# Exemplo 1
#--------------------
# Definindo uma FDP: muito simples
fy <- function(y)
  ifelse(y >= 0 & y <= 1, 
         2*y, 
         0)

# Gráfico da função
plot(fy,
     xlab='y',
     ylab='f(y)')

# Verificando se fy é uma FDP
# Primeira condição:
fy(seq(-2, 2, by=.1))
fy(seq(-2, 2, by=.1)) >= 0
all(fy(seq(-2, 2, by=.1)) >= 0)

# Segunda condição
Iteg_y <- integrate(fy, 
                    0, 
                    1)
str(Iteg_y)

integrate(fy, 
          0, 
          1)$value

# Visualizando algumas probabilidades
# Entre 0 e .5
polygon(x=c(seq(0, 
                .5, 
                l=10), 
            .5),
        y=c(fy(seq(0, 
                   .5, 
                   l=10)), 
            0),
        col=gray(.8))

# Entre .5 e 1
polygon(x=c(.5,
            seq(.5, 
                1, 
                l=10), 
            1),
        y=c(0,
            fy(seq(.5, 
                   1, 
                   l=10)), 
            0),
        col=gray(.5))

# Função text com as coordenadas de posição do texto
text(c(.3, .7), 
     c(.2, .2), 
     c(expression(P('0 < Y < .5'), 
                  P('.5 < Y < 1'))),
     cex=1.2,
     col='red')

# Calculando densidades de probabilidade
# Entre 0 e .5
integrate(fy, 
          0, 
          .5)$value

# Entre .5 e 1
integrate(fy, 
          .5, 
          1)$value

#--------------------
# Exemplo 2
#--------------------
# Definindo outra FDP: menos simples
# Disponível apenas no seminário de referência!
fx <- function(x)
  ifelse(x < 0, 
         0, 
         2*exp(-2*x))

# Gráfico da função
plot(fx,
     xlab='x',
     ylab='f(x)')

# Escolhendo limite
plot(fx, 0, 5)

# Questão A
# Verificando se fy é uma FDP
# Primeira condição
fx(seq(-2, 10))
fx(seq(-2, 10, by=.1)) >= 0
all(fx(seq(-2, 10, by=.1)) >= 0)

# Segunda condição
integrate(fx, 
          0, 
          Inf)$value

# Visualizando algumas probabilidades
polygon(x=c(1, 
            seq(1, 
                5, 
                l=20)),
        y=c(0, 
            fx(seq(1, 
                   5, 
                   l=20))),
        col=gray(.5))

polygon(x=c(0.2, 
            seq(0.2, 
                0.8, 
                l=20), 
            0.8),
        y=c(0, 
            fx(seq(0.2, 
                   0.8, 
                   l=20)), 
            0),
        col=gray(.8))

# Função text com as coordenadas de posição do texto
text(c(1.3, .5), 
     c(.08, .3), 
     c(expression(Pb,
                  Pc)),
     col='red')

# Obtendo as probabilidades pedidas
# Questão B
text(1.6, 
     .2, 
     round(integrate(fx, 
                     1, 
                     Inf)$value, 2),
     col='red')

# Questão C
text(.9, 
     .7, 
     round(integrate(fx, 
                     .2, 
                     .8)$value, 2),
     col='red')

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

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

plot(fnorm(x=tes, 
           m=10,
           s=2) ~ tes,
     type='l')

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

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

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

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

# Para o R fnorm 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: pnorm
round(integrate(fnorm,
                -Inf,
                10,
                m=10,
                s=2)$value, 
      5) ==
round(pnorm(10,
            m=10,
            s=2), 
      5)

# 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
# df: graus de liberdade (n-1)
ft <- function(x, df)
  (gamma((df+1)/2))/(sqrt(df*pi)*gamma(df/2))*((1+x^2/df))^-((df+1)/2)

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

plot(ft(x=tes, 
        df=10) ~ tes,
     type='l')

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

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

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

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

# Para o R ft 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(ft,
                -Inf,
                0,
                df=10)$value, 
      5) ==
round(pt(0,
         df=10), 
      5)

# 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)

#----------------------------------------------------------------------#
#                         Distribuição Qui-quadrado                    #
#----------------------------------------------------------------------#
# Função de densidade de probabilidades
# df: graus de liberdade (n-1)
fchisq <- function(x, df)
  (1/(2^(df/2)*gamma(df/2)))*x^(df/2-1)*exp(-x/2)

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

plot(fchisq(x=tes, 
            df=10) ~ tes,
     type='l')

curve(dchisq(x,
             df=10),
      col='red',
      lw=2,
      add=TRUE,
      n=1e3)

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

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

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

# Para o R fchisq e dchisq são uma mesma função
# é 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(fchisq,
                0,
                20,
                df=10)$value, 
      5) ==
round(pchisq(20,
             df=10), 
      5)

# 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
# df1: graus de liberdade do numerador (n-1)
# df2: graus de liberdade do numerador (n-1)
ff <- 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)))

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

plot(ff(x=tes, 
        df1=5,
        df2=15) ~ tes,
     type='l')

curve(df(x,
         df1=5,
         df2=15),
      col='red',
      lw=2,
      add=TRUE,
      n=1e3)

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

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

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

# Para o R ff 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(ff,
                0,
                3,
                df1=5,
                df2=15)$value, 
      5) ==
round(pf(3,
         df1=5,
         df2=15), 
      5)

# 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)              

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

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,    # Variando df1 e df2
             df2=j),
          col=cores[j],
          add=TRUE,
          n=1e3)