Weblog on mebius.tokaichiba.jp

吃音治療に関する興味深い論文を入手した。
"Spontaneous" late recovery from stuttering: Dimensions of reported techniques and causal attributions
成人してから吃音が「自然に」治癒した事例について、用いたテクニックと、関係したかも知れない原因やきっかけついて統計解析をしたものである。
丁度、筆者は吃音者であり、しかも今、英語と統計学の復習をしているので、読んでみるしかないと思った。

統計関連は後にして、ここではこの論文の要約を記す。
章節毎に要約する。

■ABSTRACT
目的：
(1) 青年期以降に治療なしで吃音から回復した成人が用いたテクニックと、原因やきっかけを調査する。
(2) それらのテクニックや原因・きっかけが少数の次元にまとめられるかを調べる。
(3) それらの次元が一般的な治療方法の要素と対応しているかを調べる。

方法：
吃音から回復した124人に、49の用いられそうなテクニックと、15のあり得そうな原因・きっかけについて、当てはまるかどうかを回答してもらった。

結果：
110件のアンケート結果の主成分分析により、テクニックについては6つの成分、きっかけについては3つの成分が抽出された。

議論：
テクニックの2成分（話し方の工夫(Speech Restructuring)、雄弁術）は治療方法に対応している。別の成分（リラックス／話し方観察）は、一般的に有効と考えられている対処方法に対応している。
当事者が感じる様々な原因・きっかけからの成分は、吃音が治る要因の暗黙的な理論の差を反映している。
これらの成分の特定は、単にテクニックやきっかけをリストにしただけの従来の研究と比べて新しいものであり、今後の吃音治療の道標になるかも知れない。

■1. Introduction
吃音は幼少期には5-11%の人に見られるが、成人では0.5%（女性0.2%、男性0.8%）しか見られない。幼少期の吃音の70-80%は成長の過程で無くなる。そのほとんどは治療もせず自然に無くなるものである。
青年期に持ち越した後も自然に治ることがある。この研究は青年期以降に「自然に」治った人を対象にする。
ここでは、治ったとは吃音になったことがない人と同じくらいに完全に無くなったに限らず、感情や状況によっては出るが社会生活において支障が無い程度まで無くなったのを含む。それが実際の吃音治療の現実的なゴールとも合っている。

■1.1. Validity of self-reports on recovery
(1) その人が過去に吃音だったことと治療とは無関係に治ったことをどう確認するか？
→参考人（証人）に確認する方法が開発されている
(2) 吃音から回復した人の話し方と吃音歴の無い人の話し方は同じなのか？
→多くの場合、非流暢性が残る
(3) 吃音から回復した人は喋ることを求められるなど緊張のある状況でも出ないのか？
→回復したと主張する人の60%には吃音の傾向が残っている

■1.2. How spontaneous is a "spontaneous" late recovery?
自然に治ったという人もよく何らかの対処法を口にする。それは自分で見つけた方法とは限らない。過去に受けた治療を参考にしているかも知れない。

■1.3. What can be learned from "spontaneous" late recovery?
「自然に」治った例は治療法に示唆を与えると言われてきた。
従来の研究は定性的、今回の研究は定量的。

■1.4. Aim of the study
(1) テクニックと原因・きっかけがどれくらいあったか
(2) それらはより少ない次元にまとめられるか

■2. Method
■2.1. Participants
11歳以降に吃音が治ったという124人にアンケートした。
可能な人には参考人にもアンケートしたが、参考人が得られたのは半数以下だった。
吃音以外の言語障害だったと思われる人や、規定する程度に治ったと言えない人14人は除外した。
残った110人の年齢は14-75歳、平均44歳、吃音が始まった年齢の平均は4.8歳。
吃音治療歴は42%が治療経験なし、24%は2年以下、34%は2年以上。

■2.2. Questionnaires
吃音が治ったという人へのアンケートの質問は、その人の人物特徴、吃音が始まった年齢と終わった年齢、現在の吃音の状態、49のテクニックの内使ったもの、15の原因やきっかけの内関係したと思うもの。
これら64項目は先行研究から拾ったり、一部の著者を含む数名の過去に吃音だった人や言語聴覚士から聞き取ったもの。
各項目の選択肢は
(1) "definitely not true", (2) "probably not true", (3) "don't know", (4) "probably true", (5) "definitely true"の5値だった。

■2.3. Statistical analysis
用いられたテクニックと原因・きっかけは主成分分析を用いて次元削減した。
標本サイズが十分かどうかはKaiser-Meyer-Olkin(KMO) indexが0.5以上となることで判定した。テクニックについては0.75、原因・きっかけについては0.62だった。
アンケートの値は順序尺度の5値なので、不適切な正規分布の仮定を避ける為、polychoric correlation matrixによる主成分分析を用いた。
最適な主成分の数は、平行分析によって求めた。

■3. Results
■3.1. Validity of corroborator reports
参考人の回答は31人から得られ、その内回答が揃っていたのは20人だった。
過去の状態については、吃音特有の項目の平均が2.92(N=26)、吃音に無関係の項目の平均が1.15(N=25)であり、ｔ検定ではt=8.47、p-value<.001で有意で、吃音で間違いなさそうと確認された。
現在の状態については、吃音特有の項目の平均が1.49、その他の項目の平均が1.31であり、p-value=.092で、有意な差が無かった。
参考人の回答と吃音当事者本人の回答との相関係数は、吃音関連の項目についてはr=.46、p=.017と高くて有意であり、その他の項目についてはr=.21、p=0.32と低く、有意でなかった。

■3.2. Techniques used
49のテクニックをカイザー基準化、バリマックス回転して主成分分析し、平行分析により6つが最適とされたので6つの主成分を抽出した。
それぞれの主成分について、主成分負荷量が0.5より大きいテクニックを抜き出し、主成分負荷量×テクニックの使用率を重要度とし、重要度順に並べた（表1）。

