[Yuima-commits] r56 - in pkg/yuima: . R
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Fri Jan 8 04:27:39 CET 2010
Author: hirokimasuda
Date: 2010-01-08 04:27:38 +0100 (Fri, 08 Jan 2010)
New Revision: 56
Modified:
pkg/yuima/DESCRIPTION
pkg/yuima/R/rng.R
Log:
rng file slightly updated (masuda)
Modified: pkg/yuima/DESCRIPTION
===================================================================
--- pkg/yuima/DESCRIPTION 2009-12-29 01:32:41 UTC (rev 55)
+++ pkg/yuima/DESCRIPTION 2010-01-08 03:27:38 UTC (rev 56)
@@ -1,12 +1,12 @@
-Package: yuima
-Type: Package
-Title: The YUIMA Project package
-Version: 0.0.84
-Date: 2010-12-28
-Depends: methods, zoo, adapt
-Author: YUIMA Project Team.
-Maintainer: Stefano M. Iacus <stefano.iacus at R-project.org>
-Description: The YUIMA Project for Simulation and Inference
- for Stochastic Differential Equations.
-License: GPL-2
-URL: http://R-Forge.R-project.org/projects/yuima/
+Package: yuima
+Type: Package
+Title: The YUIMA Project package
+Version: 0.0.85
+Date: 2010-1-8
+Depends: methods, zoo, adapt
+Author: YUIMA Project Team.
+Maintainer: Stefano M. Iacus <stefano.iacus at R-project.org>
+Description: The YUIMA Project for Simulation and Inference
+ for Stochastic Differential Equations.
+License: GPL-2
+URL: http://R-Forge.R-project.org/projects/yuima/
Modified: pkg/yuima/R/rng.R
===================================================================
--- pkg/yuima/R/rng.R 2009-12-29 01:32:41 UTC (rev 55)
+++ pkg/yuima/R/rng.R 2010-01-08 03:27:38 UTC (rev 56)
@@ -1,84 +1,12 @@
-## random number generators for yuima packages
+##################################################################
+###### "Random number generators and
+###### related density functions for yuima packages"
+###### (Last modified: Jan. 8, 2010)
+##################################################################
-rIG <- function(x,delta,gamma){
- if( delta <= 0 )
- stop("delta must be positive value.")
- if( gamma <= 0 )
- stop("deltat must be positive value.")
- V <- rchisq(x,df=1);
-
- z1 <- ( delta/gamma + V/(2*gamma^2) ) - sqrt( V*delta/(gamma^3) + ( V/(2*gamma^2) )^2 )
- U <- runif(x,min=0,max=1)
- idx <- which( U < (delta/(delta+gamma*z1)) )
- z2 <- (delta/gamma)^2 /z1[-idx]
- ret <- numeric(x)
- ret[idx] <- z1[idx]
- ret[-idx] <- z2
-
- return(ret)
-}
+## Bilateral gamma
-## multi-variate normal inverse Gaussian
-rNIG <- function(x,alpha=1,beta=0,delta=1,mu=0,Lambda){
- ## Error check
- if(length(mu)!=length(beta)){
- stop("Error: wrong input dimension.")
- }
-
- if(missing(Lambda)){
- ## univariate case
- gamma <- sqrt(alpha^2 - beta^2)
-
- if (gamma == 0) {
- V = rnorm(x)^2
- Z = delta * delta/V
- X = sqrt(Z) * rnorm(x)
- }else{
- Z <- rIG(x,delta,gamma)
- N <- rnorm(x,0,1)
- X <- mu + beta*Z + sqrt(Z)*N
- }
- return(X)
-
- }else{ ## multivariate case
- if( det(Lambda) < 10^(-15)){
- stop("Determinant of Lambda must be one.")
- }
- tmp <- as.numeric(alpha^2 - t(beta) %*% Lambda %*% beta)
- if(tmp <0){
- stop("gamma must be positive.")
- }
- gamma <- sqrt(tmp)
- }
-
- tau <- rIG(x,delta,gamma)
- eta <- rnorm(x*length(beta))
- sqrt.L <- svd(Lambda)
- sqrt.L <- sqrt.L$u %*% diag(sqrt(sqrt.L$d)) %*% t(sqrt.L$v)
- z <- mu + t(matrix(rep(tau,length(beta)),x,length(beta))) * matrix(rep(Lambda %*% beta,x),length(beta),x)+t(matrix(rep(tau,length(beta)),x,length(beta))) * (Lambda %*% t(matrix(eta,x,length(beta))))
-
- X <- z
- return(X)
-}
-
-
-## rNIG <- function(x,alpha=1,beta=0,delta=1,mu=0){
-## gamma <- sqrt(alpha^2 - beta^2)
-## if (gamma == 0) {
-## V = rnorm(x)^2
-## Z = delta * delta/V
-## X = sqrt(Z) * rnorm(x)
-## }else{
-## Z <- rIG(x,delta,gamma)
-## N <- rnorm(x,0,1)
-## X <- mu + beta*Z + sqrt(Z)*N
-## }
-## return(X)
-## }
-
-
-
rbgamma <- function(x,delta.plus,gamma.plus,delta.minus,gamma.minus){
if( delta.plus <= 0 )
stop("delta.plus must be positive value.")
@@ -93,7 +21,8 @@
}
-## temporaly, desined only for univariate.
+## (Multivariate) Normal gamma
+
rngamma <- function(x,lambda,alpha,beta,mu,Lambda){
## Error check
if(length(mu)!=length(beta)){
@@ -147,7 +76,20 @@
}
-## ## temporaly, desined only for univariate.
+dungamma <- function(x,lam,al,be,mu){
+ if( lam <= 0 )
+ stop("lambda must be positive value.")
+ if( al <= 0 )
+ stop("alpha must be positive value.")
+ if( al^2 - be^2 <= 0 )
+ stop("alpha^2 - beta^2 must be positive value.")
+
+ dens <- exp(be*(x-mu))*((al^2 - be^2)^(lam))*besselK(al*abs(x-mu),lam-1/2)*abs(x-mu)^(lam-1/2)/(gamma(lam)*sqrt(pi)*((2*al)^(lam-1/2)))
+ dens
+}
+
+
+## ## One-dim. normal gamma
## rngamma <- function(x,lambda,alpha,beta,mu){
## tmp <- alpha^2 - beta^2
## if( lambda <= 0 )
@@ -188,7 +130,114 @@
## }
-rstable <- function(x,alpha,beta,sigma, gamma){
+## Inverse Gaussian
+
+rIG <- function(x,delta,gamma){
+ if( delta <= 0 )
+ stop("delta must be positive value.")
+ if( gamma <= 0 )
+ stop("gamma must be positive value.")
+ V <- rchisq(x,df=1);
+
+ z1 <- ( delta/gamma + V/(2*gamma^2) ) - sqrt( V*delta/(gamma^3) + ( V/(2*gamma^2) )^2 )
+ U <- runif(x,min=0,max=1)
+ idx <- which( U < (delta/(delta+gamma*z1)) )
+ z2 <- (delta/gamma)^2 /z1[-idx]
+ ret <- numeric(x)
+ ret[idx] <- z1[idx]
+ ret[-idx] <- z2
+
+ return(ret)
+}
+
+
+dIG <- function(x,delta,gamma){
+ if( delta <= 0 )
+ stop("delta must be positive value.")
+ if( gamma <= 0 )
+ stop("gamma must be positive value.")
+
+ dens <- delta*e^(delta*gamma)*(x^(-3/2))*exp(-((delta^2)/x+x*gamma^2)/2)/sqrt(2*pi)
+ dens
+}
+
+
+## (Multivariate) Normal inverse Gaussian
+
+rNIG <- function(x,alpha=1,beta=0,delta=1,mu=0,Lambda){
+ ## Error check
+ if(length(mu)!=length(beta)){
+ stop("Error: wrong input dimension.")
+ }
+
+ if(missing(Lambda)){
+ ## univariate case
+ gamma <- sqrt(alpha^2 - beta^2)
+
+ if (gamma == 0) {
+ V = rnorm(x)^2
+ Z = delta * delta/V
+ X = sqrt(Z) * rnorm(x)
+ }else{
+ Z <- rIG(x,delta,gamma)
+ N <- rnorm(x,0,1)
+ X <- mu + beta*Z + sqrt(Z)*N
+ }
+ return(X)
+
+ }else{ ## multivariate case
+ if( det(Lambda) < 10^(-15)){
+ stop("Determinant of Lambda must be one.")
+ }
+ tmp <- as.numeric(alpha^2 - t(beta) %*% Lambda %*% beta)
+ if(tmp <0){
+ stop("gamma must be positive.")
+ }
+ gamma <- sqrt(tmp)
+ }
+
+ tau <- rIG(x,delta,gamma)
+ eta <- rnorm(x*length(beta))
+ sqrt.L <- svd(Lambda)
+ sqrt.L <- sqrt.L$u %*% diag(sqrt(sqrt.L$d)) %*% t(sqrt.L$v)
+ z <- mu + t(matrix(rep(tau,length(beta)),x,length(beta))) * matrix(rep(Lambda %*% beta,x),length(beta),x)+t(matrix(rep(tau,length(beta)),x,length(beta))) * (Lambda %*% t(matrix(eta,x,length(beta))))
+ X <- z
+ return(X)
+}
+
+
+duNIG <- function(x,al,be,de,mu){
+ if( al < 0 )
+ stop("alpha must be nonnegative.")
+ if( al^2 - be^2 < 0 )
+ stop("alpha^2 - beta^2 must be nonnegative.")
+ if( de <= 0 )
+ stop("delta must be positive.")
+
+ dens <- al*de*exp(de*sqrt(al^{2}-be^{2})+be*(x-mu))*besselK(al*sqrt(de^{2}+(x-mu)^{2}),1)/(pi*sqrt(de^{2}+(x-mu)^{2}))
+ dens
+}
+
+
+## ## One-dim. NIG
+## rNIG <- function(x,alpha=1,beta=0,delta=1,mu=0){
+## gamma <- sqrt(alpha^2 - beta^2)
+## if (gamma == 0) {
+## V = rnorm(x)^2
+## Z = delta * delta/V
+## X = sqrt(Z) * rnorm(x)
+## }else{
+## Z <- rIG(x,delta,gamma)
+## N <- rnorm(x,0,1)
+## X <- mu + beta*Z + sqrt(Z)*N
+## }
+## return(X)
+## }
+
+
+## Non-Gaussian stable
+
+rstable <- function(x,alpha,beta,sigma,gamma){
a <- (1 + (beta*tan(alpha*pi/2))^2)^(1/(2*alpha))
b <- atan(beta*tan(alpha*pi/2))/alpha
@@ -197,12 +246,12 @@
if(alpha!=1){
s <- a * (sin(alpha*(u+b))/cos(u)^(1/alpha)) * (cos(u-alpha*(u+b))/v)^((1-alpha)/alpha)
+ X <- sigma * s + gamma * rep(1,x)
}else{
s <- (2/pi) * ((pi/2 +beta*u)*tan(u) - beta * log(v*cos(u)/(beta*u + pi/2)))
+ X <- sigma * s + (gamma + (2/pi)*beta*sigma*log(sigma)) * rep(1,x)
}
-## X <- (eps^(1/alpha)) * sigma * s + gamma * eps * rep(1,x)
- X <- sigma * s + gamma * rep(1,x)
return(X)
}
@@ -222,3 +271,31 @@
## X <- (eps^(1/alpha)) * sigma * s + gamma * eps * rep(1,x)
## return(X)
## }
+
+
+## Positive exponentially tempered stable (AR method)
+## ## This must be re-coded later!! (temporarily, rather inefficient)
+rpts <- function(x,al,a,b) {
+ if( al <= 0 | al>= 1 )
+ stop("al must lie in (0,1).")
+ if( a <= 0 )
+ stop("a must be positive value.")
+ if( b <= 0 )
+ stop("b must be positive value.")
+
+ sig <- (-a*gamma(-al)*cos(pi*al/2))^(1/al)
+ y <- c()
+ i <- 1
+ while (i <= x) {
+ u <- runif(1)
+ v <- rstable(1,al,1,sig,0)
+ w <- exp(-b*v)
+ if (u < w ){
+ y <- append(y,v)
+ i <- i+1
+ }
+ }
+ return(y)
+}
+
+
More information about the Yuima-commits
mailing list