2018-10-05

Julia1.0でマクロを@docしまくった

Julia

タイトルの通りである。

「Juliaのマクロ、全然わからん」

と思ったので、とりあえず@の後にtabを打って出てくるマクロすべてについて@docしてみた。

2018-08-10

【Julia】Julia-1.0.0のProjectを試す【仮想環境？】

Julia

Julia1.0.0が遂にリリースされました。

f:id:hotoke-X:20180810011534p:plain

仮想環境っぽいものを作って環境の切り替えができるようなので、試してみます。

Projectの作成

Anaconda Pythonではconda createコマンドで仮想環境を作れました。

Julia1.0.0ではProjectと呼ぶようです。

Anacondaでは

conda create -n envname

とする必要がありましたが、Juliaではプロジェクトのディレクトリに入って

activate .

とするようです。

その後、

julia>]
(projectname)pkg> add hogehoge

で、プロジェクトローカルな環境にhogehogeパッケージが入ります。

その後バックスペースでpkg環境を抜けて同フォルダで作業すれば良いっぽい。

ただ、Julia1.0.0がリリースされたのが記事作成時点で昨日の今日なので、対応が追い付いているパッケージを探すのが結構大変です（2018/08/10現在）。

私は情弱なのでおとなしく修正を待つことにします。

(ちなみに現状遊ぶなら0.6が無難)

2018-06-23

文献検索、管理、メモを1つの画面で完結させる

研究ライフハック

文献管理ソフトの紹介です。

以前にPaperpileが凄いって紹介をしました。 hotoke-x.hatenablog.com

個人的には好きなんですが、wordとの相性が悪い...

もっと良いのはないかと思ってたらReadcubeがありました。

数か月使ったので良い部分を簡単に紹介していきます。

※例示しているのはオープンアクセスな論文です。

Readcubeの凄いところ

複数の検索エンジンで同時検索
検索結果が自動保存される（そして簡単に消せる）
毎年の引用数がわかる
引用されている文献へ容易に飛べる
どんな風に引用されているのか見れる
Recommendationsが強力

こんな感じで全体的に便利すぎ。
全部1つの画面で完結します。
関連論文が次から次へ見つかって大変（嬉しい）。

使い方

ユーザー登録は済んでいるものとします。
また、ウェブ版での説明です（デスクトップ版もあります）。

基本画面（青い四角は身バレ防止です。）

f:id:hotoke-X:20180623035811p:plain

画面にNatureやらScienceやらGoogle Scholarやらが表示されてますが、こっから同時検索してくれます。なんかもうこれだけで凄い...。

f:id:hotoke-X:20180623035441p:plain

検索してみた結果。アブストが表示される。オープンアクセスだったり、学校などジャーナルを購読している機関からアクセスするとFigureまで表示してくれます。

f:id:hotoke-X:20180623035625p:plain

検索結果の図をクリックすると拡大表示。この段階ではまだ文献をダウンロードしたりしていません。あくまでも検索結果です。

「Add to Library」をクリックするとインポートしてくれます。
（supplementがあればそれも同時にインポートする。すごい...）

f:id:hotoke-X:20180623035836p:plain

そして、画面左側には検索結果が自動保存されます。便利！青い部分には自分でタグ付けした文献が入っています。一つの文献に複数のタグをもたせることもできます。

ここからが特にお気に入りの機能。

f:id:hotoke-X:20180623040308p:plain ライブラリで文献をクリックすると、情報が表示されます。アブストがここでざっと確認出来ます。

で、上のグラフのマークをクリックすると引用情報が見れます。 f:id:hotoke-X:20180623040530p:plain

毎年どれくらい引用されているのか、どの論文に引用されているのかを見ることができます。

f:id:hotoke-X:20180623040724p:plain

どんなふうに引用されているかも見ることができる。凄すぎて拍手。

その他

デスクトップ版ではwordへのcitation追加機能も付いてきます。
もちろん文献を読みながらハイライトしたり、ノートをとることも可能です。
デスクトップアプリが重いことが唯一残念。

Recommendationsをクリックすると、ライブラリの文献を元にお勧め論文を検索して表示してくれます。
しかもタグごとのお勧めも表示可能。
関連論文がめっちゃ見つかります。

