acceptance rates vs. target

[Ben Goodrich]

If you run the Stan examples with the default target acceptance rate of 0.8, the actual acceptance rates are usually a lot higher. Is there some setting that makes it such that we prefer to have a higher-than-targeted acceptance rate over a lower-than-targeted acceptance rate? Here are the acceptance rates:

> mat
                                chain:1 chain:2 chain:3 chain:4
binormal                          0.964   0.957   0.965   0.950
inv_wishart                       0.981   0.015   0.976   0.093
normal_mixture                    0.735   0.837   0.777   0.810
normal                            0.980   0.960   0.892   0.949
ordered_pair                      0.923   0.942   0.904   0.970
sorted_pair                       0.947   0.950   0.954   0.948
triangle                          0.929   0.916   0.911   0.888
uniform                           0.898   0.951   0.960   0.934
wishart2                          0.989   0.993   0.993   0.994
wishart2x2                        0.978   0.970   0.941   0.959
wishart                           0.950   0.963   0.982   0.962
bernoulli                         0.899   0.964   0.940   0.949
negative_binomial2                1.000   0.991   0.982   0.990
negative_binomial                 0.977   0.967   0.945   0.966
normal_censored                   0.931   0.931   0.954   0.940
normal_loc                        0.989   0.948   0.929   0.966
normal_mixture_k_prop             0.996   0.993   0.997   0.992
normal_mixture_k                  0.980   0.963   0.978   0.982
normal_mixture                    0.976   0.974   0.975   0.951
normal_truncated                  0.927   0.961   0.951   0.940
blocker                           0.911   0.822   0.996   0.970
bones                             0.974   0.949   0.964   0.974
dogs                              0.945   0.947   0.965   0.989
dyes                              0.956   0.843   0.961   0.923
epil                              0.999   1.000   1.000   1.000
equiv                             0.994   0.965   0.972   0.962
inhalers                          0.998   0.998   0.994   0.995
kidney                            0.999   0.995   1.000   1.000
leukfr                            0.973   0.997   0.865   0.989
leuk                              0.981   0.988   0.990   0.978
lsat                              0.995   0.993   0.998   0.996
magnesium                         0.961   0.787   0.528   0.978
mice                              1.000   1.000   0.999   1.000
oxford                            0.983   0.983   0.985   0.995
pump                              0.995   0.975   0.978   0.990
rats                              0.980   0.976   0.990   0.981
rats_vec                          0.988   0.996   0.988   0.990
rats_vec_unit                     0.997   0.989   0.995   0.999
salm2                             1.000   1.000   0.001   1.000
salm                              0.001   0.001   0.001   0.001
seeds                             0.998   0.996   0.993   0.991
stacks_a_normal                   0.974   0.972   0.985   0.995
stacks_b_dexp                     0.982   0.986   0.985   0.965
stacks_c_t4                       0.972   0.951   0.974   0.983
stacks_d_normal_ridge             0.961   0.984   0.980   0.977
stacks_e_dexp_ridge               0.991   0.969   0.985   0.980
stacks_f_t4_ridge                 0.978   0.982   0.980   0.982
surgical_simple                   0.989   0.978   0.990   0.947
surgical                          0.909   0.963   0.983   0.980
air                               0.993   0.985   0.988   0.985
alli2                             0.998   0.998   0.995   0.998
beetles_cloglog                   0.977   0.972   0.979   0.936
beetles_logit                     0.949   0.960   0.957   0.940
beetles_probit                    0.940   0.972   0.965   0.968
birats                            0.987   0.999   0.984   1.000
cervix2                           0.997   0.990   0.994   0.989
dugongs                           0.977   0.984   0.996   0.991
endo1                             0.949   0.955   0.952   0.962
endo2                             0.937   0.908   0.927   0.961
endo3                             0.988   0.987   0.990   0.991
eyes                              0.970   0.976   0.987   0.977
hearts                            0.939   0.962   0.974   0.965
ice                               0.902   0.941   0.979   0.952
jaws                              0.988   0.995   0.997   0.996
mvn_orange                        1.000   1.000   0.993   0.994
orange                            0.969   0.977   0.997   0.973
schools                           1.000   0.001   1.000   1.000
stagnant2                         0.993   0.987   0.994   0.985
stagnant                          0.990   0.986   0.957   0.983
estdof                            0.947   0.952   0.960   0.969
simt4                             0.996   0.991   0.999   0.993
camel2                            0.954   0.869   0.873   0.954
camel                             0.987   0.976   0.973   0.995
seeds                             0.992   0.992   0.988   0.988
circle                            0.941   0.951   0.956   0.922
hsquare                           0.931   0.918   0.911   0.935
parallelogram                     0.919   0.923   0.915   0.902
ring                              0.946   0.967   0.962   0.880
squaremc                          0.966   0.962   0.984   0.953
hepatitisME                       0.980   1.000   1.000   0.960
hepatitis                         0.832   0.913   1.000   0.761
command                           0.980   0.960   0.892   0.949
corr-lda                          1.000   1.000   0.997   0.952
lda                               0.975   0.987   0.988   0.975
naive-bayes                       0.988   0.984   0.991   0.990
naive-bayes-unsup                 0.981   0.988   0.990   0.981
soft-k-means                      0.997   0.970   0.969   0.985
fx_equicorr                          NA      NA      NA      NA
fx_factor                         0.973   0.941   0.980   0.904
eight_schools                     0.765   0.965   0.967   0.898
funnel_reparam                    0.970   0.977   0.966   0.984
funnel                            0.881   0.890   0.757   0.588
arch1                                NA      NA      NA      NA
garch1_1                             NA      NA      NA      NA
koyck                             0.992   0.985   0.981   0.977
gp-fit                            0.966   0.980   0.953   0.958
gp-logit-predict                  0.995   0.994   0.997   0.996
gp-multi-sim                         NA      NA      NA      NA
gp-predict-analytic               0.990   0.998   0.997   0.990
gp-predict                        0.998   0.994   0.998   1.000
gp-sim-cholesky                   0.990   0.991   0.988   0.990
gp-sim                            0.986   0.986   0.984   0.988
hmm-analytic                      0.992   0.993   0.980   0.996
hmm-semisup                       0.997   0.996   0.994   0.994
hmm                               0.998   0.995   0.999   0.988
hmm-sufficient                    0.990   0.993   0.992   0.993
irt2_multilevel                   1.000   0.999   0.999   0.998
irt2                              0.997   0.999   1.000   0.999
irt_multilevel                    1.000   1.000   1.000   1.000
irt                               1.000   0.998   0.999   1.000
regression                           NA      NA      NA      NA
regression_std                       NA      NA      NA      NA
arma11_alt                        0.001   0.979   0.974   0.001
arma11                            1.000   0.986   0.983   0.001
ma2                               1.000   0.981   0.963   1.000
maQ                                  NA      NA      NA      NA
stochastic-volatility-optimized   0.982   0.998   0.993   0.974
stochastic-volatility             0.998   0.999   1.000   0.999
nnmf                              0.996   0.995   0.998   0.992
nnmf_vec                          0.999   0.997   0.996   0.991
logistic                          0.962   0.954   0.984   0.958
bounded_double                    0.948   0.894   0.925   0.959
correlation_matrix                0.928   0.971   0.951   0.976
simplex                           0.925   0.933   0.927   0.944
unit_vector                       0.969   0.962   0.955   0.946

