example_raw_proj | R Documentation |
Example of the workflow with add_simulation
), implementing the method described in the original publication
(Rousset et al. 2017 <doi:10.1111/1755-0998.12627>), modified to use projectors.
if (Infusion.getOption("example_maxtime")>170) { ## Normal(mu,sd) model, with inefficient raw summary statistics: ## To illustrate that case we transform normal random deviates rnorm(,mu,sd) ## so that the mean of transformed sample is not sufficient for mu, ## and the variance of transformed sample is not sufficient for sd. blurred <- function(mu,s2,sample.size) { s <- rnorm(n=sample.size,mean=mu,sd=sqrt(s2)) s <- exp(s/4) return(c(mean=mean(s),var=var(s))) } set.seed(123) dSobs <- blurred(mu=4,s2=1,sample.size=20) ## stands for the actual data to be analyzed ## Sampling design as in canonical example parsp <- init_grid(lower=c(mu=2.8,s2=0.4,sample.size=20), upper=c(mu=5.2,s2=2.4,sample.size=20)) # simulate distributions dsimuls <- add_simulation(,Simulate="blurred", par.grid=parsp) ## Use projection to construct better summary statistics for each each parameter mufit <- project("mu",stats=c("mean","var"),data=dsimuls) s2fit <- project("s2",stats=c("mean","var"),data=dsimuls) ## additional plots for some projection method if (inherits(mufit,"HLfit")) mapMM(mufit,map.asp=1, plot.title=title(main="prediction of normal mean",xlab="exp mean",ylab="exp var")) if (inherits(s2fit,"HLfit")) mapMM(s2fit,map.asp=1, plot.title=title(main="prediction of normal var",xlab="exp mean",ylab="exp var")) ## apply projections on simulated statistics corrSobs <- project(dSobs,projectors=list("MEAN"=mufit,"VAR"=s2fit)) corrSimuls <- project(dsimuls,projectors=list("MEAN"=mufit,"VAR"=s2fit)) ## Analyze 'projected' data as any data (cf canonical example) densb <- infer_logLs(corrSimuls,stat.obs=corrSobs) } else data(densb) ######### if (Infusion.getOption("example_maxtime")>10) { slik <- infer_surface(densb) ## infer a log-likelihood surface slik <- MSL(slik) ## find the maximum of the log-likelihood surface } if (Infusion.getOption("example_maxtime")>500) { slik <- refine(slik,10, update_projectors=TRUE) ## refine iteratively }
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.