a STAN function to bivariate frank copula of discrete poissons

[andre.p...@googlemail.com]

Hi,

attached my coding of the bivariate frank copula of discrete poissons.

(I put it under the MIT license, no big thing, everybody is free to 

use it)

Technical background: 

https://en.wikipedia.org/wiki/Copula_(probability_theory)

There are more copula, same easy to use.

Since every copula can be combined with every discrete probability function

and as far as I know STAN doesn't allow functions as argument to functions, for

every combination one needs to write a separate function. 

Second, for high value of y1, y2 the function produces NaN values.

  real fcop(real tau,real mu,real ni) {

    return ( -1 / tau * log(1 + (exp(-tau * mu ) - 1) * (exp(-tau * ni) - 1) / (exp(-tau)-1)));

  }

  

  real poisson_frank_copula_log(int y1, int y2, real tau, real la1, real la2) {

    if(y1==0 && y2==0) 

      return(log(fcop(tau,poisson_cdf(0,la1),poisson_cdf(0,la2))));

    else if(y1>0 && y2==0)  

      return(log(fcop(tau,poisson_cdf(y1,la1),poisson_cdf(0,la2)) - fcop(tau,poisson_cdf(y1-1,la1),poisson_cdf(0,la2))));

    else if(y1==0 && y2>0) 

      return(log(fcop(tau,poisson_cdf(0,la1),poisson_cdf(y2,la2)) - fcop(tau,poisson_cdf(0,la1),poisson_cdf(y2-1,la2))));

      

    return(log(fcop(tau,poisson_cdf(y1,la1) ,poisson_cdf(y2,la2)) 

                 - fcop(tau,poisson_cdf(y1-1,la1),poisson_cdf(y2,la2)) 

                 - fcop(tau,poisson_cdf(y1,la1) ,poisson_cdf(y2-1,la2)) 

                 + fcop(tau,poisson_cdf(y1-1,la1),poisson_cdf(y2-1,la2)))); 

    }

Andre

[Bob Carpenter]

I don't think anyone ever thanked you for making it
available. Ben's interested in building in some copula
functionality.

> On Jun 14, 2016, at 9:33 AM, andre.pfeuffer via Stan users mailing list <stan-...@googlegroups.com> wrote:
>
> Hi,
>
> attached my coding of the bivariate frank copula of discrete poissons.
> (I put it under the MIT license, no big thing, everybody is free to

> use it without.)

>
> Technical background:
> https://en.wikipedia.org/wiki/Copula_(probability_theory)
>
> There are more copula, same easy to use.
>
> Since every copula can be combined with every discrete probability function

> and as far as I know STAN down allow function as argument to functions, for

> every combination one needs to write a separate function.
>
> Second, for high value of y1, y2 the function produces NaN values.

I'm not sure what the problem is, but given that you're
taking the log of a sum of a bunch of evaluations of fcop:

> real fcop(real tau,real mu,real ni) {
> return ( -1 / tau * log(1 + (exp(-tau * mu ) - 1) * (exp(-tau * ni) - 1) / (exp(-tau)-1)));
> }

it might be more stable to work on the log scale and then
combine with log_sum_exp and log_diff_exp.

real log_fcop(...)
return log(tau) + (exp(-tau * mu ) - 1) * (exp(-tau * ni) - 1) / (exp(-tau) - 1)) + 1;

Also, there's a log1p(x) function that's a more stable
version of log(1 + x) if you stick with your version.

>
> real poisson_frank_copula_log(int y1, int y2, real tau, real la1, real la2) {
> if(y1==0 && y2==0)
> return(log(fcop(tau,poisson_cdf(0,la1),poisson_cdf(0,la2))));
> else if(y1>0 && y2==0)
> return(log(fcop(tau,poisson_cdf(y1,la1),poisson_cdf(0,la2)) - fcop(tau,poisson_cdf(y1-1,la1),poisson_cdf(0,la2))));
> else if(y1==0 && y2>0)
> return(log(fcop(tau,poisson_cdf(0,la1),poisson_cdf(y2,la2)) - fcop(tau,poisson_cdf(0,la1),poisson_cdf(y2-1,la2))));
>
> return(log(fcop(tau,poisson_cdf(y1,la1) ,poisson_cdf(y2,la2))
> - fcop(tau,poisson_cdf(y1-1,la1),poisson_cdf(y2,la2))
> - fcop(tau,poisson_cdf(y1,la1) ,poisson_cdf(y2-1,la2))
> + fcop(tau,poisson_cdf(y1-1,la1),poisson_cdf(y2-1,la2))));
> }

It'd also be possible to vectorize this. It's more usual to
code the final case as a final else, or to code the first two
"else if" cases as just "if" given that they return. I just
like to keep a consistent coding style so I have an easier time
reading my own code in the future.

By the way, Stan doesn't need parens around returns.

- Bob

[andre.p...@googlemail.com]

Hi,

the discrete copulas are based on differences in the probability space. 

To speed up the calculation we could save the values of the cdf with 

y - 1 for y. 

And we could use log1p with does better to prevent NaNs as Bob mentioned. 

Thus my code above was "unoptimized", but shown more clearly what is going 

on. 

Unfortunately we have log(1 + exp(...) * exp(...) / exp(...)) and if we take the

log again we have log(log(1 + ... )). Thus this way (to my 2cents) its the end 

of optimization possibilities.

Andre