第1主成分Aは全分散の16.2%を占めた。主な項目は「単語の最初の音を伸ばした」「音節と音節の区切りを曖昧にして繋げた」「単語の終わりと次の単語の始まりを繋げた」「単語のアクセントを変えた」「深い声で話した」「リズムよく話した」なので、この成分を「話し方の工夫」とラベル付けした。

第2主成分Bは全分散の15.1%を占めた。主な項目は「できるだけリラックスに努めた」「話の区切りを十分に取るよう気を付けた」「何をどのように話すかを考える時間をより長くした」「常に自分の話し方、言い方に集中するようにした」「話す間のリラックスした呼吸に集中するようにした」「話すテンポを下げた」なので、この成分を「リラックス／話し方観察」とラベル付けした。

第3主成分Cは全分散の10.9%を占めた。主な項目は「大きな声で読むか、誰かの前で読む」「問題のあった音や単語を分析した」「詩を読んだ」「非吃音者と一緒に大きな声で朗読した」「鏡の前で自分自身に話した」であり、この成分を「雄弁術」とラベル付けした。

第4主成分Dは全分散の7.5%を占めた。主な項目は「歌のレッスンを受けた」「演劇のレッスンを受けた」であり、この成分を「演劇」とラベル付けした。

第5主成分Eは全分散の7.1%を占めた。主な項目は「公衆の前でスピーチをした」「議論に積極的に参加した」「敢えて困難な話す場に参加した」であり、この成分を「話す必要性の探求」とラベル付けした。

第6主成分Fは全分散の6.9%を占めた。主な項目は「話し始めを楽にする為に直前にアイコンタクトした」「特に吃音が出た時や予期した時にアイコンタクトを維持するよう努めた」であり、「安心感の確保」とラベル付けした。

この6主成分で全分散の63.7%が説明される。他の8主成分は固有値>1になり得たが、平行分析により無視された。

表1. 非流暢さの低減の為に使ったテクニック、使った人の割合、主成分負荷量、重要度（割合×負荷量）

Technique used（用いたテクニック）	%	Loading	Importance
Component A: Speech Restructuring（話し方の工夫）
Stretched out the beginnings of words 単語の最初の音を伸ばした	22	.72	15.8
Slurred speech gliding smoothly from syllable to syllable 音節と音節の区切りを曖昧にして繋げた	21	.68	14.3
Connected the end of a word with the beginning of the next one 単語の終わりと次の単語の始まりを繋げた	20	.68	14.3
Changed word accent 単語のアクセントを変えた	18	.71	12.8
Spoke with a deeper voice 深い声で話した	13	.78	10.1
Spoke rhythmically without an external pacemaker メトロノームとか無しでリズムよく話した	13	.56	7.3
Spoke with a higher voice 高い声で話した	7	.80	5.6
Controlled lip muscles such that lips did not touch each other 唇がくっつかないよう唇の筋肉を動かした	10	.51	5.1
Sung utterances (rap, parlando) 歌うように話した（ラップのように）	8	.56	4.5
Reduced perception of own voice by speaking with masking noise マスクノイズを使って自分の声が聞こえにくくした	6	.74	4.4
Took relaxing or other anti-stuttering medication 緊張を緩和する薬を使った	5	.76	3.8
Spoke with an external pacemaker (e.g. metronome) メトロノームとかを使ってリズムよく話した	5	.62	3.1
Component B: Relaxed/Monitored Speech（リラックス／話し方観察）
Tried to be more relaxed at all できるだけリラックスに努めた	40	.72	28.8
Took care to make enough speech breaks 話の区切りを十分に取るよう気を付けた	38	.74	28.1
Took more time for thinking beforehand what to say and how to say it 何をどのように話すかを考える時間をより長くした	33	.79	26.1
Always tried to concentrate on my speaking manner 常に自分の話し方、言い方に集中するようにした	33	.76	25.1
Concentrated on even and relaxed breathing while speaking 話す間のリラックスした呼吸に集中するようにした	37	.67	24.8
Reduced speech tempo 話すテンポを下げた	33	.65	21.5
Concentrated on what I will say, not how I will say it どう話すかでなく何を話すかに集中した	29	.70	20.3
First thought, then spoke 考えてから話すようにした	33	.56	18.5
Took a deep breath before each utterance 都度発声前に深呼吸した	22	.64	14.1
Tried to speak clearly with precise articulation 正確な発音で明瞭に話すよう努めた	22	.63	13.9
Relaxed before certain utterances and spoke as soon as I was relaxed 吃る音の発声前にリラックスし、リラックスできたら直ちに発声した	20	.58	11.6
Component C: Elocution（雄弁術）
Read aloud or to someone 大きな声で読むか、誰かの前で読む	29	.63	18.3
Analyzed sounds and words with which I had problems 問題のあった音や単語を分析した	25	.64	16.0
Analyzed what my articulatory muscles did during stuttering 吃った時の発声器官の筋肉の動きを分析した	18	.68	12.2
Read poems 詩を読んだ	15	.76	11.4
Worked on my secondary symptoms (e.g. tics, arm movements, grimaces) 二次的な症状に対処した（ゆすり、腕の動き、しかめっ面など）	14	.59	8.3
Read texts aloud together with a fluent speaker 非吃音者と一緒に大きな声で朗読した	12	.66	7.9
Talked to myself in front of a mirror 鏡の前で自分自身に話した	10	.63	6.3
Purposely used difficult words わざと言いにくい単語を選んだ	12	.52	6.2
Component D: Stage Performance（演劇）
Took singing lessons 歌のレッスンを受けた	10	.74	7.4
Took acting lessons 演劇のレッスンを受けた	7	.74	5.2
Learned to play a wind instrument 吹奏楽器の演奏を学んだ	2	.88	1.8
Component E: Sought Speech Demands（話す必要性の探求）
Held public speeches 公衆の前でスピーチをした	25	.89	22.3
Participated actively in discussions 議論に積極的に参加した	34	.63	21.4
Participated purposely in difficult speech situations 敢えて困難な話す場に参加した	20	.79	15.8
Component F: Reassurance（安心感の確保）
Made eye contact seconds before speaking in order to ease speech onset 話し始めを楽にする為に直前にアイコンタクトした	24	.71	17.0
Tried to keep eye contact, especially when I stuttered or expected to 特に吃音が出た時や予期した時にアイコンタクトを維持するよう努めた	26	.59	15.3
Paid heed to a stress-free life with sufficient sleep and physical activity 十分な休息と運動でストレスの少ない生活を心掛けた	23	.65	15.0
Items not loading sufficiently on any main component（どの主成分にも十分な関係が無いもの）
Stopped when I stuttered a word and repeated it 吃ったら一度止まって繰り返した	30
Sought out relaxed speech situations リラックスして話せる場を探した	28
Spoke with soft voice onset 出だしは優しく発声した	24
Quit avoiding or replacing difficult words 言いにくい言葉を避けたり言い換えたりするのをやめる	23
Increased voice volume (e.g. through deep breathing) 呼吸を深くして声を大きくした	19
Consciously relaxed parts of my speech apparatus 発声器官の部位を意識的にリラックスさせた	16
Used relaxation exercises リラックスする為の訓練法を使った	15
Stuttered purposely but differently than usual わざと普段と異なる吃り方をした	11
Participated in rhetoric courses 雄弁術の講座に参加した	9