Ben

[Bob Carpenter]

This is another piece of the puzzle we need to put into the
doc when it gets refactored. I don't see any info on it at all
in either the Stan or CmdStan doc (other than the CmdStan doc
for delta, which only says it's the target).

- Bob

> --
> You received this message because you are subscribed to the Google Groups "stan development mailing list" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to stan-dev+u...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.

[Michael Betancourt]

This is just a consequence of not warming up long enough.  The dual averaging

is optimizing too strongly in the tail yielding a shorter step size than is necessary

in the bulk.  This can be remedied by tweaking the warmup adaptation window

sizes or the dual averaging parameters but it’s really hard to select something

that works better as a default.  If you warmup for longer you should see everything

converge to 0.8.

[Ben Goodrich]

On Sunday, May 18, 2014 4:32:00 PM UTC-4, Michael Betancourt wrote:

This is just a consequence of not warming up long enough.  The dual averaging

is optimizing too strongly in the tail yielding a shorter step size than is necessary

in the bulk.  This can be remedied by tweaking the warmup adaptation window

sizes or the dual averaging parameters but it’s really hard to select something

that works better as a default.  If you warmup for longer you should see everything

converge to 0.8.

If so, shouldn't we change the defaults to be something other than warmup being half the iterations and iterations defaulting to 2000? If the warmup window is almost always too short for almost all of these simple models, then it is probably way too short for complicated models.

Ben

[Bob Carpenter]

I think it's a tradeoff here. Adapting accept rate perfectly isn't
something we care about directly --- if the accept rate's too high,
then it's a bit ineffcient, but won't be wrong. And the question is
whether running more warmup iterations will lead to a greater overall
efficiency.

Yet another research project for someone to do some evals on :-)

- Bob

[Ben Goodrich]

On Monday, May 19, 2014 5:56:46 AM UTC-4, Bob Carpenter wrote:

I think it's a tradeoff here.  Adapting accept rate perfectly isn't
something we care about directly --- if the accept rate's too high,
then it's a bit ineffcient, but won't be wrong.  And the question is
whether running more warmup iterations will lead to a greater overall
efficiency.

Yet another research project for someone to do some evals on :-)

It's not much closer to the target of 0.8 with 10000 iterations and half warmup (see below).

Multiplying the number of iterations by 5 in order to get twice as many effective samples is not going to improve n_eff / second. But I continue to think that trying to make the defaults maximize n_eff / second for the largest class of models is not what defaults should be trying to accomplish. They should try to maximize the minimum n_eff / second over a wide class of models and I don't know if having an realized acceptance rate so close to 1 accomplishes that. More worrying, we have a few acceptance rates that are essentially zero, which is the sort of disaster that default arguments should be trying to prevent.

Ben

> results2
                                chain:1 chain:2 chain:3 chain:4
binormal                         0.9372  0.9802  0.9578  0.9564
inv_wishart                      0.9470  0.0098  0.9338  0.7136
normal_mixture                   0.9556  0.8246  0.9518  0.7806
normal                           0.8968  0.9484  0.9514  0.9606
ordered_pair                     0.9706  0.9474  0.9476  0.9436
sorted_pair                      0.9386  0.9464  0.9408  0.9480
triangle                         0.9472  0.9154  0.9192  0.9330
uniform                          0.8868  0.9268  0.9382  0.9344
wishart2                         0.9922  0.9946  0.9952  0.9944
wishart2x2                       0.9708  0.9608  0.9778  0.9712
wishart                          0.9672  0.9764  0.9250  0.9660
bernoulli                        0.9718  0.9436  0.9498  0.9622
negative_binomial2               0.9796  1.0000  0.9952  1.0000
negative_binomial                0.9568  0.9522  0.9646  0.9380
normal_censored                  0.9794  0.9444  0.9396  0.9720
normal_loc                       0.9578  0.9524  0.9460  0.9506
normal_mixture_k_prop            0.9980  0.9946  0.9982  0.9926
normal_mixture_k                 0.9898  0.9890  0.9894  0.9856
normal_mixture                   0.9822  0.9706  0.9696  0.9762
normal_truncated                 0.9662  0.9308  0.9754  0.9676
blocker                          0.8328  0.8232  0.9462  0.9728
bones                            0.9740  0.9820  0.9744  0.9562
dogs                             0.9768  0.9668  0.9658  0.9480
dyes                             0.8954  0.6226  0.9586  0.9630
epil                             0.9998  0.9998  0.9996  1.0000
equiv                            0.9862  0.9880  0.9740  0.9944
inhalers                         0.9914  0.9982  0.9978  0.9922
kidney                           0.9930  0.9714  0.9968  0.9774
leukfr                           0.9956  0.9854  0.8868  0.9958
leuk                             0.9856  0.9910  0.9800  0.9766
lsat                             0.9942  0.9962  0.9970  0.9946
magnesium                        0.9346  0.8206  0.8776  0.8870
mice                             0.9994  0.9994  0.9992  0.9990
oxford                           0.9900  0.9848  0.9930  0.9414
pump                             0.9838  0.9506  0.9904  0.9762
rats                             0.9958  0.9906  0.9824  0.9816
rats_vec                         0.9894  0.9826  0.9796  0.9860
rats_vec_unit                    0.9958  0.9912  0.9870  0.9914
salm2                            1.0000  1.0000  0.0002  1.0000
salm                             0.0002  0.0002  0.0002  0.0002
seeds                            0.9652  0.9684  0.9824  0.9750
stacks_a_normal                  0.9790  0.9840  0.9820  0.9862
stacks_b_dexp                    0.9896  0.9828  0.9816  0.9686
stacks_c_t4                      0.9488  0.9480  0.9694  0.9640
stacks_d_normal_ridge            0.9906  0.9920  0.9834  0.9710
stacks_e_dexp_ridge              0.9788  0.9848  0.9802  0.9750
stacks_f_t4_ridge                0.9722  0.9702  0.9684  0.9790
surgical_simple                  0.9704  0.9732  0.9902  0.9744
surgical                         0.9738  0.9680  0.9888  0.9630
air                              0.9912  0.9918  0.9882  0.9908
alli2                            0.9948  0.9974  0.9986  0.9954
beetles_cloglog                  0.9564  0.9810  0.9638  0.9752
beetles_logit                    0.9776  0.9606  0.9468  0.9760
beetles_probit                   0.9656  0.9326  0.9530  0.9576
birats                           0.9998  0.9948  1.0000  0.9998
cervix2                          0.9964  0.9924  0.9902  0.9888
dugongs                          0.9820  0.9732  0.9850  0.9888
endo1                            0.9436  0.9296  0.9370  0.9738
endo2                            0.9422  0.9224  0.9570  0.9770
endo3                            0.9896  0.9954  0.9930  0.9944
eyes                             0.9586  0.9758  0.9378  0.9764
hearts                           0.9322  0.9630  0.9664  0.9680
ice                              0.9444  0.9832  0.9776  0.7782
jaws                             0.9924  0.9956  0.9882  0.9928
mvn_orange                       0.9982  0.9778  1.0000  0.9930
orange                           0.9994  0.9906  0.9464  0.9696
schools                          0.9992  0.0002  0.9996  1.0000
stagnant2                        0.9880  0.9834  0.9942  0.9858
stagnant                         0.9872  0.9730  0.9936  0.9728
estdof                           0.9334  0.9556  0.9654  0.9614
simt4                            0.9930  0.9860  0.9908  0.9984
camel2                           0.9360  0.9488  0.9524  0.8884
camel                            0.9784  0.9808  0.9242  0.9840
seeds                            0.9876  0.9962  0.9942  0.9946
fire                                 NA      NA      NA      NA
circle                           0.9556  0.9426  0.9452  0.9398
hsquare                          0.9236  0.9484  0.9300  0.9546
parallelogram                    0.9458  0.9604  0.9342  0.9632
ring                             0.9250  0.9644  0.9290  0.9090
squaremc                         0.8818  0.9868  0.8896  0.9916
hepatitisME                      0.9184  0.5730  0.9972  1.0000
hepatitis                        0.9840  1.0000  1.0000  0.9652
command                          0.8968  0.9484  0.9514  0.9606
corr-lda                         0.8698  0.9268  0.9124  0.8322
lda                              0.9674  0.9582  0.9724  0.9824
naive-bayes                      0.9820  0.9800  0.9850  0.9750
naive-bayes-unsup                0.9884  0.9868  0.9894  0.9794
soft-k-means                     0.9816  0.9718  0.9960  0.9940