2018-06-09

JuliaでPRML：ベイズ曲線フィッティング

Julia PRML

やってみた。

コードは

qiita.com

をかなり参考にさせてもらっています。計画行列が自力ではうまく使えないので若干コードに不安ありです。

2018-06-09

JuliaでPRML：多項式曲線フィッティング

Julia PRML

Juliaの勉強も兼ねてPRMLの多項式フィッティングを実装してみた。

誤差関数書いたけど、今回は解析的にパラメータが出てくるので結局使っていない。
GAとかMCMCとかでパラメータ推定して性能見てみるのもいいかも。
計画行列をうまく使えるようになりたい。

2018-03-22

【数理統計】ラオーブラックウェル化（Rao-Blackwellization）

統計

理解したいのはこっちでした。

本質的には同じ問題だけど。

ラオーブラックウェル化（Rao-Blackwellization）は、条件付き期待値を利用して期待値を計算する方法です。

サンプルが独立なら、ブラックウェルーラオの定理より直接期待値を計算するより精度が良くなる事が保証されます（不思議）。

名前の順番がひっくり返っているのは知りません（書籍でこう書いてあったからそのまま書いてます）

ブラックウェルーラオの定理の証明はこちら

hotoke-x.hatenablog.com

ラオーブラックウェル化（Rao-Blackwellization）

確率変数 $\boldsymbol{x}$ を二つのブロックに分割し、 $\boldsymbol{x} = (\boldsymbol{x_1}, \boldsymbol{x_2}) (x_i \in \chi_i)$ と表す。このとき、 $x_1$ の統計量の期待値

f:id:hotoke-X:20180322193210p:plain

を求めることを考えます。

$\pi (\boldsymbol{x_1}, \boldsymbol{x_2})$ からサンプリングされた $\left( \boldsymbol{x^{(1)}}, \ldots, \boldsymbol{x^{(t_0 + T)}} \right)$ を用いれば、上式の期待値は以下のように計算できます。

$\displaystyle \hat I = \frac{1}{T} \sum_{t=1}^{T} g(\boldsymbol{x_1^{(t_0 + T)}})$

ここで、 $\boldsymbol{x^{(t)}} = \left(\boldsymbol{x_1^{(t)}}, \boldsymbol{x_2^{(t)}} \right)$ であることに注意。

さらに $\boldsymbol{x_2}$ の条件付き分布で同時分布を分解すると

f:id:hotoke-X:20180322194353p:plain

となるので、この式の真ん中部分

$\displaystyle h(\boldsymbol{x_2}) = \boldsymbol{\mathrm{E}} \left[g(\boldsymbol{x_1}) | \boldsymbol{x_2} \right] = \int_{\chi_1} g(\boldsymbol{x_1})\pi(\boldsymbol{x_1} | \boldsymbol{x_2}) \mathrm{d}\boldsymbol{x_1}$

を解析的に求められれば、

$\displaystyle \hat I_{RB} = \frac{1}{T} \sum_{t=1}^{T} h(\boldsymbol{x_2}^{(t_0 + t)})$

によって計算できる。

なんだ面倒くさいと思うが、実はラオーブラックウェルの定理から

$\displaystyle \boldsymbol{\mathrm{V}} (\hat I) \geq \boldsymbol{\mathrm{V}} (\hat I_{RB})$

が成立し、より良い推定値になり得ることがわかる。

参考書籍

応用を目指す数理統計学（国友直人）
ベイズ計算統計学（古澄英男）

2018-03-22

【数理統計】ブラックウェルーラオ（Blackwell-Rao）の定理

統計

はい。

私がちゃんと理解したかった定理です。

復習したら見た瞬間理解したのですが、一応メモっておきます。

初学者向けの説明がなかなかネットに落ちていなかった（気がする）ので。

初めに行っておくと、この定理は

「十分統計量で条件づけた不偏推定量の分散は、他の不偏推定量の分散より大きくなることはない。」

って定理です（違ったらコメントください）。

十分統計量に依存した不偏推定量を使っとけばとりあえず無難ですよってことですな。