■3.3. Causal attributions for the recovery
15の原因やきっかけについて、同様に主成分分析を行い、3つの主成分を抽出した。（表2）

第1主成分(I)は全分散の24.1%を占めた。主な項目は「新しい職場で仕事を始めた」「学校を卒業した」「住居を変えた」「新しいパートナーに巡り合った」「違う言語の国に移住した」なので、この成分を「生活の変化」とラベル付けした。

第2主成分(II)は全分散の17.9%を占めた。主な項目は「自信が増した」「自分の話し方を否定的に捉えなくなった」「自分の吃音に対する他人の反応が気にならなくなった」「自分が吃音者だと思わなくなった」なので、この成分を「姿勢の変化」とラベル付けした。

第3主成分(III)は全分散の16.0%を占めた。主な項目は「大事な人に動機づけられた」「自分の吃音について何か変えようと決心した」「大事な人に吃音に取り組んでいることを伝えた」「心理学者や心理療法士の指導を受けた」なので、この成分を「社会的支援」とラベル付けした。

表2. 非流暢さの低減に繋がった原因やきっかけ、該当する人の割合、主成分負荷量、重要度（割合×負荷量）

Reason or occasion to which disfluency reduction was attributed（非流暢さの低減に繋がった原因やきっかけ）	%	Loading	Importance
Component I: Life Change（生活の変化）
I started at a new workplace 新しい職場で仕事を始めた	18	.82	14.8
I finished school 学校を卒業した	18	.82	14.8
I met a new partner 新しいパートナーに巡り合った	12	.80	9.6
I changed my residence 住居を変えた	11	.90	9.9
I moved into a country with another language 違う言語の国に移住した	5	.78	3.9
Component II: Attribute Change（姿勢の変化）
My self-confidence increased 自信が増した	56	.63	35.3
I did not evaluate my speech negatively any longer 自分の話し方を否定的に捉えなくなった	28	.82	23.0
I became insensitive about reactions by others to my stuttering 自分の吃音に対する他人の反応が気にならなくなった	22	.75	16.5
I no longer accepted that I was one who stuttered 自分が吃音者だと思わなくなった	22	.67	14.7
Component III: Social Support（社会的支援）
I decided to change something about my stuttering 自分の吃音について何か変えようと決心した	35	.53	18.6
I was motivated by persons who were important to me 大事な人に動機づけられた	29	.83	24.1
I told a person important to me that I was working on my stuttering 大事な人に吃音に取り組んでいることを伝えた	15	.85	12.8
I had counseling by a psychologist/psychotherapist 心理学者や心理療法士の指導を受けた	8	.68	5.4
Items not loading sufficiently on any main component（どの主成分にも十分な関係が無いもの）
I decided to treat my stuttering my own way 自己流で吃音を直すと決めた	40
I accepted the fact that I was one who stuttered 自分が吃音者であることを受け入れた	33

■3.4. Age and gender effects
年齢の違い（中央値より大か小か）、性別の違いについて、Mann-WhitneyのU検定で確認した結果、どの項目についても有意な差は無かった。
ただ1項目、年齢が高い層には、「自分の吃音に対する他人の反応が気にならなくなった」が多かった。

■4. Discussion
■4.1. Comparison with previous findings
従来の研究と比較する。
話す速度を下げること、話し方を変えること、は従来の研究同様、よく使われていた。
リラックスは従来の研究とは異なり、よく使われていた。
雄弁術や話すことを求められる場の追求は、今回の調査では比較的多かったが、従来の研究では特筆されなかった。

原因やきっかけについては、従来の研究と同様、自信の向上が一番であり、話し方を変えたい欲求も多かった。生活環境の変化は、割合としてはさほど多くは無いものの、従来の研究では特筆されなかったものである。

加えて、今回の研究では、アルゴリズム的に項目をグループ化したので、類似の、同時に使われやすいテクニックや、暗黙的な回復の要因を示す類似の項目が得られた。

