Multivariate S-map test 1

目的

• エラーの大きさがmultivariate S-mapの係数推定に与える影響を知りたい
• Multivariate S-mapにおける埋め込みメソッドがモデルの予測精度と係数推定に与える影響を知りたい

パッケージのロード

library(rEDM); packageVersion("rEDM") #"0.7.1"

## Warning: package 'rEDM' was built under R version 3.4.4

## [1] '0.7.1'

library(dplyr); packageVersion("dplyr") #"0.8.0.1"

## Warning: package 'dplyr' was built under R version 3.4.4

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

## [1] '0.8.0.1'

2種競争系モデル

• 以下の2種競争系モデルを定義する
• Nxt+1 = rx * Nxt * (1 - (Nxt + axy * Nyt) / Kx)
• Nyt+1 = ry * Nyt * (1 - (Nyt + ayx * Nxt) / Ky)
• Nxt+1 = … + (-rxaxy/Kx)Nxt*Nyt
• なので、∂Nxt/∂Nyt = (-rxaxy/Kx)Nxt となる

#Nxt+1 = rx * Nxt * (1 - (Nxt + axy * Nyt) / Kx)
#Nyt+1 = ry * Nyt * (1 - (Nyt + ayx * Nxt) / Ky)
#Nxt+1 = ... + (-rx*axy/Kx)*Nxt*Nyt

rx=1.4; axy=0.3; Kx=2.4
x.f <- function(v){
  rx*(1-(v[2,1]+axy*v[2,2])/Kx)*v[2,1]
}

ry=1.5; ayx=0.2; Ky=2
y.f <- function(v){
  ry*(1-(v[2,2]+ayx*v[2,1])/Ky)*v[2,2]
}

関数定義

• テスト時系列を作る関数test.ts
• Multivariate S-mapをNxtとNytだけで行う関数mvar_smap2
• Simplex projectionで埋め込み次元を決定する関数simp_emb
• Multivariate S-mapをNxtの時間遅れ埋め込みとNytで行う関数mvar_smapE

##テスト時系列を作る関数
test.ts <- function(x) {
  g <- matrix(c(1.1,1.02,0.98,0.99),nrow=2,ncol=2)
  
  ex <- rnorm(301, 0, x)
  ey <- rnorm(301, 0, x)
  
  for (i in 1:299) {
    gn <- c(x.f(g[i:(i+1),])+ex[i+1], y.f(g[i:(i+1),])+ey[i+1])
    g <- rbind(g,gn)
  } # g:1列目に種x, 2列目に種yの個体数, 301 時間ステップ分
  
  g <- cbind(g,rnorm(301,1.1,0.2)) #3列目に乱数の列であるzの動態を挿入
  dimnames(g) <- NULL
  
  g[51:100,]
}

##Multivariate S-mapをNxtとNytだけで行う関数
mvar_smap2 <- function(x) {
  mlt_smap_res0 <- block_lnlp(cbind(x[,1], x[,2]), tp = 1, target_column = 1,
                              method = "s-map", theta = seq(1,10,by=0.1), silent = T)
  Tm <- mlt_smap_res0[which.min(mlt_smap_res0$rmse), "theta"]
  block_lnlp(cbind(x[,1], x[,2]), tp = 1, target_column = 1,
                              method = "s-map", theta = Tm, save_smap_coefficients = T, silent = T)
}

##simplex projectionで埋め込み次元を決定する関数
simp_emb <- function(x){
  simp_x_res0 <- 
    simplex(x[,1], E=1:10, lib=c(1, 50), pred=c(1, 50), silent=TRUE)
  simp_x_res0[which.min(simp_x_res0$rmse), "E"]
}

##multivariate S-mapをフル埋め込み(末尾入れ替え)で行う関数
mvar_smapE <- function(x, E){
  mvarts <- cbind(x[,1], x[,2])
  for (i in 1:(E-2)){
    mvarts <- cbind(mvarts, lag(x[,1], n=as.integer((1:(E-2))[i])))
  }
  mlt_smap_res0 <- block_lnlp(mvarts, tp = 1, target_column = 1,
                              method = "s-map", theta = seq(1,10,by=0.1), silent = T)
  Tm <- mlt_smap_res0[which.min(mlt_smap_res0$rmse), "theta"]
  block_lnlp(mvarts, tp = 1, target_column = 1,
                              method = "s-map", theta = Tm, save_smap_coefficients = T, silent = T)
}

エラーがmultivariate S-mapの係数推定に与える影響

