#  MindSearch 思·索:
Mimicking Human Minds Elicits Deep AI Searcher

Zehui Chen1∗, Kuikun Liu2∗, Qiuchen Wang1, Jiangning Liu2,
Wenwei Zhang2, Kai Chen2†, Feng Zhao1†

1 University of Science and Technology of China
2 Shanghai AI Laboratory

###### Abstract

Information seeking and integration is a complex cognitive task that consumes
enormous time and effort. Search engines reshape the way of seeking
information but often fail to align with complex human intentions. Inspired by
the remarkable progress of Large Language Models (LLMs), recent works attempt
to solve the information-seeking and integration task by combining LLMs and
search engines. However, these methods still obtain unsatisfying performance
due to three challenges: (1) complex requests often cannot be accurately and
completely retrieved by the search engine once; (2) corresponding information
to be integrated is spread over multiple web pages along with massive noise;
and (3) a large number of web pages with long contents may quickly exceed the
maximum context length of LLMs. Inspired by the cognitive process when humans
solve these problems, we introduce MindSearch (思·索) to mimic the human minds
in web information seeking and integration, which can be instantiated by a
simple yet effective LLM-based multi-agent framework consisting of a
WebPlanner and WebSearcher. The WebPlanner models the human mind of multi-step
information seeking as a dynamic graph construction process: it decomposes the
user query into atomic sub-questions as nodes in the graph and progressively
extends the graph based on the search result from WebSearcher. Tasked with
each sub-question, WebSearcher performs hierarchical information retrieval
with search engines and collects valuable information for WebPlanner. The
multi-agent design of MindSearch enables the whole framework to seek and
integrate information parallelly from larger-scale (e.g., more than 300) web
pages in 3 minute, which is worth 3 hours of human effort. Based on either
GPT-4o or InternLM2.5-7B models, MindSearch demonstrates significant
improvement in the response quality in terms of depth and breadth, on both
closed-set and open-set QA problems. Besides, responses from MindSearch based
on InternLM2.5-7B are preferable by humans to ChatGPT-Web (by GPT-4o) and
Perplexity.ai applications, which implies that MindSearch with open-source
models can already deliver a competitive solution to the proprietary AI search
engine. Code and models are available at
https://github.com/InternLM/MindSearch.

††Project Page: https://mindsearch.netlify.app††∗ Equal Contribution.
Alphabetical order † Corresponding author Figure 1: The overall framework of
MindSearch. It consists of two main ingredients: WebPlanner and WebSearcher.
WebPlanner acts as a high-level planner, orchestrating the reasoning steps and
multiple WebSearchers. WebSearcher conducts fine-grained web searches and
summarizes valuable information back to the planner, formalizing a simple yet
effective multi-agent framework.

##  1 Introduction

Information seeking and integration is a necessary cognitive process before
analysis and decision-making in all walks of life, which usually consumes
enormous human efforts and time. The birth of search engines (Brin & Page,
1998; Berkhin, 2005) significantly has reshaped and eased the information-
seeking process of human society, however, it still suffers in integrating web
information based on complex human intentions. Recently, Large Language Models
(LLMs) have showcased remarkable progress in reasoning, language
understanding, and information integration across a variety of domains (Achiam
et al., 2023; Team et al., 2024; Touvron et al., 2023; Cai et al., 2024),
whereas they struggling to deliver accurate knowledge in responses (Ji et al.,
2023; Gu et al., 2024).

The complementary advantages of LLMs and search engines highlights a
compelling opportunity for their combination, where the reasoning prowess of
LLMs can be complemented by the extensive web information accessible via
search engines, potentially revolutionizing the solution of web information
seeking and integration. Previous works (Asai et al., 2023; Chan et al., 2024)
simply treat the information seeking and integration task as a vanilla
retrieve-augmented generation (RAG) task (Chen et al., 2017; Lin et al.,
2023). Such a formulation, although straightforward, often results in sub-
optimal performance due to a superficial engagement with the depth and
complexity of web-based information retrieval, facing three major challenges
for more complex user queries:

(1) Real-world problems often require in-depth analysis and proper
decomposition of the question before retrieving the related information, which
cannot be done by retrieving web pages at once.
(2) The overwhelming volume of searched web pages and massive information
noise pose great challenges for LLMs for efficient information integration.
(3) The rapid proliferation of web search content can quickly exceed the
maximum context length of LLMs, which further decreases the information
integration performance.

Inspired by how human experts solve real-world problems, we propose MindSearch
思·索111The Chinese name ’思·索‘ means thinking as human and exploring by
searching, a simple yet effective LLM-based multi-agent framework, which
consists of a WebPlanner (mimic human minds for problem reasoning) and
multiple WebSearcher (manage the information seeking process). Given a user
query, the WebPlanner first decomposes the query into multiple atomic sub-
questions that can be parallelly solved and dispatches them to the respective
WebSearcher. To further enhance the reasoning ability, WebPlanner models the
complex problem-solving process as an iterative graph construction: by
predefining a list of standard code interfaces related to the construction of
the topological mind graph, WebPlanner is able to progressively decompose the
question into sequential/parallel sub-problems by adding nodes/edges in the
graph via Python code generation. Meanwhile, the WebSearcher, tasked with each
sub-problem, employs a hierarchical retrieval process to extract valuable data
for LLMs, which significantly improves the information aggregation efficiency
facing massive search pages. By distributing different aspects of the
reasoning and retrieval process to specialized agents, MindSearch effectively
reduces the load on each single agent, facilitating a more robust handling of
long contexts. It seamlessly bridges the gap between the raw data retrieval
capabilities of search engines and the context-understanding power of LLMs.

To validate the effectiveness of MindSearch, we conducted extensive
evaluations on both closed-set and open-set question-answering (QA) problems
using GPT-4o and InternLM2.5-7B-Chat models. The experimental results
demonstrate a substantial improvement in response quality, both in the
dimensions of depth and breadth. Moreover, comparative analysis shows that the
responses of MindSearch are more preferred by human evaluators over those from
existing applications like ChatGPT-Web (based on GPT-4o) and Perplexity Pro.
These findings suggest that MindSearch with open-source LLMs can offer a
highly competitive solution for AI-driven search engines.

##  2 MindSearch

To effectively synergize the web information retrieval capabilities of search
engines and the reasoning and information integration capability of LLMs,
MindSearch consists of a WebPlanner and a group of WebSearchers (Fig. 1).
WebPlanner first decomposes the user question into sequential or parallel
search tasks via reasoning on the graph and determines the next step based on
the search feedback (Sec. 2.1). WebSearcher is tasked with the query and
performs hierarchical information retrieval on the Internet to answer sub-
questions (Sec. 2.2). We also discuss the context management within the scope
of the multi-agent design in Sec. 2.3.

Figure 2: A concrete example of how WebPlanner addresses the question step by
step via planning as coding. During each turn, WebPlanner outputs a series of
thoughts along with the generated code. The code will be executed and yield
the search results to the planner. At the last turn, the WebPlanner directly
provides the final response without any code generation.

###  2.1 WebPlanner: Planning via Graph Construction

The WebPlanner functions as a high-level planner, orchestrating the reasoning
steps and coordinating other agents. However, we observed that merely
prompting the LLM to plan the entire data workflow architecture does not yield
satisfactory performance. Specifically, current LLMs struggle with decomposing
complex questions and understanding their topological relationships, leading
to coarse-grained search queries. This approach underutilizes the potential of
LLMs to serve as intermediaries between humans and search engines,
transforming human intentions into step-by-step search tasks and delivering
accurate responses.

To enhance the capability of LLM in addressing complex questions, we model the
problem-solving process as a directed acyclic graph (DAG). Given a user
question $Q$, the solution trajectory is represented as $G(Q)=\langle
V,E\rangle$, where $V$ is a set of nodes $v$, each representing an independent
web search, including an auxiliary START node (the initial question) and an
END node (the final answer). $E$ represents directed edges indicating the
reasoning topological relationships between nodes (search contents). This DAG
formalism captures the complexity of finding the optimal execution path,
providing a more formal and intuitive representation for LLMs.

Leveraging the superior performance of current LLMs on code tasks (Guo et al.,
2024; Roziere et al., 2023), we explicitly prompt the model to interact with
the graph through code writing. To achieve this, we predefined atomic code
functions to add nodes or edges to the graph (Step 1 and 2 in Figure 2). At
each turn, the LLM first reads the entire dialogue, including previously
generated code and web search results, then outputs thoughts and new code for
reasoning on the mind graph, which is executed with a Python interpreter.
During execution, once a node is added to the reasoning graph, it invokes a
WebSearcher to execute the search process and summarize the information. Since
the newly added nodes are only dependent on nodes generated in previous steps,
we can parallel them to achieve a much faster information aggregation speed.
When all information is collected, the planner produces the final response by
adding the end node (Step 3 in Figure 2).

By integrating with the Python interpreter, WebPlanner interacts with the
graph through unified code actions, dynamically constructing the reasoning
path. This ”code as planning” process enables the LLM to fully leverage its
superior code generation ability, benefiting control and data flow in long-
context scenarios and leading to better performance in solving complex
problems.

###  2.2 WebSearcher: Web Browsing with Hierarchical Retrieval

WebSearcher acts as a sophisticated RAG (Retrieve-and-Generate) agent with
internet access, summarizing valuable responses based on search results
(Figure 3). Due to the massive content available on the web, it is challenging
for LLMs to process all related pages within a limited context length (e.g. 8K
tokens). To address this, we employ a straightforward coarse-to-fine selection
strategy. Initially, the LLM generates several similar queries based on the
assigned questions from the WebPlanner to broaden the search content and thus
improve the recall of relevant information. These queries are then executed
through various search APIs, such as Google, Bing, and DuckDuckGo, which
return key contents including web URLs, titles, and summaries. The search
results are automatically merged based on the web URLs, and the LLM is
prompted to select the most valuable pages for detailed reading. The full
content of the selected web URLs is then added to the input of LLM. After
reading these results, the LLM generates a response to answer the original
question based on the search results. This hierarchical retrieval approach
significantly reduces the difficulty of navigating massive web pages and
allows to efficiently extract highly relevant information with in-depth
details.

Figure 3: A detailed working pipeline of WebSearcher. It comprises at most 4
steps: query rewrite, search content aggregation, detailed page selection, and
final summarization.

###  2.3 LLM Context Management in MindSearch

MindSearch provides a simple multi-agent solution to complex information
seeking and integration with search engines. Such a paradigm also naturally
enables long-context management among different agents, which improves the
overall efficiency of the framework, especially under circumstances that
require the model to quickly read plenty of web pages. Since the WebPlanner
distributes the search tasks into separate search agents and only relies on
the searched results from WebSearcher, WebPlanner can purely focus on the
decomposition and analysis of the user question without being distracted by
the over-length web search results. Meanwhile, each WebSearcher only needs to
search contents for its tasked sub-query, without distraction from other
contents. Thanks to the explicit role distribution, MindSearch greatly reduces
context computation during the whole process, delivering an efficient context
management solution to long-context tasks for LLM. Such a multi-agent
framework also provides a straightforward and simple long-context task
construction pipeline for training single LLMs, which is also observed in
(Team, 2024). Eventually, MindSearch collects and integrates related
information from more than 300 pages in less than 3 minute, which could take
human experts about 3 hours to finish a similar cognitive workload.

Due to the explicit context state transfer across multiple agents, we need to
carefully handle the context during the whole workflow. We empirically find
simply focusing the decomposed query from the Planner may lose useful
information during the information collection phase due to the local receptive
field inside the search agent. How to effectively handle the context between
multiple agents is non-trivial. We find that the constructed topological
relations through the directed graph edges help us easily handle the context
across different agents. More specifically, we simply prefix the response from
its father node as well as the root node when executing each search agent.
Therefore, each WebSearcher can effectively focus on its sub-task without
losing the previous related context as well as the final goal.

##  3 Experiments

We evaluate MindSearch on two primary categories of Question Answering (QA)
tasks: closed-set QA and open-set QA, which reflects both the subjective and
objective judgment of MindSearch. For a fair comparison, all models only have
access to the Internet through BING search API, and no extra reference sources
are considered.

###  3.1 Open-Set QA

####  3.1.1 Implementation Details

To better gauge the utility and search performance, we carefully curate 100
real-world human queries and collect responses from MindSearch
(InternLM2.5-7b-chat (Cai et al., 2024)), Perplexity.ai (its Pro version) ,
and ChatGPT with search plugin Achiam et al. (2023). We ask five human experts
to manually select their preferred responses, in terms of the following three
aspects:

  * •

Depth: Depth refers to the thoroughness and profundity of an answer. A
response with depth provides detailed information and delves into the
intricacies of a question.

  * •

Breadth: Breadth pertains to the scope and diversity covered by an answer. A
response with breadth touches on various aspects of the question or multiple
related fields, offering different perspectives or solutions.

  * •

Factuality: Factuality is the degree to which an answer is accurate and fact-
based. It should be grounded in reliable data and information, avoiding errors
or misleading content, and ensuring the truthfulness and credibility of the
information provided.

The final results are determined based on major votes. During the evaluation,
the correspondence between the response and its method is invisible to the
evaluators to guarantee fairness.

Figure 4: Subjective evaluation results judged by human experts on open-set QA
questions. MindSearch outperforms ChatGPT-Web and Perplexity.ai Pro by a large
margin in terms of depth, breadth, and facticity.

####  3.1.2 Results and Analysis

The evaluation results are depicted in Figure 4 and we also provide
quantitative results in Figure 5. From Figure 4, we can observe an absolute
improvement in terms of the depth and breadth of the model response, which
validates the superiority of our proposed WebPlanner. By integrating code into
the DAG construction phase, LLM is able to progressively decompose the complex
problem into executable queries while balancing the tradeoff between time
efficiency and the exploration of the search space. Besides, MindSearch goes
through more fine-grained search topics about the question, therefore
providing more compact and detailed responses compared to other models.
However, MindSearch does not yield better performance in terms of facticity.
We suspect that more detailed search results may distract the concentration of
the model on the initial problem, especially when LLM holds incomplete long-
context capability. Therefore, a natural future work of MindSearch is to
alleviate the hallucination issues during the web browsing process.

Figure 5: Solution trajectory comparison between MindSearch and Perplexity.ai (Pro) on the same question. MindSearch provides more detailed and proper responses thanks to its fine-grained searches. Table 1: Performance comparison on various closed-set QA tasks. We select two representative LLMs: GPT-4o (close-sourced) and InternLM2.5-7b-chat (open-sourced). Model | Bamboogle | Musique | HotpotQA | AVG
---|---|---|---|---
2-hop | 3-hop | 4-hop | Easy | Medium | Hard |
Closed-Souced LLM (GPT-4o)
w/o Search Engine | 70.4 | 54.0 | 22.0 | 20.0 | 73.0 | 69.0 | 66.0 | 53.5
ReAct Search | 75.2 | 48.0 | 25.0 | 13.3 | 81.0 | 73.0 | 70.0 | 55.1
MindSearch | 76.8 | 60.0 | 35.0 | 14.6 | 80.0 | 74.0 | 78.0 | 59.8
Open-Sourced LLM (InternLM2.5-7b-chat)
w/o Search Engine | 34.0 | 28.0 | 10.0 | 17.3 | 47.0 | 26.0 | 40.0 | 28.9
ReAct Search | 55.2 | 38.0 | 17.0 | 16.0 | 69.0 | 56.0 | 49.0 | 42.9
MindSearch | 67.8 | 46.0 | 20.0 | 18.6 | 69.0 | 66.0 | 57.0 | 49.2

