Obtain logProb_

[Keith Lau]

Hello everyone,

I was trying to obtain the logProb of my response data matrix (named "response"), but I found the dimension of the "logProb_response" (100 x 7) is different from the "response" (100 x 9). My purpose is to use the  "logProb_response" to test the loo package about WAIC and LOO.  I reduced my model to become a toy example that replicates this issue (see attached file). I wonder if the nimble authors could give me some directions to solve the issue.

Thank you for this great nimble package!

Keith

[Perry de Valpine]

Hi Keith,

For a non-scalar stochastic node, it can be a bit less obvious what the size of the logProbs object will be.  A non-scalar stochastic node needs only one logProb value.  In this case, response[n,item_index[k,1]:item_index[k,2]] declares stochastic vector nodes.  From item_index, I see the second index for response will be either 1:3, 4:6, or 7:9.  Normally the logProb entry will go in the first available element of an object (logProb_response) the same shape as the random variable but with unneeded extents dropped.  In this case, column 7 of logProb_response is needed for the nodes using columns 7:9 of response, but there is no need for columns 8:9 of logProb_response, so they are dropped.  (There can still be unused space in logProb_response, but that is in the service of some simplicity, but we do drop unneeded space beyond the last needed extent in the logProb, if that makes sense.)

I hope that helps and let me know if it's not clear.

Perry

--
You received this message because you are subscribed to the Google Groups "nimble-users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to nimble-users...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/nimble-users/a2e4c399-9a83-44c5-babb-65e589f9d98bn%40googlegroups.com.

[David Pleydell]

Dear Keith

I'm not an author of nimble, but perhaps this can help...

In your model, response[n, ] is a vector that holds three realisations of a length 3 mvnorm random variable. 

For each element of a given multi-variate (length 3) sub-vector getLogProb returns the same value. e.g. 

getLogProb(cmodel,"response[1,1]") # -3.034255
getLogProb(cmodel,"response[1,2]") # -3.034255
getLogProb(cmodel,"response[1,3]") # -3.034255
getLogProb(cmodel,"response[1,4]") # -4.096456
getLogProb(cmodel,"response[1,5]") # -4.096456
getLogProb(cmodel,"response[1,6]") # -4.096456
getLogProb(cmodel,"response[1,7]") # -4.414663
getLogProb(cmodel,"response[1,8]") # -4.414663
getLogProb(cmodel,"response[1,9]") # -4.414663

getLogProb(cmodel)                            # -1372.911

Now let's compare to the values of logProb returned by the MCMC. I use "#" to indicate the output on my machine.

colnames(logProb)[c(1,101,201)]   # "logProb_response[1, 1]" "logProb_response[1, 2]" "logProb_response[1, 3]"
tail(logProb, n=1)    [c(1,101,201)]   # -3.034255  0.000000  0.000000

So we have the same value as above is position 1, and zeros in the following two columns. 

colnames(logProb)[300+c(1,101,201)]  # "logProb_response[1, 4]" "logProb_response[1, 5]" "logProb_response[1, 6]"
tail(logProb, n=1)[300+c(1,101,201)]      # -4.096456  0.000000  0.000000

Again, the same value as above is position 1, and zeros in the following two columns. 

colnames(logProb)[600+c(1,101,201)]  # "logProb_response[1, 7]" NA                       NA                      
tail(logProb, n=1)[600+c(1,101,201)]     # -4.414663        NA        NA

Again, same value as above, but here nimble has truncated the two columns of zeros (so the indexing results in NA). 

sum(tail(logProb, n=1))              # -1372.86

Close to the above value, presumably computational rounding error explains the difference. 

All the best

David

 

 

[Keith Lau]

Thank you so much for the answers, Perry! I also thank David for providing the numerical demo that fully addressed my question! And thanks for the lesson that I can use the getLogProb function to inspect the values :)

Keith

David Pleydell 在 2022年5月31日 星期二下午10:04:41 [UTC+8] 的信中寫道：