test.fn <- function(nu){ sqrt(nu/2)*exp(lgamma((nu-1)/2)-lgamma(nu/2)) } ## Murnane-Willett formula mean.mw <- function(delta.mu,n1,n2){ nu <- n1+n2-2 return(delta.mu*test.fn(nu)) } mean.mw(4.899,500,500)