とげまる日記

九州の大学で数学の博士課程に通ってます。自分の日記として使っています。

Rのfor文でエラーが出ても処理を止めない方法

R

Rのfor文を回していると、途中でエラーが出て止まってしまいました。

その解決策としてエラーが出てもとりあえず最後まで回してくれるtry関数を紹介します。

以下のコードは10次元ベクトル a,bの成分をlog (rnorm(1,0,1))で生成。

log(rnorm(1,0,1))を計算する際、rnorm(1,0,1)は N(0,1)からのランダムサンプルであるためエラーを返すことがある。

以下ではエラーで止まることをtry関数で防ぎ、while関数でそれぞれ10個ずつ生成できるようにコードを書いた。

L <- 10
a <- numeric(L)
b <- numeric(L)
l <- 1

while(l <= L){
  try(a[l] <- log(rnorm(1,0,1)))
  if(class(a[l]) == "try-error"){
    a[l] <- NA
  }
  try(b[l] <- log(rnorm(1,0,1)))
  if(class(b[l]) == "try-error"){
    b[l] <- NA
  }
  if(is.na(a[l]) | is.na(b[l])){
    l <- l-1
  }
  l <- l+1
}
a
b

実行結果は以下のとおり。それぞれ10個正しく生成できました。

[1]  0.3591358  0.1832536 -3.6249145 -0.4584207  0.1969482 -1.8525264 -1.3818560 -1.3684374 -0.9294030 -2.8535764
 [1] -0.54935989 -0.45241981 -0.31944820  0.02388071 -1.25626751 -1.55962373 -1.34217871 -3.33795157 -3.33510165
[10] -0.52030770

お役に立てれば幸いです。

R言語で対数尤度関数を作り最適化を行う

R

Rで最尤法を行おうと思ったのですが、ネットで調べてもいまいち分からなかったため自分で作ってみました。 自然な設定として線形回帰モデルを考え、最尤推定量を求めます。

モデル：

\displaystyle y = X\mu + \epsilon

ここで y \in \mathbb{R}^n, X \in \mathbb{R}^n \times \mathbb{R}^q, \mu\in \mathbb{R}^q, \epsilon \in \mathbb{R}^n. 誤差項の分布として \epsilon \sim N(0, \sigma^{2})を仮定。

設定:

N <- 1000 # サンプル数
mu <- c(5,2,3) # 回帰係数
sig <- 3 # 分散
X <- matrix(runif(3*N,min=1,max=5), nr=N, nc=3) # 説明変数行列
par.init <- c(2,2,2,1) # 最適化の初期値

対数尤度関数:

\displaystyle \ell_n(\theta):= \sum_{j=1}^{n} \log \frac{1}{\sigma}f \left(\frac{y_j - X_j \mu}{\sigma}\right)

ここで f(x)は N(0,1)の確率密度関数、 y_j: y の j番目の成分、 X_j: Xの j行。

log.likelihood <- function(y){
  return(function(par){
    mu1 <- par[1]
    mu2 <- par[2]
    mu3 <- par[3]
    mu <- c(mu1, mu2, mu3)
    sig <- par[2]
    return(sum(log(dnorm((y-X%*%mu)/sig, 0,1)/sig)))
  })
}

推定:

~
t <- proc.time()
y0 <- X%*%mu + rnorm(N,0,sig)  
mle.sample <- optim(par.init, log.likelihood(y0), method = "L-BFGS-B", lower = c(-Inf,-Inf,-Inf,0), upper = c(Inf, Inf,Inf,Inf), control = list(fnscale = -1))$par
proc.time()-t
mle.sample
~

計算時間と推定結果

> proc.time()-t
   user  system elapsed 
  0.006   0.001   0.006 
> mle.sample
[1] 4.670545 2.559163 2.805950 1.000000

二項検定を具体例で考えてみた

初めましてPandaです!

先日大学の同期とテニスの試合をしたのですが、結果として7勝1敗で私の圧勝となりました。

嬉しかった私は「俺の方が強かったね」と同期に言ったのですが、彼は「まだサンプル数が少ないから本当に強いかどうか分からない」と負けを認めませんでした。

何としてでも自分が強いことを示したかった私は、二項検定を用いて数学的に私の強さを証明することにしました。

まず問題設定から。 X: Pandaの勝利数を表す確率変数, p: Pandaの勝つ確率. この時立てる仮説は

H_0: p=0.5, H_1: p=0.7.

この時臨界値kを計算する。

P(X\geq k |n=8, p=0.5) < 0.95となる最小のkは

P(X\geq k |n=8, p=0.5)= \sum_{i=1}^{8} \binom{8}{i} p^{i} (1-p)^{8-i} < 0.95

を計算するとk=6であることがわかる. つまり私が6回以上勝つと帰無仮説を棄却できる。今回は7回勝ったので p=0.7を採択できた。

次にこの検定の信頼性を調べるために検出力を計算する。検出力とは帰無仮説が正しくないときに帰無仮説を棄却できる確率である。

Power = 1 - P(X&lt;5 | n=8, p=0.5) \fallingdotseq 0.552.

以上の結果から残念ながら現時点でのサンプル数では検出力が0.552であり、多くの場合検出力として設定される0.8に遠く及ばない。

そのため今後有意に私の方が強いと示せるようにテニスの試合数を増やし勝利数を増やしていきたいと考えている。

分布収束と確率収束の関係

この記事の内容