In addition to quantitative results, we also provide a qualitative response
comparison between Perplexity.ai (Pro) and MindSearch to deliver an intuitive
understanding of their performance. From Figure 5, we can observe that
MindSearch yields more concrete and detailed responses. We empirically find
that our better responses can be attributed to the proper planning search
paths compared to Perplexity.ai, which also indicates that how to decompose
the human intention is the key step to the final problem.

###  3.2 Closed-Set QA

####  3.2.1 Implementation Details

We extensively evaluate our approach on a wide range of closed-set QA tasks,
including Bamboogle (Press et al., 2022), Musique (Trivedi et al., 2022), and
HotpotQA (Yang et al., 2018). To further validate the generalization of our
approach, we select both closed-source LLM (GPT-4o) and open-source LLM
(InternLM2.5-7b-chat) as our LLM backend. Since our approach adopts a zero-
shot experimental setting, we utilize a subjective LLM evaluator (GPT4-o) to
gauge the correctness of HotpotQA.

####  3.2.2 Results and Analysis

In Table 1, we compare our approach with two straight-forward baselines: raw
LLM without search engines (w/o Search Engine), and simply treating search
engines as an external tool and adopting a ReAct-style interaction (ReAct
Search). We can conclude that MindSearch significantly outperforms its vanilla
baselines by a large margin, validating the effectiveness of the proposed
method. These advantages are amplified when transferring from closed-sourced
LLMs to open-sourced LLMs, which further proves that MindSeach provides a
simple approach to enhance weak LLMs with broader knowledge and alleviate
hallucination issues.

##  4 Related Work

###  4.1 Tool Utilization with LLM

The Tool Learning framework empowers LLMs to seamlessly integrate with a
variety of tools (Qin et al., 2023; Hao et al., 2024; Zhuang et al., 2024;
Chen et al., 2023), such as search engines (Chan et al., 2024), databases
(Parisi et al., 2022), and APIs (Li et al., 2023; Patil et al., 2023),
offering dynamic solutions to complex problems. This integration is not only
beneficial for enhancing the interpretability and trustworthiness of LLMs but
also for improving their robustness and adaptability across diverse tasks,
including reducing hallucinations (Ji et al., 2024), code generation (Gou et
al., 2023), and question answering (Chen et al., 2024). Recent research has
focused on enhancing the tool integration component of Tool Learning systems.
Works such as (Huang et al., 2023; Shen et al., 2023; Schick et al., 2024)
have concentrated on improving the retrieval mechanisms, ensuring that LLMs
can access the most pertinent tools for a given task. Other studies, like
(Qian et al., 2023; Yuan et al., 2023), aim at refining the LLMs’ ability to
effectively utilize the retrieved information, optimizing the reading and
comprehension processes within the framework.

###  4.2 RAG with LLM

RAG demonstrates significant advantages in addressing knowledge-intensive
problems, especially in open-domain scenarios with the integration of search
engines (Chen et al., 2017; Li et al., 2017). RAG allows LLMs to integrate
with the retriever, providing timely information and offering effective
solutions. Moreover, RAG is also applied in various tasks such as reducing
hallucinations (Shuster et al., 2021; Gu et al., 2024), code generation (Zhou
et al., 2022), and question answering (Lewis et al., 2020). Recently, some
work (Karpukhin et al., 2020; Xiong et al., 2020; Qu et al., 2020) focuses on
enhancing the retrieval component of RAG systems, while others (Izacard &
Grave, 2020; Borgeaud et al., 2022; Yu et al., 2021; Lei et al., 2017)
enhances the language model’s ability as a reader to optimize the framework.

With the advancement of LLM capabilities, some researchers have begun to
reoptimize frameworks and redesign methodologies for model training. SAIL (Luo
et al., 2023) trains LLM to be more focused on credible and informative search
results. Self-RAG (Asai et al., 2023) enables LMMs to independently fetch,
introspect, and augment their text generation capabilities. RQ-RAG (Chan et
al., 2024) enhances query formulation by learning to refine queries through an
iterative process. Our work integrates web search capabilities into LLMs,
enhancing response quality by retrieving valuable information from the
Internet.

###  4.3 Web Agents

Web automation agents have evolved from question-answering tools to
sophisticated systems capable of complex web interactions. Early models like
WebGPT (Nakano et al., 2021) and WebGLM (Liu et al., 2023) primarily addressed
QA tasks, while recent advancements have shifted towards more dynamic
operations (Yao et al., 2022; He et al., 2024). MindAct (Deng et al., 2024)
and WebAgent (Gur et al., 2023) represent this progression, with the latter
showing exceptional web navigation despite deployment challenges due to its
size. AutoWebGLM (Lai et al., 2024) offers a practical alternative with robust
capabilities and a more compact model size. The incorporation of reinforcement
learning (Bai et al., 2024) and behavior cloning techniques (Zheng et al.,
2024; Patel et al., 2024) paves the way for even more autonomous and efficient
web automation, moving the field towards scalable and versatile solutions for
real-world applications. This paper mainly focuses more on the web
information-seeking and integration task with search engines instead of web
browsing, and solves the main challenges with a multi-agent framework.

##  5 Conclusion

This paper introduces MindSearch, a novel LLM-based multi-agent framework for
complex web information-seeking and integration tasks, by more comprehensively
leveraging the strengths of both search engines and LLMs. MindSearch conducts
effective and sufficient decomposition of complex queries followed by
hierarchical information retrieval to improve the precision and recall of the
retrieved relevant web information, by modeling the problem-solving process as
an iterative graph construction. The multi-agent design distributes the
cognitive load among specialized agents, facilitating robust handling of
complex and lengthy contexts. Extensive evaluations on closed-set and open-set
QA problems using GPT-4o and InternLM2.5-7B models demonstrated significant
advantages in the response quality of MindSearch. The results that human
evaluators preferred the responses from MindSearch over those from ChatGPT-Web
and Perplexity.ai indicate its competitive edge in AI-driven search solutions.
We wish this work pave the way for future research on multi-agent framework
for solving human-level complex cognitive tasks.

##  6 Acknowledgement

We would like to express our sincere gratitude to Jiaye Ge for her outstanding
technical insights, product design skills, and coordination abilities. We also
thank Zhongying Tu, Ying Zhao, Fang Fang, and Yiting Wang for their efficient
execution in developing the project demo. Our gratitude extends to Xingyuan
Liu, Shuaike Li, Zike Pan, Weijia Song, and Yuzhe Gu for their efforts in
project suggestion.

