############################### #Example ABC: # # Inference for alpha-stable ############################### library(stabledist) ##load library ############ SEMI-AUTOMATIC ABC REJECTION SAMPLER - 3unknown parameters. # # Rather than pre-specifying threshold for distance: store # (delta,dist) pairs; and post-process these for different thresholds # ###upload data -- y source("Alphadata.R") gamma=1 ###define covariates for linear model covariate=function(y){ qu=quantile(y,1:10/11) return(c(qu,qu^2,qu^3)) ##10-quantiles are there 2nd/3rd powers } m=length(covariate(y)) ##number of covariates N=50000 ##number of samples par.st=matrix(0,nr=N,nc=3) par.st[,1]=runif(N,0.5,0.9) par.st[,2]=runif(N,-0.5,0.5) par.st[,3]=rnorm(N,0,2) data.st=matrix(0,nrow=N,ncol=m) n=length(y) for(i in 1:N){ ysim=rstable(n,par.st[i,1],par.st[i,2],1,par.st[i,3]) ##simulate data-set data.st[i,]=covariate(ysim) ##store covariatese } ##FIt linear models out1=lm(par.st[,1]~data.st) out2=lm(par.st[,2]~data.st) out3=lm(par.st[,3]~data.st) ##Create matrix of co-efficients -- we can ignore the constants in the lm B=matrix(0,nr=3,nc=m) B[1,]=out1\$coe[-1] B[2,]=out2\$coe[-1] B[3,]=out3\$coe[-1] ##PLOTS TO LOOK AT INFORMATION IN SUMMARIES par(mfrow=c(1,3)) plot(par.st[,1],out1\$fitted) plot(par.st[,2],out2\$fitted) plot(par.st[,3],out3\$fitted) S=function(y){ x=covariate(y) return(B%*%x) } Sobs=S(y) ##observed summary statistics d=length(Sobs) ##dimension of summary statistics A=matrix(0,nrow=d,ncol=d); diag(A)=1 ##matrix for distance -- equal weight for all components; alternatively choose as inverse of variance from linear model. ABC.out=matrix(0,nrow=N,ncol=4) ##matrix to store (delta,dist) pairs for(i in 1:N){ ##LOOP delta.prop=rnorm(1,0,2) alpha.prop=runif(1,0.5,0.9) beta.prop=runif(1,-0.5,0.5) ysim=rstable(n,alpha.prop,beta.prop,gamma,delta.prop) ##simulated data Ssim=S(ysim) #summaries dist=matrix(Sobs-Ssim,nrow=1) %*% A %*% matrix(Sobs-Ssim,ncol=1) ABC.out[i,]=c(alpha.prop,beta.prop,delta.prop,dist) } h=quantile(ABC.out[,4],0.02) ##threshold for ABC ABCsam=ABC.out[ABC.out[,4]