■4.2. Do the components mirror stuttering treatments?
今回検出した、成人後に「自然に」治った事例の主成分は、成功している吃音治療と類似があるか、治療法に示唆があるか、という問いに対しては、2つの理由で、あると考える。
(1) 既存の治療方法は、理論や実証だけではなく、吃音者の肯定的な経験にも基づくから。
(2) 成人してからの回復は本当に「自然」に回復したとは限らない。いくらかは、過去の治療や、その他の情報源から得た既存の治療法を参考にしたと考えるべきである。

第1主成分の「話し方の工夫」は、単語の始まりを伸ばすなど、吃音治療で使われる方法で主に構成されている。
第3主成分の「雄弁術」は、20世紀半ばまで一般的な治療法だった。
第2主成分の「リラックス／話し方観察」は正式な治療法を反映しているようには見えないが、自身で行う方法としては一般的であり、治った本人の感覚としては効果があった要因として挙げられることが最も多いという報告もある。

原因やきっかけについて検出した3つの主成分に関しては、いずれも治療法に繋がり得るものであるし、自分や自分の話し方の「姿勢の変化」は吃音からの回復の主要因の1つとして報告されたことがある。

■4.3. Attribution biases
テクニックは実際に用いたからそれほど当てはまらないが、主観的な原因やきっかけはいくつかのバイアスがあり、額面通りに受け入れられないので注意が必要。
(1) 吃音が治ったというような通常でない作用は、同時期に起こったイベントを原因にされやすい。
(2) 原因やきっかけは事後に言っているものであり、よくわからないものを説明できるように、認知的不協和を減らせるように考えられやすい。
(3) 原因やきっかけは自分をひいきしやすい。成功の要因は自分の中に求め、失敗の要因は外部に求めやすい。

■4.4. Validation of former stuttering of the participants
吃音から回復したという自己申告はどこまで信用できるかについては、参考人の回答から、過去の吃音に関する回答の平均値と現在の吃音に関する回答の平均値の差の効果量が3σに相当するくらい大きいことから明確に説明される。

■4.5. Strengths and limitations of the study
この研究は主成分分析という主観的でない手法を使っていることに意味があり、従来の研究を否定するものではなく補強するものである。
1つの限界として、吃音が治ったの程度の問題がある。例えば「現在吃音の傾向がある」に「滅多にない」と回答した62人中26人は「緊張やプレッシャーがある時には吃る」に「時々ある」と回答した。従って、過去に吃音があった人は、状況によっては吃る程度なら治ったと考えている。
もう1つの限界として、「自然に」治ったの意味がある。人は情報や事実を覚えていても情報源を覚えていない傾向があることが知られている。回答者の52%は過去6ヶ月以上前に吃音治療を受けたことがあると回答している。過去の治療の影響は除外できない。
もう1の限界として、言語療法士による参加者の吃音の診断結果が無いことがある。ただ、(1)58%の人は言語療法士による治療を受けたこと、(2)信頼できる手法を使って過去の吃音を調べたこと、(3)参加者は吃音の専門家や自助組織を経由して集められたこと、からして参加者が過去に吃音が無かった人だとは考えられない。

■4.6. Conclusion and future directions
この研究で報告した吃音から回復した人が使ったテクニックや原因・きっかけは、吃音が、感覚運動系、社会性、動機、神経系に関連し、症状が状況、感情、注意、言語学的な要因に依存する、複雑な障害であるという現在の理解と合う。
今回統計学的な手法で抽出した次元は、現在の吃音治療の改良にも役立ち得る。(1)「話し方の工夫」「雄弁術」は現在の、そして驚くべきことに過去の吃音治療の構成要素に対応する。(2)「リラックス／話し方観察」は専門的ではないが一般に有効と考えられている吃音への対処法であり、吃音治療に取り入れられる可能性がある。(3)原因やきっかけの要素は、素人理論を反映しており、吃音治療で利用できる可能性がある。
今回の調査結果は吃音の脳画像解析にも示唆を与える可能性がある。最近の研究で、吃音が「自然に」治ったという人と流暢性形成法の治療を受けた吃音者の治療前後のfMRI画像を比較したものがあり、「自然に」治ったという人は脳神経回路における発話時の異常に対応する部位の挙動が独立だったとしている。これが「自然に」治ったことを示すものなら、今回の調査結果に注目した吃音治療によって同様の脳画像が得られることが示せるかも知れない。

前記事の続きである。
苦労の末に前記事のCompileErrorの問題をクリアして、「空間自己回帰モデル」に進んだら、そのすぐ後の状態空間モデルとCARモデルの比較でまたPyMCのトラブルに遭遇した。

まず、トラブルが起こらなかった方の、p.51〜p.52のIntrinsic CAR(Conditional Autoregressive Prior)モデルのサンプルコードをPyMCで書き直したものを貼る。

●データ

import pandas as pd
import matplotlib.pyplot as plt

# サンプルデータ
# https://kuboweb.github.io/-kubo/stat/iwanamibook/fig/spatial/Y.RData の中のYのみ
Y = [ 0,  3,  2,  5,  6, 16,  8, 14, 11, 10, 17, 19, 14, 19, 19, 18, 15,
       13, 13,  9, 11, 15, 18, 12, 11, 17, 14, 16, 15,  9,  6, 15, 10, 11,
       14,  7, 14, 14, 13, 17,  8,  7, 10,  4,  5,  5,  7,  4,  3,  1]
N = len(Y)

plt.scatter(range(1, N+1), Y)
plt.ylim(-1, 25)
plt.xlabel('position')
plt.ylabel('count')
plt.show()

これは隣接する50地点の植物の個体数をカウントしたという想定で作られたデータとのことである。この観測値はポアソン分布に従うとし、各地点の平均にIntrinsic CARモデルを当てはめてみる。 幸いPyMCにもpymc.ICARというIntrinsic CARのライブラリがあったので、これを使った。

