Adapted from `walker_ProbSampleReplace` Index version
| 147 | // Adapted from `walker_ProbSampleReplace` |
| 148 | // Index version |
| 149 | inline Vector<INTSXP> WalkerSample(const Vector<REALSXP>& p, int n, int nans, bool one_based) |
| 150 | { |
| 151 | Vector<INTSXP> a = no_init(n), ans = no_init(nans); |
| 152 | int i, j, k; |
| 153 | std::vector<double> q(n); |
| 154 | double rU; |
| 155 | |
| 156 | std::vector<int> HL(n); |
| 157 | std::vector<int>::iterator H, L; |
| 158 | |
| 159 | int adj = one_based ? 1 : 0; |
| 160 | |
| 161 | H = HL.begin() - 1; L = HL.begin() + n; |
| 162 | for (i = 0; i < n; i++) { |
| 163 | q[i] = p[i] * n; |
| 164 | if (q[i] < 1.0) { |
| 165 | *++H = i; |
| 166 | } else { |
| 167 | *--L = i; |
| 168 | } |
| 169 | } |
| 170 | |
| 171 | if (H >= HL.begin() && L < HL.begin() + n) { |
| 172 | for (k = 0; k < n - 1; k++) { |
| 173 | i = HL[k]; |
| 174 | j = *L; |
| 175 | a[i] = j; |
| 176 | q[j] += q[i] - 1; |
| 177 | |
| 178 | L += (q[j] < 1.0); |
| 179 | |
| 180 | if (L >= HL.begin() + n) { |
| 181 | break; |
| 182 | } |
| 183 | } |
| 184 | } |
| 185 | |
| 186 | for (i = 0; i < n; i++) { |
| 187 | q[i] += i; |
| 188 | } |
| 189 | |
| 190 | for (i = 0; i < nans; i++) { |
| 191 | rU = unif_rand() * n; |
| 192 | k = static_cast<int>(rU); |
| 193 | ans[i] = (rU < q[k]) ? k + adj : a[k] + adj; |
| 194 | } |
| 195 | |
| 196 | return ans; |
| 197 | } |
| 198 | |
| 199 | // Element version |
| 200 | template <int RTYPE> |