数学アレルギーでも、この程度の認識は持っておいた方が良いでしょう。

十分統計量についてはこちら hotoke-x.hatenablog.com

不偏推定量についてはこちら hotoke-x.hatenablog.com

証明

証明は結構簡単です。

今、未知母数の推定量は不変推定量しか考えないとします。このとき、標本 $\boldsymbol{X}$ について $\theta$ の不偏推定量 $\phi (\boldsymbol{X})$ 、十分統計量 $T(\boldsymbol{X})$ を考えます。この時十分統計量で条件付けした統計量を $\varphi (t) = \boldsymbol{\mathrm{E}} \left[\phi (\boldsymbol{X}) | T = t \right]$ とすると

$\boldsymbol{\mathrm{E}} \left[ \left( \phi (\boldsymbol{X}) - \theta \right)^2 \right] \\= \boldsymbol{\mathrm{E}} \left[ \left( \phi (\boldsymbol{X}) - \varphi (T) \right) + \left(\varphi (T) - \theta \right) \right]^2 \\= \boldsymbol{\mathrm{E}} \left[ \left( \phi (\boldsymbol{X}) - \varphi (T) \right)^2 \right] + \boldsymbol{\mathrm{E}} \left[ \left(\varphi (T) - \theta \right)^2 \right] \\ \geq \boldsymbol{\mathrm{V}} \left[ \varphi (T) \right]$

十分統計量に依存しない不偏推定量を用いた場合、下から2段目の式の第一項の分、推定量が悪くなり得ることを示しています。

まとめるとブラックウェルーラオ（Blackwell-Rao）の定理は以下のようになります[1]。

標本 $\boldsymbol{X} = (X_1, \ldots, X_n)'$ の各要素が互いに独立同分布に従うとき（independent and identically distributed, i.i.d.と書くことも多い）、 $T(\boldsymbol{X})$ を十分統計量、統計量 $\phi (\boldsymbol{X})$ は母数 $\theta$ の不偏推定量とする。このとき

$\varphi (t) = \boldsymbol{\mathrm{E}} \left[ \phi (\boldsymbol{X}) | T=t \right]$ は不変推定量
$\boldsymbol{\mathrm{V}} \left[ \varphi (T) | \theta \right] \leq \boldsymbol{\mathrm{V}} \left[\phi (\boldsymbol{X}) \right]$ となる

参考書籍

応用を目指す数理統計学（国友直人）

2018-03-22

【数理統計】最尤推定量とバイアス

統計

尤度

尤度については

hotoke-x.hatenablog.com

を参照。

例によって離散か連続かはあまり気にせず書きます。

尤度関数とは、母集団分布からの独立標本が得られたとき、確率関数 $p(x|\theta)$ から同時確率関数を母数 $\theta$ の関数と見た

$$L_n (\theta | x_1, \ldots, x_n) = \prod_{i=1}^{n} p(x_i | \theta) $$

のことでした[1]。

この尤度関数を最大にするような $\theta$ の値を最尤推定値といい、その推定量 $\hat \theta$ を最尤推定量（maximum likelihood estimator, MLE）と呼びます。私の読むレベルの文献では $\hat \theta_{ML}$ と書かれていることが多いです。

最尤推定

最尤推定すれば未知母数の推定もうまくいきそうですが本当にそうでしょうか。百聞は一見に如かずということで確かめてみます。

ここで、正規分布の平均と分散を未知として、最尤推定する問題を考えます。

母集団分布を平均 $\mu$ 、分散 $\sigma^2$ の正規分布 $\mathcal{N}(\mu, \sigma^2)$ とする。標本を独立にn個サンプリングし、観測値 $x_1, \ldots, x_n$ が得られたとするとその尤度関数は

$\displaystyle L_n (\mu, \sigma | x_1, \ldots, x_n) = \left( \frac{1}{\sqrt{2 \pi \sigma^2}} \right)^n \exp \left\{ -\frac{1}{2\sigma^{2}} \sum_{i=1}^{n} (x_i - \mu)^2 \right\}$

で与えられ、その対数尤度関数は

