Skip to content

Commit

Permalink
Add files via upload
Browse files Browse the repository at this point in the history
  • Loading branch information
@Andrew9Lawson
Andrew9Lawson committed Oct 9, 2018
1 parent b71b1b6 commit 675b7bd
Show file tree
Hide file tree
Showing 4 changed files with 605 additions and 0 deletions.
49 changes: 49 additions & 0 deletions Gamma_Poisson_toy_data_Nimble.txt
View file
@@ -0,0 +1,49 @@
#### Gamma Poisson toy data example #####
library(nimble)
GPCode<-nimbleCode({
for(i in 1:N){
log(theta[i])<-a0+v[i]
lambda[i]<-theta[i]*texp[i]
y[i]~dpois(lambda[i])
v[i]~dnorm(0,tauv)
}
a0~dnorm(0,tau0)
tauv~dgamma(2,0.5)
tau0~dgamma(2,0.5)})

GPConsts<-list(N=10,texp=c(4.3,1.7,6.9,12,3.0,30.0,2,1.05,1.05,3.1,10.5))
GPData<-list(y=c(5,1,5,14,3,19,1,1,4,22))
GPInits<-list(alpha=1,beta=1,theta=rep(1.0,GPConsts$N))

GP<-nimbleModel(code=GPCode,name="GP",constants=GPConsts,
data=GPData,inits=GPInits)
#####################################################################

######################## Compile ######################################
CGP<-compileNimble(GP)
CGP$theta


mcmc.out<-nimbleMCMC(code=GPCode,constants=GPConsts,
data=GPData,inits=GPInits,nchains=2,niter=10000,summary=TRUE,WAIC=TRUE)

mcmc.out$summary

#####################################################################

GPConf<-configureMCMC(GP,print=TRUE)
GPConf$addMonitors(c("alpha","beta","theta"))
GPMCMC<-buildMCMC(GPConf)
CGPMCMC<-compileNimble(GPMCMC,project=GP)
niter<-10000
set.seed(1)
samples<-runMCMC(CGPMCMC,niter=niter)
plot(samples[,"alpha"],type="l",xlab="iteration",ylab=expression(alpha))
plot(samples[,"beta"],type="l",xlab="iteration",ylab=expression(beta))
#par(mfrow=c(5,2))
for(k in 1:10){
plot(samples[,2+k],type="l",xlab="iteration",ylab=expression(theta))
x11()}

samplesDF<-data.frame(samples)
plot(density(samplesDF$theta.1.))
314 changes: 314 additions & 0 deletions Georgia_MultiScale_CONV_Nimble.txt
View file
@@ -0,0 +1,314 @@
#### Goergia county/Ph district joint multi-scale model ########
MultiSC<-nimbleCode({

for( i in 1:18){
yph[i]~dpois(muph[i])
log(muph[i])<-log(eph[i])+log(thph[i])
thph[i]<-exp(aph0+uph[i]+vph[i])
vph[i]~dnorm(0,0.0001)
}
for( j in 1:159){
yc[j]~dpois(muc[j])
log(muc[j])<-log(ec[j])+log(thc[j])
thc[j]<-exp(ac0+uc[j]+vc[j])
vc[j]~dnorm(0,0.0001)
}
for (k in 1 :nsumph){weiph[k]<-1}
uph[1:18]~dcar_normal(adj[1:nsumph],weiph[1:nsumph],num[1:18],tauph)
tauph<-1/pow(sdph,2)
aph0~dflat()
sdph~dunif(0,100)
for (o in 1 :nsumc){weic[o]<-1}
uc[1:159]~dcar_normal(adj2[1:nsumc],weic[1:nsumc],num2[1:159],tauc)
tauc<-1/pow(sdc,2)
ac0~dflat()
sdc~dunif(0,100)
})
#####################################################
MultiSCinits<-list(uc=c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0),vc=c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0),ac0=-6.6,sdc=2.4,
uph=c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0),vph=c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0),
aph0=-1.5,sdph=0.1)