• 期待値的には、かなりうまくαを推定できていた（マイナス値を叩き出すケースもあるにはあった）
• エラーの大きさはαの推定精度(precisionもaccuracyも)にほとんど影響しなかった
• これは一見奇異な結果だが、エラー大だとNxtとNytの状態空間の”大きさ”が大きくなるだけで形は変わらないと考えると納得のいく結果

set.seed(1)
axy_simu <- NULL
axy_simu_df <- NULL
for (j in c(0.02, 0.06, 0.08, 0.12)){
  axy_simu <- NULL
  for (i in 1:100){
  subset.g <- test.ts(j)
  if(is.na(subset.g[50,1])|is.infinite(subset.g[50,1])){
    axy_simu <- append(axy_simu, NA)
  } else {
    mlt_smap_res1 <- mvar_smap2(subset.g)
    axy_simu <- append(axy_simu, (-mean(mlt_smap_res1$smap_coefficients[[1]]$c_2/subset.g[,1]/rx*Kx, na.rm=TRUE)))
  }
  }
  axy_simu_df <- rbind(axy_simu_df, data.frame(axy_simu=axy_simu, SD=j))
}

par(mfrow=c(1, 1), las=1, family="sans", lwd=2)
boxplot(axy_simu ~ SD, axy_simu_df, 
        ylim=c(-3, 5), 
        names=c("0.02", "0.06", "0.08", "0.12"), 
        xlab="Error variability (SD)", 
        ylab=paste("Estimated",expression(α), "(omitting outliers)"))
abline(h=0.3, lty=2)

dev.off()

## null device 
##           1

Multivariate S-mapにおける埋め込みメソッドが予測精度とα推定に与える影響（エラーのSD=0.1）

• Fully multivariateの方が予測精度がこころもち良く、α推定はほぼ差が無かった
• 力学系の変数をすべて使えているというのが結果をもたらしているかもしれない、例えば3種系だが2種しか観察できないような状況を仮定して検証を進めるべきかもしれない

set.seed(1)
check.par <- NULL
for (i in 1:100){
  subset.g <- test.ts(0.1)
  if(is.na(subset.g[50,1])|is.infinite(subset.g[50,1])){
  } else {
      mlt_smap_res1 <- mvar_smap2(subset.g)
      rho_obs <- mlt_smap_res1$rho
      rmse_obs <- mlt_smap_res1$rmse
      axy_obs <- (-mean(mlt_smap_res1$smap_coefficients[[1]]$c_2/subset.g[,1]/rx*Kx, na.rm=TRUE))
      check.par <- rbind(check.par, data.frame(rho_obs, rmse_obs, axy_obs))
      }
}
check.par_2 <- check.par 

set.seed(1)
check.par <- NULL
for (i in 1:100){
  subset.g <- test.ts(0.1)
  if(is.na(subset.g[50,1])){
  } else {
    if(is.infinite(subset.g[50,1])){
    } else {
      Ex <- simp_emb(subset.g)
      if(Ex<3){
      } else {
        mlt_smap_res1 <- mvar_smapE(subset.g, Ex)
        rho_obs <- mlt_smap_res1$rho
        rmse_obs <- mlt_smap_res1$rmse
        axy_obs <- (-mean(mlt_smap_res1$smap_coefficients[[1]]$c_2/subset.g[,1]/rx*Kx, na.rm=TRUE))
        check.par <- rbind(check.par, data.frame(rho_obs, rmse_obs, axy_obs))
      }
    }
  }
}
check.par_E <- check.par 

check.par <- rbind(check.par_2, check.par_E)
check.par$method <- as.factor(c(rep("2", times=nrow(check.par_2)), rep("E", times=nrow(check.par_E))))

par(mfrow=c(1, 3), las=1, family="sans", lwd=2, ps=15)
plot(rho_obs~method, subset(check.par), ylab=expression(ρ), xlab="",
     names=c("", ""))
axis(1, at=1:2, labels=c("Fully\nmultivariate", "Embedding\nw. lags"), padj=0.5)
plot(rmse_obs~method, subset(check.par, axy_obs<10), ylab="RMSE", xlab="",
     names=c("", ""))
axis(1, at=1:2, labels=c("Fully\nmultivariate", "Embedding\nw. lags"), padj=0.5)
plot(axy_obs~method, subset(check.par, axy_obs<10), ylab=paste("Estimated",expression(α), "(omitting outliers)"), xlab="",
     names=c("", ""))
axis(1, at=1:2, labels=c("Fully\nmultivariate", "Embedding\nw. lags"), padj=0.5)
abline(h=0.3, lty=2, lwd=0.5)

dev.off()

## null device 
##           1