https://coinmarketcap.com/api

アクセス回数などの制限はありますが、無料で使用できます。
ここではRStudioを使ってデータの取得からグラフ描画までやってみました。

きっかけは、最近ブロックチェーン技術の学習を進めていく上で、各プロジェクトで使われるコインの価値について比較したくなったことからです。仮想通貨全体に影響を受けてしまうような市場の動向をできるだけ排除して比較するためには、チャートを目で見ていてはわからないので、このようなAPIを使うことにしました。

データの取得

curl -H “X-CMC_PRO_API_KEY: ******** ” -H “Accept: application/json” -d “start=1&limit=200&convert=JPY” -G https://pro-api.coinmarketcap.com/v1/cryptocurrency/listings/latest > date "+%Y%m%d-%H%M%S".json

API_KEYは登録すると取得できます。(定期的に取得。今回は10分間隔で実行した)

Rで解析するためにcsvに変換

cat *.json | jq -r ‘.data[] | [.symbol, .id, .name, .quote.JPY.price, .quote.JPY.volume_24h, .quote.JPY.market_cap, .quote.JPY.last_updated] | @csv’ > data.csv

RでBTCとETHのデータを抽出(lt-が化けるので表示上=にした)

> df = read.csv(“data.csv”, header=F, sep=”,”)
> colnames(df) = c(“a”,”b”,”c”,”d”,”e”,”f”,”g”)
> btc = df[df$a==”BTC”,][c(4,7)]
> eth = df[df$a==”ETH”,][c(4,7)]
> plot(as.POSIXct(btc$g, format = “%Y-%m-%dT%H:%M:%S”), btc$d, type=”l”, col=rgb(0,0,1))
> par(new=T)
> plot(as.POSIXct(eth$g, format = “%Y-%m-%dT%H:%M:%S”), eth$d, type=”l”, col=rgb(1,0,0))

プロット結果(4時間ほどのデータ。縦軸の値はあえてそろえず)

実際のサイト上のデータ

データ数は少ないですが、だいたいうまく取得できているようです。
個人的に興味深いデータであることと、動きが適度にあることからデータサンプルとしてもってこいの素材かもしれません。

> a1<-read.table(“a1.csv”, sep=”,”)
> a2<-read.table(“a2.csv”, sep=”,”)
> b1<-read.table(“b1.csv”, sep=”,”)
> b2<-read.table(“b2.csv”, sep=”,”)

それぞれplotしてみます。

a1.csv

a2.csv

b1.csv

b2.csv

データ数を揃えて、単純にx,y別々に相関係数を測ってみました。

> dim(a1)
[1] 63 3
> dim(a2)
[1] 106 3
> dim(b1)
[1] 75 3
> dim(b2)
[1] 106 3

> a1x<-approx(a1[,1], a1[,2],n=100)
> a1y<-approx(a1[,1], a1[,3],n=100)
> b1x<-approx(b1[,1], b1[,2],n=100)
> b1y<-approx(b1[,1], b1[,3],n=100)

> a2x<-approx(a2[,1], a2[,2],n=100)
> a2y<-approx(a2[,1], a2[,3],n=100)
> b2x<-approx(b2[,1], b2[,2],n=100)
> b2y<-approx(b2[,1], b2[,3],n=100)

> cor(a1x$y, a2x$y)
[1] 0.9905204
> cor(a1y$y, a2y$y)
[1] 0.9929844
> cor(b1x$y, b2x$y)
[1] 0.8621487
> cor(b1y$y, b2y$y)
[1] 0.9720394
> cor(a1x$y, b1x$y)
[1] 0.8606694
> cor(a2x$y, b2x$y)
[1] 0.9010704
> cor(a1y$y, b1y$y)
[1] 0.8947821
> cor(a2y$y, b2y$y)
[1] 0.9067105

だいたいa1とa2,b1とb2は近い関係になっていますが、bのxだけ期待と違いました。
相関が近い中の識別には、もうすこし判定につかう要素を増やす必要がありそうです。

