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 sampling does not work when execute

**jonah**#22

`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.

**ghjk**#23

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.

**ghjk**#24

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`

?

**jonah**#25

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:

**ghjk**#26

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.

**ghjk**#27

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

**bgoodri**#29

That error has nothing to do with Stan, but it looks like you need to install XQuartz from the Apple website.

**ghjk**#30

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?

**bgoodri**#31

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â

**ghjk**#32

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 46`24: #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.

**ghjk**#34

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

**ghjk**#35

@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)

**bgoodri**#36

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.

**ghjk**#37

@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.

**bgoodri**#38

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.

**ghjk**#39

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

**bgoodri**#40

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.