最も近い値を持つ要素のインデックスを返す

掲題の件、あると思います。

例えばこういうやつ（z(=3)の値に最も近いxの要素(=2.9)を探す）。

> z <- 3
> x <- c(1,2,2.9,3.5,4)
> index <- which.min(abs(z - x))
> x[index]
[1] 2.9

これをzがベクトルのときでも動くようにしたいと考えた２実装が以下。

１つ目は単純に sapply を使っていて、２つ目は findInterval で区間を見つけつつちょいと補正（単純にやると最も近い値にならないときがある）をしたもの。

1億個要素を持つベクトル（x）と、複数個値を持ってる適当なベクトル(y)を使って２つのやり方の時間を測ると

> x <- sort(runif(10^8, min=1, max=4))
> ys <- seq(2,3.9,by=0.3)
> system.time({ ce_1 <- sapply(ys, function(z)which.min(abs(z - x))) })
   ユーザ   システム       経過  
     2.592      1.650      4.362 
> system.time({ ce_2 <- findInterval(ys, 0.5*(x[-length(x)] + x[-1])) + 1 })
   ユーザ   システム       経過  
     1.285      0.544      1.830

findInterval で書いたほうが早いんで良さそう。

念の為結果をチェックすると

> all(ce_1 == ce_2)
[1] TRUE

ちゃんとあってる。

まぁ、パッケージにするまでもないし、自分で関数として持っておくか〜という決断をする。

2021-01-12

大きめのSparse行列をnormalizeしたい時の書き方

掲題の件、あると思います。

私は今まで行列っぽいデータを normalize するとき*1は以下のように書いていました。

最後の結果を見たらわかるように行ごとに和をとると１になるようになっています。

> x <- matrix(1:4,2)
> x
     [,1] [,2]
[1,]    1    3
[2,]    2    4
> x/rowSums(x)
          [,1]      [,2]
[1,] 0.2500000 0.7500000
[2,] 0.3333333 0.6666667
> x <- matrix(1:4,2)
> x
     [,1] [,2]
[1,]    1    3
[2,]    2    4
> x/rowSums(x)
          [,1]      [,2]
[1,] 0.2500000 0.7500000
[2,] 0.3333333 0.6666667

適当に作ったSparse行列でもこれで動くと言えば動いている。

> library("Matrix")
> x <- as(matrix(c(1,2,0,4,0,0,7,8,9),3),"sparseMatrix")
> x
3 x 3 sparse Matrix of class "dgCMatrix"
          
[1,] 1 4 7
[2,] 2 . 8
[3,] . . 9
> x/Matrix::rowSums(x)
3 x 3 sparse Matrix of class "dgCMatrix"
                                   
[1,] 0.08333333 0.3333333 0.5833333
[2,] 0.20000000 .         0.8000000
[3,] .          .         1.0000000

でも、行列のサイズがでかくなるとアウト、なんか一回Denseな行列にしてるんだろうと思われる。

ここの例は恣意的データに見えるがもともと俺が使ってハマったデータがこんなんだったのでこうなっている。

> size <- 10^5
> i <- c(sample(1:size, 3*10^7, replace=TRUE), size)
> j <- c(sample(1:size, 3*10^7, replace=TRUE), size)
> x <- sparseMatrix(i,j,x=rnorm(3*10^7 + 1), dims=rep(size + 1, 2))
> x <- as(x + t(x), "symmetricMatrix")
> str(x)
Formal class 'dsCMatrix' [package "Matrix"] with 7 slots
  ..@ i       : int [1:29910604] 2 23 14 6 13 28 18 24 0 19 ...
  ..@ p       : int [1:100002] 0 0 0 0 0 0 0 0 0 0 ...
  ..@ Dim     : int [1:2] 100001 100001
  ..@ Dimnames:List of 2
  .. ..$ : NULL
  .. ..$ : NULL
  ..@ x       : num [1:29910604] -1.511 -0.878 0.777 1.258 1.185 ...
  ..@ uplo    : chr "U"
  ..@ factors : list()
> x/Matrix::rowSums(x)
 エラー:  ベクトルのメモリを使い切りました (上限に達した？)

こういうときは等価な計算のこれを使おう、これならいける

> Matrix::Diagonal(x=1 / Matrix::rowSums(x)) %*% x
100001 x 100001 sparse Matrix of class "dgCMatrix"
                                                                                                                              
 [1,] . . . . . . . . . . . . . . . . . . . . . .  .          . . . . . . . . . . . . . . . . .  .        . . . . . . . ......
 [2,] . . . . . . . . . . . . . . . . . . . . . .  .          . . . . . . . . . . . . . . . . .  .        . . . . . . . ......
 [3,] . . . . . . . . . . . . . . . . . . . . . . -0.08169018 . . . . . . . . . . . . . . . . .  .        . . . . . . . ......

ちゃんと計算も合う（確認）

> x <- as(matrix(c(1,2,0,4,0,0,7,8,9),3),"sparseMatrix")
> x/Matrix::rowSums(x)
3 x 3 sparse Matrix of class "dgCMatrix"
                                   
[1,] 0.08333333 0.3333333 0.5833333
[2,] 0.20000000 .         0.8000000
[3,] .          .         1.0000000
> Matrix::Diagonal(x=1 / Matrix::rowSums(x)) %*% x
3 x 3 sparse Matrix of class "dgCMatrix"
                                   
[1,] 0.08333333 0.3333333 0.5833333
[2,] 0.20000000 .         0.8000000
[3,] .          .         1.0000000

元ネタ

https://stat.ethz.ch/pipermail/r-help/2017-May/446704.html

