Introduction to Bayesian Statistics 2th

谢谢楼主分享哈

谢谢分享

thank y

1. MACRO
   
2. BayesMultReg.mac kpar y Xmat;
   
3. sigma sig;
   
4. prior b0 V0.
   
5. # BayesMultReg.mac calculates the joint probability distribution for the
   
6. # multiple linear regression model
   
7. 
   
8. ###############################################################
   
9. #
   
10. # Instructions for using this macro
    
11. # are contained in  "Introduction to Bayesian Statistics, 3rd Edition" by
    
12. # William M. Bolstad and James Curran (John Wiley & Sons)
    
13. #
    
14. ###############################################################
    
15. #  
    
16. #  Neither Minitab Inc. nor the author of this MACRO makes
    
17. #  any claim or offers any warranty whatsoever with regard to
    
18. #  the accuracy of this MACRO or its suitability for use, and
    
19. #  Minitab Inc. and the author of this MACRO each disclaims     
    
20. #  any liability with respect thereto.            
    
21. #
    
22. ###############################################################
    
23. 
    
24. 
    
25. # MACRO NAME: BayesMultReg.mac
    
26. # AUTHORS NAME AND ADDRESS
    
27. # William M (Bill) Bolstad
    
28. # Statistics Department
    
29. # Waikato UNIVERSITY
    
30. # Hamilton, New Zealand
    
31. 
    
32. # DATE: October 2015
    
33. 
    
34. mconstant nobs kpar isig sig var  prb0 prb1  ma0
    
35. mconstant df ssr  xbar ybar xybar  x2bar y2bar bls als ssx
    
36. mconstant pre0 pre1 preLS k1 k2 k3 k4 kint npred kpred
    
37. mcolumn y bL b0 B1 x0 xc.1-xc.kpar
    
38. mcolumn fit res
    
39. mcolumn xp yp lp up
    
40. mmatrix Xmat X XT XTX XTXinv V0 VL V1 V0inv VLinv V1inv W W1 W2 W3 W4
    
41. default sig=0
    
42. 
    
43. # nobs is number of observations, kpred is number of predictors,
    
44. # kpar (=kpred+1) is the number of parameters (including intercept)
    
45. # xmat is matrix of predictors in columns xc.2-xc.kpar.  
    
46. # xc.1 is column of 1's for intercept.   xc.1-xc.kpar is combined into Matrix X
    
47. # b0 and VL are mean vector and covariance matrix of likelihood (least squares)
    
48. # b1 and V0 are prior meanvector and covariance matrix (kpar by kpar)
    
49. # bL and V1 are posterior mean vector and covariance matrix (kpar by kpar)
    
50. # sig is standard deviation. When it is not input, estimate
    
51. # calculated from residuals is used which gives approximation.
    
52. # ia indicator for intercept prior, 0=flat prior, 1=normal prior
    
53. # ib indicator for coefficients prior, 0=joint flat prior, 1=joint normal prior
    
54. 
    
55. 
    
56. let nobs=n(y)
    
57. 
    
58. ## prepare matrix X including column for intercept
    
59. copy Xmat xc.2-xc.kpar
    
60. set xc.1
    
61. nobs(1)
    
62. end
    
63. copy xc.1-xc.kpar X
    
64. 
    
65. ## Find least squares vector (MLE) and covariance matrix
    
66. 
    
67. trans X XT
    
68. mult XT X XTX
    
69. inve XTX XTXinv
    
70. mult XTXinv XT W
    
71. mult W y bL
    
72. 
    
73. print bL
    
74. mult X bL fit
    
75. let res=y-fit
    
76. 
    
77. if (sig=0)
    
78. let df=nobs-kpar
    
79. let ssr=sum(res**2)
    
80. let sig=sqrt(ssr/df)
    
81. print df ssr sig
    
82. let var=sum(res**2)/(df)
    
83. let sig=sqrt(var)
    
84. print "Standard deviation of residuals" sig
    
85. elseif (sig > 0)
    
86. print "known standard deviation" sig
    
87. let var=sig**2
    
88. let isig=1
    
89. endif
    
90. mult var XTXinv VL
    
91. 
    
92. print bL VL
    
93. 
    
94. inve V0 V0inv
    
95. inve VL VLinv
    
96. 
    
97. add V0inv VLinv V1inv
    
98. inve V1inv V1
    
99. print V0inv VLinv V1inv V1
    
100. 
     
101. mult V1 V0inv W1
     
102. mult V1 VLinv W2
     
103. print w1 w2
     
104. 
     
105. mult W1 b0  W3
     
106. mult W2 bL W4
     
107. add W3 W4 b1
     
108. print w3 w4
     
109. print b0 V0 bL VL b1 V1
     
110. 
     
111. endmacro

复制代码