MSCdata<-list(
yph=c(24,17,13,22,24,38,8,38,24,48,5,27,30,17,8,
29,18,24) ,
yc=c(3,1,0,0,5,4,2,2,1,0,8,0,0,0,1,2,1,1,1,0,0,9,0,1,15,1,1,5,2,0,8,0,19,4,4,3,0,4,0,2,0,2,0,24,0,1,9,5,0,0,1,0,2,1,1,5,5,2,2,38,1,1,5,4,7,1,28,2,3,0,0,2,0,2,7,6,0,1,0,0,2,0,0,2,0,0,2,1,1,1,1,4,3,1,0,0,3,0,0,0,2,1,1,1,4,
7,7,1,0,4,3,2,2,3,6,1,0,0,1,0,18,3,0,0,0,6,1,1,2,1,0,
2,0,1,0,2,3,1,1,0,7,0,0,0,4,2,2,0,1,0,1,0,0,1,4,0,0,2,3))
MSConsts<-list(nsumc = 860,
num2= c(6, 5, 5, 6, 5, 6, 6, 7, 5, 7,
6, 6, 6, 4, 5, 7, 5, 6, 6, 3,
5, 6, 2, 3, 2, 5, 3, 7, 5, 4,
5, 4, 5, 8, 6, 3, 5, 7, 5, 6,
1, 7, 5, 5, 6, 6, 6, 4, 5, 3,
4, 5, 10, 5, 5, 4, 4, 5, 4, 10,
6, 4, 4, 9, 3, 6, 7, 6, 9, 7,
3, 4, 3, 3, 6, 7, 5, 6, 6, 7,
8, 4, 6, 7, 5, 5, 7, 5, 5, 4,
4, 5, 6, 7, 4, 6, 7, 6, 7, 4,
7, 7, 5, 6, 4, 3, 6, 7, 7, 6,
5, 5, 5, 4, 4, 5, 6, 3, 2, 6,
4, 5, 4, 4, 3, 8, 3, 4, 8, 7,
6, 8, 7, 6, 6, 4, 6, 7, 4, 6,
4, 6, 6, 4, 7, 5, 6, 7, 6, 6,
7, 6, 6, 5, 4, 7, 6, 7, 7
),
adj2 = c(
151, 138, 132, 113, 80, 3,
148, 86, 34, 32, 10,
148, 113, 80, 34, 1,
101, 100, 49, 47, 43, 19,
158, 150, 117, 84, 70,
127, 97, 78, 69, 68, 59,
147, 108, 78, 69, 67, 29,
115, 112, 110, 64, 57, 33, 28,
156, 142, 134, 77, 34,
137, 92, 86, 77, 37, 34, 2,
143, 111, 102, 84, 76, 39,
158, 143, 116, 87, 76, 45,
151, 148, 113, 63, 24, 20,
136, 92, 37, 35,
89, 54, 51, 25, 16,
124, 82, 54, 53, 51, 21, 15,
124, 121, 82, 81, 53,
126, 107, 102, 85, 79, 75,
135, 120, 49, 47, 30, 4,
63, 24, 13,
138, 132, 54, 53, 16,
110, 74, 71, 60, 48, 38,
155, 146,
148, 20, 13,
51, 15,
152, 130, 128, 106, 98,
146, 64, 57,
112, 64, 60, 58, 42, 33, 8,
109, 108, 97, 78, 7,
120, 118, 49, 19,
126, 75, 60, 56, 44,
148, 86, 50, 2,
110, 60, 48, 28, 8,
148, 134, 80, 77, 10, 9, 3, 2,
159, 137, 136, 101, 37, 14,
121, 94, 90,
137, 92, 35, 14, 10,
141, 126, 99, 74, 60, 56, 22,
145, 133, 111, 102, 11,
159, 156, 142, 129, 88, 46,
146,
112, 93, 69, 61, 58, 55, 28,
125, 101, 100, 65, 4,
122, 75, 67, 60, 31,
156, 153, 134, 116, 87, 12,
156, 129, 116, 96, 76, 40,
159, 135, 101, 88, 19, 4,
110, 60, 33, 22,
125, 100, 30, 19, 4,
92, 86, 32,
124, 25, 16, 15,
157, 109, 97, 90, 73,
140, 138, 132, 103, 83, 82, 81, 21, 17, 16,
132, 89, 21, 16, 15,
144, 105, 93, 61, 42,
126, 60, 38, 31,
115, 64, 27, 8,
69, 67, 60, 42, 28,
127, 97, 73, 6,
67, 58, 56, 48, 44, 38, 33, 31, 28, 22,
112, 105, 93, 64, 55, 42,
150, 149, 81, 70,
151, 95, 20, 13,
155, 146, 112, 105, 61, 57, 28, 27, 8,
136, 101, 43,
131, 117, 109, 108, 104, 70,
147, 122, 69, 60, 58, 44, 7,
154, 139, 127, 119, 69, 6,
154, 93, 78, 68, 67, 58, 42, 7, 6,
150, 149, 131, 117, 66, 62, 5,
115, 110, 22,
141, 130, 106, 99,
97, 59, 52,
141, 38, 22,
126, 122, 107, 44, 31, 18,
143, 116, 111, 96, 46, 12, 11,
142, 137, 34, 10, 9,
108, 97, 69, 29, 7, 6,
117, 107, 104, 102, 84, 18,
153, 138, 134, 103, 34, 3, 1,
150, 149, 121, 94, 83, 62, 53, 17,
124, 53, 17, 16,
158, 150, 140, 87, 81, 53,
158, 143, 117, 102, 79, 11, 5,
145, 126, 114, 102, 18,
92, 50, 32, 10, 2,
158, 153, 143, 140, 83, 45, 12,
159, 135, 129, 47, 40,
132, 95, 91, 54, 15,
157, 94, 52, 36,
151, 132, 95, 89,
86, 50, 37, 14, 10,
154, 144, 69, 61, 55, 42,
157, 149, 131, 121, 90, 81, 36,
151, 91, 89, 63,
133, 129, 123, 111, 76, 46,
109, 78, 73, 59, 52, 29, 6,
152, 133, 130, 129, 123, 26,
145, 141, 130, 126, 114, 72, 38,
125, 49, 43, 4,
159, 136, 65, 47, 43, 35, 4,
145, 85, 84, 79, 39, 18, 11,
153, 140, 138, 80, 53,
147, 117, 108, 107, 79, 66,
155, 64, 61, 55,
130, 72, 26,
147, 122, 104, 79, 75, 18,
147, 109, 104, 78, 66, 29, 7,
157, 131, 108, 97, 66, 52, 29,
115, 71, 48, 33, 22, 8,
133, 96, 76, 39, 11,
64, 61, 42, 28, 8,
151, 148, 13, 3, 1,
145, 126, 99, 85,
110, 71, 57, 8,
156, 76, 46, 45, 12,
104, 84, 79, 70, 66, 5,
128, 120, 30,
139, 68,
152, 135, 128, 118, 30, 19,
94, 81, 36, 17,
147, 107, 75, 67, 44,
133, 129, 98, 96,
82, 51, 17, 16,
100, 49, 43,
114, 99, 85, 75, 56, 38, 31, 18,
68, 59, 6,
152, 120, 118, 26,
152, 135, 123, 98, 96, 88, 46, 40,
145, 133, 106, 99, 98, 72, 26,
157, 149, 109, 94, 70, 66,
151, 138, 91, 89, 54, 53, 21, 1,
145, 130, 123, 111, 98, 96, 39,
156, 153, 80, 45, 34, 9,
152, 129, 120, 88, 47, 19,
101, 65, 35, 14,
159, 142, 77, 37, 35, 10,
140, 132, 103, 80, 53, 21, 1,
154, 144, 119, 68,
153, 138, 103, 87, 83, 53,
99, 74, 72, 38,
159, 156, 137, 77, 40, 9,
158, 87, 84, 76, 12, 11,
154, 139, 93, 55,
133, 130, 114, 102, 99, 85, 39,
155, 64, 41, 27, 23,
122, 108, 107, 104, 67, 7,
113, 34, 32, 24, 13, 3, 2,
157, 131, 94, 81, 70, 62,
158, 83, 81, 70, 62, 5,
132, 113, 95, 91, 63, 13, 1,
135, 129, 128, 120, 98, 26,
140, 134, 103, 87, 80, 45,
144, 139, 93, 69, 68,
146, 105, 64, 23,
142, 134, 116, 46, 45, 40, 9,
149, 131, 109, 94, 90, 52,
150, 143, 87, 84, 83, 12, 5,
142, 137, 101, 88, 47, 40, 35),
phc=c(18,18,18,16,12,4,3,1,15,15,12,11,18,15,17,18,
13,10,16,17,18,10,1,18,17,14,1,2,3,14,7,18,5,18,
16,13,15,10,12,14,1,4,16,9,11,14,16,5,16,15,
17,3,13,18,2,10,1,4,4,6,2,13,17,1,16,3,8,4,4,
12,1,14,4,10,10,12,15,3,12,18,13,13,11,12,10,15,
11,16,17,13,17,15,4,13,17,14,3,14,10,16,16,12,11,
3,2,14,8,3,3,1,12,2,18,10,1,11,12,14,4,14,13,8,14,
13,16,10,4,14,14,14,13,18,14,11,16,16,15,18,4,11,
10,15,12,4,10,1,3,18,13,12,18,14,11,4,2,11,13,12,16),
eph=c(26.7313582387354,17.8569786812963,
19.000833240556,24.8322642703346,35.6928573604747,
38.1880967220473,12.4232592469938,40.2840771546551,
31.6839976247491,33.1291602142528,6.65287755667638,
23.232040109711,20.3645029329909
,16.4453885396069,10.949025543461,16.922717476408,
15.6506686820585,23.9598964049923),
ec=c(0.842405862334888,0.375626926592719,0.484362270840443,0.199184019993243,2.11970621276708,0.735452296043789,2.64537759886506,4.0780208537788,0.813194081624956,0.782106745171208,7.27574961919763,0.564870501143738,0.728747184259648
,0.767430521475849,1.29108568514923,2.82946339512059,1.08730655358364,1.04852943858025,0.286116028719108,2.11506421230113,0.477938492417873,4.76283314473956,2.80606583721648,0.501617383683548,11.1838451225867,0.633281396899421
,1.24508768053215,8.19055193324381,4.87414737813503,0.155530460055929,12.4232592469938,0.325830921594408,30.6655708558571,1.84643774089311,2.05199865041532,4.71650691786731,0.762178965393165,4.94096405150847,0.604304060657466
,1.03286854811939,0.749847186377576,0.893889867502633,1.34172569023226,31.6839976247491,0.914380314003821,0.544098721280977,4.4863762281011,5.02728650461759,0.566933612461935,0.19229135263472,2.09410487686399,0.980352987292544
,1.03591632620309,0.527406275161016,1.0134096572773,4.75139225470228,4.4079779980096,6.18300395395692,1.00590743430203,38.1880967220473,1.25451234814483,0.123364679049487,3.34585078028669,2.30116623098126,1.13846233649622
,0.733904962555141,32.8594552982921,1.82763529456135,7.54559520183913,0.460026935064432,1.31612435432918,1.25605968163347,1.09574655443081,0.529375608692023,7.47905986182727,5.80264124911106,0.465888046763857,2.32287578871593
,0.603272504998367,0.601115615892979,0.791625190571074,0.406995596408039,0.449476934005468,1.23013012347522,0.769446743900452,0.349931812902442,2.19008644205377,1.40258747411909,2.89529540172852,0.393772928414137,0.51240182921049
,4.49134645082222,1.12181677926985,1.00890832349214,0.522248496865523,0.653396732251846,1.28062946187746,0.333473811250458,1.06672232929526,0.289070029015618,1.11773744552705,1.09851299915294,0.420593375550704,0.79767385784488
,1.90162596865489,8.57363419391819,3.82257016147108,1.35696458065077,0.635672730472786,4.96722183192189,1.15651456053045,1.30215146403775,0.783982300915024,0.73850007412749,1.88807507840582,0.461246046297912,0.925868092934693
,0.115674900499842,0.750738075355888,0.343742478947849,9.20264870150043,3.60205169489193,0.189478019018996,0.719088072178996,0.434566265842132,2.85487717544929,1.17165967316176,0.233553578998668,1.54137859916146,0.308857142112875
,8.89013422568712E-02,1.07816321933254,0.421296708954635,0.605476282997351,0.513433384869588,2.06259554036788,1.88390196687583,1.25545012601673,0.47512515880215,0.330332255379566,2.86964717693184,0.440755599796724,0.489942049178295
,0.919350536724933,1.31781235449861,2.97177118718262,3.3732338941464,1.66994794540004,0.293243140545608,0.987526988012639,1.32217302160298,0.109063566502892,0.30890403100647,1.10634344438337,4.19472730993774,0.407511374237588
,0.496225160920078,0.477844714630683,1.03193077024748),
nsumph= 88,
num = c(3, 4, 5, 4, 6, 7, 4, 8, 6, 8,
5, 6, 5, 5, 4, 2, 2, 4),
adj = c(
10, 5, 2,
6, 5, 4, 1,
13, 12, 9, 8, 4,
8, 6, 3, 2,
10, 9, 8, 6, 2, 1,
10, 9, 8, 7, 5, 4, 2,
10, 9, 8, 6,
12, 10, 9, 7, 6, 5, 4, 3,
10, 8, 7, 6, 5, 3,
14, 12, 9, 8, 7, 6, 5, 1,
18, 15, 14, 13, 12,
14, 13, 11, 10, 8, 3,
18, 17, 12, 11, 3,
16, 15, 12, 11, 10,
18, 16, 14, 11,
15, 14,
18, 13,
17, 15, 13, 11))


