チラシの裏の落書き日記

統計とか,研究とか,日常についての備忘録的なもの。

自分用Julia言語勉強のためのリンク集(随時更新)

全般的な入門
rishida.hatenablog.com

https://nbviewer.jupyter.org/github/bicycle1885/Julia-Tutorial/blob/master/Julia高速チュートリアル.ipynb

The Fast Track to Julia

www.geocities.jp

Linear algebra · The Julia Language


行列の使い方
行列 |


Juliaでの確率分布の使い方
分布 |

文字列の処理(桁の0つめ)とか,地味に便利
bicycle1885.hatenablog.com

Rも使えるっぽいので,併用する。
JuliaからRを使う - りんごがでている


Rの関数っぽいものをJuliaでいろいろ
expand.grid
Function like expand.grid in R - Usage - JuliaLang

Rの中でJuliaの関数を書く
Writing Julia functions in R with examplesinsightr.wordpress.com



IJulia notebookを起動するために
using IJulia
notebook(detached =true)

問題を解決するために何を考えるのか

 問題を解く前に考えるべきこと。第1に何を目的にしているのかを明確にしなければならない。大きな目標からしたら,現在目の前にある”問題のようなもの”は本当の問題ではないかもしれない。実際には取り組んだり解決する必要がないものの可能性がある。もっとも大きい目標に照らして何が解決されれば,目標が達成できるのか逆算的に考える。そのなかでベストな道筋を探し,そこで出会う問題を一つ一つ解決していけばよい。

 問題の定義は目標との関係で変わりうることを念頭において,今何の問題に回答しようとしているのか,正確に見定める。目標達成のために解決する必要がない問題であれば,さっそうと無視する。取り組まなければならない問題に当たった場合には,その問題を5W1Hを明確にして述べる。文章にして書いてみるのがよい。固定された言葉は物理的なものと同様に取り扱うことができる。話し言葉で宣言するだけではふわふわしたまま変化してしまう。一度紙などに書いて現状を固定する事が必要となる。書いてみた問題は,抽象的すぎないか具体的すぎないかをチェックする。抽象的すぎると実際の行動に移すことができないし解決できない。具体的すぎると解決する意味がない問題になってしまう可能性もある。このあたりは非常に感覚的なものなので,人に相談するのもよいかもしれない。目指すべきは具体的な対象がイメージ出来ながらできるだけ一般性が高い状況を設定すること。個人的には自分の考え方の癖として抽象的すぎることが多いので具体性を高めに考えるように注意をしている。

 解を出すべき問題が見定められたら,それをさらに分解する。具体的に何が達成されれば,その問題が解決されたことになるのか,より細かい問題に分割していく。このサブ問題は大きな問題よりもより具体性が高いので,実感が持ちやすい。ポイントなのは,ストーリーとしてなぜそれらのサブ問題が必要になるのか,ということ。大きな目標を設定したときに解くべき問題がどのように定められるのかもストーリーに依存するが,細かいサブの問題にも同様にストーリーが必要となる。
なんとなくそれぞれのサブ問題が大事だから,というような思考ではいけない。それぞれを解くことが何を意味しているのか,位置付けと関係を述べることが必要になる。

 さらに,このサブ問題を解決できる具体的な作業に落とし込んでいく。理想的にはほとんど何も考えずにできるような形にできるのがよい。コンピュータにでもなったつもりで,考えないでもできる作業まで落とし込めたら,非常によい。あとは,ひたすら作業をしていくだけ。
このあたりは自分の得意な型に引き込めるとより速く,確実な解が得られる。どうしても必要であれば,新しい知識をみにつけることもいとわない姿勢は必要。目的がしっかりしていれば,学びは非常に早いものだ。

 もちろん途中で予期せぬ自体に直面することはあるが,全体のストーリやサブ問題の位置付けがしっかりしていれば,修正は容易になる。

 こうした目標の設定,問題の定式化,ストーリーの設定,問題の分解,作業の実行,という一連の営みを何度も繰り返してプロセス全体の作業に習熟するで,何が良い問題なのかなのかがわかってくるのだと思う。自分の頭で考えるということはこうしたことの積み重ねと繰り返しなのだと感じる今日このごろ。

