simul.data.R

N <- 200
K <- 4
timepos <- c(0.0, 2.08978153718, 2.6774000865, 3.47584796473, 3.78710127148, 4.75454597707, 6.66987539583, 7.24902152136, 8.49394193938, 10.4795896694, 10.7914520039, 11.4969222693, 13.3095734861, 14.7212655024, 15.2931492467, 15.4171748097, 15.7443055882, 16.3382608385, 19.567384084, 20.9367692864, 22.0817809742, 22.7254312132, 23.403135471, 23.8066143306, 25.4891533587, 26.6521259878, 28.8649882504, 28.9655706897, 29.1462279535, 30.8224952884, 32.4190207294, 32.4851131989, 33.7009422433, 34.200217965, 35.0046621559, 36.5020537758, 37.2170190379, 37.9581204385, 37.9809986145, 39.206561022, 39.5827104717, 39.9391764411, 39.9627692061, 41.1659966046, 41.8917759456, 44.1068011549, 44.4874876662, 45.468215782, 46.0485677517, 47.0857797338, 47.9105602445, 48.3423708895, 51.2038621099, 51.6573257642, 52.1885429832, 54.4739163581, 55.63263262, 55.7679908666, 57.0184307288, 58.4159798588, 58.6414905971, 61.4434801692, 61.9000345735, 62.64477926, 63.23209667, 65.3153604468, 65.7821474602, 66.3972814977, 67.4819307957, 67.642034874, 68.8888924684, 70.9154847995, 71.0796887623, 72.123688448, 73.8467543593, 75.0013685457, 75.5085864229, 75.5815519081, 77.405511118, 79.4107094328, 81.6429847207, 82.1224833113, 83.2130519322, 84.6526609062, 85.035363047, 85.7708309646, 86.6196880308, 86.802129892, 88.3978378714, 88.8534034396, 90.1716393001, 91.2954501387, 91.7125561618, 91.921002418, 93.2093553598, 93.8779790886, 94.0257278557, 94.3699186507, 96.1814772158, 96.8701307973, 96.9506721126, 97.0301069578, 97.871532448, 98.3398857528, 98.8892141757, 99.4014630423, 101.708311686, 105.909673998, 106.277166547, 106.745664352, 108.140768362, 108.187025859, 108.679062842, 108.886559918, 110.152402437, 111.01097894, 111.127776045, 113.122288143, 113.947880142, 115.331496313, 115.819757142, 116.489159976, 116.853951057, 117.034255824, 117.2726025, 118.216761036, 120.431479605, 121.4546924, 122.488697185, 123.325571302, 125.209467894, 125.84660617, 126.251136429, 128.238910211, 129.072096741, 129.232719802, 130.936217704, 131.463685156, 132.476958317, 133.892908841, 136.486167734, 137.314056686, 137.341791906, 137.604751146, 139.6556208, 140.097871273, 141.476112205, 143.2151502, 143.766581143, 146.191135206, 147.270173402, 147.326830305, 147.44773035, 149.844040599, 150.23441724, 152.550122746, 154.138763897, 157.324771342, 157.473017514, 159.07474021, 162.55387289, 163.017941069, 165.193258031, 166.598314112, 168.098356069, 170.56389857, 171.160160971, 172.322548607, 174.162889455, 175.451145642, 176.868860307, 177.456530806, 178.940508109, 179.882901272, 181.594775482, 183.297201771, 183.696394389, 186.212884067, 188.50415139, 189.267948099, 189.856539696, 190.205324705, 192.56711187, 193.395024254, 193.628192382, 195.052747219, 195.806198177, 196.567631896, 197.63576125, 198.825565729, 200.573751631, 201.567410286, 201.856550009, 203.149915004, 203.718799635, 203.951046746, 205.918864229, 206.804140439, 207.280199522, 208.365817821)
individ <- c(1, 1, 2, 4, 1, 3, 4, 3, 3, 2, 4, 2, 3, 3, 1, 2, 4, 3, 2, 3, 3, 2, 4, 3, 1, 1, 2, 1, 4, 2, 3, 1, 1, 2, 4, 2, 1, 3, 4, 1, 3, 4, 2, 1, 4, 4, 1, 2, 4, 1, 2, 4, 4, 1, 2, 4, 2, 1, 4, 2, 1, 4, 1, 3, 4, 3, 4, 2, 1, 3, 3, 4, 1, 2, 4, 1, 3, 2, 3, 3, 3, 1, 4, 4, 3, 1, 2, 4, 1, 3, 4, 1, 2, 3, 2, 1, 4, 2, 4, 1, 2, 3, 3, 4, 1, 3, 3, 2, 4, 1, 4, 2, 3, 1, 1, 2, 4, 3, 1, 3, 4, 1, 2, 3, 4, 4, 2, 4, 1, 2, 3, 1, 4, 3, 4, 3, 4, 3, 2, 3, 3, 4, 1, 3, 2, 4, 2, 4, 4, 2, 3, 1, 4, 2, 1, 2, 2, 2, 1, 1, 4, 1, 4, 2, 1, 2, 1, 2, 3, 1, 3, 4, 2, 3, 1, 2, 1, 3, 1, 3, 4, 3, 2, 1, 3, 4, 1, 2, 2, 4, 1, 2, 3, 3, 1, 4, 4, 3, 1, 2)

zf4vmrpcens.stan

data {
	int<lower=0>   N;  // num observations
	int<lower=1>   K;  // num individuals
	real   timepos[N]; // timepos of event n
	int<lower=1,upper=K>   individ[N]; // who is responsible for event n
}
parameters {
	real<lower=0> excesstimepos[N,K-1]; // for non-observed param, the excess >0, making it up to the time they "would have" called next. log-gaussian on gap will become truncated due to causality enforcement.

	real gap_self_mean;
	real<lower=0> gap_self_stdv;
}
model {
	vector[K] prevtimes;
	vector[K] currtimes;
	real gap_self;

	// priors
	gap_self_mean ~ normal(-0.7, 1.0);
	gap_self_stdv ~ inv_gamma(1.0, 1.0);

	// initialise the main loop
	prevtimes <- rep_vector(timepos[1] - exp(gap_self_mean), K);

	for(n in 1:N){
		// ok, so timepos[n] is correct for individ k. for the others, they're lower-bounded by it.
		// the moveable bounds are implemented by adding a zero-bounded hidden variable.
		for (k in 1:(individ[n] - 1)){
			currtimes[k] <- timepos[n] + excesstimepos[n,k];
		}
		currtimes[individ[n]] <- timepos[n];
		for (k in (individ[n] + 1):K){
			currtimes[k] <- timepos[n] + excesstimepos[n,k-1];
		}

		// now we place distributions over the time gaps, e.g. independent or not.
		// keep this as a separate loop over k, because might use distribs that depend on all of them.
		for(k in 1:K){
			gap_self <- currtimes[k] - prevtimes[k];
			if(k==individ[n]){
				//gap_self ~ lognormal(gap_self_mean, gap_self_stdv);
				increment_log_prob(lognormal_log(gap_self, gap_self_mean, gap_self_stdv));
			}else{
				increment_log_prob(lognormal_log(gap_self, gap_self_mean, gap_self_stdv));
				increment_log_prob(-lognormal_ccdf_log(timepos[n] - prevtimes[k], gap_self_mean, gap_self_stdv)); // adjust scaling factor of distrib, for censoring truncation
			}
		}//k

		// finally we prepare for the next iteration
		prevtimes[individ[n]] <- timepos[n]; // prevtimes is the previous times actually realised, so only the "winner" updates.
	}//n
}//model