an example of procedural stan code (no practical meaning).

procedural_stancode.R

library(rstan)

model <- "                                                                                                                                                                             
data {                                                                                                                                                                                 
       int<lower=0> N;  // number of data points                                                                                                                                       
       int<lower=1> D;  // number of dimensions                                                                                                                                        
       real<lower=0,upper=1> param[D];                                                                                                                                                 
}                                                                                                                                                                                      
parameters {                                                                                                                                                                           
       vector[D] phi;                                                                                                                                                                  
}                                                                                                                                                                                      
                                                                                                                                                                                       
transformed parameters {                                                                                                                                                               
  matrix[N,D] mm;                                                                                                                                                                      
  for (k in 1:N){                                                                                                                                                                      
   for (i in 1:D){                                                                                                                                                                     
    for (j in 1:i)                                                                                                                                                                     
      mm[k,j]<-1./i;                                                                                                                                                                   
    for (j in (i+1):N)                                                                                                                                                                 
      mm[k,j]<-0;                                                                                                                                                                      
    }                                                                                                                                                                                  
  }                                                                                                                                                                                    
}                                                                                                                                                                                      
model {                                                                                                                                                                                
  real q;                                                                                                                                                                              
  for (i in 1:D){                                                                                                                                                                      
    q<-0;                                                                                                                                                                              
    for (j in 1:i){                                                                                                                                                                    
    q<-q+1;                                                                                                                                                                            
    }                                                                                                                                                                                  
    phi[i] ~ normal(0,q);                                                                                                                                                              
  }                                                                                                                                                                                    
}                                                                                                                                                                                      
//generated quantities{}                                                                                                                                                               
"
N <- 3
D <- 10
param <- rep(1./D,D)
fit <- stan(model_code = model,
            pars=c("phi"),
            data = list(N=N,
              D=D,
              param=param),
           chains = 1, iter = 20)

#png("trace.png")
#traceplot(fit)
#dev.off()