############################################################
MSC<-nimbleModel(code=MultiSC,name="MSC",constants=MSConsts,
data=MSCdata,inits=MultiSCinits)
########
CMSC<-compileNimble(MSC)
################## Running with monitors and CODA output #####
LSTConf<-configureMCMC(MSC,print=TRUE)
LSTConf$addMonitors(c("ac0","uph","vph","uc","vc","thph","thc"))
LSTMCMC<-buildMCMC(LSTConf)
CLSTMCMC<-compileNimble(LSTMCMC,project=MSC)
niter<-20000
set.seed(1)
samples<-runMCMC(CLSTMCMC,niter=niter,nburnin=15000,nchains=1,summary=TRUE,WAIC=TRUE,samplesAsCodaMCMC=TRUE)
#samples$summary
samples$WAIC


######################################################################
sampMAT<-as.matrix(samples$samples) ### dim is 5000 by 535
par(mfrow=c(2,2))
plot(sampMAT[,"ac0"],type="l",xlab="iteration",ylab=expression(ac0))
plot(sampMAT[,"aph0"],type="l",xlab="iteration",ylab=expression(aph0))
plot(sampMAT[,"sdc"],type="l",xlab="iteration",ylab=expression(sdc))
plot(sampMAT[,"sdph"],type="l",xlab="iteration",ylab=expression(sdph))


setwd("C:/A_NEW_FOLDER/A_NEW_FOLDER/A_PROGRAMS/R_program/NIMBLE examples/CAR and other models on NIMBLE/Multi_scale/")
library(maptools)
library(sp)
source("fillmap.R")
GCpoly<-readSplus("Georgia_Counties_SPlus.txt")
plot(GCpoly)
GHDpoly<-readSplus("GeorgiaHDistrictsSPlus.txt")
plot(GHDpoly)

fillmap(GCpoly,"",,n.col=5);x11()
fillmap(GCpoly,"",,n.col=5);x11()
fillmap(GCpoly,"",,n.col=5)

0 comments on commit 675b7bd

Please sign in to comment.