# tests/testthat/test-localmax.R In Infusion: Inference Using Simulation

```cat(crayon::yellow("test likelihood with narrow peak among local maxima:\n"))

# Example from Wilkinson's ABC tutorial
# Where Infusion.options(mixmodGaussianModel) matters...
nips <- function(theta) {
s <- rnorm(n=1,mean=2*theta*(theta+2)*(theta-2),sd=sqrt(0.1+theta^2))
return(c(D=s))
}

lik <- function(theta) dnorm(x=2,mean=2*theta*(theta+2)*(theta-2),sd=sqrt(0.1+theta^2),log=FALSE)
liks <- sapply(seqx <- seq(-4,4,0.001),lik)

set.seed(123)
Dobs <- c(D=2) ## stands for the actual data to be analyzed

npoints <- 300
parsp <- data.frame(theta=runif(npoints,min=-4,max=4))
#
# library(abcrf)
# rfabc <- regAbcrf(theta~D, simuls)
# densityPlot(rfabc, data.frame(D=2), simuls)
# predict(rfabc,obs=data.frame(D=2), training=simuls)
#
densv <- infer_SLik_joint(simuls,stat.obs=Dobs)
currMSL <- MSL(densv) ## find the maximum of the log-likelihood surface
plot(currMSL); lines(seqx,liks/max(liks),col="red") # add true likelihood
if (inherits(currMSL\$jointdens, "dMixmod")) plot(currMSL\$jointdens, data=currMSL\$logLs) ## is revealing. Several more iterations are needed:
currMSL <- refine(currMSL,maxit=5)
plot(currMSL); lines(seqx,liks/max(liks),col="red") # better
if (inherits(currMSL\$jointdens, "dMixmod")) plot(currMSL\$jointdens, data=currMSL\$logLs) ## better

if (FALSE) { # when projections are bad
bigparsp <- data.frame(theta=runif(3000,min=-4,max=4))
prth <- project("theta", stats="D",data=bigsimuls)
prDobs <- project(Dobs,projectors=list(PRTH=prth))
prsimuls <- project(bigsimuls,projectors=list(PRTH=prth))
densv <- infer_SLik_joint(prsimuls,stat.obs=prDobs)
slik_j <- MSL(densv) ## find the maximum of the log-likelihood surface
plot(slik_j)
}
```

## Try the Infusion package in your browser

Any scripts or data that you put into this service are public.

Infusion documentation built on Sept. 29, 2022, 1:05 a.m.