student_t R function

Added dstudent and rstudent for non-standard student t.

R/distributions.r

@@ -453,3 +453,16 @@ rmvnorm2 <- function( n , Mu=rep(0,length(sigma)) , sigma=rep(1,length(Mu)) , Rh
     }
 }
 
+# student t
+dstudent <- function(x,nu=2,mu=0,sigma=1,log=FALSE) {
+    #y <- gamma((nu+1)/2)/(gamma(nu/2)*sqrt(pi*nu)*sigma) * ( 1 + (1/nu)*((x-mu)/sigma)^2 )^(-(nu+1)/2)
+    y <- lgamma((nu+1)/2) - lgamma(nu/2) - 0.5*log(pi*nu) - log(sigma) + (-(nu+1)/2)*log( 1 + (1/nu)*((x-mu)/sigma)^2 )
+    if ( log==FALSE ) y <- exp(y)
+    return(y)
+}
+
+rstudent <- function(n,nu=2,mu=0,sigma=1) {
+    y <- rt(n,df=nu)
+    y <- y*sigma + mu
+    return(y)
+}

R/ulam-function.R

@@ -922,6 +922,8 @@ ulam <- function( flist , data , pars , pars_omit , start , chains=1 , cores=1 ,
                             }
                         }
                     }#ii
+                    # if still NA at this point, then maybe no prior for symbol
+                    if ( is.na(the_dims) ) stop(concat("Unable to determine type and dimensions for ",left_symbol[1])," - Do all parameters have priors?")
                     if ( the_dims==1 )
                         the_dims <- "real"
                     else

R/ulam_templates.R

@@ -637,7 +637,7 @@ ulam_dists <- list(
         vectorized = TRUE
     ) ,
     student_t = list(
-        R_name = "dt",
+        R_name = "dstudent",
         Stan_name = "student_t",
         Stan_suffix = "lpdf",
         pars = 3, # nu, mu, sigma

R/ulam_tests.R

@@ -12,6 +12,26 @@ if ( FALSE ) {
             stop( concat( "Hashes not identical:" , paste( match.call() , collapse=" , " ) ) )
     }
 
+    # student t test
+
+    library(rethinking)
+    data(Howell1)
+    d <- Howell1
+    d2 <- d[ d$age >= 18 , ]
+
+    z <- ulam(
+      alist(
+        height ~ student_t( nu, mu , sigma ) ,
+        mu <- a + b*weight ,
+        a ~ normal( 178 , 20 ) ,
+        b ~ normal( 0 , 1 ) ,
+        sigma ~ cauchy( 0 , 2),
+        nu ~ gamma(2, 0.1)
+      ) , warmup=1000,iter=2000,chains=2,cores=2,
+    data=d2 )
+
+    sim2<-sim(t2)
+
     # dbeta2 test
 
     y <- rbeta2( 100 , 0.3 , 3 )
@@ -24,6 +44,23 @@ if ( FALSE ) {
             theta ~ exponential(1)
         ), data=list(y=y) , sample=TRUE )
 
+    # test "no prior" error
+    z <- ulam(
+        alist(
+            y ~ dbeta2( p , theta ),
+            logit(p) <- a,
+            #a ~ normal(0,1.5),
+            theta ~ exponential(1)
+        ), data=list(y=y) , sample=TRUE )
+
+    z <- ulam(
+        alist(
+            y ~ dbeta2( p , theta ),
+            #logit(p) <- a,
+            #a ~ normal(0,1.5),
+            theta ~ exponential(1)
+        ), data=list(y=y) , sample=TRUE )
+
     # binom tests
 
     data( UCBadmit )

man/AMTL.Rd

@@ -1,11 +1,13 @@
 \name{AMTL}
 \alias{AMTL}
+\alias{AMTL_short}
 %- Also NEED an '\alias' for EACH other topic documented here.
 \title{Ante-mortem Tooth Loss Data}
 \description{
   Data from four primate genera on tooth loss and its relationship to age and sex. Used for measurement error example in the textbook.
 }
 \usage{
+data(AMTL_short)
 data(AMTL)
 }
 %- maybe also 'usage' for other objects documented here.