$\displaystyle l_n (\mu, \sigma) = -\frac{n}{2} \log (2\pi) - \frac{n}{2} \log \sigma^2 - \frac{1}{2\sigma^2} \sum_{i=1}^{n} (x_i - \mu)^2$

となります。

母数 $\mu$ で微分して、 $\mu$ について解けば母数 $\mu$ についての最尤推定量が得られます。

$\displaystyle \frac{\partial l_n}{\partial \mu} = \frac{1}{\sigma^2} \sum_{i=1}^{n} (x_i - \mu) = 0$

$\displaystyle \hat \mu_{ML} = \frac{1}{n} \sum_{i=1}^{n} x_i = \bar x$

この推定量を使って、分散も推定ができます。そしてこの推定量を使って他の推定量を得るという作業に落とし穴があるので注意（後述）。

対数尤度に $\mu$ の最尤推定量 $\hat \mu_{ML}$ を代入して母数 $\sigma$ で微分すると

$\displaystyle \frac{\partial l_n}{\partial \sigma} = - \frac{n}{2\sigma^2} + \frac{\sum_{i=4}^{n} (x_i - \hat \mu_{ML} )^2}{2\sigma^4} =0$

が得られ、 $\sigma$ について解くと

$\displaystyle \hat \sigma_{ML}^2 = \frac{1}{n} \sum_{i=1}^{n} (x_i - \bar x)^2$

が得られます。

以上の手順で、正規分布の母数（平均と分散）の最尤推定量

$\displaystyle \hat \mu_{ML} = \frac{1}{n} \sum_{i=1}^{n} x_i = \bar x$

$\displaystyle \hat \sigma_{ML}^2 = \frac{1}{n} \sum_{i=1}^{n} (x_i - \bar x)^2$

が得られました。

バイアス

先ほどの最尤推定量を得る過程で推定量を使って他の推定量を得るなんて少し怪しいことをしたわけですが、怪しいという直感は当たっていて、分散は不偏推定量になっていません。不偏というのは推定量が真の母数の周りに対称に分布することを言います。つまり、不偏でない推定量は母数を過小（あるいは過大）評価することを意味します。

推定量と母数の二乗誤差の期待値（mean square error, MSE）を計算して見ると

$\boldsymbol{\mathrm{E}} \left[ |\hat \theta - \theta|^2 \right] \\= \boldsymbol{\mathrm{E}} \left[ \left( (\hat \theta - \boldsymbol{\mathrm{E}} (\hat \theta)) + (\boldsymbol{\mathrm{E}} (\hat \theta) - \theta) \right)^2 \right] \\= \boldsymbol{\mathrm{E}} \left[ \left( \hat \theta - \boldsymbol{\mathrm{E}} (\hat \theta) \right)^2 \right] + \left[ \left( \boldsymbol{\mathrm{E}} (\hat \theta) - \theta \right)^2 \right] \\= \boldsymbol{\mathrm{V}} \left( \hat \theta \right) + \left( \boldsymbol{\mathrm{E}} (\hat \theta) - \theta \right)^2$

となり、右辺第二項の分推定量が偏ることがわかります。不偏推定量であればこのバイアス項は0になります。数学アレルギー勢にために一応説明しておくと、 $\boldsymbol{\mathrm{V}}$ は分散（Variance）という意味です。

先ほどの正規分布の例で、分散の最尤推定量が不偏推定量になっていないというのはこのバイアス項が0にならないことを意味しています。

それは困ったと思うかもしれませんが、どの程度偏っているかわかっているなら補正すればよいだけです。

分散の最尤推定量の期待値を計算してみると（計算は気が向くか、コメントいただければまた書きます（mathjax疲れた））

$\displaystyle \boldsymbol{\mathrm{E}} (\hat \sigma_{ML}^2) = \frac{(n-1)}{n} \sigma^2$

となり、 $(n-1)/n$ 倍だけ過小評価されている事がわかります。

じゃあ最尤推定値をこの分だけ補正して母数を推定すれば良いのです。

統計で習う分散の式が

$\displaystyle \frac{1}{n-1} \sum_{i=1}^{n} (x_i - \mu)^2$

となっているのはこの為です。