## References

  * Achiam et al. (2023) Josh Achiam, Steven Adler, Sandhini Agarwal, Lama Ahmad, Ilge Akkaya, Florencia Leoni Aleman, Diogo Almeida, Janko Altenschmidt, Sam Altman, Shyamal Anadkat, et al.  Gpt-4 technical report.  _arXiv preprint arXiv:2303.08774_ , 2023.
  * Asai et al. (2023) Akari Asai, Zeqiu Wu, Yizhong Wang, Avirup Sil, and Hannaneh Hajishirzi.  Self-rag: Learning to retrieve, generate, and critique through self-reflection.  _arXiv preprint arXiv:2310.11511_ , 2023.
  * Bai et al. (2024) Hao Bai, Yifei Zhou, Mert Cemri, Jiayi Pan, Alane Suhr, Sergey Levine, and Aviral Kumar.  Digirl: Training in-the-wild device-control agents with autonomous reinforcement learning.  _arXiv preprint arXiv:2406.11896_ , 2024.
  * Berkhin (2005) Pavel Berkhin.  A survey on pagerank computing.  _Internet mathematics_ , 2(1):73–120, 2005.
  * Borgeaud et al. (2022) Sebastian Borgeaud, Arthur Mensch, Jordan Hoffmann, Trevor Cai, Eliza Rutherford, Katie Millican, George Bm Van Den Driessche, Jean-Baptiste Lespiau, Bogdan Damoc, Aidan Clark, et al.  Improving language models by retrieving from trillions of tokens.  In _International conference on machine learning_ , pp. 2206–2240. PMLR, 2022.
  * Brin & Page (1998) Sergey Brin and Lawrence Page.  The anatomy of a large-scale hypertextual web search engine.  _Computer networks and ISDN systems_ , 30(1-7):107–117, 1998.
  * Cai et al. (2024) Zheng Cai, Maosong Cao, Haojiong Chen, Kai Chen, Keyu Chen, Xin Chen, Xun Chen, Zehui Chen, Zhi Chen, Pei Chu, et al.  Internlm2 technical report.  _arXiv preprint arXiv:2403.17297_ , 2024.
  * Chan et al. (2024) Chi-Min Chan, Chunpu Xu, Ruibin Yuan, Hongyin Luo, Wei Xue, Yike Guo, and Jie Fu.  Rq-rag: Learning to refine queries for retrieval augmented generation.  _arXiv preprint arXiv:2404.00610_ , 2024.
  * Chen et al. (2017) Danqi Chen, Adam Fisch, Jason Weston, and Antoine Bordes.  Reading wikipedia to answer open-domain questions.  _arXiv preprint arXiv:1704.00051_ , 2017.
  * Chen et al. (2023) Zehui Chen, Weihua Du, Wenwei Zhang, Kuikun Liu, Jiangning Liu, Miao Zheng, Jingming Zhuo, Songyang Zhang, Dahua Lin, Kai Chen, et al.  T-eval: Evaluating the tool utilization capability step by step.  _arXiv preprint arXiv:2312.14033_ , 2023.
  * Chen et al. (2024) Zehui Chen, Kuikun Liu, Qiuchen Wang, Wenwei Zhang, Jiangning Liu, Dahua Lin, Kai Chen, and Feng Zhao.  Agent-flan: Designing data and methods of effective agent tuning for large language models.  _arXiv preprint arXiv:2403.12881_ , 2024.
  * Deng et al. (2024) Xiang Deng, Yu Gu, Boyuan Zheng, Shijie Chen, Sam Stevens, Boshi Wang, Huan Sun, and Yu Su.  Mind2web: Towards a generalist agent for the web.  _Advances in Neural Information Processing Systems_ , 36, 2024.
  * Gou et al. (2023) Zhibin Gou, Zhihong Shao, Yeyun Gong, Yujiu Yang, Minlie Huang, Nan Duan, Weizhu Chen, et al.  Tora: A tool-integrated reasoning agent for mathematical problem solving.  _arXiv preprint arXiv:2309.17452_ , 2023.
  * Gu et al. (2024) Yuzhe Gu, Ziwei Ji, Wenwei Zhang, Chengqi Lyu, Dahua Lin, and Kai Chen.  Anah-v2: Scaling analytical hallucination annotation of large language models.  _arXiv preprint arXiv:2407.04693_ , 2024.
  * Guo et al. (2024) Daya Guo, Qihao Zhu, Dejian Yang, Zhenda Xie, Kai Dong, Wentao Zhang, Guanting Chen, Xiao Bi, Yu Wu, YK Li, et al.  Deepseek-coder: When the large language model meets programming–the rise of code intelligence.  _arXiv preprint arXiv:2401.14196_ , 2024.
  * Gur et al. (2023) Izzeddin Gur, Hiroki Furuta, Austin Huang, Mustafa Safdari, Yutaka Matsuo, Douglas Eck, and Aleksandra Faust.  A real-world webagent with planning, long context understanding, and program synthesis.  _arXiv preprint arXiv:2307.12856_ , 2023.
  * Hao et al. (2024) Shibo Hao, Tianyang Liu, Zhen Wang, and Zhiting Hu.  Toolkengpt: Augmenting frozen language models with massive tools via tool embeddings.  _Advances in neural information processing systems_ , 36, 2024.
  * He et al. (2024) Hongliang He, Wenlin Yao, Kaixin Ma, Wenhao Yu, Yong Dai, Hongming Zhang, Zhenzhong Lan, and Dong Yu.  Webvoyager: Building an end-to-end web agent with large multimodal models.  _arXiv preprint arXiv:2401.13919_ , 2024.
  * Huang et al. (2023) Yue Huang, Jiawen Shi, Yuan Li, Chenrui Fan, Siyuan Wu, Qihui Zhang, Yixin Liu, Pan Zhou, Yao Wan, Neil Zhenqiang Gong, et al.  Metatool benchmark for large language models: Deciding whether to use tools and which to use.  _arXiv preprint arXiv:2310.03128_ , 2023.
  * Izacard & Grave (2020) Gautier Izacard and Edouard Grave.  Leveraging passage retrieval with generative models for open domain question answering.  _arXiv preprint arXiv:2007.01282_ , 2020.
  * Ji et al. (2023) Ziwei Ji, Nayeon Lee, Rita Frieske, Tiezheng Yu, Dan Su, Yan Xu, Etsuko Ishii, Ye Jin Bang, Andrea Madotto, and Pascale Fung.  Survey of hallucination in natural language generation.  _ACM Computing Surveys_ , 55(12):1–38, 2023.
  * Ji et al. (2024) Ziwei Ji, Yuzhe Gu, Wenwei Zhang, Chengqi Lyu, Dahua Lin, and Kai Chen.  Anah: Analytical annotation of hallucinations in large language models.  _arXiv preprint arXiv:2405.20315_ , 2024.
  * Karpukhin et al. (2020) Vladimir Karpukhin, Barlas Oğuz, Sewon Min, Patrick Lewis, Ledell Wu, Sergey Edunov, Danqi Chen, and Wen-tau Yih.  Dense passage retrieval for open-domain question answering.  _arXiv preprint arXiv:2004.04906_ , 2020.
  * Lai et al. (2024) Hanyu Lai, Xiao Liu, Iat Long Iong, Shuntian Yao, Yuxuan Chen, Pengbo Shen, Hao Yu, Hanchen Zhang, Xiaohan Zhang, Yuxiao Dong, et al.  Autowebglm: Bootstrap and reinforce a large language model-based web navigating agent.  _arXiv preprint arXiv:2404.03648_ , 2024.
  * Lei et al. (2017) Xiangyu Lei, Guilin Zhang, Shuaijun Li, Huihuan Qian, and Yangsheng Xu.  Dual-spring agv shock absorption system design: Dynamic analysis and simulations.  In _2017 IEEE International Conference on Robotics and Biomimetics (ROBIO)_ , pp. 1068–1074. IEEE, 2017.
  * Lewis et al. (2020) Patrick Lewis, Ethan Perez, Aleksandra Piktus, Fabio Petroni, Vladimir Karpukhin, Naman Goyal, Heinrich Küttler, Mike Lewis, Wen-tau Yih, Tim Rocktäschel, et al.  Retrieval-augmented generation for knowledge-intensive nlp tasks.  _Advances in Neural Information Processing Systems_ , 33:9459–9474, 2020.
  * Li et al. (2023) Minghao Li, Yingxiu Zhao, Bowen Yu, Feifan Song, Hangyu Li, Haiyang Yu, Zhoujun Li, Fei Huang, and Yongbin Li.  Api-bank: A comprehensive benchmark for tool-augmented llms.  _arXiv preprint arXiv:2304.08244_ , 2023.
  * Li et al. (2017) Shuaijun Li, Guilin Zhang, Xiangyu Lei, Xiao Yu, Huihuan Qian, and Yangsheng Xu.  Trajectory tracking control of a unicycle-type mobile robot with a new planning algorithm.  In _2017 IEEE International Conference on Robotics and Biomimetics (ROBIO)_ , pp. 780–786. IEEE, 2017.
  * Lin et al. (2023) Xi Victoria Lin, Xilun Chen, Mingda Chen, Weijia Shi, Maria Lomeli, Rich James, Pedro Rodriguez, Jacob Kahn, Gergely Szilvasy, Mike Lewis, et al.  Ra-dit: Retrieval-augmented dual instruction tuning.  _arXiv preprint arXiv:2310.01352_ , 2023.
  * Liu et al. (2023) Xiao Liu, Hanyu Lai, Hao Yu, Yifan Xu, Aohan Zeng, Zhengxiao Du, Peng Zhang, Yuxiao Dong, and Jie Tang.  Webglm: Towards an efficient web-enhanced question answering system with human preferences.  In _Proceedings of the 29th ACM SIGKDD Conference on Knowledge Discovery and Data Mining_ , pp. 4549–4560, 2023.
  * Luo et al. (2023) Hongyin Luo, Yung-Sung Chuang, Yuan Gong, Tianhua Zhang, Yoon Kim, Xixin Wu, Danny Fox, Helen Meng, and James Glass.  Sail: Search-augmented instruction learning.  _arXiv preprint arXiv:2305.15225_ , 2023.
  * Nakano et al. (2021) Reiichiro Nakano, Jacob Hilton, Suchir Balaji, Jeff Wu, Long Ouyang, Christina Kim, Christopher Hesse, Shantanu Jain, Vineet Kosaraju, William Saunders, et al.  Webgpt: Browser-assisted question-answering with human feedback.  _arXiv preprint arXiv:2112.09332_ , 2021.
  * Parisi et al. (2022) Aaron Parisi, Yao Zhao, and Noah Fiedel.  Talm: Tool augmented language models.  _arXiv preprint arXiv:2205.12255_ , 2022.
  * Patel et al. (2024) Ajay Patel, Markus Hofmarcher, Claudiu Leoveanu-Condrei, Marius-Constantin Dinu, Chris Callison-Burch, and Sepp Hochreiter.  Large language models can self-improve at web agent tasks.  _arXiv preprint arXiv:2405.20309_ , 2024.
  * Patil et al. (2023) Shishir G Patil, Tianjun Zhang, Xin Wang, and Joseph E Gonzalez.  Gorilla: Large language model connected with massive apis.  _arXiv preprint arXiv:2305.15334_ , 2023.
  * Press et al. (2022) Ofir Press, Muru Zhang, Sewon Min, Ludwig Schmidt, Noah A Smith, and Mike Lewis.  Measuring and narrowing the compositionality gap in language models.  _arXiv preprint arXiv:2210.03350_ , 2022.
  * Qian et al. (2023) Cheng Qian, Chi Han, Yi R Fung, Yujia Qin, Zhiyuan Liu, and Heng Ji.  Creator: Tool creation for disentangling abstract and concrete reasoning of large language models.  _arXiv preprint arXiv:2305.14318_ , 2023.
  * Qin et al. (2023) Yujia Qin, Shihao Liang, Yining Ye, Kunlun Zhu, Lan Yan, Yaxi Lu, Yankai Lin, Xin Cong, Xiangru Tang, Bill Qian, et al.  Toolllm: Facilitating large language models to master 16000+ real-world apis.  _arXiv preprint arXiv:2307.16789_ , 2023.
  * Qu et al. (2020) Yingqi Qu, Yuchen Ding, Jing Liu, Kai Liu, Ruiyang Ren, Wayne Xin Zhao, Daxiang Dong, Hua Wu, and Haifeng Wang.  Rocketqa: An optimized training approach to dense passage retrieval for open-domain question answering.  _arXiv preprint arXiv:2010.08191_ , 2020.
  * Roziere et al. (2023) Baptiste Roziere, Jonas Gehring, Fabian Gloeckle, Sten Sootla, Itai Gat, Xiaoqing Ellen Tan, Yossi Adi, Jingyu Liu, Tal Remez, Jérémy Rapin, et al.  Code llama: Open foundation models for code.  _arXiv preprint arXiv:2308.12950_ , 2023.
  * Schick et al. (2024) Timo Schick, Jane Dwivedi-Yu, Roberto Dessì, Roberta Raileanu, Maria Lomeli, Eric Hambro, Luke Zettlemoyer, Nicola Cancedda, and Thomas Scialom.  Toolformer: Language models can teach themselves to use tools.  _Advances in Neural Information Processing Systems_ , 36, 2024.
  * Shen et al. (2023) Yongliang Shen, Kaitao Song, Xu Tan, Wenqi Zhang, Kan Ren, Siyu Yuan, Weiming Lu, Dongsheng Li, and Yueting Zhuang.  Taskbench: Benchmarking large language models for task automation.  _arXiv preprint arXiv:2311.18760_ , 2023.
  * Shuster et al. (2021) Kurt Shuster, Spencer Poff, Moya Chen, Douwe Kiela, and Jason Weston.  Retrieval augmentation reduces hallucination in conversation.  _arXiv preprint arXiv:2104.07567_ , 2021.
  * Team et al. (2024) Gemma Team, Thomas Mesnard, Cassidy Hardin, Robert Dadashi, Surya Bhupatiraju, Shreya Pathak, Laurent Sifre, Morgane Rivière, Mihir Sanjay Kale, Juliette Love, et al.  Gemma: Open models based on gemini research and technology.  _arXiv preprint arXiv:2403.08295_ , 2024.
  * Team (2024) Qwen Team.  Generalizing an llm from 8k to 1m context using qwen-agent, May 2024.  URL https://qwenlm.github.io/blog/qwen-agent-2405/.
  * Touvron et al. (2023) Hugo Touvron, Louis Martin, Kevin Stone, Peter Albert, Amjad Almahairi, Yasmine Babaei, Nikolay Bashlykov, Soumya Batra, Prajjwal Bhargava, Shruti Bhosale, et al.  Llama 2: Open foundation and fine-tuned chat models.  _arXiv preprint arXiv:2307.09288_ , 2023.
  * Trivedi et al. (2022) Harsh Trivedi, Niranjan Balasubramanian, Tushar Khot, and Ashish Sabharwal.  Musique: Multihop questions via single-hop question composition.  _Transactions of the Association for Computational Linguistics_ , 10:539–554, 2022.
  * Xiong et al. (2020) Lee Xiong, Chenyan Xiong, Ye Li, Kwok-Fung Tang, Jialin Liu, Paul Bennett, Junaid Ahmed, and Arnold Overwijk.  Approximate nearest neighbor negative contrastive learning for dense text retrieval.  _arXiv preprint arXiv:2007.00808_ , 2020.
  * Yang et al. (2018) Zhilin Yang, Peng Qi, Saizheng Zhang, Yoshua Bengio, William W Cohen, Ruslan Salakhutdinov, and Christopher D Manning.  Hotpotqa: A dataset for diverse, explainable multi-hop question answering.  _arXiv preprint arXiv:1809.09600_ , 2018.
  * Yao et al. (2022) Shunyu Yao, Howard Chen, John Yang, and Karthik Narasimhan.  Webshop: Towards scalable real-world web interaction with grounded language agents.  _Advances in Neural Information Processing Systems_ , 35:20744–20757, 2022.
  * Yu et al. (2021) Donghan Yu, Chenguang Zhu, Yuwei Fang, Wenhao Yu, Shuohang Wang, Yichong Xu, Xiang Ren, Yiming Yang, and Michael Zeng.  Kg-fid: Infusing knowledge graph in fusion-in-decoder for open-domain question answering.  _arXiv preprint arXiv:2110.04330_ , 2021.
  * Yuan et al. (2023) Lifan Yuan, Yangyi Chen, Xingyao Wang, Yi R Fung, Hao Peng, and Heng Ji.  Craft: Customizing llms by creating and retrieving from specialized toolsets.  _arXiv preprint arXiv:2309.17428_ , 2023.
  * Zheng et al. (2024) Boyuan Zheng, Boyu Gou, Jihyung Kil, Huan Sun, and Yu Su.  Gpt-4v (ision) is a generalist web agent, if grounded.  _arXiv preprint arXiv:2401.01614_ , 2024.
  * Zhou et al. (2022) Shuyan Zhou, Uri Alon, Frank F Xu, Zhiruo Wang, Zhengbao Jiang, and Graham Neubig.  Docprompting: Generating code by retrieving the docs.  _arXiv preprint arXiv:2207.05987_ , 2022.
  * Zhuang et al. (2024) Yuchen Zhuang, Yue Yu, Kuan Wang, Haotian Sun, and Chao Zhang.  Toolqa: A dataset for llm question answering with external tools.  _Advances in Neural Information Processing Systems_ , 36, 2024.


# UQ_Handbook_Chapter

John Salvatier, Thomas V. Wiecki, Christopher Fonnesbeck

# Probabilistic Programming in Python using PyMC

John Salvatier, Thomas V. Wiecki, Christopher Fonnesbeck

##  1 Introduction

Probabilistic programming (PP) allows flexible specification of Bayesian
statistical models in code. PyMC3 is a new, open-source PP framework with an
intutive and readable, yet powerful, syntax that is close to the natural
syntax statisticians use to describe models. It features next-generation
Markov chain Monte Carlo (MCMC) sampling algorithms such as the No-U-Turn
Sampler (NUTS; Hoffman, 2014), a self-tuning variant of Hamiltonian Monte
Carlo (HMC; Duane, 1987). This class of samplers works well on high
dimensional and complex posterior distributions and allows many complex models
to be fit without specialized knowledge about fitting algorithms. HMC and NUTS
take advantage of gradient information from the likelihood to achieve much
faster convergence than traditional sampling methods, especially for larger
models. NUTS also has several self-tuning strategies for adaptively setting
the tunable parameters of Hamiltonian Monte Carlo, whicstatisticalh means you
usually don’t need to have specialized knowledge about how the algorithms
work. PyMC3, Stan (Stan Development Team, 2014), and the LaplacesDemon package
for R are currently the only PP packages to offer HMC.

Probabilistic programming in Python confers a number of advantages including
multi-platform compatibility, an expressive yet clean and readable syntax,
easy integration with other scientific libraries, and extensibility via C,
C++, Fortran or Cython. These features make it relatively straightforward to
write and use custom statistical distributions, samplers and transformation
functions, as required by Bayesian analysis.

While most of PyMC3’s user-facing features are written in pure Python, it
leverages Theano (Bergstra et al., 2010) to transparently transcode models to
C and compile them to machine code, thereby boosting performance. Theano is a
library that allows expressions to be defined using generalized vector data
structures called _tensors_ , which are tightly integrated with the popular
NumPy ndarray data structure, and similarly allow for broadcasting and
advanced indexing, just as NumPy arrays do. Theano also automatically
optimizes the likelihood’s computational graph for speed and provides simple
GPU integration.

Here, we present a primer on the use of PyMC3 for solving general Bayesian
statistical inference and prediction problems. We will first see the basics of
how to use PyMC3, motivated by a simple example: installation, data creation,
model definition, model fitting and posterior analysis. Then we will cover two
case studies and use them to show how to define and fit more sophisticated
models. Finally we will show how to extend PyMC3 and discuss other useful
features: the Generalized Linear Models subpackage, custom distributions,
custom transformations and alternative storage backends.

##  2 Installation

Running PyMC3 requires a working Python interpreter, either version 2.7 (or
more recent) or 3.4 (or more recent); we recommend that new users install
version 3.4. A complete Python installation for Mac OSX, Linux and Windows can
most easily be obtained by downloading and installing the free Anaconda Python
Distribution by ContinuumIO.

