search

chriselrod / ProbabilityModels.jl

Probability models with `y ~ Normal(mu, sigma)` style syntax and reverse mode AD for gradients; WIP. Prioritizes runtime performance.
MIT License
20 stars 1 forks source link

"Bit array vectors with 4 < 8 not yet supported." #2

Closed goedman closed 5 years ago

goedman commented 5 years ago

Not sure if you want to take issues yet, no problem if not, I'll just close it.

I've installed all packages. Did an extra Pkg.build("VectorizationBase"), after installing all packages in .julia/dev.

Trying:

julia> include("/Users/rob/.julia/dev/DynamicHMCModels/scripts/ProbabilityModels/LR.jl")
    Defined model: LogisticRegressionModel.
    Unknowns: μ1, β1, y, μ0, σ0, β0, X, σ1.

386
ERROR: Error while loading expression starting at /Users/rob/.julia/dev/DynamicHMCModels/scripts/ProbabilityModels/LR.jl:29
caused by [exception 1]
"Bit array vectors with 4 < 8 not yet supported."
Stacktrace:
 [1] #s43#106 at /Users/rob/.julia/dev/SIMDPirates/src/pirate.jl:96 [inlined]
 [2] #s43#106(::Any, ::Any, ::Any, ::Any, ::Any) at ./none:0
 [3] (::Core.GeneratedFunctionStub)(::Any, ::Vararg{Any,N} where N) at ./boot.jl:524
 [4] macro expansion at /Users/rob/.julia/dev/LoopVectorization/src/LoopVectorization.jl:208 [inlined]
 [5] macro expansion at ./gcutils.jl:224 [inlined]
 [6] macro expansion at /Users/rob/.julia/dev/ProbabilityDistributions/src/distribution_functions.jl:179 [inlined]
...

Did you build a special version of Julia?

chriselrod commented 5 years ago

Thank you for the issue and interest! The problem should be fixed.

I was running it on the latest Julia master, but it was not otherwise special and should work on Julia 1.0 and newer.

The issue was that I had tested and optimized the code for avx512, and used UInt8s as bitmasks. With 512 bit vectors, it was calculating vectors of 8 log(p[i]) and log(1-p[i]) at a time, and loading a UInt8 from the BitVector y to use as a mask to select between them. There is no such thing as a UInt4, so this did not work with avx2. My fix was to use vectors of length 8 anyway, even with avx2 or sse. I have not tested it thoroughly, however, and it will be a bit less efficient. But my brief tests suggest it's still much faster than Stan.

Note that

a = X*beta
# negative inverse logit
# OmP = One minus P
OmP = 1 / (1 + exp(a) )
# P = 1 - OmP
logOmP = log(OmP)
logP = a + logOmP

so calculating both log(1-p) and log(p) is cheap, and requires only one logarithm (because a = log(p / (1 - p) ) = log(p) - log(1-p)).

Please let me know if you have any more problems.

goedman commented 5 years ago

Thanks Chris! Indeed this fixed it!

Just FYI, my other mistake was to try it on Julia 1.3- which gives:

julia> include("/Users/rob/.julia/dev/DynamicHMCModels/scripts/ProbabilityModels/LR.jl")
    Defined model: LogisticRegressionModel.
    Unknowns: μ1, β1, y, μ0, σ0, β0, X, σ1.

