a kind of relaxed version of a changepoint model, with an added assumption about common cause in stacked observations

gchangepoint_lts.stan

/*
NOTE: changepoint model for vector timeseries, plus "linear transform stack" feature.
Originally designed for the lellouch data (where each vector is actually a stack of 3 related vectors).

So imagine, for example, that the "inherent" thing you're observing generates a spectrogram of size 256,
but it's actually observed from 3 different vantage points, to create an obsn of size 768. Each different vantage point we'll assume
applies a linear transform to the underlying data, (hobs * x + bgspec), where the "hobs" and "bgspec" are unchanging in time
but are independently specified per bin and per vantage point.

cd ~/dev/stan/stan
make ~/git/stored_docs/dev/stan/gchangepoint/gchangepoint_lts

cd ~/git/stored_docs/dev/stan/gchangepoint/
./gchangepoint_lts --data=data_llreshaped_lts_spectros.data.R
./gchangepoint_lts --data=data_llreshaped_lts_spectros.data.R --init=init_llreshaped_lts_spectros.data.R
mv samples.csv samples_lellouch_lts.csv
~/dev/stan/stan/bin/print samples_lellouch_lts.csv | grep changeness | cut -c 25-29 > inferred_changeness_llreshaped_lts.csv
xdg-open inferred_changeness_llreshaped_lts.csv
*/
data {
	int<lower=1> D; // num dims in an UN-STACKED observation
	int<lower=0> N; // num observations
	int<lower=1> Nstack; // how many times a vector is transformed and stacked to create an observation
	vector[D] x[N, Nstack]; // obsns

	real<lower=0> beta_a; // param for beta prior on changeness
	real<lower=0> beta_b;  // ditto
	// - NB in our case we want to encourage sparsity in changeness, and have just a few coefs at the upper end, hence asymmetric
	// beta possibilities: (0.05, 0.5) v strong bimodal poss too harsh, (0.6, 0.8) gentler but still good lop
}
parameters {
	// here's the crossfade strength per transition
	real<lower=0,upper=1> changeness[N-1];
	vector[D] newgen_means[N];
	vector<lower=0>[D] newgen_stdvs[N];

/*
	// LTS-specific: for each stack, we need to know the fixed background spectrum, and the multiplier for the obsn
	vector<lower=0>[D] bgspec[Nstack];
	vector<lower=0,upper=2>[D] hobs[Nstack];
*/
}
model {
	vector[D] gen_means[N];
	vector[D] gen_stdvs[N];
	vector[D] gen_means_stacked[N, Nstack];
	vector[D] gen_stdvs_stacked[N, Nstack];
	int stackpos;

	changeness ~ beta(beta_a, beta_b);

	// each datum is associated with its own generating gaussian
	// Common prior for the brandnew (would be nice to write this without looping but stan says no)
	for (n in 1:N){
		newgen_means[n] ~ normal(0.0, 100.0);
		newgen_stdvs[n] ~ exponential(2.0);
	}

	// each datum's actual gaussian is some crossfaded mix
	// x[1] is special, it simply comes from the brandnew
	gen_means[1] <- newgen_means[1];
	gen_stdvs[1] <- newgen_stdvs[1];
	for (whichstack in 1:Nstack){
		x[1, whichstack] ~ normal(newgen_means[1], newgen_stdvs[1]);
	}
	// x[2:] is sampled from a gaussian with parameters crossfaded from the prev and the brandnew
	for (n in 2:N) {
		gen_means[n] <- (changeness[n-1] * newgen_means[n]) + ((1.0 - changeness[n-1]) * gen_means[n-1]);
		gen_stdvs[n] <- (changeness[n-1] * newgen_stdvs[n]) + ((1.0 - changeness[n-1]) * gen_stdvs[n-1]); // to consider: better stdev crossfade
	}

/*
	// Let's have a simple prior on the projections for each stack
	for (whichstack in 1:Nstack){
		bgspec[whichstack] ~ normal(0.0, 100.0);
		hobs[  whichstack] ~ normal(1.0, 0.5);
	}
*/

	// now the obsns actually get generated, via the projections for each stack
	for (n in 1:N) {
		for (whichstack in 1:Nstack) {
/*
			gen_means_stacked[n, whichstack] <- (gen_means[n] .* hobs[whichstack]) + bgspec[whichstack];
			gen_stdvs_stacked[n, whichstack] <- (gen_stdvs[n] .* hobs[whichstack]);
*/
			gen_means_stacked[n, whichstack] <- gen_means[n];
			gen_stdvs_stacked[n, whichstack] <- gen_stdvs[n];

			x[n, whichstack] ~ normal(gen_means_stacked[n, whichstack], gen_stdvs_stacked[n, whichstack]);
		}
	}
}