●コード

import numpy as np
import pymc as pm

# adjacency matrix（隣接行列）
adj = np.array([[j == i+1 or j == i-1 for i in range(N)] for j in range(N)]).astype(int)

with pm.Model() as model:
    # beta = pm.Normal('beta', 0, 1e-4)  # 本の通りだがこれだと全く平滑化されない、それっぽい結果の時のbetaは2より大きい
    beta = pm.Normal('beta', 0, 100)
    sigma = pm.Uniform('sigma', 0, 100)

    # centered parameterization
    # これだと全く平滑化されない
    # S = pm.ICAR('S', W=adj, sigma=sigma)
    # l = pm.Deterministic('lambda', (beta + S).exp())

    # non-centered parameterization
    S = pm.ICAR('S', W=adj)
    l = pm.Deterministic('lambda', (beta + sigma * S).exp())

    Yobs = pm.Poisson('Yobs', mu=l, observed=Y)
    
with model:
    trace = pm.sample(1000, progressbar=False)

●実行結果

# 結果の可視化
import arviz as az

x = range(1, N+1)
plt.scatter(x, Y, label='observed')
plt.ylim(-1, 25)

lambda_mean_CAR = pm.summary(trace, var_names='lambda')['mean']
plt.plot(x, lambda_mean_CAR, color='r', linewidth=2, label='posterior lambda')

# 95%ベイズ信用区間
hdi_CAR = az.hdi(trace, hdi_prob=0.95, var_names='lambda')
plt.fill_between(x, hdi_CAR['lambda'][:,0], hdi_CAR['lambda'][:,1], color='pink', alpha=0.15, label='95% credible interval')

plt.xlabel('position')
plt.ylabel('count')
plt.legend()
plt.show()

理屈がわかってない為に、平滑化されるまでに苦労したが、本と大体同じ結果が得られた。

次に、問題のp.55辺りの状態空間モデルの例である。1次元1階差分のIntrinsic CARモデルはローカルレベルモデルと等価とのことで、大体同じ結果になることを確認する。

●コード

%%time

# 発見個体数からのパラメーター推定
import pymc as pm
import pandas as pd

# 不具合回避用
# N = 31  # N > 31だとAttributeError: 'Scratchpad' object has no attribute 'ufunc'

with pm.Model() as model:
    # 事前分布
    alpha = [pm.Uniform('alpha[1]', 0, 10)]
    sigma = pm.Uniform('sigma', 0, 100)
    
    # システムモデル
    for i in range(1, N):
        alpha.append(pm.Normal(f'alpha[{i+1}]', mu=alpha[i-1], sigma=sigma))

    # 観測モデル
    for i in range(N):
        lambda_tmp = pm.Deterministic(f'lambda[{i+1}]', alpha[i].exp())
        Yobs = pm.Poisson(f'Y[{i+1}]', mu=lambda_tmp, observed=Y[i])

with model:
    trace = pm.sample(1000, progressbar=False)

●実行結果

File ~/Library/Python/3.10/lib/python/site-packages/pytensor/tensor/elemwise.py:747, in Elemwise.perform(self, node, inputs, output_storage)
    745         ufunc = self.ufunc
    746     else:
--> 747         ufunc = node.tag.ufunc
    748 else:
    749     ufunc = node.tag.ufunc

File ~/Library/Python/3.10/lib/python/site-packages/pytensor/graph/utils.py:286, in Scratchpad.__getattribute__(self, name)
    285 def __getattribute__(self, name):
--> 286     return super().__getattribute__(name)

AttributeError: 'Scratchpad' object has no attribute 'ufunc'

HINT: Re-running with most PyTensor optimizations disabled could provide a back-trace showing when this node was created. This can be done by setting the PyTensor flag 'optimizer=fast_compile'. If that does not work, PyTensor optimizations can be disabled with 'optimizer=None'.
HINT: Use the PyTensor flag `exception_verbosity=high` for a debug print-out and storage map footprint of this Apply node.

これまた、前回と同じく、N<=31ではエラーが発生しなかった。しかし、前回のようにg++のCompileErrorではないので、当然同じ対策は無効だった。
今度はエラーメッセージにも回避策に繋がるヒントがほとんど無い。Webで検索しても、同じエラーが出たという情報はあるが解決方法が全く見つからなかった。

最初に発生を確認したバージョンは Python 3.10.14 + PyMC 5.15.0 だが、PyMCのドキュメントに従ってMiniforgeのcondaを用いて Python 3.12.3 + PyMC 5.15.1 をインストールしても結果は同じだった。

発生箇所は、PyTensorのelemwise.pyのperform()のこの部分の最後の行である。

            if self.ufunc:
                ufunc = self.ufunc
            elif not hasattr(node.tag, "ufunc"):
                # It happen that make_thunk isn't called, like in
                # get_underlying_scalar_constant_value
                self.prepare_node(node, None, None, "py")
                # prepare_node will add ufunc to self or the tag
                # depending if we can reuse it or not. So we need to
                # test both again.
                if self.ufunc:
                    ufunc = self.ufunc
                else:
                    ufunc = node.tag.ufunc  # <--

この部分で、self.ufuncもnode.tag.ufuncも無ければself.prepare_node()を呼び出してself.ufuncかnode.tag.ufuncかを得ることになっているが、prepare_node()のコードでは、ある条件ではどちらにもufuncが入らないようになっている。

    def prepare_node(self, node, storage_map, compute_map, impl):
        # Postpone the ufunc building to the last minutes due to:
        # - NumPy ufunc support only up to 32 operands (inputs and outputs)
        #   But our c code support more.
        # - nfunc is reused for scipy and scipy is optional
        if (len(node.inputs) + len(node.outputs)) > 32 and impl == "py":
            impl = "c"