393
ERROR: Error while loading expression starting at /Users/rob/.julia/dev/DynamicHMCModels/scripts/ProbabilityModels/LR.jl:30
caused by [exception 1]
MethodError: vmuladd(::Float64, ::ConstantFixedSizePaddedArray{Tuple{5},Float64,1,8,8}, ::ConstantFixedSizePaddedArray{Tuple{5},Float64,1,8,8}) is ambiguous. Candidates:
  vmuladd(a::T, x::Union{Adjoint{T,#s39} where #s39<:(PaddedMatrices.AbstractFixedSizePaddedArray{Tuple{P1},T,1,L1,L} where L), #s40} where #s40<:(PaddedMatrices.AbstractFixedSizePaddedArray{Tuple{P1},T,1,L1,L} where L), y::Union{Adjoint{T,#s37} where #s37<:(PaddedMatrices.AbstractFixedSizePaddedArray{Tuple{P2},T,1,L2,L} where L), #s38} where #s38<:(PaddedMatrices.AbstractFixedSizePaddedArray{Tuple{P2},T,1,L2,L} where L)) where {T, P1, L1, P2, L2} in PaddedMatrices at /Users/rob/.julia/dev/PaddedMatrices/src/blas.jl:329
  vmuladd(a::T, x::PaddedMatrices.AbstractFixedSizePaddedArray{S,T,N,P,L}, y::PaddedMatrices.AbstractFixedSizePaddedArray{S,T,N,P,L}) where {S, T<:Number, N, P, L} in PaddedMatrices at /Users/rob/.julia/dev/PaddedMatrices/src/blas.jl:320
Possible fix, define
  vmuladd(::T<:Number, ::#s40 where #s40<:PaddedMatrices.AbstractFixedSizePaddedArray{Tuple{P1},T<:Number,1,P,L}, ::#s38 where #s38<:PaddedMatrices.AbstractFixedSizePaddedArray{Tuple{P2},T<:Number,1,P,L})

Julia 1.1 and 1.2 run fine.

chriselrod commented 5 years ago

I pushed a few more changes, including a fix for that. I built the latest master as of today (1.3) and tested it locally.

goedman commented 5 years ago

Thanks Chris! I tried 1.3 again but no luck. Could it be (on OSX);

WARNING: could not import LinearAlgebra.ARPACKException into Arpack

(just before DynamicHMC starts to run).

I've successfully redone the LogisticRegression example ( and indeed see a between 2 and 3x better performance with PM/dHMC) and the CmdStan part of ITP. See similar CmdStan ITP results. Will analyze and run the PM/dHMC versions next.

chriselrod commented 5 years ago

I have been seeing that warning too on occasion, but don't recall running into subsequent errors. Seems that's because of this commit. I haven't tried this yet, but switching to Arpack master should solve the problem. None of my libraries depend on Arpack as far as I know. Manifest.toml for this project, for example.

Distributions.jl requires PDMats.jl, which requires Arpack. I'm guessing something, somewhere, is loading Distributions.

Do you mean you ran the LogisticRegression example with both PM/dHMC and CmdStan, seeing 2-3x better performance with PM/dHMC? And that you ran the ITP model with CmdStan, but not the PM/dHMC version yet?

When you do, I'd be very interested in learning the results. Unfortunately not all chains converged in both. I'm not sure if I could change that via using different options, or if a more manual intervention is needed. I'm also of course interested in runtime performance.

Especially if the Stan model runs much faster than mine. I have another version of that Stan model that takes about 2 hours to run locally, but I was told by someone else who ran my code with some modifications (using rstan and drake) on a cluster said that all chains finished in < 20 minutes. The cluster has more cores than my computes of course, but I didn't use anything like map_rect, so I don't think the chains should finish any faster for hardware reasons. I did ask for more details. This suggests to me that different Stan versions, or C++ compilers/versions could make a substantial difference. I don't have the impression that gcc 8.2 or 8.3, as I've been using, vectorizes Stan at all. That's somewhat subjectively based on the observation that Stan tends to run at least as fast on my Ryzen desktop as on my avx512 desktop, implying no benefit (Julia code can sometimes be 2x faster with avx512).

goedman commented 5 years ago

I have been seeing that warning too on occasion, but don't recall running into subsequent errors. Seems that's because of this commit https://github.com/JuliaLang/julia/commit/1925c211f5609a880f15e046673b2917e63e2b52. I haven't tried this yet, but switching to Arpack master https://github.com/JuliaLinearAlgebra/Arpack.jl/commit/6c1bfec0a7ad0abd4e18e701b595b3f6392fa58c should solve the problem. None of my libraries depend on Arpack as far as I know. Manifest.toml https://github.com/chriselrod/ProbabilityModels.jl/blob/master/Manifest.toml for this project, for example.

Distributions.jl https://github.com/JuliaStats/Distributions.jl/blob/master/REQUIRE requires PDMats.jl https://github.com/JuliaStats/PDMats.jl/blob/master/REQUIRE, which requires Arpack. I'm guessing something, somewhere, is loading Distributions.

Julia 1.3-DEV fails consistently on my system always complaining about ambiguous vmuladd definitions on lines 320 and 329in PaddedMatrices/src/blas.I’ve set Julia 1.3 aside for now. Do you mean you ran the LogisticRegression example with both PM/dHMC and CmdStan, seeing 2-3x better performance with PM/dHMC?

Yes, all in Julia release-1.2, a typical result is [here](https://github.com/StatisticalRethinkingJulia/DynamicHMCModels.jl/blob/master/scripts/ProbabilityModels/LogisticRegressionModel/results.pdf https://github.com/StatisticalRethinkingJulia/DynamicHMCModels.jl/blob/master/scripts/ProbabilityModels/LogisticRegressionModel/results.pdf). I cheat as in the PDF shown I run 4 chains in parallel in CmdStan (each 10000 draws). If I run 1 chain of 40000 draws it takes close to 15 secs. Both results are pretty consistent, I ran about 8 4-chain runs and 4 1-chain runs. The results vary somewhat with the sum(y) result, e.g. 385, 396, 421, etc. And that you ran the ITP model with CmdStan, but not the PM/dHMC version yet?

Yes, initial CmdStan results are [here](https://github.com/StatisticalRethinkingJulia/DynamicHMCModels.jl/blob/master/scripts/ProbabilityModels/ITP/itp_01_results.txt https://github.com/StatisticalRethinkingJulia/DynamicHMCModels.jl/blob/master/scripts/ProbabilityModels/ITP/itp_01_results.txt). I’ve seen 2 timing results with PM/dHMC taking ~600 secs (8 cores, 1TB) and ~1000 secs (newer machine, 4 cores, less memory - .5 TB). I’m trying to do a few more runs as I have seen some inconsistencies, e.g. some runs show AD problems (which in the StatisticalRethinking models I catch with ```∇P = LogDensityRejectErrors(ADgradient(:ForwardDiff, P));

Above timings are without compilation.

When you do, I'd be very interested in learning the results. Unfortunately not all chains converged in both. I'm not sure if I could change that via using different options, or if a more manual intervention is needed. I'm also of course interested in runtime performance.

Let me know if you want to see other performance numbers. Especially if the Stan model runs much faster than mine.

I have never seen CmdStan faster than PM/dHMC. CmdStan and RStan are nearly always identical in my experience (OSX). I have another version of that Stan model that takes about 2 hours to run locally, but I was told by someone else who ran my code with some modifications (using rstan and drake) on a cluster said that all chains finished in < 20 minutes. The cluster has more cores than my computes of course, but I didn't use anything like map_rect, so I don't think the chains should finish any faster for hardware reasons. I did ask for more details.

This suggests to me that different Stan versions, or C++ compilers/versions could make a substantial difference. I don't have the impression that gcc 8.3, as I've been using, vectorizes Stan at all.

Right now I use clang++ with O=3 on OSX 10.14.5 beta with Xcode 10.2. CmdStan on my system uses cmdstan 2.19.1

I use the latest version of CmdStan which converts the chains to MCMCChains.Chains objects by default and save those to do a similar analysis w.r.t. the chains (which converge, which don’t, etc.). But I’m still working on that. Same for the PM/dHMC arrays of NamedTuples, which are quite easily to convert to MCMCChains.Chains. Unfortunately saved .jls files are not portable. And also MCMCChains is not up to par with stansummary w.r.t. ess and Rhat. WIP!

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/chriselrod/ProbabilityModels.jl/issues/2#issuecomment-485182921, or mute the thread https://github.com/notifications/unsubscribe-auth/AAAMX2DIILGK4O2JB7WY5TTPROIGBANCNFSM4HG6WDIQ.

chriselrod commented 5 years ago

Julia 1.3-DEV fails consistently on my system always complaining about ambiguous vmuladd definitions on lines 320 and 329in PaddedMatrices/src/blas.I’ve set Julia 1.3 aside for now.

That should be fixed on master. I commented out one of those definitions.

Yes, initial CmdStan results are here

And also MCMCChains is not up to par with stansummary w.r.t. ess and Rhat.

That's in the same ballpark as it took on my computer. Is there a problem with the ess column of the mcmc chains? Seems to just indicate the total number of iterations. I've been using MCMCDiagnostics for effective sample size calculations. Otherwise, I do like the posterior summaries.

I’ve seen 2 timing results with PM/dHMC taking ~600 secs (8 cores, 1TB) and ~1000 secs (newer machine, 4 cores, less memory - .5 TB).

I'm hoping for better than that. If you're running my code, it should have said addprocs(Sys.CPU_THREADS ÷ 2) instead of hard coding 14.

Or instead run a single chain at a time to track allocations. Doing that... Depending on the tree depths favored by adaptation, on a computer with a (Ryzen) avx2 cpu takes about twice as long. Chains took at most 215 seconds, and that was with 99% reaching a depth of 10. Shallower tree depths reduce the runtime to as little as 78 seconds. I got 11.5 - 33 GB for runs of 2000 iterations on that computer, rather than 0.5-1 TB.

I’m trying to do a few more runs as I have seen some inconsistencies, e.g. some runs show AD problems

What sort of AD problems? Do you mean it is throwing errors, which you wish to convert to -Inf with LogDensityRejectErrors? Or that some of the gradients may be wrong?

Please report any errors or problems you see, and I will do my best to fix them.

I'm also on CmdStan 2.19. My make/local:

O = 3
O_STANC = 3
CXXFLAGS_OS += -O3 -march=native -mprefer-vector-width=512 -funsafe-math-optimizations -ffinite-math-only -fno-math-errno

The extra flags allow gcc to vectorize elementary functions like exp and log. Based on run time, those functions obviously aren't getting vectorized. I may try clang++.

Same for the PM/dHMC arrays of NamedTuples, which are quite easily to convert to MCMCChains.Chains.

ProbabilityModels already has a file named dynamic_hmc_interface.jl, to provide functions that interface with DynamicHMC. So it would probably be convenient for me to provide a function that runs MCMC while also conveniently returning MCMCChains.Chains objects instead of forcing the user to do it manually.

goedman commented 5 years ago

Chris, I’ll try to split my answers up so we can

On Apr 21, 2019, at 07:04, chriselrod notifications@github.com wrote:

Julia 1.3-DEV fails consistently on my system always complaining about ambiguous vmuladd definitions on lines 320 and 329in PaddedMatrices/src/blas.I’ve set Julia 1.3 aside for now.

That should be fixed on master. I commented out https://github.com/chriselrod/PaddedMatrices.jl/blob/master/src/blas.jl#L337 one of those definitions.

Below are results for LinearRegression (single chain, 40000 draws) on Julia 1.3-DEV and CmdStan using the setting you suggested in make.local + clang++.

Conclusions: Julia 1.3 issue solved, CmdStan has gotten substantially faster (~2x) but still ~ 50% slower than PM/dHMC. Here I’m using 6 additional procs on the 8 core machine. No compilation time in these numbers.

Will now switch my other machine over to Julia 1.3 and rebuild cmdstan.

julia> include("/Users/rob/.julia/dev/DynamicHMCModels/scripts/ProbabilityModels/LogisticRegressionModel/lr_pm.jl") Defined model: LogisticRegressionModel. Unknowns: μ1, β1, y, μ0, σ0, β0, X, σ1.

406

ProbabilityModels result:

MCMC, adapting ϵ (75 steps) 7.5e-5 s/step ...done MCMC, adapting ϵ (25 steps) 0.00012 s/step ...done MCMC, adapting ϵ (50 steps) 0.0002 s/step ...done MCMC, adapting ϵ (100 steps) 8.7e-5 s/step ...done MCMC, adapting ϵ (200 steps) 6.3e-5 s/step ...done MCMC, adapting ϵ (400 steps) 7.7e-5 s/step ...done MCMC, adapting ϵ (50 steps) 8.6e-5 s/step ...done MCMC (40000 steps) step 12949 (of 40000), 7.7e-5 s/step step 25968 (of 40000), 7.7e-5 s/step step 38893 (of 40000), 7.7e-5 s/step 7.7e-5 s/step ...done 4.216498 seconds (6.79 M allocations: 1.146 GiB, 3.06% gc time) 2-element Array{ChainDataFrame,1}

Summary Statistics

│ Row │ parameters │ mean │ std │ naive_se │ mcse │ ess │ │ │ Symbol │ Float64 │ Float64 │ Float64 │ Float64 │ Any │ ├─────┼────────────┼────────────┼───────────┼─────────────┼─────────────┼─────────┤ │ 1 │ β0 │ -0.0332652 │ 0.0955164 │ 0.000477582 │ 0.000414905 │ 40000.0 │ │ 2 │ β1[1] │ -1.55186 │ 0.133454 │ 0.000667271 │ 0.000811637 │ 40000.0 │ │ 3 │ β1[2] │ -1.73221 │ 0.138593 │ 0.000692963 │ 0.000855202 │ 40000.0 │ │ 4 │ β1[3] │ -0.0688642 │ 0.0969617 │ 0.000484809 │ 0.000438076 │ 40000.0 │ │ 5 │ β1[4] │ 0.512149 │ 0.102308 │ 0.000511542 │ 0.000489728 │ 40000.0 │

Quantiles

│ Row │ parameters │ 2.5% │ 25.0% │ 50.0% │ 75.0% │ 97.5% │ │ │ Symbol │ Float64 │ Float64 │ Float64 │ Float64 │ Float64 │ ├─────┼────────────┼───────────┼────────────┼────────────┼─────────────┼──────────┤ │ 1 │ β0 │ -0.218864 │ -0.0972999 │ -0.032977 │ 0.0306965 │ 0.154842 │ │ 2 │ β1[1] │ -1.82067 │ -1.64079 │ -1.54893 │ -1.46035 │ -1.29861 │ │ 3 │ β1[2] │ -2.01168 │ -1.82382 │ -1.72919 │ -1.6378 │ -1.46891 │ │ 4 │ β1[3] │ -0.261978 │ -0.13357 │ -0.0684109 │ -0.00351573 │ 0.121338 │ │ 5 │ β1[4] │ 0.314291 │ 0.442826 │ 0.511881 │ 0.580661 │ 0.715105 │

CmdStan result:

7.549346 seconds (4.01 M allocations: 168.327 MiB, 0.53% gc time) 2-element Array{ChainDataFrame,1}

Summary Statistics

│ Row │ parameters │ mean │ std │ naive_se │ mcse │ ess │ │ │ Symbol │ Float64 │ Float64 │ Float64 │ Float64 │ Any │ ├─────┼────────────┼────────────┼───────────┼─────────────┼─────────────┼─────────┤ │ 1 │ beta0 │ -0.0329488 │ 0.0962408 │ 0.000481204 │ 0.000462075 │ 40000.0 │ │ 2 │ beta1.1 │ -1.55049 │ 0.133561 │ 0.000667804 │ 0.000672946 │ 40000.0 │ │ 3 │ beta1.2 │ -1.73123 │ 0.138118 │ 0.000690589 │ 0.000712985 │ 40000.0 │ │ 4 │ beta1.3 │ -0.0689648 │ 0.0964434 │ 0.000482217 │ 0.000462029 │ 40000.0 │ │ 5 │ beta1.4 │ 0.511849 │ 0.10195 │ 0.000509748 │ 0.000484273 │ 40000.0 │

Quantiles

│ Row │ parameters │ 2.5% │ 25.0% │ 50.0% │ 75.0% │ 97.5% │ │ │ Symbol │ Float64 │ Float64 │ Float64 │ Float64 │ Float64 │ ├─────┼────────────┼───────────┼────────────┼────────────┼─────────────┼──────────┤ │ 1 │ beta0 │ -0.222728 │ -0.0975711 │ -0.0330679 │ 0.032645 │ 0.154335 │ │ 2 │ beta1.1 │ -1.82055 │ -1.63825 │ -1.54831 │ -1.45893 │ -1.29423 │ │ 3 │ beta1.2 │ -2.00893 │ -1.82302 │ -1.72947 │ -1.63622 │ -1.46714 │ │ 4 │ beta1.3 │ -0.260844 │ -0.133277 │ -0.0686811 │ -0.00436442 │ 0.119362 │ │ 5 │ beta1.4 │ 0.313982 │ 0.443209 │ 0.511489 │ 0.579326 │ 0.714491 │

goedman commented 5 years ago

And also MCMCChains is not up to par with stansummary w.r.t. ess and Rhat.

Is there a problem with the ess column of the mcmc chains? Seems to just indicate the total number of iterations. I've been using MCMCDiagnostics for effective sample size calculations. Otherwise, I do like the posterior summaries.

Yes, @cpfiffer and I are still working on Rhat and ess in MCMCChains summaries. My preference would be to implement the new proposal (https://arxiv.org/abs/1903.08008 https://arxiv.org/abs/1903.08008) but that will take some time. I’m also adding this request to the informal Stan developer focus survey (https://docs.google.com/forms/d/e/1FAIpQLSfE1C2s9z5woTrKm2Su_kNfp1DtmlNp2Rg3kIqHgiytemWL2w/viewform?usp=send_form https://docs.google.com/forms/d/e/1FAIpQLSfE1C2s9z5woTrKm2Su_kNfp1DtmlNp2Rg3kIqHgiytemWL2w/viewform?usp=send_form).

goedman commented 5 years ago

Below an example of the type of errors I occasionally. Several of the StatisticalRethinking models show errors if the errors are not rejected by encapsulating ADGradient.

What sort of AD problems? Do you mean it is throwing errors, which you wish to convert to -Inf with LogDensityRejectErrors? Or that some of the gradients may be wrong?

julia> include("/Users/rob/.julia/dev/DynamicHMCModels/scripts/ProbabilityModels/ITP/itp_pm.jl") Defined model: ITPModel. Unknowns: Y₂, domains, μh₁, μᵣ₁, μh₂, σ, σᵦ, δh, θ, μᵣ₂, ρ, Y₁, κ, L, σh, βᵣ₂, t, βᵣ₁.

1040.118877 seconds (806.01 k allocations: 99.710 MiB, 0.72% gc time) ERROR: Error while loading expression starting at /Users/rob/.julia/dev/DynamicHMCModels/scripts/ProbabilityModels/ITP/itp_pm.jl:106 caused by [exception 1] On worker 18: ArgumentError: 0 ≤ a must hold. Got 0 => 0.0 a => NaN macro expansion at /Users/rob/.julia/packages/ArgCheck/BUMkA/src/checks.jl:165 [inlined] adapt_stepsize at /Users/rob/.julia/packages/DynamicHMC/gnaxv/src/stepsize.jl:162 tune at /Users/rob/.julia/packages/DynamicHMC/gnaxv/src/sampler.jl:74 macro expansion at /Users/rob/.julia/dev/ProbabilityModels/src/dynamic_hmc_interface.jl:354 [inlined] stable_tune at /Users/rob/.julia/dev/ProbabilityModels/src/dynamic_hmc_interface.jl:348

NUTS_init_tune_mcmc_default#79 at /Users/rob/.julia/dev/ProbabilityModels/src/dynamic_hmc_interface.jl:998

chriselrod commented 5 years ago

That error is within DynamicHMC, not the logdensity. So I don't imagine such a wrapper will help in this case. ProbabilityModel's ValueGradient methods check to make sure the logdensity and gradient isfinite, returning -Inf as the density if either aren't. This behaves fairly similarly to LogDensityRejectErrors.

Although, some of my functions currently assume the arguments are within a range that produces finite answers.

julia> using SLEEFPirates, SLEEF

julia> exp(-150450.0)
0.0

julia> SLEEF.exp(-150450.0)
0.0

julia> SLEEFPirates.exp(-150450.0)
2.4825405905217762e10

Perhaps I ought to revert that change, adding back those checks. [EDIT: For now, I did revert that change for exp and log. They should now give correct results when the answer will be Inf,-Inf, and 0.] It may be the case that before the sampler converges, it may enter some badly behaved regions very far from the mode that lead to the above sort of incorrect behavior, and then this causes failure. The fix is to just compare the inputs with

julia> SLEEFPirates.min_exp(Float64)
-745.1332191019412

julia> SLEEFPirates.max_exp(Float64)
709.7827111495575

However, I am not sure yet that this is the problem. Tracking it: adapt_stepsize fails the check @argcheck 0 ≤ a ≤ 1, because isnan(A) == true. It was called by mcmc_adapting_ϵ. get_acceptance_rate(x::DivergenceStatistic) = x.∑a / x.steps So, if it isnan, either x.∑a and x.steps are 0, or isnan(x.∑a) == true (x.steps isa Int). The divergence statistic is created by a call to sample_trajectory.

This eventually leads to leaf:

function leaf(trajectory::Trajectory, z, isinitial)
    @unpack H, π₀, min_Δ = trajectory
    Δ = isinitial ? zero(π₀) : neg_energy(H, z) - π₀
    isdiv = min_Δ > Δ
    d = isinitial ? divergence_statistic() : divergence_statistic(isdiv, Δ)
    ζ = isdiv ? nothing : Proposal(z, Δ)
    τ = isdiv ? nothing : (p♯ = get_p♯(trajectory.H, z); TurnStatistic(p♯, p♯, z.p))
    ζ, τ, d
end

If neg_energy(H,z) == -Inf and π₀ == -Inf, then...

 Δ = isinitial ? zero(π₀) : neg_energy(H, z) - π₀

Δ = NaN. In this case, isdiv = false ( x > NaN == false ), and divergence_statistic:

divergence_statistic(isdivergent, Δ) =
DivergenceStatistic(isdivergent, Δ ≥ 0 ? one(Δ) : exp(Δ), 1)

NaN > 0 also returns false, therefore it returns: DivergenceStatistic(false, NaN, 1).

I'm not sure if that is what's happening. I'd need a reproducible example and some way to actually dive in and run the code. But this does seem like one possibility that I think will result in the observed error.

chriselrod commented 5 years ago

Glad the 1.3 issues is solved. Your CmdStan runtime is about the same as what I have locally, while your PM/dHMC is 2-4 times slower. I'll try clang to see if it also shrinks the gap on my machines.

The intervals aren't too different, and means within a couple mcse.

Having that new rhat would be great.

chriselrod commented 5 years ago

clang++-8.0 is slightly slower at 11 vs 10 seconds for the y ~ bernoulli_logit(beta0 + X * beta1); version of the model, but faster when using the new y ~ bernoulli_logit_glm(X, beta0, beta1); taking about 4.8 seconds instead of 7.4. I'll try clang for the other models. (EDIT: clang generated slower code [based on "Gradient evaluation took..."] for the ITP model and crashed.) PM/dHMC takes 1.2 seconds on this computer. A little slower than before, I think because I lowered the clock speed slightly because the CPU overheated enough to crash the computer. Although it's possible some of my latest changes made the code slower.

goedman commented 5 years ago

Chris, this is the first of a complete run. This is with 2 CmdStan chains at a time on my 4 core CPU:

Julia Version 1.3.0-DEV.87
Commit 773140e1cb (2019-04-22 18:06 UTC)
Platform Info:
  OS: macOS (x86_64-apple-darwin18.6.0)
  CPU: Intel(R) Core(TM) i5-6267U CPU @ 2.90GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-6.0.1 (ORCJIT, skylake)
Environment:
  JULIA_PKG3_PRECOMPILE = true
  JULIA_CMDSTAN_HOME = /Users/rob/Projects/StanSupport/cmdstan
  JULIA_EDITOR = see
  JULIA_SVG_BROWSER = Google Chrome.app
  JULIA_SPECIALFUNCTIONS_BUILD_SOURCE = true

Results:

               _
   _       _ _(_)_     |  Documentation: https://docs.julialang.org
  (_)     | (_) (_)    |
   _ _   _| |_  __ _   |  Type "?" for help, "]?" for Pkg help.
  | | | | | | |/ _` |  |
  | | |_| | | | (_| |  |  Version 1.3.0-DEV.87 (2019-04-22)
 _/ |\__'_|_|_|\__'_|  |  Commit 773140e1cb (0 days old master)
|__/                   |

julia> include("/Users/rob/.julia/dev/DynamicHMCModels/scripts/ProbabilityModels/ITP/pm_itp.jl")
    Defined model: ITPModel.
    Unknowns: Y₂, domains, μh₁, μᵣ₁, μh₂, σ, σᵦ, δh, θ, μᵣ₂, ρ, Y₁, κ, L, σh, βᵣ₂, t, βᵣ₁.

110.642647 seconds (37.99 M allocations: 1.995 GiB, 1.06% gc time)
 47.384251 seconds (58.66 k allocations: 6.979 MiB)
4×3 Array{Float64,2}:
  755.727   845.057  1758.28
  540.127   278.012  2000.0 
  953.185  1204.47   1991.65
 1028.14    936.552  2000.0 

converged = vec(sum(ess, dims=2)) .> 1000 = Bool[1, 1, 1, 1]

not_converged = .!(converged) = Bool[0, 0, 0, 0]

NUTS_statistics.(chains[not_converged]) = DynamicHMC.NUTS_Statistics{Float64,DataStructures.Accumulator{DynamicHMC.Termination,Int64},DataStructures.Accumulator{Int64,Int64}}[]

NUTS_statistics.(chains[converged]) = DynamicHMC.NUTS_Statistics{Float64,DataStructures.Accumulator{DynamicHMC.Termination,Int64},DataStructures.Accumulator{Int64,Int64}}[Hamiltonian Monte Carlo sample of length 2000
  acceptance rate mean: 0.81, min/25%/median/75%/max: 0.0 0.77 0.92 0.97 1.0
  termination: AdjacentDivergent => 4% AdjacentTurn => 1% DoubledTurn => 95%
  depth: 2 => 0% 3 => 1% 4 => 1% 5 => 4% 6 => 94% 7 => 0% 8 => 0%
, Hamiltonian Monte Carlo sample of length 2000
  acceptance rate mean: 0.86, min/25%/median/75%/max: 0.0 0.85 0.95 0.98 1.0
  termination: AdjacentDivergent => 3% AdjacentTurn => 1% DoubledTurn => 96%
  depth: 2 => 0% 3 => 1% 4 => 1% 5 => 1% 6 => 97%
, Hamiltonian Monte Carlo sample of length 2000
  acceptance rate mean: 0.91, min/25%/median/75%/max: 0.0 0.89 0.96 0.99 1.0
  termination: AdjacentDivergent => 1% AdjacentTurn => 12% DoubledTurn => 88%
  depth: 3 => 0% 4 => 0% 5 => 0% 6 => 96% 7 => 3%
, Hamiltonian Monte Carlo sample of length 2000
  acceptance rate mean: 0.89, min/25%/median/75%/max: 0.0 0.89 0.96 0.99 1.0
  termination: AdjacentDivergent => 1% AdjacentTurn => 1% DoubledTurn => 98%
  depth: 3 => 0% 4 => 0% 5 => 1% 6 => 99% 7 => 0%
]

size(poi_chain) = (3,)
"Major Quantiles for paramter with true values: -3.0:"
2×5 Array{Float64,2}:
  0.05      0.25      0.5       0.75      0.95   
 -4.63938  -3.93392  -3.51467  -3.09541  -2.41137
"Major Quantiles for paramter with true values: 9.0:"
2×5 Array{Float64,2}:
 0.05     0.25     0.5      0.75      0.95  
 7.96502  8.65698  9.07663  9.51046  10.2059
"Major Quantiles for paramter with true values: 0.7:"
2×5 Array{Float64,2}:
 0.05      0.25      0.5       0.75      0.95    
 0.695012  0.698564  0.701065  0.703691  0.707269

File /Users/rob/.julia/dev/DynamicHMCModels/scripts/ProbabilityModels/ITP/tmp/itp1.stan will be updated.

Inference for Stan model: itp1_model
2 chains: each with iter=(2000,2000); warmup=(0,0); thin=(1,1); 4000 iterations saved.

Warmup took (2660, 3371) seconds, 1.7 hours total
Sampling took (10028, 10049) seconds, 5.6 hours total

                    Mean     MCSE   StdDev        5%       50%       95%  N_Eff  N_Eff/s    R_hat
lp__            -1.0e+04  2.0e-01  7.2e+00  -1.0e+04  -1.0e+04  -1.0e+04   1339  6.7e-02  1.0e+00
accept_stat__    9.9e-01  2.6e-04  1.4e-02   9.7e-01   1.0e+00   1.0e+00   2808  1.4e-01  1.0e+00
stepsize__       9.0e-03  5.0e-04  5.0e-04   8.5e-03   9.5e-03   9.5e-03    1.0  5.0e-05  8.7e+13
treedepth__      9.0e+00      nan  8.2e-02   9.0e+00   9.0e+00   9.0e+00    nan      nan  1.0e+00
n_leapfrog__     5.1e+02      nan  1.3e+01   5.1e+02   5.1e+02   5.1e+02    nan      nan  1.0e+00
divergent__      0.0e+00      nan  0.0e+00   0.0e+00   0.0e+00   0.0e+00    nan      nan      nan
energy__         1.0e+04  2.8e-01  1.0e+01   1.0e+04   1.0e+04   1.0e+04   1273  6.3e-02  1.0e+00
muh[1]          -3.5e+00  2.0e-02  7.9e-01  -4.8e+00  -3.5e+00  -2.2e+00   1525  7.6e-02  1.0e+00
muh[2]           9.1e+00  2.2e-02  8.0e-01   7.8e+00   9.1e+00   1.0e+01   1333  6.6e-02  1.0e+00
rho              7.0e-01  8.6e-05  3.8e-03   6.9e-01   7.0e-01   7.1e-01   1997  9.9e-02  1.0e+00
kappa[1]         2.8e-01  1.6e-04  1.1e-02   2.6e-01   2.8e-01   2.9e-01   5086  2.5e-01  1.0e+00
kappa[2]         1.3e-01  8.9e-05  5.7e-03   1.2e-01   1.3e-01   1.4e-01   4053  2.0e-01  1.0e+00
kappa[3]         2.1e-01  9.0e-05  6.6e-03   2.0e-01   2.1e-01   2.2e-01   5319  2.6e-01  1.0e+00
kappa[4]         3.0e-01  1.5e-04  1.1e-02   2.8e-01   3.0e-01   3.2e-01   5307  2.6e-01  1.0e+00
kappa[5]         2.5e-01  1.7e-04  1.2e-02   2.3e-01   2.5e-01   2.7e-01   5067  2.5e-01  1.0e+00
kappa[6]         1.6e-01  1.4e-04  7.8e-03   1.4e-01   1.6e-01   1.7e-01   2918  1.5e-01  1.0e+00
kappa[7]         8.6e-02  8.1e-05  5.4e-03   7.7e-02   8.6e-02   9.5e-02   4460  2.2e-01  1.0e+00
kappa[8]         2.7e-01  1.6e-04  1.2e-02   2.5e-01   2.7e-01   2.9e-01   5557  2.8e-01  1.0e+00
kappa[9]         2.5e-01  1.7e-04  1.3e-02   2.3e-01   2.5e-01   2.7e-01   5835  2.9e-01  1.0e+00
theta[1]        -2.6e+00  1.2e-03  8.6e-02  -2.7e+00  -2.5e+00  -2.4e+00   5445  2.7e-01  1.0e+00
theta[2]         1.2e+00  7.8e-04  6.2e-02   1.1e+00   1.2e+00   1.3e+00   6259  3.1e-01  1.0e+00
theta[3]         1.8e+00  8.8e-04  7.0e-02   1.7e+00   1.8e+00   1.9e+00   6203  3.1e-01  1.0e+00
theta[4]        -7.1e-01  1.1e-03  9.1e-02  -8.6e-01  -7.1e-01  -5.6e-01   6502  3.2e-01  1.0e+00
theta[5]        -3.5e+00  1.2e-03  8.4e-02  -3.6e+00  -3.5e+00  -3.3e+00   5072  2.5e-01  1.0e+00
theta[6]        -1.2e-01  1.0e-03  7.2e-02  -2.4e-01  -1.2e-01   5.3e-04   5082  2.5e-01  1.0e+00
theta[7]        -2.4e-01  1.1e-03  9.0e-02  -4.0e-01  -2.4e-01  -9.6e-02   6181  3.1e-01  1.0e+00
theta[8]        -4.0e+00  9.7e-04  8.1e-02  -4.1e+00  -4.0e+00  -3.8e+00   7036  3.5e-01  1.0e+00
theta[9]        -1.8e+00  1.2e-03  8.9e-02  -2.0e+00  -1.8e+00  -1.7e+00   5406  2.7e-01  1.0e+00
L[1,1]           1.0e+00      nan  6.7e-16   1.0e+00   1.0e+00   1.0e+00    nan      nan  1.0e+00
L[1,2]           0.0e+00      nan  0.0e+00   0.0e+00   0.0e+00   0.0e+00    nan      nan      nan
L[1,3]           0.0e+00      nan  0.0e+00   0.0e+00   0.0e+00   0.0e+00    nan      nan      nan
L[1,4]           0.0e+00      nan  0.0e+00   0.0e+00   0.0e+00   0.0e+00    nan      nan      nan
L[1,5]           0.0e+00      nan  0.0e+00   0.0e+00   0.0e+00   0.0e+00    nan      nan      nan
L[1,6]           0.0e+00      nan  0.0e+00   0.0e+00   0.0e+00   0.0e+00    nan      nan      nan
L[1,7]           0.0e+00      nan  0.0e+00   0.0e+00   0.0e+00   0.0e+00    nan      nan      nan
L[1,8]           0.0e+00      nan  0.0e+00   0.0e+00   0.0e+00   0.0e+00    nan      nan      nan
L[1,9]           0.0e+00      nan  0.0e+00   0.0e+00   0.0e+00   0.0e+00    nan      nan      nan
L[2,1]          -2.1e-01  1.5e-04  1.2e-02  -2.3e-01  -2.1e-01  -1.9e-01   6384  3.2e-01  1.0e+00
L[2,2]           9.8e-01  3.4e-05  2.7e-03   9.7e-01   9.8e-01   9.8e-01   6306  3.1e-01  1.0e+00
L[2,3]           0.0e+00      nan  0.0e+00   0.0e+00   0.0e+00   0.0e+00    nan      nan      nan
L[2,4]           0.0e+00      nan  0.0e+00   0.0e+00   0.0e+00   0.0e+00    nan      nan      nan
L[2,5]           0.0e+00      nan  0.0e+00   0.0e+00   0.0e+00   0.0e+00    nan      nan      nan
L[2,6]           0.0e+00      nan  0.0e+00   0.0e+00   0.0e+00   0.0e+00    nan      nan      nan
L[2,7]           0.0e+00      nan  0.0e+00   0.0e+00   0.0e+00   0.0e+00    nan      nan      nan
L[2,8]           0.0e+00      nan  0.0e+00   0.0e+00   0.0e+00   0.0e+00    nan      nan      nan
L[2,9]           0.0e+00      nan  0.0e+00   0.0e+00   0.0e+00   0.0e+00    nan      nan      nan
L[3,1]           2.8e-02  1.6e-04  1.3e-02   6.8e-03   2.7e-02   4.8e-02   6806  3.4e-01  1.0e+00
L[3,2]          -4.1e-01  1.2e-04  1.1e-02  -4.3e-01  -4.1e-01  -4.0e-01   8820  4.4e-01  1.0e+00
L[3,3]           9.1e-01  5.3e-05  5.0e-03   9.0e-01   9.1e-01   9.2e-01   8809  4.4e-01  1.0e+00
L[3,4]           0.0e+00      nan  0.0e+00   0.0e+00   0.0e+00   0.0e+00    nan      nan      nan
L[3,5]           0.0e+00      nan  0.0e+00   0.0e+00   0.0e+00   0.0e+00    nan      nan      nan
L[3,6]           0.0e+00      nan  0.0e+00   0.0e+00   0.0e+00   0.0e+00    nan      nan      nan
L[3,7]           0.0e+00      nan  0.0e+00   0.0e+00   0.0e+00   0.0e+00    nan      nan      nan
L[3,8]           0.0e+00      nan  0.0e+00   0.0e+00   0.0e+00   0.0e+00    nan      nan      nan
L[3,9]           0.0e+00      nan  0.0e+00   0.0e+00   0.0e+00   0.0e+00    nan      nan      nan
L[4,1]           1.8e-02  1.5e-04  1.3e-02  -4.1e-03   1.8e-02   4.0e-02   7709  3.8e-01  1.0e+00
L[4,2]          -4.1e-02  1.4e-04  1.3e-02  -6.2e-02  -4.1e-02  -2.0e-02   8794  4.4e-01  1.0e+00
L[4,3]          -5.8e-02  1.6e-04  1.3e-02  -7.9e-02  -5.8e-02  -3.6e-02   6784  3.4e-01  1.0e+00
L[4,4]           1.0e+00  1.2e-05  9.6e-04   1.0e+00   1.0e+00   1.0e+00   6611  3.3e-01  1.0e+00
L[4,5]           0.0e+00      nan  0.0e+00   0.0e+00   0.0e+00   0.0e+00    nan      nan      nan
L[4,6]           0.0e+00      nan  0.0e+00   0.0e+00   0.0e+00   0.0e+00    nan      nan      nan
L[4,7]           0.0e+00      nan  0.0e+00   0.0e+00   0.0e+00   0.0e+00    nan      nan      nan
L[4,8]           0.0e+00      nan  0.0e+00   0.0e+00   0.0e+00   0.0e+00    nan      nan      nan
L[4,9]           0.0e+00      nan  0.0e+00   0.0e+00   0.0e+00   0.0e+00    nan      nan      nan
L[5,1]           4.7e-01  1.3e-04  1.0e-02   4.5e-01   4.7e-01   4.9e-01   6637  3.3e-01  1.0e+00
L[5,2]          -1.3e-01  1.3e-04  1.1e-02  -1.5e-01  -1.3e-01  -1.2e-01   7843  3.9e-01  1.0e+00
L[5,3]           3.5e-02  1.3e-04  1.1e-02   1.6e-02   3.5e-02   5.4e-02   7730  3.9e-01  1.0e+00
L[5,4]           6.4e-02  1.2e-04  1.1e-02   4.5e-02   6.4e-02   8.3e-02   8583  4.3e-01  1.0e+00
L[5,5]           8.7e-01  7.0e-05  5.6e-03   8.6e-01   8.7e-01   8.8e-01   6405  3.2e-01  1.0e+00
L[5,6]           0.0e+00      nan  0.0e+00   0.0e+00   0.0e+00   0.0e+00    nan      nan      nan
L[5,7]           0.0e+00      nan  0.0e+00   0.0e+00   0.0e+00   0.0e+00    nan      nan      nan
L[5,8]           0.0e+00      nan  0.0e+00   0.0e+00   0.0e+00   0.0e+00    nan      nan      nan
L[5,9]           0.0e+00      nan  0.0e+00   0.0e+00   0.0e+00   0.0e+00    nan      nan      nan
L[6,1]           2.6e-01  1.5e-04  1.2e-02   2.4e-01   2.6e-01   2.8e-01   6087  3.0e-01  1.0e+00
L[6,2]          -5.6e-02  1.4e-04  1.3e-02  -7.6e-02  -5.6e-02  -3.5e-02   8204  4.1e-01  1.0e+00
L[6,3]           2.1e-01  1.4e-04  1.2e-02   1.9e-01   2.1e-01   2.3e-01   8149  4.1e-01  1.0e+00
L[6,4]           2.1e-01  1.4e-04  1.2e-02   1.9e-01   2.1e-01   2.3e-01   7311  3.6e-01  1.0e+00
L[6,5]           4.6e-02  1.4e-04  1.2e-02   2.5e-02   4.6e-02   6.6e-02   8440  4.2e-01  1.0e+00
L[6,6]           9.2e-01  5.9e-05  4.8e-03   9.1e-01   9.2e-01   9.2e-01   6784  3.4e-01  1.0e+00
L[6,7]           0.0e+00      nan  0.0e+00   0.0e+00   0.0e+00   0.0e+00    nan      nan      nan
L[6,8]           0.0e+00      nan  0.0e+00   0.0e+00   0.0e+00   0.0e+00    nan      nan      nan
L[6,9]           0.0e+00      nan  0.0e+00   0.0e+00   0.0e+00   0.0e+00    nan      nan      nan
L[7,1]           2.5e-01  1.5e-04  1.2e-02   2.3e-01   2.5e-01   2.7e-01   7217  3.6e-01  1.0e+00
L[7,2]           1.7e-01  1.4e-04  1.2e-02   1.5e-01   1.7e-01   1.9e-01   7911  3.9e-01  1.0e+00
L[7,3]           1.5e-02  1.4e-04  1.3e-02  -4.6e-03   1.5e-02   3.7e-02   7788  3.9e-01  1.0e+00
L[7,4]          -9.1e-02  1.4e-04  1.2e-02  -1.1e-01  -9.1e-02  -7.1e-02   7321  3.6e-01  1.0e+00
L[7,5]          -1.4e-01  1.5e-04  1.3e-02  -1.6e-01  -1.4e-01  -1.2e-01   7671  3.8e-01  1.0e+00
L[7,6]          -2.8e-01  1.1e-04  1.1e-02  -3.0e-01  -2.8e-01  -2.6e-01  10436  5.2e-01  1.0e+00
L[7,7]           9.0e-01  6.0e-05  5.4e-03   8.9e-01   9.0e-01   9.0e-01   8060  4.0e-01  1.0e+00
L[7,8]           0.0e+00      nan  0.0e+00   0.0e+00   0.0e+00   0.0e+00    nan      nan      nan
L[7,9]           0.0e+00      nan  0.0e+00   0.0e+00   0.0e+00   0.0e+00    nan      nan      nan
L[8,1]          -1.5e-01  1.5e-04  1.3e-02  -1.7e-01  -1.5e-01  -1.3e-01   8040  4.0e-01  1.0e+00
L[8,2]          -1.3e-01  1.4e-04  1.3e-02  -1.5e-01  -1.3e-01  -1.1e-01   8726  4.3e-01  1.0e+00
L[8,3]          -7.4e-02  1.4e-04  1.3e-02  -9.5e-02  -7.4e-02  -5.3e-02   8621  4.3e-01  1.0e+00
L[8,4]          -7.9e-02  1.4e-04  1.3e-02  -1.0e-01  -7.9e-02  -5.8e-02   8261  4.1e-01  1.0e+00
L[8,5]          -6.4e-02  1.3e-04  1.3e-02  -8.5e-02  -6.4e-02  -4.4e-02   9220  4.6e-01  1.0e+00
L[8,6]           8.6e-02  1.3e-04  1.3e-02   6.5e-02   8.6e-02   1.1e-01   9688  4.8e-01  1.0e+00
L[8,7]          -3.9e-02  1.3e-04  1.3e-02  -6.0e-02  -3.9e-02  -1.7e-02  10696  5.3e-01  1.0e+00
L[8,8]           9.7e-01  3.2e-05  3.1e-03   9.6e-01   9.7e-01   9.7e-01   9312  4.6e-01  1.0e+00
L[8,9]           0.0e+00      nan  0.0e+00   0.0e+00   0.0e+00   0.0e+00    nan      nan      nan
L[9,1]           2.1e-01  1.5e-04  1.3e-02   1.9e-01   2.1e-01   2.3e-01   7268  3.6e-01  1.0e+00
L[9,2]          -6.8e-02  1.5e-04  1.3e-02  -8.9e-02  -6.8e-02  -4.7e-02   7453  3.7e-01  1.0e+00
L[9,3]          -3.3e-02  1.5e-04  1.3e-02  -5.4e-02  -3.3e-02  -1.2e-02   8154  4.1e-01  1.0e+00
L[9,4]           1.0e-01  1.4e-04  1.3e-02   8.0e-02   1.0e-01   1.2e-01   7962  4.0e-01  1.0e+00
L[9,5]          -2.7e-01  1.2e-04  1.2e-02  -2.9e-01  -2.7e-01  -2.5e-01  10067  5.0e-01  1.0e+00
L[9,6]           1.7e-01  1.3e-04  1.2e-02   1.5e-01   1.7e-01   1.9e-01   9177  4.6e-01  1.0e+00
L[9,7]           1.5e-02  1.3e-04  1.2e-02  -5.3e-03   1.4e-02   3.5e-02   9158  4.6e-01  1.0e+00
L[9,8]          -1.5e-01  1.1e-04  1.2e-02  -1.7e-01  -1.5e-01  -1.3e-01  10827  5.4e-01  1.0e+00
L[9,9]           9.0e-01  5.8e-05  5.2e-03   8.9e-01   9.0e-01   9.1e-01   8321  4.1e-01  1.0e+00
muraw[1,1]      -2.9e-01  1.2e-02  5.2e-01  -1.1e+00  -2.8e-01   5.5e-01   1775  8.8e-02  1.0e+00
muraw[1,2]      -1.2e+00  1.6e-02  6.4e-01  -2.3e+00  -1.2e+00  -1.9e-01   1718  8.6e-02  1.0e+00
muraw[1,3]       2.5e-01  1.3e-02  5.2e-01  -6.0e-01   2.6e-01   1.1e+00   1715  8.5e-02  1.0e+00
muraw[1,4]       1.2e+00  1.6e-02  6.5e-01   1.7e-01   1.2e+00   2.3e+00   1647  8.2e-02  1.0e+00
muraw[2,1]      -1.7e-01  1.3e-02  5.2e-01  -1.0e+00  -1.9e-01   7.0e-01   1531  7.6e-02  1.0e+00
muraw[2,2]       1.0e+00  1.5e-02  6.0e-01   4.6e-02   9.8e-01   2.0e+00   1605  8.0e-02  1.0e+00
muraw[2,3]      -1.7e-01  1.4e-02  5.2e-01  -1.0e+00  -1.7e-01   6.8e-01   1450  7.2e-02  1.0e+00
muraw[2,4]      -5.2e-01  1.4e-02  5.4e-01  -1.4e+00  -5.2e-01   3.6e-01   1458  7.3e-02  1.0e+00
betaraw[1,1]     1.4e-01  1.2e-02  9.6e-01  -1.5e+00   1.5e-01   1.7e+00   6700  3.3e-01  1.0e+00
betaraw[1,2]    -1.7e-01  1.1e-02  9.4e-01  -1.7e+00  -1.9e-01   1.4e+00   7645  3.8e-01  1.0e+00
betaraw[1,3]     1.2e-01  1.1e-02  9.1e-01  -1.4e+00   1.1e-01   1.6e+00   6312  3.1e-01  1.0e+00
betaraw[1,4]    -2.4e-01  1.1e-02  9.3e-01  -1.8e+00  -2.3e-01   1.3e+00   7261  3.6e-01  1.0e+00
betaraw[1,5]     2.5e-01  1.1e-02  9.0e-01  -1.2e+00   2.5e-01   1.8e+00   7403  3.7e-01  1.0e+00
betaraw[1,6]    -2.2e-01  1.2e-02  9.0e-01  -1.7e+00  -2.4e-01   1.3e+00   5986  3.0e-01  1.0e+00
betaraw[1,7]    -1.0e-02  1.1e-02  9.5e-01  -1.6e+00  -1.6e-02   1.5e+00   7420  3.7e-01  1.0e+00
betaraw[1,8]    -3.3e-01  1.1e-02  9.1e-01  -1.8e+00  -3.5e-01   1.1e+00   6739  3.4e-01  1.0e+00
betaraw[1,9]     4.6e-01  1.2e-02  9.1e-01  -1.0e+00   4.6e-01   2.0e+00   6164  3.1e-01  1.0e+00
betaraw[2,1]    -6.3e-02  1.2e-02  9.5e-01  -1.6e+00  -7.3e-02   1.5e+00   6217  3.1e-01  1.0e+00
betaraw[2,2]     5.2e-02  1.1e-02  9.6e-01  -1.5e+00   4.8e-02   1.6e+00   7191  3.6e-01  1.0e+00
betaraw[2,3]    -1.8e-01  1.1e-02  9.2e-01  -1.7e+00  -1.9e-01   1.3e+00   6481  3.2e-01  1.0e+00
betaraw[2,4]     2.7e-01  1.2e-02  9.5e-01  -1.3e+00   3.0e-01   1.8e+00   6649  3.3e-01  1.0e+00
betaraw[2,5]    -3.6e-01  1.2e-02  9.5e-01  -1.9e+00  -3.7e-01   1.2e+00   6741  3.4e-01  1.0e+00
betaraw[2,6]     3.5e-01  1.1e-02  9.6e-01  -1.2e+00   3.6e-01   1.9e+00   7107  3.5e-01  1.0e+00
betaraw[2,7]    -1.3e-01  1.1e-02  1.0e+00  -1.8e+00  -1.4e-01   1.5e+00   7675  3.8e-01  1.0e+00
betaraw[2,8]     1.3e-01  1.1e-02  9.3e-01  -1.4e+00   1.4e-01   1.7e+00   7029  3.5e-01  1.0e+00
betaraw[2,9]    -4.1e-02  1.0e-02  9.4e-01  -1.6e+00  -3.0e-02   1.4e+00   8566  4.3e-01  1.0e+00
sigma_beta       1.2e-01  2.7e-03  9.4e-02   8.0e-03   9.7e-02   3.0e-01   1251  6.2e-02  1.0e+00
sigma_h          1.5e+00  1.6e-02  6.1e-01   8.2e-01   1.3e+00   2.6e+00   1493  7.4e-02  1.0e+00
sigma[1]         1.6e+00  3.2e-04  1.8e-02   1.6e+00   1.6e+00   1.6e+00   3127  1.6e-01  1.0e+00
sigma[2]         1.2e+00  2.4e-04  1.4e-02   1.2e+00   1.2e+00   1.2e+00   3195  1.6e-01  1.0e+00
sigma[3]         1.2e+00  2.3e-04  1.3e-02   1.2e+00   1.2e+00   1.2e+00   3208  1.6e-01  1.0e+00
sigma[4]         1.7e+00  3.2e-04  1.9e-02   1.7e+00   1.7e+00   1.7e+00   3531  1.8e-01  1.0e+00
sigma[5]         1.5e+00  3.1e-04  1.7e-02   1.5e+00   1.5e+00   1.5e+00   3065  1.5e-01  1.0e+00
sigma[6]         1.3e+00  2.8e-04  1.5e-02   1.3e+00   1.3e+00   1.4e+00   2792  1.4e-01  1.0e+00
sigma[7]         1.9e+00  3.6e-04  2.1e-02   1.9e+00   1.9e+00   1.9e+00   3409  1.7e-01  1.0e+00
sigma[8]         1.4e+00  2.7e-04  1.6e-02   1.4e+00   1.4e+00   1.5e+00   3602  1.8e-01  1.0e+00
sigma[9]         1.6e+00  3.0e-04  1.8e-02   1.6e+00   1.6e+00   1.6e+00   3404  1.7e-01  1.0e+00

Samples were drawn using hmc with nuts.
For each parameter, N_Eff is a crude measure of effective sample size,
and R_hat is the potential scale reduction factor on split chains (at 
convergence, R_hat=1).

14050.393514 seconds (19.86 M allocations: 1013.451 MiB, 0.01% gc time)

describe(chns_itp1_01)
2-element Array{ChainDataFrame,1}

Summary Statistics

│ Row │ parameters │ mean     │ std        │ naive_se   │ mcse      │ ess    │
│     │ Symbol     │ Float64  │ Float64    │ Float64    │ Float64   │ Any    │
├─────┼────────────┼──────────┼────────────┼────────────┼───────────┼────────┤
│ 1   │ muh.1      │ -3.51205 │ 0.791016   │ 0.0125071  │ 0.0205182 │ 4000.0 │
│ 2   │ muh.2      │ 9.05276  │ 0.795115   │ 0.0125719  │ 0.0237154 │ 4000.0 │
│ 3   │ rho        │ 0.701141 │ 0.00382832 │ 6.05311e-5 │ 7.6547e-5 │ 4000.0 │

Quantiles

│ Row │ parameters │ 2.5%     │ 25.0%    │ 50.0%    │ 75.0%    │ 97.5%    │
│     │ Symbol     │ Float64  │ Float64  │ Float64  │ Float64  │ Float64  │
├─────┼────────────┼──────────┼──────────┼──────────┼──────────┼──────────┤
│ 1   │ muh.1      │ -5.10913 │ -3.95511 │ -3.5256  │ -3.05695 │ -1.97156 │
│ 2   │ muh.2      │ 7.45182  │ 8.59805  │ 9.06974  │ 9.51532  │ 10.6231  │
│ 3   │ rho        │ 0.693702 │ 0.698566 │ 0.701114 │ 0.703697 │ 0.708838 │

Inference for Stan model: itp1_model
2 chains: each with iter=(2000,2000); warmup=(0,0); thin=(1,1); 4000 iterations saved.

Warmup took (2608, 3440) seconds, 1.7 hours total
Sampling took (7812, 4637) seconds, 3.5 hours total

                    Mean     MCSE   StdDev        5%       50%       95%  N_Eff  N_Eff/s    R_hat
lp__            -1.0e+04  2.1e-01  7.3e+00  -1.0e+04  -1.0e+04  -1.0e+04   1224  9.8e-02  1.0e+00
accept_stat__    9.9e-01  4.2e-03  2.3e-02   9.5e-01   9.9e-01   1.0e+00     31  2.5e-03  1.0e+00
stepsize__       1.1e-02  2.0e-03  2.0e-03   8.9e-03   1.3e-02   1.3e-02    1.0  8.0e-05  1.6e+14
treedepth__      8.5e+00  4.9e-01  5.0e-01   8.0e+00   9.0e+00   9.0e+00    1.0  8.2e-05  5.3e+00
n_leapfrog__     4.0e+02  1.2e+02  1.3e+02   2.6e+02   5.1e+02   5.1e+02    1.2  9.4e-05  2.3e+00
divergent__      0.0e+00      nan  0.0e+00   0.0e+00   0.0e+00   0.0e+00    nan      nan      nan
energy__         1.0e+04  2.9e-01  1.0e+01   1.0e+04   1.0e+04   1.0e+04   1181  9.5e-02  1.0e+00
muh[1]          -3.5e+00  2.3e-02  7.8e-01  -4.8e+00  -3.5e+00  -2.3e+00   1161  9.3e-02  1.0e+00
muh[2]           9.0e+00  2.4e-02  8.2e-01   7.7e+00   9.1e+00   1.0e+01   1182  9.5e-02  1.0e+00
rho              7.0e-01  1.0e-04  3.9e-03   6.9e-01   7.0e-01   7.1e-01   1469  1.2e-01  1.0e+00
kappa[1]         2.8e-01  2.0e-04  1.2e-02   2.6e-01   2.8e-01   2.9e-01   3217  2.6e-01  1.0e+00
kappa[2]         1.3e-01  1.1e-04  5.7e-03   1.2e-01   1.3e-01   1.4e-01   2619  2.1e-01  1.0e+00
kappa[3]         2.1e-01  1.1e-04  6.7e-03   2.0e-01   2.1e-01   2.2e-01   3828  3.1e-01  1.0e+00
kappa[4]         3.0e-01  1.8e-04  1.1e-02   2.8e-01   3.0e-01   3.2e-01   3552  2.9e-01  1.0e+00
kappa[5]         2.5e-01  2.1e-04  1.2e-02   2.3e-01   2.5e-01   2.7e-01   3325  2.7e-01  1.0e+00
kappa[6]         1.6e-01  1.7e-04  7.7e-03   1.4e-01   1.6e-01   1.7e-01   2157  1.7e-01  1.0e+00
kappa[7]         8.6e-02  9.9e-05  5.3e-03   7.8e-02   8.6e-02   9.5e-02   2870  2.3e-01  1.0e+00
kappa[8]         2.7e-01  1.8e-04  1.2e-02   2.5e-01   2.7e-01   2.9e-01   4553  3.7e-01  1.0e+00
kappa[9]         2.5e-01  1.8e-04  1.2e-02   2.3e-01   2.4e-01   2.7e-01   4665  3.7e-01  1.0e+00
theta[1]        -2.6e+00  1.5e-03  8.6e-02  -2.7e+00  -2.6e+00  -2.4e+00   3145  2.5e-01  1.0e+00
theta[2]         1.2e+00  9.7e-04  6.2e-02   1.1e+00   1.2e+00   1.3e+00   4045  3.2e-01  1.0e+00
theta[3]         1.8e+00  1.1e-03  6.6e-02   1.7e+00   1.8e+00   1.9e+00   3497  2.8e-01  1.0e+00
theta[4]        -7.1e-01  1.6e-03  9.1e-02  -8.6e-01  -7.1e-01  -5.7e-01   3406  2.7e-01  1.0e+00
theta[5]        -3.5e+00  1.5e-03  8.5e-02  -3.6e+00  -3.5e+00  -3.3e+00   3215  2.6e-01  1.0e+00
theta[6]        -1.2e-01  1.3e-03  7.3e-02  -2.4e-01  -1.2e-01   2.4e-03   3280  2.6e-01  1.0e+00
theta[7]        -2.5e-01  1.4e-03  8.7e-02  -3.9e-01  -2.5e-01  -1.1e-01   3704  3.0e-01  1.0e+00
theta[8]        -4.0e+00  1.3e-03  7.8e-02  -4.1e+00  -4.0e+00  -3.8e+00   3551  2.9e-01  1.0e+00
theta[9]        -1.9e+00  1.5e-03  8.9e-02  -2.0e+00  -1.9e+00  -1.7e+00   3627  2.9e-01  1.0e+00
L[1,1]           1.0e+00      nan  6.7e-16   1.0e+00   1.0e+00   1.0e+00    nan      nan  1.0e+00
L[1,2]           0.0e+00      nan  0.0e+00   0.0e+00   0.0e+00   0.0e+00    nan      nan      nan
L[1,3]           0.0e+00      nan  0.0e+00   0.0e+00   0.0e+00   0.0e+00    nan      nan      nan
L[1,4]           0.0e+00      nan  0.0e+00   0.0e+00   0.0e+00   0.0e+00    nan      nan      nan
L[1,5]           0.0e+00      nan  0.0e+00   0.0e+00   0.0e+00   0.0e+00    nan      nan      nan
L[1,6]           0.0e+00      nan  0.0e+00   0.0e+00   0.0e+00   0.0e+00    nan      nan      nan
L[1,7]           0.0e+00      nan  0.0e+00   0.0e+00   0.0e+00   0.0e+00    nan      nan      nan
L[1,8]           0.0e+00      nan  0.0e+00   0.0e+00   0.0e+00   0.0e+00    nan      nan      nan
L[1,9]           0.0e+00      nan  0.0e+00   0.0e+00   0.0e+00   0.0e+00    nan      nan      nan
L[2,1]          -2.1e-01  2.0e-04  1.2e-02  -2.3e-01  -2.1e-01  -1.9e-01   3927  3.2e-01  1.0e+00
L[2,2]           9.8e-01  4.4e-05  2.7e-03   9.7e-01   9.8e-01   9.8e-01   3925  3.2e-01  1.0e+00
L[2,3]           0.0e+00      nan  0.0e+00   0.0e+00   0.0e+00   0.0e+00    nan      nan      nan
L[2,4]           0.0e+00      nan  0.0e+00   0.0e+00   0.0e+00   0.0e+00    nan      nan      nan
L[2,5]           0.0e+00      nan  0.0e+00   0.0e+00   0.0e+00   0.0e+00    nan      nan      nan
L[2,6]           0.0e+00      nan  0.0e+00   0.0e+00   0.0e+00   0.0e+00    nan      nan      nan
L[2,7]           0.0e+00      nan  0.0e+00   0.0e+00   0.0e+00   0.0e+00    nan      nan      nan
L[2,8]           0.0e+00      nan  0.0e+00   0.0e+00   0.0e+00   0.0e+00    nan      nan      nan
L[2,9]           0.0e+00      nan  0.0e+00   0.0e+00   0.0e+00   0.0e+00    nan      nan      nan
L[3,1]           2.7e-02  2.0e-04  1.3e-02   5.5e-03   2.7e-02   4.9e-02   4400  3.5e-01  1.0e+00
L[3,2]          -4.1e-01  1.5e-04  1.1e-02  -4.3e-01  -4.1e-01  -4.0e-01   5484  4.4e-01  1.0e+00
L[3,3]           9.1e-01  6.7e-05  5.0e-03   9.0e-01   9.1e-01   9.2e-01   5474  4.4e-01  1.0e+00
L[3,4]           0.0e+00      nan  0.0e+00   0.0e+00   0.0e+00   0.0e+00    nan      nan      nan
L[3,5]           0.0e+00      nan  0.0e+00   0.0e+00   0.0e+00   0.0e+00    nan      nan      nan
L[3,6]           0.0e+00      nan  0.0e+00   0.0e+00   0.0e+00   0.0e+00    nan      nan      nan
L[3,7]           0.0e+00      nan  0.0e+00   0.0e+00   0.0e+00   0.0e+00    nan      nan      nan
L[3,8]           0.0e+00      nan  0.0e+00   0.0e+00   0.0e+00   0.0e+00    nan      nan      nan
L[3,9]           0.0e+00      nan  0.0e+00   0.0e+00   0.0e+00   0.0e+00    nan      nan      nan
L[4,1]           1.8e-02  1.7e-04  1.3e-02  -3.3e-03   1.8e-02   4.0e-02   6066  4.9e-01  1.0e+00
L[4,2]          -4.1e-02  1.7e-04  1.3e-02  -6.2e-02  -4.1e-02  -1.9e-02   5817  4.7e-01  1.0e+00
L[4,3]          -5.8e-02  1.8e-04  1.3e-02  -8.0e-02  -5.8e-02  -3.6e-02   5529  4.4e-01  1.0e+00
L[4,4]           1.0e+00  1.3e-05  9.9e-04   1.0e+00   1.0e+00   1.0e+00   5723  4.6e-01  1.0e+00
L[4,5]           0.0e+00      nan  0.0e+00   0.0e+00   0.0e+00   0.0e+00    nan      nan      nan
L[4,6]           0.0e+00      nan  0.0e+00   0.0e+00   0.0e+00   0.0e+00    nan      nan      nan
L[4,7]           0.0e+00      nan  0.0e+00   0.0e+00   0.0e+00   0.0e+00    nan      nan      nan
L[4,8]           0.0e+00      nan  0.0e+00   0.0e+00   0.0e+00   0.0e+00    nan      nan      nan
L[4,9]           0.0e+00      nan  0.0e+00   0.0e+00   0.0e+00   0.0e+00    nan      nan      nan
L[5,1]           4.7e-01  1.5e-04  1.0e-02   4.5e-01   4.7e-01   4.9e-01   4774  3.8e-01  1.0e+00
L[5,2]          -1.3e-01  1.4e-04  1.1e-02  -1.5e-01  -1.3e-01  -1.2e-01   6464  5.2e-01  1.0e+00
L[5,3]           3.5e-02  1.5e-04  1.2e-02   1.6e-02   3.5e-02   5.4e-02   6169  5.0e-01  1.0e+00
L[5,4]           6.4e-02  1.6e-04  1.1e-02   4.5e-02   6.4e-02   8.3e-02   5129  4.1e-01  1.0e+00
L[5,5]           8.7e-01  8.2e-05  5.7e-03   8.6e-01   8.7e-01   8.8e-01   4826  3.9e-01  1.0e+00
L[5,6]           0.0e+00      nan  0.0e+00   0.0e+00   0.0e+00   0.0e+00    nan      nan      nan
L[5,7]           0.0e+00      nan  0.0e+00   0.0e+00   0.0e+00   0.0e+00    nan      nan      nan
L[5,8]           0.0e+00      nan  0.0e+00   0.0e+00   0.0e+00   0.0e+00    nan      nan      nan
L[5,9]           0.0e+00      nan  0.0e+00   0.0e+00   0.0e+00   0.0e+00    nan      nan      nan
L[6,1]           2.6e-01  1.8e-04  1.2e-02   2.4e-01   2.6e-01   2.8e-01   4641  3.7e-01  1.0e+00
L[6,2]          -5.5e-02  1.7e-04  1.3e-02  -7.6e-02  -5.5e-02  -3.5e-02   5357  4.3e-01  1.0e+00
L[6,3]           2.1e-01  1.7e-04  1.2e-02   1.9e-01   2.1e-01   2.3e-01   5183  4.2e-01  1.0e+00
L[6,4]           2.1e-01  1.7e-04  1.2e-02   1.9e-01   2.1e-01   2.3e-01   5359  4.3e-01  1.0e+00
L[6,5]           4.6e-02  1.5e-04  1.2e-02   2.6e-02   4.6e-02   6.7e-02   6901  5.5e-01  1.0e+00
L[6,6]           9.2e-01  6.8e-05  4.8e-03   9.1e-01   9.2e-01   9.2e-01   5024  4.0e-01  1.0e+00
L[6,7]           0.0e+00      nan  0.0e+00   0.0e+00   0.0e+00   0.0e+00    nan      nan      nan
L[6,8]           0.0e+00      nan  0.0e+00   0.0e+00   0.0e+00   0.0e+00    nan      nan      nan
L[6,9]           0.0e+00      nan  0.0e+00   0.0e+00   0.0e+00   0.0e+00    nan      nan      nan
L[7,1]           2.5e-01  1.7e-04  1.2e-02   2.3e-01   2.5e-01   2.7e-01   5061  4.1e-01  1.0e+00
L[7,2]           1.7e-01  1.7e-04  1.3e-02   1.5e-01   1.7e-01   1.9e-01   5315  4.3e-01  1.0e+00
L[7,3]           1.5e-02  1.7e-04  1.3e-02  -5.5e-03   1.5e-02   3.6e-02   5381  4.3e-01  1.0e+00
L[7,4]          -9.1e-02  1.6e-04  1.3e-02  -1.1e-01  -9.1e-02  -7.0e-02   6041  4.9e-01  1.0e+00
L[7,5]          -1.4e-01  1.5e-04  1.2e-02  -1.6e-01  -1.4e-01  -1.2e-01   6891  5.5e-01  1.0e+00
L[7,6]          -2.8e-01  1.5e-04  1.1e-02  -3.0e-01  -2.8e-01  -2.6e-01   5986  4.8e-01  1.0e+00
L[7,7]           9.0e-01  7.5e-05  5.4e-03   8.9e-01   9.0e-01   9.0e-01   5168  4.2e-01  1.0e+00
L[7,8]           0.0e+00      nan  0.0e+00   0.0e+00   0.0e+00   0.0e+00    nan      nan      nan
L[7,9]           0.0e+00      nan  0.0e+00   0.0e+00   0.0e+00   0.0e+00    nan      nan      nan
L[8,1]          -1.5e-01  1.7e-04  1.3e-02  -1.7e-01  -1.5e-01  -1.3e-01   5520  4.4e-01  1.0e+00
L[8,2]          -1.3e-01  1.7e-04  1.3e-02  -1.5e-01  -1.3e-01  -1.1e-01   5815  4.7e-01  1.0e+00
L[8,3]          -7.4e-02  1.8e-04  1.3e-02  -9.5e-02  -7.4e-02  -5.3e-02   5468  4.4e-01  1.0e+00
L[8,4]          -7.9e-02  1.6e-04  1.3e-02  -1.0e-01  -7.9e-02  -5.7e-02   6807  5.5e-01  1.0e+00
L[8,5]          -6.5e-02  1.8e-04  1.2e-02  -8.5e-02  -6.4e-02  -4.5e-02   4662  3.7e-01  1.0e+00
L[8,6]           8.6e-02  1.6e-04  1.2e-02   6.6e-02   8.6e-02   1.1e-01   6170  5.0e-01  1.0e+00
L[8,7]          -3.9e-02  1.8e-04  1.3e-02  -5.9e-02  -3.8e-02  -1.8e-02   4997  4.0e-01  1.0e+00
L[8,8]           9.7e-01  4.0e-05  3.2e-03   9.6e-01   9.7e-01   9.7e-01   6620  5.3e-01  1.0e+00
L[8,9]           0.0e+00      nan  0.0e+00   0.0e+00   0.0e+00   0.0e+00    nan      nan      nan
L[9,1]           2.1e-01  1.7e-04  1.3e-02   1.9e-01   2.1e-01   2.3e-01   5547  4.5e-01  1.0e+00
L[9,2]          -6.8e-02  1.7e-04  1.3e-02  -8.9e-02  -6.8e-02  -4.7e-02   5796  4.7e-01  1.0e+00
L[9,3]          -3.3e-02  1.9e-04  1.3e-02  -5.4e-02  -3.3e-02  -1.2e-02   4575  3.7e-01  1.0e+00
L[9,4]           1.0e-01  1.7e-04  1.3e-02   8.0e-02   1.0e-01   1.2e-01   5817  4.7e-01  1.0e+00
L[9,5]          -2.7e-01  1.7e-04  1.2e-02  -2.9e-01  -2.7e-01  -2.5e-01   5111  4.1e-01  1.0e+00
L[9,6]           1.7e-01  1.5e-04  1.2e-02   1.5e-01   1.7e-01   1.9e-01   6553  5.3e-01  1.0e+00
L[9,7]           1.5e-02  1.4e-04  1.2e-02  -4.7e-03   1.5e-02   3.4e-02   7368  5.9e-01  1.0e+00
L[9,8]          -1.5e-01  1.5e-04  1.2e-02  -1.7e-01  -1.5e-01  -1.3e-01   6371  5.1e-01  1.0e+00
L[9,9]           9.0e-01  7.1e-05  5.1e-03   8.9e-01   9.0e-01   9.1e-01   5142  4.1e-01  1.0e+00
muraw[1,1]      -2.8e-01  1.4e-02  5.1e-01  -1.1e+00  -2.8e-01   5.6e-01   1364  1.1e-01  1.0e+00
muraw[1,2]      -1.2e+00  1.7e-02  6.3e-01  -2.3e+00  -1.2e+00  -2.1e-01   1412  1.1e-01  1.0e+00
muraw[1,3]       2.5e-01  1.4e-02  5.1e-01  -5.6e-01   2.5e-01   1.1e+00   1335  1.1e-01  1.0e+00
muraw[1,4]       1.2e+00  1.8e-02  6.5e-01   1.7e-01   1.2e+00   2.3e+00   1305  1.0e-01  1.0e+00
muraw[2,1]      -1.6e-01  1.4e-02  5.2e-01  -9.9e-01  -1.7e-01   7.1e-01   1317  1.1e-01  1.0e+00
muraw[2,2]       1.0e+00  1.6e-02  6.0e-01   2.0e-02   1.0e+00   2.0e+00   1430  1.1e-01  1.0e+00
muraw[2,3]      -1.6e-01  1.5e-02  5.3e-01  -1.0e+00  -1.6e-01   7.2e-01   1297  1.0e-01  1.0e+00
muraw[2,4]      -5.1e-01  1.5e-02  5.4e-01  -1.4e+00  -5.1e-01   4.1e-01   1269  1.0e-01  1.0e+00
betaraw[1,1]     1.5e-01  1.4e-02  9.1e-01  -1.4e+00   1.5e-01   1.6e+00   4066  3.3e-01  1.0e+00
betaraw[1,2]    -1.6e-01  1.4e-02  9.3e-01  -1.7e+00  -1.7e-01   1.4e+00   4218  3.4e-01  1.0e+00
betaraw[1,3]     1.3e-01  1.4e-02  9.6e-01  -1.4e+00   1.4e-01   1.7e+00   4675  3.8e-01  1.0e+00
betaraw[1,4]    -2.3e-01  1.4e-02  9.4e-01  -1.7e+00  -2.5e-01   1.3e+00   4796  3.9e-01  1.0e+00
betaraw[1,5]     2.7e-01  1.4e-02  9.4e-01  -1.3e+00   2.7e-01   1.8e+00   4832  3.9e-01  1.0e+00
betaraw[1,6]    -2.2e-01  1.3e-02  9.1e-01  -1.7e+00  -2.2e-01   1.3e+00   4769  3.8e-01  1.0e+00
betaraw[1,7]     3.3e-02  1.3e-02  9.4e-01  -1.5e+00   4.4e-02   1.6e+00   5228  4.2e-01  1.0e+00
betaraw[1,8]    -3.3e-01  1.4e-02  9.1e-01  -1.8e+00  -3.4e-01   1.2e+00   4294  3.4e-01  1.0e+00
betaraw[1,9]     4.4e-01  1.4e-02  9.1e-01  -1.1e+00   4.6e-01   1.9e+00   4187  3.4e-01  1.0e+00
betaraw[2,1]    -7.3e-02  1.5e-02  9.7e-01  -1.7e+00  -5.9e-02   1.5e+00   4215  3.4e-01  1.0e+00
betaraw[2,2]     3.9e-02  1.6e-02  9.7e-01  -1.5e+00   4.1e-02   1.6e+00   3742  3.0e-01  1.0e+00
betaraw[2,3]    -1.7e-01  1.3e-02  9.3e-01  -1.7e+00  -1.6e-01   1.4e+00   5178  4.2e-01  1.0e+00
betaraw[2,4]     2.4e-01  1.5e-02  9.5e-01  -1.3e+00   2.4e-01   1.8e+00   3944  3.2e-01  1.0e+00
betaraw[2,5]    -3.6e-01  1.4e-02  9.6e-01  -1.8e+00  -3.9e-01   1.2e+00   4786  3.8e-01  1.0e+00
betaraw[2,6]     3.6e-01  1.5e-02  9.8e-01  -1.3e+00   3.8e-01   1.9e+00   4514  3.6e-01  1.0e+00
betaraw[2,7]    -1.3e-01  1.5e-02  1.0e+00  -1.8e+00  -1.4e-01   1.5e+00   4399  3.5e-01  1.0e+00
betaraw[2,8]     1.3e-01  1.3e-02  9.6e-01  -1.5e+00   1.4e-01   1.7e+00   5216  4.2e-01  1.0e+00
betaraw[2,9]    -6.3e-02  1.3e-02  9.3e-01  -1.6e+00  -6.5e-02   1.4e+00   5454  4.4e-01  1.0e+00
sigma_beta       1.2e-01  2.7e-03  9.0e-02   1.1e-02   9.6e-02   3.0e-01   1066  8.6e-02  1.0e+00
sigma_h          1.5e+00  1.8e-02  6.0e-01   8.4e-01   1.3e+00   2.6e+00   1133  9.1e-02  1.0e+00
sigma[1]         1.6e+00  3.9e-04  1.8e-02   1.6e+00   1.6e+00   1.6e+00   2090  1.7e-01  1.0e+00
sigma[2]         1.2e+00  2.9e-04  1.3e-02   1.2e+00   1.2e+00   1.2e+00   2172  1.7e-01  1.0e+00
sigma[3]         1.2e+00  2.8e-04  1.4e-02   1.2e+00   1.2e+00   1.2e+00   2323  1.9e-01  1.0e+00
sigma[4]         1.7e+00  3.8e-04  1.9e-02   1.7e+00   1.7e+00   1.7e+00   2602  2.1e-01  1.0e+00
sigma[5]         1.5e+00  3.6e-04  1.7e-02   1.5e+00   1.5e+00   1.5e+00   2257  1.8e-01  1.0e+00
sigma[6]         1.3e+00  3.1e-04  1.5e-02   1.3e+00   1.3e+00   1.4e+00   2362  1.9e-01  1.0e+00
sigma[7]         1.9e+00  4.0e-04  2.1e-02   1.9e+00   1.9e+00   1.9e+00   2665  2.1e-01  1.0e+00
sigma[8]         1.4e+00  3.2e-04  1.6e-02   1.4e+00   1.4e+00   1.5e+00   2378  1.9e-01  1.0e+00
sigma[9]         1.6e+00  3.8e-04  1.8e-02   1.6e+00   1.6e+00   1.6e+00   2376  1.9e-01  1.0e+00

Samples were drawn using hmc with nuts.
For each parameter, N_Eff is a crude measure of effective sample size,
and R_hat is the potential scale reduction factor on split chains (at 
convergence, R_hat=1).

10810.826194 seconds (4.30 M allocations: 200.594 MiB, 0.01% gc time)

julia> describe(chns_itp2_01)
2-element Array{ChainDataFrame,1}

Summary Statistics

│ Row │ parameters │ mean     │ std        │ naive_se   │ mcse       │ ess    │
│     │ Symbol     │ Float64  │ Float64    │ Float64    │ Float64    │ Any    │
├─────┼────────────┼──────────┼────────────┼────────────┼────────────┼────────┤
│ 1   │ muh.1      │ -3.52398 │ 0.782709   │ 0.0123757  │ 0.0216147  │ 4000.0 │
│ 2   │ muh.2      │ 9.04412  │ 0.816613   │ 0.0129118  │ 0.0244462  │ 4000.0 │
│ 3   │ rho        │ 0.701125 │ 0.00389867 │ 6.16434e-5 │ 9.65721e-5 │ 4000.0 │

Quantiles

│ Row │ parameters │ 2.5%     │ 25.0%    │ 50.0%    │ 75.0%    │ 97.5%    │
│     │ Symbol     │ Float64  │ Float64  │ Float64  │ Float64  │ Float64  │
├─────┼────────────┼──────────┼──────────┼──────────┼──────────┼──────────┤
│ 1   │ muh.1      │ -5.0916  │ -3.99384 │ -3.53361 │ -3.05392 │ -1.92731 │
│ 2   │ muh.2      │ 7.36243  │ 8.58357  │ 9.05681  │ 9.50098  │ 10.6922  │
│ 3   │ rho        │ 0.693551 │ 0.698563 │ 0.701091 │ 0.703774 │ 0.708749 │

A 2nd result will follow as soon as the run is completed. Below 2 plots with the 3 poi parms. Earlier I had a run clearly showing not convergent chains (1st plot below). .

chns_itp1_01.pdf

chns_itp2_01.pdf

chns_itp1_01.pdf

goedman commented 5 years ago

Chris, these are the promised results 2nd run on the 8-core machine.

The top part of below file 'results2.txt' is the slightly shortened 'normal' output. At the end of that file, after selecting and merging all valid chains with the script attached at the bottom, the MCMCChains results look like:

results2.txt

These are the plots for the remaining 6 chains:

itp.pdf

Let me know if this is enough info.

The 'merge' script I used, not very different from yours:

using CmdStan, StatsPlots

ProjDir = @__DIR__

stansummary = CMDSTAN_HOME*"/bin/stansummary"
set_tuple = (
  tmp1 = (name = "itp1", chains = 1:4, num_samples=2000), 
  tmp2 = (name = "itp1", chains = 1:4, num_samples=2000) 
);

# Loop over all paths where the samples were saved.
for key in keys(set_tuple)
  theset = set_tuple[key]
  resdir = joinpath(ProjDir, String(key))
  m = Stanmodel(name = theset.name,
    num_samples=theset.num_samples)

  cd(resdir) do
    #run(`$stansummary $(m.name)_samples_$[i for i ∈ theset.chains].csv`)
  end
end

# Valid chains can be obtained from above stansummary runs
set_tuple = (
  tmp1 = (name = "itp1", chains = [1, 3, 4], num_samples=2000), 
  tmp2 = (name = "itp1", chains = 1:3, num_samples=2000) 
);

# Loop over all valid sets of samples
chns_array = Vector(undef, length(set_tuple))
for (i, key) in enumerate(keys(set_tuple))
  theset = set_tuple[key]
  # Where is teset saved?
  resdir = joinpath(ProjDir, String(key))
  # Dummy Stanmodel, name and num_samples need to be correct
  m = Stanmodel(name = theset.name,
    num_samples=theset.num_samples)

  cd(resdir) do
    run(`$stansummary $(m.name)_samples_$[i for i ∈ theset.chains].csv`)
    a3d, cnames = CmdStan.read_samples(m)
    chns_array[i] = Chains(a3d[:, :, theset.chains], cnames)    
  end
end

chns = chainscat(chns_array...)

chns_itp = set_section(chns, Dict(
  :parameters => ["muh.1", "muh.2", "rho"],
  :sigma_beta => ["sigma_beta", "sigma_h"],
  :L => reshape(["L.$i.$j" for i in 1:9, j in 1:9], 81),
  :betaraw => reshape(["betaraw.$i.$j" for i in 1:2, j in 1:9], 18),
  :kappa => ["kappa.$i" for i in 1:9],
  :sigma => ["sigma.$i" for i in 1:9],
  :theta => ["theta.$i" for i in 1:9],
  :muraw => reshape(["muraw.$i.$j" for i in 1:2, j in 1:4], 8),
  :internals => ["lp__", "accept_stat__", "stepsize__", "treedepth__", "n_leapfrog__",
    "divergent__", "energy__"]
  )
)

write("itp_sections.jls", chns_itp)
pfig = plot(chns_itp)
savefig(pfig, joinpath(ProjDir, "itp.pdf"))

show(chns_itp)

describe(chns_itp, sections=[:internals])