Q&Apymc3でのモデル関数が条件分岐を含む場合の書き方を教えていただきたい
pymc3を用いて、データ解析を行っています。
モデル関数（下記参照）がifを含む条件分岐を含んでいます。
条件分岐を含むpymcでのモデル関数の書き方について教えていただきたい。
条件分岐を含むモデル関数
uの位置で、1次関数として立ち上がる関数です。
Python
1@np.vectorize 2def model_line(x, u, nor, bg): 3 4 ud = x - u 5 if ud <= 0: 6 f = bg 7 else: 8 f = nor*ud + bg 9 return f
この関数を使ってpymcを行いました。
python
1import pymc3 as pm 2import numpy as np 3import scipy.stats as stats 4import matplotlib.pyplot as plt 5import seaborn as sns 6plt.style.use('seaborn-darkgrid') 7np.set_printoptions(precision=2) 8%matplotlib inline
この関数の作図用コード。
python
1#条件分岐を含む関数の作図 2xx = np.arange(3, 6.1, 0.1) 3 4y_data = model_line(xx,4.0,1.0,1.0) 5 6plt.plot(xx,y_data) 7plt.show()
事後分布を求める
python
1# x, 2# u, gauss 3# nor guass 4# bg gauss 5# y_dataはとりあえず作図用に計算で作ったものを入れています。 6with pm.Model() as model: 7 u= pm.Normal('u', mu=4.0, sd=1.0) 8 nor = pm.Normal('nor', mu=1, sd=1) 9 bg = pm.Normal('bg', mu=1, sd=1) 10 epsilon = pm.HalfCauchy('epsilon', 1) 11 12 mu =model_line(xx,u,nor,bg) 13 14 y_pred = pm.Normal('y_pred', mu=mu, sd=epsilon, observed=y_data) 15 16 trace_model = pm.sample(2000, njobs=1) 17 18pm.traceplot(trace_model) 19plt.tight_layout() 20 21plt.figure()
発生している問題・エラーメッセージ
--------------------------------------------------------------------------- TypeError Traceback (most recent call last) <ipython-input-5-a9cdfea62e7a> in <module> 9 epsilon = pm.HalfCauchy('epsilon', 1) 10 ---> 11 mu =model_line(xx,u,nor,bg) 12 13 y_pred = pm.Normal('y_pred', mu=mu, sd=epsilon, observed=y_data) ~\Anaconda3\envs\pymc3\lib\site-packages\numpy\lib\function_base.py in __call__(self, *args, **kwargs) 1970 vargs.extend([kwargs[_n] for _n in names]) 1971 -> 1972 return self._vectorize_call(func=func, args=vargs) 1973 1974 def _get_ufunc_and_otypes(self, func, args): ~\Anaconda3\envs\pymc3\lib\site-packages\numpy\lib\function_base.py in _vectorize_call(self, func, args) 2040 res = func() 2041 else: -> 2042 ufunc, otypes = self._get_ufunc_and_otypes(func=func, args=args) 2043 2044 # Convert args to object arrays first ~\Anaconda3\envs\pymc3\lib\site-packages\numpy\lib\function_base.py in _get_ufunc_and_otypes(self, func, args) 2000 2001 inputs = [arg.flat[0] for arg in args] -> 2002 outputs = func(*inputs) 2003 2004 # Performance note: profiling indicates that -- for simple <ipython-input-2-22edf684422a> in model_line(x, u, nor, bg) 17 """ 18 ud = x - u ---> 19 if ud <= 0: 20 f = bg 21 ~\Anaconda3\envs\pymc3\lib\site-packages\theano\tensor\var.py in __bool__(self) 89 else: 90 raise TypeError( ---> 91 "Variables do not support boolean operations." 92 ) 93 TypeError: Variables do not support boolean operations.
試したこと
PyMC3で簡単なMCMCチュートリアルを試したメモ
には、以下のようなコードが載っていますが、理解できていません。
（1）itypesは、私の場合、4つありますが、一番最初のxはベクトルにすべきかどうか。Outputもベクトルにすべきか。
（2）vectorizeは関数のデコレータで行っておいた方が良いか
（3）下記のコードを参考に書き換えてもうまくいきませんでした。
python
1# 参照したHPに記載の情報 2import theano.tensor as T 3from theano.compile.ops import as_op 4 5@as_op(itypes=[T.lscalar], otypes=[T.lscalar]) 6def crazy_modulo3(value): 7 if value > 0: 8 return value % 3 9 else : 10 return (-value + 1) % 3
これを踏まえ先ほどのモデル関数を以下のように書き換えました。
python
1import theano.tensor as T 2from theano.compile.ops import as_op 3 4@as_op(itypes=[T.dvector,T.dscalar,T.dscalar,T.dscalar], otypes=[T.dvector]) 5@np.vectorize 6def model_line(x, u, nor, bg): 7 ud = x - u 8 if ud <= 0: 9 f = bg 10 else: 11 f = nor*ud + bg 12 return f
再度事後分布を実行
python
1# x, 2# u, gauss 3# nor guass 4# bg gauss 5with pm.Model() as model: 6 u= pm.Normal('u', mu=4.0, sd=1.0) 7 nor = pm.Normal('nor', mu=1, sd=1) 8 bg = pm.Normal('bg', mu=1, sd=1) 9 epsilon = pm.HalfCauchy('epsilon', 1) 10 mu =model_line(xx,u,nor,bg) 11 y_pred = pm.Normal('y_pred', mu=mu, sd=epsilon, observed=y_data) 12 trace_model = pm.sample(2000, njobs=1) 13 14pm.traceplot(trace_model) 15plt.tight_layout() 16plt.figure()
再度エラーがでました
python
1--------------------------------------------------------------------------- 2AttributeError Traceback (most recent call last) 3<ipython-input-8-a9cdfea62e7a> in <module> 4 9 epsilon = pm.HalfCauchy('epsilon', 1) 5 10 6---> 11 mu =model_line(xx,u,nor,bg) 7 12 8 13 y_pred = pm.Normal('y_pred', mu=mu, sd=epsilon, observed=y_data) 9 10~\Anaconda3\envs\pymc3\lib\site-packages\theano\gof\op.py in __call__(self, *inputs, **kwargs) 11 613 """ 12 614 return_list = kwargs.pop('return_list', False) 13--> 615 node = self.make_node(*inputs, **kwargs) 14 616 15 617 if config.compute_test_value != 'off': 16 17~\Anaconda3\envs\pymc3\lib\site-packages\theano\gof\op.py in make_node(self, *inputs) 18 981 raise ValueError("We expected %d inputs but got %d." % 19 982 (len(self.itypes), len(inputs))) 20--> 983 if not all(inp.type == it for inp, it in zip(inputs, self.itypes)): 21 984 raise TypeError( 22 985 "We expected inputs of types '%s' but got types '%s' " % 23 24~\Anaconda3\envs\pymc3\lib\site-packages\theano\gof\op.py in <genexpr>(.0) 25 981 raise ValueError("We expected %d inputs but got %d." % 26 982 (len(self.itypes), len(inputs))) 27--> 983 if not all(inp.type == it for inp, it in zip(inputs, self.itypes)): 28 984 raise TypeError( 29 985 "We expected inputs of types '%s' but got types '%s' " % 30 31AttributeError: 'numpy.ndarray' object has no attribute 'type' 32
補足情報
Pythonは、3.5を利用しています。
Ancondaを利用してpymc3はインストールしました。
回答1件
0
自己解決
質問をした後、いろいろと調べたところ以下のように書き換えたところうまく動作いたしました。
ご検討いただいた皆様ありがとうございました！
model_line関数とは別の関数を作ります。
python
1def model_line3(x, u, nor, bg): 2 3 ud = x - u 4 f1=bg 5 f2=nor*ud + bg 6 7 return pm.math.switch(ud<=0,f1,f2)
pymc3のswich関数を用いました。引数の1番目は条件、2番目は真の条件での実行、3番目は、偽の条件での実行を入力します。
そして、再度mu関数をmodel_line3に書き換え実行します。
python
1# x, 2# u, gauss 3# nor guass 4# bg gauss 5with pm.Model() as model: 6 u= pm.Normal('u', mu=4.0, sd=1.0) 7 nor = pm.Normal('nor', mu=1, sd=1) 8 bg = pm.Normal('bg', mu=1, sd=1) 9 epsilon = pm.HalfCauchy('epsilon', 1) 10 11 mu =model_line3(xx,u,nor,bg) 12 13 y_pred = pm.Normal('y_pred', mu=mu, sd=epsilon, observed=yy_data) 14 15 trace_model = pm.sample(2000, njobs=1) 16 17pm.traceplot(trace_model) 18plt.tight_layout() 19 20plt.figure()
参考にしたサイトPyMC3を使って変化点検出のベイズ推論をする
投稿2019/03/25 00:05
総合スコア7
あなたの回答
tips
プレビュー
バッドをするには、ログインかつ
こちらの条件を満たす必要があります。