NumPyの制限により、引数の数が32を超える場合はNumPyでなく"c code"を使うと書かれているが、impl = "c"の場合、この後selfにもnode.tagにもufuncが入らない。
perform()の最初の方のコメントにも

    def perform(self, node, inputs, output_storage):
        if (len(node.inputs) + len(node.outputs)) > 32:
            # Some versions of NumPy will segfault, other will raise a
            # ValueError, if the number of operands in an ufunc is more than 32.
            # In that case, the C version should be used, or Elemwise fusion
            # should be disabled.
            # FIXME: This no longer calls the C implementation!
            super().perform(node, inputs, output_storage)

と書かれており、その場合にはNumPyでなくCのufuncを使う前提の設計になっているようだが、どこを探してもCのufuncの実装が見つからなかった。FIXME:の所に書かれているように、これは直さないといけないことが認識されているが、直されていないようだ。
最新のPyMCのコードを確認しても、この辺りのコードに変化は無かった。

super().perform() は OpenMPOp.perform() のようなので、OpenMPを有効にすれば変わるかも知れないと思ったが、 ~/.pytensorrc で openmp = True にしてもWarningが出てFalseになった。
"Elemwise fusion"を停止するには、この辺りを見ると ~/.pytensorrc で optimizer = None とか optimizer_excluding = elemwise_fusion とかとすると良さそうだが、結果は同じだった。

前回と同様、システムモデルの部分にpm.GaussianRandomWalkを使えば何の問題も起こらなかったが、それでは本のサンプルコードと形が違い過ぎるし、状態空間モデルをPyMCで直接的に書きたくてこの問題に遭遇した場合に解決しないので、満足できなかった。

●コード

# 発見個体数からのパラメーター推定（pm.GaussianRandomWalk使用）
import pymc as pm
import pandas as pd

N = 50

with pm.Model() as model:
    # 事前分布
    sigma = pm.Uniform('sigma', 0, 100)

    # システムモデル
    alpha = pm.GaussianRandomWalk('alpha', mu=0, sigma=sigma, init_dist=pm.Uniform.dist(0, 10), steps=N-1)
    
    # 観測モデル
    l = pm.Deterministic(f'lambda', alpha.exp())
    for i in range(N):
        Yobs = pm.Poisson(f'Y[{i+1}]', mu=l[i], observed=Y[i])

initvals = {}
with model:
    trace = pm.sample(1000, initvals=initvals, progressbar=False)

●実行結果

# 結果の可視化
import arviz as az
import matplotlib.pyplot as plt

x = range(1, N+1)
plt.scatter(x, Y, label='observed')
plt.ylim(-1, 25)

# ICARを用いた結果（比較用）
plt.plot(x, lambda_mean_CAR, color='r', linewidth=2, label='posterior lambda')
plt.fill_between(x, hdi_CAR['lambda'][:,0], hdi_CAR['lambda'][:,1], color='pink', alpha=0.15, label='95% credible interval')

# 状態遷移モデルを用いた結果
lambda_mean_RW = pm.summary(trace, var_names='lambda')['mean']
plt.plot(x, lambda_mean_RW, color='b', linewidth=2, label='posterior lambda')

# 95%ベイズ信用区間
hdi_RW = az.hdi(trace, hdi_prob=0.95, var_names='lambda')
plt.fill_between(x, hdi_RW['lambda'][:,0], hdi_RW['lambda'][:,1], color='cyan', alpha=0.15, label='95% credible interval')

plt.xlabel('position')
plt.ylabel('count')
plt.legend()
plt.show()

一応、Intrinsic CARの結果と状態空間モデルの結果がほぼ一致することが確認できた。

これでは納得できなかったので、PyMC Discourseというコミュニティで質問をしてみたら、pytensor.scanを使えば良い（そもそもPyMCのモデルではループを使うべきでない）という回答を得た。
AttributeError: 'Scratchpad' object has no attribute 'ufunc' - Questions / v5 - PyMC Discourse
この投稿者のmaruinenは筆者で、やんわりとPyTensorの不具合について指摘したつもりだったが、それについては誰からもコメントされなかった。

このやり取りにあるように、何とかpytensor.scanを用いた状態空間モデルの書き方を理解して、上のAttributeErrorが出た例を書き直して動かすことに成功した。

●コード

# 発見個体数からのパラメーター推定（pytensor.scanを使う例）
import pymc as pm
import pandas as pd
from pymc.pytensorf import collect_default_updates
import pytensor
import pytensor.tensor as pt

N = 50

# Helper function for pm.CustomDist
n_steps = N - 1
def statespace_dist(mu_init, sigma_level, size):

    def grw_step(mu_tm1, sigma_level):
        mu_t = mu_tm1 + pm.Normal.dist(sigma=sigma_level)
        return mu_t, collect_default_updates(outputs=[mu_t])

    mu, updates = pytensor.scan(fn=grw_step, 
                                outputs_info=[{"initial": mu_init}],
                                non_sequences=[sigma_level], 
                                n_steps=N-1,
                                name='statespace',
                                strict=True)

    return mu

with pm.Model() as model:
    # 事前分布
    alpha1 = pm.Uniform('alpha[1]', 0, 10)
    sigma = pm.Uniform('sigma', 0, 100)
    
    # システムモデル
    alpha = pm.CustomDist('alpha', alpha1, sigma, dist=statespace_dist)
    alpha = pt.concatenate([[alpha1], alpha])

    # 観測モデル
    l = pm.Deterministic('lambda', alpha.exp())
    Yobs = pm.Poisson('Y', mu=l, observed=Y)

with model:
    trace = pm.sample(1000, progressbar=False)

●実行結果

# 結果の可視化
import arviz as az
import matplotlib.pyplot as plt