PyMC3 can be installed using pip
(https://pip.pypa.io/en/latest/installing.html):

    pip install git+https://github.com/pymc-devs/pymc3

PyMC3 depends on several third-party Python packages which will be
automatically installed when installing via pip. The four required
dependencies are: Theano, NumPy, SciPy, and Matplotlib.

To take full advantage of PyMC3, the optional dependencies Pandas and Patsy
should also be installed. These are _not_ automatically installed, but can be
installed by:

    pip install patsy pandas

The source code for PyMC3 is hosted on GitHub at https://github.com/pymc-
devs/pymc3 and is distributed under the liberal Apache License 2.0. On the
GitHub site, users may also report bugs and other issues, as well as
contribute code to the project, which we actively encourage.

##  3 A Motivating Example: Linear Regression

To introduce model definition, fitting and posterior analysis, we first
consider a simple Bayesian linear regression model with normal priors for the
parameters. We are interested in predicting outcomes $Y$ as normally-
distributed observations with an expected value $\mu$ that is a linear
function of two predictor variables, $X_{1}$ and $X_{2}$.

$\displaystyle Y$ $\displaystyle\sim\mathcal{N}(\mu,\sigma^{2})$
$\displaystyle\mu$ $\displaystyle=\alpha+\beta_{1}X_{1}+\beta_{2}X_{2}$

where $\alpha$ is the intercept, and $\beta_{i}$ is the coefficient for
covariate $X_{i}$, while $\sigma$ represents the observation error. Since we
are constructing a Bayesian model, the unknown variables in the model must be
assigned a prior distribution. We choose zero-mean normal priors with variance
of 100 for both regression coefficients, which corresponds to _weak_
information regarding the true parameter values. We choose a half-normal
distribution (normal distribution bounded at zero) as the prior for $\sigma$.

$\displaystyle\alpha$ $\displaystyle\sim\mathcal{N}(0,100)$
$\displaystyle\beta_{i}$ $\displaystyle\sim\mathcal{N}(0,100)$
$\displaystyle\sigma$ $\displaystyle\sim\lvert\mathcal{N}(0,1){\rvert}$

###  3.1 Generating data

We can simulate some artificial data from this model using only NumPy’s random
module, and then use PyMC3 to try to recover the corresponding parameters. We
are intentionally generating the data to closely correspond the PyMC3 model
structure.

,commandchars=
{},,

{\color{incolor}In [{\color{incolor}1}]:} \PY{k+kn}{import} \PY{n+nn}{numpy}
\PY{k+kn}{as} \PY{n+nn}{np} \PY{k+kn}{import} \PY{n+nn}{pylab} \PY{k+kn}{as}
\PY{n+nn}{pl} \PY{c}{\PYZsh{} Intialize random number generator}
\PY{n}{np}\PY{o}{.}\PY{n}{random}\PY{o}{.}\PY{n}{seed}\PY{p}{(}\PY{l+m+mi}{123}\PY{p}{)}
\PY{c}{\PYZsh{} True parameter values} \PY{n}{alpha}\PY{p}{,} \PY{n}{sigma}
\PY{o}{=} \PY{l+m+mi}{1}\PY{p}{,} \PY{l+m+mi}{1} \PY{n}{beta} \PY{o}{=}
\PY{p}{[}\PY{l+m+mi}{1}\PY{p}{,} \PY{l+m+mf}{2.5}\PY{p}{]} \PY{c}{\PYZsh{}
Size of dataset} \PY{n}{size} \PY{o}{=} \PY{l+m+mi}{100} \PY{c}{\PYZsh{}
Predictor variable} \PY{n}{X1} \PY{o}{=}
\PY{n}{np}\PY{o}{.}\PY{n}{linspace}\PY{p}{(}\PY{l+m+mi}{0}\PY{p}{,}
\PY{l+m+mi}{1}\PY{p}{,} \PY{n}{size}\PY{p}{)} \PY{n}{X2} \PY{o}{=}
\PY{n}{np}\PY{o}{.}\PY{n}{linspace}\PY{p}{(}\PY{l+m+mi}{0}\PY{p}{,}\PY{o}{.}\PY{l+m+mi}{2}\PY{p}{,}
\PY{n}{size}\PY{p}{)} \PY{c}{\PYZsh{} Simulate outcome variable} \PY{n}{Y}
\PY{o}{=} \PY{n}{alpha} \PY{o}{+}
\PY{n}{beta}\PY{p}{[}\PY{l+m+mi}{0}\PY{p}{]}\PY{o}{*}\PY{n}{X1} \PY{o}{+}
\PY{n}{beta}\PY{p}{[}\PY{l+m+mi}{1}\PY{p}{]}\PY{o}{*}\PY{n}{X2} \PY{o}{+}
\PY{n}{np}\PY{o}{.}\PY{n}{random}\PY{o}{.}\PY{n}{randn}\PY{p}{(}\PY{n}{size}\PY{p}{)}\PY{o}{*}\PY{n}{sigma}

Here is what the simulated data look like. We use the pylab module from the
plotting library matplotlib.

,commandchars=
{},,

{\color{incolor}In [{\color{incolor}2}]:} \PY{o}{\PYZpc{}}\PY{k}{pylab} inline
\PY{n}{fig}\PY{p}{,} \PY{n}{axes} \PY{o}{=}
\PY{n}{subplots}\PY{p}{(}\PY{l+m+mi}{1}\PY{p}{,} \PY{l+m+mi}{2}\PY{p}{,}
\PY{n}{sharex}\PY{o}{=}\PY{n+nb+bp}{True}\PY{p}{,}
\PY{n}{figsize}\PY{o}{=}\PY{p}{(}\PY{l+m+mi}{10}\PY{p}{,}\PY{l+m+mi}{4}\PY{p}{)}\PY{p}{)}
\PY{n}{axes}\PY{p}{[}\PY{l+m+mi}{0}\PY{p}{]}\PY{o}{.}\PY{n}{scatter}\PY{p}{(}\PY{n}{X1}\PY{p}{,}
\PY{n}{Y}\PY{p}{)}
\PY{n}{axes}\PY{p}{[}\PY{l+m+mi}{1}\PY{p}{]}\PY{o}{.}\PY{n}{scatter}\PY{p}{(}\PY{n}{X2}\PY{p}{,}
\PY{n}{Y}\PY{p}{)}
\PY{n}{axes}\PY{p}{[}\PY{l+m+mi}{0}\PY{p}{]}\PY{o}{.}\PY{n}{set\PYZus{}ylabel}\PY{p}{(}\PY{l+s}{\PYZsq{}}\PY{l+s}{Y}\PY{l+s}{\PYZsq{}}\PY{p}{)}\PY{p}{;}
\PY{n}{axes}\PY{p}{[}\PY{l+m+mi}{0}\PY{p}{]}\PY{o}{.}\PY{n}{set\PYZus{}xlabel}\PY{p}{(}\PY{l+s}{\PYZsq{}}\PY{l+s}{X1}\PY{l+s}{\PYZsq{}}\PY{p}{)}\PY{p}{;}
\PY{n}{axes}\PY{p}{[}\PY{l+m+mi}{1}\PY{p}{]}\PY{o}{.}\PY{n}{set\PYZus{}xlabel}\PY{p}{(}\PY{l+s}{\PYZsq{}}\PY{l+s}{X2}\PY{l+s}{\PYZsq{}}\PY{p}{)}\PY{p}{;}

,commandchars=
{},,

Populating the interactive namespace from numpy and matplotlib

,commandchars=
{},,

WARNING: pylab import has clobbered these variables: [’beta’,
’size’]‘\%matplotlib‘ prevents importing * from pylab and numpy

###  3.2 Model Specification

Specifiying this model in PyMC3 is straightforward because the syntax is as
close to the statistical notation. For the most part, each line of Python code
corresponds to a line in the model notation above.

First, we import the components we will need from PyMC.

,commandchars=
{},,

{\color{incolor}In [{\color{incolor}3}]:} \PY{k+kn}{from} \PY{n+nn}{pymc3}
\PY{k+kn}{import} \PY{n}{Model}\PY{p}{,} \PY{n}{Normal}\PY{p}{,}
\PY{n}{HalfNormal}

Now we build our model, which we will present in full first, then explain each
part line-by-line.

,commandchars=
{},,

{\color{incolor}In [{\color{incolor}4}]:} \PY{n}{basic\PYZus{}model} \PY{o}{=}
\PY{n}{Model}\PY{p}{(}\PY{p}{)} \PY{k}{with}
\PY{n}{basic\PYZus{}model}\PY{p}{:} \PY{c}{\PYZsh{} Priors for unknown model
parameters} \PY{n}{alpha} \PY{o}{=}
\PY{n}{Normal}\PY{p}{(}\PY{l+s}{\PYZsq{}}\PY{l+s}{alpha}\PY{l+s}{\PYZsq{}}\PY{p}{,}
\PY{n}{mu}\PY{o}{=}\PY{l+m+mi}{0}\PY{p}{,}
\PY{n}{sd}\PY{o}{=}\PY{l+m+mi}{10}\PY{p}{)} \PY{n}{beta} \PY{o}{=}
\PY{n}{Normal}\PY{p}{(}\PY{l+s}{\PYZsq{}}\PY{l+s}{beta}\PY{l+s}{\PYZsq{}}\PY{p}{,}
\PY{n}{mu}\PY{o}{=}\PY{l+m+mi}{0}\PY{p}{,}
\PY{n}{sd}\PY{o}{=}\PY{l+m+mi}{10}\PY{p}{,}
\PY{n}{shape}\PY{o}{=}\PY{l+m+mi}{2}\PY{p}{)} \PY{n}{sigma} \PY{o}{=}
\PY{n}{HalfNormal}\PY{p}{(}\PY{l+s}{\PYZsq{}}\PY{l+s}{sigma}\PY{l+s}{\PYZsq{}}\PY{p}{,}
\PY{n}{sd}\PY{o}{=}\PY{l+m+mi}{1}\PY{p}{)} \PY{c}{\PYZsh{} Expected value of
outcome} \PY{n}{mu} \PY{o}{=} \PY{n}{alpha} \PY{o}{+}
\PY{n}{beta}\PY{p}{[}\PY{l+m+mi}{0}\PY{p}{]}\PY{o}{*}\PY{n}{X1} \PY{o}{+}
\PY{n}{beta}\PY{p}{[}\PY{l+m+mi}{1}\PY{p}{]}\PY{o}{*}\PY{n}{X2}
\PY{c}{\PYZsh{} Likelihood (sampling distribution) of observations}
\PY{n}{Y\PYZus{}obs} \PY{o}{=}
\PY{n}{Normal}\PY{p}{(}\PY{l+s}{\PYZsq{}}\PY{l+s}{Y\PYZus{}obs}\PY{l+s}{\PYZsq{}}\PY{p}{,}
\PY{n}{mu}\PY{o}{=}\PY{n}{mu}\PY{p}{,}
\PY{n}{sd}\PY{o}{=}\PY{n}{sigma}\PY{p}{,}
\PY{n}{observed}\PY{o}{=}\PY{n}{Y}\PY{p}{)}

The first line,

commandchars=
{},,,

\NormalTok{basic_model = Model()}

creates a new Model object which is a container for the model random
variables.

Following instantiation of the model, the subsequent specification of the
model components is performed inside a with statement:

commandchars=
{},,,

\KeywordTok{with} \NormalTok{basic_model:}

This creates a _context manager_ , with our basic_model as the context, that
includes all statements until the indented block ends. This means all PyMC3
objects introduced in the indented code block below the with statement are
added to the model behind the scenes. Absent this context manager idiom, we
would be forced to manually associate each of the variables with basic_model
right after we create them. If you try to create a new random variable without
a with model: statement, it will raise an error since there is no obvious
model for the variable to be added to.

The first three statements in the context manager:

commandchars=
{},,,

\NormalTok{alpha = Normal(}\StringTok{’alpha’}\NormalTok{,
mu=}\DecValTok{0}\NormalTok{, sd=}\DecValTok{10}\NormalTok{)}\NormalTok{beta =
Normal(}\StringTok{’beta’}\NormalTok{, mu=}\DecValTok{0}\NormalTok{,
sd=}\DecValTok{10}\NormalTok{,
shape=}\DecValTok{2}\NormalTok{)}\NormalTok{sigma =
HalfNormal(}\StringTok{’sigma’}\NormalTok{, sd=}\DecValTok{1}\NormalTok{)}

create a stochastic random variables with a Normal prior distributions for the
regression coefficients with a mean of 0 and standard deviation of 10 for the
regression coefficients, and a half-normal distribution for the standard
deviation of the observations, $\sigma$. These are stochastic because their
values are partly determined by its parents in the dependency graph of random
variables, which for priors are simple constants, and partly random (or
stochastic).

We call the Normal constructor to create a random variable to use as a normal
prior. The first argument is always the _name_ of the random variable, which
should almost always match the name of the Python variable being assigned to,
since it sometimes used to retrieve the variable from the model for
summarizing output. The remaining required arguments for a stochastic object
are the parameters, in this case mu, the mean, and sd, the standard deviation,
which we assign hyperparameter values for the model. In general, a
distribution’s parameters are values that determine the location, shape or
scale of the random variable, depending on the parameterization of the
distribution. Most commonly used distributions, such as Beta, Exponential,
Categorical, Gamma, Binomial and many others, are available in PyMC3.

The beta variable has an additional shape argument to denote it as a vector-
valued parameter of size 2. The shape argument is available for all
distributions and specifies the length or shape of the random variable, but is
optional for scalar variables, since it defaults to a value of one. It can be
an integer, to specify an array, or a tuple, to specify a multidimensional
array (_e.g._ shape=(5,7) makes random variable that takes on 5 by 7 matrix
values).

Detailed notes about distributions, sampling methods and other PyMC3 functions
are available via the help function.

,commandchars=
{},,

{\color{incolor}In [{\color{incolor}5}]:}
\PY{n}{help}\PY{p}{(}\PY{n}{Normal}\PY{p}{)} \PY{c}{\PYZsh{}try help(Model),
help(Uniform) or help(basic\PYZus{}model)}

,commandchars=
{},,

Help on class Normal in module pymc3.distributions.continuous:class Normal(pymc3.distributions.distribution.Continuous) | Normal log-likelihood. | | .. math::ight\textbackslash{}\\} | | Parameters | ---------- | mu : float | Mean of the distribution. | tau : float | Precision of the distribution, which corresponds to | :math:‘1/\textbackslash{}sigma\^{}2‘ (tau > 0). | sd : float | Standard deviation of the distribution. Alternative parameterization. | | .. note:: | - :math:‘E(X) = \textbackslash{}mu‘ | - :math:‘Var(X) = 1/ au‘ | | Method resolution order: | Normal | pymc3.distributions.distribution.Continuous | pymc3.distributions.distribution.Distribution | \\_\\_builtin\\_\\_.object | | Methods defined here: | | \\_\\_init\\_\\_(self, mu=0.0, tau=None, sd=None, *args, **kwargs) | | logp(self, value) | | ---------------------------------------------------------------------- | Methods inherited from pymc3.distributions.distribution.Distribution: | | \\_\\_getnewargs\\_\\_(self) | | default(self) | | get\\_test\\_val(self, val, defaults) | | getattr\\_value(self, val) | | ---------------------------------------------------------------------- | Class methods inherited from pymc3.distributions.distribution.Distribution: | | dist(cls, *args, **kwargs) from \\_\\_builtin\\_\\_.type | | ---------------------------------------------------------------------- | Static methods inherited from pymc3.distributions.distribution.Distribution: | | \\_\\_new\\_\\_(cls, name, *args, **kwargs) | | ---------------------------------------------------------------------- | Data descriptors inherited from pymc3.distributions.distribution.Distribution: | | \\_\\_dict\\_\\_ | dictionary for instance variables (if defined) | | \\_\\_weakref\\_\\_ | list of weak references to the object (if defined)

Having defined the priors, the next statement creates the expected value mu of
the outcomes, specifying the linear relationship:

commandchars=
{},,,

\NormalTok{mu = alpha + beta[}\DecValTok{0}\NormalTok{]*X1 +
beta[}\DecValTok{1}\NormalTok{]*X2}

This creates a deterministic random variable, which implies that its value is
_completely_ determined by its parents’ values. That is, there is no
uncertainty beyond that which is inherent in the parents’ values. Here, mu is
just the sum of the intercept alpha and the two products of the coefficients
in beta and the predictor variables, whatever their values may be.

PyMC3 random variables and data can be arbitrarily added, subtracted, divided,
multiplied together and indexed-into to create new random variables. This
allows for great model expressivity. Many common mathematical functions like
sum, sin, exp and linear algebra functions like dot (for inner product) and
inv (for inverse) are also provided.

The final line of the model, defines Y_obs, the sampling distribution of the
outcomes in the dataset.

commandchars=
{},,,

\NormalTok{Y_obs = Normal(}\StringTok{’Y_obs’}\NormalTok{, mu=mu, sd=sigma,
observed=Y)}

This is a special case of a stochastic variable that we call an observed
stochastic, and represents the data likelihood of the model. It is identical
to a standard stochastic, except that its observed argument, which passes the
data to the variable, indicates that the values for this variable were
observed, and should not be changed by any fitting algorithm applied to the
model. The data can be passed in the form of either a numpy.ndarray or
pandas.DataFrame object.

Notice that, unlike for the priors of the model, the parameters for the normal
distribution of Y_obs are not fixed values, but rather are the deterministic
object mu and the stochastic sigma. This creates parent-child relationships
between the likelihood and these two variables.

###  3.3 Model fitting

