[Yuima-commits] r454 - in pkg/yuima: R src
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Fri Jul 8 05:18:38 CEST 2016
Author: kamatani
Date: 2016-07-08 05:18:37 +0200 (Fri, 08 Jul 2016)
New Revision: 454
Added:
pkg/yuima/R/RcppExports.R
pkg/yuima/src/RcppExports.cpp
pkg/yuima/src/qmlecpp.cpp
Modified:
pkg/yuima/R/qmle.R
Log:
add rcpp option to qmle
Added: pkg/yuima/R/RcppExports.R
===================================================================
--- pkg/yuima/R/RcppExports.R (rev 0)
+++ pkg/yuima/R/RcppExports.R 2016-07-08 03:18:37 UTC (rev 454)
@@ -0,0 +1,23 @@
+# This file was generated by Rcpp::compileAttributes
+# Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393
+
+detcpp <- function(A) {
+ .Call('yuima_detcpp', PACKAGE = 'yuima', A)
+}
+
+Smake <- function(b, d) {
+ .Call('yuima_Smake', PACKAGE = 'yuima', b, d)
+}
+
+solvecpp <- function(A) {
+ .Call('yuima_solvecpp', PACKAGE = 'yuima', A)
+}
+
+trace <- function(S, b) {
+ .Call('yuima_trace', PACKAGE = 'yuima', S, b)
+}
+
+likndim <- function(dx, b, A, h) {
+ .Call('yuima_likndim', PACKAGE = 'yuima', dx, b, A, h)
+}
+
Modified: pkg/yuima/R/qmle.R
===================================================================
--- pkg/yuima/R/qmle.R 2016-07-08 02:20:05 UTC (rev 453)
+++ pkg/yuima/R/qmle.R 2016-07-08 03:18:37 UTC (rev 454)
@@ -122,7 +122,7 @@
}
qmle <- function(yuima, start, method="BFGS", fixed = list(), print=FALSE,
- lower, upper, joint=FALSE, Est.Incr="NoIncr",aggregation=TRUE, threshold=NULL, ...){
+ lower, upper, joint=FALSE, Est.Incr="NoIncr",aggregation=TRUE, threshold=NULL,rcpp=FALSE, ...){
if(is(yuima at model, "yuima.carma")){
NoNeg.Noise<-FALSE
cat("\nStarting qmle for carma ... \n")
@@ -480,7 +480,6 @@
assign("measure.var", args[1], envir=env)
}
-
f <- function(p) {
mycoef <- as.list(p)
if(!is.CARMA(yuima)){
@@ -495,7 +494,7 @@
names(mycoef) <- nm
}
mycoef[fixed.par] <- fixed
- minusquasilogl(yuima=yuima, param=mycoef, print=print, env)
+ minusquasilogl(yuima=yuima, param=mycoef, print=print, env,rcpp=rcpp)
}
# SMI-2/9/14:
@@ -528,7 +527,7 @@
mycoef[fixed.par] <- fixed
}
mycoef[fixed.par] <- fixed
- minusquasilogl(yuima=yuima, param=mycoef, print=print, env)
+ minusquasilogl(yuima=yuima, param=mycoef, print=print, env,rcpp=rcpp)
}
oout <- NULL
@@ -765,7 +764,7 @@
names(mycoef) <- drift.par
mycoef[diff.par] <- coef[diff.par]
}
- minusquasilogl(yuima=yuima, param=mycoef, print=print, env)
+ minusquasilogl(yuima=yuima, param=mycoef, print=print, env,rcpp=rcpp)
}
fDiff <- function(p) {
@@ -774,7 +773,7 @@
names(mycoef) <- diff.par
mycoef[drift.par] <- coef[drift.par]
}
- minusquasilogl(yuima=yuima, param=mycoef, print=print, env)
+ minusquasilogl(yuima=yuima, param=mycoef, print=print, env,rcpp=rcpp)
}
# coef <- oout$par
@@ -963,10 +962,10 @@
if(length(c(diff.par,drift.par))>0 & !is.CARMA(yuima)){ # LM 04/09/14
- min.diff <- minusquasilogl(yuima=yuima, param=mycoef[c(diff.par,drift.par)], print=print, env)
+ min.diff <- minusquasilogl(yuima=yuima, param=mycoef[c(diff.par,drift.par)], print=print, env,rcpp=rcpp)
}else{
if(length(c(diff.par,drift.par))>0 & is.CARMA(yuima)){
- min.diff <- minusquasilogl(yuima=yuima, param=mycoef, print=print, env)
+ min.diff <- minusquasilogl(yuima=yuima, param=mycoef, print=print, env,rcpp=rcpp)
}
}
@@ -1499,7 +1498,7 @@
}
-quasilogl <- function(yuima, param, print=FALSE){
+quasilogl <- function(yuima, param, print=FALSE,rcpp=FALSE){
d.size <- yuima at model@equation.number
if (is(yuima at model, "yuima.carma")){
@@ -1530,11 +1529,11 @@
assign("h", deltat(yuima at data@zoo.data[[1]]), envir=env)
assign("time", as.numeric(index(yuima at data@zoo.data[[1]])), envir=env)
- -minusquasilogl(yuima=yuima, param=param, print=print, env)
+ -minusquasilogl(yuima=yuima, param=param, print=print, env,rcpp=rcpp)
}
-minusquasilogl <- function(yuima, param, print=FALSE, env){
+minusquasilogl <- function(yuima, param, print=FALSE, env,rcpp=FALSE){
diff.par <- yuima at model@parameter at diffusion
@@ -1758,7 +1757,7 @@
# yuima.warn("carma(p,q): the scale parameter is equal to 1. We will implemented as soon as possible")
# return(NULL)
# }
- } else{
+ } else if (!rcpp) {
drift <- drift.term(yuima, param, env)
diff <- diffusion.term(yuima, param, env)
@@ -1796,7 +1795,59 @@
}
}
}
- }
+ } else {
+ b <- yuima at model@drift
+ a <- yuima at model@diffusion
+ d <- d.size
+ ####data <- yuima at data@original.data
+ data <- matrix(0,length(yuima at data@zoo.data[[1]]),d.size)
+ for(i in 1:d) data[,i] <- as.numeric(yuima at data@zoo.data[[i]])
+ ####delta <- yuima at sampling@delta
+ delta <- deltat(yuima at data@zoo.data[[1]])
+ thetadim <- length(yuima at model@parameter at all)
+ ####r <- length(a[[1]])
+ r <- yuima at model@noise.number
+ xdim <- length(yuima at model@state.variable)
+
+ #if(thetadim!=length(initial)) stop("check dim of initial") #error check
+
+ d_b <- NULL
+ for(i in 1:d){
+ check_x <- NULL
+ for(k in 1:xdim) check_x <- cbind(check_x,length(grep(yuima at model@state.variable[k],b[[i]])))
+ if(sum(check_x)>0){
+ d_b[[i]] <- b[[i]] #this part of model includes "x"(state.variable)
+ }
+ else{
+ d_b[[i]] <- parse(text=paste("(",b[[i]][2],")*rep(1,length(data[,1])-1)",sep="")) #vectorization
+ }
+ }
+ #d_b <- c(d_b,b[[i]])
+
+ v_a<-matrix(list(NULL),d,r)
+ for(i in 1:d){
+ for(j in 1:r){
+ check_x <- NULL
+ for(k in 1:xdim) check_x <- cbind(check_x,length(grep(yuima at model@state.variable[k],a[[i]][[j]])))
+ if(sum(check_x)>0){
+ v_a[[i,j]] <- a[[i]][[j]] #this part of model includes "x"(state.variable)
+ }
+ else{
+ v_a[[i,j]] <- parse(text=paste("(",a[[i]][[j]][2],")*rep(1,length(data[,1])-1)",sep="")) #vectorization
+ }
+ }
+ }
+
+ for(i in 1:d) assign(yuima at model@state.variable[i], data[-length(data[,1]),i])
+ dx <- as.matrix((data-rbind(numeric(d),as.matrix(data[-length(data[,1]),])))[-1,])
+ drift <- diffusion <- NULL
+ for(i in 1:thetadim) assign(names(param)[i], param[[i]])
+ for(i in 1:d) drift <- cbind(drift,eval(d_b[[i]]))
+ for(i in 1:r){
+ for(j in 1:d) diffusion <- cbind(diffusion,eval(v_a[[j,i]]))
+ }
+ QL <- (likndim(dx,drift,diffusion,delta)*(-0.5) + (n-1)*(-0.5*d.size * log( (2*pi*h) )))
+ }
if(!is.finite(QL)){
Added: pkg/yuima/src/RcppExports.cpp
===================================================================
--- pkg/yuima/src/RcppExports.cpp (rev 0)
+++ pkg/yuima/src/RcppExports.cpp 2016-07-08 03:18:37 UTC (rev 454)
@@ -0,0 +1,67 @@
+// This file was generated by Rcpp::compileAttributes
+// Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393
+
+#include <Rcpp.h>
+
+using namespace Rcpp;
+
+// detcpp
+double detcpp(NumericMatrix A);
+RcppExport SEXP yuima_detcpp(SEXP ASEXP) {
+BEGIN_RCPP
+ Rcpp::RObject __result;
+ Rcpp::RNGScope __rngScope;
+ Rcpp::traits::input_parameter< NumericMatrix >::type A(ASEXP);
+ __result = Rcpp::wrap(detcpp(A));
+ return __result;
+END_RCPP
+}
+// Smake
+NumericMatrix Smake(NumericVector b, int d);
+RcppExport SEXP yuima_Smake(SEXP bSEXP, SEXP dSEXP) {
+BEGIN_RCPP
+ Rcpp::RObject __result;
+ Rcpp::RNGScope __rngScope;
+ Rcpp::traits::input_parameter< NumericVector >::type b(bSEXP);
+ Rcpp::traits::input_parameter< int >::type d(dSEXP);
+ __result = Rcpp::wrap(Smake(b, d));
+ return __result;
+END_RCPP
+}
+// solvecpp
+NumericMatrix solvecpp(NumericMatrix A);
+RcppExport SEXP yuima_solvecpp(SEXP ASEXP) {
+BEGIN_RCPP
+ Rcpp::RObject __result;
+ Rcpp::RNGScope __rngScope;
+ Rcpp::traits::input_parameter< NumericMatrix >::type A(ASEXP);
+ __result = Rcpp::wrap(solvecpp(A));
+ return __result;
+END_RCPP
+}
+// trace
+double trace(NumericMatrix S, NumericVector b);
+RcppExport SEXP yuima_trace(SEXP SSEXP, SEXP bSEXP) {
+BEGIN_RCPP
+ Rcpp::RObject __result;
+ Rcpp::RNGScope __rngScope;
+ Rcpp::traits::input_parameter< NumericMatrix >::type S(SSEXP);
+ Rcpp::traits::input_parameter< NumericVector >::type b(bSEXP);
+ __result = Rcpp::wrap(trace(S, b));
+ return __result;
+END_RCPP
+}
+// likndim
+double likndim(NumericMatrix dx, NumericMatrix b, NumericMatrix A, double h);
+RcppExport SEXP yuima_likndim(SEXP dxSEXP, SEXP bSEXP, SEXP ASEXP, SEXP hSEXP) {
+BEGIN_RCPP
+ Rcpp::RObject __result;
+ Rcpp::RNGScope __rngScope;
+ Rcpp::traits::input_parameter< NumericMatrix >::type dx(dxSEXP);
+ Rcpp::traits::input_parameter< NumericMatrix >::type b(bSEXP);
+ Rcpp::traits::input_parameter< NumericMatrix >::type A(ASEXP);
+ Rcpp::traits::input_parameter< double >::type h(hSEXP);
+ __result = Rcpp::wrap(likndim(dx, b, A, h));
+ return __result;
+END_RCPP
+}
Added: pkg/yuima/src/qmlecpp.cpp
===================================================================
--- pkg/yuima/src/qmlecpp.cpp (rev 0)
+++ pkg/yuima/src/qmlecpp.cpp 2016-07-08 03:18:37 UTC (rev 454)
@@ -0,0 +1,97 @@
+#include <Rcpp.h>
+using namespace Rcpp;
+
+// Below is a simple example of exporting a C++ function to R. You can
+// source this function into an R session using the Rcpp::sourceCpp
+// function (or via the Source button on the editor toolbar)
+
+// For more on using Rcpp click the Help button on the editor toolbar
+
+// [[Rcpp::export]]
+double detcpp(NumericMatrix A){ //det(A)
+ int n = A.nrow();
+ double det = 1.0,buf;
+ NumericMatrix B = clone(A);
+
+ for(int i=0;i<n;i++){
+ buf = 1/B(i,i);
+ for(int j=i+1;j<n;j++){
+ for(int k=i+1;k<n;k++){
+ B(j,k)-=B(i,k)*B(j,i)*buf;
+ }
+ }
+ det*=B(i,i);
+ }
+ return det;
+}
+
+
+// [[Rcpp::export]]
+NumericMatrix Smake(NumericVector b,int d){ //tcrossprod(matrix(b,d,r))
+ int r = b.length()/d;
+ NumericMatrix S(d,d);
+
+ for(int i=0;i<d;i++){
+ for(int j=0;j<d;j++){
+ for(int k=0;k<r;k++){
+ S(i,j) += b(d*k+i)*b(d*k+j);
+ }
+ }
+ }
+ return S;
+}
+// [[Rcpp::export]]
+NumericMatrix solvecpp(NumericMatrix A){ //solve(A)
+ int n = A.ncol();
+ double buf;
+ NumericMatrix B = clone(A);
+ NumericMatrix inv(n,n);
+
+ for(int i=0;i<n;i++){
+ inv(i,i) = 1.0;
+ buf = 1.0/B(i,i);
+ for(int j=0;j<n;j++){
+ if(j<i+1) inv(i,j)*=buf;
+ else B(i,j)*=buf;
+ }
+ for(int j=0;j<n;j++){
+ if(i!=j){
+ buf=B(j,i);
+ for(int k=0;k<n;k++){
+ if(k<i+1) inv(j,k)-=inv(i,k)*buf;
+ else B(j,k)-=B(i,k)*buf;
+ }
+ }
+ }
+ }
+ return inv;
+}
+
+// [[Rcpp::export]]
+double trace(NumericMatrix S,NumericVector b){ // tr(S%*%tcrossprod(b))
+ int n = S.nrow();
+ double tr = 0;
+
+ for(int i=0;i<n;i++){
+ for(int k=0;k<i;k++){
+ tr += S(i,k)*b(k)*b(i);
+ }
+ }
+ tr*=2;
+ for(int i=0;i<n;i++) tr += S(i,i)*b(i)*b(i);
+ return tr;
+}
+// [[Rcpp::export]]
+double likndim(NumericMatrix dx,NumericMatrix b,NumericMatrix A,double h){
+ int n = dx.nrow();
+ int m = dx.ncol();
+ double tmp1 = 0;double tmp2 = 0;
+ NumericMatrix S(m,m);
+
+ for(int i=0;i<n;i++){
+ S = Smake(A(i,_),m);
+ tmp1 += log(detcpp(S));
+ tmp2 += trace(solvecpp(S),dx(i,_)-h*b(i,_));
+ }
+ return tmp1+tmp2/h;
+}
\ No newline at end of file
More information about the Yuima-commits
mailing list