式だけでいきなり「ほぉそうか」と納得できる人はこれで良いのですが、数学アレルギー勢は釈然としないと思います。

絵で描くとわかりやすいので、気が向いたら絵でも描こうと思います。

今日は疲れたのでこの辺で。

参考書籍

応用を目指す数理統計学（国友直人）

2018-03-22

【数理統計】尤度関数と十分統計量

統計

勉強したことを自分なりにまとめたものです。
間違い等があれば指摘していただけると嬉しいです。

私は「せいきぶんぷ～？」、「のんぱらめとりっく～？」という状態から

応用を目指す数理統計学（国友直人）

で勉強しました。初学者が手を出すには敷居が高いですが、（本気で）頑張ればなんとかなるように書かれています。参考文献も充実しているのでとりあえず持っておくのはおすすめです。安いし。（アフィリエイト等はないので安心して買ってください。）

以下、応用を目指す数理統計学（国友直人）に基づいています。

自分なりの解釈も書いたりしてます。

Chapter 7 統計的推論より

統計量・推定量・推定値

統計量とは「標本の関数」のこと
推定量とは「統計量を使って母集団を表現している未知母数を推測するときの統計量のこと」
推定値とは「実際にデータを代入した場合の推定量のこと」

標本の関数でさえあれば統計量と呼ぶようです。その中でも、母数を推測する統計量を推定量と区別することに注意が必要といったところでしょうか。

「応用を目指す数理統計学」では正規分布に従うことの妥当性を標本積率を使って調べる方法を例に挙げてくれています（p108, 109）。（書きすぎるとただの本のコピペになるので、気になる方は購入してください）。母集団と標本の積率を対応づけるので積率法と呼ばれるものの一例のようです。最近では積率法を一般化した一般化積率法がよく使われているとのこと。10章で解説してくれているそうなのでまあたとで勉強しよう。

尤度関数と十分統計量

尤度関数とは「母数の関数」のこと
十分統計量とは「母数に関する情報損失が無い統計量」のこと

尤度関数

簡単のため、離散確率分布の場合を考えます（連続でも同様）。標本は大文字、確率変数は小文字で書くのがスタンダードっぽい。

$\theta$ を母数、 $x$ を標本の確率変数、 $p(x|\theta)$ を母集団が従う確率変数とする。

この時、 $n$ 個の標本 $X_1, X_2, \ldots , X_n$ の同時確率関数は

$$ p(x_1, x_2, \ldots , x_n | \theta) = \prod_{i=1}^{n} p(x_i|\theta) $$

で与えられる。これを、同時確率分布ではなく、母数 $\theta$ の関数として考えると

$$L_n (\theta | x_1, \ldots, x_n) = \prod_{i=1}^{n} p(x_i | \theta) $$

と書ける。これを尤度関数と呼ぶ。 $L$ はLikelihoodの頭文字。

得られた標本から計算された尤度関数の意味としては「起きた事象はどの程度起こりやすい事象だったのか」を定量化したことになりそうです。俗にいう「尤もらしさ」とはこのことらしい。ただし、既知にしろ未知にしろ、「ある分布に関する母数」を設定してあるということはパラメトリックな分布を想定していることになります。真の分布ではなく、統計モデルに落とし込んでいるということは意識として持っておく必要がありますね。

実際には、（生物系では）次々にデータ（標本）が増えるシチュエーションはなかなかないので、標本は固定で $\theta$ を色々動かして母数を推定することになります（例えば最尤推定）。

十分統計量

データから母数 $\theta$ を推定する事をパラメトリック推測と呼ぶ。この時使う統計量が、母数に関して情報損失が無い場合、その統計量を十分統計量という。すなわち、十分統計量は母数 $\theta$ とは無関係の関数となる。

Fisherの因子分解定理が成立すれば十分統計量ってことらしい。

Fisherの因子分解定理
標本 $X_1, \ldots, X_n$ に対して同時密度関数 $f(x_1, \ldots, x_n)$ が存在し、 $\boldsymbol{x} = (x_1, \ldots, x_n)$ とする。この時

$$f (\boldsymbol{x} | \theta) = g(T(\boldsymbol{x}) | \theta) h(\boldsymbol{x})$$