fx_equicorr                          NA      NA      NA      NA

fx_factor                        0.9112  0.9714  0.9286  0.9160
eight_schools                    0.8486  0.8582  0.8376  0.7942
funnel_reparam                   0.9700  0.9718  0.9322  0.9772
funnel                           0.8506  0.6342  0.4318  0.8662

arch1                                NA      NA      NA      NA
garch1_1                             NA      NA      NA      NA

koyck                            0.9820  0.9824  0.9906  0.9882
gp-fit                           0.9702  0.9762  0.9616  0.9868
gp-logit-predict                 0.9926  0.9898  0.9956  0.9932
gp-multi-sim                         NA      NA      NA      NA
gp-predict-analytic              0.9894  0.9932  0.9902  0.9906
gp-predict                       0.9926  0.9968  0.9974  0.9922
gp-sim-cholesky                  0.9910  0.9888  0.9950  0.9900
gp-sim                           0.9940  0.9888  0.9956  0.9920
hmm-analytic                     0.9934  0.9920  0.9816  0.9954
hmm-semisup                      0.9922  0.9976  0.9980  0.9958
hmm                              0.9958  0.9886  0.9852  0.9874
hmm-sufficient                   0.9900  0.9956  0.9928  0.9900
irt2_multilevel                  0.9998  0.9996  0.9996  0.9996
irt2                             0.9990  0.9988  1.0000  0.9968
irt_multilevel                   0.9998  0.9996  0.9998  0.9996
irt                              0.9998  0.9992  0.9998  0.9992

regression                           NA      NA      NA      NA
regression_std                       NA      NA      NA      NA

arma11_alt                       0.0002  0.9726  0.9894  0.9892
arma11                           0.0002  0.9912  0.9706  0.0002
ma2                              0.0002  0.9758  0.9574  0.0002

maQ                                  NA      NA      NA      NA

stochastic-volatility-optimized  0.9838  0.9784  0.9914  0.9832
stochastic-volatility            0.9984  0.9952  0.9978  0.9972
nnmf                             0.9980  0.9918  0.9938  0.9906
nnmf_vec                         0.9978  0.9974  0.9972  0.9972
logistic                         0.9796  0.9548  0.9736  0.9546
bounded_double                   0.9324  0.9292  0.8874  0.9462
correlation_matrix               0.9614  0.9308  0.9408  0.9198
simplex                          0.9488  0.9234  0.9648  0.9554
unit_vector                      0.9402  0.9718  0.8908  0.9744
> results2
                                chain:1 chain:2 chain:3 chain:4
