<jupyter_start><jupyter_text>Copyright 2020 The TensorFlow Authors.Licensed under the Apache License, Version 2.0 (the "License");<jupyter_code>#@title Licensed under the Apache License, Version 2.0 (the "License"); { display-mode: "form" }
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
# https://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.<jupyter_output><empty_output><jupyter_text>Chapter 3 - Sampling the Imaginary      View on TensorFlow.org        Run in Google Colab        View source on GitHub        Download notebook    Imports and utility functions<jupyter_code>#@title Install { display-mode: "form" }
TF_Installation = 'System' #@param ['TF Nightly', 'TF Stable', 'System']

if TF_Installation == 'TF Nightly':
  !pip install -q --upgrade tf-nightly
  print('Installation of `tf-nightly` complete.')
elif TF_Installation == 'TF Stable':
  !pip install -q --upgrade tensorflow
  print('Installation of `tensorflow` complete.')
elif TF_Installation == 'System':
  pass
else:
  raise ValueError('Selection Error: Please select a valid '
                   'installation option.')
#@title Install { display-mode: "form" }
TFP_Installation = "System" #@param ["Nightly", "Stable", "System"]

if TFP_Installation == "Nightly":
  !pip install -q tfp-nightly
  print("Installation of `tfp-nightly` complete.")
elif TFP_Installation == "Stable":
  !pip install -q --upgrade tensorflow-probability
  print("Installation of `tensorflow-probability` complete.")
elif TFP_Installation == "System":
  pass
else:
  raise ValueError("Selection Error: Please select a valid "
                   "installation option.")
#@title Install { display-mode: "form" }

# Install packages that are not installed in colab
try:
  import google.colab
  IN_COLAB = True
except:
  IN_COLAB = False

if IN_COLAB:
  import sys
  if sys.executable:  # make sure pip is available
    print("Installing arviz ...")
    !{sys.executable} -m pip install -q arviz
# Core
import numpy as np
import arviz as az
import pandas as pd
import tensorflow as tf
import tensorflow_probability as tfp
import scipy.stats as stats

# visualization 
import matplotlib.pyplot as plt

# aliases
tfd = tfp.distributions
az.style.use('seaborn-colorblind')<jupyter_output><empty_output><jupyter_text>Introduction We are interested in a blood test that correctly detects vampirisim 95% of time:$$\operatorname{Pr}(\text{Positive}|\text{Vampire}) = 0.95$$The test has a false positive rate of:$$\operatorname{Pr}(\text{Positive}|\text{Mortal}) = 0.01$$We also know that vampires are rare--about 0.1% of population:$$\operatorname{Pr}(\text{Vampire}) = 0.001$$To compute $\operatorname{Pr}(\text{Vampire}|\text{Positive})$ we will apply Bayes' rule:$$\operatorname{Pr}(\text{Vampire}|\text{Positive})  = \frac{\operatorname{Pr}(\text{Positive}|\text{Vampire}) \operatorname{Pr}(\text{Vampire})}{\operatorname{Pr}(\text{Positive})}$$  Code 3.1<jupyter_code>Pr_Positive_Vampire = 0.95
Pr_Positive_Mortal = 0.01
Pr_Vampire = 0.001
tmp = Pr_Positive_Vampire * Pr_Vampire
Pr_Positive = tmp + Pr_Positive_Mortal * (1 - Pr_Vampire)
Pr_Vampire_Positive = tmp / Pr_Positive
Pr_Vampire_Positive<jupyter_output><empty_output><jupyter_text>This result shows that there is only 8.7% chance that suspect is vampire even if the test is positive, because of the low incidence rate (prior probability).  3.1 Sampling from a grid-approximate posterior  Code 3.2Let's compute the posterior for the globe tossing model, the probability of *p* conditional on the data:<jupyter_code>p_grid = tf.linspace(start=0., stop=1., num=1000)
prob_p = tf.ones(1000)
prob_data = tfd.Binomial(total_count=9, probs=p_grid).prob(6)
joint_prob = prob_data * prob_p
posterior = joint_prob / tf.reduce_sum(joint_prob)<jupyter_output><empty_output><jupyter_text>We now wish to draw 10000 samples from the posterior : Code 3.3<jupyter_code>samples = tfd.Categorical(probs=posterior).sample(10_000)<jupyter_output><empty_output><jupyter_text>Now let's display the resulting samples: Code 3.4<jupyter_code>sample_rows = p_grid.numpy()[samples]
_, ax = plt.subplots()
ax.scatter(range(len(sample_rows)), sample_rows, alpha=0.2)
ax.set(ylabel='proportion water (p)',
       xlabel='sample number',
       title='Figure 3.1 (left panel)',
       ylim=(0, 1));<jupyter_output><empty_output><jupyter_text>As the plot shows, there are many samples in the 0.6 region and very few below 0.25. Here's the density estimate of these samples, which is very similar to the ideal posterior you computed via grid approximation in Chapter 2, section 2.4.3. Code 3.5<jupyter_code>ax, = az.plot_density(sample_rows, credible_interval=1.)
