Simple changepoint detection from BNPy HMM code

bnpy_changepoint.py

## Using notation from HDP-HMM paper in NeurIPS 2015 (Hughes, Stephenson, Sudderth)

## Start for a single sequence of interest,
## Compute the following as in the tutorial here:
https://bnpy.readthedocs.io/en/latest/examples/08_mocap6/plot-03-demo=interpret_hdphmm_params_and_run_viterbi.html#sphx-glr-examples-08-mocap6-plot-03-demo-interpret-hdphmm-params-and-run-viterbi-py

# start_prob_K : 1D array size K, sums to one
# trans_proba_KK : 2D array, size K x K, rows sum to one
# log_lik_seq_TK : 2D array, size T x K
#     Entry t,k represents log likelihood of observation at time t under state k
#     p( x[t] | z[t] = k)

## Get the length of the current sequence
T = log_lik_seq_TK.shape[0]

## Using the forward backward algorithm,
## Compute the pair-wise responsibilities and the single timestep responsibilities
## 
## For method details, see
## https://github.com/bnpy/bnpy/blob/master/bnpy/allocmodel/hmm/HMMUtil.py#L142

resp_TK, respPair_TKK, log_evidence = bnpy.allocmodel.hmm.HMMUtil.FwdBwdAlg(
	start_prob_K,
	trans_prob_KK,
    log_lik_seq_TK)


## Compute for each timestep, the probability that it is a change point:
## p( z_t != z_{t-1} | x_1, ... x_T)
##
## Remember that this is undefined at index 0.
## we can only start asking this question for t=1, 2, 3, ....
## By default, we'll just say the probability of changepoint at t=0 is 0.0

proba_changepoint_T = np.zeros(T)
for t in range(1, T):
	# Get joint probability table (K x K array) for timestep t and t-1
	joint_t_KK = respPair_TKK[t]

	# Compute sum of off-diagonal entries of this table
	proba_changepoint_T[t] = 1.0 - np.sum(np.diag(joint_t_KK))

## Now, use proba_changepoint_T as a per-timestep real-valued scalar signal
## All timesteps t with proba value above a threshold would be "changepoints"