binormal                         0.9372  0.9802  0.9578  0.9564
inv_wishart                      0.9470  0.0098  0.9338  0.7136
normal_mixture                   0.9556  0.8246  0.9518  0.7806
normal                           0.8968  0.9484  0.9514  0.9606
ordered_pair                     0.9706  0.9474  0.9476  0.9436
sorted_pair                      0.9386  0.9464  0.9408  0.9480
triangle                         0.9472  0.9154  0.9192  0.9330
uniform                          0.8868  0.9268  0.9382  0.9344
wishart2                         0.9922  0.9946  0.9952  0.9944
wishart2x2                       0.9708  0.9608  0.9778  0.9712
wishart                          0.9672  0.9764  0.9250  0.9660
bernoulli                        0.9718  0.9436  0.9498  0.9622
negative_binomial2               0.9796  1.0000  0.9952  1.0000
negative_binomial                0.9568  0.9522  0.9646  0.9380
normal_censored                  0.9794  0.9444  0.9396  0.9720
normal_loc                       0.9578  0.9524  0.9460  0.9506
normal_mixture_k_prop            0.9980  0.9946  0.9982  0.9926
normal_mixture_k                 0.9898  0.9890  0.9894  0.9856
normal_mixture                   0.9822  0.9706  0.9696  0.9762
normal_truncated                 0.9662  0.9308  0.9754  0.9676
blocker                          0.8328  0.8232  0.9462  0.9728
bones                            0.9740  0.9820  0.9744  0.9562
dogs                             0.9768  0.9668  0.9658  0.9480
dyes                             0.8954  0.6226  0.9586  0.9630
epil                             0.9998  0.9998  0.9996  1.0000
equiv                            0.9862  0.9880  0.9740  0.9944
inhalers                         0.9914  0.9982  0.9978  0.9922
kidney                           0.9930  0.9714  0.9968  0.9774
leukfr                           0.9956  0.9854  0.8868  0.9958
leuk                             0.9856  0.9910  0.9800  0.9766
lsat                             0.9942  0.9962  0.9970  0.9946
magnesium                        0.9346  0.8206  0.8776  0.8870
mice                             0.9994  0.9994  0.9992  0.9990
oxford                           0.9900  0.9848  0.9930  0.9414
pump                             0.9838  0.9506  0.9904  0.9762
rats                             0.9958  0.9906  0.9824  0.9816
rats_vec                         0.9894  0.9826  0.9796  0.9860
rats_vec_unit                    0.9958  0.9912  0.9870  0.9914
salm2                            1.0000  1.0000  0.0002  1.0000
salm                             0.0002  0.0002  0.0002  0.0002
seeds                            0.9652  0.9684  0.9824  0.9750
stacks_a_normal                  0.9790  0.9840  0.9820  0.9862
stacks_b_dexp                    0.9896  0.9828  0.9816  0.9686
stacks_c_t4                      0.9488  0.9480  0.9694  0.9640
stacks_d_normal_ridge            0.9906  0.9920  0.9834  0.9710
stacks_e_dexp_ridge              0.9788  0.9848  0.9802  0.9750
stacks_f_t4_ridge                0.9722  0.9702  0.9684  0.9790
surgical_simple                  0.9704  0.9732  0.9902  0.9744
surgical                         0.9738  0.9680  0.9888  0.9630
air                              0.9912  0.9918  0.9882  0.9908
alli2                            0.9948  0.9974  0.9986  0.9954
beetles_cloglog                  0.9564  0.9810  0.9638  0.9752
beetles_logit                    0.9776  0.9606  0.9468  0.9760
beetles_probit                   0.9656  0.9326  0.9530  0.9576
birats                           0.9998  0.9948  1.0000  0.9998
cervix2                          0.9964  0.9924  0.9902  0.9888
dugongs                          0.9820  0.9732  0.9850  0.9888
endo1                            0.9436  0.9296  0.9370  0.9738
endo2                            0.9422  0.9224  0.9570  0.9770
endo3                            0.9896  0.9954  0.9930  0.9944
eyes                             0.9586  0.9758  0.9378  0.9764
hearts                           0.9322  0.9630  0.9664  0.9680
ice                              0.9444  0.9832  0.9776  0.7782
jaws                             0.9924  0.9956  0.9882  0.9928
mvn_orange                       0.9982  0.9778  1.0000  0.9930
orange                           0.9994  0.9906  0.9464  0.9696
schools                          0.9992  0.0002  0.9996  1.0000
stagnant2                        0.9880  0.9834  0.9942  0.9858
stagnant                         0.9872  0.9730  0.9936  0.9728
estdof                           0.9334  0.9556  0.9654  0.9614
simt4                            0.9930  0.9860  0.9908  0.9984
camel2                           0.9360  0.9488  0.9524  0.8884
camel                            0.9784  0.9808  0.9242  0.9840
seeds                            0.9876  0.9962  0.9942  0.9946
fire                                 NA      NA      NA      NA
circle                           0.9556  0.9426  0.9452  0.9398
hsquare                          0.9236  0.9484  0.9300  0.9546
parallelogram                    0.9458  0.9604  0.9342  0.9632
ring                             0.9250  0.9644  0.9290  0.9090
squaremc                         0.8818  0.9868  0.8896  0.9916
hepatitisME                      0.9184  0.5730  0.9972  1.0000
hepatitis                        0.9840  1.0000  1.0000  0.9652
command                          0.8968  0.9484  0.9514  0.9606
corr-lda                         0.8698  0.9268  0.9124  0.8322
lda                              0.9674  0.9582  0.9724  0.9824
naive-bayes                      0.9820  0.9800  0.9850  0.9750
naive-bayes-unsup                0.9884  0.9868  0.9894  0.9794
soft-k-means                     0.9816  0.9718  0.9960  0.9940

fx_equicorr                          NA      NA      NA      NA

fx_factor                        0.9112  0.9714  0.9286  0.9160
eight_schools                    0.8486  0.8582  0.8376  0.7942
funnel_reparam                   0.9700  0.9718  0.9322  0.9772
funnel                           0.8506  0.6342  0.4318  0.8662

arch1                                NA      NA      NA      NA
garch1_1                             NA      NA      NA      NA

koyck                            0.9820  0.9824  0.9906  0.9882
gp-fit                           0.9702  0.9762  0.9616  0.9868
gp-logit-predict                 0.9926  0.9898  0.9956  0.9932
gp-multi-sim                         NA      NA      NA      NA
gp-predict-analytic              0.9894  0.9932  0.9902  0.9906
gp-predict                       0.9926  0.9968  0.9974  0.9922
gp-sim-cholesky                  0.9910  0.9888  0.9950  0.9900
gp-sim                           0.9940  0.9888  0.9956  0.9920
hmm-analytic                     0.9934  0.9920  0.9816  0.9954
hmm-semisup                      0.9922  0.9976  0.9980  0.9958
hmm                              0.9958  0.9886  0.9852  0.9874
hmm-sufficient                   0.9900  0.9956  0.9928  0.9900
irt2_multilevel                  0.9998  0.9996  0.9996  0.9996
irt2                             0.9990  0.9988  1.0000  0.9968
irt_multilevel                   0.9998  0.9996  0.9998  0.9996
irt                              0.9998  0.9992  0.9998  0.9992

regression                           NA      NA      NA      NA
regression_std                       NA      NA      NA      NA

arma11_alt                       0.0002  0.9726  0.9894  0.9892
arma11                           0.0002  0.9912  0.9706  0.0002
ma2                              0.0002  0.9758  0.9574  0.0002

maQ                                  NA      NA      NA      NA