と分解できれば $T$ は十分統計量である。ということらしい。

証明は非常に簡単で、 $f(\boldsymbol{x} | \theta)$ にこの分解定理を適用して確率密度を計算しなおしてみればよい。

$$f (\boldsymbol{x | \theta}) = f (\boldsymbol{x} | T=t)g(t|\theta) = h(\boldsymbol{x}; t)g(t|\theta)$$

ここでセミコロンは、変数かgivenな値か区別する記号（だと思われる）。これを全事象について積分して割ってやれば確率密度を計算できるので、全事象を $A_t = \{\boldsymbol{x} | T(\boldsymbol{x}) = t \}$ として

$$\int_{A_t} f(\boldsymbol{x}|\theta) \mathrm{d}\boldsymbol{x} = g(T(\boldsymbol{x}) | \theta) \int_{A_t} h(\boldsymbol{x}; t) \mathrm{d}\boldsymbol{x}$$

とすると、

$$f(\boldsymbol{x} | T=t(\boldsymbol{x})) = \frac{f (\boldsymbol{x | \theta}) }{\int_{A_t} f(\boldsymbol{x}|\theta) \mathrm{d}\boldsymbol{x} } = \frac{h(\boldsymbol{x}; t)g(t|\theta)}{g(t | \theta) \int_{A_t} h(\boldsymbol{x} ; t) \mathrm{d}\boldsymbol{x}} = \frac{h(\boldsymbol{x}; t)}{\int_{A_t} h(\boldsymbol{x} ; t) \mathrm{d}\boldsymbol{x}}$$

となって母数 $\theta$ と無関係になる。

言われてみれば当然なんだけど、十分統計量だけ調べれば母数の推測は不要という点で強力ですね。

2018-03-08

Julia tips #9: CVODE (Sundials.jl) の数値積分が終わらない時の回避方法

Julia

追記
どうやら、せっかく設定したhminは内部でフラグチェックが行われないようです。
solve関数のmaxitersオプションで対応するしかないかもしれません。また後程検証してみます。

以下は検証が十分ではありません。
ソルバ自体の動作に支障はありませんが、期待通りに動いているかは検証する必要があります。
そのことをご理解いただいたうえでご覧ください。

2018/03/08当時の環境
- Julia 0.6.2 - Sundials 1.2.0

事の発端

ODEモデルのパラメータ推定をしていたら、非線形なモデルだとどうしても途中でソルバがフリーズしてしまう...。

原因はソルバが無限に数値積分のステップを続けてしまうことでした。

数値計算に疎いのでどうしてそんなことになってしまうのかはちゃんとわかっていないのですが、stiffなODEを解こうとすると、数値積分時に細かいステップを刻まないといけなくなるのが原因です。

「DifferentialEquationsから簡単にSundialsが使える！」と喜んだのもつかの間。

DifferentialEquations.jlでもこの件について触れられていました。

Frequently Asked Questions · DifferentialEquations.jl

公式によると、現状アルゴリズムを変えたり、ソルバを変えたりするしかないようです。

しかし、stiffなODEが解けないような状況においてCVODE以外ではお話にならないのです（経験上）。

そもそもODEが解けないようなパラメータは捨ててしまえばよいので、途中で計算を止めてなんとか捨てたいわけです。

Sundials? CVODE?

Sundialsは常微分方程式のソルバ、CVODEを作っているところです（他に幾つも作ってる）。多用しないのでわかりませんが、良く名前を聞くのでシェアは大きいと思います、MATLABにもあります。

stiffな問題をCVODEに解かせたら他のソルバより圧倒的に高速で数値積分してくれます。
stiffでなくてもめっちゃ速いです。過去記事中のリンク先でJuliaのソルバで速度比較をした方がいます。

hotoke-x.hatenablog.com

ちなみにCVODEとよく比較されるソルバにLSODAっていうのがあります。自分の場合はCVODEの方が5倍くらい速かったです。問題にもよると思うので、CVODEがダメだったら試す価値はあるかもしれません。JuliaならPkg.add("LSODA")ってやればDiffEq.jlのインターフェースから使えるようになります。以下のブログ記事ではLSODAも良いぞ！って言ってますね。

