REPORT(); not working for a specific variabe

[AndreasKarp]

Hi,

I am not sure if this is the right place to ask this question, but I am unable to find a solution myself and I don't know where else I can ask the question. I am trying to implement the following hierarhical ordered probit model with the brilliant TMB package. I assume that model runs fine, because it converges. I however have a problem with the REPORT() macro. lt, ut, and si are exported but mu prints the starting values of blat instead of the product of the matrix latent and the vector blat (in R it would be latent %*% blat). That is strange because because the estimated parameter values for blat make sense and also the results I get when I do this matrix multiplication outside C++ in R.
Does anybody have an idea what the reason can be?

// Hopit Model 1.0
#include </home/ak/R/x86_64-unknown-linux-gnu-library/3.1/TMB/include/TMB.hpp>

template<class Type>
Type objective_function<Type>::operator() ()
{
  //data: 
  DATA_VECTOR(y);
  DATA_MATRIX(threshold);
  DATA_MATRIX(variance);
  DATA_MATRIX(latent);

  //parameters:
  PARAMETER(y0);
  PARAMETER_VECTOR(blow);
  PARAMETER_VECTOR(bupp);
  PARAMETER_VECTOR(bvar);
  PARAMETER_VECTOR(blat);

  //procedures: (transformed parameters)
  using namespace density;
  int i;
  Type nll = 0.0; //initialize negative log likelihood

  int n = y.size();

  vector<Type> lt(n);
  vector<Type> ut(n);
  vector<Type> si(n);
  vector<Type> mu(n);

  lt  = y0 + exp(threshold*blow);
  ut  = lt + exp(threshold*bupp);  
  si  = exp(variance*bvar);
  mu  = latent *blat;

  //reports on transformed parameters:
  REPORT(mu);
  REPORT(lt);
  REPORT(ut);
  REPORT(si);

  //get negative log likelihood 
  for(i=0;i<n;i++) {
  if(y(i)==1) nll -= log(pnorm(lt(i),mu(i),si(i))); 
  if(y(i)==2) nll -= log(pnorm(ut(i),mu(i),si(i)) - pnorm(lt(i),mu(i),si(i)));
  if(y(i)==3) nll -= log(1 - pnorm(ut(i),mu(i),si(i)));
  }

  return nll;
}

[AndreasKarp]

It was quite late already yesterday when I wrote the question, I now played around with it a bit again and it is working, although I don't understand completely the reason: Simply moving the REPORT section between the loop and the return did the trick.

I however would also have another question for you: I would like to integrate some kind of penalty when the parameters in blow and bupp (the parameters for the lower and upper thresholds) are different. In R I would write something like if(any(sign(blow)!=sign(bupp)) nll = nll + Inf. How could I achieve that in TMB?

[skaug]

Your penealty is non-differentiable, which can lead to problem for the function optimizer. Try a quadratic penalty instead, available easily via dnorm(blow-lupp,0,1,true).sum();

[AndreasKarp]

Hey, thank you very much, that is a brilliant idea!