Help language development. Donate to The Perl Foundation

Algorithm::HierarchicalPAM cpan:TITSUKI last updated on 2019-12-29

src/util.c
#include <math.h>
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include "util.h"

// This implementation is based on Eq. 31 in:
// Minka, Thomas. "Estimating a Dirichlet distribution." (2000): 4.
void hpam_update_parameter(int num_parameter, int num_observed, int** n, double* sigma_alpha, double* alpha) {
  int* sum_n = (int*)malloc(sizeof(int) * num_observed);
  double sum_alpha = 0.0;

  for (int iter_i = 0; iter_i < 100; iter_i++) {
    for (int i = 0; i < num_observed; i++) {
      sum_n[i] = 0;
    }
    sum_alpha = 0.0;

    for (int k = 0; k < num_parameter; k++) {
      sum_alpha += alpha[k];
      for (int i = 0; i < num_observed; i++) {
        sum_n[i] += n[k][i];
      }
    }
    for (int k = 0; k < num_parameter; k++) {
      double numerator = 0.0;
      double denominator = 0.0;
      double prev_alpha;
      double coeff;
      for (int i = 0; i < num_observed; i++) {
        numerator += hpam_digamma(n[k][i] + alpha[k]) - hpam_digamma(alpha[k]);
      }

      for (int i = 0; i < num_observed; i++) {
        denominator += hpam_digamma(sum_n[i] + sum_alpha) - hpam_digamma(sum_alpha);
      }

      coeff = numerator / denominator;
      if (isnan(coeff)) {
        continue;
      }
      prev_alpha = alpha[k];
      alpha[k] = alpha[k] * coeff;
      sum_alpha += alpha[k] - prev_alpha;
      *sigma_alpha = sum_alpha;
    }
  }
  free(sum_n);
}

double hpam_digamma(double x) {
  double p;
  x=x+6;
  p=1/(x*x);
  p=(((0.004166666666667*p-0.003968253986254)*p+
      0.008333333333333)*p-0.083333333333333)*p;
  p=p+log(x)-0.5/x-1/(x-1)-1/(x-2)-1/(x-3)-1/(x-4)-1/(x-5)-1/(x-6);
  return p;
}

double hpam_log_sum(double log_a, double log_b) {
  double v;

  if (log_a < log_b) {
    v = log_b + log(1 + exp(log_a - log_b));
  } else {
    v = log_a + log(1 + exp(log_b - log_a));
  }
  return v;
}