[Yuima-commits] r591 - in pkg/yuima: . R src
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu Feb 23 03:13:39 CET 2017
Author: kyuta
Date: 2017-02-23 03:13:38 +0100 (Thu, 23 Feb 2017)
New Revision: 591
Removed:
pkg/yuima/src/euler.c
Modified:
pkg/yuima/DESCRIPTION
pkg/yuima/NEWS
pkg/yuima/R/sim.euler.R
Log:
modified sim.euler.R and removed euler.c due to a memory corruption
fixed a bug in sim.euler.R
Modified: pkg/yuima/DESCRIPTION
===================================================================
--- pkg/yuima/DESCRIPTION 2017-02-22 10:28:54 UTC (rev 590)
+++ pkg/yuima/DESCRIPTION 2017-02-23 02:13:38 UTC (rev 591)
@@ -1,7 +1,7 @@
Package: yuima
Type: Package
Title: The YUIMA Project Package for SDEs
-Version: 1.5.9
+Version: 1.6.0
Depends: R(>= 2.10.0), methods, zoo, stats4, utils, expm, cubature, mvtnorm
Imports: Rcpp (>= 0.12.1)
Author: YUIMA Project Team
Modified: pkg/yuima/NEWS
===================================================================
--- pkg/yuima/NEWS 2017-02-22 10:28:54 UTC (rev 590)
+++ pkg/yuima/NEWS 2017-02-23 02:13:38 UTC (rev 591)
@@ -43,4 +43,6 @@
2016/07/08: fixed some bugs in llag.R and cce_functions.c
2016/10/04: modified setMultiModel.R, sim.euler.R and yuima.model.R to generate nts and pts process
2016/12/16: added rGIG, rGH, dGIG and dGH in rng.R and the corresponding c language file YU
-2017/01/25: modified sim.euler.R and added euler.c to implement the Euler-Maruyama scheme by the C code (only the diffusion case)
\ No newline at end of file
+2017/01/25: modified sim.euler.R and added euler.c to implement the Euler-Maruyama scheme by the C code (only the diffusion case)
+2017/02/23: modified sim.euler.R and removed euler.c due to a memory corruption
+ fix a bug in sim.euler.R
\ No newline at end of file
Modified: pkg/yuima/R/sim.euler.R
===================================================================
--- pkg/yuima/R/sim.euler.R 2017-02-22 10:28:54 UTC (rev 590)
+++ pkg/yuima/R/sim.euler.R 2017-02-23 02:13:38 UTC (rev 591)
@@ -102,12 +102,13 @@
if(!length(sdeModel at measure.type)){ ##:: Wiener Proc
##:: using Euler-Maruyama method
- if(0){ # old version (Jan 25, 2017)
+ #if(0){ # old version (Jan 25, 2017)
if(var.in.diff & (!is.Poisson(sdeModel))){ ##:: diffusions have state variables and it is not Poisson
##:: calcurate difference eq.
for( i in 1:n){
# dX <- dX + p.b(t=i*delta, X=dX) %*% dW[, i]
- dX <- dX + p.b(t=yuima at sampling@Initial+i*delta, X=dX) %*% dW[, i] # LM
+ #dX <- dX + p.b(t=yuima at sampling@Initial+i*delta, X=dX) %*% dW[, i] # LM
+ dX <- dX + p.b(t=yuima at sampling@Initial+(i-1)*delta, X=dX) %*% dW[, i] # YK
X_mat[,i+1] <- dX
}
}else{ ##:: diffusions have no state variables (not use p.b(). faster)
@@ -154,8 +155,9 @@
X_mat[, i+1] <- dX
}
}
- }
+ #}
+ if(0){ # currently ignored due to a bug (YK, Feb 23, 2017)
# new version (Jan 25, 2017)
b <- parse(text=paste("c(",paste(as.character(V0),collapse=","),")"))
vecV <- parse(text=paste("c(",paste(as.character(unlist(V)),collapse=","),")"))
@@ -163,7 +165,7 @@
X_mat <- .Call("euler", dX, Initial, as.integer(r.size),
rep(1, n) * delta, dW, modeltime, modelstate, quote(eval(b, env)),
quote(eval(vecV, env)), env, new.env())
-
+ }
#tsX <- ts(data=t(X_mat), deltat=delta , start=0)
tsX <- ts(data=t(X_mat), deltat=delta , start = yuima at sampling@Initial) #LM
}else{ ##:: Levy
@@ -286,15 +288,18 @@
}
# Pi <- p.b.j(t=i*delta,X=dX) * J #LM
- Pi <- p.b.j(t=yuima at sampling@Initial+i*delta,X=dX) * J
+ #Pi <- p.b.j(t=yuima at sampling@Initial+i*delta,X=dX) * J
+ Pi <- p.b.j(t=yuima at sampling@Initial+(i-1)*delta,X=dX) * J # YK
}else{# we add this part since we allow the user to specify the increment of CP LM 05/02/2015
# Pi <- p.b.j(t=i*delta,X=dX) %*% dL[, i] #LM
- Pi <- p.b.j(t=yuima at sampling@Initial+i*delta,X=dX) %*% dL[, i]
+ #Pi <- p.b.j(t=yuima at sampling@Initial+i*delta,X=dX) %*% dL[, i]
+ Pi <- p.b.j(t=yuima at sampling@Initial+(i - 1)*delta,X=dX) %*% dL[, i] # YK
}
##Pi <- p.b.j(t=i*delta, X=dX)
}
# dX <- dX + p.b(t=i*delta, X=dX) %*% dW[, i] + Pi # LM
- dX <- dX + p.b(t=yuima at sampling@Initial + i*delta, X=dX) %*% dW[, i] + Pi
+ #dX <- dX + p.b(t=yuima at sampling@Initial + i*delta, X=dX) %*% dW[, i] + Pi
+ dX <- dX + p.b(t=yuima at sampling@Initial + (i - 1)*delta, X=dX) %*% dW[, i] + Pi # YK
X_mat[, i+1] <- dX
}
# tsX <- ts(data=t(X_mat), deltat=delta, start=0) #LM
@@ -357,7 +362,8 @@
}
# cat("\np.b.j call\n")
# tmp.j <- p.b.j(t=i*delta, X=dX) #LM
- tmp.j <- p.b.j(t=yuima at sampling@Initial+i*delta, X=dX)
+ #tmp.j <- p.b.j(t=yuima at sampling@Initial+i*delta, X=dX)
+ tmp.j <- p.b.j(t=yuima at sampling@Initial+(i - 1)*delta, X=dX) # YK
#print(str(tmp.j))
#cat("\np.b.j cback and dZ\n")
# print(str(dZ[,i]))
@@ -367,7 +373,8 @@
#print(str(tmp.j))
#print(str(p.b(t = i * delta, X = dX) %*% dW[, i]))
# dX <- dX + p.b(t=i*delta, X=dX) %*% dW[, i] +tmp.j %*% dZ[,i] #LM
- dX <- dX + p.b(t=yuima at sampling@Initial+i*delta, X=dX) %*% dW[, i] +tmp.j %*% dZ[,i]
+ #dX <- dX + p.b(t=yuima at sampling@Initial+i*delta, X=dX) %*% dW[, i] +tmp.j %*% dZ[,i]
+ dX <- dX + p.b(t=yuima at sampling@Initial+(i - 1)*delta, X=dX) %*% dW[, i] +tmp.j %*% dZ[,i] # YK
X_mat[, i+1] <- dX
}
# tsX <- ts(data=t(X_mat), deltat=delta, start=0) #LM
Deleted: pkg/yuima/src/euler.c
===================================================================
--- pkg/yuima/src/euler.c 2017-02-22 10:28:54 UTC (rev 590)
+++ pkg/yuima/src/euler.c 2017-02-23 02:13:38 UTC (rev 591)
@@ -1,83 +0,0 @@
-//
-// euler.c
-//
-//
-// Created by Yuta Koike on 2016/01/23.
-//
-//
-
-#include <R.h>
-#include <Rmath.h>
-#include <Rdefines.h>
-#include <Rinternals.h>
-
-SEXP euler(SEXP x0, SEXP t0, SEXP R, SEXP dt, SEXP dW, SEXP modeltime, SEXP modelstate,
- SEXP drift, SEXP diffusion, SEXP env, SEXP rho){
-
- int i, j, k, n, d, r;
- double *rdt, *rdW, *rX, *rx0, *b, *sigma;
- SEXP X, xpar, tpar;
-
- PROTECT(x0 = AS_NUMERIC(x0));
- rx0 = REAL(x0);
- PROTECT(dt = AS_NUMERIC(dt));
- rdt = REAL(dt);
- PROTECT(dW = AS_NUMERIC(dW));
- rdW = REAL(dW);
-
- d = length(x0);
- n = length(dt);
-
- r = *INTEGER(R);
-
- /* PROTECT(X = allocVector(REALSXP, n + 1));
- rX = REAL(X);
- rX[0] = *REAL(x0);*/
- PROTECT(X = allocMatrix(REALSXP, d, n + 1));
- rX = REAL(X);
-
- for (j = 0; j < d; j++) {
- rX[j] = rx0[j];
- }
-
- PROTECT(tpar = allocVector(REALSXP, 1));
- PROTECT(xpar = allocVector(REALSXP, 1));
-
- PROTECT(t0 = AS_NUMERIC(t0));
- REAL(tpar)[0] = REAL(t0)[0]; /* initial time */
-
- for (i = 0; i < n; i++) {
-
- /* assign the current variables */
- defineVar(installChar(STRING_ELT(modeltime, 0)), tpar, env);
-
- for (j = 0; j < d; j++) {
- /*xpar = allocVector(REALSXP, 1);*/
- REAL(xpar)[0] = rX[j + i * d];
- defineVar(installChar(STRING_ELT(modelstate, j)), duplicate(xpar), env);
- /*defineVar(installChar(STRING_ELT(modelstate, j)), xpar, env);*/
- }
-
- /*defineVar(install("env"), env, rho);*/
-
- /* evaluate coefficients */
- b = REAL(eval(drift, rho));
- sigma = REAL(eval(diffusion, rho));
-
- for (j = 0; j < d; j++) {
- rX[j + (i + 1) * d] = rX[j + i * d] + b[j] * rdt[i];
- for (k = 0; k < r; k++) {
- rX[j + (i + 1) * d] += sigma[k + j * r] * rdW[k + i * r];
- }
- }
-
- /*defineVar(install("t"), tpar, rho);*/
- /*rX[i + 1] = rX[i] + *REAL(eval(drift, rho)) * rdt[i] + *REAL(eval(diffusion, rho)) * rdW[i];*/
- /*rX[i + 1] = rX[i] + *REAL(eval(drift, rho)) * REAL(dt)[i] + *REAL(eval(diffusion, rho)) * REAL(dW)[i];*/
-
- REAL(tpar)[0] += rdt[i];
- }
-
- UNPROTECT(7);
- return(X);
-}
More information about the Yuima-commits
mailing list