1
2 #ifndef MCMCNBUTIL_CC
3 #define MCMCNBUTIL_CC
4
5
6 #include "matrix.h"
7 #include "algorithm.h"
8 #include "distributions.h"
9 #include "la.h"
10 #include "smath.h"
11
12
13 #include "MCMCnbutil.h"
14
15 #define R_NO_REMAP_MATH
16 #include <R.h>
17 #include <R_ext/Utils.h>
18 #include <Rinternals.h>
19 #include <Rmath.h>
20
21 using namespace std;
22 using namespace scythe;
23
24
25 // component selector
component_selector(const int yt)26 Matrix<> component_selector(const int yt) {
27 const int nrange = 8;
28 int rtake = 7; // default to one component
29 Matrix<int> ranges(nrange, 2);
30 Matrix<int> ncomp(nrange + 1, 1);
31
32 ranges = 1, 5, 20, 50, 440, 1600, 10000, 30000;
33 ncomp = 10, 9, 4, 3, 2, 2, 2, 1;
34 for (int i = 0; i < nrange-1; i++) {
35 if (yt >= ranges[i] && yt < ranges[i + 1]) {
36 rtake = i;
37 }
38 }
39 int nr = yt - ranges[rtake];
40 int nc = ncomp[rtake];
41 Matrix<> rmat(nc, 3);
42
43 double mu = -::Rf_digamma(yt);
44 double sigma = sqrt(::Rf_trigamma(yt));
45 Matrix<double,Row,Concrete> P10(4, 10);
46 Matrix<double,Row,Concrete> M10(4, 10);
47 Matrix<double,Row,Concrete> V10(4, 10);
48
49 P10 = 0.00396984425, 0.0396244597, 0.16776747, 0.147036501, 0.125306271,
50 0.1014852, 0.103758531, 0.115972617, 0.107065659, 0.0880134482, // yt == 1
51 0.00396860865, 0.039627049, 0.16777003, 0.147037759, 0.125304523,
52 0.101481714, 0.103759705, 0.115973128, 0.107065554, 0.0880119305, // yt == 2
53 0.00396904734, 0.0396280642, 0.167770514, 0.14703607, 0.12530043,
54 0.10148242, 0.103759287, 0.115974323, 0.107066971, 0.0880128738, // yt == 3
55 0.00396883344, 0.039627504, 0.167771274, 0.147036659, 0.125301189,
56 0.101481755, 0.103760036, 0.115974339, 0.107065718, 0.0880126919; // yt == 4
57
58 M10 = 3.55887454, 2.11415904, 0.968124631, 0.51537638, 0.145465449,
59 -0.145346445, -0.416660312, -0.689002855, -0.974965634, -1.27310004,
60 2.78807754, 1.84979328, 0.94844169, 0.577673108, 0.223449219,
61 -0.0831666379, -0.387174155, -0.69613969, -1.01843553, -1.34112844,
62 2.43454312, 1.7315327, 0.942157786, 0.60208557, 0.251664821,
63 -0.0644746918, -0.379817508, -0.696781518, -1.0293035, -1.35705784,
64 2.30484474, 1.71231656, 0.952907078, 0.601128034, 0.252368847,
65 -0.059032783, -0.375605704, -0.699071542, -1.03734211, -1.3609072;
66
67 V10 = 2.62603032, 1.21263644, 0.66586521, 0.256650604, 0.120071142,
68 0.0649909219, 0.0473513798, 0.046639443, 0.0576541602, 0.0888536903,
69 2.39619753, 1.16995764, 0.688870128, 0.307084756, 0.155644328,
70 0.0899360571, 0.0707828448, 0.0751755614, 0.0990773728, 0.15471843,
71 2.16215586, 1.11062998, 0.682294453, 0.324750601, 0.173204837,
72 0.108063698, 0.0917073596, 0.100257256, 0.131371692, 0.200024832,
73 1.92939547, 1.00671896, 0.638983371, 0.322852776, 0.18445103,
74 0.122217472, 0.106400052, 0.116936918, 0.154113316, 0.233525098;
75
76 Matrix<double,Row,Concrete> P9(15, 9);
77 Matrix<double,Row,Concrete> M9(15, 9);
78 Matrix<double,Row,Concrete> V9(15, 9);
79
80 P9 = 0.0435820277, 0.167794347, 0.147040722, 0.125310654, 0.10147112,
81 0.10376347, 0.115973878, 0.107056197, 0.0880075845, // yt == 5
82 0.0435817033, 0.167795985, 0.1470426, 0.125311016, 0.101470666,
83 0.103763084, 0.115972864, 0.107055437, 0.0880066471, // yt == 6
84 0.0435798639, 0.167797087, 0.147042073, 0.125313291, 0.101470979,
85 0.103761847, 0.115973234, 0.107054351, 0.0880072753, // yt == 7
86 0.043578895, 0.167797426, 0.147041988, 0.125313875, 0.101470922,
87 0.103761581, 0.115973137, 0.107054001, 0.0880081751, // yt == 8
88 0.0435786725, 0.167797743, 0.1470428, 0.125313553, 0.101470946,
89 0.103761391, 0.115973188, 0.10705364, 0.0880080663, //yt == 9
90 0.0435779307, 0.167797788, 0.147042734, 0.125314068, 0.101471449,
91 0.10376142, 0.115973187, 0.107053473, 0.0880079505, // yt == 10
92 0.043576761, 0.167801375, 0.147042624, 0.125314075, 0.101470546,
93 0.103761069, 0.115973226, 0.107051966, 0.0880083593, // yt == 11
94 0.0435771819, 0.167801103, 0.147042441, 0.125313864, 0.101470305,
95 0.103761519, 0.11597319, 0.107052417, 0.0880079809, // yt == 12
96 0.0435778469, 0.167800486, 0.147041951, 0.125313914, 0.101470076,
97 0.103761707, 0.115973611, 0.107052756, 0.0880076518, // yt == 13
98 0.0435786417, 0.16779926, 0.147042119, 0.125313391, 0.101470554,
99 0.103762378, 0.115973792, 0.107052537, 0.0880073289, // yt == 14
100 0.043581505, 0.167797871, 0.147043608, 0.125312521, 0.101469081,
101 0.103762173, 0.115973414, 0.107054363, 0.0880054639, // yt == 15
102 0.0435811435, 0.167798952, 0.147043687, 0.125312616, 0.101468918,
103 0.103762052, 0.115973417, 0.107053968, 0.0880052462, //yt == 16
104 0.0435812603, 0.167798873, 0.147044518, 0.125312321, 0.101468879,
105 0.103761729, 0.115972692, 0.107054049, 0.0880056789, //yt == 17
106 0.0435808733, 0.167799002, 0.147044529, 0.125312675, 0.101468951,
107 0.103761472, 0.115972643, 0.107053883, 0.0880059719, // yt == 18
108 0.0435807283, 0.167799231, 0.14704464, 0.12531292, 0.101468814,
109 0.103761275, 0.115972628, 0.107053662, 0.088006103; // yt == 19
110
111 M9 = 1.31113348, 0.963928895, 0.659198795, 0.240742429, -0.108844644,
112 -0.252087404, -0.6546691, -1.04146524, -1.37874376,
113 1.25919247, 0.957217299, 0.66710982, 0.251658342, -0.125234491,
114 -0.240137829, -0.64912733, -1.03921002, -1.37439461,
115 1.21602216, 0.94778507, 0.671484869, 0.265435387, -0.104709908,
116 -0.24708343, -0.653441223, -1.04076324, -1.36988994,
117 1.18027937, 0.939725546, 0.67760436, 0.293497817, -0.110879079,
118 -0.257696481, -0.655613756, -1.0406543, -1.36465528,
119 1.14996911, 0.934206664, 0.686267712, 0.311595579, -0.112948479,
120 -0.274222612, -0.653808807, -1.04092104, -1.35962481,
121 1.12841748, 0.932206841, 0.69102714, 0.319038554, -0.109581301,
122 -0.302963892, -0.641448217, -1.03858769, -1.35274157,
123 1.10451126, 0.925180162, 0.689947194, 0.309587296, -0.123979787,
124 -0.239246368, -0.658582798, -1.03932069, -1.347407,
125 1.08624068, 0.918081034, 0.697616213, 0.330655882, -0.106424319,
126 -0.290644969, -0.644517493, -1.04099153, -1.34370607,
127 1.0671125, 0.915784215, 0.70024231, 0.330800476, -0.125598534,
128 -0.244656951, -0.661886313, -1.04447342, -1.33948264,
129 1.05721516, 0.918099637, 0.698999193, 0.325014717, -0.153165358,
130 -0.225909041, -0.659788653, -1.03711782, -1.33064663,
131 1.02150943, 0.896206397, 0.702224917, 0.344137939, -0.119895501,
132 -0.256590721, -0.641185291, -1.03810889, -1.32943558,
133 1.02508782, 0.902555642, 0.699256309, 0.336391119, -0.121902141,
134 -0.242730179, -0.6538063, -1.0385784, -1.32415888,
135 0.997274184, 0.88197491, 0.696155279, 0.3460138, -0.128981232,
136 -0.227346713, -0.630322077, -1.03647508, -1.32316505,
137 0.995086849, 0.891409674, 0.70171109, 0.341992158, -0.127906113,
138 -0.245952673, -0.638792902, -1.03392281, -1.31486719,
139 0.997741814, 0.892112396, 0.698155553, 0.337731787, -0.122630195,
140 -0.240878604, -0.651951415, -1.02898878, -1.3062535;
141
142 V9 = 1.5732832, 0.745075965, 0.340530976, 0.206325108, 0.206977107,
143 0.133034557, 0.123981078, 0.155417698, 0.247661591,
144 1.52550277, 0.745216293, 0.347702459, 0.213195645, 0.220928839,
145 0.147502243, 0.139478204, 0.17271313, 0.269719569,
146 1.48970429, 0.74910777, 0.35810967, 0.221696291, 0.216470192,
147 0.155837875, 0.148481868, 0.185394632, 0.28822907,
148 1.46105103, 0.752441091, 0.365198621, 0.220104509, 0.199190433,
149 0.167708126, 0.15761138, 0.197076001, 0.304425302,
150 1.43764551, 0.754190306, 0.367534375, 0.215776065, 0.185257157,
151 0.180142183, 0.165402413, 0.206954388, 0.318591695,
152 1.41468216, 0.75198881, 0.368357589, 0.215271168, 0.178178434,
153 0.198636491, 0.176790288, 0.218155881, 0.332156859,
154 1.39851898, 0.755429842, 0.377058085, 0.229287048, 0.214645547,
155 0.18489307, 0.178139004, 0.226237823, 0.343708183,
156 1.38111403, 0.759024378, 0.379809227, 0.222659694, 0.185443843,
157 0.206181273, 0.184773494, 0.231840962, 0.353714302,
158 1.36922993, 0.759197249, 0.381687395, 0.225704876, 0.199623554,
159 0.195711194, 0.18270427, 0.236837387, 0.363050264,
160 1.35398708, 0.753650144, 0.381381699, 0.231057971, 0.208319112,
161 0.210190241, 0.194957855, 0.249236388, 0.373774124,
162 1.35652837, 0.774748407, 0.400413698, 0.238592235, 0.199467639,
163 0.230239828, 0.19924794, 0.251600772, 0.380054821,
164 1.33546695, 0.763749521, 0.396745563, 0.241905327, 0.212176877,
165 0.218950701, 0.201882762, 0.257807637, 0.388524892,
166 1.33575722, 0.781739895, 0.417799104, 0.256728889, 0.211230256,
167 0.254750255, 0.208700024, 0.26087813, 0.393822372,
168 1.3227643, 0.771070524, 0.406631212, 0.249617029, 0.210958958,
169 0.249496089, 0.214362668, 0.270024593, 0.402615529,
170 1.30630549, 0.765952536, 0.407914566, 0.255018833, 0.226289944,
171 0.236887588, 0.221124118, 0.280039124, 0.411219814;
172
173
174 Matrix<double,Row,Concrete> coeff_p(4, nc);
175 Matrix<double,Row,Concrete> coeff_m(5, nc);
176 Matrix<double,Row,Concrete> coeff_v(5, nc);
177
178 // yt from 20 to 49
179 if (rtake == 2) {
180 coeff_p = -5.644536495326e-009, 7.299190941772e-009, -1.788056445701e-008, 9.259794020097e-009,
181 -1.266992312621e-006, 1.387196986613e-006, -2.391966642312e-006, 1.224613603301e-006,
182 1.000000000000e+000, 1.000000000000e+000, 1.000000000000e+000, 1.000000000000e+000,
183 4.730022618640e+000, 3.672627139064e+000, 4.871566292199e+000, 3.215154075256e+000;
184
185 coeff_m = 4.552797222246e-005, 2.284729919322e-005, -3.900177124794e-005, -2.486737015928e-005,
186 4.638009105861e-002, -1.095058888700e-002, 4.731686443506e-002, 2.978371498898e-002,
187 1.000000000000e+000, 1.000000000000e+000, 1.000000000000e+000, 1.000000000000e+000,
188 9.627160143020e-002, -1.690501643196e-002, -5.610095109269e-001, -4.643825308040e-002,
189 1.143772956136e+000, 1.944583776810e+000, -6.033854619021e+000, -1.105498133467e+000;
190
191 coeff_v = -2.191015160635e-005, 7.060864706965e-005, 1.823003483481e-004, 1.613752763707e-004,
192 9.939739739229e-002, 1.143203813438e-001, 1.675101633325e-001, 1.943336591437e-001,
193 1.000000000000e+000, 1.000000000000e+000, 1.000000000000e+000, 1.000000000000e+000,
194 9.208564449364e-002, 1.548617268518e-001, 2.735383307281e-001, 2.797653349940e-001,
195 7.148740127686e-001, 2.428636911969e+000, 4.861423133312e+000, 3.840341872065e+000
196 ;
197 }
198 // yt from 50 to 439
199 if (rtake == 3) {
200 coeff_p = -5.639545796991280e-010, 2.651836392450035e-010, -2.384482520627535e-011,
201 4.698743002874532e-007, -1.280380156002802e-007, -1.227680572544847e-007,
202 1.000000000000000e+000, 1.000000000000000e+000, 1.000000000000000e+000,
203 4.730482920811330e+000, 2.093982718501769e+000, 3.214956149674574e+000;
204
205 coeff_m = -1.653173201148335e-006, -8.298537364426537e-007, -1.431525987300163e-006,
206 1.036578627632170e-002, 5.017456263052972e-003, 8.386323466104712e-003,
207 1.000000000000000e+000, 1.000000000000000e+000, 1.000000000000000e+000,
208 2.349390607953272e-002, 5.123168011502721e-002, -1.841057020139425e-002,
209 1.432904956685477e+000, 6.453910704667408e+000, -1.410602407670769e+000;
210
211 coeff_v = -2.726183914412441e-007, 1.118379212729684e-006, 2.197737873275589e-006,
212 2.788507874891710e-002, 2.433214514397419e-002, 3.186581505796005e-002,
213 1.000000000000000e+000, 1.000000000000000e+000, 1.000000000000000e+000,
214 2.777086294607445e-002, 2.778340896223197e-002, 3.808382220884354e-002,
215 8.369406298984288e-001, 1.489387981224663e+000, 1.958805931276004e+000;
216 }
217
218 // yt from 440 to 1599
219 if (rtake == 4) {
220 coeff_p = 1.034981541036597e-010, -2.291586556531707e-010,
221 -2.445177000398938e-007, 5.414543692806514e-007,
222 1.000000000000000e+000, 1.000000000000000e+000,
223 1.451229377475864e+000, 3.216167113242079e+000;
224
225 coeff_m = -6.578325435644067e-008, -6.292364160498604e-008,
226 1.648723149067166e-003, 1.618047470065775e-003,
227 1.000000000000000e+000, 1.000000000000000e+000,
228 1.594968525045459e-002, -7.091699113800587e-003,
229 5.566082591106806e+000, -2.516741952410371e+000;
230
231 coeff_v = -2.802162650788337e-009, 3.776558110733883e-008,
232 4.051503597935380e-003, 5.022619018941299e-003,
233 1.000000000000000e+000, 1.000000000000000e+000,
234 4.018981069179972e-003, 5.525253413878772e-003,
235 9.654061278849895e-001, 1.450507513327352e+000;
236 }
237
238 // yt from 1600 to 10000
239 if (rtake == 5) {
240 coeff_p = -1.586037487404490e-013, 1.291237745205579e-013,
241 3.575996226727867e-009, -2.911316152726367e-009,
242 1.000000000000000e+000, 1.000000000000000e+000,
243 2.228310599179340e+000, 1.814126328168031e+000;
244
245 coeff_m = -2.419956255676409e-009, -2.419411092563945e-009,
246 3.245753451748892e-004, 3.245669014250788e-004,
247 1.000000000000000e+000, 1.000000000000000e+000,
248 1.895335618211674e-003, -2.327930564510444e-003,
249 3.388553853864067e+000, -4.162274939236667e+000;
250
251 coeff_v = -6.024563976875348e-011, 5.024317053887777e-010,
252 -6.540694956580495e-004, 8.898044793516080e-004,
253 1.000000000000000e+000, 1.000000000000000e+000,
254 -6.582951415419203e-004, 9.246987493760628e-004,
255 1.006399508694657e+000, 1.149073788967684e+000;
256 }
257
258 // yt from 10000 to 30000
259 if (rtake == 6) {
260 coeff_p = -1.663426552872397e-014, 1.354267905471566e-014,
261 1.141056828884990e-009, -9.289835742028532e-010,
262 1.000000000000000e+000, 1.000000000000000e+000,
263 2.228285989630589e+000, 1.814142639751731e+000;
264
265 coeff_m = -8.929405559006038e-011, -8.931137480031157e-011,
266 6.319814700961324e-005, 6.320244393309693e-005,
267 1.000000000000000e+000, 1.000000000000000e+000,
268 4.785131212048377e-004, -5.877524860249395e-004,
269 4.271922830906078e+000, -5.247218808668549e+000;
270
271 coeff_v = -1.418731402291282e-012, 1.576782807097003e-011,
272 -5.512224505288543e-006, 1.914006058179041e-004,
273 1.000000000000000e+000, 1.000000000000000e+000,
274 -5.638714069888806e-006, 1.959753272178233e-004,
275 1.006201804172733e+000, 1.087101027065273e+000;
276 }
277
278 // one component
279 if (rtake == 7) {
280 rmat(_, 0) = 1.0;
281 rmat(_, 1) = mu;
282 rmat(_, 2) = sigma * sigma;
283 }
284
285 if (rtake == 0) {
286 rmat(_, 0) = P10(nr, _);
287 rmat(_, 1) = M10(nr, _);
288 rmat(_, 2) = V10(nr, _);
289 rmat(_, 1) = rmat(_, 1) * sigma + mu;
290 rmat(_, 2) = rmat(_, 2) * (sigma * sigma);
291 }
292 if (rtake == 1) {
293 rmat(_, 0) = P9(nr, _);
294 rmat(_, 1) = M9(nr, _);
295 rmat(_, 2) = V9(nr, _);
296 rmat(_, 1) = rmat(_, 1) * sigma + mu;
297 rmat(_, 2) = rmat(_, 2) * (sigma * sigma);
298 }
299 if (rtake > 1 && rtake < 7) {
300 rmat(_, 0) = (1 + coeff_p(0, _) * yt * yt + coeff_p(1, _)*yt) / (coeff_p(3, _));
301 rmat(_, 1) = (1 + coeff_m(0, _) * yt * yt + coeff_m(1, _)*yt) / (coeff_m(3, _) * yt + coeff_m(4, _));
302 rmat(_, 2) = (1 + coeff_v(0, _) * yt * yt + coeff_v(1, _)*yt) / (coeff_v(3, _) * yt + coeff_v(4, _));
303 rmat(_, 1) = rmat(_, 1) * sigma + mu;
304 rmat(_, 2) = rmat(_, 2) * (sigma * sigma);
305 }
306
307 return rmat;
308
309 }
310
311
312
313
rho_conditional_log_density(const double rho,const Matrix<> & y,const Matrix<> & lambda,const double g,const double e,const double f)314 double rho_conditional_log_density(const double rho,
315 const Matrix<>& y,
316 const Matrix<>& lambda,
317 const double g,
318 const double e,
319 const double f) {
320 const int n = y.rows();
321 double logprior;
322 double loglike = 0;
323 double logfunval;
324
325 if (rho > 0) {
326 logprior = (e - 1) * log(rho) - (e+f) * log(rho + g);
327 for (int t = 0; t < n; t++) {
328 loglike += lngammafn(rho + y[t]) - lngammafn(rho) - lngammafn(y[t] +1);
329 loglike += rho * log(rho) + y[t]*log(lambda[t])- (rho + y[t]) * log(rho + lambda[t]);
330
331 }
332 double lpost = logprior + loglike;
333 logfunval = lpost;
334 } else {
335 logfunval = -numeric_limits<double>::infinity();
336 }
337 return logfunval;
338 }
339
340
341
342 #endif
343