- For speed up, vectorization of for loops (I am aware vectorization may provide limited performance boost Vectorising user-defined multivariate lpdf/lpmf functions - #10 by fergusjchadwick )
- To implement vectorization, replace singe real with a pseudotype
- Pseudotype reals seems to be the correct one. Based on the description of Stan Function Reference 10.8 Vectorization: “Argument positions declared as
reals
may be filled with a real, a one-dimensional array, a vector, or a row-vector.”
But it hasn’t be implemented in rstan yet? Got error message (Non-void) type expected in function argument declaration if reals used. So real (square brackets) used instead. - real works fine with an array of real
- real not compatible with a single real
Testing by re-creating a simple normal probability density function and generating pointwise likelihoods.
Step 4 implemented as the code below:
functions{
real myNormal_lpdf(real[] y, real mu, real sigma){
real LL[dims(y)[1]];
for(i in 1:dims(y)[1]){
LL[i] = -log(sigma) - square(mu-y[i]) / (2 * sigma^2);
}
return sum(LL);
}
}
data{
int N;
real y[N];
// hyperparameters
real prior_mu[2];
real prior_sigma[2];
}
parameters{
real mu;
real<lower=0> sigma;
}
model{
// priors
mu ~ normal(prior_mu[1],prior_mu[2]);
sigma ~ cauchy(prior_sigma[1],prior_sigma[2]);
// likelihood
// for(n in 1:N){ // loop
// target += myNormal_lpdf(y[n] | mu,sigma);
// }
target += myNormal_lpdf(y | mu,sigma); // vectorization
}
The sampling statement works just fine because y was specified as an array of reals real which is the same as defined in the function.
However, the input of a single real is not compatible with real. Problem in step 5 by adding a Generated Quantities block:
functions{
real myNormal_lpdf(real[] y, real mu, real sigma){
real LL[dims(y)[1]];
for(i in 1:dims(y)[1]){
LL[i] = -log(sigma) - square(mu-y[i]) / (2 * sigma^2);
}
return sum(LL);
}
}
data{
int N;
real y[N];
// hyperparameters
real prior_mu[2];
real prior_sigma[2];
}
parameters{
real mu;
real<lower=0> sigma;
}
model{
// priors
mu ~ normal(prior_mu[1],prior_mu[2]);
sigma ~ cauchy(prior_sigma[1],prior_sigma[2]);
// likelihood
// for(n in 1:N){ // loop
// target += myNormal_lpdf(y[n] | mu,sigma);
// }
target += myNormal_lpdf(y | mu,sigma); // vectorization
}
generated quantities{
real log_lik[N];
for(n in 1:N){
log_lik[n] = myNormal_lpdf(y[n] | mu,sigma);
}
}
Got an error message “Ill-typed arguments supplied to function ‘myNormal_lpdf’. Available signatures:
(real[], real, real) => real
Instead supplied arguments of incompatible type: real, real, real.”
I suppose real should be able to accept a one-dimensional array of reals with any number of elements including one of course. Why it isn’t working?
Am I using the pseudotype real correctly?
What’s the best way to use a pseudotype?
Thank you very much!