Is there any way to achieve the effect of target*=real a
?
This will be useful in simulated tempering and some general post-processing, where we want to reshape the posterior by some power transformation.
Would this work
target += target() * a - target()
edit. typo in the docs target()()
--> target()
Thanks. But target()()
does not seem to be recognized.
I think that is a typo in the docs and its just target()
.
Great, thanks! This is exactly what I want.
OK it turns out I have not fully solved my problem. I am now using a geometric bridge between two models: For two densities p_1(\theta) and p_2(\theta), I want to sample from a density that is proportional to p_1^\lambda(\theta) p_2^{(1-\lambda)}(\theta). Using the saved target()
, I can do this by:
parameters{
real lambda;
real theta;
}
model {
real log_q;//lp of the first model;
theta~foo; // model 1
log_q=target();
target+=-target();
real log_p;//lp of the alternative model
theta~foo2; // model 2
log_p=target();
target+=target()*(lambda)-target() + (1-lambda)*log_q;
lambda~foo3
}
The sampling is fine. But for some post process, I want to access my local variable log_q
and log_p
defined in the model block. For one time use, I could put them into transformed parameters, to compute and store log_p=foo1_lpdf(theta)
.
But I want an automated procedure that applied to a general model where the density is hard to compute in the transformed parameters.
In other words, in most cases, local variables in model block can be easily moved back to transformed parameters, expect when these variables depends on target().
In this structure, (before hacking into stan language), is there any easy way I can rewrite my code to
- save local variable in model block ; or
- compute transformed parameters in model block; or
- save intermediate states of target()?
which are three equivalent description my question.
Thanks.
That is going to be less trivial I believe as local variables are not saved in the results and you can’t call the local variables from model {} in generated quantities.
A “trick” you can do is use print. Example:
parameters {
real y;
}
model {
y ~ normal(0,1);
real b = target()*5.0;
print("line = ",y,b);
}
generated quantities {
print("----");
}
The print in GQ is just so we mark iterations. On the standard output you would get:
Iteration: 1001 / 2000 [ 50%] (Sampling)
line = -0.243293-0.147978
line = 1.30349-4.24769
line = 1.33544-4.45849
line = -0.184566-0.0851618
----
line = 1.33544-4.45849
line = 0.316915-0.251088
line = -1.06991-2.86175
line = -1.21335-3.68057
----
line = 0.316915-0.251088
line = 1.0784-2.90739
line = 0.586643-0.860376
line = -0.586876-0.861059
----
...
In order to know which of the printed lines was the one selected by the sampler you can compare y
with the output, which is:
lp__,accept_stat__,stepsize__,treedepth__,n_leapfrog__,divergent__,energy__,y
# Adaptation terminated
# Step size = 1.02103
# Diagonal elements of inverse mass matrix:
# 1.11476
-0.891699,0.855488,1.02103,2,3,0,1.12987,1.33544
-0.0502177,1,1.02103,1,3,0,0.672531,0.316915
-0.172212,0.929122,1.02103,2,3,0,0.470398,-0.586876
-0.324924,0.956602,1.02103,1,1,0,0.351634,-0.806132
-1.33348,0.722085,1.02103,1,3,0,2.28903,1.63308
-0.0789315,0.875911,1.02103,2,3,0,2.27032,0.39732
-0.687362,0.830736,1.02103,1,3,0,1.02702,-1.17249
-1.71255,0.742411,1.02103,1,1,0,1.78042,-1.8507
-0.564982,1,1.02103,1,1,0,1.41475,-1.063
-0.962718,0.890871,1.02103,1,1,0,1.06254,-1.3876
-0.303665,1,1.02103,2,3,0,0.788109,-0.779313
-0.00679864,0.999117,1.02103,2,3,0,0.31023,0.116607
-0.00638563,1,1.02103,1,1,0,0.00844971,0.11301
-0.166766,0.962197,1.02103,1,3,0,0.173918,0.577522
So for the first sample its the third printed line, for the second one its the second line and for the 3rd sample it was the fourth one. With some scripting you could automate this pretty easily.
But not ideal. Still beats hacking in the C++ backend I guess.
Thanks, Rok. It works. It seems the actually saved parameter is always the first printed leapfrog line after each “—”, hence the last iteration will not be saved, but that is fine. That is
--- \\ printed from the generated quantile in the i-th iter"
line = 0.316915-0.251088 \\ first-line in the i+1 iter: corresponds to the saved parameter in the i-th iter
It is also interesting that the printed value is not numerically identical to the saved parameter (they have the same precision), up to one or two digit difference (~1e-6).
By the way, I used the solution here to facilitate a model/alternative model syntax in https://github.com/yao-yl/path-tempering
It could (1) fit two models at the same time (2)ease metastability/multimodality and (3) estimate the normalization constant.
Thanks for your help!
Nice! Glad it worked!
The precision difference is because of different precision settings for stdout (print) and output stream (CSV file).