Skip to content

Commit

Permalink
version 1.0
Browse files Browse the repository at this point in the history
  • Loading branch information
Guy Nason authored and gaborcsardi committed Jan 7, 2005
0 parents commit 093aa6f
Show file tree
Hide file tree
Showing 20 changed files with 654 additions and 0 deletions.
12 changes: 12 additions & 0 deletions DESCRIPTION
@@ -0,0 +1,12 @@
Package: NORMT3
Title: Evaluates complex erf, erfc and density of sum of Gaussian and Student's t
Version: 1.0
Date: 2005-01-07
Author: Guy Nason <g.p.nason@bristol.ac.uk>
Maintainer: Guy Nason <g.p.nason@bristol.ac.uk>
Description: Evaluates the probability density function of the sum of the
Gaussian and Student's t density on 3 degrees of freedom.
Evaluates the p.d.f. of the sphered Student's t density function.
Also evaluates the erf and erfc functions on complex-valued arguments.
License: GPL version 2 or newer
Packaged: Mon Jan 10 09:49:09 2005; magpn
14 changes: 14 additions & 0 deletions R/dnormt3.R
@@ -0,0 +1,14 @@
"dnormt3" <-
function(x, mean=0, sd=1){

if (sd==0) {
f <- rep(0, length(x))
f[x==0] <- Inf
}
else {
f <- sqrt(2)*normt3ip(mu=(x-mean)*sqrt(2)/sd, sigma=1)/sd
sv <- abs(x - mean)*sqrt(2)/sd >= 37
f[sv] <- 0
}
f
}
4 changes: 4 additions & 0 deletions R/dst.R
@@ -0,0 +1,4 @@
"dst" <-
function(x,nu=3){
gamma( (nu+1)/2)/(gamma(nu/2)*sqrt(pi)*sqrt(nu-2)*( 1+ (x^2)/(nu-2))^((nu+1)/2))
}
2 changes: 2 additions & 0 deletions R/erf.R
@@ -0,0 +1,2 @@
"erf" <-
function(z) 1 - erfc(z)
17 changes: 17 additions & 0 deletions R/erfc.R
@@ -0,0 +1,17 @@
"erfc" <-
function(z){

ans <- .C("IPerfcvec",
x=as.double(Re(z)),
y=as.double(Im(z)),
ansx=as.double(Re(z)),
ansy=as.double(Im(z)),
n = as.integer(length(z)),
error = as.integer(0),
PACKAGE="NORMT3")
if (ans$error != 0)
stop(paste("Error code from TOMS 680 was ", ans$error))

return(complex(real=ans$ansx, im=ans$ansy))

}
5 changes: 5 additions & 0 deletions R/firstlib.R
@@ -0,0 +1,5 @@
.First.lib <- function(lib,pkg)
{
library.dynam("NORMT3",pkg,lib)
cat("NORMT3 loaded\n")
}
8 changes: 8 additions & 0 deletions R/ic1.R
@@ -0,0 +1,8 @@
"ic1" <-
function (p,d)
{
cc <- complex(re=0, im=p/2)

(sqrt(pi)/2)*exp(-(p^2)/4)*(1-Re(erf(cc) - erf(cc-d)))

}
8 changes: 8 additions & 0 deletions R/is1.R
@@ -0,0 +1,8 @@
"is1" <-
function (p,d)
{

cc <- complex(re=0, im=p/2)

(sqrt(pi)/2)*exp(-(p^2)/4)*Im(erf(cc-d))
}
17 changes: 17 additions & 0 deletions R/normt3ip.R
@@ -0,0 +1,17 @@
"normt3ip" <-
function (mu, sigma)
{
a <- 1/(sigma^2)
b <- sqrt(2)/sigma
d <- a/b
p <- mu*b

T1 <- ((1-a)*cos(mu*a) + p*b*sin(mu*a)/2)* ic1(p,d)

T2 <- ((1-a)*sin(mu*a) - p*b*cos(mu*a)/2)*is1(p,d)

T3 <- b*exp( - d^2)/2

total <- b*exp(a/2)*(T1+T2+T3)/pi
total
}
10 changes: 10 additions & 0 deletions README
@@ -0,0 +1,10 @@
1. Put any C/C++/Fortran code in 'src'
2. If you have compiled code, add a .First.lib() function in 'R'
to load the shared library
3. Edit the help file skeletons in 'man'
4. Run R CMD build to create the index files
5. Run R CMD check to check the package
6. Run R CMD build to make the package file


