It finally works! The resolution is extremely bad (numbers are blurry), but I was able to zoom in to see it clearly now. Thanks a ton for your help. Your help did help me learn a ton about Stan and MCMC in R, so I really appreciate it. And I did not rename any of my parameters.
mcmc_trace uses ggplot under the hood so the resolution shouldnât be any better or worse that what you would get by making the plot yourself with ggplot. The ggsave function will also let you change the resolution and other stuff when saving a plot.
If âyTâ is a scalar parameter setting pars="yT" should work fine with mcmc_trace. If the parameters are really yT[1], yT[2], etc. then you can use those names or use regex_pars="yT" to include all parameters with âyTâ in their names.
1 Like
Thank you for your great suggestion, Jonah! I tried regex_pars="yT", but got this error:
Error: No matches for 'regex_pars'
In addition: Warning message: In is.na(<S4 object of class "stanfit">) : is.na() applied to non-(list or vector) of type 'S4'
You are right, my yT is a vector, not scalar.
Btw, I have another question on this topic: when I extract the simulated output yT out of the stanfit object by using the command head(list_of_draws$yT), I donât understand why the dimension of yT is 6x24 rather than 1x24. I meant, where does the 6 come from in this case? It is not the number of iterations, is it??
Based on what it said here about the length of parameters mu and tau, it is equal to the number of post-warmup iterations times the number of chains. But my iterations = 1500, warmup = 1000, and chains = 4, so shouldnât I get 500*4 = 2000 as the length of yT?
I suggest trying something like this (replacing âthetaâ with the appropriate parameter name from your Stan program:
# assuming you have a stanfit object "stanfit" corresponding to a
# Stan program with a vector/array parameter "theta"
posterior_theta <- as.array(stanfit, pars = "theta")
# see what the parameter names are
print(dimnames(posterior_theta)[[3]])
# suppose the names are theta[1], theta[2], ..., theta[K]
library("bayesplot") # also make sure bayesplot is up-to-date (latest release is v1.4.0)
color_scheme_set("mix-brightblue-darkgray")
# trace plot for theta[1]
mcmc_trace(posterior_theta, pars = "theta[1]")
# trace plots for theta[1] and theta[2]
mcmc_trace(posterior_theta, pars = c("theta[1]", "theta[2]"))
# trace plots for theta[1] and theta[2] (alternate method)
mcmc_trace(posterior_theta[, , 1:2])
The are also many other plots in the bayesplot package. The vignettes have some good examples:
1 Like
Thank you a ton for your help, Jonah! It worked out well for me now, so helps settle what I want above:)
Unfortunately, I need to change the modeling part a little bit due to the nature of a project I am working on. When I did make a change on changing sigma_x and sigma_y from a vector to a matrix, and re-ran my program, it gave me the error message:
normal(A*yT+epsilon, sigma_x) cannot take on the form [vector, matrix].
The reason this problem arises is because each of our parameters yT and x are **vectors **, and so the variance of yT, which is sigma_y, theoretically must be a co-variance matrix with size [nrow, nrow]. Similar, sigma_x should also be a co-variance matrix with size [ncol, ncol].
Question: I wonder if you (or anyone in this forum) have encountered a similar model, where we are basically given the normal distributions of vectors, rather than random variables, and try to make Bayesian inference for the posterior distribution in the same way as my model in the original post? If anyone did, could you show me how you overcame the problem of having two vectors following normal distributions, which have different dimensions for their corresponding means and variance.
Btw, I tried installing the matlib library into R (version 3.3.3), but then I got this error message:
Error : .onLoad failed in loadNamespace() for 'rgl', details:
call: dyn.load(file, DLLpath = DLLpath, ...)
error: unable to load shared object '/Library/Frameworks/R.framework/Versions/3.3/Resources/library/rgl/libs/rgl.so':
dlopen(/Library/Frameworks/R.framework/Versions/3.3/Resources/library/rgl/libs/rgl.so, 6): Library not loaded: /opt/X11/lib/libGLU.1.dylib
Referenced from: /Library/Frameworks/R.framework/Versions/3.3/Resources/library/rgl/libs/rgl.so
Reason: image not found
Error: package or namespace load failed for âmatlibâ
Is it due to the file directory path? I am not sure how to fix it;p
Nobody wants to help me resolve this new issue?
That error has nothing to do with Stan, but it looks like you need to install XQuartz from the Apple website.
1 Like
Thank you so much for helping me once again, Sir!! Does it mean matlib library does not work in the OS environment? From what I just read here (https://www.xquartz.org/index.html), XQuart creates a virtual Window system being run on OS system. Is there a similar package to matlib specifically designed for R?
I donât even know what matlib is, but I do know it is not Stan and that this is not a general R help website. The XQuartz thing is just the first hit when I googled for
rgl âLibrary not loaded: /opt/X11/lib/libGLU.1.dylibâ
Thanks a lot for your help. I figured out that I do not actually need matlib. But I have another question:
I want to change the likelihood function in Stan above from Normal to Poisson, but when I wrote x ~ poisson(to_vector(A*yT + epsilon)) or using the for loop to simulate each entry of the column vector x, I got this error message:
Warning (non-fatal): Comments beginning with # are deprecated. Please use // in place of # for line comments.
No matches for:
real ~ poisson(real)
Available argument signatures for poisson:
int ~ poisson(real)
int ~ poisson(real)
int ~ poisson(vector)
int ~ poisson(row vector)
int ~ poisson(real)
int ~ poisson(real)
int ~ poisson(vector)
int ~ poisson(row vector)require real scalar return type for probability function.
error in âmodel25d7d6722c9_stanmodel1â at line 26, column 4624: #x ~ normal(A*yT + epsilon, sigma_x); 25: for (i in 1:nrow) 26: x[i] ~ poisson((A*yT)[i] + epsilon[i]); ^
Error in stanc(file = file, model_code = model_code, model_name = model_name, :
failed to parse Stan model âstanmodel1â due to the above error.
Could you please help me with this issue? I donât see why I got this error with the for loop.
Because you declared x as a real array or vector instead of an integer array.
1 Like
Many thanks for your help, Sir! I fixed that to become int x[nrow], and it worked out nicely. Did not know we could declare an array in this way;p
@bgoodri: I have another question, could you please help? Due to modeling purpose, I would need to transform the parameter yT in the Stan model from vector to an integer array, because Stan only accept things like yT ~ poisson(yh) where yh is a vector and yT is integer or int[]. I tried specifying yT as integer array, but when I fit the model, it said the following errors:
parameters or transformed parameters cannot be integer or integer array; found declared type int, parameter name=yT
Problem with declaration.
error in âmodel2245b774e56_stanmodel1â at line 14, column 16
12: parameters {
13: #vector[ncol] yT;
14: int yT[ncol];
^
15: }
Itâs annoying that Stan does not allow transforming from one datatype to another datatype (i.e, from vector to integer array, or from column matrix to column vector)
You canât have an int or int array in the parameters block (only data, transformed data, generated quantities, and local blocks). You cannot do HMC with discrete unknowns.
@bgoodri: Thank you for your help. But I donât understand: you mean HMC does not work with prior ~ poisson?? The current research paper that I am reading claimed to do just that, so I am not sure why we canât work with priors ~ poisson? Actually, I only need to make the part about prior ~ poisson work, not that the parameters have to be integer.
In order for yT ~ poisson to make sense as a prior for yT, yT would have to be both a parameter and an integer array and that combination is impossible for HMC.
1 Like
Thank you so much for your explanation. So do we have any way to model this situation? I read about taking poisson_log or poisson_lpdf, but even that does not work with integer parameters, as required by the default input of Stan;p
You have to marginalize discrete unknowns out of the posterior distribution (and possibly bring them back in generated quantities). With unknown counts, that is difficult because there are an infinite number of things a count could be, but truncating the sum at some huge number is probably okay.