*1:各行ごとに足したら１になる、でnormalizeを定義。他L2ノルムなどによるnormalizeも同様

2020-08-26

Pythonの純粋仮想関数は引数の数適当に変えて実装してもOK

python

掲題の件、そういいうことです。例えば適当に get() メソッドの引数を追加して実装しても

from abc import ABCMeta, abstractmethod

class Hoge(metaclass=ABCMeta):
    @abstractmethod
    def get(self):
        pass

class Moge(Hoge):
    def get(self, x = None, y = None):
        print("hi")

ちゃんと動く＆元のクラスのインスタンスは作れない

>>> x = Moge()
>>> x.get(y=123)
hi
>>> 
>>> Hoge()
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
TypeError: Can't instantiate abstract class Hoge with abstract methods get

2020-08-20

米国株式とUSDJPYの2020年1月1日からの相関構造

ドル円（USDJPY）が 0.6% でダウ平均が2.8% くらいのリスク（標準偏差、日率）

> 100*sd(df$`usdjpy=x`)
[1] 0.6186612
> 100*sd(df$dji)
[1] 2.842537

全体

library("dplyr")
library("tidyr")
library("stringr")
library("tidyquant")
library("PerformanceAnalytics")
get_from_yahoo <- function(x, range_date){
  # Why max(range_date) +1? because: from <= x < to
  df <- tq_get(x, from = min(range_date), to = max(range_date) + 1)
  dplyr::select(df, symbol, date, adjusted) %>%
    rename(price=adjusted) %>%
    mutate(symbol = str_to_lower(str_replace_all(symbol, "\\^", ""))) %>%
    arrange(date)
}
arithmetic_return <- function(x){
  c(NA, x[-1]/x[-length(x)] - 1)
}
# https://finance.yahoo.com/quote/%5EDJI/
range_date <- c(lubridate::ymd("20200101"), lubridate::ymd("20200818"))
symbols <- c("^DJI", "USDJPY=X")
df <- purrr::map_df(symbols, ~ get_from_yahoo(.x, range_date)) %>%
  tidyr::pivot_wider(names_from = symbol, values_from = price) %>%
  arrange(date) %>%
  mutate_at(
    vars(-"date"),
    arithmetic_return
  ) %>%
  na.omit()

chart.Correlation(dplyr::select(df, -date), histogram=TRUE, pch=19)
100*sd(df$`usdjpy=x`)
100*sd(df$dji)

2020-08-19

米国株式市場の2020年1月1日からの相関構造

すげぇ強かった（おしまい）

"^DJI" : ダウ平均
"^RUT": ラッセル指数
"^GSPC": S&P500
"^IXIC": ナスダック指数

再現用のCode

library("dplyr")
library("tidyr")
library("stringr")
library("tidyquant")
library("PerformanceAnalytics")
get_from_yahoo <- function(x, range_date){
  # Why max(range_date) +1? because: from <= x < to
  df <- tq_get(x, from = min(range_date), to = max(range_date) + 1)
  dplyr::select(df, symbol, date, adjusted) %>%
    rename(price=adjusted) %>%
    mutate(symbol = str_to_lower(str_replace_all(symbol, "\\^", ""))) %>%
    arrange(date)
}
arithmetic_return <- function(x){
  c(NA, x[-1]/x[-length(x)] - 1)
}
# https://finance.yahoo.com/quote/%5EDJI/
range_date <- c(lubridate::ymd("20200101"), lubridate::ymd("20200818"))
symbols <- c("^DJI", "^RUT", "^GSPC", "^IXIC")
df <- purrr::map_df(symbols, ~ get_from_yahoo(.x, range_date)) %>%
  tidyr::pivot_wider(names_from = symbol, values_from = price) %>%
  arrange(date) %>%
  mutate_at(
    vars(-"date"),
    arithmetic_return
  ) %>%
  na.omit()

chart.Correlation(dplyr::select(df, -date), histogram=TRUE, pch=19)

2020-08-12

特定の列だけ外して処理したい場合には . (dot) を名前に持つ列と lsを組み合わせると良い

株式会社ホクソエムの社長から教えていただいた。いちいち dplyr::select なぞせんでもこうするだけで .y を除いた処理を実行できる。

> df <- data.frame(
+   x = 1:3,
+   .y = 2:4
+ )
> df
  x .y
1 1  2
2 2  3
3 3  4
> ls(df)
[1] "x"

2020-04-17

dplyrでうっかり八兵衛にならないように !! を意識する

うっかりやっちまったのでメモっておく。まず適当な data.frame() を定義する。

> df <- data.frame(x=sample(1:2, 10, replace=TRUE), y=sample(c("a", "b"), 10, replace = TRUE))
> df
   x y
1  2 a
2  2 b
3  2 b
4  2 b
5  1 b
6  1 a
7  2 b
8  1 b
9  1 b
10 2 b

で、これを適当にFilteringしたいとするとこう書く。

> dplyr::filter(df, x == 1, y == "b")
  x y
1 1 b
2 1 b
3 1 b

そして「よくよく考えると x を変数にしておいたほうが楽そうだぞ〜」ということで以下のように書くと うっかり八兵衛になる （答えが狂う）。

> x <- 1
> dplyr::filter(df, x == x, y == "b")
  x y
1 2 b
2 2 b
3 2 b
4 1 b
5 2 b
6 1 b
7 1 b
8 2 b

正しくは以下のように rlang::!! をちゃんとかませておかないといけない。

> dplyr::filter(df, x == !!x, y == "b")
  x y
1 1 b
2 1 b
3 1 b

rlang は偉大なり、うっかりうっかり。