Having completely specified our model, the next step is to obtain posterior
estimates for the unknown variables in the model. Ideally, we could calculate
the posterior estimates analytically, but for most non-trivial models, this is
not feasible. We will consider two approaches, whose appropriateness depends
on the structure of the model and the goals of the analysis: finding the
_maximum a posteriori_ (MAP) point using optimization methods, and computing
summaries based on samples drawn from the posterior distribution using Markov
Chain Monte Carlo (MCMC) sampling methods.

#### Maximum a posteriori methods

The maximum a posteriori (MAP) estimate for a model, is the mode of the
posterior distribution and is generally found using numerical optimization
methods. This is often fast and easy to do, but only gives a point estimate
for the parameters and can be biased if the mode isn’t representative of the
distribution. PyMC3 provides this functionality with the find_MAP function.

Below we find the MAP for our original model. The MAP is returned as a
parameter point, which is always represented by a Python dictionary of
variable names to NumPy arrays of parameter values.

,commandchars=
{},,

{\color{incolor}In [{\color{incolor}6}]:} \PY{k+kn}{from} \PY{n+nn}{pymc3}
\PY{k+kn}{import} \PY{n}{find\PYZus{}MAP} \PY{n}{map\PYZus{}estimate}
\PY{o}{=}
\PY{n}{find\PYZus{}MAP}\PY{p}{(}\PY{n}{model}\PY{o}{=}\PY{n}{basic\PYZus{}model}\PY{p}{)}
\PY{k}{print}\PY{p}{(}\PY{n}{map\PYZus{}estimate}\PY{p}{)}

,commandchars=
{},,

\\{’alpha’: array(1.0136640995128503), ’beta’: array([ 1.46791595,
0.29358319]), ’sigma\\_log’: array(0.11928764983495602)\\}

By default, find_MAP uses the Broyden–Fletcher–Goldfarb–Shanno (BFGS)
optimization algorithm to find the maximum of the log-posterior but also
allows selection of other optimization algorithms from the scipy.optimize
module. For example, below we use Powell’s method to find the MAP.

,commandchars=
{},,

{\color{incolor}In [{\color{incolor}7}]:} \PY{k+kn}{from} \PY{n+nn}{scipy}
\PY{k+kn}{import} \PY{n}{optimize} \PY{n}{map\PYZus{}estimate} \PY{o}{=}
\PY{n}{find\PYZus{}MAP}\PY{p}{(}\PY{n}{model}\PY{o}{=}\PY{n}{basic\PYZus{}model}\PY{p}{,}
\PY{n}{fmin}\PY{o}{=}\PY{n}{optimize}\PY{o}{.}\PY{n}{fmin\PYZus{}powell}\PY{p}{)}
\PY{k}{print}\PY{p}{(}\PY{n}{map\PYZus{}estimate}\PY{p}{)}

,commandchars=
{},,

\\{’alpha’: array(1.0175522115056725), ’beta’: array([ 1.51426781,
0.03520891]), ’sigma\\_log’: array(0.1181510683418693)\\}

It is important to note that the MAP estimate is not always reasonable,
especially if the mode is at an extreme. This can be a subtle issue; with high
dimensional posteriors, one can have areas of extremely high density but low
total probability because the volume is very small. This will often occur in
hierarchical models with the variance parameter for the random effect. If the
individual group means are all the same, the posterior will have near infinite
density if the scale parameter for the group means is almost zero, even though
the probability of such a small scale parameter will be small since the group
means must be extremely close together.

Most techniques for finding the MAP estimate also only find a _local_ optimium
(which is often good enough), but can fail badly for multimodal posteriors if
the different modes are meaningfully different.

#### Sampling methods

Though finding the MAP is a fast and easy way of obtaining estimates of the
unknown model parameters, it is limited because there is no associated
estimate of uncertainty produced with the MAP estimates. Instead, a
simulation-based approach such as Markov chain Monte Carlo (MCMC) can be used
to obtain a Markov chain of values that, given the satisfaction of certain
conditions, are indistinguishable from samples from the posterior
distribution.

To conduct MCMC sampling to generate posterior samples in PyMC3, we specify a
step method object that corresponds to a particular MCMC algorithm, such as
Metropolis, Slice sampling, or the No-U-Turn Sampler (NUTS). PyMC3’s
step_methods submodule contains the following samplers: NUTS, Metropolis,
Slice, HamiltonianMC, and BinaryMetropolis.

#### Gradient-based sampling methods

PyMC3 has the standard sampling algorithms like adaptive Metropolis-Hastings
and adaptive slice sampling, but PyMC3’s most capable step method is the No-U-
Turn Sampler. NUTS is especially useful on models that have many continuous
parameters, a situatiuon where other MCMC algorithms work very slowly. It
takes advantage of information about where regions of higher probability are,
based on the gradient of the log posterior-density. This helps it achieve
dramatically faster convergence on large problems than traditional sampling
methods achieve. PyMC3 relies on Theano to analytically compute model
gradients via automatic differentation of the posterior density. NUTS also has
several self-tuning strategies for adaptively setting the tunable parameters
of Hamiltonian Monte Carlo. For random variables that are undifferentiable
(namely, discrete variables) NUTS cannot be used, but it may still be used on
the differentiable variables in a model that contains undifferentiable
variables.

NUTS requires a scaling matrix parameter, which is analogous to the variance
parameter for the jump proposal distribution in Metropolis-Hastings, althrough
NUTS uses it somewhat differently. The matrix gives the rough shape of the
distribution so that NUTS does not make jumps that are too large in some
directions and too small in other directions. It is important to set this
scaling parameter to a reasonable value to facilitate efficient sampling. This
is especially true for models that have many unobserved stochastic random
variables or models with highly non-normal posterior distributions. Poor
scaling parameters will slow down NUTS significantly, sometimes almost
stopping it completely. A reasonable starting point for sampling can also be
important for efficient sampling, but not as often.

Fortunately NUTS can often make good guesses for the scaling parameters. If
you pass a point in parameter space (as a dictionary of variable names to
parameter values, the same format as returned by find_MAP) to NUTS, it will
look at the local curvature of the log posterior-density (the diagonal of the
Hessian matrix) at that point to make a guess for a good scaling vector, which
often results in a good value. The MAP estimate is often a good point to use
to initiate sampling. It is also possible to supply your own vector or scaling
matrix to NUTS, though this is a more advanced use. If you wish to modify a
Hessian at a specific point to use as your scaling matrix or vector, you can
use find_hessian or find_hessian_diag.

For our basic linear regression example in basic_model, we will use NUTS to
sample 2000 draws from the posterior using the MAP as the starting point and
scaling point. This must also be performed inside the context of the model.

,commandchars=
{},,

{\color{incolor}In [{\color{incolor}8}]:} \PY{k+kn}{from} \PY{n+nn}{pymc3}
\PY{k+kn}{import} \PY{n}{NUTS}\PY{p}{,} \PY{n}{sample} \PY{k}{with}
\PY{n}{basic\PYZus{}model}\PY{p}{:} \PY{c}{\PYZsh{} obtain starting values via
MAP} \PY{n}{start} \PY{o}{=}
\PY{n}{find\PYZus{}MAP}\PY{p}{(}\PY{n}{fmin}\PY{o}{=}\PY{n}{optimize}\PY{o}{.}\PY{n}{fmin\PYZus{}powell}\PY{p}{)}
\PY{c}{\PYZsh{} instantiate sampler} \PY{n}{step} \PY{o}{=}
\PY{n}{NUTS}\PY{p}{(}\PY{n}{scaling}\PY{o}{=}\PY{n}{start}\PY{p}{)}
\PY{c}{\PYZsh{} draw 2000 posterior samples} \PY{n}{trace} \PY{o}{=}
\PY{n}{sample}\PY{p}{(}\PY{l+m+mi}{2000}\PY{p}{,} \PY{n}{step}\PY{p}{,}
\PY{n}{start}\PY{o}{=}\PY{n}{start}\PY{p}{)}

,commandchars=
{},,

[-----------------100\%-----------------] 2000 of 2000 complete in 10.4 sec

,commandchars=
{},,

/Library/Python/2.7/site-
packages/theano/scan\\_module/scan\\_perform\\_ext.py:133: RuntimeWarning:
numpy.ndarray size changed, may indicate binary incompatibility from
scan\\_perform.scan\\_perform import *

The sample function runs the step method(s) passed to it for the given number
of iterations and returns a Trace object containing the samples collected, in
the order they were collected. The trace object can be queried in a similar
way to a dict containing a map from variable names to numpy.arrays. The first
dimension of the array is the sampling index and the later dimensions match
the shape of the variable. We can see the last 5 values for the alpha variable
as follows

,commandchars=
{},,

{\color{incolor}In [{\color{incolor}9}]:}
\PY{n}{trace}\PY{p}{[}\PY{l+s}{\PYZsq{}}\PY{l+s}{alpha}\PY{l+s}{\PYZsq{}}\PY{p}{]}\PY{p}{[}\PY{o}{\PYZhy{}}\PY{l+m+mi}{5}\PY{p}{:}\PY{p}{]}

,commandchars=
{},,

{\color{outcolor}Out[{\color{outcolor}9}]:} array([ 0.93582185, 0.94703037,
1.0502649 , 0.91058163, 0.96962219])

###  3.4 Posterior analysis

PyMC3 provides plotting and summarization functions for inspecting the
sampling output. A simple posterior plot can be created using traceplot.

,commandchars=
{},,

{\color{incolor}In [{\color{incolor}10}]:} \PY{k+kn}{from} \PY{n+nn}{pymc3}
\PY{k+kn}{import} \PY{n}{traceplot}
\PY{n}{traceplot}\PY{p}{(}\PY{n}{trace}\PY{p}{)}\PY{p}{;}

The left column consists of a smoothed histogram (using kernel density
estimation) of the marginal posteriors of each stochastic random variable
while the right column contains the samples of the Markov chain plotted in
sequential order. The beta variable, being vector-valued, produces two
histograms and two sample traces, corresponding to both predictor
coefficients.

In addition, the summary function provides a text-based output of common
posterior statistics:

,commandchars=
{},,

{\color{incolor}In [{\color{incolor}11}]:} \PY{k+kn}{from} \PY{n+nn}{pymc3}
\PY{k+kn}{import} \PY{n}{summary}
\PY{n}{summary}\PY{p}{(}\PY{n}{trace}\PY{p}{)}

,commandchars=
{},,

alpha: Mean SD MC Error 95\% HPD interval
\------------------------------------------------------------------- 1.027
0.227 0.007 [0.592, 1.487] Posterior quantiles: 2.5 25 50 75 97.5
|--------------|==============|==============|--------------| 0.552 0.878
1.031 1.185 1.456beta: Mean SD MC Error 95\% HPD interval
\------------------------------------------------------------------- 1.340
1.861 0.099 [-1.997, 4.799] 0.838 9.188 0.488 [-16.926, 16.472] Posterior
quantiles: 2.5 25 50 75 97.5
|--------------|==============|==============|--------------| -1.994 0.027
1.413 2.540 4.805 -16.534 -5.052 0.713 7.316 16.939sigma\\_log: Mean SD MC
Error 95\% HPD interval
\------------------------------------------------------------------- 0.132
0.068 0.002 [-0.002, 0.258] Posterior quantiles: 2.5 25 50 75 97.5
|--------------|==============|==============|--------------| 0.002 0.084
0.130 0.176 0.269sigma: Mean SD MC Error 95\% HPD interval
\------------------------------------------------------------------- 1.143
0.079 0.002 [0.998, 1.294] Posterior quantiles: 2.5 25 50 75 97.5
|--------------|==============|==============|--------------| 1.002 1.088
1.139 1.193 1.308

##  4 Case study 1: Stochastic volatility

We present a case study of stochastic volatility, time varying stock market
volatility, to illustrate PyMC3’s use in addressing a more realistic problem.
The distribution of market returns is highly non-normal, which makes sampling
the volatlities significantly more difficult. This example has 400+ parameters
so using common sampling algorithms like Metropolis-Hastings would get bogged
down, generating highly autocorrelated samples. Instead, we use NUTS, which is
dramatically more efficient.

###  4.1 The Model

Asset prices have time-varying volatility (variance of day over day returns).
In some periods, returns are highly variable, while in others they are very
stable. Stochastic volatility models address this with a latent volatility
variable, which changes over time. The following model is similar to the one
described in the NUTS paper (Hoffman 2014, p. 21).

$\displaystyle\sigma$ $\displaystyle\sim exp(50)$ $\displaystyle\nu$
$\displaystyle\sim exp(.1)$ $\displaystyle s_{i}$
$\displaystyle\sim\mathcal{N}(s_{i-1},\sigma^{-2})$ $\displaystyle log(y_{i})$
$\displaystyle\sim t(\nu,0,exp(-2s_{i}))$

Here, $y$ is the daily return series which is modeled with a Student-t
distribution with an unknown degrees of freedom parameter, and a scale
parameter determined by a latent process $s$. The individual $s_{i}$ are the
individual daily log volatilities in the latent log volatility process.

###  4.2 The Data

Our data consist of the last 400 daily returns of the S&P 500.

,commandchars=
{},,

{\color{incolor}In [{\color{incolor}12}]:} \PY{n}{n} \PY{o}{=}
\PY{l+m+mi}{400} \PY{n}{returns} \PY{o}{=}
\PY{n}{np}\PY{o}{.}\PY{n}{genfromtxt}\PY{p}{(}\PY{l+s}{\PYZdq{}}\PY{l+s}{data/SP500.csv}\PY{l+s}{\PYZdq{}}\PY{p}{)}\PY{p}{[}\PY{o}{\PYZhy{}}\PY{n}{n}\PY{p}{:}\PY{p}{]}
\PY{n}{pl}\PY{o}{.}\PY{n}{plot}\PY{p}{(}\PY{n}{returns}\PY{p}{)}\PY{p}{;}

###  4.3 Model Specification

As with the linear regession example, specifying the model in PyMC3 mirrors
its statistical specification. This model employs several new distributions:
the Exponential distribution for the $\nu$ and $\sigma$ priors, the student-t
(T) distribution for distribution of returns, and the GaussianRandomWalk for
the prior for the latent volatilities.

In PyMC3, variables with purely positive priors like Exponential are
transformed with a log transform. This makes sampling more robust. Behind the
scenes, a variable in the unconstrained space (named “variableName_log”) is
added to the model for sampling. In this model his happens behind the scenes
for both the degrees of freedom, nu, and the scale parameter for the
volatility process, sigma, since they both have exponential priors. Variables
with priors that constrain them on two sides, like Beta or Uniform, are also
transformed to be unconstrained but with a log odds transform.

Although, unlike model specifiation in PyMC2, we do not typically provide
starting points for variables at the model specification stage, we can also
provide an initial value for any distribution (called a “test value”) using
the testval argument. This overrides the default test value for the
distribution (usually the mean, median or mode of the distribution), and is
most often useful if some values are illegal and we want to ensure we select a
legal one. The test values for the distributions are also used as a starting
point for sampling and optimization by default, though this is easily
overriden.

