copula functions

[DomDF]

When trying to use the copula functions (either by including a path to the functions when defining the Stan model, or by copying the definitions into the functions block of the stan file), I see the following error: "A returning function was expected but an undeclared identifier 'clayton_copula' was supplied."

Can someone please advise on what I am missing? Many thanks!

[DomDF]

I no longer receive this error when the copulas are not defined as _lpdf in the functions block.

Renaming the functions as, for example, real clayton_copula() instead of real clayton_copula_lpdf so that it has the same name as used in the model block, allows me to run the Stan file.

Renaming the function in the model block to clayton_copula_lpdf() does not work because of the arguments it expects from an _lpdf function.

Please let me know if I have misunderstood the original functions - otherwise, I hope this helps!

[spinkney]

That makes sense, I guess in order to use the _lpdf suffix it would need to take in a vector or array as the first argument. Stan does not allow two or more different arguments to be _lpdf.

Does it work if you do

real clayton_copula_lpdf(vector uv, real theta) {
  if (theta <= 0.0) 
    reject("clayton_copula: theta must > 0");
  
  return log1p(theta) - (theta + 1) * sum(log(uv))
         - (1 + 2 * theta) / theta * log(  sum(pow(uv, -theta)) - 1);
}
...

where uv will be your distributions (i.e. vector[2] uv = [log_normal_cdf(), log_normal_cdf()]';).

[DomDF]

It would be great to have an _lpdf defined defined to give us the functionality of a custom distribution.

The pow function requires 2 reals, (and not a vector and a real) as inputs, but I managed to get the below to work.
The input is a row_vector of size 2, and I just used u as the first element, and v as the second:

functions {
    
  real clayton_copula_lpdf(row_vector uv, real theta) {
    if (theta <= 0.0) 
      reject("clayton_copula: theta must > 0");
  
    return log1p(theta) - (theta + 1) * (log(uv[1]) + log(uv[2]))
           - (1 + 2 * theta) / theta * log(pow(uv[1], -theta) + pow(uv[2], -theta) - 1);
  }

}

[spinkney]

There is an elementwise operator defined for powers as .^ so I think it's just

real clayton_copula_lpdf(vector uv, real theta) {
  if (theta <= 0.0) 
    reject("clayton_copula: theta must > 0");
  
  return log1p(theta) - (theta + 1) * sum(log(uv))
         - (1 + 2 * theta) / theta * log(  sum( uv .^ -theta ) - 1);
}

I will update the repo accordingly once confirmed