# rm(list=ls())
options(OutDec = ",") 
#==============================================================================
# Exemplo para gerar de uma mistura discreta de duas normais 
#==============================================================================
set.seed(54345)  # Semente
n      <- 100000 # Tamanho da amostra
omega  <- 0.3
mu1    <- 10
sigma1 <- 1
mu2    <- 15
sigma2 <- 1.5
y      <- rbinom(n,1,omega)        # gerando da mistura
mu     <- y*mu1    + (1-y)*mu2     # gerando da mistura
sigma  <- y*sigma1 + (1-y)*sigma2  # gerando da mistura
x      <- rnorm(n,mu,sigma)        # gerando da mistura

#==============================================================================
# Histograma
#==============================================================================
par(mfrow=c(1,1),lwd=2,cex.lab=1.5,cex.axis=1.5,lab=c(14,6,0),
    mar=c(4.5,5,2,1),bty="n")
hist(x,nclass=50,prob=T,main="",ylim=c(0,0.20),xlim=c(5,23),col="darkgreen",
     ylab=expression(f(x)),xlab=expression(x))
xseq <- seq(5,23,length=1000)
yseq <- ( omega*dnorm(xseq,mu1,sigma1)
        + (1-omega)*dnorm(xseq,mu2,sigma2) )
lines(xseq,yseq,col="red",lwd=3)

#==============================================================================
# Exemplo 1: gerar de uma mistura discreta de duas normais pelo metodo da
# reamostragem ponderada 
#==============================================================================
# Proposta U(0,20)
#==============================================================================
m      <- 500000  # Tamanho da amostra (reamostra de tamanho n)
xprop  <- runif(m,0,20)
pesos  <- 20*(     omega*dnorm(xprop,mu1,sigma1)
              +(1-omega)*dnorm(xprop,mu2,sigma2)) 
prob   <- pesos/sum(pesos)	
print(sum(prob))
y      <- sample(xprop,n,replace=T,prob=prob)

[1] 1

#==============================================================================
# Histograma
#==============================================================================
par(mfrow=c(1,1),lwd=2,cex.lab=1.5,cex.axis=1.5,lab=c(14,6,0),
    mar=c(4.5,5,2,1),bty="n")
hist(y,nclass=50,prob=T,main="",ylim=c(0,0.20),xlim=c(5,23),col="darkgreen",
     ylab=expression(f(x)),xlab=expression(x))
lines(xseq,yseq,col="red",lwd=3)

#==============================================================================
# Problema de cobertura na cauda superior!
#==============================================================================
# Probabilidade de X>20
#==============================================================================
print(1-(omega*pnorm(20,mu1,sigma1)+(1-omega)*pnorm(20,mu2,sigma2)))

#==============================================================================
# Probabilidade de X<0
#==============================================================================
print(omega*pnorm(0,mu1,sigma1)+(1-omega)*pnorm(0,mu2,sigma2))

[1] 0,0003003422
[1] 7,619853e-24

#==============================================================================
# Exemplo 2: gerar de uma mistura discreta de duas normais pelo metodo da 
# reamostragem ponderada 
#==============================================================================
# Proposta N(13,5; 3^2)
#==============================================================================
mu     <-  omega*mu1+(1-omega)*mu2
print("===============================")
print(mu)

#==============================================================================
# Valor esperado da variancia condicional
#==============================================================================
aux1   <- omega*(sigma1^2)+(1-omega)*(sigma2^2)

#==============================================================================
# Variancia do valor esperado condicional
#==============================================================================
aux2   <- omega*(mu1-mu)^2+(1-omega)*(mu2-mu)^2
sig2   <- aux1+aux2          
print("===============================")
print(sig2)   # Escolhemos um valor um pouco maior: 9

#==============================================================================
# Reamostragem ponderada
#==============================================================================
sigma  <- 3  
m      <- 500000  # Tamanho da amostra (reamostra de tamanho n)
xprop  <- rnorm(m,mu,sigma)
pesos  <- ( (     omega*dnorm(xprop,mu1,sigma1)
             +(1-omega)*dnorm(xprop,mu2,sigma2))/
		 dnorm(xprop,mu,sigma) )
prob   <- pesos/sum(pesos)	
print("===============================")
print(sum(prob))
w      <- sample(xprop,n,replace=T,prob=prob)

[1] "==============================="
[1] 13,5
[1] "==============================="
[1] 7,125
[1] "==============================="
[1] 1

#==============================================================================
# Histograma
#==============================================================================
par(mfrow=c(1,1),lwd=2,cex.lab=1.5,cex.axis=1.5,lab=c(14,6,0),
    mar=c(4.5,5,2,1),bty="n")
hist(w,nclass=50,prob=T,main="",ylim=c(0,0.20),xlim=c(5,23),col="darkgreen",
     ylab=expression(f(x)),xlab=expression(x))
lines(xseq,yseq,col="red",lwd=3)

#==============================================================================
# Verificacao rapida: ha' lgum valor maior que 20?
#==============================================================================
any(w>20)
1-pnorm(20,mu,sigma)

# rm(list=ls()) 
#==============================================================================
# graphics.off()
#==============================================================================
# Fim
#==============================================================================