t distribution function (stats::pt)

[Jeremy D]

Hi,

Can anyone suggest a way to implement the t-distribution function (i.e., stats::pt) in a TMB model?

This function is needed for the skew t-distribution (sn::dst), which is currently not available in TMB.

Thank you!

[Hans Skaug]

I guess that the reason pt() is not already in TMB is that its numerical implementation is a bit tricky.

However, if seems that there are other ways of getting dst:

https://en.wikipedia.org/wiki/Skewed_generalized_t_distribution#Skewed_t_distribution

Hans

[Ben Bolker]

Slightly naive Q: since (AFAIK) the R functions for pt() and dt()
are available in TMB and since the derivative of pt() is dt() (by
definition), it would seem not too hard (for some definition of "hard")
to define a new atomic function following

http://kaskr.github.io/adcomp/_book/AtomicFunctions.html

The code below *almost* compiles, but gets an error

error: expected constructor, destructor, or type conversion before ‘(’ token
34 | VECTORIZE_1t(pt)

There's probably some tiny syntax error (see
<https://stackoverflow.com/questions/8958044/expected-constructor-destructor-or-type-conversion-before-token>:
"If the compiler can't match the definition's signature to the
declaration's signature it thinks the definition is not a constructor
and then doesn't know how to parse the code and displays this error")

Maybe someone else can see what's wrong, or can tell me why this
wouldn't work.

====
#include <TMB.hpp>

extern "C" {
double Rf_pt(double q);
double Rf_dt(double x);
}

TMB_ATOMIC_VECTOR_FUNCTION(
// ATOMIC_NAME
pt
,
// OUTPUT_DIM
1,
// ATOMIC_DOUBLE
ty[0] = Rf_pt(tx[0]); // Call the 'double' version
,
// ATOMIC_REVERSE
Type x = ty[0]; // Function value from forward pass
Type dp = Rf_dt(x); // Derivative (do we need a cast here??)
px[0] = dp * py[0]; // Reverse mode chain rule
)

template<class Type>
Type pt(Type x){
CppAD::vector<Type> tx(1);
tx[0] = x;
return pt(tx)[0];
}

// Vectorized version
VECTORIZE_1t(pt)

> --
> To post to this group, send email to us...@tmb-project.org. Before
> posting, please check the wiki and issuetracker at
> https://github.com/kaskr/adcomp/ <https://github.com/kaskr/adcomp/>.
> Please try to create a simple repeatable example to go with your
> question (e.g issues 154, 134, 51). Use the issuetracker to report bugs.
> ---
> You received this message because you are subscribed to the Google
> Groups "TMB Users" group.
> To unsubscribe from this group and stop receiving emails from it, send
> an email to tmb-users+...@googlegroups.com
> <mailto:tmb-users+...@googlegroups.com>.
> To view this discussion on the web visit
> https://groups.google.com/d/msgid/tmb-users/a5679c00-80d9-4195-904b-0373adfb9401n%40googlegroups.com
> <https://groups.google.com/d/msgid/tmb-users/a5679c00-80d9-4195-904b-0373adfb9401n%40googlegroups.com?utm_medium=email&utm_source=footer>.

--
Dr. Benjamin Bolker
Professor, Mathematics & Statistics and Biology, McMaster University
Director, School of Computational Science and Engineering
(Acting) Graduate chair, Mathematics & Statistics

[Jeremy D]

Thank you both for your replies.

I'm looking into the reformulations of skewed-t that Hank's suggested, but I'm still hoping we could get Ben's code to work so we could use pt, and also understand more generally how the syntax error he mentioned can be fixed.

[Hans Skaug]

Ben is right about using atomic functions, but the problem is that

we only have a formula for the derivative wrt "x", and not wrt to "df".

This is what prevents a straight forward implementation of pt() in TMB (I think):

Below I have written a code where "df" is hardcoded. Then there is no

need for the derivative wrt df, but of course "df" cannot be estimated as a parameter. 

 I do not know if this is  acceptable for Jeremy's application.

(I also skip the vectorization stuff, so the code only works for scalar "x".)  

My code seems to work, but I do not claim to have a full understanding of this, 

so my code should be verified numerically before being used.

Hans

------------------------------------------

// pt() with hardcoded df
#include <TMB.hpp>

extern "C" {
   double Rf_pt(double q,double df,int,int);

}

TMB_ATOMIC_VECTOR_FUNCTION(
     // ATOMIC_NAME
     pt
     ,
     // OUTPUT_DIM
     1,
     // ATOMIC_DOUBLE

     double df=2.0;  // First location of hardcoding
     ty[0] = Rf_pt(tx[0],df,1,0); // Call the 'double' version
     ,
     // ATOMIC_REVERSE
     Type df = 2.0;  // Second location of hardcoding
     Type x  = tx[0];      
     px[0] = dt(x,df,false) * py[0];   // Reverse mode chain rule

                           )
template<class Type>
Type pt(Type x){
   CppAD::vector<Type> tx(1);
   tx[0] = x;
   return pt(tx)[0];
}

template<class Type>
Type objective_function<Type>::operator() ()
{
  DATA_VECTOR(z);    
  PARAMETER_VECTOR(x);                             // Underlying latent random
  return pt(x(0));
}

[Kasper Kristensen]

I agree with the analysis of Hans.

An alternative is to look at the R source https://github.com/wch/r-source/blob/trunk/src/nmath/pt.c#L24 which gives some hints on how to express pt in terms of pbeta (which is available in TMB).

[Kasper Kristensen]

To follow up on my own suggestion here's an R function that should be straight forward to translate to TMB (using CondExpGe for ifelse). You may want to verify against stats::pt for accuracy:

pt <- function(x, df) {

    nx <- 1 + (x/df)*x
    val <- pbeta (1. / nx, df / 2., 0.5)
    val <- val / 2
    ifelse(x>0, 1-val, val)

[Jeremy D]

I can confirm that both Hans' and Kasper's functions seem to work well and, as far as I can tell, give identical result to stats::pt. 

Performance-wise, Kasper's version runs much slower in my tests, but if I fix df (to 3, as well as in Hans' pt version), then running time is similar.

Here's Kasper's version that I used:

template<class Type>
Type pt_3(Type x) {
    Type nx = Type(1) + (x / Type(3)) * x;
    Type val = pbeta(Type(1) / nx, Type(3) / Type(2), Type(0.5)) / Type(2);
    return CppAD::CondExpGt(x, Type(0), Type(1) - val, val);
}

Thanks again for all of your helpful replies!