## Diff of /src/modules/bugs/distributions/DHyper.cc[f8b988] .. [ef03d6]  Maximize  Restore

### Switch to unified view

a/src/modules/bugs/distributions/DHyper.cc b/src/modules/bugs/distributions/DHyper.cc
`...`
`...`
71
```    return false;
```
71
```    return false;
```
72
```    else
```
72
```    else
```
73
```    return true;
```
73
```    return true;
```
74
```}
```
74
```}
```
75
75
76
```/* Calculates the mode of the hypergeometric distribution
```
77
78
```   We solve the equation p(x) = p(x-1) for continuous x and then take
```
79
```   the floor of x.
```
80
81
```   This reduces to solving a quadratic equation in x, by the
```
82
```   recurrence relation embodied in the function "rfunction", i.e.
```
83
```
```
84
```   p(x) = p(x - 1) * rfunction(n1, n2, m1, psi, x)
```
85
``` */
```
76
```static int modeCompute(int n1, int n2, int m1, double psi)
```
86
```static int modeCompute(int n1, int n2, int m1, double psi)
```
77
```{
```
87
```{
```
78
```    double a =  psi - 1;
```
88
```    double a =  psi - 1;
```
79
```    double b =  -((m1 + n1 + 2) * psi + n2 - m1);
```
89
```    double b =  -((m1 + n1 + 2) * psi + n2 - m1);
```
80
```    double c = psi * (n1 + 1) * (m1 + 1);
```
90
```    double c = psi * (n1 + 1) * (m1 + 1);
```
`...`
`...`
93
```    else {
```
103
```    else {
```
94
```    return static_cast<int>(q/a);
```
104
```    return static_cast<int>(q/a);
```
95
```    }
```
105
```    }
```
96
```}
```
106
```}
```
97
107
108
```/*
```
109
```   The recurrence relation
```
110
```
```
111
```   p(x) = p(x - 1) * rfunction(n1, n2, m1, psi, x)
```
112
```
```
113
```   allows us to calculate the hypergeometric probabilities
```
114
```   recursively, avoiding combinatoric problems.
```
115
``` */
```
98
```double rfunction(int n1, int n2, int m1, double psi, int i)
```
116
```double rfunction(int n1, int n2, int m1, double psi, int i)
```
99
```{
```
117
```{
```
100
```    return psi * (n1 - i + 1) * (m1 - i + 1)/(i * (n2 - m1 + i));
```
118
```    return psi * (n1 - i + 1) * (m1 - i + 1)/(i * (n2 - m1 + i));
```
101
```}
```
119
```}
```
102
120
121
```/**
```
122
``` * Returns a vector of normalized probabilities for a hypergeometric
```
123
``` * distribution. If the returned vector p is of length N, then p[0]
```
124
``` * represents the probability of the lowest possible value ll = max(0,
```
125
``` * m1 - n2) and p[N - 1] represents the probability of the largest
```
126
``` * possible value uu = min(n1, m1)
```
127
``` */
```
103
```static vector<double> density(int n1, int n2, int m1, double psi)
```
128
```static vector<double> density(int n1, int n2, int m1, double psi)
```
104
```{
```
129
```{
```
105
```    int ll = max((int) 0, m1 - n2);
```
130
```    int ll = max(0U, m1 - n2);
```
106
```    int uu = min(n1, m1);
```
131
```    int uu = min(n1, m1);
```
107
```    int N = uu - ll + 1;
```
132
```    int N = uu - ll + 1;
```
108
```    vector<double> p(N, 1);
```
133
```    vector<double> p(N);
```
109
134
135
```    // Density at mode has reference value 1
```
110
```    int mode = modeCompute(n1, n2, m1, psi);
```
136
```    int mode = modeCompute(n1, n2, m1, psi);
```
137
```    p[mode - ll] = 1;
```
111
138
112
```    // Calculate density above the mode
```
139
```    // Calculate relative density above the mode
```
113
```    if (mode < uu) {
```
140
```    if (mode < uu) {
```
114
```    double r = 1;
```
141
```    double r = 1;
```
115
```    for (int i = mode + 1; i <= uu; ++i) {
```
142
```    for (int i = mode + 1; i <= uu; ++i) {
```
116
```        r *= rfunction(n1, n2, m1, psi, i);
```
143
```        r *= rfunction(n1, n2, m1, psi, i);
```
117
```        p[i - ll] = r;
```
144
```        p[i - ll] = r;
```
118
```    }
```
145
```    }
```
119
```    }
```
146
```    }
```
120
```    // Calculate density below the node
```
147
```    // Calculate relative density below the node
```
121
```    if (mode > ll) {
```
148
```    if (mode > ll) {
```
122
```    double r = 1;
```
149
```    double r = 1;
```
123
```    for (int i = mode; i > ll; --i) {
```
150
```    for (int i = mode - 1; i >= ll; --i) {
```
124
```        r /= rfunction(n1, n2, m1, psi, i);
```
151
```        r /= rfunction(n1, n2, m1, psi, i + 1);
```
125
```        p[i - ll - 1] = r;
```
152
```        p[i - ll] = r;
```
126
```    }
```
153
```    }
```
127
```    }
```
154
```    }
```
128
```
```
155
```
```
129
```    //Normalize
```
156
```    //Normalize
```
130
```    double sump = 0;
```
157
```    double sump = 0;
```
`...`
`...`
135
```    p[i] /= sump;
```
162
```    p[i] /= sump;
```
136
```    }
```
163
```    }
```
137
```    return p;
```
164
```    return p;
```
138
```}
```
165
```}
```
139
166
140
```static int
```
167
```/*
```
141
```sampleLowToHigh(int lower_end, double ran, vector<double> const &pi)
```
168
``` * Sample from a unimodal discrete distribution given a vector of normalized
```
142
```{
```
169
``` * probabilities pi[] and an estimate of the mode. U should be drawn from a
```
143
```    int upper = pi.size() - 1;
```
170
``` * uniform(0,1) distribution.
```
144
```    if (lower_end < 0 || lower_end > upper) {
```
171
``` */
```
145
```  throw logic_error("Internal error in Hypergeometric distribution");
```
146
```    }
```
147
```    for (int i = lower_end; i < upper; ++i) {
```
148
```  ran -= pi[i];
```
149
```  if (ran <= 0) return i;
```
150
```    }
```
151
```    return upper;
```
152
```}
```
153
154
```static int
```
155
```sampleHighToLow(int upper_end, double ran, vector<double> const &pi)
```
156
```{
```
157
```    if (upper_end < 0 || upper_end >= pi.size()) {
```
158
```  throw logic_error("Internal error in Hypergeometric distribution");
```
159
```    }
```
160
161
```    for (int i = upper_end; i > 0; --i) {
```
162
```  ran -= pi[i];
```
163
```  if (ran <= 0) return i;
```
164
```    }
```
165
```    return 0;
```
166
```}
```
167
168
```static int singleDraw(int mode, vector<double> const &pi, double ran)
```
172
```static int sampleWithMode(int mode, vector<double> const &pi, double U)
```
169
```{
```
173
```{
```
170
```    if (mode == 0) return sampleLowToHigh(0, ran, pi);
```
174
```    //We try to spend our probability allocation U as quickly as
```
171
```    int max = pi.size() - 1;
```
175
```    //possible, starting with the mode.
```
172
```    if (mode == max) return sampleHighToLow(max, ran, pi);
```
173
174
```    ran -= pi[mode];
```
176
```    U -= pi[mode];
```
175
```    if (ran < 0) return mode;
```
177
```    if (U <= 0) return mode;
```
176
178
179
```    //Then we work our way outwards from the mode, one step at a time.
```
177
```    int lower = mode - 1;
```
180
```    int lower = mode - 1;
```
178
```    int upper = mode + 1;
```
181
```    int upper = mode + 1;
```
179
```    while (true) {
```
182
```    while (true) {
```
183
```  if (lower < 0) {
```
184
```      for ( ; upper < pi.size() - 1; ++upper) {
```
185
```      U -= p[upper];
```
186
```      if (U <= 0) break;
```
187
```      }
```
188
```      return upper;
```
189
```  }
```
190
```  else if (upper >= pi.size()) {
```
191
```      for ( ; lower > 0; --lower) {
```
192
```      U -= p[lower];
```
193
```      if (U <= 0) break;
```
194
```      }
```
195
```      return lower;
```
196
```  }
```
180
```    if (pi[upper] >= pi[lower]) {
```
197
```    else if (p[upper] > p[lower]) {
```
181
```        ran -= pi[upper];
```
198
```        U -= p[upper];
```
182
```        if (ran < 0) return upper;
```
199
```        if (U <= 0) return upper;
```
183
```      else if (upper == max) return sampleHighToLow(lower, ran, pi);
```
184
```        ++upper;
```
200
```        else ++upper;
```
185
```    }
```
201
```    }
```
186
```    else {
```
202
```    else {
```
187
```        ran -= pi[lower];
```
203
```        U -= p[lower];
```
188
```        if (ran < 0) return lower;
```
204
```        if (U <= 0) return lower;
```
189
```      else if (lower == 0) return sampleLowToHigh(upper, ran, pi);
```
190
```        --lower;
```
205
```        else --lower;
```
191
```    }
```
206
```    }
```
192
```    }
```
207
```    }
```
208
209
```    return mode; //-Wall
```
193
```}
```
210
```}
```
194
211
195
```double DHyper::d(double z, PDFType type,
```
212
```double DHyper::d(double z, PDFType type,
```
196
```                 vector<double const *> const &parameters,
```
213
```                 vector<double const *> const &parameters,
```
197
```         bool give_log) const
```
214
```         bool give_log) const
```
`...`
`...`
199
```    int n1,n2,m1;
```
216
```    int n1,n2,m1;
```
200
```    double psi;
```
217
```    double psi;
```
201
```    getParameters(n1, n2, m1, psi, parameters);
```
218
```    getParameters(n1, n2, m1, psi, parameters);
```
202
219
203
```    int x = static_cast<int>(z);
```
220
```    int x = static_cast<int>(z);
```
204
```    int ll = max((int) 0, m1 - n2);
```
221
```    int ll = max(0U, m1 - n2);
```
205
```    int uu = min(n1, m1);
```
222
```    int uu = min(n1, m1);
```
206
223
207
```    double den = 0;
```
224
```    double den = 0;
```
208
```    if (x >= ll && x <= uu) {
```
225
```    if (x >= ll && x <= uu) {
```
209
```    den = density(n1, n2, m1, psi)[x - ll];
```
226
```    den = density(n1, n2, m1, psi)[x - ll];
```
`...`
`...`
282
```{
```
299
```{
```
283
```    int n1,n2,m1;
```
300
```    int n1,n2,m1;
```
284
```    double psi;
```
301
```    double psi;
```
285
```    getParameters(n1, n2, m1, psi, parameters);
```
302
```    getParameters(n1, n2, m1, psi, parameters);
```
286
303
304
```    int ll = max(0U, m1 - n2);
```
287
```    int mode = modeCompute(n1, n2, m1, psi);
```
305
```    int mode = modeCompute(n1, n2, m1, psi);
```
288
```    vector<double> pi = density(n1, n2, m1, psi);
```
306
```    vector<double> pi = density(n1, n2, m1, psi);
```
289
```    return singleDraw(mode, pi, rng->uniform());
```
307
```    return sampleWithMode(mode - ll, pi, rng->uniform()) + ll;
```
290
```}
```
308
```}
```
291
309
292
```double DHyper::l(vector<double const *> const &parameters) const
```
310
```double DHyper::l(vector<double const *> const &parameters) const
```
293
```{
```
311
```{
```
294
```    int n1,n2,m1;
```
312
```    int n1,n2,m1;
```
295
```    double psi;
```
313
```    double psi;
```
296
```    getParameters(n1, n2, m1, psi, parameters);
```
314
```    getParameters(n1, n2, m1, psi, parameters);
```
297
```
```
315
```
```
298
```    return max((int) 0, m1 - n2);
```
316
```    return max(0U, m1 - n2);
```
299
```}
```
317
```}
```
300
318
301
```double DHyper::u(vector<double const *> const &parameters) const
```
319
```double DHyper::u(vector<double const *> const &parameters) const
```
302
```{
```
320
```{
```
303
```    int n1,n2,m1;
```
321
```    int n1,n2,m1;
```
`...`
`...`
317
```{
```
335
```{
```
318
```    int n1a,n2a,m1a;
```
336
```    int n1a,n2a,m1a;
```
319
```    double psia;
```
337
```    double psia;
```
320
```    getParameters(n1a, n2a, m1a, psia, para);
```
338
```    getParameters(n1a, n2a, m1a, psia, para);
```
321
339
322
```    int lla = max((int) 0, m1a - n2a);
```
340
```    int lla = max(0U, m1a - n2a);
```
323
```    int uua = min(n1a, m1a);
```
341
```    int uua = min(n1a, m1a);
```
324
342
325
```    int n1b,n2b,m1b;
```
343
```    int n1b,n2b,m1b;
```
326
```    double psib;
```
344
```    double psib;
```
327
```    getParameters(n1b, n2b, m1b, psib, para);
```
345
```    getParameters(n1b, n2b, m1b, psib, para);
```
328
346
329
```    int llb = max((int) 0, m1b - n2b);
```
347
```    int llb = max(0U, m1b - n2b);
```
330
```    int uub = min(n1b, m1b);
```
348
```    int uub = min(n1b, m1b);
```
331
349
332
```    if (lla < llb || uua > uub)
```
350
```    if (lla < llb || uua > uub)
```
333
```    return JAGS_POSINF;
```
351
```    return JAGS_POSINF;
```
334
352