マウス操作のデータ取得には、Processingを使いました。

PrintWriter file;
void setup() {
  size(400, 400);
  background(255);
  file = createWriter("data.csv");
}
void draw() {}
int i = 0;
int cnt = 0;
void mouseDragged() {
  stroke(255-i, 255-i, 255-i);
  fill(0,0,100 + i, 20);
  ellipse(mouseX, mouseY, 40+i/10, 40+i/10);
  i += 3;
  file.println(cnt + "," + mouseX + "," + mouseY);
  cnt ++;
}
void mouseReleased(){
  file.flush();
  file.close();
}

「Stan超初心者入門 」

とてもわかりやすい説明で参考にさせていただきました。
他の書籍などで、モデルで推定したデータなのか、実際のサンプリングデータなのかどちらわかりづらかったのですが、二項分布binomial.stanのシンプルなコードを見て理解できました。
ここでは正規分布の確率モデルで推定をやってみました。

test.stan(Home Directory にrというディレクトリを作ってその下に格納した)
平均 0.3、標準偏差0.1のモデルを記述しpを推定

parameters{
  real<lower=0, upper=1> p;
}
model{
  p ~ normal(0.3, 0.1);
}

R Consloeでの操作

> library(rstan)
> fit<-stan("r/test.stan")
> stan_hist(fit, pars="p")
> pp<-rstan::extract(fit)$p
> hist(pp)

上がStanのヒストグラム、下がRのヒストグラムです。Stanの領域とRの領域は別なので変数名はpのままでよかったのですが、わかりやすくするために分けました。これでStanで生成したデータがRで使えることが確認できました。
特徴的なのはStanのモデル定義にある”〜”です。プログラムというのは具体的に物事を記述しますが、このような抽象的表現で記述されることが斬新でした。しかも動的言語でなく静的言語C++にコンパイルされているとのことです。
まだまだ入り口ですが、次に進めたいと思います。

R 4.1.0
RStudio 1.4.1717
R Commander 2.7-1

Cluster Dendrogram

前回は、コードを書きましたが、今回はR Commander を使うため、ノーコードで実行しました。

> library(Rcmdr)

で起動します。
データはエクセルで作成します。
sinカーブに乱数でノイズをのせたデータ4パターンと-0.5から0.5の乱数1バターンの5パターンを用意します。(100行:100項目の計測データがそれぞれあるという想定)

R Commanderでインポートするために、行と列を入れ替えてヘッダ(1から100)をつけます。(セルが更新されるためデータは変化します)

範囲選択して、コピーしてR commanderにクリップボード経由でインポートします。

次にクラスタ分析を実行します。

全ての変数を選びます。

オプションはデフォルトのまま。

結果は、系列1と2、系列3と4は振幅が違うたけで、近い関係なので、グループにされます。5は更新するたびにどちらのグループにつくかぶれます。

前回コードを書いたことと比べると、大量のデータを分類するときにとても便利です。

par(col.axis='firebrick2', col.lab='white', bg='gray33', fg='white')

m <- matrix(c(1,1,1,2,2,2,3,1,3,4,1,1,5,5,5,6,4,2),6,3, byrow=T)
rownames(m)<-c("A","B","C","D","E","F")
colnames(m)<-c("X","Y","Z")
plot(data.frame(m),col='yellow')

d<-dist(m)
hc = hclust(d)
plot(hc)

mがサンプルデータマトリックス、dがデータ間の距離を表しています。
距離はX,Y,Zそれぞれの距離の2乗和の平方になります。
B-Dなら、

(2-4)^2 + (2-1)^2 + (2-1) ^2 = 6
√6 = 2.449489…

この方法で導き出した距離をもとにグルーピングされます。

ついでに散布図も。

par(col.axis='firebrick2', col.lab='white', bg='gray33', fg='white')

