Biased urn problem & Wallenius' noncentral hypergeometric distribution

General
21

In the function you can precalculate y - m in the transformed data block and pass into the function. I think the way I’ve written below may be a bit faster because of the dot_product instead of the explicit multiplication.

  real approximate_wallenus_integral(
    vector weights,
    data array[] int y,
    data array[] int m,
    data vector log_n_dt
    ) 
  {
    int num_categories = size(y);
    
    vector[num_categories] vec_y = to_vector(y);  // number of balls of each category that were drawn
    vector[num_categories] vec_y_min_m = vec_y - to_vector(m); // number of balls of each category in the urn
    
    real one_over_D = 1.0 / dot_product(weights, vec_y_min_m);
   real slice_prod = log1m_exp(one_over_D  * dot_product(log_n_dt, weights));
    
    return log_sum_exp(slice_prod * vec_y);
  }
1 Like
22

Cool!
One thing: the return type for slice_prod should be matrix and not real:

If you look at the sizes of the vectors:
log_n_dt is a vector of size num_points-1 and omegas is a vector of size num_categories.

Therefore the multiplication

(log_n_dt * weights')

produces a matrix of size (num_points-1,num_categories) which is then run through the log1m_exp function (not changing the size of the matrix) and then multiplied by vec_y (which is sized num_categories) leaving a vector of size num_points-1 which is what we log_sum_exp over…

But maybe there is a better function to do vector times row_vector multiplication?

1 Like
23

Got it. In cpp we won’t want to create that full matrix. We can loop over the num_points - 1 and keep track of a streaming log_sum_exp. But it depends on what values we need to keep track of for the derivative.

How does the function behave if you replace log1m_exp with expm1 (of course multiplying all the values by -1 before), sum, and then log that value?

24

The values of

one_over_D .* (log_n_dt * weights')

I think through sampling cross -0.693147 and so we should be using the more numerically stable log1m_exp().

If we are computing pointwise, then we can just “unroll” the function and do:

if (x > -0.693147) {
    return log(-expm1(x))
} else {
  return log1p(-exp(x))
}

which is the definition anyway… (and we could nest the definitions of expm1 and log1p to get a fully unrolled pointwise stable calculation.