The vector of latent volatilities s is given a prior distribution by
GaussianRandomWalk. As its name suggests GaussianRandomWalk is a vector valued
distribution where the values of the vector form a random normal walk of
length n, as specified by the shape argument. The scale of the innovations of
the random walk, sigma, is specified in terms of the precision of the normally
distributed innovations and can be a scalar or vector.

,commandchars=
{},,

{\color{incolor}In [{\color{incolor}14}]:} \PY{k+kn}{from} \PY{n+nn}{pymc3}
\PY{k+kn}{import} \PY{n}{Exponential}\PY{p}{,} \PY{n}{T}\PY{p}{,}
\PY{n}{exp}\PY{p}{,} \PY{n}{Deterministic} \PY{k+kn}{from}
\PY{n+nn}{pymc3.distributions.timeseries} \PY{k+kn}{import}
\PY{n}{GaussianRandomWalk} \PY{k}{with} \PY{n}{Model}\PY{p}{(}\PY{p}{)}
\PY{k}{as} \PY{n}{sp500\PYZus{}model}\PY{p}{:} \PY{n}{nu} \PY{o}{=}
\PY{n}{Exponential}\PY{p}{(}\PY{l+s}{\PYZsq{}}\PY{l+s}{nu}\PY{l+s}{\PYZsq{}}\PY{p}{,}
\PY{l+m+mf}{1.}\PY{o}{/}\PY{l+m+mi}{10}\PY{p}{,}
\PY{n}{testval}\PY{o}{=}\PY{o}{.}\PY{l+m+mi}{1}\PY{p}{)} \PY{n}{sigma}
\PY{o}{=}
\PY{n}{Exponential}\PY{p}{(}\PY{l+s}{\PYZsq{}}\PY{l+s}{sigma}\PY{l+s}{\PYZsq{}}\PY{p}{,}
\PY{l+m+mf}{1.}\PY{o}{/}\PY{o}{.}\PY{l+m+mo}{02}\PY{p}{,}
\PY{n}{testval}\PY{o}{=}\PY{o}{.}\PY{l+m+mi}{1}\PY{p}{)} \PY{n}{s} \PY{o}{=}
\PY{n}{GaussianRandomWalk}\PY{p}{(}\PY{l+s}{\PYZsq{}}\PY{l+s}{s}\PY{l+s}{\PYZsq{}}\PY{p}{,}
\PY{n}{sigma}\PY{o}{*}\PY{o}{*}\PY{o}{\PYZhy{}}\PY{l+m+mi}{2}\PY{p}{,}
\PY{n}{shape}\PY{o}{=}\PY{n}{n}\PY{p}{)} \PY{n}{volatility\PYZus{}process}
\PY{o}{=}
\PY{n}{Deterministic}\PY{p}{(}\PY{l+s}{\PYZsq{}}\PY{l+s}{volatility\PYZus{}process}\PY{l+s}{\PYZsq{}}\PY{p}{,}
\PY{n}{exp}\PY{p}{(}\PY{o}{\PYZhy{}}\PY{l+m+mi}{2}\PY{o}{*}\PY{n}{s}\PY{p}{)}\PY{p}{)}
\PY{n}{r} \PY{o}{=}
\PY{n}{T}\PY{p}{(}\PY{l+s}{\PYZsq{}}\PY{l+s}{r}\PY{l+s}{\PYZsq{}}\PY{p}{,}
\PY{n}{nu}\PY{p}{,}
\PY{n}{lam}\PY{o}{=}\PY{l+m+mi}{1}\PY{o}{/}\PY{n}{volatility\PYZus{}process}\PY{p}{,}
\PY{n}{observed}\PY{o}{=}\PY{n}{returns}\PY{p}{)}

Notice that we transform the log volatility process s into the volatility
process by exp(-2*s). Here, exp is a Theano function, rather than the
corresponding function in NumPy; Theano provides a large subset of the
mathematical functions that NumPy does.

Also note that we have declared the Model name sp500_model in the first
occurrence of the context manager, rather than splitting it into two lines, as
we did for the first example.

###  4.4 Fitting

Before we draw samples from the posterior, it is prudent to find a decent
starting valuwa by finding a point of relatively high probability. For this
model, the full _maximum a posteriori_ (MAP) point over all variables is
degenerate and has infinite density. But, if we fix log_sigma and nu it is no
longer degenerate, so we find the MAP with respect only to the volatility
process s keeping log_sigma and nu constant at their default values (remember
that we set testval=.1 for sigma). We use the Limited-memory BFGS (L-BFGS)
optimizer, which is provided by the scipy.optimize package, as it is more
efficient for high dimensional functions and we have 400 stochastic random
variables (mostly from s).

To do the sampling, we do a short initial run to put us in a volume of high
probability, then start again at the new starting point. trace[-1] gives us
the last point in the sampling trace. NUTS will recalculate the scaling
parameters based on the new point, and in this case it leads to faster
sampling due to better scaling.

,commandchars=
{},,

{\color{incolor}In [{\color{incolor}14}]:} \PY{k+kn}{import} \PY{n+nn}{scipy}
\PY{k}{with} \PY{n}{sp500\PYZus{}model}\PY{p}{:} \PY{n}{start} \PY{o}{=}
\PY{n}{find\PYZus{}MAP}\PY{p}{(}\PY{n+nb}{vars}\PY{o}{=}\PY{p}{[}\PY{n}{s}\PY{p}{]}\PY{p}{,}
\PY{n}{fmin}\PY{o}{=}\PY{n}{scipy}\PY{o}{.}\PY{n}{optimize}\PY{o}{.}\PY{n}{fmin\PYZus{}l\PYZus{}bfgs\PYZus{}b}\PY{p}{)}
\PY{n}{step} \PY{o}{=}
\PY{n}{NUTS}\PY{p}{(}\PY{n}{scaling}\PY{o}{=}\PY{n}{start}\PY{p}{)}
\PY{n}{trace} \PY{o}{=} \PY{n}{sample}\PY{p}{(}\PY{l+m+mi}{50}\PY{p}{,}
\PY{n}{step}\PY{p}{,} \PY{n}{progressbar}\PY{o}{=}\PY{n+nb+bp}{False}\PY{p}{)}
\PY{c}{\PYZsh{} Start next run at the last sampled position.} \PY{n}{step}
\PY{o}{=}
\PY{n}{NUTS}\PY{p}{(}\PY{n}{scaling}\PY{o}{=}\PY{n}{trace}\PY{p}{[}\PY{o}{\PYZhy{}}\PY{l+m+mi}{1}\PY{p}{]}\PY{p}{,}
\PY{n}{gamma}\PY{o}{=}\PY{o}{.}\PY{l+m+mi}{25}\PY{p}{)} \PY{n}{trace}
\PY{o}{=} \PY{n}{sample}\PY{p}{(}\PY{l+m+mi}{2000}\PY{p}{,}
\PY{n}{step}\PY{p}{,}
\PY{n}{start}\PY{o}{=}\PY{n}{trace}\PY{p}{[}\PY{o}{\PYZhy{}}\PY{l+m+mi}{1}\PY{p}{]}\PY{p}{,}\PY{n}{progressbar}\PY{o}{=}\PY{n+nb+bp}{False}\PY{p}{,}
\PY{n}{njobs}\PY{o}{=}\PY{l+m+mi}{4}\PY{p}{)}

We can check our samples by looking at the traceplot for nu and sigma.

,commandchars=
{},,

{\color{incolor}In [{\color{incolor}16}]:} \PY{c}{\PYZsh{}figsize(12,6)}
\PY{n}{traceplot}\PY{p}{(}\PY{n}{trace}\PY{p}{,} \PY{p}{[}\PY{n}{nu}\PY{p}{,}
\PY{n}{sigma}\PY{p}{]}\PY{p}{)}\PY{p}{;}

Finally we plot the distribution of volatility paths by plotting many of our
sampled volatility paths on the same graph. Each is rendered partially
transparent (via the alpha argument in Matplotlib’s plot function) so the
regions where many paths overlap are shaded more darkly.

,commandchars=
{},,

{\color{incolor}In [{\color{incolor}16}]:}
\PY{n}{figsize}\PY{p}{(}\PY{l+m+mi}{15}\PY{p}{,}\PY{l+m+mi}{8}\PY{p}{)}
\PY{n}{pl}\PY{o}{.}\PY{n}{title}\PY{p}{(}\PY{l+s}{\PYZdq{}}\PY{l+s}{volatility\PYZus{}process}\PY{l+s}{\PYZdq{}}\PY{p}{)}\PY{p}{;}
\PY{n}{pl}\PY{o}{.}\PY{n}{plot}\PY{p}{(}\PY{n}{trace}\PY{p}{[}\PY{l+s}{\PYZsq{}}\PY{l+s}{volatility\PYZus{}process}\PY{l+s}{\PYZsq{}}\PY{p}{,}\PY{p}{:}\PY{p}{:}\PY{l+m+mi}{30}\PY{p}{]}\PY{o}{.}\PY{n}{T}\PY{p}{,}\PY{l+s}{\PYZsq{}}\PY{l+s}{b}\PY{l+s}{\PYZsq{}}\PY{p}{,}
\PY{n}{alpha}\PY{o}{=}\PY{o}{.}\PY{l+m+mo}{03}\PY{p}{)}\PY{p}{;}
\PY{n}{pl}\PY{o}{.}\PY{n}{ylim}\PY{p}{(}\PY{l+m+mi}{0}\PY{p}{,}
\PY{o}{.}\PY{l+m+mo}{001}\PY{p}{)}
\PY{n}{pl}\PY{o}{.}\PY{n}{xlabel}\PY{p}{(}\PY{l+s}{\PYZsq{}}\PY{l+s}{time}\PY{l+s}{\PYZsq{}}\PY{p}{)}\PY{p}{;}
\PY{n}{pl}\PY{o}{.}\PY{n}{ylabel}\PY{p}{(}\PY{l+s}{\PYZsq{}}\PY{l+s}{log
volatility}\PY{l+s}{\PYZsq{}}\PY{p}{)}\PY{p}{;}
\PY{n}{pl}\PY{o}{.}\PY{n}{plot}\PY{p}{(}\PY{n}{returns}\PY{o}{.}\PY{n}{cumsum}\PY{p}{(}\PY{p}{)}\PY{o}{*}\PY{o}{.}\PY{l+m+mo}{0025}\PY{o}{+}\PY{o}{.}\PY{l+m+mo}{0005}\PY{p}{)}\PY{p}{;}

##  5 Case study 2: Coal mining disasters

Consider the following time series of recorded coal mining disasters in the UK
from 1851 to 1962 (Jarrett, 1979). The number of disasters is thought to have
been affected by changes in safety regulations during this period.
Unfortunately, we also have pair of years with missing data, identified as
missing by a NumPy MaskedArray using -999 as the marker value.

Next we will build a model for this series and attempt to estimate when the
change occured. At the same time, we will see how to handle missing data, use
multiple samplers and sample from discrete random variables.

,commandchars=
{},,

{\color{incolor}In [{\color{incolor}17}]:} \PY{n}{disaster\PYZus{}data}
\PY{o}{=}
\PY{n}{np}\PY{o}{.}\PY{n}{ma}\PY{o}{.}\PY{n}{masked\PYZus{}values}\PY{p}{(}\PY{p}{[}\PY{l+m+mi}{4}\PY{p}{,}
\PY{l+m+mi}{5}\PY{p}{,} \PY{l+m+mi}{4}\PY{p}{,} \PY{l+m+mi}{0}\PY{p}{,}
\PY{l+m+mi}{1}\PY{p}{,} \PY{l+m+mi}{4}\PY{p}{,} \PY{l+m+mi}{3}\PY{p}{,}
\PY{l+m+mi}{4}\PY{p}{,} \PY{l+m+mi}{0}\PY{p}{,} \PY{l+m+mi}{6}\PY{p}{,}
\PY{l+m+mi}{3}\PY{p}{,} \PY{l+m+mi}{3}\PY{p}{,} \PY{l+m+mi}{4}\PY{p}{,}
\PY{l+m+mi}{0}\PY{p}{,} \PY{l+m+mi}{2}\PY{p}{,} \PY{l+m+mi}{6}\PY{p}{,}
\PY{l+m+mi}{3}\PY{p}{,} \PY{l+m+mi}{3}\PY{p}{,} \PY{l+m+mi}{5}\PY{p}{,}
\PY{l+m+mi}{4}\PY{p}{,} \PY{l+m+mi}{5}\PY{p}{,} \PY{l+m+mi}{3}\PY{p}{,}
\PY{l+m+mi}{1}\PY{p}{,} \PY{l+m+mi}{4}\PY{p}{,} \PY{l+m+mi}{4}\PY{p}{,}
\PY{l+m+mi}{1}\PY{p}{,} \PY{l+m+mi}{5}\PY{p}{,} \PY{l+m+mi}{5}\PY{p}{,}
\PY{l+m+mi}{3}\PY{p}{,} \PY{l+m+mi}{4}\PY{p}{,} \PY{l+m+mi}{2}\PY{p}{,}
\PY{l+m+mi}{5}\PY{p}{,} \PY{l+m+mi}{2}\PY{p}{,} \PY{l+m+mi}{2}\PY{p}{,}
\PY{l+m+mi}{3}\PY{p}{,} \PY{l+m+mi}{4}\PY{p}{,} \PY{l+m+mi}{2}\PY{p}{,}
\PY{l+m+mi}{1}\PY{p}{,} \PY{l+m+mi}{3}\PY{p}{,}
\PY{o}{\PYZhy{}}\PY{l+m+mi}{999}\PY{p}{,} \PY{l+m+mi}{2}\PY{p}{,}
\PY{l+m+mi}{1}\PY{p}{,} \PY{l+m+mi}{1}\PY{p}{,} \PY{l+m+mi}{1}\PY{p}{,}
\PY{l+m+mi}{1}\PY{p}{,} \PY{l+m+mi}{3}\PY{p}{,} \PY{l+m+mi}{0}\PY{p}{,}
\PY{l+m+mi}{0}\PY{p}{,} \PY{l+m+mi}{1}\PY{p}{,} \PY{l+m+mi}{0}\PY{p}{,}
\PY{l+m+mi}{1}\PY{p}{,} \PY{l+m+mi}{1}\PY{p}{,} \PY{l+m+mi}{0}\PY{p}{,}
\PY{l+m+mi}{0}\PY{p}{,} \PY{l+m+mi}{3}\PY{p}{,} \PY{l+m+mi}{1}\PY{p}{,}
\PY{l+m+mi}{0}\PY{p}{,} \PY{l+m+mi}{3}\PY{p}{,} \PY{l+m+mi}{2}\PY{p}{,}
\PY{l+m+mi}{2}\PY{p}{,} \PY{l+m+mi}{0}\PY{p}{,} \PY{l+m+mi}{1}\PY{p}{,}
\PY{l+m+mi}{1}\PY{p}{,} \PY{l+m+mi}{1}\PY{p}{,} \PY{l+m+mi}{0}\PY{p}{,}
\PY{l+m+mi}{1}\PY{p}{,} \PY{l+m+mi}{0}\PY{p}{,} \PY{l+m+mi}{1}\PY{p}{,}
\PY{l+m+mi}{0}\PY{p}{,} \PY{l+m+mi}{0}\PY{p}{,} \PY{l+m+mi}{0}\PY{p}{,}
\PY{l+m+mi}{2}\PY{p}{,} \PY{l+m+mi}{1}\PY{p}{,} \PY{l+m+mi}{0}\PY{p}{,}
\PY{l+m+mi}{0}\PY{p}{,} \PY{l+m+mi}{0}\PY{p}{,} \PY{l+m+mi}{1}\PY{p}{,}
\PY{l+m+mi}{1}\PY{p}{,} \PY{l+m+mi}{0}\PY{p}{,} \PY{l+m+mi}{2}\PY{p}{,}
\PY{l+m+mi}{3}\PY{p}{,} \PY{l+m+mi}{3}\PY{p}{,} \PY{l+m+mi}{1}\PY{p}{,}
\PY{o}{\PYZhy{}}\PY{l+m+mi}{999}\PY{p}{,} \PY{l+m+mi}{2}\PY{p}{,}
\PY{l+m+mi}{1}\PY{p}{,} \PY{l+m+mi}{1}\PY{p}{,} \PY{l+m+mi}{1}\PY{p}{,}
\PY{l+m+mi}{1}\PY{p}{,} \PY{l+m+mi}{2}\PY{p}{,} \PY{l+m+mi}{4}\PY{p}{,}
\PY{l+m+mi}{2}\PY{p}{,} \PY{l+m+mi}{0}\PY{p}{,} \PY{l+m+mi}{0}\PY{p}{,}
\PY{l+m+mi}{1}\PY{p}{,} \PY{l+m+mi}{4}\PY{p}{,} \PY{l+m+mi}{0}\PY{p}{,}
\PY{l+m+mi}{0}\PY{p}{,} \PY{l+m+mi}{0}\PY{p}{,} \PY{l+m+mi}{1}\PY{p}{,}
\PY{l+m+mi}{0}\PY{p}{,} \PY{l+m+mi}{0}\PY{p}{,} \PY{l+m+mi}{0}\PY{p}{,}
\PY{l+m+mi}{0}\PY{p}{,} \PY{l+m+mi}{0}\PY{p}{,} \PY{l+m+mi}{1}\PY{p}{,}
\PY{l+m+mi}{0}\PY{p}{,} \PY{l+m+mi}{0}\PY{p}{,} \PY{l+m+mi}{1}\PY{p}{,}
\PY{l+m+mi}{0}\PY{p}{,} \PY{l+m+mi}{1}\PY{p}{]}\PY{p}{,}
\PY{n}{value}\PY{o}{=}\PY{o}{\PYZhy{}}\PY{l+m+mi}{999}\PY{p}{)} \PY{n}{year}
\PY{o}{=} \PY{n}{np}\PY{o}{.}\PY{n}{arange}\PY{p}{(}\PY{l+m+mi}{1851}\PY{p}{,}
\PY{l+m+mi}{1962}\PY{p}{)} \PY{n}{plot}\PY{p}{(}\PY{n}{year}\PY{p}{,}
\PY{n}{disaster\PYZus{}data}\PY{p}{,}
\PY{l+s}{\PYZsq{}}\PY{l+s}{o}\PY{l+s}{\PYZsq{}}\PY{p}{,}
\PY{n}{markersize}\PY{o}{=}\PY{l+m+mi}{8}\PY{p}{)}\PY{p}{;}
\PY{n}{ylabel}\PY{p}{(}\PY{l+s}{\PYZdq{}}\PY{l+s}{Disaster
count}\PY{l+s}{\PYZdq{}}\PY{p}{)}
\PY{n}{xlabel}\PY{p}{(}\PY{l+s}{\PYZdq{}}\PY{l+s}{Year}\PY{l+s}{\PYZdq{}}\PY{p}{)}

