* @brief This function samples a new topic for a word in a document based on * the topic counts computed on the rest of the corpus. This is the core * function in the Gibbs sampling inference algorithm for LDA. * @param topic_num The number of topics * @param topic The current topic assignment to a word * @param count_d_z The document topic counts * @param count_w_z The w
| 63 | * called by intruders. |
| 64 | **/ |
| 65 | static int32_t __lda_gibbs_sample( |
| 66 | int32_t topic_num, int32_t topic, const int32_t * count_d_z, const int32_t * count_w_z, |
| 67 | const int64_t * count_z, double alpha, double beta) |
| 68 | { |
| 69 | /* The cumulative probability distribution of the topics */ |
| 70 | double * topic_prs = new double[topic_num]; |
| 71 | if(!topic_prs) |
| 72 | throw std::runtime_error("out of memory"); |
| 73 | |
| 74 | /* Calculate topic (unnormalised) probabilities */ |
| 75 | double total_unpr = 0; |
| 76 | for (int32_t i = 0; i < topic_num; i++) { |
| 77 | int32_t nwz = count_w_z[i]; |
| 78 | int32_t ndz = count_d_z[i]; |
| 79 | int64_t nz = count_z[i]; |
| 80 | |
| 81 | /* Adjust the counts to exclude current word's contribution */ |
| 82 | if (i == topic) { |
| 83 | nwz--; |
| 84 | ndz--; |
| 85 | nz--; |
| 86 | } |
| 87 | |
| 88 | /* Compute the probability */ |
| 89 | // Note that ndz, nwz, nz are non-negative, and topic_num, alpha, and |
| 90 | // beta are positive, so the division by zero will not occure here. |
| 91 | double unpr = |
| 92 | (ndz + alpha) * (static_cast<double>(nwz) + beta) |
| 93 | / (static_cast<double>(nz) + topic_num * beta); |
| 94 | total_unpr += unpr; |
| 95 | topic_prs[i] = total_unpr; |
| 96 | } |
| 97 | |
| 98 | /* Normalise the probabilities */ |
| 99 | // Note that the division by zero will not occure here, so no need to check |
| 100 | // whether total_unpr is zero |
| 101 | for (int32_t i = 0; i < topic_num; i++) |
| 102 | topic_prs[i] /= total_unpr; |
| 103 | |
| 104 | /* Draw a topic at random */ |
| 105 | double r = drand48(); |
| 106 | int32_t retopic = 0; |
| 107 | while (true) { |
| 108 | if (retopic == topic_num - 1 || r < topic_prs[retopic]) |
| 109 | break; |
| 110 | retopic++; |
| 111 | } |
| 112 | |
| 113 | delete[] topic_prs; |
| 114 | return retopic; |
| 115 | } |
| 116 | |
| 117 | /** |
| 118 | * @brief Get the min value of an array - for parameter checking |