-
Notifications
You must be signed in to change notification settings - Fork 0
/
plotPost_upper.R
113 lines (109 loc) · 4.9 KB
/
plotPost_upper.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
plotPost_upper = function( paramSampleVec , credMass=0.95 , compVal=NULL ,
HDItextPlace=0.7 , ROPE=NULL , yaxt=NULL , ylab=NULL ,
xlab=NULL , cex.lab=NULL , cex=NULL , xlim=lims1 , main=NULL ,
col=NULL , border=NULL , showMode=F , showCurve=F , breaks=NULL ,
... ) {
# Override defaults of hist function, if not specified by user:
# (additional arguments "..." are passed to the hist function)
if ( is.null(xlab) ) xlab="Parameter"
if ( is.null(cex.lab) ) cex.lab=1.5
if ( is.null(cex) ) cex=1.4
if ( is.null(xlim) ) xlim=range( c( compVal , paramSampleVec ) )
if ( is.null(main) ) main=""
if ( is.null(yaxt) ) yaxt="n"
if ( is.null(ylab) ) ylab=""
if ( is.null(col) ) col="skyblue"
if ( is.null(border) ) border="white"
postSummary = matrix( NA , nrow=1 , ncol=11 ,
dimnames=list( c( xlab ) ,
c("mean","median","mode",
"hdiMass","hdiLow","hdiHigh",
"compVal","pcGTcompVal",
"ROPElow","ROPEhigh","pcInROPE")))
postSummary[,"mean"] = mean(paramSampleVec)
postSummary[,"median"] = median(paramSampleVec)
mcmcDensity = density(paramSampleVec)
postSummary[,"mode"] = mcmcDensity$x[which.max(mcmcDensity$y)]
source("HDIofMCMC.R")
HDI = HDIofMCMC( paramSampleVec , credMass )
postSummary[,"hdiMass"]=credMass
postSummary[,"hdiLow"]=HDI[1]
postSummary[,"hdiHigh"]=HDI[2]
# Plot histogram.
if ( is.null(breaks) ) {
breaks = c( seq( from=min(paramSampleVec) , to=max(paramSampleVec) ,
by=(HDI[2]-HDI[1])/18 ) , max(paramSampleVec) )
}
if ( !showCurve ) {
par(xpd=NA)
histinfo = hist( paramSampleVec , xlab=xlab , yaxt=yaxt , ylab=ylab ,
freq=F , border=border , col=col ,
xlim=xlim , main=main , cex=cex , cex.lab=cex.lab ,
breaks=breaks , ... )
}
if ( showCurve ) {
par(xpd=NA)
histinfo = hist( paramSampleVec , plot=F )
densCurve = density( paramSampleVec , adjust=2 )
plot( densCurve$x , densCurve$y , type="l" , lwd=5 , col=col , bty="n" ,
xlim=xlim , xlab=xlab , yaxt=yaxt , ylab=ylab ,
main=main , cex=cex , cex.lab=cex.lab , ... )
}
cenTendHt = 0.9*max(histinfo$density)
cvHt = 0.7*max(histinfo$density)
ROPEtextHt = 0.55*max(histinfo$density)
# Display mean or mode:
if ( showMode==F ) {
meanParam = mean( paramSampleVec )
text( meanParam , cenTendHt ,
bquote(mean==.(signif(meanParam,3))) , adj=c(.5,0) , cex=cex )
} else {
dres = density( paramSampleVec )
modeParam = dres$x[which.max(dres$y)]
text( modeParam , cenTendHt ,
bquote(mode==.(signif(modeParam,3))) , adj=c(.5,0) , cex=cex )
}
# Display the comparison value.
# if ( !is.null( compVal ) ) {
# cvCol = "darkgreen"
# pcgtCompVal = round( 100 * sum( paramSampleVec > compVal )
# / length( paramSampleVec ) , 1 )
# pcltCompVal = 100 - pcgtCompVal
# lines( c(compVal,compVal) , c(0.96*cvHt,0) ,
# lty="dotted" , lwd=1 , col=cvCol )
# text( compVal , cvHt ,
# bquote( .(pcltCompVal)*"% < " *
# .(signif(compVal,3)) * " < "*.(pcgtCompVal)*"%" ) ,
# adj=c(pcltCompVal/100,0) , cex=0.8*cex , col=cvCol )
# postSummary[,"compVal"] = compVal
# postSummary[,"pcGTcompVal"] = ( sum( paramSampleVec > compVal )
# / length( paramSampleVec ) )
# }
# Display the ROPE.
if ( !is.null( ROPE ) ) {
ropeCol = "darkred"
pcInROPE = ( sum( paramSampleVec > ROPE[1] & paramSampleVec < ROPE[2] )
/ length( paramSampleVec ) )
lines( c(ROPE[1],ROPE[1]) , c(0.96*ROPEtextHt,0) , lty="dotted" , lwd=2 ,
col=ropeCol )
# lines( c(ROPE[2],ROPE[2]) , c(0.96*ROPEtextHt,0) , lty="dotted" , lwd=2 ,
# col=ropeCol)
text( mean(ROPE) , ROPEtextHt ,
bquote( .(round(100*pcInROPE))*"% in ROPE" ) ,
adj=c(.5,0) , cex=1 , col=ropeCol )
postSummary[,"ROPElow"]=ROPE[1]
postSummary[,"ROPEhigh"]=ROPE[2]
postSummary[,"pcInROPE"]=pcInROPE
}
# Display the HDI.
lines( HDI , c(0,0) , lwd=4 )
text( mean(HDI) , 0 , bquote(.(100*credMass) * "% HDI" ) ,
adj=c(.5,-1.7) , cex=cex )
text( HDI[1] , 0 , bquote(.(signif(HDI[1],3))) ,
adj=c(0.5,1.5) , cex=cex )
text( HDI[2] , 0 , bquote(.(signif(HDI[2],3))) ,
adj=c(0.5,1.5) , cex=cex )
par(xpd=F)
#
return( postSummary )
}