研究について,最近思うこと

 研究を行う上で非常に重要なのは,やはりどのような問題に決着をつけるべきなのかその問題を見極めることだと思う。問題を設定するために,はじめは指導教員や先輩の研究者などに教えを請うのが望ましい。それによって,何が研究としてあるべき形なのかを身につけることができると考えられる。
 もちろん,自分で地道に考えていくことで研究関心を見つけていくのもなしではないが,効率が悪いのみならず外した後のリスクが大きくなるように感じられる。また,自分一人で考えたことは自分にこだわりがあまりにも強く出過ぎてしまう可能性もある。
 自分の問題ではなく,その分野の発展に資する問題関心を理解することが必要だろう。ただし,あまりに王道の研究関心は競争相手も非常に多く,マンパワー金のパワーの勝負になりやすい。そういう研究室に所属していて,指導教員の研究の一端を担うならばそれでもいいかもしれないが,個人の研究としてことを進めるのであれば,それはおすすめできない。
 狙うべきはニッチな領域である。ニッチな領域を押さえてそこでてっぺんを取ることが必要。人間はやはり成果が出ないことをやり続けることは難しいが,ニッチなことであれば相対的に成果を出しやすくモチベーションも維持しやすい。モチベーションを維持し続けられるような研究関心にしつつ,成果にも繋がる部分というものを見定めるためには,他人の力も借りて自分の考えを相対化する必要があろう。
 また,ここまで未定義に使ってきた問題という言葉であるが,何が「問題」なのかは注意して考える必要がある。個人的には「今解決されないと多くの人が困ること」を問題とするのがわかりやすいように感じている。(もちろん,2つ以上の集団で決着がついていない根本問題などを扱ってもいいかもしれないが,これを扱うのは必ずしも容易ではない。)もう少しラフに,「今こういうことがわかるとみんながハッピーになること」というような点に注目すると良いかと思う。
ちなみに,「今」というのもポイントだ。ある問題について今解決する必要がないというのは緊急性が低く,人に重要性を伝えられない。ある問題が解決されればよいけど今は別の問題の方が大事,ということもある。問題は一過性のものかもしれないが,「今」だからこそ解決する価値があり,時が立てば意味がなくなる問題も往々にして存在するということには注意したい。
 ちなみに自分の専門の統計学に近づけていえば,統計を使う人の役に立つこと,というのが研究としては大切なのではないかと思っている。新しい手法を提案して,新しい知見を生み出す道具を提供し続けること,今までの道具の不備を解消してより使いやすい方法を提案することなどなどやるべきことがいくらでもある。統計の特定の分野に関していえば,こうした発想の勝負として研究テーマは数多くでてくるように思われる。
 さて,問題を「今解決するとうれしいこと」とみなすと問題設定と同程度に重要な「ストーリー」を組み立てるのが容易になる。論文の導入ではなぜその研究を行う必要があるのか,特定の問題関心についての前提知識を読者と共有し,問題が解決されないとどれほど良くないことがあるのかを示す。それが伝わればあとは妥当は方法論を使って研究を遂行すれば良い。このストーリーが伝わらないと厄介である。
 なぜその研究を行うのか,どう嬉しいのかという説得を行うことができなければ話を理解することができない。その意味で,論文を書くことは一種の説得活動であるともいえる。通常の説得と違うのは単なる感情論での説得は意味をなさず,先行研究を使って論理的に説得するということである。
 論理的な説得を行うためには問題関心に関連した大きなメッセージ,その論文を読んだ人がどのような行動をしてほしいのかといった事がいえなければならない。この大きなメッセージのために,先行研究を積上げ,必然的に疑問に至らざるを得ないように組み立てる。
 どのような「問題」であれば重要とみなされるのか,どのようなストーリーであれば説得的か,論理的飛躍が小さいか,といった点は少なからず感覚的なものにならざるをえない。(論理的な飛躍は,数学などの分野では当然非常に厳密だが。。。)なので,必要なのはそういうことが上手い人の仕事を間近で見ることが最大の学びになる。もちろん,自分で試行錯誤し査読のコメントから学ぶということも多くある。しかし,王道の展開やよくある論理のパタンなどをしっかり身につけるのも感覚がない人には難しい(自分がかなりそう)。
 いつでもどこでも「問題」,「ストーリー」,「説得」ということは頭に入れて研究を行いたいものだ。このあたりについて,もう少し考えが深まって体系的な学習方法や教育方法がわかってきたら,まとまった形にして残しておきたい。

(殴り書き終了)

JAGSで予測分布の作成その2

IRTモデルで予測分布や対数尤度の計算などをおこなってみる。
対数密度の記述は以下の頁の質問に少し書かれている。
https://sourceforge.net/p/mcmc-jags/discussion/610037/thread/c728e87b/

The general rule is that if y ~ dfoo(alpha, beta, gamma, ...) then the contribution to the log likelihood of y is given by lik[i] <- logdensity.foo(y, alpha, beta, gamma, ...).

とのことなので,これに従って書いてみよう。(binomial分布はlogdensity.binなので注意)
また,bugsモジュールも読み込んでおく。
(どうでもいいが,T(0,)という関数はサンプリングする範囲を制限するために使うっぽいので,活用すること。)