stochastic-volatility-optimized  0.9838  0.9784  0.9914  0.9832
stochastic-volatility            0.9984  0.9952  0.9978  0.9972
nnmf                             0.9980  0.9918  0.9938  0.9906
nnmf_vec                         0.9978  0.9974  0.9972  0.9972
logistic                         0.9796  0.9548  0.9736  0.9546
bounded_double                   0.9324  0.9292  0.8874  0.9462
correlation_matrix               0.9614  0.9308  0.9408  0.9198
simplex                          0.9488  0.9234  0.9648  0.9554
unit_vector                      0.9402  0.9718  0.8908  0.9744

[Bob Carpenter]

I think the main aim was and should continue to be
to maximize robustness without too much sacrifice in speed.

Smaller step sizes and higher acceptance rates are more robust,
but slower. It would be nice if the user could actually control
the tradeoff in a meaningful way, but it doesn't look like
that's happening.

Too high acceptance rates means we're erring on the side of
caution, I think. So maybe we want to aim more toward speed?

What are the NA values? Is that 0 acceptance rate?

It looks to me that the adaptation isn't that robust in the sense
that it's converging to different mass matrices and/or step sizes
in different chains leading to different acceptance rates.

I think some of the models (like maybe salm) are just poorly
translated from BUGS or in some other way problematic. We've
had trouble with lots of these before.

- Bob

[Ben Goodrich]

NA means the model didn't run for some reason. A realized acceptance rate somewhat above 0.8 is probably not too dangerous, but an acceptance rate of virtually 1 could easily not cover the tails. The near zeros on something as simple as inv_wishart are a bigger concern though.

Ben

 

[Bob Carpenter]

Agreed. This should be investigated. Maybe after this release
is out? Any ideas, Michael?

- Bob

On May 20, 2014, at 9:27 PM, Ben Goodrich <goodri...@gmail.com> wrote:

> NA means the model didn't run for some reason. A realized acceptance rate somewhat above 0.8 is probably not too dangerous, but an acceptance rate of virtually 1 could easily not cover the tails. The near zeros on something as simple as inv_wishart are a bigger concern though.
>
> Ben
>
>
>

[Michael Betancourt]

We need to separate out a few issues.

Firstly, the adaptation works fine one models we know HMC should fit, with the target
value of 0.8 quickly achieved.

Secondly, we can’t expect that success to carry over to all of the BUGS models, most
of which break every MCMC algorithm out there. With the nominal parameterizations
none of these are going to be geometrically ergodic and adaptation will be brittle as
a consequence. For example, step sizes that yield an average acceptance probability
of 0.8 in the bulk may cause way too many divergences in the tail (which we treat as
rejections, and hence skew the acceptance probability). The flip size are cases
where adaptation goes crazy during warmup (as is common without geometric
ergodicity) and the acceptance probability vanishes.

Before worrying about the adaptation we have to worry about the well-posedness
of all of these models. They need to be reparameterized and checked for heavy
tails (a la the HMC for hierarchical models paper). Once we know they fit then
we can see if the adaptation is still causing problems.

[Bob Carpenter]

Now I'm worried (and maybe I should've been all along).

I think this means we need to radically rethink our manual.
We do not want to be providing example models (translations of
BUGS models or otherwise), which should not be used.

And when you say won't fit, what's going to go wrong in terms
of sampling? Are we not going to be sampling from the posterior or
is it going to be so slow that the effect is the same?

- Bob

[Michael Betancourt]

> Now I'm worried (and maybe I should've been all along).

Most people don’t take geometric ergodicity seriously, but it’s
a big deal in practice and manifests as all of these pathologies
we keep seeing when people take careful looks at their models.

> I think this means we need to radically rethink our manual.
> We do not want to be providing example models (translations of
> BUGS models or otherwise), which should not be used.

I agree — but it’s hard to abandon them completely. Will users
be content when we say “BUGS/JAGS has never been fitting these
well, you just never had the diagnostics to see it”? Or, “we’ve
rewritten these examples so only compare these and ignore all
of the others”?

> And when you say won't fit, what's going to go wrong in terms
> of sampling? Are we not going to be sampling from the posterior or
> is it going to be so slow that the effect is the same?

It depends on your perspective. These algorithms are all ergodic,
just not geometrically ergodic. Technically this means that they will
eventually explore the entire posterior BUT we no longer have the
Monte Carlo Central Limit Theorem so the accuracy of the Monte
Carlo estimates is unknown. For example, biases estimates
of hyperparameter means, as we see in hierarchical models.

[Ben Goodrich]

On Wednesday, May 21, 2014 12:16:38 PM UTC-4, Michael Betancourt wrote:

Before worrying about the adaptation we have to worry about the well-posedness
of all of these models.  They need to be reparameterized and checked for heavy
tails (a la the HMC for hierarchical models paper).  Once we know they fit then
we can see if the adaptation is still causing problems.

I've been saying for a while that the BUGS models are only useful for demonstrating features of the language. But the slowness or inability to get close to the target acceptance rate holds just as much for the non-BUGS models where we are just sampling from a distribution (e.g. inv_wishart) or are generating the data and fitting the model (everything below funnel_reparam in that matrix). And I think funnel_reparam is a model where the chains ostensibly should converge quickly to the target and they don't with 10000 iterations

funnel_reparam: 0.9700  0.9718  0.9322  0.9772


Ben

[Michael Betancourt]

inv_wishart has notoriously heavy tails and runs into the same ergodicity problem.

funnel_reparam is literally a 2D gaussian and should be investigated.

[Marcus Brubaker]

In case anyone is interested, I found this article http://www.i-journals.org/ps/viewarticle.php?id=29 to be a relatively informative review of MCMC CLT results.

Cheers,

Marcus

On Wed, May 21, 2014 at 12:41 PM, Michael Betancourt <betan...@gmail.com> wrote:

> Now I'm worried (and maybe I should've been all along).

Most people don't take geometric ergodicity seriously, but it's
a big deal in practice and manifests as all of these pathologies
we keep seeing when people take careful looks at their models.

> I think this means we need to radically rethink our manual.
> We do not want to be providing example models (translations of
> BUGS models or otherwise), which should not be used.

I agree -- but it's hard to abandon them completely.  Will users

be content when we say "BUGS/JAGS has never been fitting these
well, you just never had the diagnostics to see it"?  Or, "we've
rewritten these examples so only compare these and ignore all
of the others"?

> And when you say won't fit, what's going to go wrong in terms
> of sampling?  Are we not going to be sampling from the posterior or
> is it going to be so slow that the effect is the same?

It depends on your perspective.  These algorithms are all ergodic,
just not geometrically ergodic.  Technically this means that they will
eventually explore the entire posterior BUT we no longer have the
Monte Carlo Central Limit Theorem so the accuracy of the Monte
Carlo estimates is unknown.  For example, biases estimates
of hyperparameter means, as we see in hierarchical models.

