Question about printing from the model block

Modeling
1

Hi all,

I’m using CmdStanR. I have a question about using the print() function from the model block to check if my transformed parameters block is doing what I want it to do. The program manages to compile and I do get the desired output from the print which is the number of participants in the control arm, but I’m getting multiple prints. How do I get it to just print one time?

R code

rm(list=ls())
set.seed(070125)
ntotal <- 1e1
kactive <- 5
theta_con <- 0.3
theta_trt <- 0.5

v.outcome <- c(rbinom(ntotal, 1, theta_con), 
               rbinom(ntotal, 1, theta_trt),
               rbinom(ntotal, 1, theta_trt),
               rbinom(ntotal, 1, theta_trt),
               rbinom(ntotal, 1, theta_trt),
               rbinom(ntotal, 1, theta_trt))

v.allocations <- c(rep(0,ntotal))
for (k in 1:kactive) {
  
  v.allocations <- c(v.allocations, rep(k, ntotal))
  
}

v.allocation.indices <- unique(v.allocations)

model.stan <- cmdstanr::cmdstan_model(file.path(getwd(), "Program.stan"))

list.data <- list(
  N = length(v.outcome),
  K = kactive,
  K_indices = v.allocation.indices,
  outcome_i = v.outcome,
  K_arm_i = v.allocations
)

model.fit <- model.stan$sample(
  data = list.data,
  chains = 1,
  iter_warmup = 1,
  iter_sampling = 0, 
  parallel_chains = 1
)

R output

Running MCMC with 1 chain...

Chain 1 Num control = 10 
Chain 1 Num control = 10 
Chain 1 WARNING: No variance estimation is 
Chain 1          performed for num_warmup < 20 
Chain 1 Num control = 10 
Chain 1 Num control = 10 
Chain 1 Num control = 10 
Chain 1 Num control = 10 
Chain 1 Iteration: 1 / 1 [100%]  (Warmup) 
Chain 1 finished in 0.0 seconds.

Stan code

data {
  
  // Total sample size
  int<lower=0> N;
  
  // Number of experimental arms active at any given time
  int<lower=0> K;
  
  // Vector of the arm indices. 0 is for control, and anything greater than 0 is a treatment arm index
  vector<lower=0>[K+1] K_indices;

  // Outcome array
  array[N] int<lower=0, upper=1> outcome_i;
  
  // Vector of arm allocations
  vector<lower=0>[N] K_arm_i;
  
}

transformed data {

  // Create indicator variables for each observation
  // vector of 1s and 0s to indicate if an individual is a control in arm i
  // Control arm indicator
  vector<lower=0>[N] control_indicator;
  for (n in 1:N) {
      
      control_indicator[n] = (K_arm_i[n]==0);
  }
  // count the number of control participants
  real<lower=0> N_0 = sum(control_indicator);

  
}

parameters {
  
  real<lower=0, upper=1> theta_con;
  vector<lower=0, upper=1>[K] theta_trt;
}

model {
  
  print("Num control = ", N_0);
}
1 Like
2

This is because in each iteration, the Hamiltonian Monte Carlo algorithm in Stan performs multiple evaluations of the target density (usually called steps), see MCMC Sampling in the Stan manual for more details.

1 Like
3

Thank you for this. I’m guessing it wouldn’t be possible to get it to print once?

4

Managed to find a workaround to get a single print using the stringr package and sink function.

sink(file.path(getwd(), "Program_sink.txt"))
model.fit <- model.stan$sample(
  data = list.data,
  chains = 1,
  iter_warmup = 1,
  iter_sampling = 0, 
  parallel_chains = 1
)
sink()

output.readlines <- readr::read_lines(file.path(getwd(), "Program_sink.txt"))
output.readlines <- output.readlines[max(str_which(output.readlines, "Num control")):length(output.readlines)]
output.readlines
5

In this case you could put the print statement itself in transformed data, and then it would only run once

2 Likes