options(dig=12) require(statmod) ls(2) #?matvec ?gauss.quad ## Exemplo 1 f1 <- function(x) x^2 ## \int_{-1}^1 f1(x) dx ## a) fazer via r-yacas/sympy! # R: 2/3 I <- 2/3 ## b) Ia <- integrate(f1, -1, 1) c(I, Ia$val, (Ia$val-I)/I) ## c) N <- 100 nos <- ((1-(-1))*(1:N-0.5))/N -1 nos Ia <- (1-(-1)) * sum(f1(nos))/N Ia c(I, Ia, (Ia-I)/I) ## Nota: para ajustar os nós a limites de integração {a,b} ## (interpolacao linear/regra de tres) a foirma geral é # nos <- ((b-a)*(1:N-0.5))/N + a ## d) xw <- gauss.quad(5, kind="legendre") #xw <- gauss.quad(5, kind="jacobi") #xw <- gauss.quad(200, kind="chebyshev1") ## requer muitos pontos!! #f1 <- function(x) x^2 * sqrt(1-x^2) #xw <- gauss.quad(5, kind="chebyshev2") ## requer ajustes names(xw) xw Ia <- sum(f1(xw$n) * xw$w) Ia c(I, Ia, (Ia-I)/I) ## Exemplo 1a ## \int_{0}^{3} f1(x) dx ## a) fazer via r-yacas/sympy! # R: 9 I <- 9 ## b) Ia <- integrate(f1, 0, 3) Ia c(I, Ia$val, (Ia$val-I)/I) ## c) N <- 100; a <- 0 ; b <- 3 #nos <- ((3-0)*(1:N-0.5))/N -0 nos <- ((b-a)*(1:N-0.5))/N + a nos Ia <- (3-0) * sum(f1(nos))/N Ia c(I, Ia, (Ia-I)/I) ## d) xw <- gauss.quad(5, kind="legendre") #xw <- gauss.quad(5, kind="jacobi") #xw <- gauss.quad(5, kind="chebychev1") Ia <- ((3-0)/2) * sum(f1(((3-0)/2)*xw$n + (0+3)/2) * xw$w) Ia c(I, Ia, (Ia-I)/I) ## ## \frac{b-a}{2} \int_{-1}^{1} f1(((b-a)/2)x + (a+b)/2) dx f1a <- function(x, a, b){ x <- (((b-a)/2)*x) + ((a+b)/2) x^2 } Ia <- ((3-0)/2) * sum(f1a(xw$n, a=0, b=3) * xw$w) Ia ## Exemplo 2 f2 <- function(x) exp(x) ## \int_{0}^1 f2(x) dx ## a) fazer via r-yacas/sympy! # I <- exp(1) - 1 I ## b) Ia <- integrate(f2, 0, 1) Ia c(I, Ia$val, (Ia$val-I)/I) ## c) N <- 50; a <- 0 ; b <- 1 #nos <- ((1-0)*(1:N-0.5))/N -0 nos <- ((b-a)*(1:N-0.5))/N + a nos Ia <- (1-(0)) * sum(f2(nos))/N Ia c(I, Ia, (Ia-I)/I) ## d) f2a <- function(x, a, b){ x <- (((b-a)/2)*x) + ((a+b)/2) exp(x) } xw <- gauss.quad(5, kind="legendre") Ia <- ((1-0)/2) *sum(f2a(xw$n, a=1, b=0) * xw$w) Ia c(I, Ia, (Ia-I)/I) ## \int_{1}^4 f2(x) dx ## a) fazer via r-yacas/sympy! # 51.87986820468519 I <- exp(4) - exp(1) I ## b) Ia <- integrate(f2, 1, 4) Ia c(I, Ia$val, (Ia$val-I)/I) ## c) N <- 50; a <- 1 ; b <- 4 nos <- ((b-a)*(1:N-0.5))/N + a #nos <- ((4-1)*(1:N-0.5))/N +1 nos Ia <- (4-(1)) * sum(f2(nos))/N Ia c(I, Ia, (Ia-I)/I) ## d) xw <- gauss.quad(5, kind="legendre") Ia <- ((4-1)/2) *sum(f2a(xw$n, a=4, b=1) * xw$w) Ia c(I, Ia, (Ia-I)/I) #### ## Exemplo 3 f3 <- function(x) exp(-x/5) ## \int_{0}^{\Infty} f3(x) dx ## a) fazer via r-yacas/sympy! I <- 5 ## b) Ia <- integrate(f3, 0, Inf) Ia c(I, Ia$val, (Ia$val-I)/I) ## c) ## d) #f3a <- function(x) exp(-x/5) * exp(x) [x^alpha * exp(-x)] f3a <- function(x){ exp(4*x/5)} xw <- gauss.quad(10, kind="laguerre") xw Ia <- sum(f3a(xw$n) * xw$w) Ia c(I, Ia, (Ia-I)/I) ## \int_{1}^{3} f3(x) dx ## a) fazer via r-yacas/sympy! # R (maxima): 1.349595584919777 I <- 1.349595584919777 ## b) Ia <- integrate(f3, 1, 3) Ia c(I, Ia$val, (Ia$val-I)/I) ## c) ## d) f3b <- function(x, a, b){ x <- (((b-a)/2)*x) + ((a+b)/2) exp(-x/5) } xw <- gauss.quad(10, kind="legendre") xw Ia <- ((3-1)/2)*sum(f3b(xw$n, a=1, b=3) * xw$w) Ia c(I, Ia, (Ia-I)/I) #f3c <- function(x) exp(-x/5) * sqrt(1-x^2) [1/sqrt(1-x^2)] f3c <- function(x, a, b){ y <- (((b-a)/2)*x) + ((a+b)/2) exp(-y/5) * sqrt(1-x^2) } xw <- gauss.quad(100, kind="chebyshev1") xw Ia <- ((3-1)/2)*sum(f3c(xw$n, a=1, b=3) * xw$w) Ia c(I, Ia, (Ia-I)/I) ## note a necessidade de muitos pontos para uma aproximacao razoável! ## Exemplo 4 # f(x) = e^{-5 x} f4 <- function(x) exp(-5*x) ## \int_{0}^{\Infty} f3(x) dx ## a) fazer via r-yacas/sympy! # R: 1/5 I <- 0.2 Ia <- integrate(f4, 0, Inf) Ia c(I, Ia$val, (Ia$val-I)/I) #f4a <- function(x) exp(-5*x) * x^(-alpha) exp(x) [x^alpha * exp(-x)] #f4a <- function(x) exp(-5*x) * exp(-alpha *log(x)) * exp(x) [x^alpha * exp(-x)] f4a <- function(x, alpha=1) exp(-4*x - alpha*log(x)) f4a <- function(x, alpha) exp(-4*x)/(x^alpha) alfa <- 0 xw <- gauss.quad(100, kind="laguerre", alpha=alfa) xw Ia <- sum(f4a(xw$n, alpha=alfa) * xw$w) Ia c(I, Ia, (Ia-I)/I) ## Exemplo 5 ## OBS: f(x) é o núcleo da chi^2 #f(x) = x^{\frac{k}{2}-1} exp{-x/2} f5 <- function(x) x^((k/2)-1) * exp(-x/2) #I = 2^{k/2} \Gamma(k/2) k <- 6 2^(k/2) * gamma(k/2) I = 16 Ia <- integrate(f5, 0, Inf) Ia c(I, Ia$val, (Ia$val-I)/I) ## OBS: Note uma forma mais adequada de escrever f_5(x)!!!!! curve(f5(x), 0, 10) f5 <- function(x) exp(((k/2)-1)*log(x) - (x/2)) curve(f5(x), 0, 10, add=T, col=2) ## d) #f5 <- function(x) x^((k/2)-1) * exp(-x/2) * * x^(-alpha) exp(x) [x^alpha * exp(-x)] f5 <- function(x, alpha) x^((k/2)-1-alpha) * exp(x/2) alfa <- 0 xw <- gauss.quad(5, kind="laguerre", alpha=alfa) xw Ia <- sum(f5(xw$n, alpha=alfa) * xw$w) Ia c(I, Ia, (Ia-I)/I) ## cauchy #f(x) = \frac{1}{pi} \[\frac{\gamma}{(x-x_0)^2 + \gamma^2} \] x0 <- -1 g <-1.5 f5 <- function(x) 1/((x-x0)^2 + g^2) curve(f5, from=-5, to=4) I <- pi/g I Ia <- integrate(f5, -Inf, Inf) Ia c(I, Ia$val, (Ia$val-I)/I) ## adicionando termo para Hermite #f5a <- function(x) 1/((x-x0)^2 + g^2) * exp(x^2) [exp(-x^2)] f5a <- function(x) exp(x^2)/((x-x0)^2 + g^2) xw <- gauss.quad(10, kind="hermite") xw Ia <- sum(f5a(xw$n) * xw$w) Ia c(I, Ia, (Ia-I)/I) ## a aproximação é ruim... e nam melhora muito adicionando mais pontos... #### f6 <- function(x) exp(-x^2/2) I <- sqrt(2*pi) I Ia <- integrate(f6, -Inf, Inf) Ia c(I, Ia$val, (Ia$val-I)/I) #f6a <- function(x) exp(-x^2/2)*exp(x^2) [exp(-x^2)] f6a <- function(x) exp(x^2/2) xw <- gauss.quad(10, kind="hermite") xw Ia <- sum(f6a(xw$n) *xw$w) Ia c(I, Ia, (Ia-I)/I) ## Note que o seguinte fornece o mesmo resultado !!!! f6b <- function(x) sqrt(2) Ia <- sum(f6b(xw$n) *xw$w) Ia c(I, Ia, (Ia-I)/I) ## de forma geral, para ## funcoes do tipo f(x) = e^{-k x^2} ka <- 3 f7 <- function(x, k) exp(-k * x^2) integrate(f7, -Inf, Inf, k=ka) xw <- gauss.quad(10, kind="hermite") f7a <- function(x, k) exp((1-k) * x^2) sum(f7a(xw$n, k=ka) *xw$w) ## mas, por substituição de variável, podemos fazer f7b <- function(x, k) sqrt(1/k) sum(f7b(xw$n, k=ka) * xw$w) ########################################### ########################################### ## Integração em 2D ## \int f(x,y) w(x,y) dx dy \approx \sum_i \sum_j w_i w_j f(x_i, y_j) ## estratégia simplificada ## supondo mesmos nós e pesos para ambas dimensões ## construindo grid de nós e pesos Q2D <- function(x,w) { if(length(x) != length(w)) stop("x e w devem ter mesmo comprimento") X <- as.matrix(expand.grid(x, x)) W <- as.matrix(expand.grid(w, w)) dimnames(X) <- dimnames(W) <- NULL w <- exp(rowSums(log(W))) list(nodes=X,weigths=w) } f11 <- function(x) x[1]^2 + x[2]^2 ## Exemplo: \int f11 em [-1, 1] x [-1, 1] ## solucao analitica #I <- ## pelo método do ponto médio N <- 100; a <- -1 ; b <- 1 nos <- ((b-a)*(1:N-0.5))/N + a nos pesos <- rep((b-a)/N, N) xw2d <- Q2D(nos, pesos) Ia <- sum(xw2d$w * f11(xw2d$n)) Ia ## pelo método da quadratura N <- 30 xw <- gauss.quad(N) xw2d <- Q2D(xw$n, xw$w) Ia <- sum(xw2d$w * f11(xw2d$n)) Ia ## fazer aqui tb pela adapt() do R (pacote adapt) e/ou outras opcoes! ## ## Exemplo ??? (acrescentar outro exemplo aqui com limites ## de integracao que nao sejam entre -1 e 1 ## ## Exemplo ??? (acrescentar outro exemplo com outra funcao aqui com limites ## de integracao que nao sejam infinitos e/ou diferentes entre as dimensoes ## por exemplo de (0, \Infty) x (-Infty, Infty) ## a funcao Q2D pode ser facilmente generalizada para ## diferentes pontos de integracao em diferentes dimensões Q2Da <- function(x1, w1, x2, w2) { if(length(x1) != length(w1)) stop("x1 e w1 devem ter mesmo comprimento") if(length(x2) != length(w2)) stop("x2 e w2 devem ter mesmo comprimento") X <- as.matrix(expand.grid(x1, x2)) W <- as.matrix(expand.grid(w1, w2)) dimnames(X) <- dimnames(W) <- NULL w <- exp(rowSums(log(W))) list(nodes=X,weigths=w) } ## conferir, verificando que usando esta temos os mesmos resultados dos exemplos anteriores! ######################################################## ## n dimensões ## Exemplo: dimensão genérica ## considera a seguinta função (quen é beeeemmmm suave!) ## f(x) = exp(\sum_{i=1}^d x_i) f16 <- function(x) {exp(rowSums(x))} # resultado analítico #I <- (exp(1)-exp(-1))^d ## grid com mesmos nós em todas dimensoes!!! QnD <- function(x, w, d){ nx <- length(x) if(nx != length(w)) stop("x e w devem ter mesmo comprimento") X <- expand.grid(as.data.frame(matrix(x, nx, d))) W <- expand.grid(as.data.frame(matrix(w, nx, d))) # dimnames(X) <- dimnames(W) <- NULL w <- exp(rowSums(log(W))) list(nodes=X,weigths=w) } ## aproximação por ponto médio in 2D N <- 10; d <- 2 ; a <- -1; b <- 1 I <- (exp(1)-exp(-1))^d I nos <- ((b-a)*(1:N-0.5))/N + a nos pesos <- rep((b-a)/N, N) xw <- QnD(nos, pesos, d=d) Ia <- sum(xw$w*f16(xw$n)) c(Ia,I,(I-Ia)/I) ## note o problema com o aumento de dimensão, # por exemplo fazer com d = 5 ## agora por quadratura gaussiana N <- 4; d <- 2 xw <- gauss.quad(N) xw <- QnD(xw$n, xw$w, d=d) Ia <- sum(xw$w*f16(xw$n)) c(Ia,I,(I-Ia)/I) ## fazer agora: ## achar um N para cada um dos dois metodos que forneca um erro comparável ## comparar os tempos computacionais para cada um deles com d=5 ## Desafio # Considere o modelo de poisson com 2 covariáveis e efeitos aleatórios # cf o esquema de simulacao abaixo. set.seed(0) n <- 50 x <- runif(n) z1 <- sample(c(0,1),n,replace=TRUE) z2 <- sample(c(0,1),n,replace=TRUE) b <- rnorm(4) X <- model.matrix(~x) Z <- cbind(z1,1-z1,z2,1-z2) eta <- X%*%c(0,3) + Z%*%b mu <- exp(eta) y <- rpois(mu,mu) # obter: ## um função que avalie log{f(y,b)}... ## Obs: depois de montar a função considere obte-la agora de uma ## forma vetorizada e menos sensivel a "overflow" ## Diferentes abordagems ## 1. quadratura por pontos médios ## 2. quadratura gaussiana ## 3. aproximação de Laplace ## Dica: escrever uma função de gradiente e Hessiano de log{f(y,b)} em relacao a b ## encontrar modas a posteriori ## e avaliar a paroximnação de Laplace nestas modas ### 4. Monte Carlo ("força bruta") (com Uniforme) ### 5. Laplace om amostragem por importância