確率変数の分布収束について考えた時,X_n \xrightarrow{L} Xが成り立つ時、Xを移項したような形のX_n - X　\xrightarrow{L} 0が成り立つのか気になったので、確率変数の収束の性質について書いてみました.また,分布収束の部分を確率収束に変えた時の性質も書きました。最後に分布収束と確率収束の関係について書いています。

命題1　X_n \xrightarrow{L} X　 \nRightarrow　X_n - X　\xrightarrow{L} 0

証明

反例をあげれば良い。
X_n ~ N(0,1) , X~N(0,1)と仮定する.この時 , X_n  \xrightarrow{L} Xは明らかに成立し,仮定を満たす.
しかし,
Y_n := X_n-Xと定義しておくと正規分布の再生性から , Y_n ~ N(0,2)であり , Y_n  \xrightarrow{L} N(0,2)となり0に分布収束しない.

命題2　X_n-X \xrightarrow{L} 0 \Rightarrow X_n \xrightarrow{L} X

証明

　
X_n-X \xrightarrow{L}  0  \Leftrightarrow X_n \xrightarrow{p}であることから , XとX_nは漸近同値である.
漸近同値な分布は同じ分布に分布収束するため、X_n  \xrightarrow{L} Xである.

注意　命題2の証明中に漸近同値な分布は同じ分布に分布収束するという定理を用いたがその定理は必携統計的大標本論の定理6(b) にのっています。

命題3　X_n \xrightarrow{p} X \Leftrightarrow X_n - X \xrightarrow{L} 0

証明

これは確率収束の定義から明らか.



以上の性質をまとめた関係図を載せてこの記事を終わりにする.

f:id:napton:20210113020622j:plain
分布収束と確率収束の関係

正規分布の偶数次モーメントの求め方

正規分布の偶数次モーメントは公式で覚えておくと便利なので、その公式と導出を書いてみました.

X~N(0, \sigma) \Rightarrow E[X2k]=(2k-1)!! \sigma^{2k} (kは自然数)

証明


E[X^{2k}] \\\\
=\int_{-\infty}^{\infty} x^{2k}\frac{1}{\sqrt{2\pi\sigma^2}}exp(-\frac{x^2}{2\sigma^2})dx\\\\
=\frac{1}{\sqrt{2\pi\sigma^2}}\int_{-\infty}^{\infty}x^{2k} exp(-\frac{x^2}{2\sigma^2})dx\\\\
=\frac{1}{\sqrt{2\pi\sigma^2}}\int_{0}^{\infty}x^{2k} exp(-\frac{x^2}{2\sigma^2})dx　(\because非積分関数が偶関数)\\\\
=\frac{(2\sigma^2)^k}{\sqrt{\pi}}\int_{0}^{\infty} u^{k-\frac{1}{2}} e^{-u}du　(\because u=\frac{x^2}{2\sigma^2}と変数変換)\\\\
=\frac{(2\sigma^2)^k}{\sqrt{\pi}}\Gamma(k+\frac{1}{2})　(\becauseガンマ関数の定義から)\\\\
=\frac{(2\sigma^2)^k}{\sqrt{\pi}}\frac{(2k-1)!!}{2^{k}}\Gamma(\frac{1}{2})　(\because ガンマ関数の性質より)\\\\
=(2k-1)!!\sigma^{2k}

参考文献

1.ガンマ関数の計算に関しては以下のサイトを参考にした. https://ja.wikipedia.org/wiki/%E3%82%AC%E3%83%B3%E3%83%9E%E9%96%A2%E6%95%B0

最尤推定量の漸近正規性から一致性を導く方法

この記事の内容

漸近正規性を持つ最尤推定量は一致性を持つという記述を渡辺澄夫さんのサイト http://watanabe-www.math.dis.titech.ac.jp/users/swatanab/ML_cons2.pdf で見て、なぜ成り立つのかピンとこなかったのでその証明を書いてみました .

本題を示す前に , まず補題を示します .

補題の主張

実数列 (a_n) が a_n \xrightarrow{n} a で , X_n \xrightarrow{L} Xであるとき ,

a_nX_n \xrightarrow{L} aX

補題の証明

分布収束と特性関数の収束は連続定理より , 同値であるため特性関数の収束を示す . \varphi(t),\varphi_{n}(t) をそれぞれX,X_n の特性関数とする.


\varphi_{a_n X_n}(t)\\\\
= E[exp(ita_n X_n)] \\\\
= exp(a_n)E[exp(itX_n)] \\\\
= exp(a_n)\varphi \_{n}(t) \\\\
\xrightarrow{n} exp(a)\varphi(t)\\\\
= \varphi\_{aX}(t)\\\\

ここから本題を示していきます.\

主張

\sqrt{n}(\hat{\theta_{n}}-\theta) \xrightarrow{L} N(0,1) \Rightarrow (\hat{\theta_{n}}-\theta) \xrightarrow{L} 0

証明

X_n=\sqrt{n}(\hat{\theta_{n}}-\theta)とおく . このとき仮定より , X_n \xrightarrow{L} N(0,1)が成立. また,nを無限大に飛ばすと, \frac{1}{\sqrt{n}} \xrightarrow{} 0が成立 .

よって補題より \frac{1}{\sqrt{n}}X_n \xrightarrow{L} N(0,1)　つまり , (\hat{\theta_{n}}-\theta) \xrightarrow{L} 0が成立 .

補足　

補題の定理は確率収束や概収束の場合も成り立つ . 平均r次収束の場合は確認していない .