## HPP com N \sim P(lambda=50) ## simulação no intervalo [0, 100] x <- numeric() while(sum(x) <= 100) x <- c(x, rexp(1, rate=1/2)) x <- x[-length(x)] n <- length(x) pontos <- cumsum(x) plot(pontos, rep(1.2, n), type="n", ylim=c(0, 1.2)) arrows(pontos, rep(1, n), pontos, rep(0, n), length=0.1) plot(sort(pontos));segments(0,0,n,100) plot(density(pontos)) ## a density "engana" (não possui correcao de borda...) rug(pontos) ## função G(d) = 1 - exp( - lambda * pi * d^2) #distEV <- sort(as.vector(dist(pontos))) #plot(distEV, (1:length(distEV))/(length(distEV)+1), type="s") #G <- function(d, lambda=100) 1 - exp(-lambda * pi * d^2) #curve(G(x, lambda=0.001), 0, 100, add=T) ## outras opcoes descritivas em resumos de PP (funções K, J, G, etc) ## ## simulacao de IPP ## via amostragem por rejeição ## Para simular IPP com intensidade λ(x): ## 1. encontrar λ "grande" tal que λ(x) = λ p(x) ## (se necessário maximizar λ(x) ) λ = sup λ(x) ## 2. Simular um HPP X* com taxa λ sobre a área ## 3. Obter IPP X aceitando cada ponto de X* com prob λ(t)/λ ## ## alternativa: ## 1. simular no de pontos de N ~ P( \int \lambda(x) dx ) ## 2. discretizar o dominio em conjunto finito de pontos x_i ## \lambda(x_i)/\int \lambda(x) dx = \lambda(x_i)/\sum \lambda(x_i) lambda.x <- function(x, par, log=FALSE){ eta <- par[1] + par[2] * x if(log) return(eta) else return(exp(eta)) } b0 <- 1.2 ; b1 <- 0.05 (LX <- round(optimize(lambda.x, c(0, 100), par=c(b0, b1), maximum=TRUE)$objective)) X <- runif(LX, 0, 100) p <- lambda.x(X, c(b0, b1))/LX pontos <- sort(X[runif(LX) <= p]) (n <- length(pontos)) plot(pontos, rep(1.2, n), type="n", ylim=c(0, 1.2)) arrows(pontos, rep(1, n), pontos, rep(0, n), length=0.1) plot(sort(pontos)); segments(0,0,n,100) plot(density(pontos)) ## a density "engana" (não possui correcao de borda...) rug(pontos) ## outras opcoes descritivas em resumos de PP (funções K, J, G, etc) nlikIPP <- function(par, X, interval, numerico=FALSE){ a <- interval[1] ; b = interval[2] if(numerico) intLx <- integrate(lambda.x, low=a, upp=b, par=par)$value else intLx <- exp(par[1])/par[2] * (exp(b*par[2]) - exp(a*par[2])) slxi <- sum(lambda.x(x=X, par = par, log=TRUE)) return(-(slxi - intLx)) } nlikIPP(c(1, 0.1), X=pontos, interval=c(0, 100)) nlikIPP(c(1, 0.1), X=pontos, interval=c(0, 100), num=T) optim(c(1.2, 0.05), nlikIPP, X = pontos, interval=c(0, 100))[1:2] optim(c(1.2, 0.05), nlikIPP, X = pontos, interval=c(0, 100), num=T)[1:2] ## o seguinte tb não resolve lambdaR.x <- function(x, par, log=FALSE){ eta <- par[1] * par[2]^x if(log) return(log(eta)) else return(eta) } lambda.x(50, c(b0, b1), log=T) lambdaR.x(50, exp(c(b0, b1)), log=T) lambda.x(50, c(b0, b1), log=F) lambdaR.x(50, exp(c(b0, b1)), log=F) nlikIPPR <- function(par, X, interval, numerico=FALSE){ a <- interval[1] ; b = interval[2] intLx <- integrate(lambdaR.x, low=a, upp=b, par=par)$value # if(numerico) intLx <- integrate(lambdaR.x, low=a, upp=b, par=par)$value # else intLx <- exp(par[1])/par[2] * (exp(b*par[2]) - exp(a*par[2])) slxi <- sum(lambdaR.x(x=X, par = par, log=TRUE)) return(-(slxi - intLx)) } nlikIPPR(exp(c(1, 0.1)), X=pontos, interval=c(0, 100)) #nlikIPP(c(1, 0.1), X=pontos, interval=c(0, 100), num=T) ap <- optim(exp(c(1.2, 0.05)), nlikIPPR, X = pontos, interval=c(0, 100))[1:2] log(ap[[1]]) ## ## Reescalonando para uma janela de observação unitária nlikIPP <- function(par, X, interval){ X <- (X - interval[1])/interval[2] intLx <- integrate(lambda.x, low=0, upp=1, par=par)$value slxi <- sum(lambda.x(x=X, par = par, log=T)) return(-(slxi - intLx)) } nlikIPP(c(1, 0.1), X=pontos, interval=c(0, 100)) optim(c(1, 0.1), nlikIPP, X = pontos, interval=c(0, 100))[1:2] ### loglambda<-function(x, alpha, beta) { l<-alpha+beta*x return(l) } L <- function(alphabeta, x){ l<-loglambda(x, alpha=alphabeta[1], beta=alphabeta[-1]) l<-sum(l) intL<-integrate(f=function(x, alpha, beta){ exp(loglambda(x, alpha, beta))}, low=0, up=100, alpha=alphabeta[1], beta=alphabeta[-1]) l<-l-intL$value return(l)#Optim minimises } optim(c(1, 0.03), L, x = pontos, control=list(fnscale=-1,maxit=10000, reltol=.00001))[1:2] require(adapt) loglambda<-function(x, alpha, beta){alpha+beta*x} L<-function(alphabeta, x){ l<-sum(loglambda(x, alpha=alphabeta[1], beta=alphabeta[-1])) intL<-adapt(1, c(0,100), functn=function(x, alpha, beta){ exp(loglambda(x, alpha, beta))}, alpha=alphabeta[1], beta=alphabeta[-1]) l<-l-intL$value return(l)#Optim minimises } optim(c(1, 0.03), L, x = pontos, control=list(fnscale=-1,maxit=10000, reltol=.00001))[1:2] ## ## IPP em 2-D ## ### loglambda<-function(x, alpha, beta) { l<-alpha+sum(beta*c(x, x*x, prod(x))) return(l) } L<-function(alphabeta, x) { l<-apply(x,1,loglambda, alpha=alphabeta[1], beta=alphabeta[-1]) l<-sum(l) intL<-adapt(2, c(0,0), c(1,1), functn=function(x, alpha, beta) { exp(loglambda(x, alpha, beta)) }, alpha=alphabeta[1], beta=alphabeta[-1]) l<-l-intL$value # print(l) return(l)#Optim minimises } IPPfit <- function(formula, ini, intervalo, data){ ## só permite uma variavel na formula mf <- model.frame(formula, data) pts <- model.response(mf) X <- model.matrix(mf) ll.ipp <- function(par, ...){ lam.x <- function(x) exp(X%*%par) integrate ## precisa montar exp(X%*%beta) no grid sobre a area para integrate integrate(, area...)