TOST for TERGM sampling

[AdrienLeGuillou]

Following up on this discussion I played around the idea of the Two One Sided Student Test (TOST) of the Hamming distance between two time steps to check if the model as converged.

Below is a simple implementation I came up to.

The main change compared to the code in tergm/src/MCMCDyn.c is that here I do not estimate the variance.
The below implementation only counts when a transition actually occurs. Therefore, the value passed to ham_advance would be -1 or +1 but never 0. The reasoning being that otherwise a model could theoretically be seen as converging just because no transition ever occurs.
The side benefit is that the variance that we need to calculate the TOST is 1 when the model has converged and <1 otherwise. It then makes sense to just use 1 always I think.

Are my assumptions there correct or am I missing an important point?

// hamming_test.h
#ifndef INCLUDE_INCLUDE_HAMMING_TEST_H_
#define INCLUDE_INCLUDE_HAMMING_TEST_H_

#include <stddef.h>
#include <stdint.h>
#include <stdbool.h>

typedef struct {
    double sum_x;           // sum of increments
    double sum_weights;     // sum of weights
    double step_decay;      // per step decay
    double tost_thresh;     // threshold for the TOST
} hamming_t;

hamming_t ham_new(size_t min_mh_steps);
void ham_advance(hamming_t *ham, int32_t new_discord);
bool ham_has_converged(const hamming_t *ham);

#endif  // INCLUDE_INCLUDE_HAMMING_TEST_H_

// hamming_test.c
#include "hamming_test.h"
#include <math.h>
#include <stdio.h>

#define SQUARE(x) x * x

hamming_t ham_new(size_t min_mh_steps) {
    hamming_t ham = {
        .sum_x = 0,
        .sum_weights = 0,
        .step_decay = 1.0 - 1.0 / min_mh_steps,
        .tost_thresh = 0.5,
    };
    return ham;
}

void ham_advance(hamming_t *ham, int32_t x) {
    ham->sum_x *= ham->step_decay;
    ham->sum_x += x;
    ham->sum_weights *= ham->step_decay;
    ham->sum_weights++;
}

bool ham_has_converged(const hamming_t *ham) {
    double mean_x = ham->sum_x / ham->sum_weights;
    // at equilibrium (what we want), the `true` variance is 1
    //   same number of -1 and +1 (don't advance on no transition)
    // before equilibrium the variance is < 1 (inbalance of -1 or +1)
    // so testing with var == 1 is ok
    double sem_x = sqrt(1 / ham->sum_weights);
    double zi_lo = (mean_x - (0.0 - ham->tost_thresh)) / sem_x;
    double zi_hi = (mean_x - (0.0 + ham->tost_thresh)) / sem_x;
    double t_crit = 1.64; // qnorm(0.95)
    bool equivalence = zi_hi < -t_crit && zi_lo > t_crit;

    return equivalence;
}

Then, with such an implementation, we would need to define the p-value and threshold to use it within the MCMCDyn1Step