ax.set(title="Figure 3.1 (right panel)",
       xlabel="proportion water (p)",
       ylabel="Density",
       xlim=(0, 1));<jupyter_output><empty_output><jupyter_text>3.2 Sampling to summarizeOnce we have the posterior distribution the next step is to summarize it. Some of the common tasks are:>* How much posterior probability lies below some parameter value?* How much posterior probability lies between two parameter values?* Which parameter value marks the lower 5% of the posterior probability?* Which range of parameter values contains 90% of the posterior probability? * Which parameter value has highest posterior probability?>These simple questions can be usefully divided into questions about (1) intervals of defined boundaries, (2) questions about intervals of defined probability mass, and (3) questions about point estimates. We’ll see how to approach these questions using samples from the posterior.  3.2.1 Intervals of defined boundariesWhat is the posterior probability that proportion of water is less than 0.5?  We can simply sum all of the probabilities where the parameter value is less than 0.5. Let's add up the posterior probability where p < 0.5: Code 3.6<jupyter_code>tf.reduce_sum(posterior[p_grid < 0.5])<jupyter_output><empty_output><jupyter_text>This means that about 17% of posterior probability is below 0.5. Now let's find the frequency of parameter values below 0.5: Code 3.7<jupyter_code>tf.where(sample_rows < 0.5).shape[0] / 10_000<jupyter_output><empty_output><jupyter_text>This is nearly the same answer as that of grid approximation, in Code 3.6, Using the same approach, how much posterior probability lies between 0.5 and 0.75? Code 3.8<jupyter_code>condition = (sample_rows > 0.5) & (sample_rows < 0.75)
tf.reduce_sum(tf.cast(condition, float))/ 10_000<jupyter_output><empty_output><jupyter_text>3.2.2 Intervals of defined massScientific journals will commonly report on an "interval of defined mass" also known as a *Confidence Interval*. The text will use the term *Compatiblity Interval* instead, since the interval indicates a range of parameter values *compatible* with the model and data. If we want to know what interval of parameter values contains 80% of the probability mass, we can simply examine the samples from the posterior, with an interval starting at *p*=0 until we have reached the 80th percentile. Code 3.9<jupyter_code>tfp.stats.percentile(sample_rows, q=80.)<jupyter_output><empty_output><jupyter_text>Similarly, the middle 80% interval lies between the 10th percentile and the 90th percentile: Code 3.10<jupyter_code>tfp.stats.percentile(sample_rows, q=[10.,90.])<jupyter_output><empty_output><jupyter_text>The text refers to intervals like this, which assign equal probability mass to each tail, as *Percentile Intervals*. They can be useful to characterize a distribution if it is fairly symmetrical. By contrast, consider a highly skewed distribution with a maximum value at *p*=1, like the posterior for observing three waters in three tosses with a uniform prior. We can compute this posterior using grid approximation: Code 3.11<jupyter_code>p_grid = tf.linspace(start=0., stop=1., num=1000)
prior = tf.ones(1000)
likelihood = tfd.Binomial(total_count=3, probs=p_grid).prob(3)
joint_prob = likelihood * prior
posterior = joint_prob / tf.reduce_sum(joint_prob)

samples = tfd.Categorical(probs=posterior).sample(10_000)
sample_rows = p_grid.numpy()[samples]<jupyter_output><empty_output><jupyter_text>Let's compute a 50% percentile compatibility interval that provides the central 50% probability by assigning 25% of the probability below the interval, and 25% above it: Code 3.12<jupyter_code>tfp.stats.percentile(sample_rows, q=[25.,75.])<jupyter_output><empty_output><jupyter_text>Given the assymetric shape of the distribution, this interval is misleading because it fails to contain the most probable parameter values at *p*=1. The *Highest Posterior Density Interval* (HPDI) is the narrowest interval containing the specified probability mass, which always contains the most probable parameter value.We can use arviz to compute it: Code 3.13<jupyter_code>az.hpd(sample_rows, credible_interval=0.5)<jupyter_output><empty_output><jupyter_text>3.2.3 Point EstimatesHow can we use the posterior distribution to create a point estimate to summarize the distribution for getting three waters from three tosses? Let's compare three different types of point estimates.  Code 3.14The *maximum a posteriori* (MAP) estimate is the parameter value with the high posterior probability.<jupyter_code>p_grid[posterior == max(posterior)]<jupyter_output><empty_output><jupyter_text>With our samples from the posterior, we can approximate the MAP: Code 3.15<jupyter_code>tf.reduce_max(sample_rows)<jupyter_output><empty_output><jupyter_text>We can also compute the posterior mean and median. Code 3.16<jupyter_code>median = tfp.stats.percentile(sample_rows, q=50.)
mean = tf.math.reduce_mean(sample_rows)