[Michael Betancourt]

The hardest part about all of this stuff is the dearth of research on it.

As Marcus notes Jones’ stuff is really good and probably the most

readable.  Outside of that you’d want to look at Gareth Roberts’ work,

but the level of math increases exponentially.

[Bob Carpenter]

On May 21, 2014, at 7:11 PM, Marcus Brubaker <marcus.a...@gmail.com> wrote:

> In case anyone is interested, I found this article http://www.i-journals.org/ps/viewarticle.php?id=29 to be a relatively informative review of MCMC CLT results.
>
> Cheers,
> Marcus
>
>
> On Wed, May 21, 2014 at 12:41 PM, Michael Betancourt <betan...@gmail.com> wrote:
>
> > Now I'm worried (and maybe I should've been all along).
>
> Most people don't take geometric ergodicity seriously, but it's
> a big deal in practice and manifests as all of these pathologies
> we keep seeing when people take careful looks at their models.
>
> > I think this means we need to radically rethink our manual.
> > We do not want to be providing example models (translations of
> > BUGS models or otherwise), which should not be used.
>
> I agree -- but it's hard to abandon them completely. Will users
> be content when we say "BUGS/JAGS has never been fitting these
> well, you just never had the diagnostics to see it"?

Yes, if we can show them those diagnostics. If they don't
accept it, there's not much we can do.

Is there a good diagnostic that would show these not working?
In such a way that we can quantify the differences in inferences
vs. a working sampler?

> Or, "we've
> rewritten these examples so only compare these and ignore all
> of the others"?

The diagnostics will then motivate the rewriting of the examples
so they do work, if that's possible.

If not, maybe we just give up MCMC and spend all of our time
on optimization :-)

- Bob

[Ben Goodrich]

On Wednesday, May 21, 2014 1:40:22 PM UTC-4, Bob Carpenter wrote:

> I agree -- but it's hard to abandon them completely.  Will users
> be content when we say "BUGS/JAGS has never been fitting these
> well, you just never had the diagnostics to see it"?

Yes, if we can show them those diagnostics.  If they don't
accept it, there's not much we can do.

Is there a good diagnostic that would show these not working?
In such a way that we can quantify the differences in inferences
vs. a working sampler?

I don't think that many of them just don't work

https://groups.google.com/d/msg/stan-dev/jEZvT9sOwnA/A4uCIQGsH5gJ

The ones that I classify as having bad mixing at a minimum are ones whose dependence plot looks like the leukfr example at the end of

https://groups.google.com/d/msg/stan-dev/_SJBiCaUbt4/I_G38zlHGmoJ

When the accepted draws closer together than they ideally would be, you should get a hockey stick for a dependence plot like

https://groups.google.com/d/msg/stan-dev/89GSQ9wrIQI/KBMKxvYYW8MJ

where it takes a while for the chains to get to the same distribution but there is not that much dependence between the draws thereafter. But the hockey stick examples aren't of very much concern, compared to the situation where you can't get the within-chain dependence to go to zero for any reasonable number of lags.

Ben

[Aki Vehtari]

Bob:

> > I think this means we need to radically rethink our manual.
> > We do not want to be providing example models (translations of
> > BUGS models or otherwise), which should not be used.

Mike:

> I agree — but it’s hard to abandon them completely. Will users
> be content when we say “BUGS/JAGS has never been fitting these
> well, you just never had the diagnostics to see it”? Or, “we’ve
> rewritten these examples so only compare these and ignore all
> of the others”?

I think manual should tell what you want people to do. Technical report could give details (including sampler comparisons) why the model was rewritten. There are lot of people who have never used BUGS/JAGS and when they start using Stan they will suffer if the manual is full of bad legacy examples. When I was teaching Stan to some students, I had a problem that I had to explain why the example models were not good and the current knowledge on how models should be constructed is different. You could put BUGS/JAGS examples in a separate chapter with instructions how to port them to Stan.

Aki

[Bob Carpenter]

Absolutely. It's echoing what Andrew's been saying all along and what
I've been trying to update the doc to reflect.

Part of the issue wasn't so much pedagogical as we didn't have the better
ways to do things implemented or figured out when the manual was written!

At this point, I think it could use a thorough rewrite and also the BUGS
models could use a thorough rewrite into best practices for both computation
and statistical modeling.

- Bob

[Michael Betancourt]

I agree with Bob — much of the confusing arises from the fact

that we’re learning on the job.  Speaking for myself my current

research into the more theoretical perspectives like ergodicity

brings new understanding into common pathologies and motivates

new strategies.  The very construction of Stan is fueling new

understanding of both HMC and MCMC!

The manual is due for a big update and we’ve made the first 

steps towards that.  Hopefully out outreach efforts (speaking,

meet ups, etc) will also contribute to getting the most recent

principles out sooner.

On May 23, 2014, at 2:40 PM, Bob Carpenter <ca...@alias-i.com> wrote:

On May 23, 2014, at 3:08 PM, Aki Vehtari <aki.v...@aalto.fi> wrote:

Bob:

I think this means we need to radically rethink our manual.
We do not want to be providing example models (translations of
BUGS models or otherwise), which should not be used.

Mike:

I agree -- but it's hard to abandon them completely.  Will users

be content when we say "BUGS/JAGS has never been fitting these
well, you just never had the diagnostics to see it"?  Or, "we've
rewritten these examples so only compare these and ignore all
of the others"?

I think manual should tell what you want people to do. Technical report could give details (including sampler comparisons) why the model was rewritten. There are lot of people who have never used BUGS/JAGS and when they start using Stan they will suffer if the manual is full of bad legacy examples. When I was teaching Stan to some students, I had a problem that I had to explain why the example models were not good and the current knowledge on how models should be constructed is different. You could put BUGS/JAGS examples in a separate chapter with instructions how to port them to Stan.

Absolutely.  It's echoing what Andrew's been saying all along and what
I've been trying to update the doc to reflect. 

Part of the issue wasn't so much pedagogical as we didn't have the better
ways to do things implemented or figured out when the manual was written!

At this point, I think it could use a thorough rewrite and also the BUGS
models could use a thorough rewrite into best practices for both computation
and statistical modeling. 

- Bob

[Ben Goodrich]

On Wednesday, May 21, 2014 1:09:12 PM UTC-4, Michael Betancourt wrote:

inv_wishart has notoriously heavy tails and runs into the same ergodicity problem.

funnel_reparam is literally a 2D gaussian and should be investigated.