Read "Writing R Extensions" for more information.
39 changes: 39 additions & 0 deletions man/dnormt3.Rd
@@ -0,0 +1,39 @@
\name{dnormt3}
\alias{dnormt3}
\title{Density function of sum of Gaussian and Student's t on 3 df}
\description{
Computes the probability density function of the sum of the
Gaussian distribution and the Student's t distribution on
3 degrees of freedom.
}
\usage{
dnormt3(x, mean = 0, sd = 1)
}
\arguments{
\item{x}{ Where to evaluate the density function}
\item{mean}{ The mean of the Gaussian}
\item{sd}{ The standard deviation of the Gaussian}
}
\details{
Computes the probability density function of the sum of the
Gaussian distribution and the Student's t distribution on
3 degrees of freedom.
}
\value{The appropriate pdf value.
}
\references{ Nason, G.P. (2005) On the sum of the Gaussian and Student's t
random variables. \emph{Technical Report}, 05:01,
Statistics Group, Department of Mathematics, University of Bristol.}
\author{ Guy Nason }

\seealso{ \code{\link{normt3ip}},\code{\link{dst}},\code{\link{dnorm}}}

\examples{
dnormt3(0)
#
# Should be 0.4501582 = sqrt(2)/pi
#
x <- seq(from=-5, to=5, length=100)
plot(x, dnormt3(x), type="l") # Density plot
}
\keyword{ distribution }% at least one, from doc/KEYWORDS
34 changes: 34 additions & 0 deletions man/dst.Rd
@@ -0,0 +1,34 @@
\name{dst}
\alias{dst}
\title{ Density function of sphered Student's t distribution}
\description{
Evaluates probability density function of sphered Student's
t distribution on nu degrees of freedom. This is just the
standard Student's t distribution rescaled to have unit variance.
}
\usage{
dst(x, nu = 3)
}
\arguments{
\item{x}{ Where to evaluate the density}
\item{nu}{ The degrees of freedom }
}
\details{Description says it all.
}
\value{
The appropriate probability density.
}
\references{Nason, G.P. (2001) Robust projection indices. \emph{Journal of
the Royal Statistical Society}, Series B, \bold{63}, 551--567. }
\author{ Guy Nason}
\seealso{ \code{\link{dnorm}}, \code{\link{dnormt3}}}
\examples{
dst(0)
#
# Should be 2/pi = 0.6366198
#
x <- seq(from=-5, to=5, length=100)
plot(x, dst(x), type="l") # Produces a density plot
}
\keyword{ distribution }% at least one, from doc/KEYWORDS
45 changes: 45 additions & 0 deletions man/erf.Rd
@@ -0,0 +1,45 @@
\name{erf}
\alias{erf}
\title{Error function}
\description{
Computes the error function of a (possibly) complex
valued argument. This function is 1-erfc(z).
}
\usage{
erf(z)
}
\arguments{
\item{z}{Argument of error function}
}
\details{
Computes the error function of a (possibly) complex
valued argument by computing the complementary error function
and subtracting the answer from 1.
}
\value{
The error function of z
}
\references{
Poppe, G.P.M. and Wijers, C.M.J. (1990) More efficient computation of the
complex error function. \emph{ACM Transactions on Mathematical
Software}, \bold{16}, 38--46.
}
\author{Guy P. Nason, Department of Mathematics, University of Bristol}


\seealso{ \code{\link{erfc}}}
\examples{
erf(0)
#
# Should give 0
#
erf(1)
#
# Should give 0.8427008+0i
#
erf(complex(re=1, im=1))
#
# Should give 1.316151+0.1904535i
#
}
\keyword{math}
49 changes: 49 additions & 0 deletions man/erfc.Rd
@@ -0,0 +1,49 @@
\name{erfc}
\alias{erfc}
\title{Complementary error function}
\description{
Computes the complementary error function of a (possibly) complex
valued argument. This function is
$2/\sqrt{\pi} \int_{z}^{\infty} \exp^{-t^2} dt$.
}
\usage{
erfc(z)
}
\arguments{
\item{z}{Argument of complementary error function}
}
\details{
Computes the complementary error function of a (possibly) complex
valued argument. This function is
$2/\sqrt{\pi} \int_{z}^{\infty} \exp^{-t^2} dt$. This function
actually calls FORTRAN code (algorithm TOMS 680) which computes
the Faddeeva's function and then with a slight modification
obtains the erfc function of a complex-valued argument.
}
\value{
The complementary error function of z
}
\references{
Poppe, G.P.M. and Wijers, C.M.J. (1990) More efficient computation of the
complex error function. \emph{ACM Transactions on Mathematical
Software}, \bold{16}, 38--46.
}
\author{Guy P. Nason, Department of Mathematics, University of Bristol}
\seealso{\code{\link{erf}}}
\examples{
erfc(0)
#
# Should give 1
#
erfc(1)
#
# Should give 0.1572992+0i
#
erfc(complex(re=1, im=1))
#
# Should give -0.3161513-0.1904535i
#
}
\keyword{math}
34 changes: 34 additions & 0 deletions man/ic1.Rd
@@ -0,0 +1,34 @@
\name{ic1}
\alias{ic1}
\title{Compute IC1 formula from Nason (2005) }
\description{
Computes $I_{C1}(p,d) = \int_{d}^{\infty} \cos (px) \exp (-x^2) dx$}
}
\usage{
ic1(p, d)
}
%- maybe also 'usage' for other objects documented here.
\arguments{
\item{p}{ Argument for IC1 function}
\item{d}{ Argument for IC1 function}
}
\details{
An intermediate function for the computation of the inner product
between Gaussian and Student's t distribution.
}
\value{
Value of the function IC1
}
\references{ Nason, G.P. (2005) On the sum of the Gaussian and Student's t
random variables. \emph{Technical Report}, 05:01,
Statistics Group, Department of Mathematics, University of Bristol.}
\author{ Guy Nason }

\seealso{ \code{\link{normt3ip}}}
\examples{
ic1(1,1)
#
# Should give 0.03401986.
#
}
\keyword{ math }
34 changes: 34 additions & 0 deletions man/is1.Rd
@@ -0,0 +1,34 @@
\name{is1}
\alias{is1}
\title{Compute IS1 formula from Nason (2005) }
\description{
Computes $I_{S1}(p,d) = \int_{d}^{\infty} \cos (px) \exp (-x^2) dx$}
}
\usage{
is1(p, d)
}
%- maybe also 'usage' for other objects documented here.
\arguments{
\item{p}{ Argument for IS1 function}
\item{d}{ Argument for IS1 function}
}
\details{
An intermediate function for the computation of the inner product
between Gaussian and Student's t distribution.
}
\value{
Value of the function IS1
}
\references{ Nason, G.P. (2005) On the sum of the Gaussian and Student's t
random variables. \emph{Technical Report}, 05:01,
Statistics Group, Department of Mathematics, University of Bristol.}
\author{ Guy Nason }