print(f"mean={mean:.6f}, median={median:.6f}")<jupyter_output>mean=0.799774, median=0.841842<jupyter_text>We can use a loss function to compute the cost of using any specific point estimate. Suppose our loss function is proportional to the difference between our decision (e.g., our point estimate) and the true value of the parameter.If we chose *p*=0.5 as our decision for the parameter value, we can use the posterior distribution to compute the expected loss, by computing the weighted average loss: We weight each loss by its corresponding posterior probability: Code 3.17<jupyter_code>tf.reduce_sum(posterior * abs(0.5 - p_grid))<jupyter_output><empty_output><jupyter_text>We can repeat this loss computation for every possible decision: Code 3.18<jupyter_code>loss = tf.map_fn(lambda d: tf.reduce_sum(posterior * np.abs(d - p_grid)), p_grid)<jupyter_output><empty_output><jupyter_text>Code 3.19Now, we can find the parameter value that minimizes the loss.<jupyter_code>p_grid[tf.math.argmin(loss)]<jupyter_output><empty_output><jupyter_text>It turns out that the value that minimizes our loss will be the same as that for the posterior median, which splits the posterior density so that half the mass is above it, with half below it.  3.3 Sampling to simulate predictionWe can use our Bayesian models to produce simulated observations, since all Bayesian models are generative.  3.3.1 Dummy data>We will call such simulated data *dummy data*, to indicate that it is a stand-in for actual data.Recall from the globe tossing model that the probability of observing *W* counts of water in *N* tosses with proportion of water *p* is given by the binomial likelihood:$$\operatorname{Pr}(W|N,p) = \frac{N!}{W!(N-W)!}p^W(1-p)^{N-W}$$For two tosses of the globe and proportion of water at 0.7 we can compute the probability of observing 0, 1, and 2 counts of water:  Code 3.20<jupyter_code>tfd.Binomial(total_count=2, probs=0.7).prob(np.arange(3))<jupyter_output><empty_output><jupyter_text>This means that there’s a 9% chance of observing *w* = 0, a 42% chance of *w* = 1, and a 49% chance of *w* = 2.We can simulate a dummy observation of *W* from our model by sampling from the binomial distribution:  Code 3.21<jupyter_code>tfd.Binomial(total_count=2, probs=0.7).sample()<jupyter_output><empty_output><jupyter_text>The result is the number of water observations in 2 tosses of the globe.A set of 10 simulations can be made by:  Code 3.22<jupyter_code>tfd.Binomial(total_count=2, probs=0.7).sample(10)<jupyter_output><empty_output><jupyter_text>When we generate 100,000 dummy observations, note that each value of *w* appears in proportion to its likelihood Code 3.23<jupyter_code>dummy_w = tfd.Binomial(total_count=2, probs=0.7).sample(100_000)

tf.unique_with_counts(tf.sort(dummy_w)).count / int(1e5)<jupyter_output><empty_output><jupyter_text>>Only two tosses of the globe isn’t much of a sample, though. So now let’s simulate the same sample size as before, 9 tosses. Code 3.24<jupyter_code>dummy_w = tfd.Binomial(total_count=9, probs=0.7).sample((100000,))
_, ax = plt.subplots()
ax.hist(dummy_w.numpy(), bins=np.arange(11), rwidth=0.2)
ax.set(xlabel="dummy water count", 
       ylabel="Frequency");<jupyter_output><empty_output><jupyter_text>3.3.2. Model Checking>We’d like to propagate the parameter uncertainty—carry it forward—as we evaluate the implied predictions. All that is required is averaging over the posterior density for *p*, while computing the predictions. For each possible value of the parameter *p*, there is an implied distribution of outcomes. So if you were to compute the sampling distribution of outcomes at each value of *p*, then you could average all of these prediction distributions together, using the posterior probabilities of each value of *p*, to get a POSTERIOR PREDICTIVE DISTRIBUTION.To simulate predicted observations for nine globe tosses, for a single value of *p*=0.6, we can use `Binomial` to generate random binomial samples:  Code 3.25<jupyter_code>w = tfd.Binomial(total_count=9, probs=0.6).sample(1e4)<jupyter_output><empty_output><jupyter_text>>All you need to propagate parameter uncertainty into these predictions is replace the value 0.6 with samples from the posterior.... Since the sampled values appear in proportion to their posterior probabilities, the resulting simulated observations are averaged over the posterior. Code 3.26<jupyter_code>w = tfd.Binomial(total_count=9, probs=sample_rows).sample()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.0929573571083005, 'fraction_numerical': 0.02467343976777939, 'mean_word_length': 5.689782709200185, 'pattern_counts': {'":': 3, '<': 116, '<?xml version=': 0, '>': 119, 'https://': 1, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '28', 'dir': 'jupyter-structured-clean-dedup', 'id': '365709', 'license': 'Apache-2.0', 'n_tokens_mistral': 4637, 'n_tokens_neox': 4283, 'n_words': 1579, 'path': 'tensorflow_probability/examples/statistical_rethinking/notebooks/03_sampling_the_imaginary.ipynb', 'repo_name': 'Leengit/probability', 'size': '14469'}