NEFI_pheno/plotDiurnalParametersExample.R at master · EcoForecast/NEFI_pheno · GitHub

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
114
library("ecoforecastR")
library("rjags")

diurnalExp <- function(a,c,k,xseq){
  k <- round(k,digits=1)
  #print(k)
  bk <- which(round(xseq,digits=1)==k)
  #print(bk)
  left <- -a*exp(-1*(xseq[1:bk]-k))+c
  right.xseq <- xseq[(bk+1):length(xseq)]
  right <- -a*exp((right.xseq-k))+c
  #print(length(c(left,right)))
  return(c(left,right))
}

siteName <- "russellSage"
xseq <- seq(0,25,0.1)

outputFileName <- paste(siteName,"_diurnalPosteriorExampleALL.pdf",sep="")
pdf(file=outputFileName,width=10,height=5)
par(mfrow=c(2,4),mai=c(0.4,0.4,0.2,0.2))

#fileName <- "GOES_Diurnal_russellSage_2017186.csv"

load("russellSage_231_varBurn_ALLParamaters.RData")
out.mat <- as.matrix(var.burn)
colnames(out.mat)
a <- out.mat[,1]
alp <- out.mat[,2]
bet <- out.mat[,3]
c <- out.mat[,4]
k <- out.mat[,5]
p.cloud <- out.mat[,6]
prec <- out.mat[,7]

h = hist(a,plot=FALSE)
h$density = h$counts/sum(h$counts)
plot(h,freq=FALSE,main="Posterior Density of a",xlab="",ylab="",cex.axis=1.5,cex.main=1.5)
qnts <- quantile(a,c(0.025,0.975))
polygon(c(qnts[1],qnts[1],qnts[2],qnts[2]),c(0,1,1,0),border=NA,col=rgb(0,0,1,0.3))
abline(v=mean(a),col="red",lwd=2)
print("a")
print(qnts)
print(mean(a))
print(qnts[2]-qnts[1])


h = hist(c,plot=FALSE)
h$density = h$counts/sum(h$counts)
plot(h,freq=FALSE,main="Posterior Density of c",xlab="",ylab="",cex.axis=1.5,cex.main=1.5)
qnts <- quantile(c,c(0.025,0.975))
polygon(c(qnts[1],qnts[1],qnts[2],qnts[2]),c(0,1,1,0),border=NA,col=rgb(0,0,1,0.3))
abline(v=mean(c),col="red",lwd=2)
print("c")
print(qnts)
print(mean(c))
print(qnts[2]-qnts[1])

h = hist(k,plot=FALSE)
h$density = h$counts/sum(h$counts)
plot(h,freq=FALSE,main="Posterior Density of k",xlab="",ylab="",cex.axis=1.5,cex.main=1.5)
qnts <- quantile(k,c(0.025,0.975))
polygon(c(qnts[1],qnts[1],qnts[2],qnts[2]),c(0,1,1,0),border=NA,col=rgb(0,0,1,0.3))
abline(v=mean(k),col="red",lwd=2)
print("k")
print(qnts)
print(mean(k))
print(qnts[2]-qnts[1])

h = hist(prec,plot=FALSE)
h$density = h$counts/sum(h$counts)
plot(h,freq=FALSE,main="Posterior Density of prec",xlab="",ylab="",cex.axis=1.5,cex.main=1.5)
qnts <- quantile(prec,c(0.025,0.975))
polygon(c(qnts[1],qnts[1],qnts[2],qnts[2]),c(0,1,1,0),border=NA,col=rgb(0,0,1,0.3))
abline(v=mean(prec),col="red",lwd=2)
print("prec")
print(qnts)
print(mean(prec))
print(qnts[2]-qnts[1])

h = hist(alp,plot=FALSE)
h$density = h$counts/sum(h$counts)
plot(h,freq=FALSE,main="Posterior Density of alp",xlab="",ylab="",cex.axis=1.5,cex.main=1.5)
qnts <- quantile(alp,c(0.025,0.975))
polygon(c(qnts[1],qnts[1],qnts[2],qnts[2]),c(0,1,1,0),border=NA,col=rgb(0,0,1,0.3))
abline(v=mean(alp),col="red",lwd=2)
print("alp")
print(qnts)
print(mean(alp))
print(qnts[2]-qnts[1])

h = hist(bet,plot=FALSE)
h$density = h$counts/sum(h$counts)
plot(h,freq=FALSE,main="Posterior Density of bet",xlab="",ylab="",cex.axis=1.5,cex.main=1.5)
qnts <- quantile(bet,c(0.025,0.975))
polygon(c(qnts[1],qnts[1],qnts[2],qnts[2]),c(0,1,1,0),border=NA,col=rgb(0,0,1,0.3))
abline(v=mean(bet),col="red",lwd=2)
print("bet")
print(qnts)
print(mean(bet))
print(qnts[2]-qnts[1])

h = hist(p.cloud,plot=FALSE)
h$density = h$counts/sum(h$counts)
plot(h,freq=FALSE,main="Posterior Density of p.cloud",xlab="",ylab="",cex.axis=1.5,cex.main=1.5)
qnts <- quantile(p.cloud,c(0.025,0.975))
polygon(c(qnts[1],qnts[1],qnts[2],qnts[2]),c(0,1,1,0),border=NA,col=rgb(0,0,1,0.3))
abline(v=mean(p.cloud),col="red",lwd=2)
print("p.cloud")
print(qnts)
print(mean(p.cloud))
print(qnts[2]-qnts[1])

dev.off()