x = range(1, N+1)
plt.scatter(x, Y, label='observed')
plt.ylim(-1, 25)

# ICARを用いた結果（比較用）
plt.plot(x, lambda_mean_CAR, color='r', linewidth=2, label='posterior lambda')
plt.fill_between(x, hdi_CAR['lambda'][:,0], hdi_CAR['lambda'][:,1], color='pink', alpha=0.15, label='95% credible interval')

# 状態遷移モデルを用いた結果
lambda_mean_SS = pm.summary(trace, var_names='lambda')['mean']
plt.plot(x, lambda_mean_SS, color='b', linewidth=2, label='posterior lambda')

# 95%ベイズ信用区間
hdi_SS = az.hdi(trace, hdi_prob=0.95, var_names='lambda')
plt.fill_between(x, hdi_SS['lambda'][:,0], hdi_SS['lambda'][:,1], color='cyan', alpha=0.15, label='95% credible interval')

plt.xlabel('position')
plt.ylabel('count')
plt.legend()
plt.show()

「岩波データサイエンス Vol.1 [特集] ベイズ推論とMCMCのフリーソフト」という本を借りたので読んでみると、PyMCというライブラリを使えばPythonでもできると書かれていたので、サンプルコードがBUGSやRで書かれている例題を含めて全てPyMCでやってみることにした。
BUGSやRのコードを読み解きながら、p.39〜p.46辺りの、時系列データに状態空間モデルのトレンドモデルや、観測モデルに対数正規分布を用いたモデルを当てはめてパラメーター推定する例題まで順調に進められたが、p.47〜p.48の「観測値が二項分布のモデル」でPyMCのトラブルに遭遇した。その回避方法がわかったので、ここに記す。

まず成功例として、p.45辺りの、観測モデルに対数正規分布を用いた例を紹介する。
なお、今回はなるべくこの本のサンプルコードと比較できるように、サンプルコードに近い形で書くことを目指している為、PyMCを上手く使う書き方にはしていない。例えば、PyMCではモデル作成時にループを使うことが推奨されていないようである。そもそも筆者はPyMCを使い始めたばかりで、ではどう書くのが良いのかがわからない。

import matplotlib.pyplot as plt
# 年輪データ
# https://github.com/iwanami-datascience/vol1/blob/master/ito/nenrin.data より借用

nenrin_data = [
    4.71, 7.7 , 7.97, 8.35, 5.7 , 7.33, 3.1 , 4.98, 3.75, 3.35, 1.84, 3.28, 2.77, 2.72, 2.54,
    3.23, 2.45, 1.9 , 2.56, 2.12, 1.78, 3.18, 2.64, 1.86, 1.69, 0.81, 1.02, 1.4 , 1.31, 1.57]
index = range(1960, 1990)

plt.plot(index, nenrin_data, color='k')
plt.ylabel('Tree-ring width')
plt.show()

こういう、京都市内で採取された杉の木の年輪幅の推移のデータがあり、これに次のような状態空間モデルを当てはめる。

α_tが年輪幅で、1行目はα_tのシステムモデルとしてトレンドモデルを仮定するもの、2行目は1行目を変形したもの、3行目は観測モデルである。 このα_tを次のコードで推定する。

●コード

import pymc as pm
import pandas as pd

N = nenrin_data.shape[0]

with pm.Model() as model:
    # 事前分布
    alpha = [pm.Uniform('alpha1', 0, 10), pm.Uniform('alpha2', 0, 10)]
    sigma = [pm.Uniform('sigma1', 0, 10), pm.Uniform('sigma2', 0, 10)]  # σ_ε, σ_η
             
    # システムモデル
    for i in range(2, N):
        alpha.append(pm.Normal(f'alpha{i+1}', mu=2 * alpha[i-1] - alpha[i-2], sigma=sigma[1]))

    # 観測モデル
    y = []
    for i in range(N):
        y.append(pm.LogNormal(f'y{i+1}', mu=alpha[i], sigma=sigma[0], observed=nenrin_data['width'].iloc[i]))

with model:
    trace = pm.sample(10000, progressbar=False)

●結果

# 結果の可視化
# alphaはyの対数に対応しているのでexp(alpha)に変換している
import numpy as np

plt.plot(index, nenrin_data, color='k', label='observed')

# 事後分布の平均
alpha_mean = [trace['posterior'][f'alpha{i+1}'].mean() for i in range(N)]
alpha_mean_exp = np.exp(alpha_mean)
plt.plot(index, alpha_mean_exp, color='b', linewidth=2, label='posterior')

# 95%ベイズ信用区間
alpha_025 = [np.percentile(trace['posterior'][f'alpha{i+1}'], 2.5) for i in range(N)]
alpha_975 = [np.percentile(trace['posterior'][f'alpha{i+1}'], 97.5) for i in range(N)]
alpha_025_exp = np.exp(alpha_025)
alpha_975_exp = np.exp(alpha_975)
plt.fill_between(index, alpha_025_exp, alpha_975_exp, color='cyan', alpha=0.15, label='95% credible interval')

plt.ylabel('Tree-ring width')
plt.legend()
plt.show()

この結果は本に記載されている図とほぼ同じである。

次に、1つ目の問題が発生した、p.47〜p.48の「観測値が二項分布のモデル」の例である。

ある林にいる鳥の数をカウントするとする。 鳥の個体を発見できる確率をp、1回の観測につき5回(i=1〜5)繰り返しカウントする、それを50回(t=1〜50)観測するとする。 実際には観測できない真の個体数N_tは期待値λ_tのポアソン分布に従うとする。
t回目の観測におけるi回目のカウントでの発見個体数をN_ti^obs とすると、観測モデルは次のようになる。

