[Yuima-commits] r340 - in pkg/yuima: . R src
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Fri Oct 17 23:15:54 CEST 2014
Author: lorenzo
Date: 2014-10-17 23:15:53 +0200 (Fri, 17 Oct 2014)
New Revision: 340
Modified:
pkg/yuima/DESCRIPTION
pkg/yuima/R/qmle.R
pkg/yuima/src/carmafilter.c
Log:
updated carmafilter.c
Modified: pkg/yuima/DESCRIPTION
===================================================================
--- pkg/yuima/DESCRIPTION 2014-10-14 11:30:29 UTC (rev 339)
+++ pkg/yuima/DESCRIPTION 2014-10-17 21:15:53 UTC (rev 340)
@@ -1,8 +1,8 @@
Package: yuima
Type: Package
Title: The YUIMA Project package for SDEs
-Version: 1.0.39
-Date: 2014-10-14
+Version: 1.0.40
+Date: 2014-10-17
Depends: methods, zoo, stats4, utils, expm, cubature, mvtnorm
Author: YUIMA Project Team
Maintainer: Stefano M. Iacus <stefano.iacus at unimi.it>
Modified: pkg/yuima/R/qmle.R
===================================================================
--- pkg/yuima/R/qmle.R 2014-10-14 11:30:29 UTC (rev 339)
+++ pkg/yuima/R/qmle.R 2014-10-17 21:15:53 UTC (rev 340)
@@ -1831,7 +1831,7 @@
# u<-diff(tt)[1]
-# Amatx<-yuima.carma.eigen(A)
+# Amatx<-yuima.ca2rma.eigen(A)
# expA<-Amatx$vectors%*%expm(diag(Amatx$values*u),
# method="Pade",
# order=6,
@@ -1887,7 +1887,7 @@
# input: V_inf, p, expA
# output: SigMatr (pre-alloc in R)
#
- SigMatr <- .Call("carma_tmp", V_inf, as.integer(p), matrix(expA,p,p), PACKAGE="yuima")
+ SigMatr <- .Call("carma_tmp", V_inf, as.integer(p), expA, PACKAGE="yuima")
#
# SigMatr <- V_inf - expA %*% V_inf %*% expAT
@@ -1897,45 +1897,184 @@
# stop("")
#
statevar<-matrix(rep(0, p),p,1)
- Qmatr<-SigMatr
+ Qmatr<-SigMatr+0
# set
#statevar<-statevar0
- # SigMatr<-expA%*%V_inf%*%t(expA)+Qmatr
+ # SigMatCheckr<-expA%*%V_inf%*%t(expA)+Qmatr
#SigMatr<-Qmatr
+
#SigMatr <- V_inf
zc <- matrix(bvector,1,p)
loglstar <- 0
- loglstar1 <- 0
+
# zcT<-matrix(bvector,p,1)
zcT <- t(zc)
- ### statevar, expA, times.obs, Qmatr, SigMatr, zc
- for(t in 1:times.obs){
- # prediction
- statevar <- expA %*% statevar
- SigMatr <- expA %*% SigMatr %*% expAT + Qmatr
- # forecast
- Uobs <- y[t] - zc %*% statevar
- dum.zc <- zc %*% SigMatr
- sd_2 <- dum.zc %*% zcT
- Inv_sd_2 <- 1/sd_2
- #correction
- Kgain <- SigMatr %*% zcT %*% Inv_sd_2
- statevar <- statevar+Kgain %*% Uobs
- SigMatr <- SigMatr - Kgain %*% dum.zc
- term_int<- -0.5 * (log(sd_2)+ Uobs %*% Uobs %*% Inv_sd_2)
- loglstar <- loglstar + term_int
- }
- return(list(loglstar=(loglstar-0.5*log(2*pi)*times.obs),s2hat=sd_2))
+# Cycle_Carma(SEXP StateVar, SEXP ExpA, SEXP Times.Obs, SEXP P,
+# SEXP Qmatr, SEXP SigMatr, SEXP Zc, SEXP Logstar)
+ # sd_2<-0
+ loglstar<- .Call("Cycle_Carma", y, statevar, expA, as.integer(times.obs),
+ as.integer(p), Qmatr, SigMatr, zc,
+ PACKAGE="yuima")
+ return(list(loglstar=as.numeric(loglstar[1,1])-0.5*log(2*pi)*times.obs,s2hat=loglstar[2,1]))
+
+ # statevar, expA, times.obs, Qmatr, SigMatr, zc, logstar
+
+ #loglstar<-CycleCarma(y, statevar, expA, as.integer(times.obs), as.integer(p), Qmatr, SigMatr, zc, loglstar)
+#return(list(loglstar=loglstar[1,1]-0.5*log(2*pi)*times.obs,s2hat=loglstar[2,1]))
+# for(t in 1:times.obs){
+# # prediction
+# statevar <- expA %*% statevar
+# SigMatr <- expA %*% SigMatr %*% expAT + Qmatr
+# # forecast
+# Uobs <- y[t] - zc %*% statevar
+# dum.zc <- zc %*% SigMatr
+# sd_2 <- dum.zc %*% zcT
+# Inv_sd_2 <- 1/sd_2
+# #correction
+# Kgain <- SigMatr %*% zcT %*% Inv_sd_2
+# statevar <- statevar+Kgain %*% Uobs
+# SigMatr <- SigMatr - Kgain %*% dum.zc
+# term_int<- -0.5 * (log(sd_2)+ Uobs %*% Uobs %*% Inv_sd_2)
+# loglstar <- loglstar + term_int
+# }
+# return(list(loglstar=loglstar-0.5*log(2*pi)*times.obs,s2hat=sd_2))
}
+# CycleCarma<-function(y, statevar, expA, times.obs=integer(),
+# p=integer(), Qmatr, SigMatr, zc, loglstar){
+# # expAT=t(expA)
+# # zcT=t(zc)
+# # for(t in 1:times.obs){
+# # t=1
+# # # # prediction
+# # statevar <- expA %*% statevar
+# # SigMatr <- expA %*% SigMatr %*% t(expA) + Qmatr
+# # # forecast
+# # Uobs <- y[t] - zc %*% statevar # 1 - 1Xp px1
+# # dum.zc <- zc %*% SigMatr # 1xp pxp
+# # sd_2 <- dum.zc %*% t(zc) # 1xp px1
+# # Inv_sd_2 <- 1/sd_2
+# # #correction
+# # Kgain <- SigMatr %*% t(zc) %*% Inv_sd_2 # pxp px1*1
+# # statevar <- statevar+Kgain %*% Uobs # px1+px1
+# # SigMatr <- SigMatr - Kgain %*% dum.zc # pxp-px1 1x+
+# # term_int<- -0.5 * (log(sd_2)+ Uobs %*% Uobs %*% Inv_sd_2) # every entries are scalars
+# # loglstar <- loglstar + term_int # every entries are scalars
+# # }
+# # expA=matrix(c(1:16),nrow=4,ncol=4)
+# # SigMatr=matrix(c(1:16),nrow=4,ncol=4)+1
+# # Qmatr=matrix(c(1:16),nrow=4,ncol=4)+2
+# # vvvvv<-expA%*%SigMatr
+# # ppppp<-expA%*%SigMatr%*%t(expA)+Qmatr
+# rY=as.numeric(y)
+# rStateVar=as.numeric(statevar)
+# rExpA=as.numeric(expA)
+# rtime_obs=times.obs
+# p=p
+# rQmatr=as.numeric(Qmatr)
+# rSigMatr=as.numeric(SigMatr)
+# rZc=as.numeric(zc)
+# rLogstar=loglstar
+# In_dum=0
+# sd_2=0
+# rMat21=numeric(length=p)
+# rdum_zc=numeric(length=p)
+# rMat22int=numeric(length=p*p)
+# rMat22est=numeric(length=p*p)
+# rKgain=numeric(length=p)
+# for(t in 1:rtime_obs){
+# # prediction
+# for(i in 1:p){
+# rMat21[(i-1)+1] = 0
+# for(j in 1:p){
+# # statevar <- expA %*% statevar: px1=pxp px1
+# rMat21[(i-1)+1] = rMat21[(i-1)+1]+rExpA[(i-1)+(j-1)*p+1]*rStateVar[(j-1)+1]
+# }
+# rStateVar[(i-1)+1] = rMat21[(i-1)+1] # statevar <- expA %*% statevar
+# }
+#
+# # SigMatr <- expA %*% SigMatr %*% expAT + Qmatr: pxp = pxp pxp pxp
+# # First We compute rMat22int <- expA %*% SigMatr : pxp = pxp pxp
+# for(i in 1:p){
+# for(j in 1:p){
+# rMat22int[(i-1)+(j-1)*p+1]=0
+# for(h in 1:p){
+# rMat22int[(i-1)+(j-1)*p+1]=rMat22int[(i-1)+(j-1)*p+1]+rExpA[(i-1)+(h-1)*p+1]*
+# rSigMatr[(h-1)+(j-1)*p+1]
+# }
+# }
+# }
+# # Second We compute rMat22est <- rMat22int %*% t(expA) + Qmatr: pxp = pxp pxp + pxp
+# for(i in 1:p){
+# for(j in 1:p){
+# rMat22est[(i-1)+(j-1)*p+1]=0
+# for(h in 1:p){
+# rMat22est[(i-1)+(j-1)*p+1]=rMat22est[(i-1)+(j-1)*p+1]+rMat22int[(i-1)+(h-1)*p+1]*rExpA[(j-1)+(h-1)*p+1]
+#
+# }
+# rSigMatr[(i-1)+(j-1)*p+1]=rMat22est[(i-1)+(j-1)*p+1]+rQmatr[(i-1)+(j-1)*p+1]
+# }
+# }
+# # # forecast
+#
+# # Uobs <- y[t] - zc %*% statevar # 1 - 1Xp px1
+# rMat22est[1]=0
+# for(i in c(1:p)){
+# rMat22est[1]=rMat22est[1]+rZc[i]*rStateVar[i]
+# }
+# Uobs=rY[t]-rMat22est[1]
+#
+# # dum.zc <- zc %*% SigMatr # 1xp pxp
+#
+#
+# for(i in c(1:p)){
+# rdum_zc[i]=0
+# for(j in c(1:p)){
+# rdum_zc[i]=rdum_zc[i]+rZc[j]*rSigMatr[(i-1)*h+j-1+1]
+# }
+# }
+# # sd_2 <- dum.zc %*% zcT # 1xp px1
+# sd_2=0
+# for(i in c(1:p)){
+# sd_2=sd_2+rdum_zc[i]*rZc[i]
+# }
+# # #correction
+# # Kgain <- SigMatr %*% zcT %*% 1/sd_2 # pxp px1*1
+# for(i in c(1:p)){
+# rMat21[i]=0
+# for(j in c(1:p)){
+# rMat21[i]=rMat21[i]+rSigMatr[(i-1)+(j-1)*p+1]*rZc[j]
+# }
+# rKgain[i]=rMat21[i]/sd_2
+# }
+#
+#
+# # statevar <- statevar+Kgain %*% Uobs # px1+px1
+# for(i in c(1:p)){
+# rStateVar[i] = rStateVar[i] + rKgain[i]*Uobs
+# }
+# # SigMatr <- SigMatr - Kgain %*% dum.zc # pxp-px1 1xp
+# for(i in c(1:p)){
+# for(j in c(1:p)){
+# rSigMatr[(i-1)+(j-1)*p+1] =rSigMatr[(i-1)+(j-1)*p+1]-rKgain[i]*rdum_zc[j]
+# }
+# }
+#
+# term_int = -0.5 * (log(sd_2)+ Uobs * Uobs * 1/sd_2) # every entries are scalars
+# loglstar = loglstar + term_int # every entries are scalars
+#
+#
+# }
+# Res<-matrix(c(loglstar,sd_2),nrow=2,ncol=1)
+# return(Res)
+# }
-
yuima.PhamBreton.Alg<-function(a){
p<-length(a)
gamma<-a[p:1]
Modified: pkg/yuima/src/carmafilter.c
===================================================================
--- pkg/yuima/src/carmafilter.c 2014-10-14 11:30:29 UTC (rev 339)
+++ pkg/yuima/src/carmafilter.c 2014-10-17 21:15:53 UTC (rev 340)
@@ -43,6 +43,11 @@
SEXP carma_tmp(SEXP V, SEXP P, SEXP A);
+SEXP Cycle_Carma(SEXP Y, SEXP StateVar, SEXP ExpA, SEXP Times_Obs, SEXP P,
+ SEXP Qmatr, SEXP SigMatr, SEXP Zc);
+
+
+
SEXP carma_tmp(SEXP V, SEXP P, SEXP A){
int p;
@@ -73,8 +78,6 @@
PROTECT(Sigma = allocMatrix(REALSXP, p, p));
rSigma = REAL(Sigma);
-
-
/* B = A %*% V */
for(i=0; i<p; i++){
for(j=0; j<p; j++){
@@ -83,27 +86,365 @@
rB[i+j*p] = rB[i+j*p] + rA[i+h*p] * rV[h+j*p];
}
}
- }
-
-
- /* C = B %*% A^T */
- for(i=0; i<p; i++){
+ /*}
+ C = B %*% A^T
+ for(i=0; i<p; i++){*/
for(j=0; j<p; j++){
rC[i+j*p] = 0;
for(h=0; h<p; h++){
rC[i+j*p] = rC[i+j*p] + rB[i+h*p] * rA[j+h*p];
}
+ rSigma[i+j*p] = rV[i+j*p] - rC[i+j*p];
}
}
-
- for(i=0; i<p*p; i++){
- rSigma[i] = rV[i] - rC[i];
- }
-
+ /*for(i=0; i<p*p; i++){
+
+ }*/
UNPROTECT(5);
- return(Sigma);
+ return Sigma;
}
+SEXP Cycle_Carma(SEXP Y, SEXP StateVar, SEXP ExpA, SEXP Times_Obs, SEXP P,
+ SEXP Qmatr, SEXP SigMatr, SEXP Zc){
+
+ /* Declare Integer Variable */
+
+ int times_obs, p;
+ int i, j, h, k, t;
+
+ /* Declare pointer */
+
+ double *rY, *rStateVar, *rExpA, *rQmatr, *rSigMatr, *rZc;
+ double *rKgain, *rdum_zc, *rMat21, *rMat22int, *rMat22est, *rResult;
+ double ppp=0;
+ double Uobs=0;
+ double sd_2=0;
+ double LogL=0;
+
+ /* Declare SEXP */
+
+ SEXP Kgain, dum_zc, Mat21, Mat22int, Mat22est, Result;
+
+ /* Check the type of variables*/
+
+ if(!isInteger(P)) error("`P' must be integer");
+ if(!isInteger(Times_Obs)) error("`Times_Obs' must be integer");
+ if(!isNumeric(Y)) error("`Y' must be numeric");
+ if(!isNumeric(StateVar)) error("`StateVar' must be numeric");
+ if(!isNumeric(ExpA)) error("`ExpA' must be numeric");
+ if(!isNumeric(Qmatr)) error("`Qmatr' must be numeric");
+ if(!isNumeric(SigMatr)) error("`SigMatr' must be numeric");
+ if(!isNumeric(Zc)) error("`Zc' must be numeric");
+
+
+ /* Protect Objects */
+
+ PROTECT(Y = AS_NUMERIC(Y));
+ rY = REAL(Y);
+
+ PROTECT(StateVar = AS_NUMERIC(StateVar));
+ rStateVar = REAL(StateVar);
+
+ PROTECT(ExpA = AS_NUMERIC(ExpA));
+ rExpA = REAL(ExpA);
+
+ PROTECT(Qmatr = AS_NUMERIC(Qmatr));
+ rQmatr = REAL(Qmatr);
+
+ PROTECT(SigMatr = AS_NUMERIC(SigMatr));
+ rSigMatr = REAL(SigMatr);
+
+ PROTECT(Zc = AS_NUMERIC(Zc));
+ rZc = REAL(Zc);
+
+ /*PROTECT(Logstar = AS_NUMERIC(Logstar));
+ rLoglstar = REAL(Logstar);*/
+
+
+ /* Declare dimensions for the state variables and observations */
+ times_obs = *INTEGER(Times_Obs);
+ p = *INTEGER(P);
+
+ PROTECT(Mat21 = allocMatrix(REALSXP, p, 1));
+ rMat21 = REAL(Mat21);
+
+ PROTECT(Kgain = allocMatrix(REALSXP, p, 1));
+ rKgain = REAL(Kgain);
+
+ PROTECT(dum_zc = allocMatrix(REALSXP, p,1));
+ rdum_zc = REAL(dum_zc);
+
+ PROTECT(Mat22int = allocMatrix(REALSXP, p, p));
+ rMat22int = REAL(Mat22int);
+
+ PROTECT(Mat22est = allocMatrix(REALSXP, p, p));
+ rMat22est = REAL(Mat22est);
+
+ PROTECT(Result = allocMatrix(REALSXP, 2, 1));
+ rResult = REAL(Result);
+
+
+ /* Main Code */
+ /* Dimension of Inputs:
+ Y = Vector p dimension;
+ StateVar = matrix p x 1;
+ ExpA = matrix p x p;
+ Times_Obs = integer;
+ P = integer;
+ Qmatr = matrix p x p;
+ SigMatr = matrix p x p;
+ Zc = matrix 1 x p;
+ Logstar = scalar.
+
+ for(t in 1:rtime_obs){*/
+
+ /* rResult[0]=0;*/
+ /* for(i=0; i<p; i++){
+ rdum_zc[i]=0;
+ }*/
+
+ for(t=0; t<times_obs; t++){
+ /*t=0;*/
+ /* prediction */
+ for(i=0; i<p; i++){
+ rMat21[i] = 0;
+ for(j=0; j<p; j++){
+ /* statevar <- expA %*% statevar: px1=pxp px1 */
+ rMat21[i] = rMat21[i]+rExpA[i+j*p]*rStateVar[j];
+ }
+ rStateVar[i] = rMat21[i];
+ /* rMat21[i]=0; statevar <- expA %*% statevar */
+ /*}
+ SigMatr <- expA %*% SigMatr %*% expAT + Qmatr: pxp = pxp pxp pxp
+ First We compute rMat22int <- expA %*% SigMatr : pxp = pxp pxp
+ for(i=0; i<p; i++){*/
+ for(j=0; j<p; j++){
+ rMat22int[i+j*p]=0;
+ for(h=0; h<p; h++){
+ rMat22int[i+j*p]=rMat22int[i+j*p]+rExpA[i+h*p]*
+ rSigMatr[h+j*p];
+ }
+ }
+ }
+ /* Second We compute rMat22est <- rMat22int %*% t(expA)
+ + Qmatr: pxp = pxp pxp + pxp */
+ for(i=0; i<p; i++){
+ for(j=0; j<p; j++){
+ rMat22est[i+j*p]=0;
+ for(h=0; h<p; h++){
+ rMat22est[i+j*p]=rMat22est[i+j*p]+rMat22int[i+h*p]*rExpA[j+h*p];
+
+ }
+ rSigMatr[i+j*p]=rMat22est[i+j*p]+rQmatr[i+j*p];
+ }
+
+ }
+
+ /*forecast*/
+
+ /*Uobs <- y[t] - zc %*% statevar # 1 - 1Xp px1*/
+ rMat22est[0]=0;
+ for(i=0; i<p; i++){
+ rMat22est[0]=rMat22est[0]+rZc[i]*rStateVar[i];
+ }
+ Uobs=rY[t]-rMat22est[0];
+ /*rMat22est[0]=0;*/
+ /* dum.zc <- zc %*% SigMatr # 1xp pxp */
+ for(i=0; i<p; i++){
+ rMat21[i] = 0;
+ for(j=0; j<p; j++){
+ rMat21[i] =rMat21[i]+rSigMatr[i+j*p]*rZc[j];
+ }
+ rdum_zc[i]=rMat21[i];
+ }
+
+
+ /*rResult[1]=rSigMatr[0];*/
+ /* Sd_2 <- dum.zc %*% zcT # 1xp px1 */
+ rResult[1] =0;
+ for(i=0; i<p; i++){
+ rResult[1] = rResult[1]+rdum_zc[i]*rZc[i];
+ }
+
+ /* correction */
+ /* Kgain <- SigMatr %*% zcT %*% 1/Sd_2 # pxp px1*1 */
+ for(i=0; i<p; i++){
+ rMat21[i]=0;
+ for(j=0; j<p; j++){
+ rMat21[i]=rMat21[i]+rSigMatr[i+j*p]*rZc[j];
+ }
+ rKgain[i]=rMat21[i]/rResult[1];
+ /*rMat21[i]=0;*/
+ }
+ /* statevar <- statevar+Kgain %*% Uobs # px1+px1 */
+ for(i=0; i<p; i++){
+ rStateVar[i] = rStateVar[i] + rKgain[i]*Uobs;
+ }
+
+
+ /* SigMatr <- SigMatr - Kgain %*% dum.zc # pxp-px1 1xp */
+ for(i=0; i<p; i++){
+ for(j=0; j<p; j++){
+ rSigMatr[i+j*p]=rSigMatr[i+j*p]-rKgain[i]*rdum_zc[j];
+ }
+ }
+ /*for(i=0; i<p*p; i++){
+ rMat22int[i];
+ }*/
+ /*term_int = -0.5 * (log(Sd_2)+ Uobs * Uobs * 1/Sd_2) every entries are scalars*/
+ rResult[0] = rResult[0] -0.5 * (log(rResult[1])+ Uobs * Uobs /rResult[1]);
+ /* manual debug */
+
+
+ }
+
+ UNPROTECT(12);
+ return Result;
+
+ }
+
+ /* t=1;*/
+
+
+ /* next step of the loop */
+
+
+ /* for(i=0; i<p; i++){
+ rMat21[i] = 0;
+ for(j=0; j<p; j++){
+ /* statevar <- expA %*% statevar: px1=pxp px1 */
+ /* rMat21[i] = rMat21[i]+rExpA[i+j*p]*rStateVar[j];
+ }
+ rStateVar[i] = rMat21[i];
+ rMat21[i]=0;/* statevar <- expA %*% statevar */
+ /* }
+
+
+
+ /* SigMatr <- expA %*% SigMatr %*% expAT + Qmatr: pxp = pxp pxp pxp */
+ /* First We compute rMat22int <- expA %*% SigMatr : pxp = pxp pxp */
+ /* for(i=0; i<p; i++){
+ for(j=0; j<p; j++){
+ rMat22int[i+j*p]=0;
+ for(h=0; h<p; h++){
+ rMat22int[i+j*p]=rMat22int[i+j*p]+rExpA[i+h*p]*
+ rSigMatr[h+j*p];
+ }
+ }
+ }
+
+ /* print(rMat22int) */
+ /* Second We compute rMat22est <- rMat22int %*% t(expA) + Qmatr: pxp = pxp pxp + pxp */
+ /* for(i=0; i<p; i++){
+ for(j=0; j<p; j++){
+ rMat22est[i+j*p]=0;
+ for(h=0; h<p; h++){
+ rMat22est[i+j*p]=rMat22est[i+j*p]+rMat22int[i+h*p]*rExpA[j+h*p];
+
+ }
+ rSigMatr[i+j*p]=rMat22est[i+j*p];
+
+ }
+
+ }
+
+
+
+
+ for(i=0; i<p*p; i++){
+ rSigMatr[i]=rSigMatr[i]+rQmatr[i];
+ }
+
+
+ for(i=0; i<p*p; i++){
+ rMat22est[i]=0;
+ rMat22int[i]=0;
+ }
+
+ /*forecast*/
+
+ /*Uobs <- y[t] - zc %*% statevar # 1 - 1Xp px1*/
+ /* rMat22est[0]=0;
+ for(i=0; i<p; i++){
+ rMat22est[0]=rMat22est[0]+rZc[i]*rStateVar[i];
+ }
+ Uobs=rY[t]-rMat22est[0];
+ rMat22est[0]=0;
+ /* dum.zc <- zc %*% SigMatr # 1xp pxp */
+
+
+
+ /*
+
+
+
+ for(i=0; i<p; i++){
+ rResult[1] = rResult[1]+rdum_zc[i]*rZc[i];
+ }*/
+
+
+
+ /* for(i=0; i<p; i++){
+ rMat21[i] = 0;
+ for(j=0; j<p; j++){
+ rMat21[i] =rMat21[i]+rSigMatr[i+j*p]*rZc[j];
+ }
+ rdum_zc[i]=rMat21[i];
+ rMat21[i]=0;
+ }
+
+ /*rResult[1]=rSigMatr[0];*/
+ /* Sd_2 <- dum.zc %*% zcT # 1xp px1 */
+ /* rResult[1] =0;
+ for(i=0; i<p; i++){
+ rResult[1] = rResult[1]+rdum_zc[i]*rZc[i];
+ }
+
+ /* correction */
+ /* Kgain <- SigMatr %*% zcT %*% 1/Sd_2 # pxp px1*1 */
+ /* for(i=0; i<p; i++){
+ rMat21[i]=0;
+ for(j=0; j<p; j++){
+ rMat21[i]=rMat21[i]+rSigMatr[i+j*p]*rZc[j];
+ }
+ rKgain[i]=rMat21[i]/rResult[1];
+ rMat21[i]=0;
+ }
+
+
+ /* statevar <- statevar+Kgain %*% Uobs # px1+px1 */
+ /* for(i=0; i<p; i++){
+ rStateVar[i] = rStateVar[i] + rKgain[i]*Uobs;
+ }
+
+
+ /* SigMatr <- SigMatr - Kgain %*% dum.zc # pxp-px1 1xp */
+ /* for(i=0; i<p; i++){
+ for(j=0; j<p; j++){
+ rSigMatr[i+j*p]=rSigMatr[i+j*p]-rKgain[i]*rdum_zc[j];
+ }
+ }
+ /*for(i=0; i<p*p; i++){
+ rMat22int[i];
+ }*/
+ /*term_int = -0.5 * (log(Sd_2)+ Uobs * Uobs * 1/Sd_2) every entries are scalars*/
+ /* rResult[0] = rResult[0] -0.5 * (log(rResult[1])+ Uobs * Uobs /rResult[1]);
+
+
+
+
+ /* Unprotect steps */
+ /* UNPROTECT(12);
+ return(Result);
+
+
+
+
+
+
+}
+
+
/*
SEXP sde_sim_euler(SEXP x0, SEXP t0, SEXP delta, SEXP N, SEXP M,
SEXP d,888 SEXP s, SEXP sx,
More information about the Yuima-commits
mailing list