,commandchars=
{},,

{\color{outcolor}Out[{\color{outcolor}17}]:} <matplotlib.text.Text at
0x1129c3650>

Occurrences of disasters in the time series is thought to follow a Poisson
process with a large rate parameter in the early part of the time series, and
from one with a smaller rate in the later part. We are interested in locating
the change point in the series, which perhaps is related to changes in mining
safety regulations.

In our model,

$\displaystyle D_{t}$
$\displaystyle\sim\text{Pois}(r_{t}),r_{t}=\begin{cases}l,&\text{if }t<s\\\
r,&\text{if }t\geq s\end{cases}$ $\displaystyle s$
$\displaystyle\sim\text{Unif}(t_{l},t_{h})$ $\displaystyle e$
$\displaystyle\sim\text{exp}(1)$ $\displaystyle l$
$\displaystyle\sim\text{exp}(1)$

the parameters are defined as follows: * $D_{t}$: The number of disasters in
year $t$ * $r_{t}$: The rate parameter of the Poisson distribution of
disasters in year $t$. * $s$: The year in which the rate parameter changes
(the switchpoint). * $e$: The rate parameter before the switchpoint $s$. *
$l$: The rate parameter after the switchpoint $s$. * $t_{l}$, $t_{h}$: The
lower and upper boundaries of year $t$.

This model is built much like our previous models. The major differences are
the introduction of discrete variables with the Poisson and discrete-uniform
priors and the novel form of the deterministic random variable rate.

,commandchars=
{},,

{\color{incolor}In [{\color{incolor}18}]:} \PY{k+kn}{from} \PY{n+nn}{pymc3}
\PY{k+kn}{import} \PY{n}{DiscreteUniform}\PY{p}{,} \PY{n}{Poisson}\PY{p}{,}
\PY{n}{switch} \PY{k}{with} \PY{n}{Model}\PY{p}{(}\PY{p}{)} \PY{k}{as}
\PY{n}{disaster\PYZus{}model}\PY{p}{:} \PY{n}{switchpoint} \PY{o}{=}
\PY{n}{DiscreteUniform}\PY{p}{(}\PY{l+s}{\PYZsq{}}\PY{l+s}{switchpoint}\PY{l+s}{\PYZsq{}}\PY{p}{,}
\PY{n}{lower}\PY{o}{=}\PY{n}{year}\PY{o}{.}\PY{n}{min}\PY{p}{(}\PY{p}{)}\PY{p}{,}
\PY{n}{upper}\PY{o}{=}\PY{n}{year}\PY{o}{.}\PY{n}{max}\PY{p}{(}\PY{p}{)}\PY{p}{,}
\PY{n}{testval}\PY{o}{=}\PY{l+m+mi}{1900}\PY{p}{)} \PY{c}{\PYZsh{} Priors for
pre\PYZhy{} and post\PYZhy{}switch rates number of disasters}
\PY{n}{early\PYZus{}rate} \PY{o}{=}
\PY{n}{Exponential}\PY{p}{(}\PY{l+s}{\PYZsq{}}\PY{l+s}{early\PYZus{}rate}\PY{l+s}{\PYZsq{}}\PY{p}{,}
\PY{l+m+mi}{1}\PY{p}{)} \PY{n}{late\PYZus{}rate} \PY{o}{=}
\PY{n}{Exponential}\PY{p}{(}\PY{l+s}{\PYZsq{}}\PY{l+s}{late\PYZus{}rate}\PY{l+s}{\PYZsq{}}\PY{p}{,}
\PY{l+m+mi}{1}\PY{p}{)} \PY{c}{\PYZsh{} Allocate appropriate Poisson rates to
years before and after current} \PY{n}{rate} \PY{o}{=}
\PY{n}{switch}\PY{p}{(}\PY{n}{switchpoint} \PY{o}{\PYZgt{}}\PY{o}{=}
\PY{n}{year}\PY{p}{,} \PY{n}{early\PYZus{}rate}\PY{p}{,}
\PY{n}{late\PYZus{}rate}\PY{p}{)} \PY{n}{disasters} \PY{o}{=}
\PY{n}{Poisson}\PY{p}{(}\PY{l+s}{\PYZsq{}}\PY{l+s}{disasters}\PY{l+s}{\PYZsq{}}\PY{p}{,}
\PY{n}{rate}\PY{p}{,}
\PY{n}{observed}\PY{o}{=}\PY{n}{disaster\PYZus{}data}\PY{p}{)}

The logic for the rate random variable,

commandchars=
{},,,

\NormalTok{rate = switch(switchpoint >= year, early_rate, late_rate)}

is implemented using switch, a Theano function that works like an if
statement. It uses the first argument to switch between the next two
arguments.

Missing values are handled transparently by passing a MaskedArray or a
pandas.DataFrame with NaN values to the observed argument when creating an
observed stochastic random variable. Behind the scenes, another random
variable, disasters.missing_values is created to model the missing values. All
we need to do to handle the missing values is ensure we sample this random
variable as well.

Unfortunately because they are discrete variables and thus have no meaningful
gradient, we cannot use NUTS for sampling switchpoint or the missing disaster
observations. Instead, we will sample using a Metroplis step method, which
implements adaptive Metropolis-Hastings, because it is designed to handle
discrete values.

We sample with both samplers at once by passing them to the sample function in
a list. Each new sample is generated by first applying step1 then step2.

,commandchars=
{},,

{\color{incolor}In [{\color{incolor}19}]:} \PY{k+kn}{from} \PY{n+nn}{pymc3}
\PY{k+kn}{import} \PY{n}{Metropolis} \PY{k}{with}
\PY{n}{disaster\PYZus{}model}\PY{p}{:} \PY{n}{step1} \PY{o}{=}
\PY{n}{NUTS}\PY{p}{(}\PY{p}{[}\PY{n}{early\PYZus{}rate}\PY{p}{,}
\PY{n}{late\PYZus{}rate}\PY{p}{]}\PY{p}{)} \PY{c}{\PYZsh{} Use Metropolis for
switchpoint, and missing values since it accomodates discrete variables}
\PY{n}{step2} \PY{o}{=}
\PY{n}{Metropolis}\PY{p}{(}\PY{p}{[}\PY{n}{switchpoint}\PY{p}{,}
\PY{n}{disasters}\PY{o}{.}\PY{n}{missing\PYZus{}values}\PY{p}{[}\PY{l+m+mi}{0}\PY{p}{]}\PY{p}{]}
\PY{p}{)} \PY{n}{trace} \PY{o}{=}
\PY{n}{sample}\PY{p}{(}\PY{l+m+mi}{10000}\PY{p}{,}
\PY{n}{step}\PY{o}{=}\PY{p}{[}\PY{n}{step1}\PY{p}{,}
\PY{n}{step2}\PY{p}{]}\PY{p}{)}

,commandchars=
{},,

[-----------------100\%-----------------] 10000 of 10000 complete in 10.3 sec

In the trace plot below we can see that there’s about a 10 year span that’s
plausible for a significant change in safety, but a 5 year span that contains
most of the probability mass. The distribution is jagged because of the jumpy
relationship beween the year switchpoint and the likelihood and not due to
sampling error.

,commandchars=
{},,

{\color{incolor}In [{\color{incolor}21}]:}
\PY{n}{traceplot}\PY{p}{(}\PY{n}{trace}\PY{p}{)}\PY{p}{;}

##  6 Arbitrary deterministics

Due to its reliance on Theano, PyMC3 provides many mathematical functions and
operators for transforming random variables into new random variables.
However, the library of functions in Theano is not exhaustive, therefore
Theano and PyMC3 provide functionality for creating arbitrary Theano functions
in pure Python, and including these functions in PyMC models. This is
supported with the as_op function decorator.

Theano needs to know the types of the inputs and outputs of a function, which
are specified for as_op by itypes for inputs and otypes for outputs. The
Theano documentation includes an overview of the available types.

,commandchars=
{},,

{\color{incolor}In [{\color{incolor}22}]:} \PY{k+kn}{import}
\PY{n+nn}{theano.tensor} \PY{k+kn}{as} \PY{n+nn}{T} \PY{k+kn}{from}
\PY{n+nn}{theano.compile.ops} \PY{k+kn}{import} \PY{n}{as\PYZus{}op}
\PY{n+nd}{@as\PYZus{}op}\PY{p}{(}\PY{n}{itypes}\PY{o}{=}\PY{p}{[}\PY{n}{T}\PY{o}{.}\PY{n}{lscalar}\PY{p}{]}\PY{p}{,}
\PY{n}{otypes}\PY{o}{=}\PY{p}{[}\PY{n}{T}\PY{o}{.}\PY{n}{lscalar}\PY{p}{]}\PY{p}{)}
\PY{k}{def}
\PY{n+nf}{crazy\PYZus{}modulo3}\PY{p}{(}\PY{n}{value}\PY{p}{)}\PY{p}{:}
\PY{k}{if} \PY{n}{value} \PY{o}{\PYZgt{}} \PY{l+m+mi}{0}\PY{p}{:}
\PY{k}{return} \PY{n}{value} \PY{o}{\PYZpc{}} \PY{l+m+mi}{3} \PY{k}{else}
\PY{p}{:} \PY{k}{return} \PY{p}{(}\PY{o}{\PYZhy{}}\PY{n}{value} \PY{o}{+}
\PY{l+m+mi}{1}\PY{p}{)} \PY{o}{\PYZpc{}} \PY{l+m+mi}{3} \PY{k}{with}
\PY{n}{Model}\PY{p}{(}\PY{p}{)} \PY{k}{as}
\PY{n}{model\PYZus{}deterministic}\PY{p}{:} \PY{n}{a} \PY{o}{=}
\PY{n}{Poisson}\PY{p}{(}\PY{l+s}{\PYZsq{}}\PY{l+s}{a}\PY{l+s}{\PYZsq{}}\PY{p}{,}
\PY{l+m+mi}{1}\PY{p}{)} \PY{n}{b} \PY{o}{=}
\PY{n}{crazy\PYZus{}modulo3}\PY{p}{(}\PY{n}{a}\PY{p}{)}

An important drawback of this approach is that it is not possible for theano
to inspect these functions in order to compute the gradient required for the
Hamiltonian-based samplers. Therefore, it is not possible to use the HMC or
NUTS samplers for a model that uses such an operator. However, it is possible
to add a gradient if we inherit from theano.Op instead of using as_op. The
PyMC example set includes a more elaborate example of the usage of as_op.

##  7 Arbitrary distributions

Similarly, the library of statistical distributions in PyMC3 is not
exhaustive, but PyMC allows for the creation of user-defined functions for an
arbitrary probability distribution. For simple statistical distributions, the
DensityDist function takes as an argument any function that calculates a log-
probability $log(p(x))$. This function may employ other random variables in
its calculation. Here is an example inspired by a blog post by Jake Vanderplas
on which priors to use for a linear regression (Vanderplas, 2014).

commandchars=
{},,,