λ_tは対数スケールで、分散σ²の正規分布に従ってランダムウォークするとすると、システムモデルは次のようになる。

●データ

import numpy as np

# 模擬データ作成
np.random.seed(5)

Nt = 50
Nr = 5
p = 0.6
sigma = 0.1

# log λ_t
log_lambda = [np.log(30)]
for i in range(1, Nt):
    log_lambda.append(log_lambda[i-1] + np.random.normal(0, sigma))

# 真の個体数
N_true = [np.random.poisson(np.exp(log_lambda[i])) for i in range(Nt)] 

# 発見個体数
Nobs = np.array([np.random.binomial(N_true[t], p, Nr) for t in range(50)])

# 可視化
import matplotlib.pyplot as plt

plt.plot(range(1, Nt+1), N_true, color='k', linewidth=0.5, label='true')
plt.scatter([[t+1]*5 for t in range(50)], Nobs, marker='.', alpha=0.5, label='observed')

plt.xlabel('time')
plt.title('Bird count')
plt.legend(loc='upper center')
plt.show()

このデータから真の個体数を推定するとする。
上の成功例と同じ要領で、次のようなコードを実行すると、CompileErrorが発生した。

●コード

# 発見個体数からのパラメーター推定
import pymc as pm

# CompileError対策
# import pytensor
# pytensor.config.gcc__cxxflags += ' -fbracket-depth=384'

Nt = 50  # Nt >= 32だとCompileErrorになる
Nr = 5

p = None
log_lambda = None
sigma = None
N = None

with pm.Model() as model:
    # 事前分布
    p = pm.Uniform('p', 0, 1)
    log_lambda = [pm.Normal('loglambda[1]', 0, 100)]  # sigma=100は書籍だと1e-4だが、それだと結果がかなりずれる
    sigma = pm.Uniform('sigma', 0, 100)

    # システムモデル
    for t in range(1, Nt):
        log_lambda.append(pm.Normal(f'loglambda[{t+1}]', mu=log_lambda[t-1], sigma=sigma))

    # 観測モデル
    l = [pm.Deterministic(f'lambda[{t+1}]', log_lambda[t].exp()) for t in range(Nt)]
    N = [pm.Poisson(f'N[{t+1}]', mu=l[t]) for t in range(Nt)]
    Y = [pm.Binomial(f'Nobs[{t+1}]', n=N[t], p=p, observed=Nobs[t]) for t in range(Nt)]

# 初期値を与えないとSamplingErrorになる
initvals = {f'N[{t+1}]': 100 for t in range(Nt)}

with model:
    trace = pm.sample(10000, initvals=initvals)

●結果

CompileError: Compilation failed (return status=1):

/Users/ynomura/.pytensor/compiledir_macOS-14.5-arm64-arm-64bit-arm-3.10.14-64/tmp7dc45x_s/mod.cpp:26982:32: fatal error: bracket nesting level exceeded maximum of 256
        if (!PyErr_Occurred()) {
                               ^
/Users/ynomura/.pytensor/compiledir_macOS-14.5-arm64-arm-64bit-arm-3.10.14-64/tmp7dc45x_s/mod.cpp:26982:32: note: use -fbracket-depth=N to increase maximum nesting level

HINT: Use a linker other than the C linker to print the inputs' shapes and strides.
HINT: Re-running with most PyTensor optimizations disabled could provide a back-trace showing when this node was created. This can be done by setting the PyTensor flag 'optimizer=fast_compile'. If that does not work, PyTensor optimizations can be disabled with 'optimizer=None'.
HINT: Use the PyTensor flag `exception_verbosity=high` for a debug print-out and storage map footprint of this Apply node.

筆者の環境は、macOS 14.5 + python 3.10.14 + PyMC 5.15.0 + PyTensor 2.20.0 だった。PyMCはpipでインストールしたものである。
コードのコメントにも書いたが、色々試した結果、 Nt = 31 なら発生しないことがわかったので、Pythonのコードが間違っている訳ではなさそうだった。

Webで検索すると、このCompileErrorが発生したと書かれているページがいくつか見つかったが、具体的な解決方法が見つけられなかった。
PyMCのサイトにAnacondaやMiniforgeを使ってインストールすることを推奨すると書かれているので、Miniforgeを使ったインストールも試したが、結果は同じだった。その環境は python 3.12.3 + PyMC 5.15.1 + PyTensor 2.22.1だった。

上のエラーメッセージから、Cコンパイラの{}のネストが最大ネスト数の256を超えたというエラーが出ており、コンパイラの引数に-fbracket-depth=Nを加えると最大ネスト数が変えられると書いてある。色々調べた結果、次の2行（上のコードではコメントアウトしている）を加えると、このCompileErrorを回避できた。

import pytensor
pytensor.config.gcc__cxxflags += ' -fbracket-depth=384'

●結果

import matplotlib.pyplot as plt
import numpy as np

x = range(1, Nt+1)
plt.plot(x, N_true[:Nt], color='k', linewidth=0.5, label='true')
plt.scatter([[t+1]*5 for t in range(50)], Nobs, marker='.', alpha=0.5, label='observed')

lambda_mean = [trace['posterior'][f'lambda[{t+1}]'].mean() for t in range(Nt)]
plt.plot(x, lambda_mean, color='b', linewidth=2, label='posterior lambda')

# 95%ベイズ信用区間
lambda_025 = [np.percentile(trace['posterior'][f'lambda[{t+1}]'], 2.5)  for t in range(Nt)]
lambda_975 = [np.percentile(trace['posterior'][f'lambda[{t+1}]'], 97.5)  for t in range(Nt)]
plt.fill_between(x, lambda_025, lambda_975, color='cyan', alpha=0.15, label='95% credible interval')

plt.xlabel('time')
plt.title('Bird count')
plt.legend()
plt.show()