Skip to content

Commit 173f463

Browse files
committed
Center and Identity 1D works
Basic check and it looks like it's working as expected.
1 parent 4290e7d commit 173f463

1 file changed

Lines changed: 13 additions & 6 deletions

File tree

nimbleQuad/R/Laplace.R

Lines changed: 13 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -758,7 +758,8 @@ buildOneAGHQuad1D <- nimbleFunction(
758758
}
759759
}
760760
## Given all the saved values, weights and log density, do quadrature sum.
761-
res <- log(ans) + saved_inner_max_value - 0.5 * saved_inner_logdetNegHess
761+
res <- log(ans) + saved_inner_max_value
762+
if(quadTransform_ != "center") res <- res - 0.5 * saved_inner_logdetNegHess
762763
quadrature_previous_p <<- p ## Cache this to make sure you have it for
763764
return(res)
764765
returnType(double())
@@ -807,7 +808,10 @@ buildOneAGHQuad1D <- nimbleFunction(
807808
}
808809
## Sum gradient of each node.
809810
grp_AGHQuad_sum <- gr_AGHQuad_nodes(p = p, method = 2)
810-
AGHQuad_saved_gr <<- grp_AGHQuad_sum - 0.5 * (grlogdetNegHesswrtp + grlogdetNegHesswrtre * gr_rehatwrtp)
811+
if(quadTransform_ == "center")
812+
AGHQuad_saved_gr <<- grp_AGHQuad_sum
813+
else
814+
AGHQuad_saved_gr <<- grp_AGHQuad_sum - 0.5 * (grlogdetNegHesswrtp + grlogdetNegHesswrtre * gr_rehatwrtp)
811815
}
812816
# N.B. An extra negation is built into gr_logdet because this is gradient of hessian, but the uptri_Omega_invNegHess is from the negative Hessian.
813817

@@ -859,7 +863,7 @@ buildOneAGHQuad1D <- nimbleFunction(
859863
nodes <<- quadGrid$nodes(0)
860864
wgts <<- quadGrid$weights(0)
861865
logDensity_quad <<- numeric(value = 0, length = nQ)
862-
for(i in 1:nQ) logDensity_quad[i] <- logLik_RE(reTransform = nodes[i,])
866+
for(i in 1:nQ) logDensity_quad[i] <<- logLik_RE(reTransform = nodes[i,])
863867
maxLogDens <- max(logDensity_quad)
864868
res <- log(sum(wgts*exp(logDensity_quad - maxLogDens))) + maxLogDens
865869
quadrature_previous_p <<- p ## Cache this to make sure you have it for later
@@ -872,8 +876,12 @@ buildOneAGHQuad1D <- nimbleFunction(
872876
if(any(p != quadrature_previous_p)){
873877
margLogLik_saved_value <<- calcLogLik_identity(p)
874878
}
879+
nQ <- quadGrid$gridSize()
875880
gr_wgted_wrt_p <- numeric(value = 0, length = dim(p)[1])
876-
for(i in 1:nQ) gr_wgted_wrt_p <- gr_wgted_wrt_p + gr_P_RE_b(p, nodes[i,])*exp(logDensity_quad[i])*wgts[i]
881+
for(i in 1:nQ) {
882+
gr_jointlogLikwrtp <- gr_P_RE_b(p, nodes[i,])[p_indices]
883+
gr_wgted_wrt_p <- gr_wgted_wrt_p + gr_jointlogLikwrtp*exp(logDensity_quad[i])*wgts[i]
884+
}
877885
return(gr_wgted_wrt_p/exp(margLogLik_saved_value))
878886
returnType(double(1))
879887
},
@@ -2233,8 +2241,7 @@ buildAGHQ <- nimbleFunction(
22332241
## stop("updateSettings: `computeMethod` must be 1, 2, or 3")
22342242
## }
22352243
if(quadTransform != "NULL") {
2236-
if(quadTransform == "centre") quadTransform <- "center" ## British spelling okay.
2237-
if(!any(quadTransform == c("spectral", "cholesky", "identity", "center")))
2244+
if(quadTransform != "spectral" & quadTransform != "cholesky" & quadTransform != "identity" & quadTransform != "center")
22382245
stop("`quadTransform` must be either cholesky or spectral.")
22392246
}
22402247
# actions

0 commit comments

Comments
 (0)