Any chance we can try to fix this before 2.4.0? I don't know of any model that gets close to the targeted acceptance rate of 0.8 in a reasonable time. At the JSM presentation, I am supposed to be telling them how efficient Stan is for hierarchical models, so it would be nice if the acceptance rates weren't systematically high. Even for sampling a standard normal with 1 million warmup iterations, it is way above target.

library(rstan)
library(parallel)

mc <- '
parameters {
  real theta;
}
model {
  theta ~ normal(0,1);
}
'

fit <- stan(model_code = mc, chains = 0)
AR <- unlist(mclapply(1:8, mc.cores = 8L, FUN = function(i) {
  length(unique(extract(stan(fit = fit, warmup = 10^6, iter = 10^6 + 1000,
                             chains = 1, seed = 12345, chain_id = i, refresh = 0),
          pars = "theta", permuted = FALSE))) / 1000
}))
AR
[1] 0.956 0.953 0.894 0.951 0.944 0.947 0.962 0.945

Ben

[Michael Betancourt]

Well, it would help if you understood how adaptation works.

As in documented in the manual, adaptation proceeds in three stages: an initial window where the metric is not updated, an intermediate window where the metric and step size are updated together, and a terminal window where the metric is set and the step size is given some time to find a good value conditioned on the final metric.  What's happening here is that for a simple model the terminal window is too small to compute a reasonable estimate of the average acceptance probability.  In particular, the estimated average is biased towards initial values and smaller acceptance probabilities (both because of who the Nesterov dual averaging works and because the distribution of acceptance probabilities is highly skewed) which is why the final step size is too conservative.

There are two simple fixes for a simple model: either increase the terminal window with adapt term_buffer=100 (the default is 75) or mess around with the dual averaging parameters.  One way to deemphasize the early values, for example, is to increase the gamma parameter with adapt gamma=0.15.  Both of these work for the Gaussian case, but without having an ensemble of models I'm not going to advocate changing any of the defaults.  We're on the safe side right now.

[Ben Goodrich]

On Thursday, July 10, 2014 6:18:50 AM UTC-4, Michael Betancourt wrote:

Well, it would help if you understood how adaptation works.

I indeed do not understand how the estimated average acceptance probability could be biased towards initial values when there are 1 million warmup iterations.
 

There are two simple fixes for a simple model: either increase the terminal window with adapt term_buffer=100 (the default is 75) or mess around with the dual averaging parameters.  One way to deemphasize the early values, for example, is to increase the gamma parameter with adapt gamma=0.15.  Both of these work for the Gaussian case, but without having an ensemble of models I'm not going to advocate changing any of the defaults.  We're on the safe side right now.

In the cmdstan documentation, it says adapt term_buffer defaults to 50 and in the rstan documentation it says it defaults to 25. Is it supposed to default to 75? In any event, neither changing term_buffer to 100 nor gamma to 0.15 separately or together works for the this Gaussian case:

> AR <- unlist(mclapply(1:8, mc.cores = 8L, FUN = function(i) {

+   length(unique(extract(stan(fit = fit, warmup = 10^6, iter = 10^6 + 1000,
+                              control = list(adapt_term_buffer = 100),
+                              chains = 1, seed = 12345, chain_id = i, refresh = 0),
+           pars = "theta", permuted = FALSE))) / 1000
+ }))
> AR
[1] 0.927 0.932 0.924 0.914 0.946 0.923 0.917 0.925

> AR <- unlist(mclapply(1:8, mc.cores = 8L, FUN = function(i) {

+   length(unique(extract(stan(fit = fit, warmup = 10^6, iter = 10^6 + 1000,
+                              control = list(adapt_gamma = 0.15),
+                              chains = 1, seed = 12345, chain_id = i, refresh = 0),
+           pars = "theta", permuted = FALSE))) / 1000
+ }))
> AR
[1] 0.884 0.898 0.832 0.903 0.916 0.901 0.837 0.897

> AR <- unlist(mclapply(1:8, mc.cores = 8L, FUN = function(i) {

+   length(unique(extract(stan(fit = fit, warmup = 10^6, iter = 10^6 + 1000,
+                              control = list(adapt_term_buffer = 100,
+                                             adapt_gamma = 0.15),
+                              chains = 1, seed = 12345, chain_id = i, refresh = 0),
+           pars = "theta", permuted = FALSE))) / 1000
+ }))
> AR
[1] 0.856 0.822 0.852 0.879 0.907 0.894 0.869 0.905

It gets a little closer if I tweak it to 500000 term_buffer and 0.75 gamma

> AR <- unlist(mclapply(1:8, mc.cores = 8L, FUN = function(i) {

+   length(unique(extract(stan(fit = fit, warmup = 10^6, iter = 10^6 + 1000,
+                              control = list(adapt_term_buffer = 500000,
+                                             adapt_gamma = 0.75),
+                              chains = 1, seed = 12345, chain_id = i, refresh = 0),
+           pars = "theta", permuted = FALSE))) / 1000
+ }))
> AR
[1] 0.853 0.865 0.844 0.859 0.858 0.842 0.834 0.852

Perhaps it would be easier to say that if you want to target an acceptance rate of 0.8, then you should target an acceptance rate of 0.6-ish

> AR <- unlist(mclapply(1:8, mc.cores = 8L, FUN = function(i) {

+   length(unique(extract(stan(fit = fit, warmup = 10^6, iter = 10^6 + 1000,
+                              control = list(adapt_delta = 0.63),
+                              chains = 1, seed = 12345, chain_id = i, refresh = 0),
+           pars = "theta", permuted = FALSE))) / 1000
+ }))
> AR
[1] 0.758 0.887 0.905 0.831 0.737 0.759 0.642 0.725
> mean(AR)
[1] 0.7805

Ben

[Michael Betancourt]

I indeed do not understand how the estimated average acceptance probability could be biased towards initial values when there are 1 million warmup iterations.

Because you are not understanding how adaptation works.  As I have explained multiple times, only the term_buffer iterations are used to estimate the final step size.

If you run 1 million warmup iterations you still only have a term_buffer of 100.

 

In the cmdstan documentation, it says adapt term_buffer defaults to 50 and in the rstan documentation it says it defaults to 25. Is it supposed to default to 75?

I misspoke -- it's 50.  Rstan sets it to whatever it wants to.

 

In any event, neither changing term_buffer to 100 nor gamma to 0.15 separately or together works for the this Gaussian case:

For a 10D IID Gaussian

./normal sample adapt term_buffer=100 gamma=0.15

works fine,

                    Mean     MCSE   StdDev    5%       50%   95%  N_Eff  N_Eff/s    R_hat