f2 <- function(x,y){predict(lm(y~x + I(x^2)))}
f3 <- function(x,y){predict(lm(y~x + I(x^2) +  I(x^3)))}
f4 <- function(x,y){predict(lm(y~x + I(x^2) +  I(x^3) + I(x^4)))}

x0 <- 1:100/100
y0 <- 3 * sin(x0 * 10) + rnorm(x0)
plot(x0, y0)

lines(x0, f2(x0,y0), col="cyan")
lines(x0, f3(x0,y0), col="green")
lines(x0, f4(x0,y0), col="yellow")

対象データは、sin関数を正規分布で散らしたものを使いました。
これを2次、3次、4次曲線でフィッテングします。

以下は、f3と同じになります。係数を確認するときはこのように実行する必要があります。
確認のため、下に0.2ずらした位置にPlotしました。

ce <- coef(lm(y0~x0 + I(x0^2) +  I(x0^3)))
lines(x0, ce[1] + ce[2]*x0 + ce[3]*x0^2 + ce[4]*x0^3 - 0.2, col="blue")

その前にある区間の最大値、最小値を求めるRの関数についてのテストです。

par(col.axis='firebrick2', col.lab='white', bg='gray33', fg='white')
f <- function(x){x^3+(x-1)^2}
x = seq(-2,2,length=201)
t1 <- optimize(f, c(-2,0), maximum=TRUE)
t2 <- optimize(f, c(-1,2), maximum=FALSE)
plot(x,f(x), col=7, type='l')
abline(v=t1$maximum, col=4)
abline(v=t2$minimum, col=4)

三次関数の極大値、極小値を図示しています。optimizeを使うと、指定区間のこのような値を求めてくれます。

そして、尤度関数についても、最大値を求めてみます。

L <- function(p){p^3*(1-p)^7}
p <- seq(0,1, length=101)
t3 <- optimize(L, c(0, 1), maximum=TRUE)
plot(p, L(p), col=7, type='l')
abline(v=t3$maximum, col=4)

この試行は、コインを10回投げて、3回表がでたと想定します。
尤度関数が最大となる0.3が表がでる確率を表しています。
最大値は、モード(最尤値)として以下のようにも計算できます。
((3+1)-1)/((3+1)+(7+1)-2) = 0.3

このp^n(1-p)^nは、ベータ分布(ベルヌーイ分布)と呼ばれ、関数にしてみました。(係数は簡略化しています)

Beta <- function(x, a, b){x^a*(1-x)^b}
par(mfrow = c(3,1))
plot(p, Beta(p, 0, 0), col=7, type='l')
plot(p, Beta(p, 1, 1), col=7, type='l')
plot(p, Beta(p, 2, 0), col=7, type='l')

一番上は、一様分布を表します。
Beta関数の第二引数、第三引数は、上の例で言えば、表が出た回数と裏が出た回数に相当します。

そして三番目の分布を事前分布として、尤度をかけて事後分布をもとめると以下のようになります。

par(mfrow = c(2,1))
plot(p, Beta(p, 2, 0) * L(p), col=7, type='l')
plot(p, Beta(p, 2, 0) * L(p) * L(p), col=7, type='l')

事後分布を事前分布に更新(ベイズ更新)して繰り返し尤度をかけていくやり方は、マルコフの状態遷移行列にも似てるように感じました。

環境: R 3.0.0 / Windows 7

状態遷移行列の内容や掛け合わせる回数を変えながら、グラフの特性を確認できます。
各行の数値は遷移先への確率を表しています。自分自身への遷移は0で、他を合計すると1.0になるようにします。

以下実行コマンド