\CharTok{import} \NormalTok{theano.tensor }\CharTok{as}
\NormalTok{T}\CharTok{from} \NormalTok{pymc3 }\CharTok{import}
\NormalTok{DensityDist, Uniform}\KeywordTok{with} \NormalTok{Model()
}\CharTok{as} \NormalTok{model:} \NormalTok{alpha =
Uniform(}\StringTok{’intercept’}\NormalTok{, -}\DecValTok{100}\NormalTok{,
}\DecValTok{100}\NormalTok{)} \CommentTok{# Create custom densities}
\NormalTok{beta = DensityDist(}\StringTok{’beta’}\NormalTok{,
}\KeywordTok{lambda} \NormalTok{value: -}\FloatTok{1.5} \NormalTok{*
T.log(}\DecValTok{1} \NormalTok{+ value**}\DecValTok{2}\NormalTok{),
testval=}\DecValTok{0}\NormalTok{)} \NormalTok{eps =
DensityDist(}\StringTok{’eps’}\NormalTok{, }\KeywordTok{lambda}
\NormalTok{value: -T.log(T.abs_(value)), testval=}\DecValTok{1}\NormalTok{)}
\CommentTok{# Create likelihood} \NormalTok{like =
Normal(}\StringTok{’y_est’}\NormalTok{, mu=alpha + beta * X, sd=eps,
observed=Y)}

For more complex distributions, one can create a subclass of Continuous or
Discrete and provide the custom logp function, as required. This is how the
built-in distributions in PyMC are specified. As an example, fields like
psychology and astrophysics have complex likelihood functions for a particular
process that may require numerical approximation. In these cases, it is
impossible to write the function in terms of predefined theano operators and
we must use a custom theano operator using as_op or inheriting from theano.Op.

Implementing the beta variable above as a Continuous subclass is shown below,
along with a sub-function using the as_op decorator, though this is not
strictly necessary.

,commandchars=
{},,

{\color{incolor}In [{\color{incolor}23}]:} \PY{k+kn}{from}
\PY{n+nn}{pymc3.distributions} \PY{k+kn}{import} \PY{n}{Continuous}
\PY{k}{class} \PY{n+nc}{Beta}\PY{p}{(}\PY{n}{Continuous}\PY{p}{)}\PY{p}{:}
\PY{k}{def}
\PY{n+nf}{\PYZus{}\PYZus{}init\PYZus{}\PYZus{}}\PY{p}{(}\PY{n+nb+bp}{self}\PY{p}{,}
\PY{n}{mu}\PY{p}{,} \PY{o}{*}\PY{n}{args}\PY{p}{,}
\PY{o}{*}\PY{o}{*}\PY{n}{kwargs}\PY{p}{)}\PY{p}{:}
\PY{n+nb}{super}\PY{p}{(}\PY{n}{Beta}\PY{p}{,}
\PY{n+nb+bp}{self}\PY{p}{)}\PY{o}{.}\PY{n}{\PYZus{}\PYZus{}init\PYZus{}\PYZus{}}\PY{p}{(}\PY{o}{*}\PY{n}{args}\PY{p}{,}
\PY{o}{*}\PY{o}{*}\PY{n}{kwargs}\PY{p}{)}
\PY{n+nb+bp}{self}\PY{o}{.}\PY{n}{mu} \PY{o}{=} \PY{n}{mu}
\PY{n+nb+bp}{self}\PY{o}{.}\PY{n}{mode} \PY{o}{=} \PY{n}{mu} \PY{k}{def}
\PY{n+nf}{logp}\PY{p}{(}\PY{n+nb+bp}{self}\PY{p}{,}
\PY{n}{value}\PY{p}{)}\PY{p}{:} \PY{n}{mu} \PY{o}{=}
\PY{n+nb+bp}{self}\PY{o}{.}\PY{n}{mu} \PY{k}{return}
\PY{n}{beta\PYZus{}logp}\PY{p}{(}\PY{n}{value} \PY{o}{\PYZhy{}}
\PY{n}{mu}\PY{p}{)}
\PY{n+nd}{@as\PYZus{}op}\PY{p}{(}\PY{n}{itypes}\PY{o}{=}\PY{p}{[}\PY{n}{T}\PY{o}{.}\PY{n}{dscalar}\PY{p}{]}\PY{p}{,}
\PY{n}{otypes}\PY{o}{=}\PY{p}{[}\PY{n}{T}\PY{o}{.}\PY{n}{dscalar}\PY{p}{]}\PY{p}{)}
\PY{k}{def}
\PY{n+nf}{beta\PYZus{}logp}\PY{p}{(}\PY{n}{value}\PY{p}{)}\PY{p}{:}
\PY{k}{return} \PY{o}{\PYZhy{}}\PY{l+m+mf}{1.5} \PY{o}{*}
\PY{n}{np}\PY{o}{.}\PY{n}{log}\PY{p}{(}\PY{l+m+mi}{1} \PY{o}{+}
\PY{p}{(}\PY{n}{value}\PY{p}{)}\PY{o}{*}\PY{o}{*}\PY{l+m+mi}{2}\PY{p}{)}
\PY{k}{with} \PY{n}{Model}\PY{p}{(}\PY{p}{)} \PY{k}{as} \PY{n}{model}\PY{p}{:}
\PY{n}{beta} \PY{o}{=}
\PY{n}{Beta}\PY{p}{(}\PY{l+s}{\PYZsq{}}\PY{l+s}{slope}\PY{l+s}{\PYZsq{}}\PY{p}{,}
\PY{n}{mu}\PY{o}{=}\PY{l+m+mi}{0}\PY{p}{,}
\PY{n}{testval}\PY{o}{=}\PY{l+m+mi}{0}\PY{p}{)}

##  8 Generalized Linear Models

Generalized Linear Models (GLMs) are a class of flexible models that are
widely used to estimate regression relationships between a single outcome
variable and one or multiple predictors. Because these models are so common,
PyMC3 offers a glm submodule that allows flexible creation of various GLMs
with an intuitive R-like syntax that is implemented via the patsy module.

The glm submodule requires data to be included as a pandas DataFrame. Hence,
for our linear regression example:

,commandchars=
{},,

{\color{incolor}In [{\color{incolor}24}]:} \PY{c}{\PYZsh{} Convert X and Y to
a pandas DataFrame} \PY{k+kn}{import} \PY{n+nn}{pandas} \PY{n}{df} \PY{o}{=}
\PY{n}{pandas}\PY{o}{.}\PY{n}{DataFrame}\PY{p}{(}\PY{p}{\PYZob{}}\PY{l+s}{\PYZsq{}}\PY{l+s}{x1}\PY{l+s}{\PYZsq{}}\PY{p}{:}
\PY{n}{X1}\PY{p}{,} \PY{l+s}{\PYZsq{}}\PY{l+s}{x2}\PY{l+s}{\PYZsq{}}\PY{p}{:}
\PY{n}{X2}\PY{p}{,} \PY{l+s}{\PYZsq{}}\PY{l+s}{y}\PY{l+s}{\PYZsq{}}\PY{p}{:}
\PY{n}{Y}\PY{p}{\PYZcb{}}\PY{p}{)}

The model can then be very concisely specified in one line of code.

,commandchars=
{},,

{\color{incolor}In [{\color{incolor}25}]:} \PY{k+kn}{from}
\PY{n+nn}{pymc3.glm} \PY{k+kn}{import} \PY{n}{glm} \PY{k}{with}
\PY{n}{Model}\PY{p}{(}\PY{p}{)} \PY{k}{as} \PY{n}{model\PYZus{}glm}\PY{p}{:}
\PY{n}{glm}\PY{p}{(}\PY{l+s}{\PYZsq{}}\PY{l+s}{y \PYZti{} x1 +
x2}\PY{l+s}{\PYZsq{}}\PY{p}{,} \PY{n}{df}\PY{p}{)}

The error distribution, if not specified via the family argument, is assumed
to be normal. In the case of logistic regression, this can be modified by
passing in a Binomial family object.

,commandchars=
{},,

{\color{incolor}In [{\color{incolor}26}]:} \PY{k+kn}{from}
\PY{n+nn}{pymc3.glm.families} \PY{k+kn}{import} \PY{n}{Binomial}
\PY{n}{df\PYZus{}logistic} \PY{o}{=}
\PY{n}{pandas}\PY{o}{.}\PY{n}{DataFrame}\PY{p}{(}\PY{p}{\PYZob{}}\PY{l+s}{\PYZsq{}}\PY{l+s}{x1}\PY{l+s}{\PYZsq{}}\PY{p}{:}
\PY{n}{X1}\PY{p}{,} \PY{l+s}{\PYZsq{}}\PY{l+s}{x2}\PY{l+s}{\PYZsq{}}\PY{p}{:}
\PY{n}{X2}\PY{p}{,} \PY{l+s}{\PYZsq{}}\PY{l+s}{y}\PY{l+s}{\PYZsq{}}\PY{p}{:}
\PY{n}{Y} \PY{o}{\PYZgt{}} \PY{l+m+mi}{0}\PY{p}{\PYZcb{}}\PY{p}{)}
\PY{k}{with} \PY{n}{Model}\PY{p}{(}\PY{p}{)} \PY{k}{as}
\PY{n}{model\PYZus{}glm\PYZus{}logistic}\PY{p}{:}
\PY{n}{glm}\PY{p}{(}\PY{l+s}{\PYZsq{}}\PY{l+s}{y \PYZti{} x1 +
x2}\PY{l+s}{\PYZsq{}}\PY{p}{,} \PY{n}{df\PYZus{}logistic}\PY{p}{,}
\PY{n}{family}\PY{o}{=}\PY{n}{Binomial}\PY{p}{(}\PY{p}{)}\PY{p}{)}

##  9 Backends

PyMC3 has support for different ways to store samples during and after
sampling, called backends, including in-memory (default), text file, and
SQLite. These can be found in pymc.backends:

By default, an in-memory ndarray is used but if the samples would get too
large to be held in memory we could use the sqlite backend:

,commandchars=
{},,

{\color{incolor}In [{\color{incolor}27}]:} \PY{k+kn}{from}
\PY{n+nn}{pymc3.backends} \PY{k+kn}{import} \PY{n}{SQLite} \PY{k}{with}
\PY{n}{model\PYZus{}glm\PYZus{}logistic}\PY{p}{:} \PY{n}{backend} \PY{o}{=}
\PY{n}{SQLite}\PY{p}{(}\PY{l+s}{\PYZsq{}}\PY{l+s}{trace.sqlite}\PY{l+s}{\PYZsq{}}\PY{p}{)}
\PY{n}{trace} \PY{o}{=} \PY{n}{sample}\PY{p}{(}\PY{l+m+mi}{5000}\PY{p}{,}
\PY{n}{Metropolis}\PY{p}{(}\PY{p}{)}\PY{p}{,}
\PY{n}{trace}\PY{o}{=}\PY{n}{backend}\PY{p}{)}

,commandchars=
{},,

[-----------------100\%-----------------] 5000 of 5000 complete in 1.7 sec

,commandchars=
{},,

{\color{incolor}In [{\color{incolor}28}]:}
\PY{n}{summary}\PY{p}{(}\PY{n}{trace}\PY{p}{,}
\PY{n+nb}{vars}\PY{o}{=}\PY{p}{[}\PY{l+s}{\PYZsq{}}\PY{l+s}{x1}\PY{l+s}{\PYZsq{}}\PY{p}{,}
\PY{l+s}{\PYZsq{}}\PY{l+s}{x2}\PY{l+s}{\PYZsq{}}\PY{p}{]}\PY{p}{)}

,commandchars=
{},,

x1: Mean SD MC Error 95\% HPD interval
\------------------------------------------------------------------- 0.141
5.371 0.522 [-10.872, 11.579] Posterior quantiles: 2.5 25 50 75 97.5
|--------------|==============|==============|--------------| -11.702 -2.722
-0.043 2.959 11.026x2: Mean SD MC Error 95\% HPD interval
\------------------------------------------------------------------- -0.627
26.921 2.620 [-56.304, 57.165] Posterior quantiles: 2.5 25 50 75 97.5
|--------------|==============|==============|--------------| -55.090 -14.627
0.022 13.840 59.016

The stored trace can then later be loaded using the load command:

,commandchars=
{},,

{\color{incolor}In [{\color{incolor}29}]:} \PY{k+kn}{from}
\PY{n+nn}{pymc3.backends.sqlite} \PY{k+kn}{import} \PY{n}{load} \PY{k}{with}
\PY{n}{basic\PYZus{}model}\PY{p}{:} \PY{n}{trace\PYZus{}loaded} \PY{o}{=}
\PY{n}{load}\PY{p}{(}\PY{l+s}{\PYZsq{}}\PY{l+s}{trace.sqlite}\PY{l+s}{\PYZsq{}}\PY{p}{)}

More information about backends can be found in the docstring of
pymc.backends.

##  10 References

Patil, A., D. Huard and C.J. Fonnesbeck. (2010) PyMC: Bayesian Stochastic
Modelling in Python. Journal of Statistical Software, 35(4), pp. 1-81

Bastien, F., Lamblin, P., Pascanu, R., Bergstra, J., Goodfellow, I., Bergeron,
A., Bouchard, N., Warde-Farley, D., and Bengio, Y. (2012) “Theano: new
features and speed improvements”. NIPS 2012 deep learning workshop.

Bergstra, J., Breuleux, O., Bastien, F., Lamblin, P., Pascanu, R., Desjardins,
G., Turian, J., Warde-Farley, D., and Bengio, Y. (2010) “Theano: A CPU and GPU
Math Expression Compiler”. Proceedings of the Python for Scientific Computing
Conference (SciPy) 2010. June 30 - July 3, Austin, TX

Lunn, D.J., Thomas, A., Best, N., and Spiegelhalter, D. (2000) WinBUGS – a
Bayesian modelling framework: concepts, structure, and extensibility.
Statistics and Computing, 10:325–337.

Neal, R.M. Slice sampling. Annals of Statistics. (2003). doi:10.2307/3448413.

van Rossum, G. The Python Library Reference Release 2.6.5., (2010). URL
http://docs.python.org/library/.

Duane, S., Kennedy, A. D., Pendleton, B. J., and Roweth, D. (1987) “Hybrid
Monte Carlo”, Physics Letters, vol. 195, pp. 216-222.

Stan Development Team. (2014). Stan: A C++ Library for Probability and
Sampling, Version 2.5.0. http://mc-stan.org.

Gamerman, D. Markov Chain Monte Carlo: statistical simulation for Bayesian
inference. Chapman and Hall, 1997.

Hoffman, M. D., & Gelman, A. (2014). The No-U-Turn Sampler: Adaptively Setting
Path Lengths in Hamiltonian Monte Carlo. The Journal of Machine Learning
Research, 30.

Vanderplas, Jake. “Frequentism and Bayesianism IV: How to be a Bayesian in
Python.” Pythonic Perambulations. N.p., 14 Jun 2014. Web. 27 May. 2015\.
https://jakevdp.github.io/blog/2014/06/14/frequentism-and-
bayesianism-4-bayesian-in-python/.

R.G. Jarrett. A note on the intervals between coal mining disasters.
Biometrika, 66:191–193, 1979.