|
Dear Bugs users,
I am fitting a model with many observations. To avoid some computational problems, I fit the model on the first 500 observations with vague priors, then I fit the 500 next observations using the posteriors from the first adjustment as priors for the second one, and so on. The problem is that my model contains a multivariate normal distribution, for which I use a Wishart distribution as prior distribution for the variance-covariance matrix.
So I have to estimate the degrees of freedom from estimates for specifying the wishart prior distribution in the next adjustment.
I have a practical solution for that :
Note Vhat_i the estimate of the variance covariance matrix at adjustment i, for example :
[,1] [,2] [,3] [,4]
[1,] 0.4362 0.3054 0.2037 0.1481
[2,] 0.3054 0.6048 0.2859 0.3447
[3,] 0.2037 0.2859 0.4172 0.2531
[4,] 0.1481 0.3447 0.2531 0.6900
I simulate some variance covariance matrix from wishart distributions with different df (n):
model {
P1[1:4,1:4] ~ dwish(X1[1:4,1:4], n[1])
P2[1:4,1:4] ~ dwish(X2[1:4,1:4], n[2])
P3[1:4,1:4] ~ dwish(X3[1:4,1:4], n[3])
P4[1:4,1:4] ~ dwish(X4[1:4,1:4], n[4])
Xest1[1:4,1:4]<-inverse(P1[1:4,1:4])
Xest2[1:4,1:4]<-inverse(P2[1:4,1:4])
Xest3[1:4,1:4]<-inverse(P3[1:4,1:4])
Xest4[1:4,1:4]<-inverse(P4[1:4,1:4])
}
data (X_j = n_j * Vhat_i) :
list(X1= structure(.Data= c(43.62, 30.54, 20.37, 14.81, 30.54, 60.48, 28.59, 34.47, 20.37, 28.59, 41.72, 25.31, 14.81, 34.47, 25.31, 69), .Dim=c(4, 4)),
X2= structure(.Data= c(87.24, 61.08, 40.74, 29.62, 61.08, 120.96, 57.18, 68.94, 40.74, 57.18, 83.44, 50.62, 29.62, 68.94, 50.62, 138), .Dim=c(4, 4)),
X3= structure(.Data= c(130.86, 91.62, 61.11, 44.43, 91.62, 181.44, 85.77, 103.41, 61.11, 85.77, 125.16, 75.93, 44.43, 103.41, 75.93, 207), .Dim=c(4, 4)),
X4= structure(.Data= c(174.48, 122.16, 81.48, 59.24, 122.16, 241.92, 114.36, 137.88, 81.48, 114.36, 166.88, 101.24, 59.24, 137.88, 101.24, 276), .Dim=c(4, 4)),
n=c(100, 200, 300, 400))
Then I fit a linear regression of the n_j on the posterior precision of each element of the Xest_j (j=1,..,4).
Finally I am using the regression parameters to predict n from the precision of the estimates of the actual Vhat_i.
The procedure works fine, I have written a R function to do this, so if you are interested in you can email me (I have found some messages on that subject in the list but no answer).
From the first fit, I get the following results:
i j mean sd var prec intercept slope n
1 1 1 0.4362 0.07505 Vhat[1,1] 177.54098 18.657363 0.3793318 86
2 1 2 0.3054 0.05031 Vhat[1,2] 395.08575 25.971650 0.3523914 165
3 1 3 0.2037 0.03925 Vhat[1,3] 649.11355 8.838502 0.2244511 155
4 1 4 0.1481 0.04344 Vhat[1,4] 529.93227 17.135482 0.3207961 187
5 2 2 0.6048 0.10250 Vhat[2,2] 95.18144 24.107450 0.7260863 93
6 2 3 0.2859 0.04988 Vhat[2,3] 401.92693 14.076226 0.3342414 148
7 2 4 0.3447 0.06759 Vhat[2,4] 218.89463 11.448776 0.5384537 129
8 3 3 0.4172 0.07354 Vhat[3,3] 184.90674 14.683878 0.3491776 79
9 3 4 0.2531 0.05085 Vhat[3,4] 386.73910 19.092368 0.3511461 155
10 4 4 0.6900 0.12100 Vhat[4,4] 68.30135 17.196565 0.9527488 82
(n = round ( intercept + slope * prec ))
As you can see, n varies according to the element of Vhat_1 considered. In particular, n is quite lower for the variance elements.
I have choose to take the minimum n as degree of freedom for the next fit, but this is not very satisfactory because some other elements have a much bigger precision.
I would be interested in your comments and any idea on how I can circumvent this problem.
Thanks in advance,
Regards,
Etienne
-------------------------------------------------------------------
This list is for discussion of modelling issues and the BUGS software.
For help with crashes and error messages, first mail [log in to unmask]
To mail the BUGS list, mail to [log in to unmask]
Before mailing, please check the archive at www.jiscmail.ac.uk/lists/bugs.html
Please do not mail attachments to the list.
To leave the BUGS list, send LEAVE BUGS to [log in to unmask]
If this fails, mail [log in to unmask], NOT the whole list
|