\seealso{ \code{\link{normt3ip}}}
\examples{
is1(1,1)
#
# Should give 0.1297382.
#
}
\keyword{ math }
34 changes: 34 additions & 0 deletions man/normt3ip.Rd
@@ -0,0 +1,34 @@
\name{normt3ip}
\alias{normt3ip}
\title{Compute inner product function of Gaussian and t3 distribution}
\description{This function evaluates the inner product function of
the (sphered) Student's t distribution on 3 df and the Gaussian
distribution as defined by Theorem~1 of Nason, 2005
}
\usage{
normt3ip(mu, sigma)
}
\arguments{
\item{mu}{ Mean of the Gaussian}
\item{sigma}{ Standard deviation of the Gaussian}
}
\details{No extra details }
\value{
The evaluated function.
}
\references{ Nason, G.P. (2005) On the sum of the Gaussian and Student's t
random variables. \emph{Technical Report}, 05:01,
Statistics Group, Department of Mathematics, University of Bristol.}
\author{ Guy Nason }


%\seealso{\code{\link{ic1}},\code{\link{is1}},\code{\link{dnormt3}}}
\examples{
#normt3ip(0,1)
#
# Answer should be 0.3183099 which is 1/pi
#
x <- 3
}
\keyword{ math }

0 comments on commit 093aa6f

Please sign in to comment.