Tolle Frage für einen Neuling !!!
Ihr ABC-Algorithmus liefert Ihnen ein Beispiel aus der ABC-posterioren Verteilung. Für jede Komponente des Vektors Sie somit eine Stichprobe der Größe vom marginalen ABC-posterior. Zum Beispiel ist hier ein Spielzeugbeispiel über den normalen posterioren Mittelwert der Varianz, wenn Median und Mad als Zusammenfassungen verwendet werden: θ M.θ1, … , ΘM.θM.
#normal data with 100 observations
x=rnorm(100)
#observed summaries
sumx=c(median(x),mad(x))
#normal x gamma prior
priori=function(N){
return(cbind(rnorm(N,sd=10),1/sqrt(rgamma(N,shape=2,scale=5))))
}
ABC=function(N){
prior=priori(N) #reference table
#pseudo-data
summ=matrix(0,N,2)
for (i in 1:N){
xi=rnorm(100)*prior[i,2]+prior[i,1]
summ[i,]=c(median(xi),mad(xi)) #summaries
}
#normalisation factor for the distance
mads=c(mad(summ[,1]),mad(summ[,2]))
#distance
dist=(abs(sumx[1]-summ[,1])/mads[1])+(abs(sumx[2]-summ[,2])/mads[2])
#selection
posterior=prior[dist<quantile(dist,.05),]
return(posterior)
}
Wenn Sie planen
res=ABC(10^5);hist(res[,1])
Sie erhalten das marginale ABC-posterior für den normalen Mittelwert.
Wenn Sie jedoch eine posteriore Vorhersageprüfung durchführen möchten, können Sie nicht jeweils eine Komponente Ihres posterioren Tests generieren, um Pseudodaten und die entsprechenden Zusammenfassungen abzurufen. Sie benötigen sowohl Mittelwert als auch Varianz, um eine neue normale Stichprobe zu erhalten! Also wäre mein R-Code dann
postsample=res[sample(1:length(res[,1]),10^3),]
Um eine Probe aus dem ABC-Posterior zu ziehen, würden dann die Pseudodaten wie zuvor erzeugt:
#pseudo-data
summ=matrix(0,M,2)
for (i in 1:M){
xi=rnorm(100)*postsample[i,2]+postsample[i,1]
summ[i,]=c(median(xi),mad(xi)) #summaries
}