text:
```python
import astropy
from astropy import units as u
from astropy import cosmology as cosmo
```

```python
from astropy.cosmology import WMAP5
WMAP5.H(0)
```

$70.2 \; \mathrm{\frac{km}{Mpc\,s}}$

```python
WMAP5.critical_density(0)# critical density at z = 0 in g/cm^3
```

$9.2565426 \times 10^{-30} \; \mathrm{\frac{g}{cm^{3}}}$

```python
WMAP5.H(7)  # Hubble parameter at z = 7 in km/s/Mpc
```

$839.15905 \; \mathrm{\frac{km}{Mpc\,s}}$

```python
from astropy.cosmology import FlatLambdaCDM
import astropy.units as u
cosmo = FlatLambdaCDM(H0=70.2 * u.km / u.s / u.Mpc, Om0=0.3)
```

```python
cosmo
```

    FlatLambdaCDM(H0=70.2 km / (Mpc s), Om0=0.3, Tcmb0=0 K, Neff=3.04, m_nu=None, Ob0=None)

```python
FlatLambdaCDM(H0=70.2, Om0=0.3, Tcmb0=2.725 ,
              Neff=3.04, m_nu=[0], Ob0=None)
```

    FlatLambdaCDM(H0=70.2 km / (Mpc s), Om0=0.3, Tcmb0=2.725 K, Neff=3.04, m_nu=[0. 0. 0.] eV, Ob0=None)

```python
from astropy.cosmology import WMAP5   # WMAP 5-year cosmology
WMAP5.Om([0, 1.0, 2.0]), WMAP5.Ode([0., 1.0, 2.0])
```

    (array([0.277     , 0.7536719 , 0.91109749]),
     array([0.72291523, 0.24586683, 0.08806609]))

```python
from astropy.cosmology import WMAP5
print(WMAP5.__doc__)
```

    WMAP5 instance of FlatLambdaCDM cosmology
    (from Komatsu et al. 2009, ApJS, 180, 330, doi: 10.1088/0067-0049/180/2/330. Table 1 (WMAP + BAO + SN ML).)

```python

```

```python
import pysm

import numpy as np
import healpy as hp
import matplotlib.pyplot as plt
%matplotlib inline
```

```python
from pysm.nominal import models
```

```python
import os
import json
import numpy as np
from scipy.spatial import cKDTree
from astropy import units as u
from astropy.coordinates import Angle, Distance, ICRSCoordinates
from astropy.cosmology import Planck13, WMAP5
from vpf_calc import vpf
from dat.cartesian_cosmo import mk_coords
from dat.h5_funcs import h5_arr, arr2h5

def mk_mock_srch(radecfile, nzdictfile, Nsph, simul_cosmo):

    if simul_cosmo == "Planck":
        # First make h free
        Planck13.__init__(100.0, Planck13.Om0)
        cosmo = Planck13
    elif simul_cosmo == "WMAP":
        WMAP5.__init__(100.0, WMAP5.Om0)
        cosmo = WMAP5
    comv = cosmo.comoving_distance

    radecarr = h5_arr(radecfile, "good_pts")
    nzdict = json.load(open(nzdictfile))

    Nrands = radecarr.shape[0]
    Narrs = Nsph / Nrands
    remain = Nsph % Nrands

    radecz = np.zeros((Nsph, 3))

    for i in range(Narrs):

        start = Nrands * i
        stop = Nrands * (i + 1)
        radecz[start:stop, :2] = radecarr[:, :]

    endchunk = Nrands * (Narrs)
    radecz[endchunk:, :2] = radecarr[:remain, :]

    rad = np.arange(1.0, 67.0, 5.0)
    zlo = nzdict["zlo"]
    zhi = nzdict["zhi"]

    radeczlist = len(rad) * [radecz]

    for r_i, r in enumerate(rad):

        dis_near = Distance(comv(zlo).value + r, u.Mpc)
        dis_far = Distance(comv(zhi).value - r, u.Mpc)

        z_a = dis_near.compute_z(cosmology=cosmo)

        z_b = dis_far.compute_z(cosmology=cosmo)

        randz = (z_a ** 3 + \
                 (z_b ** 3 - z_a ** 3) * np.random.rand(Nsph)) ** (1. / 3.)

        radeczlist[r_i][:, 2] = randz[:]

        arr2h5(radeczlist[r_i], "{0}/{1}/mocks/mock_srch_pts.hdf5".format(os.path.dirname(radecfile), simul_cosmo), "radecz_{0}".format(str(r_i * 5 + 1)))

def mk_mock_coords(radeczfile, outfile, simul_cosmo):

    if simul_cosmo == "Planck":
        Planck13.__init__(100.0, Planck13.Om0)
        cosmo = Planck13
    elif simul_cosmo == "WMAP":
        WMAP5.__init__(100.0, WMAP5.Om0)
        cosmo = WMAP5

    rad = np.arange(1.0, 67.0, 5.0)

    radecz = h5_arr(radeczfile, "radecz")

    cart = np.zeros(radecz.shape)

    for i, rdz in enumerate(radecz):

        ra = Angle(rdz[0], u.deg)
        dec = Angle(rdz[1], u.deg)

        losd = cosmo.comoving_distance(rdz[2])
        dis = Distance(losd)

        coord = ICRSCoordinates(ra, dec, distance=dis)

        cart[i, :] = np.array([coord.x.value, coord.y.value, coord.z.value])

    np.savetxt(outfile, cart)

def mock_vpf(mock_cart_coords, spheresfile, simul_cosmo, rad):

    gals = h5_arr(mock_cart_coords, "coords")

    print gals

    name = mock_cart_coords.split("/")[-1].split(".")[0]

    if simul_cosmo == "Planck":
        # First make h free
        Planck13.__init__(100.0, Planck13.Om0)
        cosmo = Planck13
    elif simul_cosmo == "WMAP":
        WMAP5.__init__(100.0, WMAP5.Om0)
        cosmo = WMAP5
    comv = cosmo.comoving_distance

    gal_baum = cKDTree(gals)

    spheres = h5_arr(spheresfile, "radecz_{0}".format(str(int(rad))))

    print spheres

    for i, sphere in enumerate(spheres):

        rang = Angle(sphere[0], u.deg)
        decang = Angle(sphere[1], u.deg)

        dis = Distance(comv(sphere[2]), u.Mpc)

        coord = ICRSCoordinates(rang, decang, distance=dis)

        sph_cen = np.array([coord.x, coord.y, coord.z])

        nn = gal_baum.query(sph_cen)

        print "rad: ", rad, ", sphere: ", i

        f = open("{0}/vpf_out/ascii/{1}_{2}.dat".format(os.path.dirname(spheresfile), name, str(int(rad))), 'a')

        if not nn[0] < rad:
            f.write("1\n")
        else:
            f.write("0\n")

        f.close()
```

```python
from nbodykit.lab import cosmology

cosmo = cosmology.Cosmology()
print("original sigma8 = %.4f" % cosmo.sigma8)

new_cosmo = cosmo.match(sigma8=0.82)
print("new sigma8 = %.4f" % new_cosmo.sigma8)
```

```python
$ sudo port install autoconf automake libtool
```

```python
install autoconf automake libtool
```

```python
from nbodykit.lab import *
from nbodykit import setup_logging
```

```python
import nbodykit; print(nbodykit, nbodykit.__version__);
```

