{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# CHAPTER 16. Metric-Predicted Variable on One or Two Groups"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* 싸이그래머 / 베이지안R [1]\n",
    "* 김무성"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Contents"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* 16.1. Estimating the Mean and Standard Deviation of a Normal Distribution\n",
    "* 16.2. Outliers and Robust Estimation : The t Distribution\n",
    "* 16.3. Two Groups\n",
    "* 16.4. Other Noise Distributions and Transforming Data\n",
    "* 16.5. EXERCISES"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"figures/tbl15.3.png\" width=600 />"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this chapter, we consider a situation in which we have a <font color=\"red\">metric-predicted variable</font> that is observed for items <font color=\"blue\">from one or two groups</font>."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* For example, \n",
    "    - we could measure \n",
    "        - the <font color=\"red\">blood pressure</font> (i.e., a metric variable) for people \n",
    "        - randomly sampled <font color=\"blue\">from first-year university students</font> (i.e., a single group).\n",
    "    - In this case, we might be interested in \n",
    "        - how much the <font color=\"orange\">group’s typical blood pressure differs</font> \n",
    "        - from the recommended value for people of that age as published by a federal agency."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* As another example, \n",
    "    - we could measure \n",
    "        - the <font color=\"red\">IQ</font> (i.e., a metric variable) of people \n",
    "        - randomly sampled from everyone self-described as <font color=\"blue\">vegetarian</font> (i.e., a single group). \n",
    "    - In this case, we could be interested in \n",
    "        - how much this <font color=\"orange\">group’s IQ differs</font> \n",
    "        - from the general population’s average IQ of 100."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"figures/tbl15.1.png\" width=600 />"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* In the context of the generalized linear model (GLM) introduced in the previous chapter, this chapter’s situation involves the most trivial cases of the linear core of the GLM, as indicated in the left cells of Table 15.1 (p. 434), with a link function that is the identity along with a normal distribution for describing noise in the data, as indicated in the first row of Table 15.2 (p. 443). "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"figures/tbl15.2.png\" width=600 />"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* We will explore options for the prior distribution on parameters of the normal distribution, and methods for Bayesian estimation of the parameters. \n",
    "* We will also consider alternative noise distributions for describing data that have outliers."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 16.1. Estimating the Mean and Standard Deviation of a Normal Distribution"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* 16.1.1 Solution by mathematical analysis\n",
    "* 16.1.2 Approximation by MCMC in JAGS"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The normal distribution specifies the probability density of a value y, given the values of two\n",
    "parameters, the mean μ and standard deviation σ :"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"figures/eq16.1.png\" width=600 />"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* To get an intuition for the <font color=\"red\">normal distribution as a likelihood function</font>, consider three data values y1 = 85, y2 = 100, and y3 = 115, which are plotted as large dots in Figure 16.1. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* Figure 16.1 <font color=\"red\">shows p(D|μ, σ ) for different values of μ and σ</font> . As you can see, there are values of μ and σ that make the data most probable, but other nearby values also accommodate the data reasonably well"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"figures/fig16.1.png\" width=600 />"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<font color=\"red\">The question is, given the data, how should we allocate credibility to combinations of μ and σ?</font>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* likelhood \n",
    "    - Figure 16.1 shows examples of p(D|μ, σ ) for a particular data set at different values of μ and σ \n",
    "* prior \n",
    "    -The prior, p(μ,σ), specifies the credibility of each combination of μ,σ values in the two-dimensional joint parameter space, without the data.\n",
    "*  posterior\n",
    "    - Bayes’ rule says that the posterior credibility of each combination of μ, σ values is the prior credibility times the likelihood, normalized by the marginal likelihood."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"figures/eq16.2.png\" width=600 />"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<font color=\"red\">Our goal now is to evaluate Equation 16.2 for reasonable choices of the prior distribution, p(μ, σ ).</font>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 16.1.1 Solution by mathematical analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* We take a short algebraic tour before moving on to MCMC implementations."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### When σ is fixed, "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* 참고 : Conjugate prior\n",
    "    - https://en.wikipedia.org/wiki/Conjugate_prior"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "* It is convenient first to consider the case in which the standard deviation of the likelihood function is fixed at a specific value. In other words, <font color=\"red\">the prior distribution on σ is a spike over that specific value</font>. We’ll denote that <font color=\"red\">fixed value as σ = Sy</font>.  \n",
    "* With this simplifying assumption, <font color=\"red\">we are only estimating μ</font> because we are <font color=\"red\">assuming perfectly certain prior knowledge about σ</font> .\n",
    "* proir : When σ is fixed, then the <font color=\"red\">prior distribution on μ</font> in Equation 16.2 can be easily chosen to be <font color=\"red\">conjugate to the normal likelihood</font>.\n",
    "    - The term “conjugate prior” was defined in Section 6.2, p. 126.\n",
    "    - It turns out that the product of normal distributions is again a normal distribution; \n",
    "    - in other words, if the prior on μ is normal, then the posterior on μ is normal.\n",
    "    - <font color=\"red\">Let the prior distribution on μ be normal with mean Mμ and standard deviation Sμ</font> .\n",
    "* likelihood * prior : "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"figures/eq16.3.png\" width=600 />"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"figures/eq16.4.png\" width=600 />"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"figures/cap16.3.png\" width=600 />"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### precision & posterior precision"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"figures/cap16.4.png\" width=600 />"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"figures/eq16.5.png\" width=600 />"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* <font color=\"red\">Thus, the posterior precision is the sum of the prior precision and the likelihood precision.</font>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### posterior mean"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"figures/cap16.5.png\" width=600 />"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"figures/eq16.6.png\" width=600 />"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* <font color=\"red\">In other words, the posterior mean is a weighted average of the prior mean and the datum, with the weighting corresponding to the relative precisions of the prior and the likelihood.</font>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### N value case"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"figures/cap16.6.png\" width=600 />"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### posterior mean"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"figures/cap16.1.png\" width=600 />"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### posterior precision"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"figures/cap16.2.png\" />"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* Notice that as the sample size N increases, the posterior mean is dominated by the data mean.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### When μ is fixed, "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* We can instead <font color=\"red\">estimate the σ parameter</font> when μ is fixed. \n",
    "* It turns out that when μ is fixed, a <font color=\"red\">conjugate prior for the precision</font> is <font color=\"red\">the gamma distribution</font> (e.g., Gelman et al., 2013, p. 43). \n",
    "    - It is important to <font color=\"red\">understand the meaning of a gamma prior on precision</font>. \n",
    "    - Consider a gamma distribution \n",
    "        - that is loaded <font color=\"orange\">heavily over very small values, but has a long shallow tail extending over large values</font>. \n",
    "            - This sort of  gamma distribution on precision indicates that <font color=\"orange\">we believe most strongly in small precisions</font>, but we <font color=\"orange\">admit that large precisions are possible</font>. \n",
    "            - If this is a <font color=\"orange\">belief about the precision of a normal likelihood function</font>, then this sort of <font color=\"orange\">gamma distribution expresses a belief that the data will be more spread out</font>, because <font color=\"orange\">small precisions imply large standard deviations</font>. \n",
    "        -  If the gamma distribution is instead loaded over <font color=\"blue\">large values of precision</font>, \n",
    "            - it expresses a <font color=\"blue\">belief that the data will be tightly clustered</font>."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### 참고: gamma distribution ?\n",
    "\n",
    "* https://en.wikipedia.org/wiki/Gamma_distribution"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"https://upload.wikimedia.org/wikipedia/commons/f/fc/Gamma_distribution_pdf.png\" width=600 />"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### conjugate priors & gamma distiribution & precision"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* Because of its role in <font color=\"red\">conjugate priors for the normal likelihood function</font>, the <font color=\"red\">gamma distribution is routinely used</font> as a <font color=\"red\">prior on precision</font>. \n",
    "* But there is no logical necessity to do so, and <font color=\"red\">modern MCMC methods permit more flexible</font> specification of priors. \n",
    "* Indeed, <font color=\"red\">because precision is less intuitive than standard deviation</font>, it can be <font color=\"red\">more useful instead to give the standard deviation a uniform prior that spans a wide range</font>."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Summary"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We have assumed that the data are generated by a normal likelihood function,\n",
    "parameterized by a mean μ and standard deviation σ , and denoted y ∼ normal(y|μ, σ ).\n",
    "For purposes of mathematical derivation, we made unrealistic assumptions that the prior\n",
    "distribution is either a spike on σ or a spike on μ, in order to make three main\n",
    "points:\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1. <font color=\"red\">A natural way to express a prior on μ</font> is with a <font color=\"red\">normal distribution</font>, because this is <font color=\"red\">conjugate</font> with the <font color=\"red\">normal likelihood</font> when its <font color=\"red\">standard deviation is fixed</font>.\n",
    "2. <font color=\"red\">A way to express a prior on the precision</font> 1/σ 2 is with a <font color=\"red\">gamma distribution</font>, because this is <font color=\"red\">conjugate with the normal likelihood when its mean is fixed</font>. However <font color=\"blue\">in practice</font> the <font color=\"blue\">standard deviation</font> can instead be <font color=\"blue\">given a uniform prior</font> (or anything else that reflects prior beliefs, of course).\n",
    "3. The <font color=\"red\">formulas for Bayesian updating</font> of the parameter distribution are <font color=\"red\">more conveniently expressed</font> in terms of <font color=\"red\">precision than standard deviation</font>. <font color=\"red\">Normal distributions are described</font> <font color=\"orange\">sometimes</font> in terms of <font color=\"orange\">standard deviation</font> and <font color=\"green\">sometimes</font> in terms of <font color=\"green\">precision</font>, so it is important to glean from context which is being referred to. <font color=\"red\">In R and Stan</font>, the <font color=\"red\">normal distribution</font> is parameterized <font color=\"red\">by mean and standard deviation</font>. <font color=\"blue\">In JAGS and BUGS</font>, the normal distribution is parameterized <font color=\"blue\">by mean and precision</font>."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 16.1.2 Approximation by MCMC in JAGS"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is easy to estimate the mean and standard deviation in JAGS."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* The right panel of Figure 16.2 instead puts a broad uniform distribution directly\n",
    "on σ . The low and high values of the uniform distribution are set to be far outside any realistic\n",
    "value for the data, so that the prior has minimal influence on the posterior. \n",
    "* The uniform prior on σ is easier to intuit than a gamma prior on precision, but the priors are not equivalent."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"figures/fig16.2.png\" width=600 />"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "dataList = list(\n",
    "    y = y ,\n",
    "    Ntotal = length(y) ,\n",
    "    meanY = mean(y) ,\n",
    "    sdY = sd(y)\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "model {\n",
    "    for ( i in 1:Ntotal ) {\n",
    "        y[i]  ̃ dnorm( mu , 1/sigmaˆ2 ) # JAGS uses precision\n",
    "    }\n",
    "    mu  ̃ dnorm( meanY , 1/(100*sdY)ˆ2 ) # JAGS uses precision\n",
    "    sigma  ̃ dunif( sdY/1000 , sdY*1000 )\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* For purposes of illustration, we use fictitious data.\n",
    "* The data are \n",
    "    - <font color=\"red\">IQ</font> (intelligence quotient) scores \n",
    "        - from <font color=\"red\">a group</font> of people who have consumed a <font color=\"red\">“smart drug.”</font> \n",
    "    - We know that IQ tests have been normed to the general population so that they have an <font color=\"blue\">average score of 100 and a standard deviation of 15</font>. \n",
    "* Therefore, we would like to know <font color=\"red\">how differently the smart-drug group</font> has performed relative to the <font color=\"red\">general population average</font>."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"figures/fig16.3.png\" width=600 />"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Jags-Ymet-Xnom1grp-Mnormal-Example.R"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# 책에는 없는 코드. 실습을 위해 data 폴더에 넣은 코드를 실행시키기 위한 작업\n",
    "cur_dir = getwd()\n",
    "setwd(sprintf(\"%s/%s\", cur_dir, 'data'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Example for Jags-Ymet-Xnom1grp-Mnormal.R \n",
    "#------------------------------------------------------------------------------- \n",
    "# Optional generic preliminaries:\n",
    "#graphics.off() # This closes all of R's graphics windows.\n",
    "#rm(list=ls())  # Careful! This clears all of R's memory!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#------------------------------------------------------------------------------- \n",
    "# Load The data file \n",
    "myDataFrame = read.csv( file=\"TwoGroupIQ.csv\" )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table>\n",
       "<thead><tr><th></th><th scope=col>Score</th><th scope=col>Group</th></tr></thead>\n",
       "<tbody>\n",
       "\t<tr><th scope=row>1</th><td>102</td><td>Smart Drug</td></tr>\n",
       "\t<tr><th scope=row>2</th><td>107</td><td>Smart Drug</td></tr>\n",
       "\t<tr><th scope=row>3</th><td>92</td><td>Smart Drug</td></tr>\n",
       "\t<tr><th scope=row>4</th><td>101</td><td>Smart Drug</td></tr>\n",
       "\t<tr><th scope=row>5</th><td>110</td><td>Smart Drug</td></tr>\n",
       "\t<tr><th scope=row>6</th><td>68</td><td>Smart Drug</td></tr>\n",
       "\t<tr><th scope=row>7</th><td>119</td><td>Smart Drug</td></tr>\n",
       "\t<tr><th scope=row>8</th><td>106</td><td>Smart Drug</td></tr>\n",
       "\t<tr><th scope=row>9</th><td>99</td><td>Smart Drug</td></tr>\n",
       "\t<tr><th scope=row>10</th><td>103</td><td>Smart Drug</td></tr>\n",
       "\t<tr><th scope=row>11</th><td>90</td><td>Smart Drug</td></tr>\n",
       "\t<tr><th scope=row>12</th><td>93</td><td>Smart Drug</td></tr>\n",
       "\t<tr><th scope=row>13</th><td>79</td><td>Smart Drug</td></tr>\n",
       "\t<tr><th scope=row>14</th><td>89</td><td>Smart Drug</td></tr>\n",
       "\t<tr><th scope=row>15</th><td>137</td><td>Smart Drug</td></tr>\n",
       "\t<tr><th scope=row>16</th><td>119</td><td>Smart Drug</td></tr>\n",
       "\t<tr><th scope=row>17</th><td>126</td><td>Smart Drug</td></tr>\n",
       "\t<tr><th scope=row>18</th><td>110</td><td>Smart Drug</td></tr>\n",
       "\t<tr><th scope=row>19</th><td>71</td><td>Smart Drug</td></tr>\n",
       "\t<tr><th scope=row>20</th><td>114</td><td>Smart Drug</td></tr>\n",
       "\t<tr><th scope=row>21</th><td>100</td><td>Smart Drug</td></tr>\n",
       "\t<tr><th scope=row>22</th><td>95</td><td>Smart Drug</td></tr>\n",
       "\t<tr><th scope=row>23</th><td>91</td><td>Smart Drug</td></tr>\n",
       "\t<tr><th scope=row>24</th><td>99</td><td>Smart Drug</td></tr>\n",
       "\t<tr><th scope=row>25</th><td>97</td><td>Smart Drug</td></tr>\n",
       "\t<tr><th scope=row>26</th><td>106</td><td>Smart Drug</td></tr>\n",
       "\t<tr><th scope=row>27</th><td>106</td><td>Smart Drug</td></tr>\n",
       "\t<tr><th scope=row>28</th><td>129</td><td>Smart Drug</td></tr>\n",
       "\t<tr><th scope=row>29</th><td>115</td><td>Smart Drug</td></tr>\n",
       "\t<tr><th scope=row>30</th><td>124</td><td>Smart Drug</td></tr>\n",
       "\t<tr><th scope=row>31</th><td>137</td><td>Smart Drug</td></tr>\n",
       "\t<tr><th scope=row>32</th><td>73</td><td>Smart Drug</td></tr>\n",
       "\t<tr><th scope=row>33</th><td>69</td><td>Smart Drug</td></tr>\n",
       "\t<tr><th scope=row>34</th><td>95</td><td>Smart Drug</td></tr>\n",
       "\t<tr><th scope=row>35</th><td>102</td><td>Smart Drug</td></tr>\n",
       "\t<tr><th scope=row>36</th><td>116</td><td>Smart Drug</td></tr>\n",
       "\t<tr><th scope=row>37</th><td>111</td><td>Smart Drug</td></tr>\n",
       "\t<tr><th scope=row>38</th><td>134</td><td>Smart Drug</td></tr>\n",
       "\t<tr><th scope=row>39</th><td>102</td><td>Smart Drug</td></tr>\n",
       "\t<tr><th scope=row>40</th><td>110</td><td>Smart Drug</td></tr>\n",
       "\t<tr><th scope=row>41</th><td>139</td><td>Smart Drug</td></tr>\n",
       "\t<tr><th scope=row>42</th><td>112</td><td>Smart Drug</td></tr>\n",
       "\t<tr><th scope=row>43</th><td>122</td><td>Smart Drug</td></tr>\n",
       "\t<tr><th scope=row>44</th><td>84</td><td>Smart Drug</td></tr>\n",
       "\t<tr><th scope=row>45</th><td>129</td><td>Smart Drug</td></tr>\n",
       "\t<tr><th scope=row>46</th><td>112</td><td>Smart Drug</td></tr>\n",
       "\t<tr><th scope=row>47</th><td>127</td><td>Smart Drug</td></tr>\n",
       "\t<tr><th scope=row>48</th><td>106</td><td>Smart Drug</td></tr>\n",
       "\t<tr><th scope=row>49</th><td>113</td><td>Smart Drug</td></tr>\n",
       "\t<tr><th scope=row>50</th><td>109</td><td>Smart Drug</td></tr>\n",
       "\t<tr><th scope=row>51</th><td>208</td><td>Smart Drug</td></tr>\n",
       "\t<tr><th scope=row>52</th><td>114</td><td>Smart Drug</td></tr>\n",
       "\t<tr><th scope=row>53</th><td>107</td><td>Smart Drug</td></tr>\n",
       "\t<tr><th scope=row>54</th><td>50</td><td>Smart Drug</td></tr>\n",
       "\t<tr><th scope=row>55</th><td>169</td><td>Smart Drug</td></tr>\n",
       "\t<tr><th scope=row>56</th><td>133</td><td>Smart Drug</td></tr>\n",
       "\t<tr><th scope=row>57</th><td>50</td><td>Smart Drug</td></tr>\n",
       "\t<tr><th scope=row>58</th><td>97</td><td>Smart Drug</td></tr>\n",
       "\t<tr><th scope=row>59</th><td>139</td><td>Smart Drug</td></tr>\n",
       "\t<tr><th scope=row>60</th><td>72</td><td>Smart Drug</td></tr>\n",
       "\t<tr><th scope=row>61</th><td>100</td><td>Smart Drug</td></tr>\n",
       "\t<tr><th scope=row>62</th><td>144</td><td>Smart Drug</td></tr>\n",
       "\t<tr><th scope=row>63</th><td>112</td><td>Smart Drug</td></tr>\n",
       "\t<tr><th scope=row>64</th><td>109</td><td>Placebo</td></tr>\n",
       "\t<tr><th scope=row>65</th><td>98</td><td>Placebo</td></tr>\n",
       "\t<tr><th scope=row>66</th><td>106</td><td>Placebo</td></tr>\n",
       "\t<tr><th scope=row>67</th><td>101</td><td>Placebo</td></tr>\n",
       "\t<tr><th scope=row>68</th><td>100</td><td>Placebo</td></tr>\n",
       "\t<tr><th scope=row>69</th><td>111</td><td>Placebo</td></tr>\n",
       "\t<tr><th scope=row>70</th><td>117</td><td>Placebo</td></tr>\n",
       "\t<tr><th scope=row>71</th><td>104</td><td>Placebo</td></tr>\n",
       "\t<tr><th scope=row>72</th><td>106</td><td>Placebo</td></tr>\n",
       "\t<tr><th scope=row>73</th><td>89</td><td>Placebo</td></tr>\n",
       "\t<tr><th scope=row>74</th><td>84</td><td>Placebo</td></tr>\n",
       "\t<tr><th scope=row>75</th><td>88</td><td>Placebo</td></tr>\n",
       "\t<tr><th scope=row>76</th><td>94</td><td>Placebo</td></tr>\n",
       "\t<tr><th scope=row>77</th><td>78</td><td>Placebo</td></tr>\n",
       "\t<tr><th scope=row>78</th><td>108</td><td>Placebo</td></tr>\n",
       "\t<tr><th scope=row>79</th><td>102</td><td>Placebo</td></tr>\n",
       "\t<tr><th scope=row>80</th><td>95</td><td>Placebo</td></tr>\n",
       "\t<tr><th scope=row>81</th><td>99</td><td>Placebo</td></tr>\n",
       "\t<tr><th scope=row>82</th><td>90</td><td>Placebo</td></tr>\n",
       "\t<tr><th scope=row>83</th><td>116</td><td>Placebo</td></tr>\n",
       "\t<tr><th scope=row>84</th><td>97</td><td>Placebo</td></tr>\n",
       "\t<tr><th scope=row>85</th><td>107</td><td>Placebo</td></tr>\n",
       "\t<tr><th scope=row>86</th><td>102</td><td>Placebo</td></tr>\n",
       "\t<tr><th scope=row>87</th><td>91</td><td>Placebo</td></tr>\n",
       "\t<tr><th scope=row>88</th><td>94</td><td>Placebo</td></tr>\n",
       "\t<tr><th scope=row>89</th><td>95</td><td>Placebo</td></tr>\n",
       "\t<tr><th scope=row>90</th><td>86</td><td>Placebo</td></tr>\n",
       "\t<tr><th scope=row>91</th><td>108</td><td>Placebo</td></tr>\n",
       "\t<tr><th scope=row>92</th><td>115</td><td>Placebo</td></tr>\n",
       "\t<tr><th scope=row>93</th><td>108</td><td>Placebo</td></tr>\n",
       "\t<tr><th scope=row>94</th><td>88</td><td>Placebo</td></tr>\n",
       "\t<tr><th scope=row>95</th><td>102</td><td>Placebo</td></tr>\n",
       "\t<tr><th scope=row>96</th><td>102</td><td>Placebo</td></tr>\n",
       "\t<tr><th scope=row>97</th><td>120</td><td>Placebo</td></tr>\n",
       "\t<tr><th scope=row>98</th><td>112</td><td>Placebo</td></tr>\n",
       "\t<tr><th scope=row>99</th><td>100</td><td>Placebo</td></tr>\n",
       "\t<tr><th scope=row>100</th><td>105</td><td>Placebo</td></tr>\n",
       "\t<tr><th scope=row>101</th><td>105</td><td>Placebo</td></tr>\n",
       "\t<tr><th scope=row>102</th><td>88</td><td>Placebo</td></tr>\n",
       "\t<tr><th scope=row>103</th><td>82</td><td>Placebo</td></tr>\n",
       "\t<tr><th scope=row>104</th><td>111</td><td>Placebo</td></tr>\n",
       "\t<tr><th scope=row>105</th><td>96</td><td>Placebo</td></tr>\n",
       "\t<tr><th scope=row>106</th><td>92</td><td>Placebo</td></tr>\n",
       "\t<tr><th scope=row>107</th><td>109</td><td>Placebo</td></tr>\n",
       "\t<tr><th scope=row>108</th><td>91</td><td>Placebo</td></tr>\n",
       "\t<tr><th scope=row>109</th><td>92</td><td>Placebo</td></tr>\n",
       "\t<tr><th scope=row>110</th><td>123</td><td>Placebo</td></tr>\n",
       "\t<tr><th scope=row>111</th><td>61</td><td>Placebo</td></tr>\n",
       "\t<tr><th scope=row>112</th><td>59</td><td>Placebo</td></tr>\n",
       "\t<tr><th scope=row>113</th><td>105</td><td>Placebo</td></tr>\n",
       "\t<tr><th scope=row>114</th><td>184</td><td>Placebo</td></tr>\n",
       "\t<tr><th scope=row>115</th><td>82</td><td>Placebo</td></tr>\n",
       "\t<tr><th scope=row>116</th><td>138</td><td>Placebo</td></tr>\n",
       "\t<tr><th scope=row>117</th><td>99</td><td>Placebo</td></tr>\n",
       "\t<tr><th scope=row>118</th><td>93</td><td>Placebo</td></tr>\n",
       "\t<tr><th scope=row>119</th><td>93</td><td>Placebo</td></tr>\n",
       "\t<tr><th scope=row>120</th><td>72</td><td>Placebo</td></tr>\n",
       "</tbody>\n",
       "</table>\n"
      ],
      "text/latex": [
       "\\begin{tabular}{r|ll}\n",
       "  & Score & Group\\\\\n",
       "\\hline\n",
       "\t1 & 102 & Smart Drug\\\\\n",
       "\t2 & 107 & Smart Drug\\\\\n",
       "\t3 & 92 & Smart Drug\\\\\n",
       "\t4 & 101 & Smart Drug\\\\\n",
       "\t5 & 110 & Smart Drug\\\\\n",
       "\t6 & 68 & Smart Drug\\\\\n",
       "\t7 & 119 & Smart Drug\\\\\n",
       "\t8 & 106 & Smart Drug\\\\\n",
       "\t9 & 99 & Smart Drug\\\\\n",
       "\t10 & 103 & Smart Drug\\\\\n",
       "\t11 & 90 & Smart Drug\\\\\n",
       "\t12 & 93 & Smart Drug\\\\\n",
       "\t13 & 79 & Smart Drug\\\\\n",
       "\t14 & 89 & Smart Drug\\\\\n",
       "\t15 & 137 & Smart Drug\\\\\n",
       "\t16 & 119 & Smart Drug\\\\\n",
       "\t17 & 126 & Smart Drug\\\\\n",
       "\t18 & 110 & Smart Drug\\\\\n",
       "\t19 & 71 & Smart Drug\\\\\n",
       "\t20 & 114 & Smart Drug\\\\\n",
       "\t21 & 100 & Smart Drug\\\\\n",
       "\t22 & 95 & Smart Drug\\\\\n",
       "\t23 & 91 & Smart Drug\\\\\n",
       "\t24 & 99 & Smart Drug\\\\\n",
       "\t25 & 97 & Smart Drug\\\\\n",
       "\t26 & 106 & Smart Drug\\\\\n",
       "\t27 & 106 & Smart Drug\\\\\n",
       "\t28 & 129 & Smart Drug\\\\\n",
       "\t29 & 115 & Smart Drug\\\\\n",
       "\t30 & 124 & Smart Drug\\\\\n",
       "\t31 & 137 & Smart Drug\\\\\n",
       "\t32 & 73 & Smart Drug\\\\\n",
       "\t33 & 69 & Smart Drug\\\\\n",
       "\t34 & 95 & Smart Drug\\\\\n",
       "\t35 & 102 & Smart Drug\\\\\n",
       "\t36 & 116 & Smart Drug\\\\\n",
       "\t37 & 111 & Smart Drug\\\\\n",
       "\t38 & 134 & Smart Drug\\\\\n",
       "\t39 & 102 & Smart Drug\\\\\n",
       "\t40 & 110 & Smart Drug\\\\\n",
       "\t41 & 139 & Smart Drug\\\\\n",
       "\t42 & 112 & Smart Drug\\\\\n",
       "\t43 & 122 & Smart Drug\\\\\n",
       "\t44 & 84 & Smart Drug\\\\\n",
       "\t45 & 129 & Smart Drug\\\\\n",
       "\t46 & 112 & Smart Drug\\\\\n",
       "\t47 & 127 & Smart Drug\\\\\n",
       "\t48 & 106 & Smart Drug\\\\\n",
       "\t49 & 113 & Smart Drug\\\\\n",
       "\t50 & 109 & Smart Drug\\\\\n",
       "\t51 & 208 & Smart Drug\\\\\n",
       "\t52 & 114 & Smart Drug\\\\\n",
       "\t53 & 107 & Smart Drug\\\\\n",
       "\t54 & 50 & Smart Drug\\\\\n",
       "\t55 & 169 & Smart Drug\\\\\n",
       "\t56 & 133 & Smart Drug\\\\\n",
       "\t57 & 50 & Smart Drug\\\\\n",
       "\t58 & 97 & Smart Drug\\\\\n",
       "\t59 & 139 & Smart Drug\\\\\n",
       "\t60 & 72 & Smart Drug\\\\\n",
       "\t61 & 100 & Smart Drug\\\\\n",
       "\t62 & 144 & Smart Drug\\\\\n",
       "\t63 & 112 & Smart Drug\\\\\n",
       "\t64 & 109 & Placebo\\\\\n",
       "\t65 & 98 & Placebo\\\\\n",
       "\t66 & 106 & Placebo\\\\\n",
       "\t67 & 101 & Placebo\\\\\n",
       "\t68 & 100 & Placebo\\\\\n",
       "\t69 & 111 & Placebo\\\\\n",
       "\t70 & 117 & Placebo\\\\\n",
       "\t71 & 104 & Placebo\\\\\n",
       "\t72 & 106 & Placebo\\\\\n",
       "\t73 & 89 & Placebo\\\\\n",
       "\t74 & 84 & Placebo\\\\\n",
       "\t75 & 88 & Placebo\\\\\n",
       "\t76 & 94 & Placebo\\\\\n",
       "\t77 & 78 & Placebo\\\\\n",
       "\t78 & 108 & Placebo\\\\\n",
       "\t79 & 102 & Placebo\\\\\n",
       "\t80 & 95 & Placebo\\\\\n",
       "\t81 & 99 & Placebo\\\\\n",
       "\t82 & 90 & Placebo\\\\\n",
       "\t83 & 116 & Placebo\\\\\n",
       "\t84 & 97 & Placebo\\\\\n",
       "\t85 & 107 & Placebo\\\\\n",
       "\t86 & 102 & Placebo\\\\\n",
       "\t87 & 91 & Placebo\\\\\n",
       "\t88 & 94 & Placebo\\\\\n",
       "\t89 & 95 & Placebo\\\\\n",
       "\t90 & 86 & Placebo\\\\\n",
       "\t91 & 108 & Placebo\\\\\n",
       "\t92 & 115 & Placebo\\\\\n",
       "\t93 & 108 & Placebo\\\\\n",
       "\t94 & 88 & Placebo\\\\\n",
       "\t95 & 102 & Placebo\\\\\n",
       "\t96 & 102 & Placebo\\\\\n",
       "\t97 & 120 & Placebo\\\\\n",
       "\t98 & 112 & Placebo\\\\\n",
       "\t99 & 100 & Placebo\\\\\n",
       "\t100 & 105 & Placebo\\\\\n",
       "\t101 & 105 & Placebo\\\\\n",
       "\t102 & 88 & Placebo\\\\\n",
       "\t103 & 82 & Placebo\\\\\n",
       "\t104 & 111 & Placebo\\\\\n",
       "\t105 & 96 & Placebo\\\\\n",
       "\t106 & 92 & Placebo\\\\\n",
       "\t107 & 109 & Placebo\\\\\n",
       "\t108 & 91 & Placebo\\\\\n",
       "\t109 & 92 & Placebo\\\\\n",
       "\t110 & 123 & Placebo\\\\\n",
       "\t111 & 61 & Placebo\\\\\n",
       "\t112 & 59 & Placebo\\\\\n",
       "\t113 & 105 & Placebo\\\\\n",
       "\t114 & 184 & Placebo\\\\\n",
       "\t115 & 82 & Placebo\\\\\n",
       "\t116 & 138 & Placebo\\\\\n",
       "\t117 & 99 & Placebo\\\\\n",
       "\t118 & 93 & Placebo\\\\\n",
       "\t119 & 93 & Placebo\\\\\n",
       "\t120 & 72 & Placebo\\\\\n",
       "\\end{tabular}\n"
      ],
      "text/plain": [
       "    Score      Group\n",
       "1     102 Smart Drug\n",
       "2     107 Smart Drug\n",
       "3      92 Smart Drug\n",
       "4     101 Smart Drug\n",
       "5     110 Smart Drug\n",
       "6      68 Smart Drug\n",
       "7     119 Smart Drug\n",
       "8     106 Smart Drug\n",
       "9      99 Smart Drug\n",
       "10    103 Smart Drug\n",
       "11     90 Smart Drug\n",
       "12     93 Smart Drug\n",
       "13     79 Smart Drug\n",
       "14     89 Smart Drug\n",
       "15    137 Smart Drug\n",
       "16    119 Smart Drug\n",
       "17    126 Smart Drug\n",
       "18    110 Smart Drug\n",
       "19     71 Smart Drug\n",
       "20    114 Smart Drug\n",
       "21    100 Smart Drug\n",
       "22     95 Smart Drug\n",
       "23     91 Smart Drug\n",
       "24     99 Smart Drug\n",
       "25     97 Smart Drug\n",
       "26    106 Smart Drug\n",
       "27    106 Smart Drug\n",
       "28    129 Smart Drug\n",
       "29    115 Smart Drug\n",
       "30    124 Smart Drug\n",
       "31    137 Smart Drug\n",
       "32     73 Smart Drug\n",
       "33     69 Smart Drug\n",
       "34     95 Smart Drug\n",
       "35    102 Smart Drug\n",
       "36    116 Smart Drug\n",
       "37    111 Smart Drug\n",
       "38    134 Smart Drug\n",
       "39    102 Smart Drug\n",
       "40    110 Smart Drug\n",
       "41    139 Smart Drug\n",
       "42    112 Smart Drug\n",
       "43    122 Smart Drug\n",
       "44     84 Smart Drug\n",
       "45    129 Smart Drug\n",
       "46    112 Smart Drug\n",
       "47    127 Smart Drug\n",
       "48    106 Smart Drug\n",
       "49    113 Smart Drug\n",
       "50    109 Smart Drug\n",
       "51    208 Smart Drug\n",
       "52    114 Smart Drug\n",
       "53    107 Smart Drug\n",
       "54     50 Smart Drug\n",
       "55    169 Smart Drug\n",
       "56    133 Smart Drug\n",
       "57     50 Smart Drug\n",
       "58     97 Smart Drug\n",
       "59    139 Smart Drug\n",
       "60     72 Smart Drug\n",
       "61    100 Smart Drug\n",
       "62    144 Smart Drug\n",
       "63    112 Smart Drug\n",
       "64    109    Placebo\n",
       "65     98    Placebo\n",
       "66    106    Placebo\n",
       "67    101    Placebo\n",
       "68    100    Placebo\n",
       "69    111    Placebo\n",
       "70    117    Placebo\n",
       "71    104    Placebo\n",
       "72    106    Placebo\n",
       "73     89    Placebo\n",
       "74     84    Placebo\n",
       "75     88    Placebo\n",
       "76     94    Placebo\n",
       "77     78    Placebo\n",
       "78    108    Placebo\n",
       "79    102    Placebo\n",
       "80     95    Placebo\n",
       "81     99    Placebo\n",
       "82     90    Placebo\n",
       "83    116    Placebo\n",
       "84     97    Placebo\n",
       "85    107    Placebo\n",
       "86    102    Placebo\n",
       "87     91    Placebo\n",
       "88     94    Placebo\n",
       "89     95    Placebo\n",
       "90     86    Placebo\n",
       "91    108    Placebo\n",
       "92    115    Placebo\n",
       "93    108    Placebo\n",
       "94     88    Placebo\n",
       "95    102    Placebo\n",
       "96    102    Placebo\n",
       "97    120    Placebo\n",
       "98    112    Placebo\n",
       "99    100    Placebo\n",
       "100   105    Placebo\n",
       "101   105    Placebo\n",
       "102    88    Placebo\n",
       "103    82    Placebo\n",
       "104   111    Placebo\n",
       "105    96    Placebo\n",
       "106    92    Placebo\n",
       "107   109    Placebo\n",
       "108    91    Placebo\n",
       "109    92    Placebo\n",
       "110   123    Placebo\n",
       "111    61    Placebo\n",
       "112    59    Placebo\n",
       "113   105    Placebo\n",
       "114   184    Placebo\n",
       "115    82    Placebo\n",
       "116   138    Placebo\n",
       "117    99    Placebo\n",
       "118    93    Placebo\n",
       "119    93    Placebo\n",
       "120    72    Placebo"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "myDataFrame"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "     Score               Group   \n",
       " Min.   : 50.00   Placebo   :57  \n",
       " 1st Qu.: 92.75   Smart Drug:63  \n",
       " Median :102.00                  \n",
       " Mean   :104.13                  \n",
       " 3rd Qu.:112.00                  \n",
       " Max.   :208.00                  "
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "summary(myDataFrame)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# For purposes of this one-group example, use data from Smart Drug group:\n",
    "myData = myDataFrame$Score[myDataFrame$Group==\"Smart Drug\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#------------------------------------------------------------------------------- \n",
    "# Optional: Specify filename root and graphical format for saving output.\n",
    "# Otherwise specify as NULL or leave saveName and saveType arguments \n",
    "# out of function calls.\n",
    "fileNameRoot = \"OneGroupIQnormal-\" \n",
    "graphFileType = \"png\"  #\"eps\" "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "*********************************************************************\n",
      "Kruschke, J. K. (2015). Doing Bayesian Data Analysis, Second Edition:\n",
      "A Tutorial with R, JAGS, and Stan. Academic Press / Elsevier.\n",
      "*********************************************************************\n",
      "\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Loading required package: coda\n"
     ]
    }
   ],
   "source": [
    "#------------------------------------------------------------------------------- \n",
    "# Load the relevant model into R's working memory:\n",
    "source(\"Jags-Ymet-Xnom1grp-Mnormal.R\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Loading required package: rjags\n",
      "Linked to JAGS 3.4.0\n",
      "Loaded modules: basemod,bugs\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Compiling model graph\n",
      "   Resolving undeclared variables\n",
      "   Allocating nodes\n",
      "   Graph Size: 79\n",
      "\n",
      "Initializing model\n",
      "\n",
      "Burning in the MCMC chain...\n",
      "Sampling final MCMC chain...\n"
     ]
    }
   ],
   "source": [
    "#------------------------------------------------------------------------------- \n",
    "# Generate the MCMC chain:\n",
    "mcmcCoda = genMCMC( data=myData , numSavedSteps=20000 , saveName=fileNameRoot )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#------------------------------------------------------------------------------- \n",
    "# Display diagnostics of chain, for specified parameters:\n",
    "parameterNames = varnames(mcmcCoda) # get all parameter names\n",
    "for ( parName in parameterNames ) {\n",
    "  diagMCMC( codaObject=mcmcCoda , parName=parName , \n",
    "            saveName=fileNameRoot , saveType=graphFileType )\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# 창으로 뜨는 것들을 없앤다.\n",
    "graphics.off()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"data/OneGroupIQnormal-Diagmu.png\" width=600 />\n",
    "<img src=\"data/OneGroupIQnormal-Diagsigma.png\" width=600 />"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "             Mean      Median        Mode     ESS HDImass       HDIlow\n",
      "mu    107.8278271 107.8114667 107.0894331 20097.3    0.95 101.43788497\n",
      "sigma  25.9588549  25.7791320  25.0124458 11475.4    0.95  21.46449101\n",
      "effSz   0.3039062   0.3034737   0.2807922 20000.0    0.95   0.04810492\n",
      "          HDIhigh CompVal PcntGtCompVal ROPElow ROPEhigh PcntLtROPE PcntInROPE\n",
      "mu    114.4334683     100         98.95    99.0    101.0       0.41      1.625\n",
      "sigma  30.6368277      15        100.00    14.0     16.0       0.00      0.000\n",
      "effSz   0.5540163       0         98.95    -0.1      0.1       0.07      5.650\n",
      "      PcntGtROPE\n",
      "mu        97.965\n",
      "sigma    100.000\n",
      "effSz     94.280\n"
     ]
    }
   ],
   "source": [
    "#------------------------------------------------------------------------------- \n",
    "# Get summary statistics of chain:\n",
    "summaryInfo = smryMCMC( mcmcCoda , \n",
    "                        compValMu=100.0 , ropeMu=c(99.0,101.0) ,\n",
    "                        compValSigma=15.0 , ropeSigma=c(14,16) ,\n",
    "                        compValEff=0.0 , ropeEff=c(-0.1,0.1) ,\n",
    "                        saveName=fileNameRoot )\n",
    "show(summaryInfo)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Display posterior information:\n",
    "plotMCMC( mcmcCoda , data=myData , \n",
    "          compValMu=100.0 , ropeMu=c(99.0,101.0) ,\n",
    "          compValSigma=15.0 , ropeSigma=c(14,16) ,\n",
    "          compValEff=0.0 , ropeEff=c(-0.1,0.1) ,\n",
    "          pairsPlot=TRUE , showCurve=FALSE ,\n",
    "          saveName=fileNameRoot , saveType=graphFileType )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "graphics.off()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"data/OneGroupIQnormal-Post.png\" width=600 />\n",
    "<img src=\"data/OneGroupIQnormal-PostPairs.png\" width=600 />"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 16.2. Outliers and Robust Estimation : The t Distribution"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* 16.2.1 Using the t distribution in JAGS\n",
    "* 16.2.2 Using the t distribution in Stan"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### When data appear to have outliers \n",
    "beyond what would be accommodated by a normal distribution, it would be useful to be able to describe the data with a more appropriate distribution that has taller or heavier tails than the normal. \n",
    "* A well-known distribution with heavy tails is the t distribution[2-6].\n",
    "    - Figure 16.4 shows examples of the t distribution. \n",
    "    - Like the normal distribution, it has two parameters \n",
    "        - that control its <font color=\"red\">mean</font> and its width. \n",
    "        - The standard deviation is controlled indirectly via the t distribution’s <font color=\"red\">“scale” parameter</font>. \n",
    "        - In Figure 16.4, the mean is set to zero and the scale is set to one.\n",
    "    - The t distribution has a third parameter \n",
    "        - that controls the heaviness of its tails, \n",
    "            - which I will refer to as the <font color=\"red\">“normality” parameter</font>, \n",
    "            - ν (Greek letter nu). \n",
    "            - Many people might be familiar with this parameter as the “degrees of freedom” from its use in NHST. \n",
    "               - But because \n",
    "                   - we will not be using the t distribution as a sampling distribution, \n",
    "                   - and instead we will be using it only as a descriptive shape, \n",
    "               - I prefer to name the parameter by its effect on the distribution’s shape. \n",
    "            - The normality parameter can range continuously from 1 to ∞. As can be seen in Figure 16.4, \n",
    "                - when ν = 1 the t distribution has very heavy tails, and \n",
    "                - when ν approaches ∞ the t distribution becomes normal."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"figures/fig16.4.png\" width=600 />"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* Figure 16.5 shows examples of how the t distribution is robust against outliers. \n",
    "    - D = {yi}\n",
    "    - p(D|μ, σ ) <- normal\n",
    "    - p(D|μ, σ , ν) <- t"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"figures/fig16.5.png\" width=600 />"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* It is important to understand that <font color=\"red\">the scale parameter σ</font> in the t distribution <font color=\"red\">is not the standard deviation</font> of the distribution. \n",
    "* The standard deviation is actually larger than σ because of the heavy tails. \n",
    "    - In fact, when ν drops below 2 (but is still ≥ 1), the standard deviation of the mathematical t distribution goes to infinity. \n",
    "        - For example, in the upper panel of Figure 16.5, ν is only 1.14, which means that the standard deviation of the mathematical t distribution is infinity, even though the sample standard deviation of the data is finite. \n",
    "            - At the same time, the scale parameter of the t distribution has value σ = 1.47.\n",
    "* While this value of the scale parameter is not the standard deviation of the distribution, it does have an <font color=\"red\">intuitive relation to the spread of the data</font>. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"figures/fig16.6.png\" width=600 />"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* The lower panel of Figure 16.5 uses realistic data that indicates levels of inorganic phosphorous, measured in milligrams per deciliter, in 177 human subjects aged 65 or older. \n",
    "    - The authors of the data (Holcomb & Spalsbury, 2005) intentionally altered a few data points to reflect typical transcription errors and to illustrate methods for detecting and correcting such errors.\n",
    "    - We instead assume that we no longer have access to records of the original individual measurements, and must model the uncorrected data set.\n",
    "* The t distribution accommodates the outliers and fits the distribution of data much better than the normal."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### robust estimation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* The use of a heavy-tailed distribution is often called robust estimation because the estimated value of the central tendency is stable, that is, “robust,” against outliers.\n",
    "* The t distribution is useful \n",
    "    - as a <font color=\"red\">likelihood function</font> for modeling outliers at the level of observed data. But the t distribution is also useful\n",
    "    - for modeling outliers at higher levels in a <font color=\"red\">hierarchical prior</font>. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 16.2.1 Using the t distribution in JAGS"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### using normal "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"figures/fig16.2.png\" width=600 />"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "dataList = list(\n",
    "    y = y ,\n",
    "    Ntotal = length(y) ,\n",
    "    meanY = mean(y) ,\n",
    "    sdY = sd(y)\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "model {\n",
    "    for ( i in 1:Ntotal ) {\n",
    "        y[i]  ̃ dnorm( mu , 1/sigmaˆ2 ) # JAGS uses precision\n",
    "    }\n",
    "    mu  ̃ dnorm( meanY , 1/(100*sdY)ˆ2 ) # JAGS uses precision\n",
    "    sigma  ̃ dunif( sdY/1000 , sdY*1000 )\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### using t "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "model {\n",
    "    for ( i in 1:Ntotal ) {\n",
    "        y[i]  ̃ dt( mu , 1/sigmaˆ2 , nu ) # JAGS: dt uses precision\n",
    "    }\n",
    "    mu  ̃ dnorm( meanY , 1/(100*sdY)ˆ2 )\n",
    "    sigma  ̃ dunif( sdY/1000 , sdY*1000 )\n",
    "    nu <- nuMinusOne+1       # nu must be >= 1\n",
    "    nuMinusOne  ̃ dexp(1/29)  # prior on nu-1\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"figures/fig16.7.png\" width=600 />"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"figures/fig16.8.png\" width=600 />"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Jags-Ymet-Xnom1grp-Mrobust-Example.R"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"figures/fig16.9.png\" width=600 />"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# 책에는 없는 코드. 실습을 위해 data 폴더에 넣은 코드를 실행시키기 위한 작업\n",
    "cur_dir = getwd()\n",
    "setwd(sprintf(\"%s/%s\", cur_dir, 'data'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#------------------------------------------------------------------------------- \n",
    "# Optional generic preliminaries:\n",
    "graphics.off() # This closes all of R's graphics windows.\n",
    "rm(list=ls())  # Careful! This clears all of R's memory!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#------------------------------------------------------------------------------- \n",
    "# Load The data file \n",
    "myDataFrame = read.csv( file=\"TwoGroupIQ.csv\" )\n",
    "# For purposes of this one-group example, use data from Smart Drug group:\n",
    "myData = myDataFrame$Score[myDataFrame$Group==\"Smart Drug\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#------------------------------------------------------------------------------- \n",
    "# Optional: Specify filename root and graphical format for saving output.\n",
    "# Otherwise specify as NULL or leave saveName and saveType arguments \n",
    "# out of function calls.\n",
    "fileNameRoot = \"OneGroupIQrobust-Jags-\" \n",
    "graphFileType = \"png\" #\"eps\" "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "*********************************************************************\n",
      "Kruschke, J. K. (2015). Doing Bayesian Data Analysis, Second Edition:\n",
      "A Tutorial with R, JAGS, and Stan. Academic Press / Elsevier.\n",
      "*********************************************************************\n",
      "\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Installing packages into ‘/home/moosung/R/x86_64-pc-linux-gnu-library/3.2’\n",
      "(as ‘lib’ is unspecified)\n"
     ]
    },
    {
     "ename": "ERROR",
     "evalue": "Error in contrib.url(repos, type): trying to use CRAN without setting a mirror\n",
     "output_type": "error",
     "traceback": [
      "Error in contrib.url(repos, type): trying to use CRAN without setting a mirror\n"
     ]
    }
   ],
   "source": [
    "#------------------------------------------------------------------------------- \n",
    "# Load the relevant model into R's working memory:\n",
    "source(\"Jags-Ymet-Xnom1grp-Mrobust.R\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#------------------------------------------------------------------------------- \n",
    "# Generate the MCMC chain:\n",
    "mcmcCoda = genMCMC( data=myData , numSavedSteps=20000 , saveName=fileNameRoot )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#------------------------------------------------------------------------------- \n",
    "# Display diagnostics of chain, for specified parameters:\n",
    "parameterNames = varnames(mcmcCoda) # get all parameter names\n",
    "for ( parName in parameterNames ) {\n",
    "  diagMCMC( codaObject=mcmcCoda , parName=parName , \n",
    "            saveName=fileNameRoot , saveType=graphFileType )\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# 창으로 뜨는 것들을 없앤다.\n",
    "graphics.off()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#------------------------------------------------------------------------------- \n",
    "# Get summary statistics of chain:\n",
    "summaryInfo = smryMCMC( mcmcCoda , \n",
    "                        compValMu=100.0 , ropeMu=c(99.0,101.0) ,\n",
    "                        compValSigma=15.0 , ropeSigma=c(14,16) ,\n",
    "                        compValEff=0.0 , ropeEff=c(-0.1,0.1) ,\n",
    "                        saveName=fileNameRoot )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "show(summaryInfo)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Display posterior information:\n",
    "plotMCMC( mcmcCoda , data=myData , \n",
    "          compValMu=100.0 , ropeMu=c(99.0,101.0) ,\n",
    "          compValSigma=15.0 , ropeSigma=c(14,16) ,\n",
    "          compValEff=0.0 , ropeEff=c(-0.1,0.1) ,\n",
    "          pairsPlot=TRUE , showCurve=FALSE ,\n",
    "          saveName=fileNameRoot , saveType=graphFileType )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "graphics.off()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 16.2.2 Using the t distribution in Stan"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* When you run the JAGS program yourself, you will see that it uses many steps to produce a posterior sample for σ that has an ESS exceeding 10,000. \n",
    "* You will also see that the ESS for ν is less than 10,000 despite the long chain. \n",
    "    - In other words, there is high autocorrelation in the chains in JAGS.\n",
    "* We do not care that the chain for ν has a relatively small ESS because \n",
    "    - (a) we do not care about the exact value of ν when interpreting the posterior, as explained above, and \n",
    "    - (b) the exact value of ν has relatively little effect on the estimates of the other parameters. \n",
    "        - To be sure, the posterior sample of ν must be converged and truly representative of the posterior, but it does not need to be as finely detailed as the other parameters. \n",
    "* Nevertheless, it would be less worrying if ν had a larger ESS.\n",
    "    - The autocorrelation of the MCMC sampling in JAGS requires a long chain, \n",
    "        - which requires us to have patience while the computer chugs along. \n",
    "        - We have discussed two options for improving the efficiency of the sampling. \n",
    "            - One option is to run parallel chains on multiple cores using runjags (Section 8.7, p. 215).\n",
    "            - Another option is to implement the model in Stan, \n",
    "                - which may explore the parameter space more efficiently with its HMC sampling (Chapter 14). \n",
    "* This section shows how to run the model in Stan. Its results do indeed show a posterior sample of ν with higher ESS than JAGS."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"figures/fig16.2.png\" width=600 />"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "data {\n",
    "    int<lower=1> Ntotal ;\n",
    "    real y[Ntotal] ;\n",
    "    real meanY ;\n",
    "    real sdY ;\n",
    "}\n",
    "transformed data { // compute the constants for the priors\n",
    "    real unifLo ;\n",
    "    real unifHi ;\n",
    "    real normalSigma ;\n",
    "    real expLambda ;\n",
    "    unifLo <- sdY/1000 ;\n",
    "    unifHi <- sdY*1000 ;\n",
    "    normalSigma <- sdY*100 ;\n",
    "    expLambda <- 1/29.0 ;\n",
    "}\n",
    "parameters {\n",
    "    real<lower=0> nuMinusOne ;\n",
    "    real mu ;\n",
    "    real<lower=0> sigma ;\n",
    "}\n",
    "transformed parameters {\n",
    "    real<lower=0> nu ;\n",
    "    nu <- nuMinusOne + 1 ;\n",
    "}\n",
    "model {\n",
    "    sigma  ̃ uniform( unifLo , unifHi ) ;\n",
    "    mu  ̃ normal( meanY , normalSigma ) ;\n",
    "    nuMinusOne  ̃ exponential( expLambda ) ;\n",
    "    y  ̃ student_t( nu , mu , sigma ) ; // vectorized\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"figures/fig16.10-1.png\" width=600 />"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"figures/fig16.10-2.png\" width=600 />"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Stan-Ymet-Xnom1grp-Mrobust-Example.R"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# 책에는 없는 코드. 실습을 위해 data 폴더에 넣은 코드를 실행시키기 위한 작업\n",
    "cur_dir = getwd()\n",
    "setwd(sprintf(\"%s/%s\", cur_dir, 'data'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#------------------------------------------------------------------------------- \n",
    "# Optional generic preliminaries:\n",
    "graphics.off() # This closes all of R's graphics windows.\n",
    "rm(list=ls())  # Careful! This clears all of R's memory!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#------------------------------------------------------------------------------- \n",
    "# Load The data file \n",
    "myDataFrame = read.csv( file=\"TwoGroupIQ.csv\" )\n",
    "# For purposes of this one-group example, use data from Smart Drug group:\n",
    "myData = myDataFrame$Score[myDataFrame$Group==\"Smart Drug\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#------------------------------------------------------------------------------- \n",
    "# Load the relevant model into R's working memory:\n",
    "source(\"Stan-Ymet-Xnom1grp-Mrobust.R\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#------------------------------------------------------------------------------- \n",
    "# Optional: Specify filename root and graphical format for saving output.\n",
    "# Otherwise specify as NULL or leave saveName and saveType arguments \n",
    "# out of function calls.\n",
    "fileNameRoot = \"OneGroupIQrobust-Stan-\" \n",
    "graphFileType = \"png\" #\"eps\" "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#------------------------------------------------------------------------------- \n",
    "# Generate the MCMC chain:\n",
    "mcmcCoda = genMCMC( data=myData , numSavedSteps=20000 , saveName=fileNameRoot )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#------------------------------------------------------------------------------- \n",
    "# Display diagnostics of chain, for specified parameters:\n",
    "parameterNames = varnames(mcmcCoda) # get all parameter names\n",
    "for ( parName in parameterNames ) {\n",
    "  diagMCMC( codaObject=mcmcCoda , parName=parName , \n",
    "            saveName=fileNameRoot , saveType=graphFileType )\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# 창으로 뜨는 것들을 없앤다.\n",
    "graphics.off()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#------------------------------------------------------------------------------- \n",
    "# Get summary statistics of chain:\n",
    "summaryInfo = smryMCMC( mcmcCoda , \n",
    "                        compValMu=100.0 , ropeMu=c(99.0,101.0) ,\n",
    "                        compValSigma=15.0 , ropeSigma=c(14,16) ,\n",
    "                        compValEff=0.0 , ropeEff=c(-0.1,0.1) ,\n",
    "                        saveName=fileNameRoot )\n",
    "show(summaryInfo)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Display posterior information:\n",
    "plotMCMC( mcmcCoda , data=myData , \n",
    "          compValMu=100.0 , ropeMu=c(99.0,101.0) ,\n",
    "          compValSigma=15.0 , ropeSigma=c(14,16) ,\n",
    "          compValEff=0.0 , ropeEff=c(-0.1,0.1) ,\n",
    "          pairsPlot=TRUE , showCurve=FALSE ,\n",
    "          saveName=fileNameRoot , saveType=graphFileType )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "graphics.off()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 16.3. Two Groups"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* 16.3.1 Analysis by NHST"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### An often-used research design is a comparison of two groups. \n",
    "\n",
    "* For example, in the context of assaying the effect of a “smart drug” on IQ, instead of comparing the mean of the treatment group against an assumed landmark such as 100 (see Figure 16.3), it would make more sense to <font color=\"red\">compare the treatment group against an identically handled placebo group</font>.\n",
    "* When there are two groups, we <font color=\"red\">estimate the mean and scale for each group</font>. \n",
    "* When using t distributions for robust estimation, \n",
    "    - we could also <font color=\"blue\">estimate the normality of each group separately</font>. \n",
    "    - But because there usually are <font color=\"red\">relatively few outliers</font>, \n",
    "        - we will use a <font color=\"red\">single normality parameter to describe both groups</font>, \n",
    "            - so that the estimate of the normality is more stably estimated."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"figures/fig16.11.png\" width=600 />"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "data {\n",
    "    int<lower=1> Ntotal ;\n",
    "    int x[Ntotal] ;                                                 // group identifier     \n",
    "    real y[Ntotal] ;\n",
    "    real meanY ;\n",
    "    real sdY ;\n",
    "}\n",
    "transformed data {\n",
    "    real unifLo ;\n",
    "    real unifHi ;\n",
    "    real normalSigma ;\n",
    "    real expLambda ;\n",
    "    unifLo <- sdY/1000 ;\n",
    "    unifHi <- sdY*1000 ;\n",
    "    normalSigma <- sdY*100 ;\n",
    "    expLambda <- 1/29.0 ;\n",
    "}\n",
    "parameters {\n",
    "    real<lower=0> nuMinusOne ;\n",
    "    real mu[2] ;                                                   // 2 groups    \n",
    "    real<lower=0> sigma[2] ;                               // 2 groups\n",
    "}\n",
    "transformed parameters {\n",
    "    real<lower=0> nu ;\n",
    "    nu <- nuMinusOne + 1 ;\n",
    "}\n",
    "model {\n",
    "    sigma  ̃ uniform( unifLo , unifHi ) ;                    // vectorized 2 groups                     \n",
    "    mu  ̃ normal( meanY , normalSigma ) ;             // vectorized 2 groups\n",
    "    nuMinusOne  ̃ exponential( expLambda ) ;\n",
    "    for ( i in 1:Ntotal ) {\n",
    "        y[i]  ̃ student_t( nu , mu[x[i]] , sigma[x[i]] ) ;\n",
    "    }\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"figures/fig16.12.png\" width=600 />"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Stan-Ymet-Xnom1grp-Mrobust-Example.R"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#------------------------------------------------------------------------------- \n",
    "# Optional generic preliminaries:\n",
    "graphics.off() # This closes all of R's graphics windows.\n",
    "rm(list=ls())  # Careful! This clears all of R's memory!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#------------------------------------------------------------------------------- \n",
    "# Load The data file \n",
    "myDataFrame = read.csv( file=\"TwoGroupIQ.csv\" )\n",
    "# For purposes of this one-group example, use data from Smart Drug group:\n",
    "myData = myDataFrame$Score[myDataFrame$Group==\"Smart Drug\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#------------------------------------------------------------------------------- \n",
    "# Load the relevant model into R's working memory:\n",
    "source(\"Stan-Ymet-Xnom1grp-Mrobust.R\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#------------------------------------------------------------------------------- \n",
    "# Optional: Specify filename root and graphical format for saving output.\n",
    "# Otherwise specify as NULL or leave saveName and saveType arguments \n",
    "# out of function calls.\n",
    "fileNameRoot = \"OneGroupIQrobust-Stan-\" \n",
    "graphFileType = \"png\" #\"eps\" "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#------------------------------------------------------------------------------- \n",
    "# Generate the MCMC chain:\n",
    "mcmcCoda = genMCMC( data=myData , numSavedSteps=20000 , saveName=fileNameRoot )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#------------------------------------------------------------------------------- \n",
    "# Display diagnostics of chain, for specified parameters:\n",
    "parameterNames = varnames(mcmcCoda) # get all parameter names\n",
    "for ( parName in parameterNames ) {\n",
    "  diagMCMC( codaObject=mcmcCoda , parName=parName , \n",
    "            saveName=fileNameRoot , saveType=graphFileType )\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# 창으로 뜨는 것들을 없앤다.\n",
    "graphics.off()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#------------------------------------------------------------------------------- \n",
    "# Get summary statistics of chain:\n",
    "summaryInfo = smryMCMC( mcmcCoda , \n",
    "                        compValMu=100.0 , ropeMu=c(99.0,101.0) ,\n",
    "                        compValSigma=15.0 , ropeSigma=c(14,16) ,\n",
    "                        compValEff=0.0 , ropeEff=c(-0.1,0.1) ,\n",
    "                        saveName=fileNameRoot )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "show(summaryInfo)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "graphics.off()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Display posterior information:\n",
    "plotMCMC( mcmcCoda , data=myData , \n",
    "          compValMu=100.0 , ropeMu=c(99.0,101.0) ,\n",
    "          compValSigma=15.0 , ropeSigma=c(14,16) ,\n",
    "          compValEff=0.0 , ropeEff=c(-0.1,0.1) ,\n",
    "          pairsPlot=TRUE , showCurve=FALSE ,\n",
    "          saveName=fileNameRoot , saveType=graphFileType )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "graphics.off()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 16.3.1 Analysis by NHST"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* In traditional NHST, metric data from two groups would be submitted to a t-test."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "myDataFrame = read.csv( file=\"TwoGroupIQ.csv\" )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\n",
       "\tWelch Two Sample t-test\n",
       "\n",
       "data:  Score by Group\n",
       "t = -1.958, df = 111.44, p-value = 0.05273\n",
       "alternative hypothesis: true difference in means is not equal to 0\n",
       "95 percent confidence interval:\n",
       " -15.70602585   0.09366161\n",
       "sample estimates:\n",
       "   mean in group Placebo mean in group Smart Drug \n",
       "                100.0351                 107.8413 \n"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "t.test(Score~Group , data=myDataFrame)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* Notice that the p value is greater than 0.05, which means that the conventional decision would be not to reject the null hypothesis. \n",
    "    - This conclusion conflicts with the Bayesian analysis in Figure 16.12, unless we use a conservatively wide ROPE. \n",
    "    - The reason that \n",
    "        - the t test is less sensitive than the Bayesian estimation \n",
    "            - in this example is that \n",
    "                - the t test assumes normality and therefore \n",
    "                - its estimate of the within-group variances is too large \n",
    "                    - when there are outliers.\n",
    "    - The t test has other problems. Unlike the Bayesian analysis, \n",
    "        - the t test provides only a test of the equality of means, \n",
    "            - without a test of the equality of variances. \n",
    "            - To test equality of variances, \n",
    "                - we need to run an additional test, \n",
    "                    - namely an F test of the ratio of variances, which would inflate the p values of both tests. \n",
    "        - Moreover, \n",
    "            - both tests compute p values based on hypothetical normally distributed data, \n",
    "            - and the F test is particularly sensitive to violations of this assumption. \n",
    "    - Therefore it would be better to use resampling methods to compute the p values (and correcting them for multiple tests)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 16.4. Other Noise Distributions and Transforming Data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 16.5. EXERCISES"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 16.1."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### [Purpose: Practice using different data files in the high-level script, with an interesting real example about alcohol preference of sexually frustrated males.]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"figures/fig16.13.png\" width=600 />"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 참고자료"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* [1] Doing Bayesian Data Analysis, Second Edition: A Tutorial with R, JAGS, and Stan - http://www.amazon.com/Doing-Bayesian-Analysis-Second-Edition/dp/0124058884\n",
    "* [2] 스튜던트 t 분포(한글위키피디아) - https://ko.wikipedia.org/wiki/%EC%8A%A4%ED%8A%9C%EB%8D%98%ED%8A%B8_t_%EB%B6%84%ED%8F%AC\n",
    "* [3] t분포(티분포) 개념정리! - http://math7.tistory.com/55\n",
    "* [4] t분포표(티분포표) 보는 법! - http://math7.tistory.com/56\n",
    "* [5] 표준정규분포와 t분포의 차이 - http://blog.daum.net/_blog/BlogTypeView.do?blogid=0Isuz&articleno=1172552\n",
    "* [6] t 분포 (Student's t-distribution) - http://godrag77.blogspot.kr/2011/07/t-students-t-distribution.html"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "3.2.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}