Julia quits unexpectedly when compiling Stan program in Stan.jl

[maxwell]

Hi folks,

I decided to get some practice in Stan.jl by fitting a model I previously fitted in Rstan, without problem. However, Julia just quits unexpectedly when compiling the code. I have been scratching my head trying to find out why. 

I am attaching the code an the data set. Any help would be greatly appreciated it.

Thanks,

[maxwell]

Oops, I am also attaching dat.csv. 

[Bob Carpenter]

I'm not sure how much Rob Goedman monitors the stan-users list.
You might want to raise this as an issue on the Stan.jl issue
tracker. I don't know that anyone else has any Julia or Stan.jl
experience.

- Bob

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

[Daniel Lee]

what version of julia / stan.jl / cmdstan are you using?

[Rob J. Goedman]

HI Maxwell,

I’ll have a look today. Thanks for bringing this up.

Are you using Julia 0.3.6?

Regards,

Rob J. Goedman
goe...@mac.com

[Rob J. Goedman]

Maxwell,

Do you at any point in time see an error message like:

Exception: mismatch in dimension declared and found in context; processing stage=data initialization; variable name=Qj; position=0; dims declared=(23,4); dims found=(22,4)

Diagnostic information: 

Dynamic exception type: St13runtime_error

std::exception::what: mismatch in dimension declared and found in context; processing stage=data initialization; variable name=Qj; position=0; dims declared=(23,4); dims found=(22,4)

It should certainly not crash Julia (and doesn’t on my system, it looks you are on WIndows though). 

Above message to me looks like maybe input lines are being truncated. I’m looking into that possibility.

Please feel free to use my direct email address below, once we have figured out what the problem is we can report back to the list and/or fix Stan.jl.

Regards,

Rob J. Goedman
goe...@mac.com

[Rob J. Goedman]

Maxwell,

Fixed the dimension issue, now Stan runs, but all 4 chains only produce messages (a boatload of them) like:

Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:

INDEX OPERATOR [] OUT OF BOUNDS; index=4; lower bound=1; upper bound=3; index position=2; bijd

If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,

but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.

Rejecting proposed initial value with zero density.

[maxwell]

Hi Rob,

I run my code again correcting Qj and now it is able to translate, link but does not seem to be able to simulate. I get the following error:

