# rm(list=ls())
options(OutDec = ",") 
#==============================================================================
# Exemplo do calculo do valor esperado de uma variavel aleatoria qui-quadrada 
# com 'nu' graus de liberdade pela regra do Simpson
#==============================================================================
# Funcao para a integracao numerica pela regra de Simpson
#==============================================================================
simpson <- function(a,b,f,nsub){
	# (a,b) e' o intervalo de integracao
	# f e' a funcao que se deseja integrar
	# nsub e' o numero de subintervalos
	xseq  <- seq(a,b,length=(nsub+1))
	inc   <- xseq[2]-xseq[1] 
	yseq  <- f(xseq)
	aux   <- inc/6
	resp1 <- aux*(yseq[1]+yseq[nsub+1]+2*sum(yseq[2:nsub]))
	resp2 <- 4*aux*sum(f(xseq[1:nsub]+0.5*inc)) # 4f((a+b)/2)
	return(resp1+resp2)
}

#==============================================================================
# Funcao objetivo - qui-quadrado 2 graus de liberdade
#==============================================================================
fx <- function(x){
	x*dgamma(x,2/2,1/2)
}

#==============================================================================
# Avaliacoes da integral
#==============================================================================
print("===============================")
simpson(0.0001,30,fx,4)
print("===============================")
simpson(0.0001,30,fx,20)
print("===============================")
simpson(0.0001,30,fx,100)
print("===============================")
simpson(0.0001,30,fx,1000)

[1] "==============================="

1,77448318994852

[1] "==============================="

1,99934909842815

[1] "==============================="

1,99998915519745

[1] "==============================="

1,99999020852037

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