```python
pip install camb
```

    Collecting camb
      Using cached camb-1.3.5.tar.gz (1.7 MB)
    Requirement already satisfied: scipy>=1.0 in /opt/anaconda3/lib/python3.8/site-packages (from camb) (1.8.0)
    Requirement already satisfied: sympy>=1.0 in /opt/anaconda3/lib/python3.8/site-packages (from camb) (1.8)
    Requirement already satisfied: numpy<1.25.0,>=1.17.3 in /opt/anaconda3/lib/python3.8/site-packages (from scipy>=1.0->camb) (1.22.2)
    Requirement already satisfied: mpmath>=0.19 in /opt/anaconda3/lib/python3.8/site-packages (from sympy>=1.0->camb) (1.2.1)
    Building wheels for collected packages: camb
      Building wheel for camb (setup.py) ... error
      ERROR: Command errored out with exit status 1:
       command: /opt/anaconda3/bin/python -u -c 'import sys, setuptools, tokenize; sys.argv[0] = '"'"'/private/var/folders/pr/2pb42ldd0dzc3mwbz3chdrpc0000gn/T/pip-install-z3bc_mrp/camb_a39c83dcdb3e448ba60a0623d907c9aa/setup.py'"'"'; __file__='"'"'/private/var/folders/pr/2pb42ldd0dzc3mwbz3chdrpc0000gn/T/pip-install-z3bc_mrp/camb_a39c83dcdb3e448ba60a0623d907c9aa/setup.py'"'"';f=getattr(tokenize, '"'"'open'"'"', open)(__file__);code=f.read().replace('"'"'\r\n'"'"', '"'"'\n'"'"');f.close();exec(compile(code, __file__, '"'"'exec'"'"'))' bdist_wheel -d /private/var/folders/pr/2pb42ldd0dzc3mwbz3chdrpc0000gn/T/pip-wheel-_dnc0ofq
           cwd: /private/var/folders/pr/2pb42ldd0dzc3mwbz3chdrpc0000gn/T/pip-install-z3bc_mrp/camb_a39c83dcdb3e448ba60a0623d907c9aa/
      Complete output (35 lines):
      running bdist_wheel
      running build
      running build_py
      Traceback (most recent call last):
        File "<string>", line 1, in <module>
        File "/private/var/folders/pr/2pb42ldd0dzc3mwbz3chdrpc0000gn/T/pip-install-z3bc_mrp/camb_a39c83dcdb3e448ba60a0623d907c9aa/setup.py", line 268, in <module>
          setup(name=os.getenv('CAMB_PACKAGE_NAME', 'camb'),
        File "/opt/anaconda3/lib/python3.8/site-packages/setuptools/__init__.py", line 153, in setup
          return distutils.core.setup(**attrs)
        File "/opt/anaconda3/lib/python3.8/distutils/core.py", line 148, in setup
          dist.run_commands()
        File "/opt/anaconda3/lib/python3.8/distutils/dist.py", line 966, in run_commands
          self.run_command(cmd)
        File "/opt/anaconda3/lib/python3.8/distutils/dist.py", line 985, in run_command
          cmd_obj.run()
        File "/opt/anaconda3/lib/python3.8/site-packages/wheel/bdist_wheel.py", line 299, in run
          self.run_command('build')
        File "/opt/anaconda3/lib/python3.8/distutils/cmd.py", line 313, in run_command
          self.distribution.run_command(command)
        File "/opt/anaconda3/lib/python3.8/distutils/dist.py", line 985, in run_command
          cmd_obj.run()
        File "/opt/anaconda3/lib/python3.8/distutils/command/build.py", line 135, in run
          self.run_command(cmd_name)
        File "/opt/anaconda3/lib/python3.8/distutils/cmd.py", line 313, in run_command
          self.distribution.run_command(command)
        File "/opt/anaconda3/lib/python3.8/distutils/dist.py", line 985, in run_command
          cmd_obj.run()
        File "/private/var/folders/pr/2pb42ldd0dzc3mwbz3chdrpc0000gn/T/pip-install-z3bc_mrp/camb_a39c83dcdb3e448ba60a0623d907c9aa/setup.py", line 232, in run
          make_library(False)
        File "/private/var/folders/pr/2pb42ldd0dzc3mwbz3chdrpc0000gn/T/pip-install-z3bc_mrp/camb_a39c83dcdb3e448ba60a0623d907c9aa/setup.py", line 111, in make_library
          ok, gfortran_version = _compile.check_gfortran(msg=not _compile.is_windows)
        File "/private/var/folders/pr/2pb42ldd0dzc3mwbz3chdrpc0000gn/T/pip-install-z3bc_mrp/camb_a39c83dcdb3e448ba60a0623d907c9aa/camb/_compilers.py", line 71, in check_gfortran
          raise Exception(
      Exception: You need ifort or gfortran 6 or higher to compile (found: None).
      See https://camb.readthedocs.io/en/latest/fortran_compilers.html
      ----------------------------------------
      ERROR: Failed building wheel for camb
      Running setup.py clean for camb
    Failed to build camb
    Installing collected packages: camb
        Running setup.py install for camb ... error
        ERROR: Command errored out with exit status 1:
         command: /opt/anaconda3/bin/python -u -c 'import sys, setuptools, tokenize; sys.argv[0] = '"'"'/private/var/folders/pr/2pb42ldd0dzc3mwbz3chdrpc0000gn/T/pip-install-z3bc_mrp/camb_a39c83dcdb3e448ba60a0623d907c9aa/setup.py'"'"'; __file__='"'"'/private/var/folders/pr/2pb42ldd0dzc3mwbz3chdrpc0000gn/T/pip-install-z3bc_mrp/camb_a39c83dcdb3e448ba60a0623d907c9aa/setup.py'"'"';f=getattr(tokenize, '"'"'open'"'"', open)(__file__);code=f.read().replace('"'"'\r\n'"'"', '"'"'\n'"'"');f.close();exec(compile(code, __file__, '"'"'exec'"'"'))' install --record /private/var/folders/pr/2pb42ldd0dzc3mwbz3chdrpc0000gn/T/pip-record-i3yzpu7m/install-record.txt --single-version-externally-managed --compile --install-headers /opt/anaconda3/include/python3.8/camb
             cwd: /private/var/folders/pr/2pb42ldd0dzc3mwbz3chdrpc0000gn/T/pip-install-z3bc_mrp/camb_a39c83dcdb3e448ba60a0623d907c9aa/
        Complete output (37 lines):
        running install
        running build
        running build_py
        Traceback (most recent call last):
          File "<string>", line 1, in <module>
          File "/private/var/folders/pr/2pb42ldd0dzc3mwbz3chdrpc0000gn/T/pip-install-z3bc_mrp/camb_a39c83dcdb3e448ba60a0623d907c9aa/setup.py", line 268, in <module>
            setup(name=os.getenv('CAMB_PACKAGE_NAME', 'camb'),
          File "/opt/anaconda3/lib/python3.8/site-packages/setuptools/__init__.py", line 153, in setup
            return distutils.core.setup(**attrs)
          File "/opt/anaconda3/lib/python3.8/distutils/core.py", line 148, in setup
            dist.run_commands()
          File "/opt/anaconda3/lib/python3.8/distutils/dist.py", line 966, in run_commands
            self.run_command(cmd)
          File "/opt/anaconda3/lib/python3.8/distutils/dist.py", line 985, in run_command
            cmd_obj.run()
          File "/opt/anaconda3/lib/python3.8/site-packages/setuptools/command/install.py", line 61, in run
            return orig.install.run(self)
          File "/opt/anaconda3/lib/python3.8/distutils/command/install.py", line 545, in run
            self.run_command('build')
          File "/opt/anaconda3/lib/python3.8/distutils/cmd.py", line 313, in run_command
            self.distribution.run_command(command)
          File "/opt/anaconda3/lib/python3.8/distutils/dist.py", line 985, in run_command
            cmd_obj.run()
          File "/opt/anaconda3/lib/python3.8/distutils/command/build.py", line 135, in run
            self.run_command(cmd_name)
          File "/opt/anaconda3/lib/python3.8/distutils/cmd.py", line 313, in run_command
            self.distribution.run_command(command)
          File "/opt/anaconda3/lib/python3.8/distutils/dist.py", line 985, in run_command
            cmd_obj.run()
          File "/private/var/folders/pr/2pb42ldd0dzc3mwbz3chdrpc0000gn/T/pip-install-z3bc_mrp/camb_a39c83dcdb3e448ba60a0623d907c9aa/setup.py", line 232, in run
            make_library(False)
          File "/private/var/folders/pr/2pb42ldd0dzc3mwbz3chdrpc0000gn/T/pip-install-z3bc_mrp/camb_a39c83dcdb3e448ba60a0623d907c9aa/setup.py", line 111, in make_library
            ok, gfortran_version = _compile.check_gfortran(msg=not _compile.is_windows)
          File "/private/var/folders/pr/2pb42ldd0dzc3mwbz3chdrpc0000gn/T/pip-install-z3bc_mrp/camb_a39c83dcdb3e448ba60a0623d907c9aa/camb/_compilers.py", line 71, in check_gfortran
            raise Exception(
        Exception: You need ifort or gfortran 6 or higher to compile (found: None).
        See https://camb.readthedocs.io/en/latest/fortran_compilers.html
        ----------------------------------------
    ERROR: Command errored out with exit status 1: /opt/anaconda3/bin/python -u -c 'import sys, setuptools, tokenize; sys.argv[0] = '"'"'/private/var/folders/pr/2pb42ldd0dzc3mwbz3chdrpc0000gn/T/pip-install-z3bc_mrp/camb_a39c83dcdb3e448ba60a0623d907c9aa/setup.py'"'"'; __file__='"'"'/private/var/folders/pr/2pb42ldd0dzc3mwbz3chdrpc0000gn/T/pip-install-z3bc_mrp/camb_a39c83dcdb3e448ba60a0623d907c9aa/setup.py'"'"';f=getattr(tokenize, '"'"'open'"'"', open)(__file__);code=f.read().replace('"'"'\r\n'"'"', '"'"'\n'"'"');f.close();exec(compile(code, __file__, '"'"'exec'"'"'))' install --record /private/var/folders/pr/2pb42ldd0dzc3mwbz3chdrpc0000gn/T/pip-record-i3yzpu7m/install-record.txt --single-version-externally-managed --compile --install-headers /opt/anaconda3/include/python3.8/camb Check the logs for full command output.
    Note: you may need to restart the kernel to use updated packages.

```python

```

max_stars_repo_name: DhruvKumarPHY/solutions

max_issues_repo_name: DhruvKumarPHY/solutions

max_forks_repo_name: DhruvKumarPHY/solutions


text:
Probabilistic Programming
=====
and Bayesian Methods for Hackers
========

##### Version 0.1

`Original content created by Cam Davidson-Pilon`

`Ported to Python 3 and PyMC3 by Max Margenot (@clean_utensils) and Thomas Wiecki (@twiecki) at Quantopian (@quantopian)`
___

Welcome to *Bayesian Methods for Hackers*. The full Github repository is available at [github/Probabilistic-Programming-and-Bayesian-Methods-for-Hackers](https://github.com/CamDavidsonPilon/Probabilistic-Programming-and-Bayesian-Methods-for-Hackers). The other chapters can be found on the project's [homepage](https://camdavidsonpilon.github.io/Probabilistic-Programming-and-Bayesian-Methods-for-Hackers/). We hope you enjoy the book, and we encourage any contributions!

Chapter 1
======
***

The Philosophy of Bayesian Inference
------

> You are a skilled programmer, but bugs still slip into your code. After a particularly difficult implementation of an algorithm, you decide to test your code on a trivial example. It passes. You test the code on a harder problem. It passes once again. And it passes the next, *even more difficult*, test too! You are starting to believe that there may be no bugs in this code...

If you think this way, then congratulations, you already are thinking Bayesian! Bayesian inference is simply updating your beliefs after considering new evidence. A Bayesian can rarely be certain about a result, but he or she can be very confident. Just like in the example above, we can never be 100% sure that our code is bug-free unless we test it on every possible problem; something rarely possible in practice. Instead, we can test it on a large number of problems, and if it succeeds we can feel more *confident* about our code, but still not certain.  Bayesian inference works identically: we update our beliefs about an outcome; rarely can we be absolutely sure unless we rule out all other alternatives.

### The Bayesian state of mind

Bayesian inference differs from more traditional statistical inference by preserving *uncertainty*. At first, this sounds like a bad statistical technique. Isn't statistics all about deriving *certainty* from randomness? To reconcile this, we need to start thinking like Bayesians.

The Bayesian world-view interprets probability as measure of *believability in an event*, that is, how confident we are in an event occurring. In fact, we will see in a moment that this is the natural interpretation of probability.

For this to be clearer, we consider an alternative interpretation of probability: *Frequentist*, known as the more *classical* version of statistics, assume that probability is the long-run frequency of events (hence the bestowed title). For example, the *probability of plane accidents* under a frequentist philosophy is interpreted as the *long-term frequency of plane accidents*. This makes logical sense for many probabilities of events, but becomes more difficult to understand when events have no long-term frequency of occurrences. Consider: we often assign probabilities to outcomes of presidential elections, but the election itself only happens once! Frequentists get around this by invoking alternative realities and saying across all these realities, the frequency of occurrences defines the probability.

Bayesians, on the other hand, have a more intuitive approach. Bayesians interpret a probability as measure of *belief*, or confidence, of an event occurring. Simply, a probability is a summary of an opinion. An individual who assigns a belief of 0 to an event has no confidence that the event will occur; conversely, assigning a belief of 1 implies that the individual is absolutely certain of an event occurring. Beliefs between 0 and 1 allow for weightings of other outcomes. This definition agrees with the probability of a plane accident example, for having observed the frequency of plane accidents, an individual's belief should be equal to that frequency, excluding any outside information. Similarly, under this definition of probability being equal to beliefs, it is meaningful to speak about probabilities (beliefs) of presidential election outcomes: how confident are you candidate *A* will win?

Notice in the paragraph above, I assigned the belief (probability) measure to an *individual*, not to Nature. This is very interesting, as this definition leaves room for conflicting beliefs between individuals. Again, this is appropriate for what naturally occurs: different individuals have different beliefs of events occurring, because they possess different *information* about the world. The existence of different beliefs does not imply that anyone is wrong. Consider the following examples demonstrating the relationship between individual beliefs and probabilities:

- I flip a coin, and we both guess the result. We would both agree, assuming the coin is fair, that the probability of Heads is 1/2. Assume, then, that I peek at the coin. Now I know for certain what the result is: I assign probability 1.0 to either Heads or Tails (whichever it is). Now what is *your* belief that the coin is Heads? My knowledge of the outcome has not changed the coin's results. Thus we assign different probabilities to the result.

-  Your code either has a bug in it or not, but we do not know for certain which is true, though we have a belief about the presence or absence of a bug.

-  A medical patient is exhibiting symptoms $x$, $y$ and $z$. There are a number of diseases that could be causing all of them, but only a single disease is present. A doctor has beliefs about which disease, but a second doctor may have slightly different beliefs.

This philosophy of treating beliefs as probability is natural to humans. We employ it constantly as we interact with the world and only see partial truths, but gather evidence to form beliefs. Alternatively, you have to be *trained* to think like a frequentist.

To align ourselves with traditional probability notation, we denote our belief about event $A$ as $P(A)$. We call this quantity the *prior probability*.

John Maynard Keynes, a great economist and thinker, said "When the facts change, I change my mind. What do you do, sir?" This quote reflects the way a Bayesian updates his or her beliefs after seeing evidence. Even — especially — if the evidence is counter to what was initially believed, the evidence cannot be ignored. We denote our updated belief as $P(A |X )$, interpreted as the probability of $A$ given the evidence $X$. We call the updated belief the *posterior probability* so as to contrast it with the prior probability. For example, consider the posterior probabilities (read: posterior beliefs) of the above examples, after observing some evidence $X$:

1\. $P(A): \;\;$ the coin has a 50 percent chance of being Heads. $P(A | X):\;\;$ You look at the coin, observe a Heads has landed, denote this information $X$, and trivially assign probability 1.0 to Heads and 0.0 to Tails.

2\.   $P(A): \;\;$  This big, complex code likely has a bug in it. $P(A | X): \;\;$ The code passed all $X$ tests; there still might be a bug, but its presence is less likely now.

3\.  $P(A):\;\;$ The patient could have any number of diseases. $P(A | X):\;\;$ Performing a blood test generated evidence $X$, ruling out some of the possible diseases from consideration.

It's clear that in each example we did not completely discard the prior belief after seeing new evidence $X$, but we *re-weighted the prior* to incorporate the new evidence (i.e. we put more weight, or confidence, on some beliefs versus others).

By introducing prior uncertainty about events, we are already admitting that any guess we make is potentially very wrong. After observing data, evidence, or other information, we update our beliefs, and our guess becomes *less wrong*. This is the alternative side of the prediction coin, where typically we try to be *more right*.

### Bayesian Inference in Practice

 If frequentist and Bayesian inference were programming functions, with inputs being statistical problems, then the two would be different in what they return to the user. The frequentist inference function would return a number, representing an estimate (typically a summary statistic like the sample average etc.), whereas the Bayesian function would return *probabilities*.

For example, in our debugging problem above, calling the frequentist function with the argument "My code passed all $X$ tests; is my code bug-free?" would return a *YES*. On the other hand, asking our Bayesian function "Often my code has bugs. My code passed all $X$ tests; is my code bug-free?" would return something very different: probabilities of *YES* and *NO*. The function might return:

>    *YES*, with probability 0.8; *NO*, with probability 0.2

This is very different from the answer the frequentist function returned. Notice that the Bayesian function accepted an additional argument:  *"Often my code has bugs"*. This parameter is the *prior*. By including the prior parameter, we are telling the Bayesian function to include our belief about the situation. Technically this parameter in the Bayesian function is optional, but we will see excluding it has its own consequences.

#### Incorporating evidence

As we acquire more and more instances of evidence, our prior belief is *washed out* by the new evidence. This is to be expected. For example, if your prior belief is something ridiculous, like "I expect the sun to explode today", and each day you are proved wrong, you would hope that any inference would correct you, or at least align your beliefs better. Bayesian inference will correct this belief.

Denote $N$ as the number of instances of evidence we possess. As we gather an *infinite* amount of evidence, say as $N \rightarrow \infty$, our Bayesian results (often) align with frequentist results. Hence for large $N$, statistical inference is more or less objective. On the other hand, for small $N$, inference is much more *unstable*: frequentist estimates have more variance and larger confidence intervals. This is where Bayesian analysis excels. By introducing a prior, and returning probabilities (instead of a scalar estimate), we *preserve the uncertainty* that reflects the instability of statistical inference of a small $N$ dataset.

One may think that for large $N$, one can be indifferent between the two techniques since they offer similar inference, and might lean towards the computationally-simpler, frequentist methods. An individual in this position should consider the following quote by Andrew Gelman (2005)[1], before making such a decision:

> Sample sizes are never large. If $N$ is too small to get a sufficiently-precise estimate, you need to get more data (or make more assumptions). But once $N$ is "large enough," you can start subdividing the data to learn more (for example, in a public opinion poll, once you have a good estimate for the entire country, you can estimate among men and women, northerners and southerners, different age groups, etc.). $N$ is never enough because if it were "enough" you'd already be on to the next problem for which you need more data.

### Are frequentist methods incorrect then?

**No.**

Frequentist methods are still useful or state-of-the-art in many areas. Tools such as least squares linear regression, LASSO regression, and expectation-maximization algorithms are all powerful and fast. Bayesian methods complement these techniques by solving problems that these approaches cannot, or by illuminating the underlying system with more flexible modeling.

#### A note on *Big Data*
Paradoxically, big data's predictive analytic problems are actually solved by relatively simple algorithms [2][4]. Thus we can argue that big data's prediction difficulty does not lie in the algorithm used, but instead on the computational difficulties of storage and execution on big data. (One should also consider Gelman's quote from above and ask "Do I really have big data?")

The much more difficult analytic problems involve *medium data* and, especially troublesome, *really small data*. Using a similar argument as  Gelman's above, if big data problems are *big enough* to be readily solved, then we should be more interested in the *not-quite-big enough* datasets.

### Our Bayesian framework

We are interested in beliefs, which can be interpreted as probabilities by thinking Bayesian. We have a *prior* belief in event $A$, beliefs formed by previous information, e.g., our prior belief about bugs being in our code before performing tests.

Secondly, we observe our evidence. To continue our buggy-code example: if our code passes $X$ tests, we want to update our belief to incorporate this. We call this new belief the *posterior* probability. Updating our belief is done via the following equation, known as Bayes' Theorem, after its discoverer Thomas Bayes:

\begin{align}
 P( A | X ) = & \frac{ P(X | A) P(A) } {P(X) } \\\\[5pt]
& \propto P(X | A) P(A)\;\; (\propto \text{is proportional to })
\end{align}

The above formula is not unique to Bayesian inference: it is a mathematical fact with uses outside Bayesian inference. Bayesian inference merely uses it to connect prior probabilities $P(A)$ with an updated posterior probabilities $P(A | X )$.

##### Example: Mandatory coin-flip example

Every statistics text must contain a coin-flipping example, I'll use it here to get it out of the way. Suppose, naively, that you are unsure about the probability of heads in a coin flip (spoiler alert: it's 50%). You believe there is some true underlying ratio, call it $p$, but have no prior opinion on what $p$ might be.

We begin to flip a coin, and record the observations: either $H$ or $T$. This is our observed data. An interesting question to ask is how our inference changes as we observe more and more data? More specifically, what do our posterior probabilities look like when we have little data, versus when we have lots of data.

Below we plot a sequence of updating posterior probabilities as we observe increasing amounts of data (coin flips).

```python
"""
The book uses a custom matplotlibrc file, which provides the unique styles for
matplotlib plots. If executing this book, and you wish to use the book's
styling, provided are two options:
    1. Overwrite your own matplotlibrc file with the rc-file provided in the
       book's styles/ dir. See http://matplotlib.org/users/customizing.html
    2. Also in the styles is  bmh_matplotlibrc.json file. This can be used to
       update the styles in only this notebook. Try running the following code:

        import json
        s = json.load(open("../styles/bmh_matplotlibrc.json"))
        matplotlib.rcParams.update(s)

"""

# The code below can be passed over, as it is currently not important, plus it
# uses advanced topics we have not covered yet. LOOK AT PICTURE, MICHAEL!
%matplotlib inline
from IPython.core.pylabtools import figsize
import numpy as np
from matplotlib import pyplot as plt
figsize(11, 9)

import scipy.stats as stats

dist = stats.beta
n_trials = [0, 1, 2, 3, 4, 5, 8, 15, 50, 500]
data = stats.bernoulli.rvs(0.5, size=n_trials[-1])
x = np.linspace(0, 1, 100)

# For the already prepared, I'm using Binomial's conj. prior.
for k, N in enumerate(n_trials):
    sx = plt.subplot(len(n_trials)/2, 2, k+1)
    plt.xlabel("$p$, probability of heads") \
        if k in [0, len(n_trials)-1] else None
    plt.setp(sx.get_yticklabels(), visible=False)
    heads = data[:N].sum()
    y = dist.pdf(x, 1 + heads, 1 + N - heads)
    plt.plot(x, y, label="observe %d tosses,\n %d heads" % (N, heads))
    plt.fill_between(x, 0, y, color="#348ABD", alpha=0.4)
    plt.vlines(0.5, 0, 4, color="k", linestyles="--", lw=1)

    leg = plt.legend()
    leg.get_frame().set_alpha(0.4)
    plt.autoscale(tight=True)

plt.suptitle("Bayesian updating of posterior probabilities",
             y=1.02,
             fontsize=14)

plt.tight_layout()
```

The posterior probabilities are represented by the curves, and our uncertainty is proportional to the width of the curve. As the plot above shows, as we start to observe data our posterior probabilities start to shift and move around. Eventually, as we observe more and more data (coin-flips), our probabilities will tighten closer and closer around the true value of $p=0.5$ (marked by a dashed line).

Notice that the plots are not always *peaked* at 0.5. There is no reason it should be: recall we assumed we did not have a prior opinion of what $p$ is. In fact, if we observe quite extreme data, say 8 flips and only 1 observed heads, our distribution would look very biased *away* from lumping around 0.5 (with no prior opinion, how confident would you feel betting on a fair coin after observing 8 tails and 1 head?). As more data accumulates, we would see more and more probability being assigned at $p=0.5$, though never all of it.

The next example is a simple demonstration of the mathematics of Bayesian inference.

##### Example: Bug, or just sweet, unintended feature?

Let $A$ denote the event that our code has **no bugs** in it. Let $X$ denote the event that the code passes all debugging tests. For now, we will leave the prior probability of no bugs as a variable, i.e. $P(A) = p$.

We are interested in $P(A|X)$, i.e. the probability of no bugs, given our debugging tests $X$. To use the formula above, we need to compute some quantities.

What is $P(X | A)$, i.e., the probability that the code passes $X$ tests *given* there are no bugs? Well, it is equal to 1, for a code with no bugs will pass all tests.

$P(X)$ is a little bit trickier: The event $X$ can be divided into two possibilities, event $X$ occurring even though our code *indeed has* bugs (denoted $\sim A\;$, spoken *not $A$*), or event $X$ without bugs ($A$). $P(X)$ can be represented as:

\begin{align}
P(X ) & = P(X \text{ and } A) + P(X \text{ and } \sim A) \\\\[5pt]
 & = P(X|A)P(A) + P(X | \sim A)P(\sim A)\\\\[5pt]
& = P(X|A)p + P(X | \sim A)(1-p)
\end{align}

We have already computed $P(X|A)$ above. On the other hand, $P(X | \sim A)$ is subjective: our code can pass tests but still have a bug in it, though the probability there is a bug present is reduced. Note this is dependent on the number of tests performed, the degree of complication in the tests, etc. Let's be conservative and assign $P(X|\sim A) = 0.5$. Then

\begin{align}
P(A | X) & = \frac{1\cdot p}{ 1\cdot p +0.5 (1-p) } \\\\
& = \frac{ 2 p}{1+p}
\end{align}
This is the posterior probability. What does it look like as a function of our prior, $p \in [0,1]$?

```python
figsize(12.5, 4)
p = np.linspace(0, 1, 50)
plt.plot(p, 2*p/(1+p), color="#348ABD", lw=3)
#plt.fill_between(p, 2*p/(1+p), alpha=.5, facecolor=["#A60628"])
plt.scatter(0.2, 2*(0.2)/1.2, s=140, c="#348ABD")
plt.xlim(0, 1)
plt.ylim(0, 1)
plt.xlabel("Prior, $P(A) = p$")
plt.ylabel("Posterior, $P(A|X)$, with $P(A) = p$")
plt.title("Are there bugs in my code?");
```

We can see the biggest gains if we observe the $X$ tests passed when the prior probability, $p$, is low. Let's settle on a specific value for the prior. I'm a strong programmer (I think), so I'm going to give myself a realistic prior of 0.20, that is, there is a 20% chance that I write code bug-free. To be more realistic, this prior should be a function of how complicated and large the code is, but let's pin it at 0.20. Then my updated belief that my code is bug-free is 0.33.

Recall that the prior is a probability: $p$ is the prior probability that there *are no bugs*, so $1-p$ is the prior probability that there *are bugs*.

Similarly, our posterior is also a probability, with $P(A | X)$ the probability there is no bug *given we saw all tests pass*, hence $1-P(A|X)$ is the probability there is a bug *given all tests passed*. What does our posterior probability look like? Below is a chart of both the prior and the posterior probabilities.

```python
figsize(12.5, 4)
colours = ["#348ABD", "#A60628"]

prior = [0.20, 0.80]
posterior = [1./3, 2./3]
plt.bar([0, .7], prior, alpha=0.70, width=0.25,
        color=colours[0], label="prior distribution",
        lw="3", edgecolor=colours[0])

plt.bar([0+0.25, .7+0.25], posterior, alpha=0.7,
        width=0.25, color=colours[1],
        label="posterior distribution",
        lw="3", edgecolor=colours[1])

plt.xticks([0.20, .95], ["Bugs Absent", "Bugs Present"])
plt.title("Prior and Posterior probability of bugs present")
plt.ylabel("Probability")
plt.legend(loc="upper left");
```

Notice that after we observed $X$ occur, the probability of bugs being absent increased. By increasing the number of tests, we can approach confidence (probability 1) that there are no bugs present.

This was a very simple example of Bayesian inference and Bayes rule. Unfortunately, the mathematics necessary to perform more complicated Bayesian inference only becomes more difficult, except for artificially constructed cases. We will later see that this type of mathematical analysis is actually unnecessary. First we must broaden our modeling tools. The next section deals with *probability distributions*. If you are already familiar, feel free to skip (or at least skim), but for the less familiar the next section is essential.

_______

## Probability Distributions

**Let's quickly recall what a probability distribution is:** Let $Z$ be some random variable. Then associated with $Z$ is a *probability distribution function* that assigns probabilities to the different outcomes $Z$ can take. Graphically, a probability distribution is a curve where the probability of an outcome is proportional to the height of the curve. You can see examples in the first figure of this chapter.

We can divide random variables into three classifications:

-   **$Z$ is discrete**: Discrete random variables may only assume values on a specified list. Things like populations, movie ratings, and number of votes are all discrete random variables. Discrete random variables become more clear when we contrast them with...

-   **$Z$ is continuous**: Continuous random variable can take on arbitrarily exact values. For example, temperature, speed, time, color are all modeled as continuous variables because you can progressively make the values more and more precise.

- **$Z$ is mixed**: Mixed random variables assign probabilities to both discrete and continuous random variables, i.e. it is a combination of the above two categories.

### Discrete Case
If $Z$ is discrete, then its distribution is called a *probability mass function*, which measures the probability $Z$ takes on the value $k$, denoted $P(Z=k)$. Note that the probability mass function completely describes the random variable $Z$, that is, if we know the mass function, we know how $Z$ should behave. There are popular probability mass functions that consistently appear: we will introduce them as needed, but let's introduce the first very useful probability mass function. We say $Z$ is *Poisson*-distributed if:

$$P(Z = k) =\frac{ \lambda^k e^{-\lambda} }{k!}, \; \; k=0,1,2, \dots $$

$\lambda$ is called a parameter of the distribution, and it controls the distribution's shape. For the Poisson distribution, $\lambda$ can be any positive number. By increasing $\lambda$, we add more probability to larger values, and conversely by decreasing $\lambda$ we add more probability to smaller values. One can describe $\lambda$ as the *intensity* of the Poisson distribution.

Unlike $\lambda$, which can be any positive number, the value $k$ in the above formula must be a non-negative integer, i.e., $k$ must take on values 0,1,2, and so on. This is very important, because if you wanted to model a population you could not make sense of populations with 4.25 or 5.612 members.

If a random variable $Z$ has a Poisson mass distribution, we denote this by writing

$$Z \sim \text{Poi}(\lambda) $$

One useful property of the Poisson distribution is that its expected value is equal to its parameter, i.e.:

$$E\large[ \;Z\; | \; \lambda \;\large] = \lambda $$

We will use this property often, so it's useful to remember. Below, we plot the probability mass distribution for different $\lambda$ values. The first thing to notice is that by increasing $\lambda$, we add more probability of larger values occurring. Second, notice that although the graph ends at 15, the distributions do not. They assign positive probability to every non-negative integer.

```python
figsize(12.5, 4)

import scipy.stats as stats
a = np.arange(16)
poi = stats.poisson
lambda_ = [1.5, 4.25]
colours = ["#348ABD", "#A60628"]

plt.bar(a, poi.pmf(a, lambda_[0]), color=colours[0],
        label="$\lambda = %.1f$" % lambda_[0], alpha=0.60,
        edgecolor=colours[0], lw="3")

plt.bar(a, poi.pmf(a, lambda_[1]), color=colours[1],
        label="$\lambda = %.1f$" % lambda_[1], alpha=0.60,
        edgecolor=colours[1], lw="3")

plt.xticks(a + 0.4, a)
plt.legend()
plt.ylabel("probability of $k$")
plt.xlabel("$k$")
plt.title("Probability mass function of a Poisson random variable; differing \
$\lambda$ values");
```

### Continuous Case
Instead of a probability mass function, a continuous random variable has a *probability density function*. This might seem like unnecessary nomenclature, but the density function and the mass function are very different creatures. An example of continuous random variable is a random variable with *exponential density*. The density function for an exponential random variable looks like this:

$$f_Z(z | \lambda) = \lambda e^{-\lambda z }, \;\; z\ge 0$$

Like a Poisson random variable, an exponential random variable can take on only non-negative values. But unlike a Poisson variable, the exponential can take on *any* non-negative values, including non-integral values such as 4.25 or 5.612401. This property makes it a poor choice for count data, which must be an integer, but a great choice for time data, temperature data (measured in Kelvins, of course), or any other precise *and positive* variable. The graph below shows two probability density functions with different $\lambda$ values.

When a random variable $Z$ has an exponential distribution with parameter $\lambda$, we say *$Z$ is exponential* and write

$$Z \sim \text{Exp}(\lambda)$$

Given a specific $\lambda$, the expected value of an exponential random variable is equal to the inverse of $\lambda$, that is:

$$E[\; Z \;|\; \lambda \;] = \frac{1}{\lambda}$$

```python
a = np.linspace(0, 4, 100)
expo = stats.expon
lambda_ = [0.5, 1]

for l, c in zip(lambda_, colours):
    plt.plot(a, expo.pdf(a, scale=1./l), lw=3,
             color=c, label="$\lambda = %.1f$" % l)
    plt.fill_between(a, expo.pdf(a, scale=1./l), color=c, alpha=.33)

plt.legend()
plt.ylabel("PDF at $z$")
plt.xlabel("$z$")
plt.ylim(0,1.2)
plt.title("Probability density function of an Exponential random variable;\
 differing $\lambda$");
```

### But what is $\lambda \;$?

**This question is what motivates statistics**. In the real world, $\lambda$ is hidden from us. We see only $Z$, and must go backwards to try and determine $\lambda$. The problem is difficult because there is no one-to-one mapping from $Z$ to $\lambda$. Many different methods have been created to solve the problem of estimating $\lambda$, but since $\lambda$ is never actually observed, no one can say for certain which method is best!

Bayesian inference is concerned with *beliefs* about what $\lambda$ might be. Rather than try to guess $\lambda$ exactly, we can only talk about what $\lambda$ is likely to be by assigning a probability distribution to $\lambda$.

This might seem odd at first. After all, $\lambda$ is fixed; it is not (necessarily) random! How can we assign probabilities to values of a non-random variable? Ah, we have fallen for our old, frequentist way of thinking. Recall that under Bayesian philosophy, we *can* assign probabilities if we interpret them as beliefs. And it is entirely acceptable to have *beliefs* about the parameter $\lambda$.

##### Example: Inferring behaviour from text-message data

Let's try to model a more interesting example, one that concerns the rate at which a user sends and receives text messages:

>  You are given a series of daily text-message counts from a user of your system. The data, plotted over time, appears in the chart below. You are curious to know if the user's text-messaging habits have changed over time, either gradually or suddenly. How can you model this? (This is in fact my own text-message data. Judge my popularity as you wish.)

```python
figsize(12.5, 3.5)
count_data = np.loadtxt("data/txtdata.csv")
n_count_data = len(count_data)
plt.bar(np.arange(n_count_data), count_data, color="#348ABD")
plt.xlabel("Time (days)")
plt.ylabel("count of text-msgs received")
plt.title("Did the user's texting habits change over time?")
plt.xlim(0, n_count_data);
```

Before we start modeling, see what you can figure out just by looking at the chart above. Would you say there was a change in behaviour during this time period?

How can we start to model this? Well, as we have conveniently already seen, a Poisson random variable is a very appropriate model for this type of *count* data. Denoting day $i$'s text-message count by $C_i$,

$$ C_i \sim \text{Poisson}(\lambda)  $$

We are not sure what the value of the $\lambda$ parameter really is, however. Looking at the chart above, it appears that the rate might become higher late in the observation period, which is equivalent to saying that $\lambda$ increases at some point during the observations. (Recall that a higher value of $\lambda$ assigns more probability to larger outcomes. That is, there is a higher probability of many text messages having been sent on a given day.)

How can we represent this observation mathematically? Let's assume that on some day during the observation period (call it $\tau$), the parameter $\lambda$ suddenly jumps to a higher value. So we really have two $\lambda$ parameters: one for the period before $\tau$, and one for the rest of the observation period. In the literature, a sudden transition like this would be called a *switchpoint*:

$$
\lambda =
\begin{cases}
\lambda_1  & \text{if } t \lt \tau \cr
\lambda_2 & \text{if } t \ge \tau
\end{cases}
$$

If, in reality, no sudden change occurred and indeed $\lambda_1 = \lambda_2$, then the $\lambda$s posterior distributions should look about equal.

We are interested in inferring the unknown $\lambda$s. To use Bayesian inference, we need to assign prior probabilities to the different possible values of $\lambda$. What would be good prior probability distributions for $\lambda_1$ and $\lambda_2$? Recall that $\lambda$ can be any positive number. As we saw earlier, the *exponential* distribution provides a continuous density function for positive numbers, so it might be a good choice for modeling $\lambda_i$. But recall that the exponential distribution takes a parameter of its own, so we'll need to include that parameter in our model. Let's call that parameter $\alpha$.

\begin{align}
&\lambda_1 \sim \text{Exp}( \alpha ) \\\
&\lambda_2 \sim \text{Exp}( \alpha )
\end{align}

$\alpha$ is called a *hyper-parameter* or *parent variable*. In literal terms, it is a parameter that influences other parameters. Our initial guess at $\alpha$ does not influence the model too strongly, so we have some flexibility in our choice.  A good rule of thumb is to set the exponential parameter equal to the inverse of the average of the count data. Since we're modeling $\lambda$ using an exponential distribution, we can use the expected value identity shown earlier to get:

$$\frac{1}{N}\sum_{i=0}^N \;C_i \approx E[\; \lambda \; |\; \alpha ] = \frac{1}{\alpha}$$

An alternative, and something I encourage the reader to try, would be to have two priors: one for each $\lambda_i$. Creating two exponential distributions with different $\alpha$ values reflects our prior belief that the rate changed at some point during the observations.

What about $\tau$? Because of the noisiness of the data, it's difficult to pick out a priori when $\tau$ might have occurred. Instead, we can assign a *uniform prior belief* to every possible day. This is equivalent to saying

\begin{align}
& \tau \sim \text{DiscreteUniform(1,70) }\\\\
& \Rightarrow P( \tau = k ) = \frac{1}{70}
\end{align}

So after all this, what does our overall prior distribution for the unknown variables look like? Frankly, *it doesn't matter*. What we should understand is that it's an ugly, complicated mess involving symbols only a mathematician could love. And things will only get uglier the more complicated our models become. Regardless, all we really care about is the posterior distribution.

We next turn to PyMC3, a Python library for performing Bayesian analysis that is undaunted by the mathematical monster we have created.

Introducing our first hammer: PyMC3
-----

PyMC3 is a Python library for programming Bayesian analysis [3]. It is a fast, well-maintained library. The only unfortunate part is that its documentation is lacking in certain areas, especially those that bridge the gap between beginner and hacker. One of this book's main goals is to solve that problem, and also to demonstrate why PyMC3 is so cool.

We will model the problem above using PyMC3. This type of programming is called *probabilistic programming*, an unfortunate misnomer that invokes ideas of randomly-generated code and has likely confused and frightened users away from this field. The code is not random; it is probabilistic in the sense that we create probability models using programming variables as the model's components. Model components are first-class primitives within the PyMC3 framework.

B. Cronin [5] has a very motivating description of probabilistic programming:

>   Another way of thinking about this: unlike a traditional program, which only runs in the forward directions, a probabilistic program is run in both the forward and backward direction. It runs forward to compute the consequences of the assumptions it contains about the world (i.e., the model space it represents), but it also runs backward from the data to constrain the possible explanations. In practice, many probabilistic programming systems will cleverly interleave these forward and backward operations to efficiently home in on the best explanations.

Because of the confusion engendered by the term *probabilistic programming*, I'll refrain from using it. Instead, I'll simply say *programming*, since that's what it really is.

PyMC3 code is easy to read. The only novel thing should be the syntax. Simply remember that we are representing the model's components ($\tau, \lambda_1, \lambda_2$ ) as variables.

```python
import pymc3 as pm
import theano.tensor as tt

with pm.Model() as model:
    alpha = 1.0/count_data.mean()  # Recall count_data is the
                                   # variable that holds our txt counts
    lambda_1 = pm.Exponential("lambda_1", alpha)
    lambda_2 = pm.Exponential("lambda_2", alpha)

    tau = pm.DiscreteUniform("tau", lower=0, upper=n_count_data - 1)
```

In the code above, we create the PyMC3 variables corresponding to $\lambda_1$ and $\lambda_2$. We assign them to PyMC3's *stochastic variables*, so-called because they are treated by the back end as random number generators.

```python
with model:
    idx = np.arange(n_count_data) # Index
    lambda_ = pm.math.switch(tau > idx, lambda_1, lambda_2)
```

This code creates a new function `lambda_`, but really we can think of it as a random variable: the random variable $\lambda$ from above. The `switch()` function assigns `lambda_1` or `lambda_2` as the value of `lambda_`, depending on what side of `tau` we are on. The values of `lambda_` up until `tau` are `lambda_1` and the values afterwards are `lambda_2`.

Note that because `lambda_1`, `lambda_2` and `tau` are random, `lambda_` will be random. We are **not** fixing any variables yet.

```python
with model:
    observation = pm.Poisson("obs", lambda_, observed=count_data)
```

The variable `observation` combines our data, `count_data`, with our proposed data-generation scheme, given by the variable `lambda_`, through the `observed` keyword.

The code below will be explained in Chapter 3, but I show it here so you can see where our results come from. One can think of it as a *learning* step. The machinery being employed is called *Markov Chain Monte Carlo* (MCMC), which I also delay explaining until Chapter 3. This technique returns thousands of random variables from the posterior distributions of $\lambda_1, \lambda_2$ and $\tau$. We can plot a histogram of the random variables to see what the posterior distributions look like. Below, we collect the samples (called *traces* in the MCMC literature) into histograms.

```python
### Mysterious code to be explained in Chapter 3.
with model:
    step = pm.Metropolis()
    trace = pm.sample(10000, tune=5000,step=step)
```

    Multiprocess sampling (4 chains in 4 jobs)
    CompoundStep
    >Metropolis: [tau]
    >Metropolis: [lambda_2]
    >Metropolis: [lambda_1]
    Sampling 4 chains: 100%|██████████| 60000/60000 [00:28<00:00, 2118.35draws/s]
    The number of effective samples is smaller than 25% for some parameters.

```python
lambda_1_samples = trace['lambda_1']
lambda_2_samples = trace['lambda_2']
tau_samples = trace['tau']
```

```python
figsize(12.5, 10)
#histogram of the samples:

ax = plt.subplot(311)
ax.set_autoscaley_on(False)

plt.hist(lambda_1_samples, histtype='stepfilled', bins=30, alpha=0.85,
         label="posterior of $\lambda_1$", color="#A60628", normed=True)
plt.legend(loc="upper left")
plt.title(r"""Posterior distributions of the variables
    $\lambda_1,\;\lambda_2,\;\tau$""")
plt.xlim([15, 30])
plt.xlabel("$\lambda_1$ value")

ax = plt.subplot(312)
ax.set_autoscaley_on(False)
plt.hist(lambda_2_samples, histtype='stepfilled', bins=30, alpha=0.85,
         label="posterior of $\lambda_2$", color="#7A68A6", normed=True)
plt.legend(loc="upper left")
plt.xlim([15, 30])
plt.xlabel("$\lambda_2$ value")

plt.subplot(313)
w = 1.0 / tau_samples.shape[0] * np.ones_like(tau_samples)
plt.hist(tau_samples, bins=n_count_data, alpha=1,
         label=r"posterior of $\tau$",
         color="#467821", weights=w, rwidth=2.)
plt.xticks(np.arange(n_count_data))

plt.legend(loc="upper left")
plt.ylim([0, .75])
plt.xlim([35, len(count_data)-20])
plt.xlabel(r"$\tau$ (in days)")
plt.ylabel("probability");
```

### Interpretation

Recall that Bayesian methodology returns a *distribution*. Hence we now have distributions to describe the unknown $\lambda$s and $\tau$. What have we gained? Immediately, we can see the uncertainty in our estimates: the wider the distribution, the less certain our posterior belief should be. We can also see what the plausible values for the parameters are: $\lambda_1$ is around 18 and $\lambda_2$ is around 23. The posterior distributions of the two $\lambda$s are clearly distinct, indicating that it is indeed likely that there was a change in the user's text-message behaviour.

What other observations can you make? If you look at the original data again, do these results seem reasonable?

Notice also that the posterior distributions for the $\lambda$s do not look like exponential distributions, even though our priors for these variables were exponential. In fact, the posterior distributions are not really of any form that we recognize from the original model. But that's OK! This is one of the benefits of taking a computational point of view. If we had instead done this analysis using mathematical approaches, we would have been stuck with an analytically intractable (and messy) distribution. Our use of a computational approach makes us indifferent to mathematical tractability.

Our analysis also returned a distribution for $\tau$. Its posterior distribution looks a little different from the other two because it is a discrete random variable, so it doesn't assign probabilities to intervals. We can see that near day 45, there was a 50% chance that the user's behaviour changed. Had no change occurred, or had the change been gradual over time, the posterior distribution of $\tau$ would have been more spread out, reflecting that many days were plausible candidates for $\tau$. By contrast, in the actual results we see that only three or four days make any sense as potential transition points.

### Why would I want samples from the posterior, anyways?

We will deal with this question for the remainder of the book, and it is an understatement to say that it will lead us to some amazing results. For now, let's end this chapter with one more example.

We'll use the posterior samples to answer the following question: what is the expected number of texts at day $t, \; 0 \le t \le 70$ ? Recall that the expected value of a Poisson variable is equal to its parameter $\lambda$. Therefore, the question is equivalent to *what is the expected value of $\lambda$ at time $t$*?

In the code below, let $i$ index samples from the posterior distributions. Given a day $t$, we average over all possible $\lambda_i$ for that day $t$, using $\lambda_i = \lambda_{1,i}$ if $t \lt \tau_i$ (that is, if the behaviour change has not yet occurred), else we use $\lambda_i = \lambda_{2,i}$.

```python
figsize(12.5, 5)
# tau_samples, lambda_1_samples, lambda_2_samples contain
# N samples from the corresponding posterior distribution
N = tau_samples.shape[0]
expected_texts_per_day = np.zeros(n_count_data)
for day in range(0, n_count_data):
    # ix is a bool index of all tau samples corresponding to
    # the switchpoint occurring prior to value of 'day'
    ix = day < tau_samples
    # Each posterior sample corresponds to a value for tau.
    # for each day, that value of tau indicates whether we're "before"
    # (in the lambda1 "regime") or
    #  "after" (in the lambda2 "regime") the switchpoint.
    # by taking the posterior sample of lambda1/2 accordingly, we can average
    # over all samples to get an expected value for lambda on that day.
    # As explained, the "message count" random variable is Poisson distributed,
    # and therefore lambda (the poisson parameter) is the expected value of
    # "message count".
    expected_texts_per_day[day] = (lambda_1_samples[ix].sum()
                                   + lambda_2_samples[~ix].sum()) / N

plt.plot(range(n_count_data), expected_texts_per_day, lw=4, color="#E24A33",
         label="expected number of text-messages received")
plt.xlim(0, n_count_data)
plt.xlabel("Day")
plt.ylabel("Expected # text-messages")
plt.title("Expected number of text-messages received")
plt.ylim(0, 60)
plt.bar(np.arange(len(count_data)), count_data, color="#348ABD", alpha=0.65,
        label="observed texts per day")

plt.legend(loc="upper left");
```

Our analysis shows strong support for believing the user's behavior did change ($\lambda_1$ would have been close in value to $\lambda_2$ had this not been true), and that the change was sudden rather than gradual (as demonstrated by $\tau$'s strongly peaked posterior distribution). We can speculate what might have caused this: a cheaper text-message rate, a recent weather-to-text subscription, or perhaps a new relationship. (In fact, the 45th day corresponds to Christmas, and I moved away to Toronto the next month, leaving a girlfriend behind.)

##### Exercises

1\.  Using `lambda_1_samples` and `lambda_2_samples`, what is the mean of the posterior distributions of $\lambda_1$ and $\lambda_2$?

```python
#type your code here.
```

2\.  What is the expected percentage increase in text-message rates? `hint:` compute the mean of `lambda_1_samples/lambda_2_samples`. Note that this quantity is very different from `lambda_1_samples.mean()/lambda_2_samples.mean()`.

```python
#type your code here.
```

3\. What is the mean of $\lambda_1$ **given** that we know $\tau$ is less than 45.  That is, suppose we have been given new information that the change in behaviour occurred prior to day 45. What is the expected value of $\lambda_1$ now? (You do not need to redo the PyMC3 part. Just consider all instances where `tau_samples < 45`.)

```python
#type your code here.
```

### References

-  [1] Gelman, Andrew. N.p.. Web. 22 Jan 2013. [N is never large enough](http://andrewgelman.com/2005/07/31/n_is_never_larg).
-  [2] Norvig, Peter. 2009. [The Unreasonable Effectiveness of Data](http://static.googleusercontent.com/media/research.google.com/en//pubs/archive/35179.pdf).
- [3] Salvatier, J, Wiecki TV, and Fonnesbeck C. (2016) Probabilistic programming in Python using PyMC3. *PeerJ Computer Science* 2:e55 <https://doi.org/10.7717/peerj-cs.55>
- [4] Jimmy Lin and Alek Kolcz. Large-Scale Machine Learning at Twitter. Proceedings of the 2012 ACM SIGMOD International Conference on Management of Data (SIGMOD 2012), pages 793-804, May 2012, Scottsdale, Arizona.
- [5] Cronin, Beau. "Why Probabilistic Programming Matters." 24 Mar 2013. Google, Online Posting to Google . Web. 24 Mar. 2013. <https://plus.google.com/u/0/107971134877020469960/posts/KpeRdJKR6Z1>.

```python
from IPython.core.display import HTML
def css_styling():
    styles = open("../styles/custom.css", "r").read()
    return HTML(styles)
css_styling()
```

<style>
    @font-face {
        font-family: "Computer Modern";
        src: url('http://9dbb143991406a7c655e-aa5fcb0a5a4ec34cff238a2d56ca4144.r56.cf5.rackcdn.com/cmunss.otf');
    }
    div.cell{
        width:800px;
        margin-left:16% !important;
        margin-right:auto;
    }
    h1 {
        font-family: Helvetica, serif;
    }
    h4{
        margin-top:12px;
        margin-bottom: 3px;
       }
    div.text_cell_render{
        font-family: Computer Modern, "Helvetica Neue", Arial, Helvetica, Geneva, sans-serif;
        line-height: 145%;
        font-size: 130%;
        width:800px;
        margin-left:auto;
        margin-right:auto;
    }
    .CodeMirror{
            font-family: "Source Code Pro", source-code-pro,Consolas, monospace;
    }
    .prompt{
        display: None;
    }
    .text_cell_render h5 {
        font-weight: 300;
        font-size: 22pt;
        color: #4057A1;
        font-style: italic;
        margin-bottom: .5em;
        margin-top: 0.5em;
        display: block;
    }

    .warning{
        color: rgb( 240, 20, 20 )
        }
</style>

```python

```

ext: ipynb


text:
# Simple nonlinear complementarity problem

**Randall Romero Aguilar, PhD**

This demo is based on the original Matlab demo accompanying the  Computational Economics and Finance 2001 textbook by Mario Miranda and Paul Fackler.

Original (Matlab) CompEcon file: **demslv13.m**

Running this file requires the Python version of CompEcon. This can be installed with pip by running

    !pip install compecon --upgrade

Last updated: 2021-Oct-01

## About

The problem is to solve

\begin{equation}
f(x, y) = \begin{bmatrix}1+xy -2x^3-x\\ 2x^2-y\end{bmatrix}
\end{equation}

subject to $0 \leq x, y \leq 1$

```python
from compecon import MCP, jacobian
import numpy as np

x0 = [0.5, 0.5]
```

### Solving the problem without the Jacobian
To solve this problem we create a **MCP** object using a lambda function.

```python
def func(z):
    x, y = z
    return np.array([1 + x*y - 2*x**3 - x, 2*x**2 - y])

F = MCP(func, [0, 0], [1,1])
```

Solve for initial guess $x_0 = [0.5, 0.5]$

```python
x = F.zero(x0, transform='minmax', show=True)
#x = F.zero(x0, transform='ssmooth', show=True) # FIXME: generates error
print(f'Solution is {x=}.')
```

    Using the MINMAX transformation
    Solving nonlinear equations by Broyden's method
    it    bstep  change
    --------------------
       0     1  1.56e-01
       1     1  1.14e-01
       2     1  1.01e-01
       3     0  3.66e-02
       4     0  1.33e-03
       5     0  1.86e-05
       6     0  2.02e-07
       7     0  1.13e-09
    Solution is x=array([0.7937, 1.    ]).

### Solving the problem with the Jacobian

```python
def func2(z):
    x, y = z
    f = [1 + x*y - 2*x**3 - x, 2*x**2 - y]
    J = [[y-6*x**2-1, x],[4*x, -1]]
    return np.array(f), np.array(J)

F2 = MCP(func2, [0, 0], [1,1])
```

```python
x = F2.zero(x0, transform='minmax', show=True)
print(f'Solution is {x=}.')
```

    Using the MINMAX transformation
    Solving nonlinear equations by Newton's method
    it    bstep  change
    --------------------
       0     1  1.56e-01
       1     0  9.84e-03
       2     0  3.19e-05
       3     0  3.40e-10
    Solution is x=array([0.7937, 1.    ]).

max_stars_repo_path: _build/jupyter_execute/notebooks/slv/13 Simple nonlinear complementarity problem.ipynb

max_stars_repo_name: randall-romero/CompEcon-python

max_issues_repo_name: randall-romero/CompEcon

max_forks_repo_name: randall-romero/CompEcon-python


text:
# メモ
やさしく学べる離散数学
石村園子

を読む。

# このノートの読み方
1. colab で開く
1. 上の方のメニューでドライブにコピーを押す
1. セルを実行しながら読む

目次
1. 集合と論理
1. 関係と写像 => 有向グラフ、関係行列
1. 代数系 => 2項演算、半群と群、環と体
1. 重所集合と束
1. グラフ => ツリー、オイラーグラフ、ハミルトングラフ、有限オートマトン

# 集合
set

要素 element または元と言う。

$a \in A$ または $A \owns a$

$a \notin A$ または $A \not \owns a$

メモ $\quad$ latex の `\not \owns` で $\not \owns$ になった。 $\not \in$ は `\not \in`  $\neq$ は `\neq` と $\not =$ `\not =` のどちらでもオッケー。

集合は普通 `{}` 波括弧で表す。 latex の中では波括弧は `\` でエスケープすれば使える。

$A=\{1,2,3,4,5\}$

$A=\{n \ |\  1 \leq n \leq5,\  n は自然数\}$

$n$ は自然数は、$n \in \mathbb N$ とか $n \in \boldsymbol N$ と書く。

数の集合の表記 $\boldsymbol N$ 自然数、$\boldsymbol Z$ 整数、$\boldsymbol Q$ 有理数、$\boldsymbol R$ 実数、$\boldsymbol C$ 複素数。

黒板ボールドという書体もある。 $\mathbb N$ 自然数、$\mathbb Z$ 整数、$\mathbb Q$ 有理数、$\mathbb R$ 実数、$\mathbb C$ 複素数。

**全体集合** $\quad U$

**空集合** $\quad \emptyset$ $\phi$ $\varnothing \quad$ `\emptyset$ $\phi$ $\varnothing`

集合 $A$ の要素が $B$ の要素であるとき、$A$ は $B$ の部分集合である、と言い、

$A \subseteq B$ または $B \supseteq A$

で表す。

$\quad \subseteq \quad$ `\subseteq`

$\quad \supseteq \quad$ `\supseteq`

$A \subseteq B \overset{\text{def}}{\Longleftrightarrow}$ $[x \in A \implies x \in B]$

$A \subseteq B \overset{\text{def}}{\iff}$ $[x \in A \implies x \in B]$

メモ $\quad$ $\subset$ は $\subseteq$ を含まない真部分集合とする、というのが教育上の建前のようだが、数学的に無理があると思う。 実際いくつかのテキストは、$A\subset B \land  A\neq B$ のとき真部分集合と言う定義をしている。
教育上は日本もアメリカも Subset はそれ自身を含まない。含む場合は $\subseteq$ で明確にする、というルールのようである。

集合はある全体集合 $U$ の部分集合として考え、特に $U \subseteq U$ である。

空集合 $\phi$ はすべての集合の部分集合と定めておく。

**定義** $\quad A \subseteq B \land B \subseteq A$ が成立するとき、$A$ と $B$ は等しいと言い、$A=B$ と書く。

**定義** $\quad A \subseteq B \land A \neq B$ のとき $A$ は $B$ の真部分集合であると言い、 $A \subset B$と書く。

メモ $\quad$ これは別の流儀もあるとのことなので注意して使う。

この本では記号 $A \subseteq B$ を次の意味で使う。

$\quad \quad A \subseteq B \Longleftrightarrow [A \subset B \lor A=B]$

**数の集合**

$\quad\quad \boldsymbol N \subset \boldsymbol Z \subset \boldsymbol Q \subset \boldsymbol R \subset \boldsymbol C$

定義 $\quad$ 冪集合 べき集合

集合 $A$ のすべての部分集合からなる集合を $A$ の冪集合と言い、 $\mathcal P(A)$ または $2^A$ などで表す。

例題 2 p.7

(1) 次の集合の冪集合を集める。

1. $A=\{a,b\}\quad$  #=> $\mathcal P (A) = \{\phi,\{a\},\{b\},\{a,b\}\}$
2. $B=\{\phi\}\quad$ #=> $\mathcal P (B) = \{\phi,\{\phi\}\}$

(2) $A$ の部分集合の中で $\{a\} \subseteq X$ となる集合 $X$ をすべて求める。

=>  $\mathcal P (A)$ の要素の中で $\{a\} \subseteq X$ となる集合は $\{a\},\{a,b\}$

# いまここ

練習問題 2 p.7

(1) 次の集合の冪集合を求める。

1. $C=\{a,b,c\}\quad$ #=> $\mathcal P (C) = \{\phi,\{a\},\{b\},\{c\},\{a,b\},\{a,c\},\{b,c\},\{a,b,c\}\}$

2. $D=\{\phi,\{\phi\}\}\quad$ #=> $\mathcal P (D) = \{\phi,\{\phi\},\{\{\phi\}\},\{\phi,\{\phi\}\}\}$

(2) $C$ の部分集合のなkで $\{a\}\subseteq X$ となる集合 $X$ をすべて求める。

=> $\{a\},\{a,b\},\{a,c\},\{a,b,c\}$

```python
# べき集合のことを power set と言うらしい
from itertools import *
def powerset(iterable):
    "powerset([1,2,3]) --> () (1,) (2,) (3,) (1,2) (1,3) (2,3) (1,2,3)"
    s = list(iterable)
    return chain.from_iterable(combinations(s, r) for r in range(len(s)+1))

list(powerset(['a','b','c']))
```

    [(),
     ('a',),
     ('b',),
     ('c',),
     ('a', 'b'),
     ('a', 'c'),
     ('b', 'c'),
     ('a', 'b', 'c')]

全称記号 $\forall$ と 存在記号 $\exists$ の使い方

$\forall x \in \mathbf Z, \ x \in \mathbf R$

日本語訳 $\quad \mathbf Z$ に属するすべての (任意の) 要素 $x$ は $\mathbf R$ に属する。 整数はすべて実数である。

$\exists x \in \mathbf Q, \ x \in \mathbf Z$

日本語訳 $\quad \mathbf Q$ には $\mathbf Z$ に属する要素 $x$ が存在する。有理数の中には整数であるものが存在する。

$\forall$ と $\exists$ の順番で意味が変わる。

$\forall x \in \mathbf Z, \ \exists y \in \mathbf Z, \ x+y=1$

日本語訳 $\quad \mathbf Z$ の任意の要素 $x$ に対し、$\mathbf Z$ のある要素 $y$ が存在して、$x+y=1$ が成立する。
$\mathbf Z$ の任意の要素 $x$ に対し、$x+y=1$ となるような要素 $y$ が $\mathbf Z$ に存在する。

$\exists x \in \mathbf R,\ \forall y \in \mathbf R, xy=y$

日本語訳 $\quad \mathbf R$
にある要素 $x$ が存在して、$\mathbf R$ の任意の要素 $y$ に対して $xy=y$ が成立する。
$\mathbf R$ の中には、$\mathbf R$ の任意の要素 $y$ に対して $xy=y$ が成立するような要素 $x$ が存在する。 (この $x$ は 1 のことである。)

例題 3 p.8 日本語訳、論理記号訳を作る。

(1) $\forall x \in \mathbf R, \ x^2 \geq 0 \quad$

日本語訳: すべての実数 $x$ の 2乗は 0 以上である。

(2) $\exists x \in \mathbf R, \ x^2+x \leq 0 \quad$

日本語訳: $x^2+x \leq 0$ を満たす $x$ が $\mathbf R$ に存在する。

(3) $\forall n \in \mathbf Z,\ \exists m \in \mathbf Z,\ n+m=0 \quad$

日本語訳: $\mathbf Z$ の任意の $n$ について、$\mathbf Z$ に $n+m=0$ を満足させる $m$ が存在する。

(4) $\exists a \in \mathbf Z,\ \forall x \in \mathbf R,\ x^2-2x \geq a \quad $

日本語訳: すべての実数 $x$ について $x^2-2x \geq a$ が成立する整数 $a$ が存在する。

(5) 複素数の中には、2 乗すると実数になる数が存在する。

論理記号訳: $\exists z \in \mathbf C,\ z^2 \in \mathbf R$

(6) すべての実数 $x$ について、$x^2-2x+2 \geq 0$ が成立する。

論理記号訳: $\forall x \in \mathbf R,\ x^2-2x+2 \geq 0$

(7) 0 と異なる任意の実数 $x$ に対して、$xy=1$ となる実数 $y$ が存在する。

論理記号訳: $\forall x \in \mathbf R^*,\ \exists y \in \mathbf R,\ xy = 1 \quad$ ただし、$\mathbf R^* = \{x \ |\ x \in \mathbf R,\ x \neq 0\}$

(8) 任意の整数 $n$ に対し、 $n+a=n$ となる定数 $a$ が整数の中に存在する。

論理記号訳: $\forall n \in \mathbf Z,\ \exists a \in \mathbf Z,\ n+a=n$ ($a$ は 0 のこと)

メモ

$\mathrm R ^* = \mathrm R\backslash \{0\}$

$\mathrm R ^* = \{x \ |\ x \in \mathbf R,\ x \neq 0\}$

練習問題 3 p.9

(1) $\exists n \in \mathbf N,\ \sqrt n \in \mathbf N$

日本語訳: その平方根をとると自然数であるような自然数が存在する。

(2) $\forall x \in \mathbf R^*, \frac1x \in \mathbf R$

日本語訳: 0 を除く任意の実数を分母とし、分子を 1 とする分数は実数である。

(3) $\forall a \in \mathbf R^*,\ \forall x \in \mathbf R,\ \exists x \in \mathbf R,\ ax+b=0$

日本語訳: 実数 $a,b$ ただし $a \neg 0$ について、$ax+b=0$ となる実数 $x$ が存在する。

(4) $\exists a \in \mathbf R,\ \forall x \in \mathbf R,\ |\sin x| \leq a$

日本語訳: 任意の実数 $x$ について $\sin x$ の絶対値が $a$ 以下であるような実数 $a$ が存在する。

(5) 複素数の中には絶対値が 1 となる数が存在する。

論理記号訳: $\exists z \in \mathbf C,\ |z| = 1$

(6)  任意の実数 $a$ に対し、$x^2+x+2 \gt a$ となる有理数 $x$ が存在する。

論理記号訳: $\forall a \in \mathbf R ,\ \exists x \in \mathbf Q ,\ x^2+x+2 \gt a$

集合の演算

定義

和集合

$A \cup B \overset{\text{def}}{\Longleftrightarrow} \{x \ |\ x \in \lor \ x \in  B\}$

積集合

$A \cap B \overset{\text{def}}{\Longleftrightarrow} \{x \ |\ x \in \land \ x \in  B\}$

補集合

$\overline A \overset{\text{def}}{\Longleftrightarrow} \{x \ |\ x \not \in A \}$

# ベン図 Venn diagram

```python
# ベン図 venn diagram
from matplotlib_venn import venn2
import matplotlib.pyplot as plt
from sympy import *

A = FiniteSet(1, 2)
B = FiniteSet(2, 3, 4)

venn2(subsets=(A,B))
plt.show()
print()
venn2(subsets=(15,10,5), set_labels=('X','Y'))
plt.show()
```

定理 1.1

冪等律

$A \cup A = A, \quad A \cap A = A$

交換律

$A \cup B = B\cup A,\quad A\cap B=B\cap A$

結合律

$(A \cup B)\cup C = A \cup (B \cup C)$

$(A \cap B)\cap C = A \cap (B \cap C)$

分配律

$A \cup (B \cap C) = (A \cup B) \cap (A \cup C)$

$A \cap (B \cup C) = (A \cap B) \cup (A \cap C)$

吸収律

$A \cup (A \cap B) = A$

$A \cap (A \cup B) = A$

定理 1.2 p.11

1. $A \cup \phi = A,\quad A\cap \phi = \phi$
1. $A \cup U = U,\quad A \cap U = A$
1. $A \cup \overline A = U,\quad A \cap \overline A = \phi$
1. $\overline{\overline A} = A$

定理 1.3 p. 11

ド・モルガンの法則

1. $\overline{A \cup B}= \overline A \cap \overline B$

1. $\overline{A \cap B}= \overline A \cup \overline B$

例題 4 p.11

次の性質をベン図で確認する。

1. $A \cup (B \cap C) = (A \cup B) \cap (A \cup C)$

1. $\overline{A \cup B}= \overline A \cap \overline B$

```python
# ベン図 venn diagram
from matplotlib_venn import venn3
import matplotlib.pyplot as plt
venn3(subsets = (4, 4, 2, 4, 2, 2, 1), set_labels = ('A', 'B', 'C'))
plt.show()
```

```python
# ベン図 venn diagram
from matplotlib_venn import venn2
import matplotlib.pyplot as plt
venn2(subsets = (4, 4, 2), set_labels = ('A', 'B'))
plt.show()
```

練習問題 4 p.11 次の性質をベン図で確かめる。
1. $A\cap(B \cup C) = (A\cap B) \cup (A\cap C)$

1. $\overline{(A \cap B)}=\overline A \cup \overline B$

例題 5 p.12

$U=\{n \ |\ 1 \leq n \leq 15,\ n \in \mathbf Z\}$

を全体集合として、部分集合

$A = \{a \ |\ a は素数\}$

$B = \{b \ |\ b は奇数\}$

$C = \{c \ |\ c は 3 の倍数\}$

を考える。

1. $A,B,C$ の要素を列挙して、ベン図に表す。

2. $B \cup C, B \cap C, \overline{A \cup C}, \overline B \cap C$ を求めてみる。

```python
# ベン図 venn diagram
from matplotlib_venn import venn3
import matplotlib.pyplot as plt

A = set([2,3,5,7,11,13])
B = set([1,3,5,7,9,11,13,15])
C = set([3,6,9,12,15])
U = set([1,2,3,4,5,6,7,8,9,10,11,12,13,14,15])

venn3(subsets=(A,B,C))
plt.show()
```

$A \cup B \overset{\text{def}}{\Longleftrightarrow} \{x \ |\ x \in A \lor x \in B\}$

$A \cap B \overset{\text{def}}{\Longleftrightarrow} \{x \ |\ x \in A \land x \in B\}$

$\overline A \overset{\text{def}}{\Longleftrightarrow} \{x \ |\ x \not \in A \}$

練習問題 5 p.12

$U = \{ n \ |\ n は90の正の約数\}$ を全体集合とし、$U$ の部分集合を

$A=\{a \ |\ a は 2 の倍数\}$、$B=\{b \ |\ b は 3 の倍数\}$、$C=\{c \ |\ c は 5 の倍数\}$

とするとき、次の集合の要素を求める。

1. $U,A,B,C$
1. $A\cup B, A\cap C,\overline B, B\cup \overline C, A \cup (\overline B \cap \overline C)$

```python
# 問題おかしくないか?
from matplotlib_venn import venn3
import matplotlib.pyplot as plt

A = set([2,6,10,18,30,90])
B = set([3,6,9,15,18,30,45,90])
C = set([5,10,15,30,45,90])
U = set([1,2,3,5,6,9,10,15,18,30,45,90])

print(A.union(B))
print(A.intersection(C))
print(U.difference(B))
print(C.union(U.difference(C)))
print(A.union((U.difference(B)).intersection(C)))
print(sorted(A.union((U.difference(B)).intersection(C))))
venn3(subsets=(A,B,C))
plt.show()
```

要素の個数 p.13

定義

要素の数が有限である集合を有限集合と言う。

有限集合でない集合を無限集合と言う。

説明 p.13

$A_1=\{2,4,6,\cdots,98,100\}$

$B_1=\{1,3,5,\cdots,99\}$

などは有限集合であるが、

$A_2=\{2,4,6,\cdots\} = \{m \ |\ m = 2k, \ k \in \mathbf N\}$

$B_2=\{1,3,5,\cdots\} = \{m \ |\ m = 2k-1, \ k \in \mathbf N\}$

や、$\mathbf{N,Z,Q,R,C}$ はすべて無限集合である。

これらの数の集合の中で $\mathbf Z$ と $\mathbf Q$ の要素はともに $\mathbf N$ の要素と 1 対 1 に対応をつけることができるが、$\mathbf R$ と $\mathbf C$ の要素は $\mathbf N$ の要素と 1 対 1 に対応をつけることができない。

有限集合および $\mathbf {N,Z,Q}$ は離散集合、可算集合、可付番集合などと呼ばれる。

$A$ が有限集合の場合、要素の数を

$n(A), |A|, \#(A)$

などで表す。

有限集合の要素の数については次の定理が成立する。

定理 1.4 p.13

$U$ を有限な全体集合とする。

$A,B$ を $U$ の部分集合とするとき、次の式が成立する。

1. $n(A\cup B) = n(A)+n(B)-n(A\cap B)$

1. $n(\overline A) =n(U) - n(A)$

例題 6 p.14

$U=\{m\ |\ 1 \leq m \leq 100, \ m \in \mathbf N\}$ を全体集合とし

$A=\{m\ |\ m=5k,\ k\in \mathbf N\}$、$B=\{m\ |\ m=7k,\ k\in \mathbf N\}$

とするとき、次の値を求める。

1. $n(A),\ n(B),\ n(A \cap B)$
2. $n(A\cup B),\ n(\overline A),\ n(B)$
3. $n(\overline A \cup \overline B),\ n(\overline A \cap \overline B)$

ド・モルガンの法則をつかってもよい。

# 命題

定義

真か偽か 真か偽か 真か偽かどちらか一方に明確に定まる主張を命題と言う。

```python
# python では 真偽は True , False
0 == 0
```

    True

命題を $p,q,r,\cdots$ などの英小文字で表し、その中身は

$p=$ 犬は動物である

$p=[5+3=8]$

などと書く。

命題の真偽を真理値と言い、$\mathrm{T, F}$ の記号を使う。

命題 $p$ が真のときは、$\quad p = \mathrm T$

などと書く。

定義

変数にある特定の要素を代入すると真偽が定まる主張を命題関数または述語 predicate と言う。

命題関数は $x$ を変数として $p(x)$ などで表す。

例題 7 p.16 命題かどうか。命題ならば真理値を求める。

1. 100 は自然数である。
1. 100 は大きい数である。
1. 人間は優れた生物である。
1. 地球は平らである。

1. $p(x)=x$ は動物である $\quad$ について、$x=$ 石、$x=$ 馬
1. $q(x)=[x^2-1=0]$ $\quad$ について $q(0),\ q(1)$

練習問題 7 p.16 命題かどうか。命題ならば真理値を求める。

1. $\sqrt 2$ は有理数である。
1. $\sqrt 2$ は実数である。
1. 数学は美しい学問である。
1. スイカは果物である。

1. $p(x)=x$ は植物である $\quad$ について $x=$ ウサギ、$x=$ ユリ
1. $q(x)=[x^2-x-2=0]$ $\quad$ について $q(0),\ q(-1)$

命題論理と述語論理
p. 17

わかりません。

命題論理はただの真偽で、述語論理は関係を表す。

命題論理は $P(x)$ で、述語論理は $Q(x,y)$ 。

述語論理は $\forall, \exists$ を含む。

# いまここ

論理演算

p.18

定義

命題 $p,q$ の真偽 $\mathrm{T,F}$ により右表の真理値をもつ命題 $p\lor q, p\land q$ について

$p$ $\ $ $q$ | $p$ $\lor$ $q$ | $p$ $\land$ $q$
--- | :-: | :-:
$\mathrm{T}$ $\ \ $ $\mathrm{T}$ | $\mathrm{T}$ | $\mathrm{T}$
$\mathrm{T}$ $\ \ $ $\mathrm{F}$ | $\mathrm{T}$ | $\mathrm{F}$
$\mathrm{F}$ $\ \ $ $\mathrm{T}$ | $\mathrm{T}$ | $\mathrm{F}$
$\mathrm{F}$ $\ \ $ $\mathrm{F}$ | $\mathrm{F}$ | $\mathrm{F}$

# 命題と真理値

真偽が明確に定まる式や文を命題 proposition と言う。

$x+1=2$ のように変数 $x$ が含まれている場合、命題にはならない。 $x$ に要素を代入したときに真理値が定まるものを, $x$ についての述語 predicate あるいは条件 condition と言い、$c(x),p(x)$ などと表す。

**合成命題**

compound proposition

George Bool 1815-1864 イギリスの数学者

**構成規則**

でない、かつ、または、ならば =>  論理結合子 logical connective =>
$\lnot , \land, \lor, \Rightarrow $

これらの文字や記号が構成規則に従って組合されてできた記号列を論理式 well-formed formula と言う。

1. 命題を表す文字は論理式である。
1. $p$ と $q$ が論旨式ならば、$(\lnot p), (p \land q), (p \lor q), (p \Rightarrow q), (p \Leftrightarrow q)$ は論理式である。

$r,s$ を命題とするとき、構成規則により、次の論理式が得られる。

$\quad r, (\lnot r), (r \land s), ((\lnot r) \land s), ((r \land s) \lor (\lnot r)), (s \Rightarrow (\lnot r)) $

# sympy の logic

https://docs.sympy.org/latest/modules/logic.html

論理式を表示してくれる部分しかわからないが、とりあえず colab 上で動くことを確認してみる。

```python
from sympy.abc import *
from sympy import *
display( y | (x & y))
display(x | y)
display(~x)
```

$\displaystyle y \vee \left(x \wedge y\right)$

$\displaystyle x \vee y$

$\displaystyle \neg x$

```python
from sympy.abc import *
from sympy import *
display(x  >> y)
display(Implies(x, y))
display(x << y)
display(Implies(y,x))
```

$\displaystyle x \Rightarrow y$

$\displaystyle x \Rightarrow y$

$\displaystyle y \Rightarrow x$

$\displaystyle y \Rightarrow x$

```python
(y & x).subs({x: True, y: True})
```

$\displaystyle \text{True}$

```python
(x | y).atoms()
```

    {x, y}

```python
# SOPform
minterms = [[0, 0, 0, 1], [0, 0, 1, 1],
            [0, 1, 1, 1], [1, 0, 1, 1], [1, 1, 1, 1]]
dontcares = [[0, 0, 0, 0], [0, 0, 1, 0], [0, 1, 0, 1]]
SOPform([w, x, y, z], minterms, dontcares)
```

$\displaystyle \left(y \wedge z\right) \vee \left(z \wedge \neg w\right)$

```python
# 整数で示すことができる
minterms = [1,3,7,11,15]
dontcares = {0,2,5}
SOPform([w, x, y, z], minterms, dontcares)
```

```python
# 辞書で示すことができる
minterms = [{w:0,x:1},{y:1,z:1,x:0}]
SOPform([w, x, y, z], minterms)
```

```python
# 混ぜて使うこともできる
minterms = [4,7,11,[1,1,1,1]]
dontcares =[{w:0,x:0,y:0},5]
SOPform([w,x,y,z],minterms,dontcares)
```

```python
# POSform
minterms = [[0, 0, 0, 1], [0, 0, 1, 1],
            [0, 1, 1, 1], [1, 0, 1, 1], [1, 1, 1, 1]]
dontcares = [[0, 0, 0, 0], [0, 0, 1, 0], [0, 1, 0, 1]]
POSform([w, x, y, z], minterms, dontcares)
```

$\displaystyle z \wedge \left(y \vee \neg w\right)$

```python
# 整数で示すことができる
minterms = [1,3,7,11,15]
dontcares = {0,2,5}
POSform([w, x, y, z], minterms, dontcares)
```

$\displaystyle z \wedge \left(y \vee \neg w\right)$

```python
# 辞書で示すことができる
minterms = [{w:0,x:1},{y:1,z:1,x:0}]
POSform([w, x, y, z], minterms)
```

$\displaystyle \left(x \vee y\right) \wedge \left(x \vee z\right) \wedge \left(\neg w \vee \neg x\right)$

```python
# 混ぜて使うこともできる
minterms = [4,7,11,[1,1,1,1]]
dontcares =[{w:0,x:0,y:0},5]
POSform([w,x,y,z],minterms,dontcares)
```

$\displaystyle \left(w \vee x\right) \wedge \left(y \vee \neg w\right) \wedge \left(z \vee \neg y\right)$

```python
from sympy import sympify, true, false, Or
display(sympify(True))
display((_ is True, _ is true))
display(Or(true,false))
display(_ is true)
display((true >> true, True >> True))
```

$\displaystyle \text{True}$

    (False, False)

$\displaystyle \text{True}$

    False

    (True, 0)

```python
display(true.as_set()) #=> U, UniversalSet
```

$\displaystyle \mathbb{U}$

```python
from sympy import sympify, true, false, Or
display(sympify(False))
display((_ is False, _ is false))
display(Or(true,false))
display(_ is true)
display((~false, ~False))
display((false >> false, False >> False))
```

$\displaystyle \text{False}$

    (False, False)

$\displaystyle \text{True}$

    False

    (True, -1)

    (True, 0)

```python
false.as_set() #=> EmptySet
```

$\displaystyle \emptyset$

```python
display(x & y)
display(And(x,y).subs(x,1))
display(x | y)
display(Or(x,y).subs(x,0))
display(Not(True))
display(Not(False))
display(And(True,False))
display(Or(True,False))
display(Not(And(And(True,x), Or(x,False))))
display(Not(And(Or(A,B), Or(Not(A), Not(B)))))

```

$\displaystyle x \wedge y$

$\displaystyle y$

$\displaystyle x \vee y$

$\displaystyle y$

$\displaystyle \text{False}$

$\displaystyle \text{True}$

$\displaystyle \text{False}$

$\displaystyle \text{True}$

$\displaystyle \neg x$

$\displaystyle \neg \left(\left(A \vee B\right) \wedge \left(\neg A \vee \neg B\right)\right)$

```python
display(~True)
display(~true)
```

    -2

$\displaystyle \text{False}$

```python
display(Xor(True,False))
display(Xor(True,True))
display(Xor(True,False,True,True,False)) #=> True の数が奇数だと Truelay(Xor(True,False))
display(Xor(True,True))
display(Xor(True,False,True,True,False)) #=> True の数が奇数だと True
display(Xor(True,False,True,False))
display(x ^ y)
```

$\displaystyle \text{True}$

$\displaystyle \text{False}$

$\displaystyle \text{True}$

$\displaystyle \text{False}$

$\displaystyle x \veebar y$

```python
Xor(x,y).subs(y,0)
```

$\displaystyle x$

```python
# Nand は頭から評価していって、False があれば即座に True、最後までFalseがなければ False
# その結果、False がひとつでもあれば True, なければ False になる
display(Nand(False,True))
display(Nand(True,True))
display(Nand(x,y))
```

$\displaystyle \text{True}$

$\displaystyle \text{False}$

$\displaystyle \neg \left(x \wedge y\right)$

```python
# Nor は頭から評価していって、True があれば即座に False、最後までTrueなければ True
# その結果、True がひとつでもあれば Fale, なければ True になる
display(Nor(False,True))
display(Nor(True,True))
display(Nor(False,False))
display(Nor(x,y))
```

$\displaystyle \text{False}$

$\displaystyle \text{False}$

$\displaystyle \text{True}$

$\displaystyle \neg \left(x \vee y\right)$

```python
# Implies は A が True で B が False のときのみ False を返し、それ以外では True である
from sympy.abc import *
from sympy import *
display(Implies(True, False))
display(Implies(False, False))
display(Implies(True, True))
display(Implies(False, True))
display(x >> y)
display(y << x)

```

$\displaystyle \text{False}$

$\displaystyle \text{True}$

$\displaystyle \text{True}$

$\displaystyle \text{True}$

$\displaystyle x \Rightarrow y$

$\displaystyle x \Rightarrow y$

```python
display(True >> False)
display(true >> false)

```

    1

$\displaystyle \text{False}$

```python
# Equivalent は A と B の両方が True または False のときに True
# Equivalent is True iff A and B are both True or both False
from sympy.abc import *
from sympy import *
Equivalent (False, False, False)
Equivalent (True, False, False)
Equivalent (x, And(x, True)) #=> なるほど!!!!
```

$\displaystyle \text{True}$

```python
# ITE if then else clause
# A が真なら B の結果、偽なら C の結果を返す
display(ITE(True, False, True))
display(ITE(Or(True,False),And(True,True),Xor(True,True)))
display(ITE(x,y,z))
display(ITE(True, x, y))
display(ITE(False, x, y))
display(ITE(x,y,y))

```

$\displaystyle \text{False}$

$\displaystyle \text{True}$

$\displaystyle ITE\left(x, y, z\right)$

$\displaystyle x$

$\displaystyle y$

$\displaystyle y$

```python
# ITE は if then else と言っても引数に bool 値以外はエラー
# display(ITE(True, [],())) #=> TypeError
```

```python
# to_cnf(expr, simplify=False, force=False)
# Convert a propositional logical sentence s to
# conjunctive normal form: ((A | ~B | …) & (B | C | …) & …).
from sympy.abc import *
from sympy import *
display(to_cnf(~(A|B)|D))
display(to_cnf((A|B)&(A|~A)))
display(to_cnf((A|B)&(A|~A),True))

```

$\displaystyle \left(D \vee \neg A\right) \wedge \left(D \vee \neg B\right)$

$\displaystyle \left(A \vee B\right) \wedge \left(A \vee \neg A\right)$

$\displaystyle A \vee B$

```python
# to_dnf(expr, simplify=False, force=False)
# Convert a propositional logical sentence s to
# disjunctive normal form: ((A & ~B & …) | (B & C & …) | …).
from sympy.abc import *
from sympy import *
from sympy.abc import C
display(to_dnf(B & (A | C)))
display(to_dnf((A&B)|(A&~B)|(B&C)|(~B&C)))
display(to_dnf((A&B)|(A&~B)|(B&C)|(~B&C),True))
```

$\displaystyle \left(A \wedge B\right) \vee \left(B \wedge C\right)$

$\displaystyle \left(A \wedge B\right) \vee \left(A \wedge \neg B\right) \vee \left(B \wedge C\right) \vee \left(C \wedge \neg B\right)$

$\displaystyle A \vee C$

```python
from sympy.abc import *
from sympy import *
from sympy.abc import C
from sympy.logic.boolalg import is_cnf
display(is_cnf(A|B|C))
display(is_cnf(A&B&C))
display(is_cnf((A&B)|C))
```

    True

    True

    False

```python
from sympy.abc import *
from sympy import *
from sympy.abc import C
from sympy.logic.boolalg import is_dnf
display(is_dnf(A|B|C))
display(is_dnf(A&B&C))
display(is_dnf((A&B)|C))
display(is_dnf(A&(B|C)))
```

    True

    True

    True

    False

```python
# simplify_logic(expr,form=None,deep=True,force=False)
from sympy.abc import *
from sympy import *
b=(~x&~y&~z) | (~x&~y& z)
display(b)
display(simplify_logic(b))
display(S(b)) #=> SingletonRegistry
display(simplify_logic(S(b)))
```

$\displaystyle \left(z \wedge \neg x \wedge \neg y\right) \vee \left(\neg x \wedge \neg y \wedge \neg z\right)$

$\displaystyle \neg x \wedge \neg y$

$\displaystyle \left(z \wedge \neg x \wedge \neg y\right) \vee \left(\neg x \wedge \neg y \wedge \neg z\right)$

$\displaystyle \neg x \wedge \neg y$

$\displaystyle B \wedge \left(A \vee C\right)$

```python
# bool_map(bool1, bool2)
from sympy.abc import *
from sympy import *
function1 = SOPform([x, z, y],[[1, 0, 1], [0, 0, 1]])
function2 = SOPform([a, b, c],[[1, 0, 1], [1, 0, 0]])
display(bool_map(function1, function2))
```

    (y & ~z, {y: a, z: b})

```python
eq =  Or(And(Not(y), w), And(Not(y), z), And(x, y))
eq2 = Or(And(Not(c), a), And(Not(c), d), And(b, c))
display(bool_map(eq, eq2))
eq = And(Xor(a, b), c, And(c,d))
display(bool_map(eq, eq.subs(c, x)))
```

    ((x & y) | (w & ~y) | (z & ~y), {w: a, x: b, y: c, z: d})

    (c & d & (a | b) & (~a | ~b), {a: a, b: b, c: d, d: x})

```python
# inference, satisfiable
from sympy.abc import *
from sympy import *
display(satisfiable(x & ~x))
display(satisfiable((x|y)&(x|~y)&(~x|y)))
display(satisfiable(x))
```

    False

    {x: True, y: True}

    {x: True}

```python
display(satisfiable(True))
```

    {True: True}

```python
display(next(satisfiable(A & ~A, all_models=True)))
models=satisfiable((A >> B)&B, all_models=True)
display(next(models))
def use_models(models):
    for model in models:
        if model:
            # do something with the model
            print(model)
        else:
            # Given expr is unsatisfiable
            print("UNSAT")
use_models(satisfiable(A>>~A,all_models=True))
use_models(satisfiable(A^A,all_models=True))
```

    False

    {A: False, B: True}

    {A: False}
    UNSAT

# 逆

```python

```

論理式 $p \Rightarrow q$ において、$p$ を仮定 hypothesis または前提 antecedent, premise 、 $q$ を結論 conclusion, consequent と言う。

仮定と結論を入れ替えた $q \Rightarrow p$ を $p \Rightarrow q$ の逆 converse と言う。

反例 counter example

# 必要十分条件

論理式 $p \Rightarrow q$ が成り立つとき

$q$ を $p$ (が成り立つため) の必要条件 necessary condition と言う。

$p$ を $q$ (が成り立つため) の十分条件 sufficient condition と言う。

$p \Longleftrightarrow q$ が成り立つとき

$p$ を $q$ (が成り立つため) の必要十分条件と言う。このとき $p$ と $q$ は同値であると言う。

恒真 トートロジー tautology

恒偽 contradiction

# 充足可能性問題 satisfiable problem SAT

ある論理式が真となるための各命題への真理値の割り当て方(組合せ)が少なくとも 1通りあるとき、すなわち恒真でない論理式を充足可能 sastifiable と言う。


text:
# Lab 6 - Gaussian Process

In general, Gaussian Process (GP) is a stochastic process used for modeling data, which were observed over time, space or both. Main thing that can characterise GP is that is a kind of generalization of normal probability distributions, where each of them describes a random variable (scalar or vector if we deal with multivariate distribution).

Following definition captures the essence of Gaussian Process:

A Gaussian process is a collection of random variables, any Gaussian process finite number of which have a joint Gaussian distribution.

If we define mean function $m(x)$ and the covariance function mean function $k(x,x')$ of a real process $f(x)$ as:
\begin{equation}
m(x) = E[f(x)]
\end{equation}
\begin{equation}
k(x,x') = E[f(x)-m(x))(f(x')-m(x'))]
\end{equation}
Then we can define GP can as:
\begin{equation}
\label{eq:gpmain}
f(x) \sim GP(m(x),k(x,x'))
\end{equation}

GP can be fully determined by only declaring mean and covariance functions. Mean in most cases is set to value "0", because such a setting can be useful, simplifies matters and is not a difficult requirement to fulfill. Of course, there are examples
where we would like to change the mean e.g. for better model interpretability or the
specification of our prior. Covariance function (also called kernel function) represents
a similarity between data points. Note that usually covariance is chosen from
the set of already defined functions. The one should at least pick one, which represents
prior beliefs of the problem. But in fact, covariance function can be any function, which
has the property of generating a positive definite covariance matrix. Anyway, creating
and defining new covariance functions, which will simultaneously be correct and have a
practical usage, can be really difficult

The most basic and common kernel function is Radial Basis Function (RBF), which is defined with the formula:
\begin{equation}
\label{eq:rbf}
k\left(x_{i}, x_{j}\right)=\exp \left(-\frac{d\left(x_{i}, x_{j}\right)^{2}}{2 l^{2}}\right)
\end{equation}
Its main property is that its value is usually only dependent on the distance from the specified point. The parameter which RBF kernel uses is $l$ as a characteristic length scale. RBF is infinitely differentiable. That means the GP with this kernel has a mean square derivatives for all orders. Other worth noting kernel, in terms of our classess is one from exponentiated quadratic family, where $\rho$ is length-scale and $\alpha$ is marginal deviation:
$$
k\left(x_{1}, x_{2}\right)=\alpha^{2} \exp \left(-\frac{1}{2}\left(\frac{\left|x_{1}-x_{2}\right|}{\rho}\right)^{2}\right)
$$
Gaussian Process is growing in popularity of usage in domains connected with statistics and machine learning. It is especially efficient solution, while we are dealing with a problem of obtaining data e.g. for faulty cases or just an overall lack of comprehensive data and if we want to obtain not only prediction but also confidence interval. Main areas where GP is used for are: regression, prediction, classification and identification. On this classes we will focuse on a simple task of regression.

### Task 1: simulating from Gaussian Process model

To make any simulations, first you need a set of observations:

```python
#Observed data
x_obs = [-10, -8, -6, -4, -2, 0, 2, 4, 6, 8, 10]
y_obs = [0.328572824089476, 4.20607004111644, 1.35507551134795,
0.161608755204364, -5.42320349780782, -3.05851276224202, -0.0764172642034502,
-4.55218472276499, -0.902226297922731, -5.8609833528976, -1.05854090910473]
```

Model also needs sampling space and the indexes of where our observed data resides there.
For your convenience, we provide them below:

```python
x_predict = np.linspace(-11,11,551,endpoint=True)
observed_idx = [26, 76, 126, 176, 226, 276, 326, 376, 426, 476, 526]
```

In this example we are using kernel from exponentiated quadratic family, so we need to provide its hyperparameters (alpha and rho). Moreover, for sample generation, we would like to set a sigma value. Let's start with some arbitrarily selected parameters:

```python
alpha = 3
rho = 5.5
sigma = 2
```

Now we can create a model and sample from it. Do not forget to define proper data value based on requirements from stan model.

```python
model = CmdStanModel(stan_file='gaussian_process.stan')
```

```python
result = model.sample(data=data, chains = 4)
```

Your task now is to:
- plot few (about 5-10) samples generated from model
- plot mean value from model with its confidence interval (standar deviation, using errorbar function)

Also on each plot place points with obsserved data. Repeat the steps for two other parameters sets of your choice

### Task 2: optimize hyperparameters

To achive better results of GP fitting for our observed data, we need to oprimize its parameters. We can do that by using maximum marginal likelihood estimation. It can be done by placing parameters in "parameter" block in stan model and using "optimize" method.

```python
model2 = CmdStanModel(stan_file='gaussian_process_optimize.stan')
```

```python
result2 = model2.optimize(data=data, seed=5838298)
```

Optimize parameters with mentioned method and generate new samples with model from task 1. Of course use just aquired parameters. Do not forget to specify data variable based on stan model.
What can your say about results in comparison to ones from task 1?
Did the GP perform well in the task of regression?

```python

```

max_stars_repo_path: Data Analytics 2022/Lab 6 - Gaussian Process/Lab6.ipynb


text:
# Penalised Regression

## YouTube Videos
1. **Scikit Learn Linear Regression:** https://www.youtube.com/watch?v=EvnpoUTXA0E
2. **Scikit Learn Linear Penalise Regression:** https://www.youtube.com/watch?v=RhsEAyDBkTQ

## Introduction
We often do not want the coefficients/ weights to be too large. Hence we append the loss function with a penalty function to discourage large values of $w$.

\begin{align}
\mathcal{L} & = \sum_{i=1}^N (y_i-f(x_i|w,b))^2 + \alpha \sum_{j=1}^D w_j^2 + \beta \sum_{j=1}^D |w_j|
\end{align}
where, $f(x_i|w,b) = wx_i+b$. The values of $\alpha$ and $\beta$ are positive (or zero), with higher values enforcing the weights to be closer to zero.

## Lesson Structure
1. The task of this lesson is to infer the weights given the data (observations, $y$ and inputs $x$).
2. We will be using the module `sklearn.linear_model`.

```python
import numpy as np
from sklearn.linear_model import LinearRegression

import matplotlib.pyplot as plt
%matplotlib inline

# In order to reproduce the exact same number we need to set the seed for random number generators:
np.random.seed(1)
```

A normally distributed random looks as follows:

```python
e = np.random.randn(10000,1)
plt.hist(e,100) #histogram with 100 bins
plt.ylabel('y')
plt.xlabel('x')
plt.title('Histogram of Normally Distributed Numbers')
plt.show()
```

Generate observations $y$ given feature (design) matrix $X$ according to:
$$
y = Xw + \xi\\
\xi_i \sim \mathcal{N}(0,\sigma^2)
$$

In this particular case, $w$ is a 100 dimensional vector where 90% of the numbers are zero. i.e. only 10 of the numbers are non-zero.

```python
# Generate the data
N = 40 # Number of observations
D = 100 # Dimensionality

x = np.random.randn(N,D) # get random observations of x
w_true = np.zeros((D,1)) # create a weight vector of zeros
idx = np.random.choice(100,10,replace=False) # randomly choose 10 of those weights
w_true[idx] = np.random.randn(10,1) # populate then with 10 random weights

e = np.random.randn(N,1) # have a noise vector
y = np.matmul(x,w_true) + e # generate observations

# create validation set:
N_test = 50
x_test = np.random.randn(50,D)
y_test_true = np.matmul(x_test,w_true)
```

```python
model = LinearRegression()
model.fit(x,y)

# plot the true vs estimated coeffiecients
plt.plot(np.arange(100),np.squeeze(model.coef_))
plt.plot(np.arange(100),w_true)
plt.legend(["Estimated","True"])
plt.title('Estimated Weights')
plt.show()
```

One way of testing how good your model is to look at metrics. In the case of regression Mean Squared Error (MSE) is a common metric which is defined as:
$$ \frac{1}{N}\sum_{i=1}^N \xi_i^2$$ where, $\xi_i = y_i-f(x_i|w,b)$. Furthermore it is best to look at the MSE on a validation set, rather than on the training dataset that we used to train the model.

```python
y_est = model.predict(x_test)
mse = np.mean(np.square(y_test_true-y_est))
print(mse)
```

    6.808967330386871

Ridge regression is where you penalise the weights by setting the $\alpha$ parameter right at the top. It penalises it so that the higher **the square of the weights** the higher the loss.

```python
from sklearn.linear_model import Ridge

model = Ridge(alpha=5.0,fit_intercept = False)
# TODO:
# Train the model: see how model was trained above, same code!
model.fit(x,y)
# plot the true vs estimated coeffiecients
plt.plot(np.arange(100),np.squeeze(model.coef_))
plt.plot(np.arange(100),w_true)
plt.legend(["Estimated","True"])
plt.show()
```

This model is slightly better than without any penalty on the weights.

```python
y_est = model.predict(x_test)
mse = np.mean(np.square(y_test_true-y_est))
print(mse)
```

    6.422880725012181

Lasso is a model that encourages weights to go to zero exactly, as opposed to Ridge regression which encourages small weights.

```python
from sklearn.linear_model import Lasso

# TODO:
# 1. Set the model to a Lasso below, set the alpha parameter and set `fit_intercept=False`.
# 2. Train the model (same as before)
model = Lasso(fit_intercept=False)

model.fit(x,y)

# plot the true vs estimated coeffiecients
plt.plot(np.arange(100),np.squeeze(model.coef_))
plt.plot(np.arange(100),w_true)
plt.legend(["Estimated","True"])
plt.title('Lasso regression weight inference')
plt.show()
```

The MSE is significantly better than both the above models.

```python
y_est = model.predict(x_test)[:,None]
mse = np.mean(np.square(y_test_true-y_est))
print(mse)
```

    4.795231099516383

Automated Relevance Determination (ARD) regression is similar to lasso in that it encourages zero weights. However, the advantage is that you do not need to set a penalisation parameter, $\alpha$, $\beta$ in this model.

```python
from sklearn.linear_model import ARDRegression

# TODO
# 1. Set model to ARDRegression and `fit_intercept=False`.
# This model is not sensitive to its parameters which makes it powerful
# 2. Train the model, same as above
model = ARDRegression(fit_intercept=False)
model.fit(x,y)
# plot the true vs estimated coeffiecients
plt.plot(np.arange(100),np.squeeze(model.coef_))
plt.plot(np.arange(100),w_true)
plt.legend(["Estimated","True"])
plt.show()
```

```python
y_est = model.predict(x_test)[:,None]
mse = np.mean(np.square(y_test_true-y_est))
print(mse)
```

    2.870211080521882

### Note:
Rerun the above with setting N=400

## Inverse Problems (Optional)
The following section is optional and you may skip it. It is not necessary for understanding Deep Learning.

Inverse problems are where given the outputs you are required to infer the inputs. A typical example is X-rays. Given the x-ray sensor readings, the algorithm needs to build an image of an individuals bone structure.

See [compressed_sensing](http://scikit-learn.org/stable/auto_examples/applications/plot_tomography_l1_reconstruction.html#sphx-glr-auto-examples-applications-plot-tomography-l1-reconstruction-py) for an example of l1 reguralisation applied to a compressed sensing problem (has a resemblance to the x-ray problem).

```python

```

max_stars_repo_path: jupyter/Keras_TensorFlow_Course/Lesson 01 - PenalisedRegression.ipynb

max_stars_repo_name: multivacplatform/multivac-dl

max_issues_repo_name: multivacplatform/multivac-dl

max_forks_repo_name: multivacplatform/multivac-dl