Point estimates

Modeling
1

Hi,
I am try to do a multivariate Bayesian regression as follows
r = S*w + e
where r is a 9x1 vector, S is a 9x13 matrix, w is a 13x1 vector and e ~ t(v, 0, Sigma), i.e. an multivariate t distribution. I am using Rstan. The data samples for the simulation ranges from 1000-5000. i.e. an expanding window approach where the likelihood will include 1000 to 5000 samples. However, I can only get results for the first 10 runs or so(with samples 1000-1010). Then I cannot get any result back. The console tracing of the simulation shows no error message. From my previous experience, it is a typical sign of running out of memory. The stan code I used is the following :

		data {
			int<lower=1> K;  
			int<lower=1> M;   
			int<lower=1> N;   // # of samples points, usually 5000 
			vector[K] r[N];   
			matrix[K,M] S[N]; 
			int<lower=1> s_i; // usually 0
			int<lower=s_i+1> e_i; // start from 1000 then increase w each run
		}

		parameters {
			vector[M] w;
			real<lower=2, upper=20> nu;
			cholesky_factor_corr[K] L_Omega;
			vector<lower=0, upper=pi()/2>[K] sigma_unif;
		}

		transformed parameters {
			vector<lower=0>[K] sigma;
			for (k in 1:K) sigma[k] = 2.5*tan(sigma_unif[k]);
		}

		model {
			matrix[K, K] Sigma;
			w ~ normal(0, 1);
			L_Omega ~ lkj_corr_cholesky(2);
			Sigma = tcrossprod(diag_pre_multiply(sigma, L_Omega));
			for (i in s_i:e_i) {
				r[i] ~ multi_student_t(nu, S[i]*w, Sigma);
			}
		} 

		generated quantities {
			real total_err;
			total_err = sqrt(quad_form(tcrossprod(L_Omega), sigma)*nu/(nu-2));
		}

in order to get some simple result back I will try optimizing(…) to get an MAP estimation. But I previous experience with optimizing is the result often contain jumps which is not surprising considering it is an estimate of local minimal. Then according to the stan documentation. I can also get other point estimates like posterior mean and median. The document refers to each interface docs for details. However, I can not find anything on rstan document indicating how can I get these posterior mean and median. Can someone shed some light here. Also, about the memory consumption. A problem like I have what size of RAM is needed on a Linux machine?

2

You should not run out of memory at 1010 observations for this model. What is the console output?

3

Hi Ben,
The console does not show any error message, it actually showed the simulation had finished. But when I try to the get the result back via mean(extract(fit,‘w[1]’)$‘w[1]’) for example, I get nothing.

Jicun

4

Try it with just one chain, but it seems odd that it would have enough memory to run and then exhaust its memory trying to return the output.

5

How do you call sampling to get a fit? How long does the sampling take?

6

One thing I need to mention is I am running rstan via Rserve, which might have some extra cost. I will try to set chains =1. Currently I am using optimizing to get the MAP estimate. Do you have any idea how to get posterior mean/median as indicated in stan doc 31.2 and 31.3 ?

Jicun

7

fit<-stan(model_code=stan_code, data=stan_data, iter=3000)
It takes at least 30min.

jicun

8

You get summary statistic table by

print(fit)

Posterior mean / median are just mean / median of the mcmc sample.

9

There is no way to get the posterior mean or median from the posterior mode unless you just assume all three are equal, which is rarely the case in non-trivial models. We need to get the Rserve thing figured out so that you can use a genuine posterior distribution.

10

OK, that is perfectly reasonable, no magic. Thanks for clarification!

jicun

11

Yes samples from the posterior solve everything. I have encountered similar problems before. After I increase the memory for the machine(virtual), the problem is gone. At that time, I had multiple Rserve sessions which correspond to multiple processes in the Linux machine. Now I have cut down to 1 session only, the next step is to cut down chains.

jicun