library('igraph')
g <- matrix(0,10,10)
g[1,] <- c(0, 0.3, 0.7, 0, 0, 0, 0, 0, 0, 0)
g[2,] <- c(0.1, 0, 0, 0.4, 0.5, 0, 0, 0, 0, 0)
g[3,] <- c(0.7, 0.3, 0, 0, 0, 0, 0, 0, 0, 0)
g[4,] <- c(0, 0, 0, 0, 0, 0.3, 0.3, 0.4, 0, 0)
g[5,] <- c(0, 0.2, 0.2, 0.2, 0, 0, 0, 0.4, 0, 0)
g[6,] <- c(0, 0, 0, 0, 0, 0, 0.3, 0.7, 0, 0)
g[7,] <- c(0.3, 0.1, 0, 0, 0, 0.5, 0, 0.1, 0, 0)
g[8,] <- c(0, 0, 0, 0, 0.3, 0.3, 0.3, 0, 0, 0.1)
g[9,] <- c(0, 0, 0, 0.2, 0, 0, 0.2, 0.4, 0, 0.2)
g[10,] <- c(0, 0, 0, 0.0, 0, 0, 0.2, 0.4, 0.4, 0)

plot(graph.adjacency(g>0),layout=layout.circle)

colors = c("red", "yellow", "blue", "green", "orange", "cyan", "hotpink", "green4", "coral", "azure")
barplot(t(g), legend=c('1','2','3','4','5','6','7','8','9','10'), col=colors)
barplot(t(g%*%g%*%g%*%g%*%g), legend=c('1','2','3','4','5','6','7','8','9','10'), col=colors)
image(g, col=topo.colors(100))
image(g%*%g%*%g%*%g%*%g, col=topo.colors(100))

※「%*%」演算子は、行列の積をとるものです。

使用データ

初期状態(左)と5回遷移後(右)の確率を並べました。
上の「使用データ」と見比べるとイメージがしやすいです。

image関数で表示

おもしろいデータ作成になかなか時間がかかりそうです。もっと効率のいい方法をまた考えたいと思います。

環境 : R version 3.0.0 / Windows 7
初期設定

par(mfrow = c(3,1), col.axis='firebrick2', col.lab='white', bg='gray33', fg='white', col.main='yellow', col.sub='yellow')
s <- sin((1:1000)*0.1)
r <- runif(1000, 0, 1)

ランダム

plot(r)
acf(r)
spectrum(r)

サイン波

plot(s)
acf(s)
spectrum(s)

ランダム+サイン波

plot(0.8*r+0.2*s)
acf(0.8*r+0.2*s)
spectrum(0.8*r+0.2*s)

データの周期性を調べるときや、その他Plot関連のメモでした。

環境: R version 3.0 / Windows 7

対象のデータは、三次曲線にノイズを負荷して、散らしたものを使いました。

x=seq(1,10,0.1)
y=jitter(x^3,amount=50)
df = data.frame(x,y)
f1 <- function(data,par) with(data, sum((par[1] + par[2] * x - y) ^ 2))
f2 <- function(data,par) with(data, sum((par[1] + par[2] * x ^ 2 - y) ^ 2))
f3 <- function(data,par) with(data, sum((par[1] + par[2] * x ^ 3 - y) ^ 2))
res1 <- optim(par = c(0,1), f1, data = df)
res2 <- optim(par = c(0,1), f2, data = df)
res3 <- optim(par = c(0,1), f3, data = df)
plot(df, xlim=c(0, 10), ylim=c(-100,1000), ylab='y');
par(new='T');plot(x,res1$par[1] + res1$par[2] * x, type='l', col='red', xlim=c(0, 10), ylim=c(-100,1000), ylab='y')
par(new='T');plot(x,res2$par[1] + res2$par[2] * x ^ 2, type='l', col='blue', xlim=c(0, 10), ylim=c(-100,1000), ylab='y')
par(new='T');plot(x,res3$par[1] + res3$par[2] * x ^ 3, type='l', col='green', xlim=c(0, 10), ylim=c(-100,1000), ylab='y')

一次、二次、三次それぞれについてデータにフィッテングさせてみます。
これは「関数」に対して最適な係数を探し出すものでとても強力で興味深い機能です。
(元のデータをグラフで見て、適用できそうな関数を選ぶ)

以下Rのコンソール。係数を確認。

参考サイト:
How to use optim in R