Reverse Random-Walk Election Forecast: Sampling Error: Rejecting initial value: normal_lpdf: Random variable is nan, but must not be nan!


#1

Dear all,

this is my first post here, so if I have a too complex or too little explained minimal example model or made any other mistake, please tell me and I will retry.

Essentially, what I am trying to do is apply Linzer2013s reverse-random walk election forecasting model (as implemented in STAN by Kremp2016 here(see Andrew Gelman’s explanation of the model here) to the German electoral context (at this stage only one national level, but multinomial (multi-party) dependent variable, with inspiration credits to Stoetzer et al. 2017.

The key set-up of the model is as follows: a reverse random-walk starting with a structural/random prior on election day is estimated, with the log party support vector of today being a function of tomorrows value and some random variance. Polls on a given day are treated as multinomial draws from the latent (unobserved) party support vector at that moment.

Attached, I have provided a full minimal example with data-csv, R and STAN-script. Instead of describing every step at length and making this post very long, I have opted for providing the code along with detailed comments at every step. Please tell me if you had prefered a more detailed description at this place.

After a week of trial and error, the model finally builds, but does fails to compile with the following error message:

SAMPLING FOR MODEL ‘minimal_model’ NOW (CHAIN 2).
Rejecting initial value:
Error evaluating the log probability at the initial value.
Exception: normal_lpdf: Random variable is nan, but must not be nan! (in ‘modele4b066842876_minimal_model’ at line 78)

polls.csv (6.4 KB)

minimal_model.stan (3.2 KB)

minimal_example.R (1.8 KB)

Thank you so much for your help and warm welcome to your community.
Marcel
.

STAN MODEL

data {
  int<lower=0> nParties; // number of parties
  int<lower=0> nPolls; // number of polls
  int<lower=0> nPeriods; // number of 
  vector[nPolls] samplesize; //number of Observations in a given sample
  int date[nPolls]; //days until election
  matrix[nParties, 2] normal_priors; // normal-prior matrix (alternatively, had tried it with beta priors, but did not work either)
  int y[nPolls, nParties]; // outcome matrix: absolute counts of votes per party (column) per poll (row)
}

transformed data{
}

parameters {
  real alphanorm[nParties]; //vector of structural priors
  simplex[nParties] theta[nPolls];// theta to be used to as input for multinomial, obtained with dirichlet of vote shares
  real s[nParties-1]; //evolution variance

}

transformed parameters {
  matrix[nPeriods,nParties] alphastar; //daily party support vector: determined by alphanorm for election day, otherwise value for tomorrow plus random evolution variance
  matrix[nPeriods,nParties] ea; // see below
  matrix[nPeriods,nParties] alpha; // see below
  real t[nParties-1]; //transformed evolution variance
  real alrforcast[nParties]; //final forecast
  real ef[nParties];
  real forcast[nParties];
  vector[nParties] shares; 


  //log transformation of structural priors
	for (j in 1:(nParties)) {
  alphastar[nPeriods,j] = log(alphanorm[j]/alphanorm[nParties]);
	}
	
	//alphastar for [,nParties] (other parties) set to 0
	for (i in 1:(nPeriods-1)) {
			alphastar[i,nParties] = 0; 
	}
	
	//convert log shares from alphastar back to vote shares
	for (i in 1:(nPeriods)) {
		for (j in 1:(nParties)) {
		  ea[i,j] = exp(alphastar[i,j]);
    	alpha[i,j] = ea[i,j]/sum(ea[i,]);
    }
  }
    
  //transposition from row vector to vector
  for (k in 1:nPolls){
    shares = alpha[date[k],]'; 
  }
		
	//transformation of parameter for evolution variance s	
	for (j in 1:(nParties-1)) {
	  t[j] = pow(s[j], -2);
	}
		
	//final forecast
	for (j in 1:(nParties)) {
	  alrforcast[j] = alphastar[nPeriods,j];
	  ef[j] = exp(alrforcast[j]);
    forcast[j] = ef[j]/sum(ef[]);
	}
}

model {
	// Walk starts on election day with forcast based on structural prior
	for (j in 1:(nParties)) {
	  target += normal_lpdf(alphanorm[j] | normal_priors[j,1],normal_priors[j,1]);
	}
		
	// Backward-Random-Walk 
	for (i in 1:(nPeriods-1)) {
		for (j in 1:(nParties-1)) {
		  //alphastar realization of today is result of tomorrow plus random variance t
  		target += normal_lpdf(alphastar[i,j] | alphastar[i+1,j],t[j]); 
  		  //this is where the error seems to appear:
  			//("Rejecting initial value: Error evaluating the log probability at the initial value.Exception: normal_lpdf:   
  			// Random variable is nan, but must not be nan!")
	  }	
	}
		
  // Likelihood Model
	for (k in 1:nPolls) {
    theta[k] ~ dirichlet(shares); //shares is the party support vote share vector for the day in which poll k was conducted
		y[k,1:nParties] ~ multinomial(theta[k]); //y[k,] is the outcome matrix with absolute counts for poll responses
		}
		
  // evolution variance with uniform prior (see transformed parameters: t[j] = pow(s[j], -2);)
	for (j in 1:(nParties-1)) {
	  s[j] ~ uniform(0,2);
	}
}

R-Script

load("stanData.rdata")

library(rstan)
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())

polls <- read.csv("polls.csv")


party_names <- grep("cdu|spd|fdp|lin|afd|gru|other|oth", names(polls),value = T)
Y <- round(as.matrix(polls[,grep("cdu|spd|fdp|lin|afd|gru|other|oth", names(polls))]/100) * polls$sample_size) 
NObs <- apply(Y,1,sum)
nParties <-ncol(Y)
time_window <- c(365)

# Normal Priors (from external structural forecast)
normal_priors <- 
  data.frame(
    prior_mean = c(0.09096958,
                   0.35017978,
                   0.06393053,
                   0.07916824,
                   0.08401857,
                   0.05717943,
                   0.27455386),
    prior_sigma = c(0.02276087,
                    0.04404920,
                    0.02241541,
                    0.02136944,
                    0.02218496,
                    0.02231709,
                    0.02953791)
  )
rownames(normal_priors) <- 
  c("afd","cdu","fdp","gru","lin","oth","spd")


stanData <- list(y = Y,
                 nParties = nParties,
                 nPeriods =  time_window[1]+1,
                 nPolls = nrow(Y),
                 date = polls$t,
                 samplesize = NObs, 
                 normal_priors = normal_priors[party_names,] # Order others last
)


fit1 <- stan(
  file = "minimal_model.stan",  # Stan program
  data = stanData,    # named list of data
  chains = 4,             # number of Markov chains
  warmup = 1000,          # number of warmup iterations per chain
  iter = 2000,            # total number of iterations per chain
  cores = 4,              # number of cores (using 2 just for the vignette)
  refresh = 100         # show progress every 'refresh' iterations
)

summary(fit1)

#2

Sorry this one’s taken a long time to answer. Everyone’s been out over the holidays and working toward StanCon and I’m only now getting caught up.

It actually compiled OK, but failed at runtime during execution.

That’s just saying that on line 78 of your model, you are providing a not-a-number value as a y in normal_lpdf(y | mu, sigma) or y ~ normal(mu, sigma).

I don’t want to count down to line 78 of your program, but it’s usually a matter of using variables before they’re defined. It can also come from numerical errors along the way, like dividing zero by zero.