Let’s say now I define a function that needs a vector of positions, and mean and standard deviation of a univariate Gaussian distribution, then returns the probability density of Gaussian at those positions as a vector.

Below is a naive way to implement this function (using for loop)

functions{

vector gaussian_pdf(real mu, real sigma, vector z){

vector[num_elements(z)] z_pdf;

for (i in 1:num_elements(z)){

z_pdf[i] = exp( -(z[i] - mu)^2 / (2 * pow(sigma,2) )) * (1/sqrt(2 * pi() * pow(sigma,2)));

}

return z_pdf;

}

}

I was wondering that if there is a way for me to write a more efficient function to achieve the goal?