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

Get latest updates about Open Source Projects, Conferences and News.

Sign up for the SourceForge newsletter:





No, thanks