Solving Stiff ODEs

CVODEが途中で止まる

最強のCVODEといえど限界はあります。積分のステップが小さすぎると計算が終わらなくなります。途中でエラー吐けよって思うんですが、デフォルト設定だと無限にステップを刻めることになっています。こいつをちょちょいといじれば良いのです。

以下のhminというオプションをいじる（CVODE公式ドキュメント, p37より）↓ f:id:hotoke-X:20180308031520p:plain

The default value is 0.0となっていますね。

これがいかんわけです。

Sundials.jlを書き換える（自己責任で）

DiffEq.jlのインターフェースには残念なことに、hminを設定するオプションがありません。しかし、Sundials.jl自体ではちゃんとラップしてくれているのです。じゃあ直接Sundials.jlを使えばいいんですが、exampleを見てるとわかるんですがjuliaの可読性が損なわれます。何より面倒です。

なので、DiffEq.jlのインターフェースから使えるようにSundials.jlを少しだけいじりました。

自分でやる場合は自己責任でお願いします。 ...といっても、ミスっても書き換えたところを戻せば元通りです（2敗）

DiffEq.jlのインターフェースとなるsolve.jlをSundialsのパッケージディレクトリから探します。 f:id:hotoke-X:20180308033027p:plain

hminを引数としてsolveに追加（後ろにkwargs...があるので多分これやらなくても動くけど）。

f:id:hotoke-X:20180308033609p:plain

以下のようにhminをCVODEの設定として反映させる（ハイライト部分）。

f:id:hotoke-X:20180308033746p:plain

これで完了。

SundialsはCのライブラリなんですが、Sundilas.jlの作者さんがきれいにラップしてくれているおかげでインターフェース部分を少し書き換えるだけで良いです。作者さんありがとうございます。

後は動作確認したらよいです。

パッケージを書き換えたので、初めにプリコンパイルが走ります。

（以下ではCVODEが固まるほどのstiffなODEがどんなものかわからないので、ソルバが動くかどうかだけ確かめています。）

using DifferentialEquations
using Sundials
using Plots

function parameterized_lorenz(du,u,p,t)
 du[1] = p[1]*(u[2]-u[1])
 du[2] = u[1]*(p[2]-u[3]) - u[2]
 du[3] = u[1]*u[2] - p[3]*u[3]
end

u0 = [1.0,0.0,0.0]
tspan = (0.0,100.0)
p = [10.0,28.0,8/3]
alg = CVODE_BDF()
prob = ODEProblem(parameterized_lorenz, u0, tspan, p)
sol = solve(prob, hmin=1.0e-8)

plot(sol,vars=(1,2,3))

以下のようにグラフがプロットされればとりあえず正常に動作しています。

f:id:hotoke-X:20180308043930p:plain

しかし、本当に設定が反映されたのだろうか...

（どなたか、CVODEをフリーズさせる程のめちゃめちゃ硬いODEを教えてください）

元バイオ系

元バイオウェット系がデータサイエンスやらを勉強していくブログ。基本自分用のまとめ。

Julia1.0でマクロを@docしまくった

【Julia】Julia-1.0.0のProjectを試す【仮想環境？】

Projectの作成

文献検索、管理、メモを1つの画面で完結させる

Readcubeの凄いところ

使い方

その他

JuliaでPRML：ベイズ曲線フィッティング

JuliaでPRML：多項式曲線フィッティング

【数理統計】ラオーブラックウェル化（Rao-Blackwellization）

ラオーブラックウェル化（Rao-Blackwellization）

参考書籍

【数理統計】ブラックウェルーラオ（Blackwell-Rao）の定理

証明

参考書籍

【数理統計】最尤推定量とバイアス

尤度

最尤推定

バイアス

参考書籍

【数理統計】尤度関数と十分統計量

Chapter 7 統計的推論より

統計量・推定量・推定値

尤度関数と十分統計量

尤度関数

十分統計量

Julia tips #9: CVODE (Sundials.jl) の数値積分が終わらない時の回避方法

事の発端

Sundials? CVODE?

CVODEが途中で止まる

Sundials.jlを書き換える（自己責任で）