accept_stat__    7.8e-01  6.4e-03  2.0e-01  0.38   8.2e-01   1.0   1000    26413  1.0e+00

The dual averaging is very sensitive to the parameters, and pushing gamma up to huge values will not result in better results than a small increase.

 

[Ben Goodrich]

On Thursday, July 10, 2014 9:38:42 AM UTC-4, Michael Betancourt wrote:

I indeed do not understand how the estimated average acceptance probability could be biased towards initial values when there are 1 million warmup iterations.

Because you are not understanding how adaptation works.  As I have explained multiple times, only the term_buffer iterations are used to estimate the final step size.

If you run 1 million warmup iterations you still only have a term_buffer of 100.

Do you mean initial values within the term_buffer or initial values of the Markov Chain? I can see how increasing the term_buffer would yield a more precise estimate and if it is heavily skewed, then an imprecise estimate could exhibit the behavior we are seeing in the examples where many acceptance rates close to 1 and a few close to zero. But I don't see how the initial values of the Markov Chain would entail bias for a reasonable warmup size.
 

 In the cmdstan documentation, it says adapt term_buffer defaults to 50 and in the rstan documentation it says it defaults to 25. Is it supposed to default to 75?

I misspoke -- it's 50.  Rstan sets it to whatever it wants to.

I'll change it to 50 by default.
 

 In any event, neither changing term_buffer to 100 nor gamma to 0.15 separately or together works for the this Gaussian case:

For a 10D IID Gaussian

./normal sample adapt term_buffer=100 gamma=0.15

works fine,

                    Mean     MCSE   StdDev    5%       50%   95%  N_Eff  N_Eff/s    R_hat

accept_stat__    7.8e-01  6.4e-03  2.0e-01  0.38   8.2e-01   1.0   1000    26413  1.0e+00

Yours is working differently than mine. What commit are you working from?

mc <- '
parameters {
  real theta[10];
}
model {
  theta ~ normal(0,1);
}
'

> AR <- unlist(mclapply(1:8, mc.cores = 8L, FUN = function(i) {

+   nrow(unique(as.data.frame(
+     extract(stan(fit = fit, warmup = 10^6, iter = 10^6 + 1000,
+             control = list(adapt_term_buffer = 100, adapt_gamma = 0.15),

+             chains = 1, seed = 12345, chain_id = i, refresh = 0),

+           pars = "theta", permuted = FALSE)))) / 1000
+ }))
> AR
[1] 0.953 0.919 0.929 0.937 0.947 0.948 0.919 0.927

Ben

[Ben Goodrich]

On Thursday, July 10, 2014 10:07:00 AM UTC-4, Ben Goodrich wrote:

For a 10D IID Gaussian

./normal sample adapt term_buffer=100 gamma=0.15

works fine,

                    Mean     MCSE   StdDev    5%       50%   95%  N_Eff  N_Eff/s    R_hat

accept_stat__    7.8e-01  6.4e-03  2.0e-01  0.38   8.2e-01   1.0   1000    26413  1.0e+00

Yours is working differently than mine. What commit are you working from?

Upon further review (and using develop), the mean accept_stat__ is close to 0.8 for all 10 chains

goodrich@CYBERPOWERPC:/tmp/SU/AR$ find . -type f -name "*.csv" -exec /opt/stan/bin/print '{}' \; | grep "accept_stat__"
accept_stat__    8.1e-01  5.5e-03  1.7e-01  0.48   8.4e-01   1.0   1000    13579  1.0e+00
accept_stat__    7.9e-01  6.0e-03  1.9e-01  0.43   8.3e-01   1.0   1000    24491  1.0e+00
accept_stat__    8.2e-01  5.5e-03  1.7e-01  0.48   8.6e-01   1.0   1000    13700  1.0e+00
accept_stat__    8.1e-01  5.8e-03  1.8e-01  0.44   8.6e-01   1.0   1000    14148  1.0e+00
accept_stat__    8.2e-01  5.2e-03  1.7e-01  0.51   8.6e-01   1.0   1000    13257  1.0e+00
accept_stat__    7.9e-01  6.1e-03  1.9e-01  0.43   8.3e-01   1.0   1000    13747  1.0e+00
accept_stat__    8.2e-01  5.4e-03  1.7e-01  0.49   8.6e-01   1.0   1000    13388  1.0e+00
accept_stat__    7.9e-01  6.1e-03  1.8e-01  0.44   8.3e-01   1.0    925    12311  1.0e+00
accept_stat__    7.7e-01  6.4e-03  2.0e-01  0.40   8.0e-01   1.0   1000    15021  1.0e+00
accept_stat__    8.2e-01  5.5e-03  1.7e-01  0.48   8.6e-01   1.0   1000    13422  1.0e+00

but the proportion of unique draws (which is what I was calculating in R previously) is closer to 0.93

goodrich@CYBERPOWERPC:/tmp/SU/AR$ grep -v -e "^#" -e "^lp__" -e "^$" output*.csv | awk  -F "\"*,\"*" '{print $1}' | sort | uniq | wc -l
9287

I gather from the documentation that accept_stat__ is min(1, ratio) for the proposal that gets sliced. So, why does a process with a mean accept_stat__ of 0.8 yield an acceptance rate of 0.93?

Ben

[Michael Betancourt]

Do you mean initial values within the term_buffer or initial values of the Markov Chain? I can see how increasing the term_buffer would yield a more precise estimate and if it is heavily skewed, then an imprecise estimate could exhibit the behavior we are seeing in the examples where many acceptance rates close to 1 and a few close to zero. But I don't see how the initial values of the Markov Chain would entail bias for a reasonable warmup size.

No you misunderstanding.  After (n_warmup - term_buffer) iterations the metric is updated for the last time and then warmup continues of term_buffer iterations to learn the step size for the new metric.

The empirical average acceptance probability in this window is noisy, and tends to underestimate the true expectations because the distribution of acceptance probabilities is no where near normal, more

of a Beta (big number, small number).

Yours is working differently than mine. What commit are you working from?

Develop version of CmdStan and Stan. 

[Jiqiang Guo]

On Thu, Jul 10, 2014 at 10:07 AM, Ben Goodrich <goodri...@gmail.com> wrote:

 In the cmdstan documentation, it says adapt term_buffer defaults to 50 and in the rstan documentation it says it defaults to 25. Is it supposed to default to 75?

I misspoke -- it's 50.  Rstan sets it to whatever it wants to.

I'll change it to 50 by default.
 

I think it is just a documentation error.  For the sake of keeping a log, I will create an issue and fix it. 

Jiqiang