統計科学

データ増やした時の標準誤差の減り具合

掲題の件、これはふつー「1/√データ数」ですが、適当にデータとして0だけ突っ込んでいったらどうなるんじゃいと思ってやってみたら大体似たようなもんだった。 理論的な話はあとで考える。 > size_unit <- 10^3 > x <- c() > y <- c() > for (i in 0:9){ + …

N個の独立な一様乱数(-1, 1)の和の標準偏差は1/3√N

手計算で確認&ボケ防止にRでチェック。 sizes <- 2**(1:15) y <- numeric(length(sizes)) for(i in seq_along(sizes)){ x <- numeric(10**3) for(j in seq_len(10**3)){ x[j] <- sum(runif(sizes[i], min=-1, max=1)) } y[i] <- sd(x) } plot(sizes, y) lin…

purrr::map_dfr = lapply + dplyr::bind_rows

そういうことなんだよな〜シミュレーション系でよく使うのでメモ。 > dplyr::bind_rows(lapply(1:3, function(x){head(iris, 1)})) Sepal.Length Sepal.Width Petal.Length Petal.Width Species 1 5.1 3.5 1.4 0.2 setosa 2 5.1 3.5 1.4 0.2 setosa 3 5.1 3.…

確率変数を変数変換した場合の確率分布

昔やった気もするが、すぐ忘れる&頭の体操もかねてもう一度 算数 適当な確率分布に従う確率変数x(面倒なので[0, 1]区間の一様分布想定)とそれを2乗した変数yを考える。 このときyの従う確率分布はより となる。 Rでやる 0, 1区間の一様な乱数x, 及びそれを…

金融市場が混乱する間隔を日経平均ボラティリティー・インデックス(日経平均VI)から推定したい

日経平均ボラティリティー・インデックス(日経平均VI)の半減期が知りたい - My Life as a Mock Quant 日経平均ボラティリティー・インデックス(日経平均VI)を高速フーリエ変換したが、なんもなかった - My Life as a Mock Quant 日経ボラティリティー・…

確率変数の変換について

ロジックメモ ある確率変数とその写像で定義される2つの確率変数、特に確率変数の確率分布について考えたい。 の確率分布関数は と確率変数の確率分布関数を用いてあらわすことができる。ここで、は単調増加な関数であると仮定している。単調減少の場合は不…

売上げを気温で回帰してみる

あたまだし 掲題の件、最近出たわかりやすいデータ分析の入門書、 において 売上げ v.s. 気温 として回帰分析をしていたケースがあった。これが時系列データに対するやってはいけない回帰、すなわち、みせかけの回帰に相当するケースではないか?という話を…

標準誤差と(ノンパラメトリック)ブートストラップのイケない関係

標準誤差と(ノンパラメトリック)ブートストラップについて、特に理論的なお話はおいておいて、直感的にどういうことなのかということをメモりたい。 それぞれ (回帰分析の係数の)標準誤差: ある確率分布からのサンプルデータに対し推定されたある統計量(ここ…

可視化で理解しない「負の二項分布」

匿名知的集団であるホクソエムの親分が可視化で理解する「負の二項分布」 - ほくそ笑むで言及している負の二項分布ですが、シンクロニシティというか、同じ匿名知的集団ホクソエムなので当たり前と言えば当たり前なのですが、私も負の二項分布、特に負の二項…

層別サンプリングした結果を綺麗に持ちたい

層別サンプリングをするには、samplingパッケージのstrata関数使えばいいんだけど、出力とインターフェイスがイケてないのでラップする。 下記の`stratified_sampling`関数を使えば、各層(クラス)の割合を一定に保ったままサンプリングしてくれる。 library(…

(g)lmのNA省略(デフォルトの挙動)でうっかりうっかり

データにうっかり欠損(NA)が入ってると、デフォルトでその行がモデル構築&予測結果からごっそり削られるうっかりミスを防止するための備忘録うっかりうっかり。 > library(dplyr) > #デフォルト(自分で設定したオプション値だけど)のNAなし時の挙動 > get…

可視化で理解する中心極限定理(ベルヌーイ分布編)

こういう話がある。 可視化で理解する中心極限定理 #rstatsj 非常に素晴らしいので、指数分布じゃなくて、ベルヌーイ分布版をアニメーションにしてみた。 下図は ベルヌーイ分布(コイン投げでいうところの表が出る確率p=0.2)からサンプリングしたデータのヒ…

ガウシアンを一様乱数でモンテカルロ積分した際の一様乱数の幅依存性

ガウシアンを、一様乱数でモンテカルロ積分した際の結果の、一様乱数の範囲依存性がみたい。こんな感じか。 f <- function(x){exp(-0.5*x^2)} N <- 10^3 L <- c(1:10, 20, 50, 100) y <- numeric(length(L)) for(i in 1:length(L)){ y[i] <- L[i]/N*sum(f(ru…

時系列データに対するブートストラップ法(ブロック・ブートストラップ法)について

あたまだし 検定やクロスバリデーション等への応用を企図した、サンプル数を水増しするための手法としてブートストラップ法がある。これをRで実行するにはsample関数を使って自分でリサンプリングするコードを実装するか、あるいはbootパッケージのboot関数…

リサンプリング(復元抽出)で積分値を評価する

復元抽出のアルゴリズム - My Life as a Mock Quant の応用かつ、「これで間違ってねぇよな?」という自分への備忘録。インポータンスサンプリングや測度変換とは微妙に分母の形式が違ってくるので間違いやすいんだ。 数式の整理 以下の積分I を計算すること…

復元抽出のアルゴリズム

もくてき 粒子フィルタ(パーティクル・フィルタ)を実行する際には、粒子のウェイト(weight)に比例する確率でリサンプリングを実行する必要がある。そのためのアルゴリズムとコードを考えたい。ここでは手元にある各粒子のウェイトはK個の要素からなるベク…

非同次/確率的な強度λを持つ指数分布に従う確率変数τ(死亡/消滅時刻)の生成

非同次な指数分布 非同次なポアソン分布ってのがある。日本語で定義がちゃんと載ってるページがなかったので、英語版のwikipediaを参照すると Poisson point process - Wikipedia が該当する。要するにこれは、今まで、一定だと思っていたポアソン分布の強度…

やってみよう分析!F#で強化学習(Q-learning, ε-greedy行動選択)

やってみよう分析!Rで強化学習(Q-learning, ε-greedy行動選択) - My Life as a Mock Quant やってみよう分析!おまけ 2 - 1: Excel VBAで強化学習(Q-learning, ε-greedy / softmax 行動選択) のF#実装版。強化学習自体の解説は上の記事読んどいたらいい。…

(正規)分布の確率密度関数を乱数から推定する

すごい基本的なお話な気がするけど、数回は書き直しているのでメモっておく。ここでは 正規分布に従う乱数を生成して、それがある幅(dx)の中に入っているか否かをカウントして確率にする 正規分布の確率密度関数に幅(dx)をかけて確率を計算する の2つの方法…

リカレント・ニューラル・ネットワーク(Recurrent Neural Network, RNN)っていう、時系列に対するDeep Learningやりたい

はじめに 前々からDeep Learning系、特にその時系列への応用に興味があって、それにはリカレント・ニューラル・ネットワーク(Recurrent Neural Network, 以下RNN)ってのを理解する必要があるのは知っていたが、@yamano氏から@teramonagi こちらです。まずは…

二次元ヒートカーネル(熱核)の時間発展の可視化

二次元ヒートカーネルなんていうと難しく聞こえるけど、要するにラプラス方程式の解である、時間発展するガウシアンの可視化です。物理系だと微分方程式の基本解という意味で、熱核じゃなくてグリーン関数って呼んでるものと等価(なはず)。 熱核 - Wikiped…

円内に一様分布する乱数を生成する時は、俺、絶対忘れないよヤコビアンのこと

頭出し ある半径Rの円内に一様分布する乱数を生成する時には注意しないといけないことがありますよというお話。所謂「一度はやってしまうミス」系でもある。この手の話は円に限ったわけではなく、円の高次元版である球、あるいは超球(次元>3)、あるいは任意…

逆ガンマ分布(事後分布)〜正規分布(尤度関数)×逆ガンマ分布(事前分布)

一般に、ある確率測度の下で、逆ガンマ分布に従う確率変数の確率密度関数は である。この時、確率変数を分散に持つ正規分布を尤度関数とし、逆ガンマ分布を事前分布とすると、 確率変数の事後確率密度関数は 従って、 となり、のパラメーターが決定され、か…

正規分布(事後分布)〜正規分布(尤度関数)×正規分布(事前分布)

一般に、ある確率測度の下で、正規分布に従う確率変数に対し、 尤度関数: 事前分布: であることを仮定したとき、事前分布が共役事前分布となっていることから、事後分布も正規分布となり、 が成立する。ただしここで、は平均, 分散を持つ正規分布を表す。…

やってみよう分析!Rで強化学習(Q-learning, ε-greedy行動選択)

やってみよう分析!おまけ 2 - 1: Excel VBAで強化学習(Q-learning, ε-greedy / softmax 行動選択) のR実装版。強化学習自体の解説は上の記事読んどいたらいい。めんどいのでとりあえずε-greedyのみやった。計算結果は > Qlearning() [,1] [,2] [1,] 5 20.0…

逐次モンテカルロ/(粒子|パーティクル|モンテカルロ)フィルタを実装してみた

このポストの概要 逐次モンテカルロ法 (Sequential Monte Carlo: SMC)/粒子フィルタ(Particle Filter)は、状態空間モデルの状態変数の系列を観測値の系列 のみから推定するアルゴリズム パッケージを使ってもいいんだど、お勉強も兼ねて自分でコード書いてみ…

ブラウン橋(Brownian bridge)作りたい

簡単なんだけど、やってなかったなって。ブラウン橋言ってる人もいるし、ブラウニアン橋いってる人もいてどっちが正しいのかわかんねぇけど、ここではブラウン橋でいく。ブラウン橋は、を標準ブラウン運動として と書かれる確率過程だと思っておけばいい。こ…

AR(1)モデルのメモ

今更感あるが、諸事情によりメモ。やりたいことは(定常性を仮定した)AR(1)モデルの自己相関係数がちゃんと理論値と一致するのかを見たいって話。 AR(1)モデルとは AR(1)は、適当な確率過程 が、 に従うようなモデル。ここで、 は独立な平均0、標準偏差の正規…

scikit-learnでサポートベクトル回帰、及びそのパラメーター推計 with クロスバリデーションやってみる

サポートベクトル回帰(Support Vector Regression, SVR)の理論が大体脳内整理出来たので、実践もしたいぞと、そしてちょいとpythonを使う別件があるので、慣れの意味も込めてR言語ではなくpythonとその機械学習ライブラリであるscikit-learnを使ってやるぞと…

R言語でデータ拡大(Data Augmentation)サンプリング法を実装してみた

データ拡大 (Data Augmentation)サンプリング法ってのをお勉強したのでそのまとめ。 概要&アルゴリズム 説明は2つ目の参考LINKに書いてある方法でずばり良いと思う。 要するに 分布から直にサンプリングしたいんだけど、それが難しいという状況の時、適当…