lars/FYS-STK4155
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity
Files
c464ea1051ca45a1ebb0eb151ee2e70ebe6489e1
FYS-STK4155/doc/LectureNotes/_build/html/statistics.html
T
KarlHenrik c464ea1051 week41 exercises small changes
2025-10-08 15:40:16 +02:00

1621 lines
83 KiB
HTML
Raw Blame History

This file contains ambiguous Unicode characters
This file contains Unicode characters that might be confused with other characters. If you think that this is intentional, you can safely ignore this warning. Use the Escape button to reveal them.
<!DOCTYPE html>
<html lang="en" data-content_root="./" >
<head>
<meta charset="utf-8" />
<meta name="viewport" content="width=device-width, initial-scale=1.0" /><meta name="viewport" content="width=device-width, initial-scale=1" />
<title>1. Elements of Probability Theory and Statistical Data Analysis &#8212; Applied Data Analysis and Machine Learning</title>
<script data-cfasync="false">
document.documentElement.dataset.mode = localStorage.getItem("mode") || "";
document.documentElement.dataset.theme = localStorage.getItem("theme") || "";
</script>
<!-- Loaded before other Sphinx assets -->
<link href="_static/styles/theme.css?digest=dfe6caa3a7d634c4db9b" rel="stylesheet" />
<link href="_static/styles/bootstrap.css?digest=dfe6caa3a7d634c4db9b" rel="stylesheet" />
<link href="_static/styles/pydata-sphinx-theme.css?digest=dfe6caa3a7d634c4db9b" rel="stylesheet" />
<link href="_static/vendor/fontawesome/6.5.2/css/all.min.css?digest=dfe6caa3a7d634c4db9b" rel="stylesheet" />
<link rel="preload" as="font" type="font/woff2" crossorigin href="_static/vendor/fontawesome/6.5.2/webfonts/fa-solid-900.woff2" />
<link rel="preload" as="font" type="font/woff2" crossorigin href="_static/vendor/fontawesome/6.5.2/webfonts/fa-brands-400.woff2" />
<link rel="preload" as="font" type="font/woff2" crossorigin href="_static/vendor/fontawesome/6.5.2/webfonts/fa-regular-400.woff2" />
<link rel="stylesheet" type="text/css" href="_static/pygments.css?v=03e43079" />
<link rel="stylesheet" type="text/css" href="_static/styles/sphinx-book-theme.css?v=eba8b062" />
<link rel="stylesheet" type="text/css" href="_static/togglebutton.css?v=13237357" />
<link rel="stylesheet" type="text/css" href="_static/copybutton.css?v=76b2166b" />
<link rel="stylesheet" type="text/css" href="_static/mystnb.8ecb98da25f57f5357bf6f572d296f466b2cfe2517ffebfabe82451661e28f02.css?v=6644e6bb" />
<link rel="stylesheet" type="text/css" href="_static/sphinx-thebe.css?v=4fa983c6" />
<link rel="stylesheet" type="text/css" href="_static/sphinx-design.min.css?v=95c83b7e" />
<!-- Pre-loaded scripts that we'll load fully later -->
<link rel="preload" as="script" href="_static/scripts/bootstrap.js?digest=dfe6caa3a7d634c4db9b" />
<link rel="preload" as="script" href="_static/scripts/pydata-sphinx-theme.js?digest=dfe6caa3a7d634c4db9b" />
<script src="_static/vendor/fontawesome/6.5.2/js/all.min.js?digest=dfe6caa3a7d634c4db9b"></script>
<script src="_static/documentation_options.js?v=9eb32ce0"></script>
<script src="_static/doctools.js?v=9a2dae69"></script>
<script src="_static/sphinx_highlight.js?v=dc90522c"></script>
<script src="_static/clipboard.min.js?v=a7894cd8"></script>
<script src="_static/copybutton.js?v=f281be69"></script>
<script src="_static/scripts/sphinx-book-theme.js?v=887ef09a"></script>
<script>let toggleHintShow = 'Click to show';</script>
<script>let toggleHintHide = 'Click to hide';</script>
<script>let toggleOpenOnPrint = 'true';</script>
<script src="_static/togglebutton.js?v=4a39c7ea"></script>
<script>var togglebuttonSelector = '.toggle, .admonition.dropdown';</script>
<script src="_static/design-tabs.js?v=f930bc37"></script>
<script>const THEBE_JS_URL = "https://unpkg.com/thebe@0.8.2/lib/index.js"; const thebe_selector = ".thebe,.cell"; const thebe_selector_input = "pre"; const thebe_selector_output = ".output, .cell_output"</script>
<script async="async" src="_static/sphinx-thebe.js?v=c100c467"></script>
<script>var togglebuttonSelector = '.toggle, .admonition.dropdown';</script>
<script>const THEBE_JS_URL = "https://unpkg.com/thebe@0.8.2/lib/index.js"; const thebe_selector = ".thebe,.cell"; const thebe_selector_input = "pre"; const thebe_selector_output = ".output, .cell_output"</script>
<script>window.MathJax = {"options": {"processHtmlClass": "tex2jax_process|mathjax_process|math|output_area"}}</script>
<script defer="defer" src="https://cdn.jsdelivr.net/npm/mathjax@3/es5/tex-mml-chtml.js"></script>
<script>DOCUMENTATION_OPTIONS.pagename = 'statistics';</script>
<link rel="index" title="Index" href="genindex.html" />
<link rel="search" title="Search" href="search.html" />
<link rel="next" title="2. Linear Algebra, Handling of Arrays and more Python Features" href="linalg.html" />
<link rel="prev" title="Textbooks" href="textbooks.html" />
<meta name="viewport" content="width=device-width, initial-scale=1"/>
<meta name="docsearch:language" content="en"/>
</head>
<body data-bs-spy="scroll" data-bs-target=".bd-toc-nav" data-offset="180" data-bs-root-margin="0px 0px -60%" data-default-mode="">
<div id="pst-skip-link" class="skip-link d-print-none"><a href="#main-content">Skip to main content</a></div>
<div id="pst-scroll-pixel-helper"></div>
<button type="button" class="btn rounded-pill" id="pst-back-to-top">
<i class="fa-solid fa-arrow-up"></i>Back to top</button>
<input type="checkbox"
class="sidebar-toggle"
id="pst-primary-sidebar-checkbox"/>
<label class="overlay overlay-primary" for="pst-primary-sidebar-checkbox"></label>
<input type="checkbox"
class="sidebar-toggle"
id="pst-secondary-sidebar-checkbox"/>
<label class="overlay overlay-secondary" for="pst-secondary-sidebar-checkbox"></label>
<div class="search-button__wrapper">
<div class="search-button__overlay"></div>
<div class="search-button__search-container">
<form class="bd-search d-flex align-items-center"
action="search.html"
method="get">
<i class="fa-solid fa-magnifying-glass"></i>
<input type="search"
class="form-control"
name="q"
id="search-input"
placeholder="Search this book..."
aria-label="Search this book..."
autocomplete="off"
autocorrect="off"
autocapitalize="off"
spellcheck="false"/>
<span class="search-button__kbd-shortcut"><kbd class="kbd-shortcut__modifier">Ctrl</kbd>+<kbd>K</kbd></span>
</form></div>
</div>
<div class="pst-async-banner-revealer d-none">
<aside id="bd-header-version-warning" class="d-none d-print-none" aria-label="Version warning"></aside>
</div>
<header class="bd-header navbar navbar-expand-lg bd-navbar d-print-none">
</header>
<div class="bd-container">
<div class="bd-container__inner bd-page-width">
<div class="bd-sidebar-primary bd-sidebar">
<div class="sidebar-header-items sidebar-primary__section">
</div>
<div class="sidebar-primary-items__start sidebar-primary__section">
<div class="sidebar-primary-item">
<a class="navbar-brand logo" href="intro.html">
<img src="_static/logo.png" class="logo__image only-light" alt="Applied Data Analysis and Machine Learning - Home"/>
<script>document.write(`<img src="_static/logo.png" class="logo__image only-dark" alt="Applied Data Analysis and Machine Learning - Home"/>`);</script>
</a></div>
<div class="sidebar-primary-item">
<script>
document.write(`
<button class="btn search-button-field search-button__button" title="Search" aria-label="Search" data-bs-placement="bottom" data-bs-toggle="tooltip">
<i class="fa-solid fa-magnifying-glass"></i>
<span class="search-button__default-text">Search</span>
<span class="search-button__kbd-shortcut"><kbd class="kbd-shortcut__modifier">Ctrl</kbd>+<kbd class="kbd-shortcut__modifier">K</kbd></span>
</button>
`);
</script></div>
<div class="sidebar-primary-item"><nav class="bd-links bd-docs-nav" aria-label="Main">
<div class="bd-toc-item navbar-nav active">
<ul class="nav bd-sidenav bd-sidenav__home-link">
<li class="toctree-l1">
<a class="reference internal" href="intro.html">
Applied Data Analysis and Machine Learning
</a>
</li>
</ul>
<p aria-level="2" class="caption" role="heading"><span class="caption-text">About the course</span></p>
<ul class="nav bd-sidenav">
<li class="toctree-l1"><a class="reference internal" href="schedule.html">Course setting</a></li>
<li class="toctree-l1"><a class="reference internal" href="teachers.html">Teachers and Grading</a></li>
<li class="toctree-l1"><a class="reference internal" href="textbooks.html">Textbooks</a></li>
</ul>
<p aria-level="2" class="caption" role="heading"><span class="caption-text">Review of Statistics with Resampling Techniques and Linear Algebra</span></p>
<ul class="current nav bd-sidenav">
<li class="toctree-l1 current active"><a class="current reference internal" href="#">1. Elements of Probability Theory and Statistical Data Analysis</a></li>
<li class="toctree-l1"><a class="reference internal" href="linalg.html">2. Linear Algebra, Handling of Arrays and more Python Features</a></li>
</ul>
<p aria-level="2" class="caption" role="heading"><span class="caption-text">From Regression to Support Vector Machines</span></p>
<ul class="nav bd-sidenav">
<li class="toctree-l1"><a class="reference internal" href="chapter1.html">3. Linear Regression</a></li>
<li class="toctree-l1"><a class="reference internal" href="chapter2.html">4. Ridge and Lasso Regression</a></li>
<li class="toctree-l1"><a class="reference internal" href="chapter3.html">5. Resampling Methods</a></li>
<li class="toctree-l1"><a class="reference internal" href="chapter4.html">6. Logistic Regression</a></li>
<li class="toctree-l1"><a class="reference internal" href="chapteroptimization.html">7. Optimization, the central part of any Machine Learning algortithm</a></li>
<li class="toctree-l1"><a class="reference internal" href="chapter5.html">8. Support Vector Machines, overarching aims</a></li>
</ul>
<p aria-level="2" class="caption" role="heading"><span class="caption-text">Decision Trees, Ensemble Methods and Boosting</span></p>
<ul class="nav bd-sidenav">
<li class="toctree-l1"><a class="reference internal" href="chapter6.html">9. Decision trees, overarching aims</a></li>
<li class="toctree-l1"><a class="reference internal" href="chapter7.html">10. Ensemble Methods: From a Single Tree to Many Trees and Extreme Boosting, Meet the Jungle of Methods</a></li>
</ul>
<p aria-level="2" class="caption" role="heading"><span class="caption-text">Dimensionality Reduction</span></p>
<ul class="nav bd-sidenav">
<li class="toctree-l1"><a class="reference internal" href="chapter8.html">11. Basic ideas of the Principal Component Analysis (PCA)</a></li>
<li class="toctree-l1"><a class="reference internal" href="clustering.html">12. Clustering and Unsupervised Learning</a></li>
</ul>
<p aria-level="2" class="caption" role="heading"><span class="caption-text">Deep Learning Methods</span></p>
<ul class="nav bd-sidenav">
<li class="toctree-l1"><a class="reference internal" href="chapter9.html">13. Neural networks</a></li>
<li class="toctree-l1"><a class="reference internal" href="chapter10.html">14. Building a Feed Forward Neural Network</a></li>
<li class="toctree-l1"><a class="reference internal" href="chapter11.html">15. Solving Differential Equations with Deep Learning</a></li>
<li class="toctree-l1"><a class="reference internal" href="chapter12.html">16. Convolutional Neural Networks</a></li>
<li class="toctree-l1"><a class="reference internal" href="chapter13.html">17. Recurrent neural networks: Overarching view</a></li>
</ul>
<p aria-level="2" class="caption" role="heading"><span class="caption-text">Weekly material, notes and exercises</span></p>
<ul class="nav bd-sidenav">
<li class="toctree-l1"><a class="reference internal" href="exercisesweek34.html">Exercises week 34</a></li>
<li class="toctree-l1"><a class="reference internal" href="week34.html">Week 34: Introduction to the course, Logistics and Practicalities</a></li>
<li class="toctree-l1"><a class="reference internal" href="exercisesweek35.html">Exercises week 35</a></li>
<li class="toctree-l1"><a class="reference internal" href="week35.html">Week 35: From Ordinary Linear Regression to Ridge and Lasso Regression</a></li>
<li class="toctree-l1"><a class="reference internal" href="exercisesweek36.html">Exercises week 36</a></li>
<li class="toctree-l1"><a class="reference internal" href="week36.html">Week 36: Linear Regression and Gradient descent</a></li>
<li class="toctree-l1"><a class="reference internal" href="exercisesweek37.html">Exercises week 37</a></li>
<li class="toctree-l1"><a class="reference internal" href="week37.html">Week 37: Gradient descent methods</a></li>
<li class="toctree-l1"><a class="reference internal" href="exercisesweek38.html">Exercises week 38</a></li>
<li class="toctree-l1"><a class="reference internal" href="week38.html">Week 38: Statistical analysis, bias-variance tradeoff and resampling methods</a></li>
<li class="toctree-l1"><a class="reference internal" href="exercisesweek39.html">Exercises week 39</a></li>
<li class="toctree-l1"><a class="reference internal" href="week39.html">Week 39: Resampling methods and logistic regression</a></li>
<li class="toctree-l1"><a class="reference internal" href="week40.html">Week 40: Gradient descent methods (continued) and start Neural networks</a></li>
<li class="toctree-l1"><a class="reference internal" href="week41.html">Week 41 Neural networks and constructing a neural network code</a></li>
<li class="toctree-l1"><a class="reference internal" href="exercisesweek41.html">Exercises week 41</a></li>
</ul>
<p aria-level="2" class="caption" role="heading"><span class="caption-text">Projects</span></p>
<ul class="nav bd-sidenav">
<li class="toctree-l1"><a class="reference internal" href="project1.html">Project 1 on Machine Learning, deadline October 6 (midnight), 2025</a></li>
</ul>
</div>
</nav></div>
</div>
<div class="sidebar-primary-items__end sidebar-primary__section">
</div>
<div id="rtd-footer-container"></div>
</div>
<main id="main-content" class="bd-main" role="main">
<div class="sbt-scroll-pixel-helper"></div>
<div class="bd-content">
<div class="bd-article-container">
<div class="bd-header-article d-print-none">
<div class="header-article-items header-article__inner">
<div class="header-article-items__start">
<div class="header-article-item"><button class="sidebar-toggle primary-toggle btn btn-sm" title="Toggle primary sidebar" data-bs-placement="bottom" data-bs-toggle="tooltip">
<span class="fa-solid fa-bars"></span>
</button></div>
</div>
<div class="header-article-items__end">
<div class="header-article-item">
<div class="article-header-buttons">
<div class="dropdown dropdown-download-buttons">
<button class="btn dropdown-toggle" type="button" data-bs-toggle="dropdown" aria-expanded="false" aria-label="Download this page">
<i class="fas fa-download"></i>
</button>
<ul class="dropdown-menu">
<li><a href="_sources/statistics.ipynb" target="_blank"
class="btn btn-sm btn-download-source-button dropdown-item"
title="Download source file"
data-bs-placement="left" data-bs-toggle="tooltip"
>
<span class="btn__icon-container">
<i class="fas fa-file"></i>
</span>
<span class="btn__text-container">.ipynb</span>
</a>
</li>
<li>
<button onclick="window.print()"
class="btn btn-sm btn-download-pdf-button dropdown-item"
title="Print to PDF"
data-bs-placement="left" data-bs-toggle="tooltip"
>
<span class="btn__icon-container">
<i class="fas fa-file-pdf"></i>
</span>
<span class="btn__text-container">.pdf</span>
</button>
</li>
</ul>
</div>
<button onclick="toggleFullScreen()"
class="btn btn-sm btn-fullscreen-button"
title="Fullscreen mode"
data-bs-placement="bottom" data-bs-toggle="tooltip"
>
<span class="btn__icon-container">
<i class="fas fa-expand"></i>
</span>
</button>
<script>
document.write(`
<button class="btn btn-sm nav-link pst-navbar-icon theme-switch-button" title="light/dark" aria-label="light/dark" data-bs-placement="bottom" data-bs-toggle="tooltip">
<i class="theme-switch fa-solid fa-sun fa-lg" data-mode="light"></i>
<i class="theme-switch fa-solid fa-moon fa-lg" data-mode="dark"></i>
<i class="theme-switch fa-solid fa-circle-half-stroke fa-lg" data-mode="auto"></i>
</button>
`);
</script>
<script>
document.write(`
<button class="btn btn-sm pst-navbar-icon search-button search-button__button" title="Search" aria-label="Search" data-bs-placement="bottom" data-bs-toggle="tooltip">
<i class="fa-solid fa-magnifying-glass fa-lg"></i>
</button>
`);
</script>
<button class="sidebar-toggle secondary-toggle btn btn-sm" title="Toggle secondary sidebar" data-bs-placement="bottom" data-bs-toggle="tooltip">
<span class="fa-solid fa-list"></span>
</button>
</div></div>
</div>
</div>
</div>
<div id="jb-print-docs-body" class="onlyprint">
<h1>Elements of Probability Theory and Statistical Data Analysis</h1>
<!-- Table of contents -->
<div id="print-main-content">
<div id="jb-print-toc">
<div>
<h2> Contents </h2>
</div>
<nav aria-label="Page">
<ul class="visible nav section-nav flex-column">
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#domains-and-probabilities">1.1. Domains and probabilities</a><ul class="nav section-nav flex-column">
<li class="toc-h3 nav-item toc-entry"><a class="reference internal nav-link" href="#stochastic-variables-and-the-main-concepts-the-discrete-case">1.1.1. Stochastic variables and the main concepts, the discrete case</a></li>
<li class="toc-h3 nav-item toc-entry"><a class="reference internal nav-link" href="#properties-of-pdfs">1.1.2. Properties of PDFs</a></li>
<li class="toc-h3 nav-item toc-entry"><a class="reference internal nav-link" href="#expectation-values">1.1.3. Expectation values</a></li>
<li class="toc-h3 nav-item toc-entry"><a class="reference internal nav-link" href="#probability-distribution-functions">1.1.4. Probability Distribution Functions</a></li>
<li class="toc-h3 nav-item toc-entry"><a class="reference internal nav-link" href="#meet-the-covariance">1.1.5. Meet the covariance!</a></li>
<li class="toc-h3 nav-item toc-entry"><a class="reference internal nav-link" href="#numerical-experiments-and-the-covariance-central-limit-theorem">1.1.6. Numerical experiments and the covariance, central limit theorem</a></li>
<li class="toc-h3 nav-item toc-entry"><a class="reference internal nav-link" href="#random-numbers">1.1.7. Random Numbers</a></li>
<li class="toc-h3 nav-item toc-entry"><a class="reference internal nav-link" href="#autocorrelation-function">1.1.8. Autocorrelation function</a></li>
</ul>
</li>
</ul>
</nav>
</div>
</div>
</div>
<div id="searchbox"></div>
<article class="bd-article">
<!-- HTML file automatically generated from DocOnce source (https://github.com/doconce/doconce/)
doconce format html statistics.do.txt --><section class="tex2jax_ignore mathjax_ignore" id="elements-of-probability-theory-and-statistical-data-analysis">
<h1><span class="section-number">1. </span>Elements of Probability Theory and Statistical Data Analysis<a class="headerlink" href="#elements-of-probability-theory-and-statistical-data-analysis" title="Link to this heading">#</a></h1>
<section id="domains-and-probabilities">
<h2><span class="section-number">1.1. </span>Domains and probabilities<a class="headerlink" href="#domains-and-probabilities" title="Link to this heading">#</a></h2>
<p>Consider the following simple example, namely the tossing of two dice, resulting in the following possible values</p>
<div class="math notranslate nohighlight">
\[
\{2,3,4,5,6,7,8,9,10,11,12\}.
\]</div>
<p>These values are called the <em>domain</em>.
To this domain we have the corresponding <em>probabilities</em></p>
<div class="math notranslate nohighlight">
\[
\{1/36,2/36/,3/36,4/36,5/36,6/36,5/36,4/36,3/36,2/36,1/36\}.
\]</div>
<p>The numbers in the domain are the outcomes of the physical process of tossing say two dice.
We cannot tell beforehand whether the outcome is 3 or 5 or any other number in this domain.
This defines the randomness of the outcome, or unexpectedness or any other synonimous word which
encompasses the uncertitude of the final outcome.</p>
<p>The only thing we can tell beforehand
is that say the outcome 2 has a certain probability.<br />
If our favorite hobby is to spend an hour every evening throwing dice and
registering the sequence of outcomes, we will note that the numbers in the above domain</p>
<div class="math notranslate nohighlight">
\[
\{2,3,4,5,6,7,8,9,10,11,12\},
\]</div>
<p>appear in a random order. After 11 throws the results may look like</p>
<div class="math notranslate nohighlight">
\[
\{10,8,6,3,6,9,11,8,12,4,5\}.
\]</div>
<p><strong>Random variables are characterized by a domain which contains all possible values that the random value may take. This domain has a corresponding probability distribution function(PDF)</strong>.</p>
<section id="stochastic-variables-and-the-main-concepts-the-discrete-case">
<h3><span class="section-number">1.1.1. </span>Stochastic variables and the main concepts, the discrete case<a class="headerlink" href="#stochastic-variables-and-the-main-concepts-the-discrete-case" title="Link to this heading">#</a></h3>
<p>There are two main concepts associated with a stochastic variable. The
<em>domain</em> is the set <span class="math notranslate nohighlight">\(\mathbb D = \{x\}\)</span> of all accessible values
the variable can assume, so that <span class="math notranslate nohighlight">\(X \in \mathbb D\)</span>. An example of a
discrete domain is the set of six different numbers that we may get by
throwing of a dice, <span class="math notranslate nohighlight">\(x\in\{1,\,2,\,3,\,4,\,5,\,6\}\)</span>.</p>
<p>The <em>probability distribution function (PDF)</em> is a function
<span class="math notranslate nohighlight">\(p(x)\)</span> on the domain which, in the discrete case, gives us the
probability or relative frequency with which these values of <span class="math notranslate nohighlight">\(X\)</span>
occur</p>
<div class="math notranslate nohighlight">
\[
p(x) = \mathrm{Prob}(X=x).
\]</div>
<p>In the continuous case, the PDF does not directly depict the
actual probability. Instead we define the probability for the
stochastic variable to assume any value on an infinitesimal interval
around <span class="math notranslate nohighlight">\(x\)</span> to be <span class="math notranslate nohighlight">\(p(x)dx\)</span>. The continuous function <span class="math notranslate nohighlight">\(p(x)\)</span> then gives us
the <em>density</em> of the probability rather than the probability
itself. The probability for a stochastic variable to assume any value
on a non-infinitesimal interval <span class="math notranslate nohighlight">\([a,\,b]\)</span> is then just the integral</p>
<div class="math notranslate nohighlight">
\[
\mathrm{Prob}(a\leq X\leq b) = \int_a^b p(x)dx.
\]</div>
<p>Qualitatively speaking, a stochastic variable represents the values of
numbers chosen as if by chance from some specified PDF so that the
selection of a large set of these numbers reproduces this PDF.</p>
<p>Of interest to us is the <em>cumulative probability
distribution function</em> (<strong>CDF</strong>), <span class="math notranslate nohighlight">\(P(x)\)</span>, which is just the probability
for a stochastic variable <span class="math notranslate nohighlight">\(X\)</span> to assume any value less than <span class="math notranslate nohighlight">\(x\)</span></p>
<div class="math notranslate nohighlight">
\[
P(x)=\mathrm{Prob(}X\leq x\mathrm{)} =
\int_{-\infty}^x p(x^{\prime})dx^{\prime}.
\]</div>
<p>The relation between a CDF and its corresponding PDF is then</p>
<div class="math notranslate nohighlight">
\[
p(x) = \frac{d}{dx}P(x).
\]</div>
</section>
<section id="properties-of-pdfs">
<h3><span class="section-number">1.1.2. </span>Properties of PDFs<a class="headerlink" href="#properties-of-pdfs" title="Link to this heading">#</a></h3>
<p>There are two properties that all PDFs must satisfy. The first one is
positivity (assuming that the PDF is normalized)</p>
<div class="math notranslate nohighlight">
\[
0 \leq p(x) \leq 1.
\]</div>
<p>Naturally, it would be nonsensical for any of the values of the domain
to occur with a probability greater than <span class="math notranslate nohighlight">\(1\)</span> or less than <span class="math notranslate nohighlight">\(0\)</span>. Also,
the PDF must be normalized. That is, all the probabilities must add up
to unity. The probability of “anything” to happen is always unity. For
both discrete and continuous PDFs, this condition is</p>
<div class="math notranslate nohighlight">
\[\begin{split}
\begin{align*}
\sum_{x_i\in\mathbb D} p(x_i) &amp; = 1,\\
\int_{x\in\mathbb D} p(x)\,dx &amp; = 1.
\end{align*}
\end{split}\]</div>
<p>The first one
is the most basic PDF; namely the uniform distribution</p>
<!-- Equation labels as ordinary links -->
<div id="eq:unifromPDF"></div>
<div class="math notranslate nohighlight">
\[
\begin{equation}
p(x) = \frac{1}{b-a}\theta(x-a)\theta(b-x).
\label{eq:unifromPDF} \tag{1}
\end{equation}
\]</div>
<p>For <span class="math notranslate nohighlight">\(a=0\)</span> and <span class="math notranslate nohighlight">\(b=1\)</span> we have</p>
<div class="math notranslate nohighlight">
\[
\begin{array}{ll}
p(x)dx = dx &amp; \in [0,1].
\end{array}
\]</div>
<p>The latter distribution is used to generate random numbers. For other PDFs, one needs normally a mapping from this distribution to say for example the exponential distribution.</p>
<p>The second one is the Gaussian Distribution</p>
<div class="math notranslate nohighlight">
\[
p(x) = \frac{1}{\sigma\sqrt{2\pi}} \exp{(-\frac{(x-\mu)^2}{2\sigma^2})},
\]</div>
<p>with mean value <span class="math notranslate nohighlight">\(\mu\)</span> and standard deviation <span class="math notranslate nohighlight">\(\sigma\)</span>. If <span class="math notranslate nohighlight">\(\mu=0\)</span> and <span class="math notranslate nohighlight">\(\sigma=1\)</span>, it is normally called the <strong>standard normal distribution</strong></p>
<div class="math notranslate nohighlight">
\[
p(x) = \frac{1}{\sqrt{2\pi}} \exp{(-\frac{x^2}{2})},
\]</div>
<p>The following simple Python code plots the above distribution for different values of <span class="math notranslate nohighlight">\(\mu\)</span> and <span class="math notranslate nohighlight">\(\sigma\)</span>.</p>
<div class="cell docutils container">
<div class="cell_input docutils container">
<div class="highlight-none notranslate"><div class="highlight"><pre><span></span>%matplotlib inline
import numpy as np
from math import acos, exp, sqrt
from matplotlib import pyplot as plt
from matplotlib import rc, rcParams
import matplotlib.units as units
import matplotlib.ticker as ticker
rc(&#39;text&#39;,usetex=True)
rc(&#39;font&#39;,**{&#39;family&#39;:&#39;serif&#39;,&#39;serif&#39;:[&#39;Gaussian distribution&#39;]})
font = {&#39;family&#39; : &#39;serif&#39;,
&#39;color&#39; : &#39;darkred&#39;,
&#39;weight&#39; : &#39;normal&#39;,
&#39;size&#39; : 16,
}
pi = acos(-1.0)
mu0 = 0.0
sigma0 = 1.0
mu1= 1.0
sigma1 = 2.0
mu2 = 2.0
sigma2 = 4.0
x = np.linspace(-20.0, 20.0)
v0 = np.exp(-(x*x-2*x*mu0+mu0*mu0)/(2*sigma0*sigma0))/sqrt(2*pi*sigma0*sigma0)
v1 = np.exp(-(x*x-2*x*mu1+mu1*mu1)/(2*sigma1*sigma1))/sqrt(2*pi*sigma1*sigma1)
v2 = np.exp(-(x*x-2*x*mu2+mu2*mu2)/(2*sigma2*sigma2))/sqrt(2*pi*sigma2*sigma2)
plt.plot(x, v0, &#39;b-&#39;, x, v1, &#39;r-&#39;, x, v2, &#39;g-&#39;)
plt.title(r&#39;{\bf Gaussian distributions}&#39;, fontsize=20)
plt.text(-19, 0.3, r&#39;Parameters: $\mu = 0$, $\sigma = 1$&#39;, fontdict=font)
plt.text(-19, 0.18, r&#39;Parameters: $\mu = 1$, $\sigma = 2$&#39;, fontdict=font)
plt.text(-19, 0.08, r&#39;Parameters: $\mu = 2$, $\sigma = 4$&#39;, fontdict=font)
plt.xlabel(r&#39;$x$&#39;,fontsize=20)
plt.ylabel(r&#39;$p(x)$ [MeV]&#39;,fontsize=20)
# Tweak spacing to prevent clipping of ylabel
plt.subplots_adjust(left=0.15)
plt.savefig(&#39;gaussian.pdf&#39;, format=&#39;pdf&#39;)
plt.show()
</pre></div>
</div>
</div>
</div>
<p>Another important distribution in science is the exponential distribution</p>
<div class="math notranslate nohighlight">
\[
p(x) = \alpha\exp{-(\alpha x)}.
\]</div>
</section>
<section id="expectation-values">
<h3><span class="section-number">1.1.3. </span>Expectation values<a class="headerlink" href="#expectation-values" title="Link to this heading">#</a></h3>
<p>Let <span class="math notranslate nohighlight">\(h(x)\)</span> be an arbitrary continuous function on the domain of the stochastic
variable <span class="math notranslate nohighlight">\(X\)</span> whose PDF is <span class="math notranslate nohighlight">\(p(x)\)</span>. We define the <em>expectation value</em>
of <span class="math notranslate nohighlight">\(h\)</span> with respect to <span class="math notranslate nohighlight">\(p\)</span> as follows</p>
<!-- Equation labels as ordinary links -->
<div id="eq:expectation_value_of_h_wrt_p"></div>
<div class="math notranslate nohighlight">
\[
\begin{equation}
\langle h \rangle_X \equiv \int\! h(x)p(x)\,dx
\label{eq:expectation_value_of_h_wrt_p} \tag{2}
\end{equation}
\]</div>
<p>Whenever the PDF is known implicitly, like in this case, we will drop
the index <span class="math notranslate nohighlight">\(X\)</span> for clarity.<br />
A particularly useful class of special expectation values are the
<em>moments</em>. The <span class="math notranslate nohighlight">\(n\)</span>-th moment of the PDF <span class="math notranslate nohighlight">\(p\)</span> is defined as
follows</p>
<div class="math notranslate nohighlight">
\[
\langle x^n \rangle \equiv \int\! x^n p(x)\,dx
\]</div>
<p>The zero-th moment <span class="math notranslate nohighlight">\(\langle 1\rangle\)</span> is just the normalization condition of
<span class="math notranslate nohighlight">\(p\)</span>. The first moment, <span class="math notranslate nohighlight">\(\langle x\rangle\)</span>, is called the <em>mean</em> of <span class="math notranslate nohighlight">\(p\)</span>
and often denoted by the letter <span class="math notranslate nohighlight">\(\mu\)</span></p>
<div class="math notranslate nohighlight">
\[
\langle x\rangle = \mu \equiv \int x p(x)dx,
\]</div>
<p>for a continuous distribution and</p>
<div class="math notranslate nohighlight">
\[
\langle x\rangle = \mu \equiv \sum_{i=1}^N x_i p(x_i),
\]</div>
<p>for a discrete distribution.
Qualitatively it represents the centroid or the average value of the
PDF and is therefore simply called the expectation value of <span class="math notranslate nohighlight">\(p(x)\)</span>.</p>
<p>A special version of the moments is the set of <em>central moments</em>, the n-th central moment defined as</p>
<div class="math notranslate nohighlight">
\[
\langle (x-\langle x\rangle )^n\rangle \equiv \int\! (x-\langle x\rangle)^n p(x)\,dx
\]</div>
<p>The zero-th and first central moments are both trivial, equal <span class="math notranslate nohighlight">\(1\)</span> and
<span class="math notranslate nohighlight">\(0\)</span>, respectively. But the second central moment, known as the
<em>variance</em> of <span class="math notranslate nohighlight">\(p\)</span>, is of particular interest. For the stochastic
variable <span class="math notranslate nohighlight">\(X\)</span>, the variance is denoted as <span class="math notranslate nohighlight">\(\sigma^2_X\)</span> or <span class="math notranslate nohighlight">\(\mathrm{Var}(X)\)</span></p>
<div class="math notranslate nohighlight">
\[\begin{split}
\begin{align*}
\sigma^2_X &amp;=\mathrm{Var}(X) = \langle (x-\langle x\rangle)^2\rangle =
\int (x-\langle x\rangle)^2 p(x)dx\\
&amp; = \int\left(x^2 - 2 x \langle x\rangle^{2} +\langle x\rangle^2\right)p(x)dx\\
&amp; = \langle x^2\rangle - 2 \langle x\rangle\langle x\rangle + \langle x\rangle^2\\
&amp; = \langle x^2 \rangle - \langle x\rangle^2
\end{align*}
\end{split}\]</div>
<p>The square root of the variance, <span class="math notranslate nohighlight">\(\sigma =\sqrt{\langle (x-\langle x\rangle)^2\rangle}\)</span> is called the
<strong>standard deviation</strong> of <span class="math notranslate nohighlight">\(p\)</span>. It is the RMS (root-mean-square)
value of the deviation of the PDF from its mean value, interpreted
qualitatively as the “spread” of <span class="math notranslate nohighlight">\(p\)</span> around its mean.</p>
</section>
<section id="probability-distribution-functions">
<h3><span class="section-number">1.1.4. </span>Probability Distribution Functions<a class="headerlink" href="#probability-distribution-functions" title="Link to this heading">#</a></h3>
<p>The following table collects properties of probability distribution functions.
In our notation we reserve the label <span class="math notranslate nohighlight">\(p(x)\)</span> for the probability of a certain event,
while <span class="math notranslate nohighlight">\(P(x)\)</span> is the cumulative probability.</p>
<table class="dotable" border="1">
<thead>
<tr><th align="center"> </th> <th align="center"> Discrete PDF </th> <th align="center"> Continuous PDF </th> </tr>
</thead>
<tbody>
<tr><td align="left"> Domain </td> <td align="center"> $\left\{x_1, x_2, x_3, \dots, x_N\right\}$ </td> <td align="center"> $[a,b]$ </td> </tr>
<tr><td align="left"> Probability </td> <td align="center"> $p(x_i)$ </td> <td align="center"> $p(x)dx$ </td> </tr>
<tr><td align="left"> Cumulative </td> <td align="center"> $P_i=\sum_{l=1}^ip(x_l)$ </td> <td align="center"> $P(x)=\int_a^xp(t)dt$ </td> </tr>
<tr><td align="left"> Positivity </td> <td align="center"> $0 \le p(x_i) \le 1$ </td> <td align="center"> $p(x) \ge 0$ </td> </tr>
<tr><td align="left"> Positivity </td> <td align="center"> $0 \le P_i \le 1$ </td> <td align="center"> $0 \le P(x) \le 1$ </td> </tr>
<tr><td align="left"> Monotonic </td> <td align="center"> $P_i \ge P_j$ if $x_i \ge x_j$ </td> <td align="center"> $P(x_i) \ge P(x_j)$ if $x_i \ge x_j$ </td> </tr>
<tr><td align="left"> Normalization </td> <td align="center"> $P_N=1$ </td> <td align="center"> $P(b)=1$ </td> </tr>
</tbody>
</table>
<p>With a PDF we can compute expectation values of selected quantities such as</p>
<div class="math notranslate nohighlight">
\[
\langle x^k\rangle=\sum_{i=1}^{N}x_i^kp(x_i),
\]</div>
<p>if we have a discrete PDF or</p>
<div class="math notranslate nohighlight">
\[
\langle x^k\rangle=\int_a^b x^kp(x)dx,
\]</div>
<p>in the case of a continuous PDF. We have already defined the mean value <span class="math notranslate nohighlight">\(\mu\)</span>
and the variance <span class="math notranslate nohighlight">\(\sigma^2\)</span>.</p>
<p>There are at least three PDFs which one may encounter. These are the</p>
<p><strong>Uniform distribution</strong></p>
<div class="math notranslate nohighlight">
\[
p(x)=\frac{1}{b-a}\Theta(x-a)\Theta(b-x),
\]</div>
<p>yielding probabilities different from zero in the interval <span class="math notranslate nohighlight">\([a,b]\)</span>.</p>
<p><strong>The exponential distribution</strong></p>
<div class="math notranslate nohighlight">
\[
p(x)=\alpha \exp{(-\alpha x)},
\]</div>
<p>yielding probabilities different from zero in the interval <span class="math notranslate nohighlight">\([0,\infty)\)</span> and with mean value</p>
<div class="math notranslate nohighlight">
\[
\mu = \int_0^{\infty}xp(x)dx=\int_0^{\infty}x\alpha \exp{(-\alpha x)}dx=\frac{1}{\alpha},
\]</div>
<p>with variance</p>
<div class="math notranslate nohighlight">
\[
\sigma^2=\int_0^{\infty}x^2p(x)dx-\mu^2 = \frac{1}{\alpha^2}.
\]</div>
<p>Finally, we have the so-called univariate normal distribution, or just the <strong>normal distribution</strong></p>
<div class="math notranslate nohighlight">
\[
p(x)=\frac{1}{b\sqrt{2\pi}}\exp{\left(-\frac{(x-a)^2}{2b^2}\right)}
\]</div>
<p>with probabilities different from zero in the interval <span class="math notranslate nohighlight">\((-\infty,\infty)\)</span>.
The integral <span class="math notranslate nohighlight">\(\int_{-\infty}^{\infty}\exp{\left(-(x^2\right)}dx\)</span> appears in many calculations, its value
is <span class="math notranslate nohighlight">\(\sqrt{\pi}\)</span>, a result we will need when we compute the mean value and the variance.
The mean value is</p>
<div class="math notranslate nohighlight">
\[
\mu = \int_0^{\infty}xp(x)dx=\frac{1}{b\sqrt{2\pi}}\int_{-\infty}^{\infty}x \exp{\left(-\frac{(x-a)^2}{2b^2}\right)}dx,
\]</div>
<p>which becomes with a suitable change of variables</p>
<div class="math notranslate nohighlight">
\[
\mu =\frac{1}{b\sqrt{2\pi}}\int_{-\infty}^{\infty}b\sqrt{2}(a+b\sqrt{2}y)\exp{-y^2}dy=a.
\]</div>
<p>Similarly, the variance becomes</p>
<div class="math notranslate nohighlight">
\[
\sigma^2 = \frac{1}{b\sqrt{2\pi}}\int_{-\infty}^{\infty}(x-\mu)^2 \exp{\left(-\frac{(x-a)^2}{2b^2}\right)}dx,
\]</div>
<p>and inserting the mean value and performing a variable change we obtain</p>
<div class="math notranslate nohighlight">
\[
\sigma^2 = \frac{1}{b\sqrt{2\pi}}\int_{-\infty}^{\infty}b\sqrt{2}(b\sqrt{2}y)^2\exp{\left(-y^2\right)}dy=
\frac{2b^2}{\sqrt{\pi}}\int_{-\infty}^{\infty}y^2\exp{\left(-y^2\right)}dy,
\]</div>
<p>and performing a final integration by parts we obtain the well-known result <span class="math notranslate nohighlight">\(\sigma^2=b^2\)</span>.
It is useful to introduce the standard normal distribution as well, defined by <span class="math notranslate nohighlight">\(\mu=a=0\)</span>, viz. a distribution
centered around zero and with a variance <span class="math notranslate nohighlight">\(\sigma^2=1\)</span>, leading to</p>
<!-- Equation labels as ordinary links -->
<div id="_auto1"></div>
<div class="math notranslate nohighlight">
\[
\begin{equation}
p(x)=\frac{1}{\sqrt{2\pi}}\exp{\left(-\frac{x^2}{2}\right)}.
\label{_auto1} \tag{3}
\end{equation}
\]</div>
<p>The exponential and uniform distributions have simple cumulative functions,
whereas the normal distribution does not, being proportional to the so-called
error function <span class="math notranslate nohighlight">\(erf(x)\)</span>, given by</p>
<div class="math notranslate nohighlight">
\[
P(x) = \frac{1}{\sqrt{2\pi}}\int_{-\infty}^x\exp{\left(-\frac{t^2}{2}\right)}dt,
\]</div>
<p>which is difficult to evaluate in a quick way.</p>
<p>Some other PDFs which one encounters often in the natural sciences are the binomial distribution</p>
<div class="math notranslate nohighlight">
\[\begin{split}
p(x) = \left(\begin{array}{c} n \\ x\end{array}\right)y^x(1-y)^{n-x} \hspace{0.5cm}x=0,1,\dots,n,
\end{split}\]</div>
<p>where <span class="math notranslate nohighlight">\(y\)</span> is the probability for a specific event, such as the tossing of a coin or moving left or right
in case of a random walker. Note that <span class="math notranslate nohighlight">\(x\)</span> is a discrete stochastic variable.</p>
<p>The sequence of binomial trials is characterized by the following definitions</p>
<ul class="simple">
<li><p>Every experiment is thought to consist of <span class="math notranslate nohighlight">\(N\)</span> independent trials.</p></li>
<li><p>In every independent trial one registers if a specific situation happens or not, such as the jump to the left or right of a random walker.</p></li>
<li><p>The probability for every outcome in a single trial has the same value, for example the outcome of tossing (either heads or tails) a coin is always <span class="math notranslate nohighlight">\(1/2\)</span>.</p></li>
</ul>
<p>In order to compute the mean and variance we need to recall Newton’s binomial
formula</p>
<div class="math notranslate nohighlight">
\[\begin{split}
(a+b)^m=\sum_{n=0}^m \left(\begin{array}{c} m \\ n\end{array}\right)a^nb^{m-n},
\end{split}\]</div>
<p>which can be used to show that</p>
<div class="math notranslate nohighlight">
\[\begin{split}
\sum_{x=0}^n\left(\begin{array}{c} n \\ x\end{array}\right)y^x(1-y)^{n-x} = (y+1-y)^n = 1,
\end{split}\]</div>
<p>the PDF is normalized to one.
The mean value is</p>
<div class="math notranslate nohighlight">
\[\begin{split}
\mu = \sum_{x=0}^n x\left(\begin{array}{c} n \\ x\end{array}\right)y^x(1-y)^{n-x} =
\sum_{x=0}^n x\frac{n!}{x!(n-x)!}y^x(1-y)^{n-x},
\end{split}\]</div>
<p>resulting in</p>
<div class="math notranslate nohighlight">
\[
\mu =
\sum_{x=0}^n x\frac{(n-1)!}{(x-1)!(n-1-(x-1))!}y^{x-1}(1-y)^{n-1-(x-1)},
\]</div>
<p>which we rewrite as</p>
<div class="math notranslate nohighlight">
\[\begin{split}
\mu=ny\sum_{\nu=0}^n\left(\begin{array}{c} n-1 \\ \nu\end{array}\right)y^{\nu}(1-y)^{n-1-\nu} =ny(y+1-y)^{n-1}=ny.
\end{split}\]</div>
<p>The variance is slightly trickier to get. It reads <span class="math notranslate nohighlight">\(\sigma^2=ny(1-y)\)</span>.</p>
<p>Another important distribution with discrete stochastic variables <span class="math notranslate nohighlight">\(x\)</span> is<br />
the Poisson model, which resembles the exponential distribution and reads</p>
<div class="math notranslate nohighlight">
\[
p(x) = \frac{\lambda^x}{x!} e^{-\lambda} \hspace{0.5cm}x=0,1,\dots,;\lambda &gt; 0.
\]</div>
<p>In this case both the mean value and the variance are easier to calculate,</p>
<div class="math notranslate nohighlight">
\[
\mu = \sum_{x=0}^{\infty} x \frac{\lambda^x}{x!} e^{-\lambda} = \lambda e^{-\lambda}\sum_{x=1}^{\infty}
\frac{\lambda^{x-1}}{(x-1)!}=\lambda,
\]</div>
<p>and the variance is <span class="math notranslate nohighlight">\(\sigma^2=\lambda\)</span>.</p>
<p>An example of applications of the Poisson distribution could be the counting
of the number of <span class="math notranslate nohighlight">\(\alpha\)</span>-particles emitted from a radioactive source in a given time interval.
In the limit of <span class="math notranslate nohighlight">\(n\rightarrow \infty\)</span> and for small probabilities <span class="math notranslate nohighlight">\(y\)</span>, the binomial distribution
approaches the Poisson distribution. Setting <span class="math notranslate nohighlight">\(\lambda = ny\)</span>, with <span class="math notranslate nohighlight">\(y\)</span> the probability for an event in
the binomial distribution we can show that</p>
<div class="math notranslate nohighlight">
\[\begin{split}
\lim_{n\rightarrow \infty}\left(\begin{array}{c} n \\ x\end{array}\right)y^x(1-y)^{n-x} e^{-\lambda}=\sum_{x=1}^{\infty}\frac{\lambda^x}{x!} e^{-\lambda}.
\end{split}\]</div>
</section>
<section id="meet-the-covariance">
<h3><span class="section-number">1.1.5. </span>Meet the covariance!<a class="headerlink" href="#meet-the-covariance" title="Link to this heading">#</a></h3>
<p>An important quantity in a statistical analysis is the so-called covariance.</p>
<p>Consider the set <span class="math notranslate nohighlight">\(\{X_i\}\)</span> of <span class="math notranslate nohighlight">\(n\)</span>
stochastic variables (not necessarily uncorrelated) with the
multivariate PDF <span class="math notranslate nohighlight">\(P(x_1,\dots,x_n)\)</span>. The <em>covariance</em> of two
of the stochastic variables, <span class="math notranslate nohighlight">\(X_i\)</span> and <span class="math notranslate nohighlight">\(X_j\)</span>, is defined as follows</p>
<!-- Equation labels as ordinary links -->
<div id="_auto2"></div>
<div class="math notranslate nohighlight">
\[
\begin{equation}
\mathrm{Cov}(X_i,\,X_j) = \langle (x_i-\langle x_i\rangle)(x_j-\langle x_j\rangle)\rangle
\label{_auto2} \tag{4}
\end{equation}
\]</div>
<!-- Equation labels as ordinary links -->
<div id="eq:def_covariance"></div>
<div class="math notranslate nohighlight">
\[
\begin{equation}
=\int\cdots\int (x_i-\langle x_i\rangle)(x_j-\langle x_j\rangle)P(x_1,\dots,x_n)\,dx_1\dots dx_n,
\label{eq:def_covariance} \tag{5}
\end{equation}
\]</div>
<p>with</p>
<div class="math notranslate nohighlight">
\[
\langle x_i\rangle =
\int\cdots\int x_i P(x_1,\dots,x_n)\,dx_1\dots dx_n.
\]</div>
<p>If we consider the above covariance as a matrix</p>
<div class="math notranslate nohighlight">
\[
C_{ij} =\mathrm{Cov}(X_i,\,X_j),
\]</div>
<p>then the diagonal elements are just the familiar
variances, <span class="math notranslate nohighlight">\(C_{ii} = \mathrm{Cov}(X_i,\,X_i) = \mathrm{Var}(X_i)\)</span>. It turns out that
all the off-diagonal elements are zero if the stochastic variables are
uncorrelated.</p>
<div class="cell docutils container">
<div class="cell_input docutils container">
<div class="highlight-none notranslate"><div class="highlight"><pre><span></span># Importing various packages
from math import exp, sqrt
from random import random, seed
import numpy as np
import matplotlib.pyplot as plt
def covariance(x, y, n):
sum = 0.0
mean_x = np.mean(x)
mean_y = np.mean(y)
for i in range(0, n):
sum += (x[(i)]-mean_x)*(y[i]-mean_y)
return sum/n
n = 10
x=np.random.normal(size=n)
y = 4+3*x+np.random.normal(size=n)
covxy = covariance(x,y,n)
print(covxy)
z = np.vstack((x, y))
c = np.cov(z.T)
print(c)
</pre></div>
</div>
</div>
</div>
<p>Consider the stochastic variables <span class="math notranslate nohighlight">\(X_i\)</span> and <span class="math notranslate nohighlight">\(X_j\)</span>, (<span class="math notranslate nohighlight">\(i\neq j\)</span>). We have</p>
<div class="math notranslate nohighlight">
\[\begin{split}
\begin{align*}
Cov(X_i,\,X_j) &amp;= \langle (x_i-\langle x_i\rangle)(x_j-\langle x_j\rangle)\rangle\\
&amp;=\langle x_i x_j - x_i\langle x_j\rangle - \langle x_i\rangle x_j + \langle x_i\rangle\langle x_j\rangle\rangle\\
&amp;=\langle x_i x_j\rangle - \langle x_i\langle x_j\rangle\rangle - \langle \langle x_i\rangle x_j \rangle +
\langle \langle x_i\rangle\langle x_j\rangle\rangle \\
&amp;=\langle x_i x_j\rangle - \langle x_i\rangle\langle x_j\rangle - \langle x_i\rangle\langle x_j\rangle +
\langle x_i\rangle\langle x_j\rangle \\
&amp;=\langle x_i x_j\rangle - \langle x_i\rangle\langle x_j\rangle
\end{align*}
\end{split}\]</div>
<p>If <span class="math notranslate nohighlight">\(X_i\)</span> and <span class="math notranslate nohighlight">\(X_j\)</span> are independent (assuming <span class="math notranslate nohighlight">\(i \neq j\)</span>), we have that</p>
<div class="math notranslate nohighlight">
\[
\langle x_i x_j\rangle = \langle x_i\rangle\langle x_j\rangle,
\]</div>
<p>leading to</p>
<div class="math notranslate nohighlight">
\[
Cov(X_i, X_j) = 0 \hspace{0.1cm} (i\neq j).
\]</div>
<p>Now that we have constructed an idealized mathematical framework, let
us try to apply it to empirical observations. Examples of relevant
physical phenomena may be spontaneous decays of nuclei, or a purely
mathematical set of numbers produced by some deterministic
mechanism. It is the latter we will deal with, using so-called pseudo-random
number generators. In general our observations will contain only a limited set of
observables. We remind the reader that
a <em>stochastic process</em> is a process that produces sequentially a
chain of values</p>
<div class="math notranslate nohighlight">
\[
\{x_1, x_2,\dots\,x_k,\dots\}.
\]</div>
<p>We will call these
values our <em>measurements</em> and the entire set as our measured
<em>sample</em>. The action of measuring all the elements of a sample
we will call a stochastic <em>experiment</em> (since, operationally,
they are often associated with results of empirical observation of
some physical or mathematical phenomena; precisely an experiment). We
assume that these values are distributed according to some
PDF <span class="math notranslate nohighlight">\(p_X^{\phantom X}(x)\)</span>, where <span class="math notranslate nohighlight">\(X\)</span> is just the formal symbol for the
stochastic variable whose PDF is <span class="math notranslate nohighlight">\(p_X^{\phantom X}(x)\)</span>. Instead of
trying to determine the full distribution <span class="math notranslate nohighlight">\(p\)</span> we are often only
interested in finding the few lowest moments, like the mean
<span class="math notranslate nohighlight">\(\mu_X^{\phantom X}\)</span> and the variance <span class="math notranslate nohighlight">\(\sigma_X^{\phantom X}\)</span>.</p>
<p>In practical situations however, a sample is always of finite size. Let that
size be <span class="math notranslate nohighlight">\(n\)</span>. The expectation value of a sample <span class="math notranslate nohighlight">\(\alpha\)</span>, the <strong>sample mean</strong>, is then defined as follows</p>
<div class="math notranslate nohighlight">
\[
\langle x_{\alpha} \rangle \equiv \frac{1}{n}\sum_{k=1}^n x_{\alpha,k}.
\]</div>
<p>The <em>sample variance</em> is:</p>
<div class="math notranslate nohighlight">
\[
\mathrm{Var}(x) \equiv \frac{1}{n}\sum_{k=1}^n (x_{\alpha,k} - \langle x_{\alpha} \rangle)^2,
\]</div>
<p>with its square root being the <em>standard deviation of the sample</em>.</p>
<p>You can think of the above observables as a set of quantities which define
a given experiment. This experiment is then repeated several times, say <span class="math notranslate nohighlight">\(m\)</span> times.
The total average is then</p>
<!-- Equation labels as ordinary links -->
<div id="eq:exptmean"></div>
<div class="math notranslate nohighlight">
\[
\begin{equation}
\langle X_m \rangle= \frac{1}{m}\sum_{\alpha=1}^mx_{\alpha}=\frac{1}{mn}\sum_{\alpha, k} x_{\alpha,k},
\label{eq:exptmean} \tag{6}
\end{equation}
\]</div>
<p>where the last sums end at <span class="math notranslate nohighlight">\(m\)</span> and <span class="math notranslate nohighlight">\(n\)</span>.
The total variance is</p>
<div class="math notranslate nohighlight">
\[
\sigma^2_m= \frac{1}{mn^2}\sum_{\alpha=1}^m(\langle x_{\alpha} \rangle-\langle X_m \rangle)^2,
\]</div>
<p>which we rewrite as</p>
<!-- Equation labels as ordinary links -->
<div id="eq:exptvariance"></div>
<div class="math notranslate nohighlight">
\[
\begin{equation}
\sigma^2_m=\frac{1}{m}\sum_{\alpha=1}^m\sum_{kl=1}^n (x_{\alpha,k}-\langle X_m \rangle)(x_{\alpha,l}-\langle X_m \rangle).
\label{eq:exptvariance} \tag{7}
\end{equation}
\]</div>
<p>We define also the sample variance <span class="math notranslate nohighlight">\(\sigma^2\)</span> of all <span class="math notranslate nohighlight">\(mn\)</span> individual experiments as</p>
<!-- Equation labels as ordinary links -->
<div id="eq:sampleexptvariance"></div>
<div class="math notranslate nohighlight">
\[
\begin{equation}
\sigma^2=\frac{1}{mn}\sum_{\alpha=1}^m\sum_{k=1}^n (x_{\alpha,k}-\langle X_m \rangle)^2.
\label{eq:sampleexptvariance} \tag{8}
\end{equation}
\]</div>
<p>These quantities, being known experimental values or the results from our calculations,
may differ, in some cases
significantly, from the similarly named
exact values for the mean value <span class="math notranslate nohighlight">\(\mu_X\)</span>, the variance <span class="math notranslate nohighlight">\(\mathrm{Var}(X)\)</span>
and the covariance <span class="math notranslate nohighlight">\(\mathrm{Cov}(X,Y)\)</span>.</p>
</section>
<section id="numerical-experiments-and-the-covariance-central-limit-theorem">
<h3><span class="section-number">1.1.6. </span>Numerical experiments and the covariance, central limit theorem<a class="headerlink" href="#numerical-experiments-and-the-covariance-central-limit-theorem" title="Link to this heading">#</a></h3>
<p>The central limit theorem states that the PDF <span class="math notranslate nohighlight">\(\tilde{p}(z)\)</span> of
the average of <span class="math notranslate nohighlight">\(m\)</span> random values corresponding to a PDF <span class="math notranslate nohighlight">\(p(x)\)</span>
is a normal distribution whose mean is the
mean value of the PDF <span class="math notranslate nohighlight">\(p(x)\)</span> and whose variance is the variance
of the PDF <span class="math notranslate nohighlight">\(p(x)\)</span> divided by <span class="math notranslate nohighlight">\(m\)</span>, the number of values used to compute <span class="math notranslate nohighlight">\(z\)</span>.</p>
<p>The central limit theorem leads then to the well-known expression for the
standard deviation, given by</p>
<div class="math notranslate nohighlight">
\[
\sigma_m=
\frac{\sigma}{\sqrt{m}}.
\]</div>
<p>In many cases the above estimate for the standard deviation, in particular if correlations are strong, may be too simplistic. We need therefore a more precise defintion of the error and the variance in our results.</p>
<p>Our estimate of the true average <span class="math notranslate nohighlight">\(\mu_{X}\)</span> is the sample mean <span class="math notranslate nohighlight">\(\langle X_m \rangle\)</span></p>
<div class="math notranslate nohighlight">
\[
\mu_{X}^{\phantom X} \approx X_m=\frac{1}{mn}\sum_{\alpha=1}^m\sum_{k=1}^n x_{\alpha,k}.
\]</div>
<p>We can then use Eq. (<a class="reference internal" href="#eq:exptvariance"><span class="xref myst">7</span></a>)</p>
<div class="math notranslate nohighlight">
\[
\sigma^2_m=\frac{1}{mn^2}\sum_{\alpha=1}^m\sum_{kl=1}^n (x_{\alpha,k}-\langle X_m \rangle)(x_{\alpha,l}-\langle X_m \rangle),
\]</div>
<p>and rewrite it as</p>
<div class="math notranslate nohighlight">
\[
\sigma^2_m=\frac{\sigma^2}{n}+\frac{2}{mn^2}\sum_{\alpha=1}^m\sum_{k&lt;l}^n (x_{\alpha,k}-\langle X_m \rangle)(x_{\alpha,l}-\langle X_m \rangle),
\]</div>
<p>where the first term is the sample variance of all <span class="math notranslate nohighlight">\(mn\)</span> experiments divided by <span class="math notranslate nohighlight">\(n\)</span>
and the last term is nothing but the covariance which arises when <span class="math notranslate nohighlight">\(k\ne l\)</span>.</p>
<p>Our estimate of the true average <span class="math notranslate nohighlight">\(\mu_{X}\)</span> is the sample mean <span class="math notranslate nohighlight">\(\langle X_m \rangle\)</span></p>
<p>If the
observables are uncorrelated, then the covariance is zero and we obtain a total variance
which agrees with the central limit theorem. Correlations may often be present in our data set, resulting in a non-zero covariance. The first term is normally called the uncorrelated
contribution.
Computationally the uncorrelated first term is much easier to treat
efficiently than the second.
We just accumulate separately the values <span class="math notranslate nohighlight">\(x^2\)</span> and <span class="math notranslate nohighlight">\(x\)</span> for every
measurement <span class="math notranslate nohighlight">\(x\)</span> we receive. The correlation term, though, has to be
calculated at the end of the experiment since we need all the
measurements to calculate the cross terms. Therefore, all measurements
have to be stored throughout the experiment.</p>
<p>Let us analyze the problem by splitting up the correlation term into
partial sums of the form</p>
<div class="math notranslate nohighlight">
\[
f_d = \frac{1}{nm}\sum_{\alpha=1}^m\sum_{k=1}^{n-d}(x_{\alpha,k}-\langle X_m \rangle)(x_{\alpha,k+d}-\langle X_m \rangle),
\]</div>
<p>The correlation term of the total variance can now be rewritten in terms of
<span class="math notranslate nohighlight">\(f_d\)</span></p>
<div class="math notranslate nohighlight">
\[
\frac{2}{mn^2}\sum_{\alpha=1}^m\sum_{k&lt;l}^n (x_{\alpha,k}-\langle X_m \rangle)(x_{\alpha,l}-\langle X_m \rangle)=
\frac{2}{n}\sum_{d=1}^{n-1} f_d
\]</div>
<p>The value of <span class="math notranslate nohighlight">\(f_d\)</span> reflects the correlation between measurements
separated by the distance <span class="math notranslate nohighlight">\(d\)</span> in the samples. Notice that for
<span class="math notranslate nohighlight">\(d=0\)</span>, <span class="math notranslate nohighlight">\(f\)</span> is just the sample variance, <span class="math notranslate nohighlight">\(\sigma^2\)</span>. If we divide <span class="math notranslate nohighlight">\(f_d\)</span>
by <span class="math notranslate nohighlight">\(\sigma^2\)</span>, we arrive at the so called <strong>autocorrelation function</strong></p>
<!-- Equation labels as ordinary links -->
<div id="eq:autocorrelformal"></div>
<div class="math notranslate nohighlight">
\[
\begin{equation}
\kappa_d = \frac{f_d}{\sigma^2}
\label{eq:autocorrelformal} \tag{9}
\end{equation}
\]</div>
<p>which gives us a useful measure of the correlation pair correlation
starting always at <span class="math notranslate nohighlight">\(1\)</span> for <span class="math notranslate nohighlight">\(d=0\)</span>.</p>
<p>The sample variance of the <span class="math notranslate nohighlight">\(mn\)</span> experiments can now be
written in terms of the autocorrelation function</p>
<!-- Equation labels as ordinary links -->
<div id="eq:error_estimate_corr_time"></div>
<div class="math notranslate nohighlight">
\[
\begin{equation}
\sigma_m^2=\frac{\sigma^2}{n}+\frac{2}{n}\cdot\sigma^2\sum_{d=1}^{n-1}
\frac{f_d}{\sigma^2}=\left(1+2\sum_{d=1}^{n-1}\kappa_d\right)\frac{1}{n}\sigma^2=\frac{\tau}{n}\cdot\sigma^2
\label{eq:error_estimate_corr_time} \tag{10}
\end{equation}
\]</div>
<p>and we see that <span class="math notranslate nohighlight">\(\sigma_m\)</span> can be expressed in terms of the
uncorrelated sample variance times a correction factor <span class="math notranslate nohighlight">\(\tau\)</span> which
accounts for the correlation between measurements. We call this
correction factor the <em>autocorrelation time</em></p>
<!-- Equation labels as ordinary links -->
<div id="eq:autocorrelation_time"></div>
<div class="math notranslate nohighlight">
\[
\begin{equation}
\tau = 1+2\sum_{d=1}^{n-1}\kappa_d
\label{eq:autocorrelation_time} \tag{11}
\end{equation}
\]</div>
<!-- It is closely related to the area under the graph of the -->
<!-- autocorrelation function. -->
<p>For a correlation free experiment, <span class="math notranslate nohighlight">\(\tau\)</span>
equals 1.</p>
<p>From the point of view of
Eq. (<a class="reference internal" href="#eq:error_estimate_corr_time"><span class="xref myst">10</span></a>) we can interpret a sequential
correlation as an effective reduction of the number of measurements by
a factor <span class="math notranslate nohighlight">\(\tau\)</span>. The effective number of measurements becomes</p>
<div class="math notranslate nohighlight">
\[
n_\mathrm{eff} = \frac{n}{\tau}
\]</div>
<p>To neglect the autocorrelation time <span class="math notranslate nohighlight">\(\tau\)</span> will always cause our
simple uncorrelated estimate of <span class="math notranslate nohighlight">\(\sigma_m^2\approx \sigma^2/n\)</span> to
be less than the true sample error. The estimate of the error will be
too “good”. On the other hand, the calculation of the full
autocorrelation time poses an efficiency problem if the set of
measurements is very large. The solution to this problem is given by
more practically oriented methods like the blocking technique.</p>
<!-- add ref here to flybjerg --><div class="cell docutils container">
<div class="cell_input docutils container">
<div class="highlight-none notranslate"><div class="highlight"><pre><span></span># Importing various packages
from math import exp, sqrt
from random import random, seed
import numpy as np
import matplotlib.pyplot as plt
# Sample covariance, note the factor 1/(n-1)
def covariance(x, y, n):
sum = 0.0
mean_x = np.mean(x)
mean_y = np.mean(y)
for i in range(0, n):
sum += (x[(i)]-mean_x)*(y[i]-mean_y)
return sum/(n-1.)
n = 100
x = np.random.normal(size=n)
print(np.mean(x))
y = 4+3*x+np.random.normal(size=n)
print(np.mean(y))
z = x**3+np.random.normal(size=n)
print(np.mean(z))
covxx = covariance(x,x,n)
covyy = covariance(y,y,n)
covzz = covariance(z,z,n)
covxy = covariance(x,y,n)
covxz = covariance(x,z,n)
covyz = covariance(y,z,n)
print(covxx,covyy, covzz)
print(covxy,covxz, covyz)
w = np.vstack((x, y, z))
#print(w)
c = np.cov(w)
print(c)
#eigen = np.zeros(n)
Eigvals, Eigvecs = np.linalg.eig(c)
print(Eigvals)
</pre></div>
</div>
</div>
</div>
</section>
<section id="random-numbers">
<h3><span class="section-number">1.1.7. </span>Random Numbers<a class="headerlink" href="#random-numbers" title="Link to this heading">#</a></h3>
<p>Uniform deviates are just random numbers that lie within a specified range
(typically 0 to 1), with any one number in the range just as likely as any other. They
are, in other words, what you probably think random numbers are. However,
we want to distinguish uniform deviates from other sorts of random numbers, for
example numbers drawn from a normal (Gaussian) distribution of specified mean
and standard deviation. These other sorts of deviates are almost always generated by
performing appropriate operations on one or more uniform deviates, as we will see
in subsequent sections. So, a reliable source of random uniform deviates, the subject
of this section, is an essential building block for any sort of stochastic modeling
or Monte Carlo computer work.</p>
<p>A disclaimer is however appropriate. It should be fairly obvious that
something as deterministic as a computer cannot generate purely random numbers.</p>
<p>Numbers generated by any of the standard algorithms are in reality pseudo random
numbers, hopefully abiding to the following criteria:</p>
<ul class="simple">
<li><p>they produce a uniform distribution in the interval [0,1].</p></li>
<li><p>correlations between random numbers are negligible</p></li>
<li><p>the period before the same sequence of random numbers is repeated is as large as possible and finally</p></li>
<li><p>the algorithm should be fast.</p></li>
</ul>
<p>The most common random number generators are based on so-called
Linear congruential relations of the type</p>
<div class="math notranslate nohighlight">
\[
N_i=(aN_{i-1}+c) \mathrm{MOD} (M),
\]</div>
<p>which yield a number in the interval [0,1] through</p>
<div class="math notranslate nohighlight">
\[
x_i=N_i/M
\]</div>
<p>The number
<span class="math notranslate nohighlight">\(M\)</span> is called the period and it should be as large as possible
and
<span class="math notranslate nohighlight">\(N_0\)</span> is the starting value, or seed. The function <span class="math notranslate nohighlight">\(\mathrm{MOD}\)</span> means the remainder,
that is if we were to evaluate <span class="math notranslate nohighlight">\((13)\mathrm{MOD}(9)\)</span>, the outcome is the remainder
of the division <span class="math notranslate nohighlight">\(13/9\)</span>, namely <span class="math notranslate nohighlight">\(4\)</span>.</p>
<p>The problem with such generators is that their outputs are periodic;
they
will start to repeat themselves with a period that is at most <span class="math notranslate nohighlight">\(M\)</span>. If however
the parameters <span class="math notranslate nohighlight">\(a\)</span> and <span class="math notranslate nohighlight">\(c\)</span> are badly chosen, the period may be even shorter.</p>
<p>Consider the following example</p>
<div class="math notranslate nohighlight">
\[
N_i=(6N_{i-1}+7) \mathrm{MOD} (5),
\]</div>
<p>with a seed <span class="math notranslate nohighlight">\(N_0=2\)</span>. This generator produces the sequence
<span class="math notranslate nohighlight">\(4,1,3,0,2,4,1,3,0,2,...\dots\)</span>, i.e., a sequence with period <span class="math notranslate nohighlight">\(5\)</span>.
However, increasing <span class="math notranslate nohighlight">\(M\)</span> may not guarantee a larger period as the following
example shows</p>
<div class="math notranslate nohighlight">
\[
N_i=(27N_{i-1}+11) \mathrm{MOD} (54),
\]</div>
<p>which still, with <span class="math notranslate nohighlight">\(N_0=2\)</span>, results in <span class="math notranslate nohighlight">\(11,38,11,38,11,38,\dots\)</span>, a period of
just <span class="math notranslate nohighlight">\(2\)</span>.</p>
<p>Typical periods for the random generators provided in the program library
are of the order of <span class="math notranslate nohighlight">\(\sim 10^9\)</span> or larger. Other random number generators which have
become increasingly popular are so-called shift-register generators.
In these generators each successive number depends on many preceding
values (rather than the last values as in the linear congruential
generator).
For example, you could make a shift register generator whose <span class="math notranslate nohighlight">\(l\)</span>th
number is the sum of the <span class="math notranslate nohighlight">\(l-i\)</span>th and <span class="math notranslate nohighlight">\(l-j\)</span>th values with modulo <span class="math notranslate nohighlight">\(M\)</span>,</p>
<div class="math notranslate nohighlight">
\[
N_l=(aN_{l-i}+cN_{l-j})\mathrm{MOD}(M).
\]</div>
<p>Such a generator again produces a sequence of pseudorandom numbers
but this time with a period much larger than <span class="math notranslate nohighlight">\(M\)</span>.
It is also possible to construct more elaborate algorithms by including
more than two past terms in the sum of each iteration.
One example is the generator of <a class="reference external" href="http://dl.acm.org/citation.cfm?id=187154">Marsaglia and Zaman</a>
which consists of two congruential relations</p>
<!-- Equation labels as ordinary links -->
<div id="eq:mz1"></div>
<div class="math notranslate nohighlight">
\[
\begin{equation}
N_l=(N_{l-3}-N_{l-1})\mathrm{MOD}(2^{31}-69),
\label{eq:mz1} \tag{12}
\end{equation}
\]</div>
<p>followed by</p>
<!-- Equation labels as ordinary links -->
<div id="eq:mz2"></div>
<div class="math notranslate nohighlight">
\[
\begin{equation}
N_l=(69069N_{l-1}+1013904243)\mathrm{MOD}(2^{32}),
\label{eq:mz2} \tag{13}
\end{equation}
\]</div>
<p>which according to the authors has a period larger than <span class="math notranslate nohighlight">\(2^{94}\)</span>.</p>
<p>Instead of using modular addition, we could use the bitwise
exclusive-OR (<span class="math notranslate nohighlight">\(\oplus\)</span>) operation so that</p>
<div class="math notranslate nohighlight">
\[
N_l=(N_{l-i})\oplus (N_{l-j})
\]</div>
<p>where the bitwise action of <span class="math notranslate nohighlight">\(\oplus\)</span> means that if <span class="math notranslate nohighlight">\(N_{l-i}=N_{l-j}\)</span> the result is
<span class="math notranslate nohighlight">\(0\)</span> whereas if <span class="math notranslate nohighlight">\(N_{l-i}\ne N_{l-j}\)</span> the result is
<span class="math notranslate nohighlight">\(1\)</span>. As an example, consider the case where <span class="math notranslate nohighlight">\(N_{l-i}=6\)</span> and <span class="math notranslate nohighlight">\(N_{l-j}=11\)</span>. The first
one has a bit representation (using 4 bits only) which reads <span class="math notranslate nohighlight">\(0110\)</span> whereas the
second number is <span class="math notranslate nohighlight">\(1011\)</span>. Employing the <span class="math notranslate nohighlight">\(\oplus\)</span> operator yields
<span class="math notranslate nohighlight">\(1101\)</span>, or <span class="math notranslate nohighlight">\(2^3+2^2+2^0=13\)</span>.</p>
<p>In Fortran90, the bitwise <span class="math notranslate nohighlight">\(\oplus\)</span> operation is coded through the intrinsic
function <span class="math notranslate nohighlight">\(\mathrm{IEOR}(m,n)\)</span> where <span class="math notranslate nohighlight">\(m\)</span> and <span class="math notranslate nohighlight">\(n\)</span> are the input numbers, while in <span class="math notranslate nohighlight">\(C\)</span>
it is given by <span class="math notranslate nohighlight">\(m\wedge n\)</span>.</p>
<p>We show here how the linear congruential algorithm can be implemented, namely</p>
<div class="math notranslate nohighlight">
\[
N_i=(aN_{i-1}) \mathrm{MOD} (M).
\]</div>
<p>However, since <span class="math notranslate nohighlight">\(a\)</span> and <span class="math notranslate nohighlight">\(N_{i-1}\)</span> are integers and their multiplication
could become greater than the standard 32 bit integer, there is a trick via
Schrage’s algorithm which approximates the multiplication
of large integers through the factorization</p>
<div class="math notranslate nohighlight">
\[
M=aq+r,
\]</div>
<p>where we have defined</p>
<div class="math notranslate nohighlight">
\[
q=[M/a],
\]</div>
<p>and</p>
<div class="math notranslate nohighlight">
\[
r = M\hspace{0.1cm}\mathrm{MOD} \hspace{0.1cm}a.
\]</div>
<p>where the brackets denote integer division. In the code below the numbers
<span class="math notranslate nohighlight">\(q\)</span> and <span class="math notranslate nohighlight">\(r\)</span> are chosen so that <span class="math notranslate nohighlight">\(r &lt; q\)</span>.</p>
<p>To see how this works we note first that</p>
<!-- Equation labels as ordinary links -->
<div id="eq:rntrick1"></div>
<div class="math notranslate nohighlight">
\[
\begin{equation}
(aN_{i-1}) \mathrm{MOD} (M)= (aN_{i-1}-[N_{i-1}/q]M)\mathrm{MOD} (M),
\label{eq:rntrick1} \tag{14}
\end{equation}
\]</div>
<p>since we can add or subtract any integer multiple of <span class="math notranslate nohighlight">\(M\)</span> from <span class="math notranslate nohighlight">\(aN_{i-1}\)</span>.
The last term <span class="math notranslate nohighlight">\([N_{i-1}/q]M\mathrm{MOD}(M)\)</span> is zero since the integer division
<span class="math notranslate nohighlight">\([N_{i-1}/q]\)</span> just yields a constant which is multiplied with <span class="math notranslate nohighlight">\(M\)</span>.</p>
<p>We can now rewrite Eq. (<a class="reference internal" href="#eq:rntrick1"><span class="xref myst">14</span></a>) as</p>
<!-- Equation labels as ordinary links -->
<div id="eq:rntrick2"></div>
<div class="math notranslate nohighlight">
\[
\begin{equation}
(aN_{i-1}) \mathrm{MOD} (M)= (aN_{i-1}-[N_{i-1}/q](aq+r))\mathrm{MOD} (M),
\label{eq:rntrick2} \tag{15}
\end{equation}
\]</div>
<p>which results</p>
<!-- Equation labels as ordinary links -->
<div id="eq:rntrick3"></div>
<div class="math notranslate nohighlight">
\[
\begin{equation}
(aN_{i-1}) \mathrm{MOD} (M)= \left(a(N_{i-1}-[N_{i-1}/q]q)-[N_{i-1}/q]r)\right)\mathrm{MOD} (M),
\label{eq:rntrick3} \tag{16}
\end{equation}
\]</div>
<p>yielding</p>
<!-- Equation labels as ordinary links -->
<div id="eq:rntrick4"></div>
<div class="math notranslate nohighlight">
\[
\begin{equation}
(aN_{i-1}) \mathrm{MOD} (M)= \left(a(N_{i-1}\mathrm{MOD} (q)) -[N_{i-1}/q]r)\right)\mathrm{MOD} (M).
\label{eq:rntrick4} \tag{17}
\end{equation}
\]</div>
<p>The term <span class="math notranslate nohighlight">\([N_{i-1}/q]r\)</span> is always smaller or equal <span class="math notranslate nohighlight">\(N_{i-1}(r/q)\)</span> and with <span class="math notranslate nohighlight">\(r &lt; q\)</span> we obtain always a
number smaller than <span class="math notranslate nohighlight">\(N_{i-1}\)</span>, which is smaller than <span class="math notranslate nohighlight">\(M\)</span>.
And since the number <span class="math notranslate nohighlight">\(N_{i-1}\mathrm{MOD} (q)\)</span> is between zero and <span class="math notranslate nohighlight">\(q-1\)</span> then
<span class="math notranslate nohighlight">\(a(N_{i-1}\mathrm{MOD} (q))&lt; aq\)</span>. Combined with our definition of <span class="math notranslate nohighlight">\(q=[M/a]\)</span> ensures that
this term is also smaller than <span class="math notranslate nohighlight">\(M\)</span> meaning that both terms fit into a
32-bit signed integer. None of these two terms can be negative, but their difference could.
The algorithm below adds <span class="math notranslate nohighlight">\(M\)</span> if their difference is negative.
Note that the program uses the bitwise <span class="math notranslate nohighlight">\(\oplus\)</span> operator to generate
the starting point for each generation of a random number. The period
of <span class="math notranslate nohighlight">\(ran0\)</span> is <span class="math notranslate nohighlight">\(\sim 2.1\times 10^{9}\)</span>. A special feature of this
algorithm is that is should never be called with the initial seed
set to <span class="math notranslate nohighlight">\(0\)</span>.</p>
<p>As mentioned previously, the underlying PDF for the generation of
random numbers is the uniform distribution, meaning that the
probability for finding a number <span class="math notranslate nohighlight">\(x\)</span> in the interval [0,1] is <span class="math notranslate nohighlight">\(p(x)=1\)</span>.</p>
<p>A random number generator should produce numbers which are uniformly distributed
in this interval. The table shows the distribution of <span class="math notranslate nohighlight">\(N=10000\)</span> random
numbers generated by the functions in the program library.
We note in this table that the number of points in the various
intervals <span class="math notranslate nohighlight">\(0.0-0.1\)</span>, <span class="math notranslate nohighlight">\(0.1-0.2\)</span> etc are fairly close to <span class="math notranslate nohighlight">\(1000\)</span>, with some minor
deviations.</p>
<p>Two additional measures are the standard deviation <span class="math notranslate nohighlight">\(\sigma\)</span> and the mean
<span class="math notranslate nohighlight">\(\mu=\langle x\rangle\)</span>.</p>
<p>For the uniform distribution, the mean value <span class="math notranslate nohighlight">\(\mu\)</span> is then</p>
<div class="math notranslate nohighlight">
\[
\mu=\langle x\rangle=\frac{1}{2}
\]</div>
<p>while the standard deviation is</p>
<div class="math notranslate nohighlight">
\[
\sigma=\sqrt{\langle x^2\rangle-\mu^2}=\frac{1}{\sqrt{12}}=0.2886.
\]</div>
<p>The various random number generators produce results which agree rather well with
these limiting values.</p>
<table class="dotable" border="1">
<thead>
<tr><th align="center">$x$-bin </th> <th align="center"> ran0 </th> <th align="center"> ran1 </th> <th align="center"> ran2 </th> <th align="center"> ran3 </th> </tr>
</thead>
<tbody>
<tr><td align="center"> 0.0-0.1 </td> <td align="right"> 1013 </td> <td align="right"> 991 </td> <td align="right"> 938 </td> <td align="right"> 1047 </td> </tr>
<tr><td align="center"> 0.1-0.2 </td> <td align="right"> 1002 </td> <td align="right"> 1009 </td> <td align="right"> 1040 </td> <td align="right"> 1030 </td> </tr>
<tr><td align="center"> 0.2-0.3 </td> <td align="right"> 989 </td> <td align="right"> 999 </td> <td align="right"> 1030 </td> <td align="right"> 993 </td> </tr>
<tr><td align="center"> 0.3-0.4 </td> <td align="right"> 939 </td> <td align="right"> 960 </td> <td align="right"> 1023 </td> <td align="right"> 937 </td> </tr>
<tr><td align="center"> 0.4-0.5 </td> <td align="right"> 1038 </td> <td align="right"> 1001 </td> <td align="right"> 1002 </td> <td align="right"> 992 </td> </tr>
<tr><td align="center"> 0.5-0.6 </td> <td align="right"> 1037 </td> <td align="right"> 1047 </td> <td align="right"> 1009 </td> <td align="right"> 1009 </td> </tr>
<tr><td align="center"> 0.6-0.7 </td> <td align="right"> 1005 </td> <td align="right"> 989 </td> <td align="right"> 1003 </td> <td align="right"> 989 </td> </tr>
<tr><td align="center"> 0.7-0.8 </td> <td align="right"> 986 </td> <td align="right"> 962 </td> <td align="right"> 985 </td> <td align="right"> 954 </td> </tr>
<tr><td align="center"> 0.8-0.9 </td> <td align="right"> 1000 </td> <td align="right"> 1027 </td> <td align="right"> 1009 </td> <td align="right"> 1023 </td> </tr>
<tr><td align="center"> 0.9-1.0 </td> <td align="right"> 991 </td> <td align="right"> 1015 </td> <td align="right"> 961 </td> <td align="right"> 1026 </td> </tr>
<tr><td align="center"> $\mu$ </td> <td align="right"> 0.4997 </td> <td align="right"> 0.5018 </td> <td align="right"> 0.4992 </td> <td align="right"> 0.4990 </td> </tr>
<tr><td align="center"> $\sigma$ </td> <td align="right"> 0.2882 </td> <td align="right"> 0.2892 </td> <td align="right"> 0.2861 </td> <td align="right"> 0.2915 </td> </tr>
</tbody>
</table>
<p>The following simple Python code plots the distribution of the produced random numbers using the linear congruential RNG employed by Python. The trend displayed in the previous table is seen rather clearly.</p>
<div class="cell docutils container">
<div class="cell_input docutils container">
<div class="highlight-none notranslate"><div class="highlight"><pre><span></span>#!/usr/bin/env python
import numpy as np
import matplotlib.mlab as mlab
import matplotlib.pyplot as plt
import random
# initialize the rng with a seed
random.seed()
counts = 10000
values = np.zeros(counts)
for i in range (1, counts, 1):
values[i] = random.random()
# the histogram of the data
n, bins, patches = plt.hist(values, 10, facecolor=&#39;green&#39;)
plt.xlabel(&#39;$x$&#39;)
plt.ylabel(&#39;Number of counts&#39;)
plt.title(r&#39;Test of uniform distribution&#39;)
plt.axis([0, 1, 0, 1100])
plt.grid(True)
plt.show()
</pre></div>
</div>
</div>
</div>
<p>Since our random numbers, which are typically generated via a linear congruential algorithm,
are never fully independent, we can then define
an important test which measures the degree of correlation, namely the so-called<br />
auto-correlation function defined previously, see again Eq. (<a class="reference internal" href="#eq:autocorrelformal"><span class="xref myst">9</span></a>).
We rewrite it here as</p>
<div class="math notranslate nohighlight">
\[
C_k=\frac{f_d}
{\sigma^2},
\]</div>
<p>with <span class="math notranslate nohighlight">\(C_0=1\)</span>. Recall that
<span class="math notranslate nohighlight">\(\sigma^2=\langle x_i^2\rangle-\langle x_i\rangle^2\)</span> and that</p>
<div class="math notranslate nohighlight">
\[
f_d = \frac{1}{nm}\sum_{\alpha=1}^m\sum_{k=1}^{n-d}(x_{\alpha,k}-\langle X_m \rangle)(x_{\alpha,k+d}-\langle X_m \rangle),
\]</div>
<p>The non-vanishing of <span class="math notranslate nohighlight">\(C_k\)</span> for <span class="math notranslate nohighlight">\(k\ne 0\)</span> means that the random
numbers are not independent. The independence of the random numbers is crucial
in the evaluation of other expectation values. If they are not independent, our
assumption for approximating <span class="math notranslate nohighlight">\(\sigma_N\)</span> is no longer valid.</p>
</section>
<section id="autocorrelation-function">
<h3><span class="section-number">1.1.8. </span>Autocorrelation function<a class="headerlink" href="#autocorrelation-function" title="Link to this heading">#</a></h3>
<p>This program computes the autocorrelation function as discussed in the equation on the previous slide for random numbers generated with the normal distribution <span class="math notranslate nohighlight">\(N(0,1)\)</span>.</p>
<div class="cell docutils container">
<div class="cell_input docutils container">
<div class="highlight-none notranslate"><div class="highlight"><pre><span></span># Importing various packages
from math import exp, sqrt
from random import random, seed
import numpy as np
import matplotlib.pyplot as plt
def autocovariance(x, n, k, mean_x):
sum = 0.0
for i in range(0, n-k):
sum += (x[(i+k)]-mean_x)*(x[i]-mean_x)
return sum/n
n = 1000
x=np.random.normal(size=n)
autocor = np.zeros(n)
figaxis = np.zeros(n)
mean_x=np.mean(x)
var_x = np.var(x)
print(mean_x, var_x)
for i in range (0, n):
figaxis[i] = i
autocor[i]=(autocovariance(x, n, i, mean_x))/var_x
plt.plot(figaxis, autocor, &quot;r-&quot;)
plt.axis([0,n,-0.1, 1.0])
plt.xlabel(r&#39;$i$&#39;)
plt.ylabel(r&#39;$\gamma_i$&#39;)
plt.title(r&#39;Autocorrelation function&#39;)
plt.show()
</pre></div>
</div>
</div>
</div>
<p>As can be seen from the plot, the first point gives back the variance and a value of one.
For the remaining values we notice that there are still non-zero values for the auto-correlation function.</p>
</section>
</section>
</section>
<script type="text/x-thebe-config">
{
requestKernel: true,
binderOptions: {
repo: "binder-examples/jupyter-stacks-datascience",
ref: "master",
},
codeMirrorConfig: {
theme: "abcdef",
mode: "python"
},
kernelOptions: {
name: "python3",
path: "./."
},
predefinedOutput: true
}
</script>
<script>kernelName = 'python3'</script>
</article>
<footer class="prev-next-footer d-print-none">
<div class="prev-next-area">
<a class="left-prev"
href="textbooks.html"
title="previous page">
<i class="fa-solid fa-angle-left"></i>
<div class="prev-next-info">
<p class="prev-next-subtitle">previous</p>
<p class="prev-next-title">Textbooks</p>
</div>
</a>
<a class="right-next"
href="linalg.html"
title="next page">
<div class="prev-next-info">
<p class="prev-next-subtitle">next</p>
<p class="prev-next-title"><span class="section-number">2. </span>Linear Algebra, Handling of Arrays and more Python Features</p>
</div>
<i class="fa-solid fa-angle-right"></i>
</a>
</div>
</footer>
</div>
<div class="bd-sidebar-secondary bd-toc"><div class="sidebar-secondary-items sidebar-secondary__inner">
<div class="sidebar-secondary-item">
<div class="page-toc tocsection onthispage">
<i class="fa-solid fa-list"></i> Contents
</div>
<nav class="bd-toc-nav page-toc">
<ul class="visible nav section-nav flex-column">
<li class="toc-h2 nav-item toc-entry"><a class="reference internal nav-link" href="#domains-and-probabilities">1.1. Domains and probabilities</a><ul class="nav section-nav flex-column">
<li class="toc-h3 nav-item toc-entry"><a class="reference internal nav-link" href="#stochastic-variables-and-the-main-concepts-the-discrete-case">1.1.1. Stochastic variables and the main concepts, the discrete case</a></li>
<li class="toc-h3 nav-item toc-entry"><a class="reference internal nav-link" href="#properties-of-pdfs">1.1.2. Properties of PDFs</a></li>
<li class="toc-h3 nav-item toc-entry"><a class="reference internal nav-link" href="#expectation-values">1.1.3. Expectation values</a></li>
<li class="toc-h3 nav-item toc-entry"><a class="reference internal nav-link" href="#probability-distribution-functions">1.1.4. Probability Distribution Functions</a></li>
<li class="toc-h3 nav-item toc-entry"><a class="reference internal nav-link" href="#meet-the-covariance">1.1.5. Meet the covariance!</a></li>
<li class="toc-h3 nav-item toc-entry"><a class="reference internal nav-link" href="#numerical-experiments-and-the-covariance-central-limit-theorem">1.1.6. Numerical experiments and the covariance, central limit theorem</a></li>
<li class="toc-h3 nav-item toc-entry"><a class="reference internal nav-link" href="#random-numbers">1.1.7. Random Numbers</a></li>
<li class="toc-h3 nav-item toc-entry"><a class="reference internal nav-link" href="#autocorrelation-function">1.1.8. Autocorrelation function</a></li>
</ul>
</li>
</ul>
</nav></div>
</div></div>
</div>
<footer class="bd-footer-content">
<div class="bd-footer-content__inner container">
<div class="footer-item">
<p class="component-author">
By Morten Hjorth-Jensen
</p>
</div>
<div class="footer-item">
<p class="copyright">
© Copyright 2023.
<br/>
</p>
</div>
<div class="footer-item">
</div>
<div class="footer-item">
</div>
</div>
</footer>
</main>
</div>
</div>
<!-- Scripts loaded after <body> so the DOM is not blocked -->
<script src="_static/scripts/bootstrap.js?digest=dfe6caa3a7d634c4db9b"></script>
<script src="_static/scripts/pydata-sphinx-theme.js?digest=dfe6caa3a7d634c4db9b"></script>
<footer class="bd-footer">
</footer>
</body>
</html>
Reference in New Issue View Git Blame Copy Permalink