Difficulties using categorical() distribution
reuben.a...@gmail.com
Hi Stan users,
I am using CmdStan 2.14.0 on Windows 7 and also on Ubuntu 14.04.
##########################
Script
##########################
#!/bin/bash
/C/sandbox/mcmcBayes/code/sampler/projectTemplates/Stan/modeltypeid8/modeltypeid8.exe sample num_samples=5000 num_warmup=5000 thin=1 data file=/C/sandbox/mcmcBayes/data/samplerData/modeltypeid8priordatatypeid40datatypeid101datasetid3tempid26samplertypeid3.data.R init=/C/sandbox/mcmcBayes/data/samplerInits/modeltypeid8priordatatypeid40datatypeid101datasetid3tempid26samplertypeid3.init.R output file=/C/sandbox/mcmcBayes/data/samplerCoda/modeltypeid8priordatatypeid40datatypeid101datasetid3tempid26samplertypeid3output.csv
##########################
Error
##########################
$ ../mcmcBayes/temp/modeltypeid8priordatatypeid40datatypeid101datasetid3tempid26samplertypeid3script.txt
method = sample (Default)
sample
num_samples = 5000
num_warmup = 5000
save_warmup = 0 (Default)
thin = 1 (Default)
adapt
engaged = 1 (Default)
gamma = 0.050000000000000003 (Default)
delta = 0.80000000000000004 (Default)
kappa = 0.75 (Default)
t0 = 10 (Default)
init_buffer = 75 (Default)
term_buffer = 50 (Default)
window = 25 (Default)
algorithm = hmc (Default)
hmc
engine = nuts (Default)
nuts
max_depth = 10 (Default)
metric = diag_e (Default)
stepsize = 1 (Default)
stepsize_jitter = 0 (Default)
id = 0 (Default)
data
file = C:/sandbox/mcmcBayes/data/samplerData/modeltypeid8priordatatypeid40datatypeid101datasetid3tempid26samplertypeid3.data.R
init = C:/sandbox/mcmcBayes/data/samplerInits/modeltypeid8priordatatypeid40datatypeid101datasetid3tempid26samplertypeid3.init.R
random
seed = 1708780682
output
file = C:/sandbox/mcmcBayes/data/samplerCoda/modeltypeid8priordatatypeid40datatypeid101datasetid3tempid26samplertypeid3output.csv
diagnostic_file = (Default)
refresh = 100 (Default)
pr_r[1,:]=[0.1,0.1,0.1,0.5,0.1,0.1]
Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
Exception thrown at line 38: categorical_log: Number of categories is -2147483648, but must be in the interval [1, 6]
If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
Rejecting initial value:
Error evaluating the log probability at the initial value.
##########################
Model (Stan)
##########################
functions {
// Define log likelihood function for temperature parametere phi (Note: phi=0 generates the prior's loglikelihood, phi=1 generates the same for the posterior, and the corresponding temperature values in between are generated by the temperature step index, the discretization or k number of steps, and the choice for the temperature parameter c)
real nhpp_loglikelihood_phi_i(int[] N, int[] theta, real phi) {
real retval=0;
int Taux=size(N)-1;
for (i in 1:Taux) {
retval=-log_rising_factorial(1,N[i+1]-N[i])+(N[i+1]-N[i])*log(theta[i+1]-theta[i])-(theta[i+1]-theta[i]);
}
return phi*retval;
}
}
data {
real<lower=0, upper=1> phi;
int<lower=0> Tau;
int<lower=0> Rho;
real<lower=0> t[Tau];
int<lower=0> N[Tau];
vector[Rho] pr_r[Tau-1];//for each interval's PMF, contains an array of vectors with point estimates
int Nindices[Rho];//array of possible values for N(t)=n, that is 0:Rho-1
//int<lower=0> Lambdap[Tau];//MATLAB uses this to generate pr_r and Nindices above
//real<lower=0> c;//MATLAB uses this to generate pr_r and Nindices above
}
parameters {
}
model {
int r_[Tau-1];
int theta[Tau];
theta[1]=0;
for (i in 1:Tau-1) {
print("pr_r[", i, ",:]=", pr_r[i]);
r_[i]~categorical(pr_r[i]);//Using the categorical distribution to generate values from the gamma prior's PMF
print("r_[", i, "]=", r_[i]);
theta[i+1]=r_[i]-1+theta[i];//Subtracting one as my values are from 0:Rho-1
}
print("theta=", theta);
if (phi>0)//we don't need the loglikelihood for the prior step
target+=nhpp_loglikelihood_phi_i(N, theta, phi);//Model assumes a NHPP with a Gamma process cumulative intensity function for likelihood of the data
}
##########################
Inits:
##########################
r_<-c(5, 5, 5, 5, 5)
##########################
Data:
##########################
phi<-1
Tau<-6
t<-c(0, 10, 20, 30, 40, 50)
N<-c(0, 4, 8, 12, 16, 20)
Lambdap<-c(0, 2, 6, 11, 12, 13)
c<-0.468
Nindices<-c(0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99, 100, 101, 102, 103, 104, 105, 106, 107, 108, 109, 110, 111, 112, 113, 114, 115, 116, 117, 118, 119, 120, 121, 122, 123, 124, 125, 126, 127, 128, 129, 130, 131, 132, 133, 134, 135, 136, 137, 138, 139, 140, 141, 142, 143, 144, 145, 146, 147, 148, 149, 150, 151, 152, 153, 154, 155, 156, 157, 158, 159, 160, 161, 162, 163, 164, 165, 166, 167, 168, 169, 170, 171, 172, 173, 174, 175, 176, 177, 178, 179, 180, 181, 182, 183, 184, 185, 186, 187, 188, 189, 190, 191, 192, 193, 194, 195, 196, 197, 198, 199, 200)
Rho<-201
pr_r<-structure(
c(
0.0579043, 0.0129129, 0.00809342, 0.220137, 0.220137,
0.280925, 0.0637826, 0.0368704, 0.57605, 0.57605,
0.366089, 0.143729, 0.087054, 0.187957, 0.187957,
0.211525, 0.203888, 0.140096, 0.0153937, 0.0153937,
0.0677782, 0.208021, 0.171838, 0.000456263, 0.000456263,
0.0136839, 0.16428, 0.170768, 6.09983E-6, 6.09983E-6,
0.00189044, 0.105212, 0.142905, 4.2501E-8, 4.2501E-8,
0.000189328, 0.0564352, 0.103418, 1.70766E-10, 1.70766E-10,
1.43436E-5, 0.0259601, 0.0659947, 4.26257E-13, 4.26257E-13,
8.49382E-7, 0.0104275, 0.0376906, 6.9976E-16, 6.9976E-16,
4.0345E-8, 0.00370998, 0.019492, 7.90214E-19, 7.90214E-19,
1.56979E-9, 0.00118275, 0.00921505, 6.36536E-22, 6.36536E-22,
5.09095E-11, 0.00034111, 0.00401375, 3.7686E-25, 3.7686E-25,
1.39638E-12, 8.97125E-5, 0.00162131, 1.68154E-28, 1.68154E-28,
3.27993E-14, 2.16639E-5, 0.000610769, 5.77654E-32, 5.77654E-32,
6.66864E-16, 4.8317E-6, 0.000215612, 1.55606E-35, 1.55606E-35,
1.18463E-17, 1.0004E-6, 7.16269E-5, 3.33973E-39, 3.33973E-39,
1.85384E-19, 1.93157E-7, 2.24744E-5, 5.79159E-43, 5.79159E-43,
2.57432E-21, 3.4918E-8, 6.68234E-6, 8.21599E-47, 8.21599E-47,
3.19282E-23, 5.93108E-9, 1.88825E-6, 9.64014E-51, 9.64014E-51,
3.55741E-25, 9.49623E-10, 5.08413E-7, 9.44837E-55, 9.44837E-55,
3.57951E-27, 1.43732E-10, 1.30744E-7, 7.80451E-59, 7.80451E-59,
3.26823E-29, 2.06196E-11, 3.21816E-8, 5.47711E-63, 5.47711E-63,
2.71948E-31, 2.81035E-12, 7.59653E-9, 3.28973E-67, 3.28973E-67,
2.07047E-33, 3.64709E-13, 1.72274E-9, 1.7025E-71, 1.7025E-71,
1.44757E-35, 4.51548E-14, 3.75954E-10, 7.63841E-76, 7.63841E-76,
9.32531E-38, 5.34364E-15, 7.90698E-11, 2.98792E-80, 2.98792E-80,
5.55243E-40, 6.05459E-16, 1.60493E-11, 1.02437E-84, 1.02437E-84,
3.06443E-42, 6.57861E-17, 3.14796E-12, 3.09295E-89, 3.09295E-89,
1.57191E-44, 6.86476E-18, 5.97386E-13, 8.26182E-94, 8.26182E-94,
7.51267E-47, 6.88894E-19, 1.09804E-13, 1.96059E-98, 1.96059E-98,
3.35325E-49, 6.65691E-20, 1.95692E-14, 4.14962E-103, 4.14962E-103,
1.40084E-51, 6.20164E-21, 3.38492E-15, 7.86201E-108, 7.86201E-108,
5.48843E-54, 5.57627E-22, 5.68773E-16, 1.338E-112, 1.338E-112,
2.02061E-56, 4.84444E-23, 9.29232E-17, 2.05202E-117, 2.05202E-117,
7.00294E-59, 4.07042E-24, 1.47726E-17, 2.84467E-122, 2.84467E-122,
2.28866E-61, 3.31085E-25, 2.28701E-18, 3.57481E-127, 3.57481E-127,
7.06459E-64, 2.60934E-26, 3.45047E-19, 4.08341E-132, 4.08341E-132,
2.06283E-66, 1.99424E-27, 5.07672E-20, 4.25068E-137, 4.25068E-137,
5.70607E-69, 1.47921E-28, 7.28897E-21, 4.04216E-142, 4.04216E-142,
1.4973E-71, 1.06565E-29, 1.02187E-21, 3.51959E-147, 3.51959E-147,
3.73201E-74, 7.4618E-31, 1.39969E-22, 2.81221E-152, 2.81221E-152,
8.84674E-77, 5.08179E-32, 1.87418E-23, 2.06626E-157, 2.06626E-157,
1.99684E-79, 3.36835E-33, 2.45454E-24, 1.39884E-162, 1.39884E-162,
4.29648E-82, 2.17428E-34, 3.14579E-25, 8.74219E-168, 8.74219E-168,
8.82192E-85, 1.36763E-35, 3.94728E-26, 5.05276E-173, 5.05276E-173,
1.73038E-87, 8.38742E-37, 4.85155E-27, 2.70547E-178, 2.70547E-178,
3.24547E-90, 5.01796E-38, 5.84343E-28, 1.34426E-183, 1.34426E-183,
5.82613E-93, 2.93017E-39, 6.89996E-29, 6.20781E-189, 6.20781E-189,
1.00195E-95, 1.67087E-40, 7.99087E-30, 2.6685E-194, 2.6685E-194,
1.65216E-98, 9.30865E-42, 9.07989E-31, 1.06931E-199, 1.06931E-199,
2.61432E-101, 5.069E-43, 1.01267E-31, 3.99994E-205, 3.99994E-205,
3.973E-104, 2.69925E-44, 1.10896E-32, 1.39863E-210, 1.39863E-210,
5.80316E-107, 1.40616E-45, 1.19283E-33, 4.57733E-216, 4.57733E-216,
8.15307E-110, 7.16929E-47, 1.26065E-34, 1.40386E-221, 1.40386E-221,
1.10255E-112, 3.5788E-48, 1.3095E-35, 4.03976E-227, 4.03976E-227,
1.43614E-115, 1.74978E-49, 1.33737E-36, 1.09197E-232, 1.09197E-232,
1.80303E-118, 8.38257E-51, 1.34325E-37, 2.7757E-238, 2.7757E-238,
2.18323E-121, 3.93613E-52, 1.32725E-38, 6.64216E-244, 6.64216E-244,
2.55126E-124, 1.81221E-53, 1.2905E-39, 1.49786E-249, 1.49786E-249,
2.8789E-127, 8.18356E-55, 1.23507E-40, 3.18634E-255, 3.18634E-255,
3.13885E-130, 3.62582E-56, 1.16376E-41, 6.40024E-261, 6.40024E-261,
3.30847E-133, 1.57665E-57, 1.07991E-42, 1.21503E-266, 1.21503E-266,
3.37314E-136, 6.73071E-59, 9.87115E-44, 2.18205E-272, 2.18205E-272,
3.32827E-139, 2.82168E-60, 8.89015E-45, 3.71027E-278, 3.71027E-278,
3.17982E-142, 1.16198E-61, 7.8906E-46, 5.97835E-284, 5.97835E-284,
2.94307E-145, 4.70169E-63, 6.90347E-47, 9.13589E-290, 9.13589E-290,
2.64009E-148, 1.86976E-64, 5.95492E-48, 1.32514E-295, 1.32514E-295,
2.29647E-151, 7.30981E-66, 5.06554E-49, 1.82579E-301, 1.82579E-301,
1.93785E-154, 2.81012E-67, 4.25017E-50, 2.39137E-307, 2.39137E-307,
1.58705E-157, 1.06254E-68, 3.51806E-51, 2.97966E-313, 2.97966E-313,
1.26199E-160, 3.95248E-70, 2.87342E-52, 3.53445E-319, 3.53445E-319,
9.74748E-164, 1.44677E-71, 2.31619E-53, 0, 0,
7.31608E-167, 5.21228E-73, 1.84294E-54, 0, 0,
5.33805E-170, 1.84864E-74, 1.44771E-55, 0, 0,
3.78766E-173, 6.45595E-76, 1.12296E-56, 0, 0,
2.61458E-176, 2.22047E-77, 8.60266E-58, 0, 0,
1.75645E-179, 7.52296E-79, 6.50959E-59, 0, 0,
1.14874E-182, 2.51119E-80, 4.86629E-60, 0, 0,
7.31666E-186, 8.26034E-82, 3.59446E-61, 0, 0,
4.53994E-189, 2.67809E-83, 2.62376E-62, 0, 0,
2.74521E-192, 8.55929E-85, 1.89293E-63, 0, 0,
1.61819E-195, 2.6972E-86, 1.34998E-64, 0, 0,
9.30131E-199, 8.38157E-88, 9.51843E-66, 0, 0,
5.21496E-202, 2.5689E-89, 6.63597E-67, 0, 0,
2.85285E-205, 7.76695E-91, 4.57513E-68, 0, 0,
1.52319E-208, 2.31688E-92, 3.11975E-69, 0, 0,
7.93957E-212, 6.81982E-94, 2.10431E-70, 0, 0,
4.04137E-215, 1.98119E-95, 1.40418E-71, 0, 0,
2.00939E-218, 5.68104E-97, 9.27083E-73, 0, 0,
9.76156E-222, 1.6082E-98, 6.0568E-74, 0, 0,
4.63452E-225, 4.49497E-100, 3.91606E-75, 0, 0,
2.15095E-228, 1.24064E-101, 2.50603E-76, 0, 0,
9.76118E-232, 3.38189E-103, 1.58746E-77, 0, 0,
4.33237E-235, 9.10589E-105, 9.95508E-79, 0, 0,
1.88106E-238, 2.4221E-106, 6.18101E-80, 0, 0,
7.9916E-242, 6.36539E-108, 3.80007E-81, 0, 0,
3.3229E-245, 1.65301E-109, 2.31359E-82, 0, 0,
1.35253E-248, 4.24223E-111, 1.39504E-83, 0, 0,
5.3904E-252, 1.07606E-112, 8.33173E-85, 0, 0,
2.10392E-255, 2.69806E-114, 4.92917E-86, 0, 0,
8.04379E-259, 6.68791E-116, 2.88897E-87, 0, 0,
3.01305E-262, 1.63908E-117, 1.67758E-88, 0, 0,
1.10599E-265, 3.97221E-119, 9.65236E-90, 0, 0,
3.97909E-269, 9.51987E-121, 5.50342E-91, 0, 0,
1.40341E-272, 2.25654E-122, 3.1097E-92, 0, 0,
4.85329E-276, 5.29073E-124, 1.74152E-93, 0, 0,
1.64596E-279, 1.22713E-125, 9.66713E-95, 0, 0,
5.47536E-283, 2.81588E-127, 5.3194E-96, 0, 0,
1.78688E-286, 6.39333E-129, 2.90174E-97, 0, 0,
5.7219E-290, 1.43639E-130, 1.56935E-98, 0, 0,
1.79814E-293, 3.19369E-132, 8.41552E-100, 0, 0,
5.54652E-297, 7.02792E-134, 4.4748E-101, 0, 0,
1.67958E-300, 1.53079E-135, 2.35956E-102, 0, 0,
4.99385E-304, 3.30063E-137, 1.23392E-103, 0, 0,
1.45813E-307, 7.04548E-139, 6.39983E-105, 0, 0,
4.18164E-311, 1.489E-140, 3.29238E-106, 0, 0,
1.17804E-314, 3.11591E-142, 1.68013E-107, 0, 0,
3.2606E-318, 6.45682E-144, 8.50536E-109, 0, 0,
8.84378E-322, 1.32505E-145, 4.27162E-110, 0, 0,
0, 2.69315E-147, 2.12848E-111, 0, 0,
0, 5.42174E-149, 1.05234E-112, 0, 0,
0, 1.08119E-150, 5.16265E-114, 0, 0,
0, 2.13589E-152, 2.51334E-115, 0, 0,
0, 4.18031E-154, 1.21428E-116, 0, 0,
0, 8.10627E-156, 5.82242E-118, 0, 0,
0, 1.55757E-157, 2.77094E-119, 0, 0,
0, 2.96566E-159, 1.30894E-120, 0, 0,
0, 5.59594E-161, 6.13767E-122, 0, 0,
0, 1.04648E-162, 2.85697E-123, 0, 0,
0, 1.93968E-164, 1.32023E-124, 0, 0,
0, 3.56364E-166, 6.05702E-126, 0, 0,
0, 6.49017E-168, 2.75904E-127, 0, 0,
0, 1.17177E-169, 1.24788E-128, 0, 0,
0, 2.09742E-171, 5.6043E-130, 0, 0,
0, 3.72227E-173, 2.49937E-131, 0, 0,
0, 6.54997E-175, 1.10693E-132, 0, 0,
0, 1.1429E-176, 4.86875E-134, 0, 0,
0, 1.9776E-178, 2.12685E-135, 0, 0,
0, 3.39358E-180, 9.22796E-137, 0, 0,
0, 5.77556E-182, 3.97688E-138, 0, 0,
0, 9.74921E-184, 1.70242E-139, 0, 0,
0, 1.63234E-185, 7.2394E-141, 0, 0,
0, 2.71108E-187, 3.05821E-142, 0, 0,
0, 4.46673E-189, 1.28346E-143, 0, 0,
0, 7.30088E-191, 5.35141E-145, 0, 0,
0, 1.18392E-192, 2.21688E-146, 0, 0,
0, 1.90484E-194, 9.12485E-148, 0, 0,
0, 3.04091E-196, 3.73195E-149, 0, 0,
0, 4.81704E-198, 1.51667E-150, 0, 0,
0, 7.57203E-200, 6.12509E-152, 0, 0,
0, 1.18119E-201, 2.45819E-153, 0, 0,
0, 1.82864E-203, 9.80436E-155, 0, 0,
0, 2.80969E-205, 3.88635E-156, 0, 0,
0, 4.28478E-207, 1.53108E-157, 0, 0,
0, 6.48578E-209, 5.99528E-159, 0, 0,
0, 9.74496E-211, 2.33341E-160, 0, 0,
0, 1.45346E-212, 9.02729E-162, 0, 0,
0, 2.15204E-214, 3.47159E-163, 0, 0,
0, 3.16331E-216, 1.32715E-164, 0, 0,
0, 4.61634E-218, 5.04365E-166, 0, 0,
0, 6.68866E-220, 1.90556E-167, 0, 0,
0, 9.6224E-222, 7.15763E-169, 0, 0,
0, 1.37452E-223, 2.673E-170, 0, 0,
0, 1.94966E-225, 9.92491E-172, 0, 0,
0, 2.74616E-227, 3.66411E-173, 0, 0,
0, 3.84126E-229, 1.34505E-174, 0, 0,
0, 5.33601E-231, 4.90968E-176, 0, 0,
0, 7.36163E-233, 1.78207E-177, 0, 0,
0, 1.0087E-234, 6.43237E-179, 0, 0,
0, 1.37279E-236, 2.30889E-180, 0, 0,
0, 1.85571E-238, 8.2421E-182, 0, 0,
0, 2.49172E-240, 2.92609E-183, 0, 0,
0, 3.32346E-242, 1.03316E-184, 0, 0,
0, 4.40352E-244, 3.62816E-186, 0, 0,
0, 5.7962E-246, 1.26726E-187, 0, 0,
0, 7.57944E-248, 4.40262E-189, 0, 0,
0, 9.84687E-250, 1.52139E-190, 0, 0,
0, 1.27099E-251, 5.22959E-192, 0, 0,
0, 1.62998E-253, 1.78814E-193, 0, 0,
0, 2.07701E-255, 6.08218E-195, 0, 0,
0, 2.62981E-257, 2.05802E-196, 0, 0,
0, 3.30868E-259, 6.92767E-198, 0, 0,
0, 4.13661E-261, 2.31998E-199, 0, 0,
0, 5.13936E-263, 7.72949E-201, 0, 0,
0, 6.34547E-265, 2.56213E-202, 0, 0,
0, 7.78615E-267, 8.44983E-204, 0, 0,
0, 9.49511E-269, 2.77269E-205, 0, 0,
0, 1.15083E-270, 9.0526E-207, 0, 0,
0, 1.38633E-272, 2.94087E-208, 0, 0,
0, 1.65992E-274, 9.50649E-210, 0, 0,
0, 1.97551E-276, 3.05786E-211, 0, 0,
0, 2.33702E-278, 9.78767E-213, 0, 0,
0, 2.74819E-280, 3.11757E-214, 0, 0,
0, 3.21253E-282, 9.88189E-216, 0, 0,
0, 3.73315E-284, 3.11717E-217, 0, 0,
0, 4.31268E-286, 9.78566E-219, 0, 0,
0, 4.95306E-288, 3.0573E-220, 0, 0,
0, 5.65546E-290, 9.50636E-222, 0, 0,
0, 6.42012E-292, 2.94192E-223, 0, 0,
0, 7.24623E-294, 9.06144E-225, 0, 0),
.Dim=c(5,201))
##########################
Model (JAGS)
##########################
#Description: Doubly stochastic NHPP, with a Gamma process for the cumulative intensity function (Kuo-Ghosh model).
#Note, Theta(t) must be constructed from the individual r intervals. Otherwise, the resulting random variables for
#Theta(t) fail to meet the non-decreasing requirement.
#unlike OpenBUGS, JAGS forbids stochastic nodes to be declared twice
data {
##INDICES,START##
#JAGS appears to not support expressions for indices in the model, create a variable for this index maximum
T1 <- length(N)
T2 <- dim(pr_r)
##INDICES,STOP##
##LIKELIHOODDATA,START##
for (i in 2:length(N)) {
#zero the dummy variables for the custom sampling distribution
dummylogN[i-1] <- 0
}
##LIKELIHOODDATA,STOP##
}
model {
#terms for simplifying the joint loglikelihood below (these are not dependent on independent variables and are located outside the area preprocessed)
C <- 1.0e+5 #this just has to be large enough to ensure all dummylog[i]'s >0
for (i in 2:T1) {
#terms for simplifying the joint loglikelihood below
A[i-1] <- Theta[i]-Theta[i-1]
}
##LIKELIHOOD,START##
#joint likelihood function
for (i in 2:T1) {
#terms for simplifying the joint loglikelihood below
BN[i-1] <- N[i]-N[i-1]
#Using the zeros trick for custom loglikelihood
dummylogN[i-1] ~ dpois(logN[i-1])
logN[i-1] <- -loglikelihoodN[i-1] + C
loglikelihoodN[i-1] <- phi*(-1*logfact(BN[i-1])+BN[i-1]*log(A[i-1])-A[i-1])
}
##LIKELIHOOD,STOP##
##VARIABLE,1##
#prior distribution for Theta (gamma process)
#translation of the stochastic node intervals to cumulative form
Theta[1] <- 0
for (i in 2:T1) {
Theta[i] <- Theta[i-1]+r[i-1]
#Using the dcat distribution to generate values from the gamma prior's PMF
#pr_r[interval,PMF] matrix contains an array of point estimates for each interval's PMF
#Nindices[] array of possible values for N(t)=n
r_[i-1] ~ dcat(pr_r[i-1,])#sample the dcat rv to get the index
r[i-1] <- Nindices[r_[i-1]]#transform the index to the interval count
}
}
#FN1: "A Poisson(phi) observation of zero has likelihood exp(-phi), so if our observed data is a set of 0's, and
#phi[i] is set to - log(L[i]), we will obtain the correct likelihood contribution. (Note that phi[i] should always
#be > 0 as it is a Poisson mean, and so we may need to add a suitable constant to ensure that it is positive.)"
#Source: "Tricks: Advanced Use of the BUGS Language", http://users.aims.ac.za/~mackay/BUGS/Manuals/Tricks.html
Bob Carpenter
y ~ foo(theta);
It just increments the log density.
You have to marginalize out discrete parameters in Stan.
There are chapters in the manual for latent discrete parameters
and mixture models that have examples (and also HMMs in the
time-series chapter).
I didn't see where the variable r_ was being used in our
Stan program, so it's hard to provide concrete advice.
- Bob
> You received this message because you are subscribed to the Google Groups "Stan users mailing list" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.
> To post to this group, send email to stan-...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.
reuben.a...@gmail.com
Bob Carpenter
in generated quantities using (categorical)_rng functions.
Stan doesn't allow discrete parameters. If you have categorical
data, that's fine, you just write it like everything else. If
you want categorical parameters, you need to marginalize them out.
If you need categorical predictions, put them in the generated quantities
block.
- Bob
Reuben Johnston
> To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+unsubscribe@googlegroups.com.
> To post to this group, send email to stan-...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.
--
You received this message because you are subscribed to a topic in the Google Groups "Stan users mailing list" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/stan-users/8HrULFanPok/unsubscribe.
To unsubscribe from this group and all its topics, send an email to stan-users+unsubscribe@googlegroups.com.
Reuben Johnston
Bob Carpenter
to write your own discrete sampler. Or if you think you do, you'll
find it won't work in inference because the step functions involved
break derivatives.
You will be able to marginalize out, though.
You might write a simple normal mixture model in JAGS model like this:
z[n] ~ dbern(theta);
y[n] ~ dnorm(mu[z[n]], sigma[z[n]]);
In Stan, you need to marginalize out that z parameter. So you
write the same thing with no discrete parameters as:
target += log_sum_exp(log(theta) + normal_lpdf(y[n] | mu[1], sigma[1]),
log1m(theta) + normal_lpdf(y[n] | mu[2], sigma[2]));
Works the same way with more than two values for z.
It gives you the exact same inference on the parameters mu and sigma,
except it's a lot faster and a lot more robust (as in it can be more
than two orders of magnitude faster). No surprise---just a consequence
of the Rao-Blackwell theorem and the fact that discrete sampling tends to
be intractable in practice---it just mixes terribly if there are any
correlations and has a really hard time with tails.
You can reconstruct anything you need to do with z.
You can see this implicitly in the mixture formulation above, where
mu[1], sigma[1] and mu[2], sigma[2] are both used.
You can reconstruct everything you need in expectation. Specifically,
you can calculate Pr[z[n] == 1 | y]
Check out the manual chapter on latent discrete parameters to see how all
that works.
Or you could do the same marginalization in JAGS using the zeros trick.
Others have found this hugely sped up their state-space models with latent
discrete parameters.
- Bob
> To unsubscribe from this topic, visit https://groups.google.com/d/topic/stan-users/8HrULFanPok/unsubscribe.
Reuben Johnston
The parameter for this NHPP's model is a discrete stochastic process, the gamma process, that for this case is also the mean-value function for the NHPP. I don't see how the marginalizing approach would work for this application.
> > To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+unsubscribe@googlegroups.com.
> > To post to this group, send email to stan-...@googlegroups.com.
> > For more options, visit https://groups.google.com/d/optout.
>
> --
> You received this message because you are subscribed to a topic in the Google Groups "Stan users mailing list" group.
> To unsubscribe from this topic, visit https://groups.google.com/d/topic/stan-users/8HrULFanPok/unsubscribe.
> To unsubscribe from this group and all its topics, send an email to stan-users+unsubscribe@googlegroups.com.
> To post to this group, send email to stan-...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.
>
>
>
> --
> You received this message because you are subscribed to the Google Groups "Stan users mailing list" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+unsubscribe@googlegroups.com.
> To post to this group, send email to stan-...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.
--
You received this message because you are subscribed to a topic in the Google Groups "Stan users mailing list" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/stan-users/8HrULFanPok/unsubscribe.
To unsubscribe from this group and all its topics, send an email to stan-users+unsubscribe@googlegroups.com.
Bob Carpenter
> On Feb 22, 2017, at 7:22 PM, Reuben Johnston <reuben.a...@gmail.com> wrote:
>
> The parameter for this NHPP's model is a discrete stochastic process, the gamma process, that for this case is also the mean-value function for the NHPP. I don't see how the marginalizing approach would work for this application.
state like a Gaussian process or are the outputs conditionally
independent based on the discrete states like an HMM?
If it's the latter, it'd work like the forward algorithm in HMMs.
If it's the former, you're probably out of luck because the
margninalization will be intractable.
I had a hard time seeing the overall structure of the process from
your code. But then I didn't understand the Wikipedia page on
Gamma processes either.
> I do believe that it would be trivial to code a user-defined function for an inverse transform sampler that generates discrete values from the gamma process. All that would need to be called from this custom function would be functions to query from a gamma variable's cdf and to generate random uniform numbers. It would also necessitate performing the call to this custom function outside the bounds of the time steps in the main Stan model's sampling progression. However, I have no idea on the level of work involved with making this capability available in Stan.
Stan's model block---only after the fact in generated quantities. The
reason we do this is because they don't play nicely with inference.
They basically perform something like the "cut" in BUGS implicitly.
We lose the information connection statistically and we lose the
derivative information computationally.
- Bob
Reuben Johnston
I am using the counting process, non-stationary form of the Gamma process and I think of it as being very similar to the non-homogeneous Poisson process, but the interarrival times are gamma distributed, rather than exponential (as in the NHPP). A more formal definition is provided by van Noortwijk:
"Let v(t) be a non-decreasing, right-continuous,real-valued function for t>=0, with v(0)=0. The gamma process with shape function v(t)>0 and scale parameter u>0 is a continuous-time stochastic process {X(t),t>=0} with the following properties:
(1) X(0)=0 with probability one;
(2) X(tau)-X(t)~Ga(v(tau)-v(t),u), for all tau>t>=0;
(3) X(t) has independent increments
Pages 8-10 of 2007, Merrick and Soyer provide the best description of how I am using the gamma process. It will also be described in my dissertation, but that is not finished yet. A very concise description is as follows.
I estimated the values for the Gamma process's probability mass function and passed in the results for a subset of the realizations (e.g., pr(Θ(t)=n, for n=1:250) as an array to BUGS/JAGS. I then used the dcat() distribution with this pmf, to simulate values from the Gamma process prior (i.e., dcat provides samples from Θ(t) in the range 1<=N<=250, according to the probabilities in the pmf passed to dcat), in the sampling steps in BUGS/JAGS. In this model, the parameter Θ(t) is used as the cumulative intensity (or mean-value) function for a non-homogeneous Poisson process. For this alternative, the prior for Θ(t)~GammaProcess(c*Λ(t),c), where the hyper-parameter Λ(t) is described as a best guess for the NHPP mean value function and c is a constant describing the measure of certainty for the best guess given the prior history, where lower values indicate higher confidence in the best guess. I am using the counting form of the NHPP, so my data are grouped-interval counts over time. An example dataset for the intervals t=(0-10], (10-20], (20-30], would be the respective cumulative counts 10, 25, and 50.
References:
* 2008, Merrick, Soyer, Computational Issues in Semiparametric Bayesian Replacement Models, Technical Report TR-2008-4
* 2007, Merrick, Soyer, Semiparametric Bayesian Decision Models for Optimal Replacement, Technical Report TR-2007-4
* 2007, van Noortwijk, A survey of the application of gamma processes in maintenance
Bob Carpenter
answer it. Would you mind reposting the question with
"gamma process" in the title and hopefully some of
our devs or users who know more about these processes than
me can help out.
If there are unknown count variables, they usually can't
be marginalized out. If there are just categorical
variables that are data-item related, they usually can
be marginalized out.
- Bob
> For more options, visit https://groups.google.com/d/optout.
>
>
> You received this message because you are subscribed to the Google Groups "Stan users mailing list" group.