#
# IRTモデルでの予測分布の作成など。 
#
library("irtoys")
library("jagsUI")
# help(package="jagsUI")
# help(package="irtoys")
#
# 適当な真値
#
Scored2pl

# 適当なデータ
data <- list(x = Scored,
             I = nrow(Scored),
             J = ncol(Scored)
)


# JAGS code 
cat("model{
    # likelihood
    for(i in 1:I){
      for(j in 1:J){
        x[i, j] ~ dbern(prob[i, j])
        logit(prob[i, j]) <- alpha[j] *(theta[i] - delta[j])
        # 予測分布
        x_pred[i, j] ~ dbern(prob[i, j])
        
        # 個々人の個々の項目の尤度
        logLik[i,j]  <- logdensity.bern(x[i,j],prob[i,j])
        #logLik_temp[i,j]  <- logdensity.bern(x[i,j],prob[i,j])
      }
      #logLik[i] <- sum(logLik_temp[i,])

      # prior
      theta[i] ~ dnorm(0.0, 1.0)
    }
    
  # prior
    for(j in 1:J){
      delta[j] ~ dnorm(0 , pr.delta)
      alpha[j] ~ dnorm(0 , pr.alpha)T(0,)
  }
    pr.delta <- pow(100 , -2)
    pr.alpha <- pow(100 , -2)
    }", file="irt_2pl.jag")

fit_jag <- jags(data=data,
                parameters.to.save=c("alpha", "delta","x_pred","logLik","theta"), 
                model.file="irt_2pl.jag",
                modules=c("glm","bugs","basemod"),
                n.chains = 1,
                n.adapt = 0,
                n.iter = 2000,
                n.burnin = 1000,
                n.thin = 2)
fit_jag
Scored2pl$est

JAGSでの予測分布の作成

やり方は簡単で,予測分布を未知の変数としてモデルに書けばいい。
Stanではgenerated quantitiyブロックでrng系の関数を使えばいいが,JAGSではそのまま"x_pred ~ dnorm(mu, sigma)"などとするとRでいうところのx_pred <- rnom(1,mu, sigma)のように乱数を生成してくれる。
以下では母相関係数の推定を行ったあとに,2変量のデータの予測分布を作成してみる。

#
# 予測分布の作成
#

library("jagsUI")

#
# 相関係数については,以下の頁のものを参考にした。
# https://kyougokumakoto.blogspot.jp/2016/09/jags_91.html
#
library(MASS)
set.seed(1234) 
mu <- c(0,0) 
Sig <- matrix(c(1,0.75,0.75,1),2,2) 
x <- mvrnorm(100,mu,Sig) 
cor.test(x[,1],x[,2])

data = list(x = x, N = nrow(x)) 

# JAGS code 
cat("model{
    # likelihood
    for(i in 1:N){
    x[i,1:2] ~ dmnorm(mu[], prec[ , ]) 
    
    # 対数尤度
    logLik[i] <- logdensity.mnorm(x[i,1:2], mu, prec) 
    } 
    
    # 予測分布
    x_pred[1:2] ~ dmnorm(mu[], prec[ , ]) 

    # corresponding the precision matrix
    prec[1:2, 1:2] <- inverse(cov[,]) 
    
    # constructing the covariance matrix
    cov[1,1] <- sigma[1]^2
    cov[1,2] <- sigma[1]*sigma[2]*rho
    cov[2,1] <- sigma[1]*sigma[2]*rho
    cov[2,2] <- sigma[2]^2 
    
    #priors
    rho ~ dunif(-1,1) 
    mu[1] ~ dnorm(0, 1.0e-4) 
    mu[2] ~ dnorm(0, 1.0e-4)
    sigma[1] ~ dt(0, 1, 1)T(0,)
    sigma[2] ~ dt(0, 1, 1)T(0,) 
    
    }", file="cor_pred.jag")

fit_jag <- jags(data=data,
                parameters.to.save=c("mu", "sigma","rho","x_pred","logLik"), 
                model.file="cor_pred.jag",
                n.chains = 1,
                n.adapt = 100,
                n.iter = 2000,
                n.burnin = 1000,
                n.thin = 2)

# 予測分布の相関
cor(fit_jag$samples[[1]][,c("x_pred[1]","x_pred[2]")])
# 予測分布のヒストグラムの重ね合わせ
hist(fit_jag$samples[[1]][,c("x_pred[1]")],freq = F,col=rgb(1,0,1,alpha=0.5))
hist(x[,1],freq = F,add=T,col=rgb(1,0,0,alpha=0.5))