ErrorException("failed processes:\n  Process(`cmd /c fmimicme.exe sample num_sam

ples=1000 num_warmup=1000 save_warmup=0 thin=1 adapt engaged=1 gamma=0.05 delta=

0.8 kappa=0.75 t0=10.0 init_buffer=75 term_buffer=50 window=25 algorithm=hmc eng

ine=nuts max_depth=10 metric=diag_e stepsize=1.0 stepsize_jitter=1.0 random seed

=-1 init=2 id=1 data file=fmimicme_1.data.R output file=fmimicme_samples_1.csv r

efresh=100`, ProcessExited(70)) [70]\n  Process(`cmd /c fmimicme.exe sample num_

I am using Julia 0.3.3

julia> versioninfo()

Julia Version 0.3.3

Commit b24213b* (2014-11-23 20:19 UTC)

Platform Info:

  System: Windows (x86_64-w64-mingw32)

  CPU: Intel(R) Core(TM) i5-4570 CPU @ 3.20GHz

  WORD_SIZE: 64

  BLAS: libopenblas (USE64BITINT DYNAMIC_ARCH NO_AFFINITY Haswell)

  LAPACK: libopenblas

  LIBM: libopenlibm

  LLVM: libLLVM-3.3

Thank you,

[Rob J. Goedman]

Hi Maxwell,

Yes, that’s correct because I expect Stan returns an error code (and I should do a better job in Stan.jl to handle that!).

In the tmp subdirectory you should find fmimicme_run.log and several fmimicme_samples_[1-4].csv files showing the ‘fixed’ output of the 4 Stan chains as shown below but no simulation result lines. The xxx_run.log should also contain the error messages I referred to and included below.

Would you mind checking that?

Could it be that this is caused by the simulated input instead your real data?

Regards,

Rob J. Goedman
goe...@mac.com

# stan_version_major = 2

# stan_version_minor = 6

# stan_version_patch = 0

# model = fmimicme_model

# method = sample (Default)

#   sample

#     num_samples = 1000 (Default)

#     num_warmup = 1000 (Default)

#     save_warmup = 0 (Default)

#     thin = 1 (Default)

#     adapt

#       engaged = 1 (Default)

#       gamma = 0.050000000000000003 (Default)

#       delta = 0.80000000000000004 (Default)

#       kappa = 0.75 (Default)

#       t0 = 10 (Default)

#       init_buffer = 75 (Default)

#       term_buffer = 50 (Default)

#       window = 25 (Default)

#     algorithm = hmc (Default)

#       hmc

#         engine = nuts (Default)

#           nuts

#             max_depth = 10 (Default)

#         metric = diag_e (Default)

#         stepsize = 1 (Default)

#         stepsize_jitter = 1

# id = 1

# data

#   file = fmimicme_1.data.R

# init = 2 (Default)

# random

#   seed = 798457501

# output

#   file = fmimicme_samples_1.csv

#   diagnostic_file =  (Default)

#   refresh = 100 (Default)

And in fmimicme_run.log many additional lines:

Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:

INDEX OPERATOR [] OUT OF BOUNDS; index=4; lower bound=1; upper bound=3; index position=2; bijd

If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,

but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.

Rejecting proposed initial value with zero density.

….

[Roger Zoh]

Rob,

I found the problem once I open the xxx_run.log: the second column of dat had some issues. The values were off by 1. The code seems running fine now, kind of slow but it is running!

Thank you so much!!

Roger Zoh, Ph.D.

Department of Statistics

and Nutrition

Blocker 525B and Kleberg 318

Texas A&M University

[Rob Goedman]

Hi Maxwell,

Thanks for the update. Can you send me the updated dat.csv? Would like to see where the slowness occurs. It should run the 4 chains in parallel. 

Regards,

Rob

Sent from Rob Goedman's iPad Mini

[maxwell]

Hi Rob,

Sure please find attached the updated version of the data set dat. Also, I made a small change in the model(see the red part):

// likelihood

  for(ii in 1:Ns){

    i <- dat[ii,1]; // animal id

    d <- dat[ii,2]; // find the group 

    j <- dat[ii,3]; // finds what attribute (I changed this)

    W[i] ~ normal(X[i],sigu);

    X[i] ~ normal(al0d[d]+Zid[i]*al3d[d], signu[d]); // measurement error part

    eta[i] ~ normal(0,1);

    bijd[i,j] ~ multi_normal(mu0,lamb[j]);

    Y[ii] ~ multi_normal(rep_vector(X[i] * kap1jd[j,d] + eta[i] * bet2jd[j,d] + kap2jd[j,d] * Zid[i],tm) +  Qj * Cjd[j,d] + Qj*bijd[i,j], sigj[j]*eyetm);

  }

Thanks.

[Rob J. Goedman]

Hi Maxwell,

Thanks for the updated dat.csv. It sure is not fast, roughly 100 mins warm-up + 100 mins sampling for each chain. Unfortunately that’s not something I can help with.

Plotting a few of the traces (see below pdf file) do seem to indicate some other problems with the model (I think). Also, I wonder if 1000 iterations (vs. say 200?) add a lot in this case?

I asked for the updated dat.csv file to see how long describe() took, but that is negligible in this case.

[maxwell]

Hi Rob,

I agree!  I need to constrain some parameters. But I am not able to run "describe" on the MCMC output. I get the following error

julia> typeof(sim1)

Chains (constructor with 6 methods)

julia> describe(sim1)

ERROR: `describe` has no method matching describe(::Chains)

Do you have any idea what's wrong?

Thank you

[Rob J. Goedman]

Hi Maxwell,

Can you maybe check which package version of Stan you are running. I made a very late change, for Windows, to handle CRLF issues. The master branch will have that update included. You can install the ‘master’ version by running Pkg.checkout(“Stan”) in Julia’s REPL.

I should have realesed an updated version but wanted to wait on some updates to Stan itself and the Docile documentation on Julia 0.3.x before tagging Stan.jl master as Stan-0.2.1. I’ll try to complete the online docs early next week. In the meantime I believe Stan v2.6.3 might also have been released. 

With Stan-0.2.0+ (the master version), describe() and plotting should work.

Any of the examples, like test_bernoulli.jl in the Test subdirectory, should confirm that.

Regards,

Rob J. Goedman
goe...@mac.com

[Bob Carpenter]

Maxwell --- if you want to send the model along to the list,
there are many of us who can offer help on the statistics side.
The problem you're having is that each chain appears to be stuck
near a different posterior mode.

Rob --- those graphics look great --- are they built
into Stan.jl?

- Bob

[Rob J. Goedman]

Hi Bob,

All kudo’s for the graphics should go to Brian Smith and the Mamba package! Stan.jl ‘borrows’ posterior tools and graphics from Mamba.jl.

[maxwell]

Hi Rob and Bob,

Rob, I have been able to run the test_bernouilli.jl example, the code runs perfectly except the functions describe and autocor. I even run an example provided in the mamba tutorial and I still get the same error:

julia> describe(sim)

ERROR: `describe` has no method matching describe(::Chains)

So,I think my issues seems way beyond Stan.jl. Rob, do you have any idea what maybe the culprit?  I am showing my Julia version information below

julia> versioninfo()

Julia Version 0.3.6

Commit 0c24dca (2015-02-17 22:12 UTC)

Platform Info:

  System: Windows (x86_64-w64-mingw32)

  CPU: Intel(R) Core(TM) i7-2600 CPU @ 3.40GHz

  WORD_SIZE: 64

  BLAS: libopenblas (USE64BITINT DYNAMIC_ARCH NO_AFFINITY Sandybridge)

  LAPACK: libopenblas

  LIBM: libopenlibm

  LLVM: libLLVM-3.3

Bob, I am attaching the Stan file, the data and the Rcode. Let me know I am missing anything.

Thanks,

[Rob J. Goedman]

Hi Maxwell,

Hmm, that’s weird, did you checkout the latest, master version of Stan.jl, i.e. Pkg.checkout(“Stan”)?

I just ran Pkg.test(“Mamba”) on Parallel Desktop (under Windows 10 Tech Preview) without any problems. I’ll try to install Stan itself. I did this a few months ago for testing on Windows 8.1 and had no problems (except finding time, aren’t we all?).

As a short term alternative, have you ever used JuliaBox? If you are in a hurry to get your model running that might be a quick way to get going while I try to sort out the describe()-on-Windows problem. For JuliaBox you just need a Google account.

<R_code_fmi.R><fmi.rstan><dat.csv><Zid.csv><W.csv><Y.csv><Basis.csv>

[maxwell]

Hi Rob,

YES! it worked! I did Pkg.test("Mamba") and I restarted Julia and after that everything was just ok! 

Many thanks Rob!

[Robert J Goedman]

That's great! I am still embarrased I didn't tag the current master and will take care of that next week.

Now the real work, to get model convergence, starts!

Regards,

Rob

Sent from my iPhone