Morten Hjorth-Jensen
2423 lines
153 KiB
HTML
2423 lines
153 KiB
HTML
|
||
<!DOCTYPE html>
|
||
|
||
<html>
|
||
<head>
|
||
<meta charset="utf-8" />
|
||
<meta name="viewport" content="width=device-width, initial-scale=1.0" />
|
||
<title>1. Elements of Probability Theory and Statistical Data Analysis — Applied Data Analysis and Machine Learning</title>
|
||
|
||
<link href="_static/css/theme.css" rel="stylesheet" />
|
||
<link href="_static/css/index.c5995385ac14fb8791e8eb36b4908be2.css" rel="stylesheet" />
|
||
|
||
|
||
<link rel="stylesheet"
|
||
href="_static/vendor/fontawesome/5.13.0/css/all.min.css">
|
||
<link rel="preload" as="font" type="font/woff2" crossorigin
|
||
href="_static/vendor/fontawesome/5.13.0/webfonts/fa-solid-900.woff2">
|
||
<link rel="preload" as="font" type="font/woff2" crossorigin
|
||
href="_static/vendor/fontawesome/5.13.0/webfonts/fa-brands-400.woff2">
|
||
|
||
|
||
|
||
|
||
|
||
<link rel="stylesheet" href="_static/pygments.css" type="text/css" />
|
||
<link rel="stylesheet" href="_static/sphinx-book-theme.5f77b4aec8189eecf79907ce328c390d.css" type="text/css" />
|
||
<link rel="stylesheet" type="text/css" href="_static/togglebutton.css" />
|
||
<link rel="stylesheet" type="text/css" href="_static/copybutton.css" />
|
||
<link rel="stylesheet" type="text/css" href="_static/mystnb.css" />
|
||
<link rel="stylesheet" type="text/css" href="_static/sphinx-thebe.css" />
|
||
<link rel="stylesheet" type="text/css" href="_static/panels-main.c949a650a448cc0ae9fd3441c0e17fb0.css" />
|
||
<link rel="stylesheet" type="text/css" href="_static/panels-variables.06eb56fa6e07937060861dad626602ad.css" />
|
||
|
||
<link rel="preload" as="script" href="_static/js/index.1c5a1a01449ed65a7b51.js">
|
||
|
||
<script id="documentation_options" data-url_root="./" src="_static/documentation_options.js"></script>
|
||
<script src="_static/jquery.js"></script>
|
||
<script src="_static/underscore.js"></script>
|
||
<script src="_static/doctools.js"></script>
|
||
<script src="_static/language_data.js"></script>
|
||
<script src="_static/togglebutton.js"></script>
|
||
<script src="_static/clipboard.min.js"></script>
|
||
<script src="_static/copybutton.js"></script>
|
||
<script >var togglebuttonSelector = '.toggle, .admonition.dropdown, .tag_hide_input div.cell_input, .tag_hide-input div.cell_input, .tag_hide_output div.cell_output, .tag_hide-output div.cell_output, .tag_hide_cell.cell, .tag_hide-cell.cell';</script>
|
||
<script src="_static/sphinx-book-theme.12a9622fbb08dcb3a2a40b2c02b83a57.js"></script>
|
||
<script async="async" src="https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.7/latest.js?config=TeX-AMS-MML_HTMLorMML"></script>
|
||
<script type="text/x-mathjax-config">MathJax.Hub.Config({"tex2jax": {"inlineMath": [["\\(", "\\)"]], "displayMath": [["\\[", "\\]"]], "processRefs": false, "processEnvironments": false}})</script>
|
||
<script async="async" src="https://unpkg.com/thebelab@latest/lib/index.js"></script>
|
||
<script >
|
||
const thebe_selector = ".thebe"
|
||
const thebe_selector_input = "pre"
|
||
const thebe_selector_output = ".output"
|
||
</script>
|
||
<script async="async" src="_static/sphinx-thebe.js"></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-spy="scroll" data-target="#bd-toc-nav" data-offset="80">
|
||
|
||
<div class="container-fluid" id="banner"></div>
|
||
|
||
|
||
|
||
<div class="container-xl">
|
||
<div class="row">
|
||
|
||
<div class="col-12 col-md-3 bd-sidebar site-navigation show" id="site-navigation">
|
||
|
||
<div class="navbar-brand-box">
|
||
<a class="navbar-brand text-wrap" href="index.html">
|
||
|
||
<img src="_static/logo.png" class="logo" alt="logo">
|
||
|
||
|
||
<h1 class="site-logo" id="site-title">Applied Data Analysis and Machine Learning</h1>
|
||
|
||
</a>
|
||
</div><form class="bd-search d-flex align-items-center" action="search.html" method="get">
|
||
<i class="icon fas fa-search"></i>
|
||
<input type="search" class="form-control" name="q" id="search-input" placeholder="Search this book..." aria-label="Search this book..." autocomplete="off" >
|
||
</form><nav class="bd-links" id="bd-docs-nav" aria-label="Main navigation">
|
||
<div class="bd-toc-item active">
|
||
<ul class="nav bd-sidenav">
|
||
<li class="toctree-l1">
|
||
<a class="reference internal" href="intro.html">
|
||
Applied Data Analysis and Machine Learning, FYS-STK3155/4155 at the University of Oslo, Norway
|
||
</a>
|
||
</li>
|
||
</ul>
|
||
<p class="caption">
|
||
<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">
|
||
Teaching schedule with links to material
|
||
</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 class="caption">
|
||
<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 class="caption">
|
||
<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 class="caption">
|
||
<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 class="caption">
|
||
<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 Analysis
|
||
</a>
|
||
</li>
|
||
</ul>
|
||
<p class="caption">
|
||
<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>
|
||
</ul>
|
||
|
||
</div>
|
||
</nav> <!-- To handle the deprecated key -->
|
||
|
||
<div class="navbar_extra_footer">
|
||
Powered by <a href="https://jupyterbook.org">Jupyter Book</a>
|
||
</div>
|
||
|
||
</div>
|
||
|
||
|
||
|
||
|
||
|
||
|
||
<main class="col py-md-3 pl-md-4 bd-content overflow-auto" role="main">
|
||
|
||
<div class="topbar container-xl fixed-top">
|
||
<div class="topbar-contents row">
|
||
<div class="col-12 col-md-3 bd-topbar-whitespace site-navigation show"></div>
|
||
<div class="col pl-md-4 topbar-main">
|
||
|
||
<button id="navbar-toggler" class="navbar-toggler ml-0" type="button" data-toggle="collapse"
|
||
data-toggle="tooltip" data-placement="bottom" data-target=".site-navigation" aria-controls="navbar-menu"
|
||
aria-expanded="true" aria-label="Toggle navigation" aria-controls="site-navigation"
|
||
title="Toggle navigation" data-toggle="tooltip" data-placement="left">
|
||
<i class="fas fa-bars"></i>
|
||
<i class="fas fa-arrow-left"></i>
|
||
<i class="fas fa-arrow-up"></i>
|
||
</button>
|
||
|
||
|
||
<div class="dropdown-buttons-trigger">
|
||
<button id="dropdown-buttons-trigger" class="btn btn-secondary topbarbtn" aria-label="Download this page"><i
|
||
class="fas fa-download"></i></button>
|
||
|
||
<div class="dropdown-buttons">
|
||
<!-- ipynb file if we had a myst markdown file -->
|
||
|
||
<!-- Download raw file -->
|
||
<a class="dropdown-buttons" href="_sources/statistics.ipynb"><button type="button"
|
||
class="btn btn-secondary topbarbtn" title="Download source file" data-toggle="tooltip"
|
||
data-placement="left">.ipynb</button></a>
|
||
<!-- Download PDF via print -->
|
||
<button type="button" id="download-print" class="btn btn-secondary topbarbtn" title="Print to PDF"
|
||
onClick="window.print()" data-toggle="tooltip" data-placement="left">.pdf</button>
|
||
</div>
|
||
</div>
|
||
|
||
<!-- Source interaction buttons -->
|
||
|
||
<!-- Full screen (wrap in <a> to have style consistency -->
|
||
|
||
<a class="full-screen-button"><button type="button" class="btn btn-secondary topbarbtn" data-toggle="tooltip"
|
||
data-placement="bottom" onclick="toggleFullScreen()" aria-label="Fullscreen mode"
|
||
title="Fullscreen mode"><i
|
||
class="fas fa-expand"></i></button></a>
|
||
|
||
<!-- Launch buttons -->
|
||
|
||
</div>
|
||
|
||
<!-- Table of contents -->
|
||
<div class="d-none d-md-block col-md-2 bd-toc show">
|
||
|
||
<div class="tocsection onthispage pt-5 pb-3">
|
||
<i class="fas fa-list"></i> Contents
|
||
</div>
|
||
<nav id="bd-toc-nav">
|
||
<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>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#tossing-the-dice">
|
||
1.2. Tossing the dice
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#stochastic-variables">
|
||
1.3. Stochastic variables
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#stochastic-variables-and-the-main-concepts-the-discrete-case">
|
||
1.4. Stochastic variables and the main concepts, the discrete case
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#stochastic-variables-and-the-main-concepts-the-continuous-case">
|
||
1.5. Stochastic variables and the main concepts, the continuous case
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#the-cumulative-probability">
|
||
1.6. The cumulative probability
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#properties-of-pdfs">
|
||
1.7. Properties of PDFs
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#important-distributions-the-uniform-distribution">
|
||
1.8. Important distributions, the uniform distribution
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#gaussian-distribution">
|
||
1.9. Gaussian distribution
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#exponential-distribution">
|
||
1.10. Exponential distribution
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#expectation-values">
|
||
1.11. Expectation values
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#stochastic-variables-and-the-main-concepts-mean-values">
|
||
1.12. Stochastic variables and the main concepts, mean values
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#stochastic-variables-and-the-main-concepts-central-moments-the-variance">
|
||
1.13. Stochastic variables and the main concepts, central moments, the variance
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#probability-distribution-functions">
|
||
1.14. Probability Distribution Functions
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#id1">
|
||
1.15. Probability Distribution Functions
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#the-three-famous-probability-distribution-functions">
|
||
1.16. The three famous Probability Distribution Functions
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#probability-distribution-functions-the-normal-distribution">
|
||
1.17. Probability Distribution Functions, the normal distribution
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#id2">
|
||
1.18. Probability Distribution Functions, the normal distribution
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#probability-distribution-functions-the-cumulative-distribution">
|
||
1.19. Probability Distribution Functions, the cumulative distribution
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#probability-distribution-functions-other-important-distribution">
|
||
1.20. Probability Distribution Functions, other important distribution
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#probability-distribution-functions-the-binomial-distribution">
|
||
1.21. Probability Distribution Functions, the binomial distribution
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#probability-distribution-functions-poisson-s-distribution">
|
||
1.22. Probability Distribution Functions, Poisson’s distribution
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#id3">
|
||
1.23. Probability Distribution Functions, Poisson’s distribution
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#meet-the-covariance">
|
||
1.24. Meet the covariance!
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#meet-the-covariance-in-matrix-disguise">
|
||
1.25. Meet the covariance in matrix disguise
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#covariance">
|
||
1.26. Covariance
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#meet-the-covariance-uncorrelated-events">
|
||
1.27. Meet the covariance, uncorrelated events
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#numerical-experiments-and-the-covariance">
|
||
1.28. Numerical experiments and the covariance
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#id4">
|
||
1.29. Numerical experiments and the covariance
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#numerical-experiments-and-the-covariance-actual-situations">
|
||
1.30. Numerical experiments and the covariance, actual situations
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#numerical-experiments-and-the-covariance-our-observables">
|
||
1.31. Numerical experiments and the covariance, our observables
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#numerical-experiments-and-the-covariance-the-sample-variance">
|
||
1.32. Numerical experiments and the covariance, the sample variance
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#numerical-experiments-and-the-covariance-central-limit-theorem">
|
||
1.33. Numerical experiments and the covariance, central limit theorem
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#definition-of-correlation-functions-and-standard-deviation">
|
||
1.34. Definition of Correlation Functions and Standard Deviation
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#id5">
|
||
1.35. Definition of Correlation Functions and Standard Deviation
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#id6">
|
||
1.36. Definition of Correlation Functions and Standard Deviation
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#id7">
|
||
1.37. Definition of Correlation Functions and Standard Deviation
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#definition-of-correlation-functions-and-standard-deviation-sample-variance">
|
||
1.38. Definition of Correlation Functions and Standard Deviation, sample variance
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#id8">
|
||
1.39. Definition of Correlation Functions and Standard Deviation
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#code-to-compute-the-covariance-matrix-and-the-covariance">
|
||
1.40. Code to compute the Covariance matrix and the Covariance
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#random-numbers">
|
||
1.41. Random Numbers
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#random-numbers-better-name-pseudo-random-numbers">
|
||
1.42. Random Numbers, better name: pseudo random numbers
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#random-number-generator-rng">
|
||
1.43. Random number generator RNG
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#random-number-generator-rng-and-periodic-outputs">
|
||
1.44. Random number generator RNG and periodic outputs
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#random-number-generator-rng-and-its-period">
|
||
1.45. Random number generator RNG and its period
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#random-number-generator-rng-other-examples">
|
||
1.46. Random number generator RNG, other examples
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#id9">
|
||
1.47. Random number generator RNG, other examples
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#random-number-generator-rng-ran0">
|
||
1.48. Random number generator RNG, RAN0
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#id10">
|
||
1.49. Random number generator RNG, RAN0
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#id11">
|
||
1.50. Random number generator RNG, RAN0
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#id12">
|
||
1.51. Random number generator RNG, RAN0
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#random-number-generator-rng-ran0-code">
|
||
1.52. Random number generator RNG, RAN0 code
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#properties-of-selected-random-number-generators">
|
||
1.53. Properties of Selected Random Number Generators
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#id13">
|
||
1.54. Properties of Selected Random Number Generators
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#id14">
|
||
1.55. Properties of Selected Random Number Generators
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#simple-demonstration-of-rngs-using-python">
|
||
1.56. Simple demonstration of RNGs using python
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#id15">
|
||
1.57. Properties of Selected Random Number Generators
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#autocorrelation-function">
|
||
1.58. Autocorrelation function
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#correlation-function-and-which-random-number-generators-should-i-use">
|
||
1.59. Correlation function and which random number generators should I use
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#which-rng-should-i-use">
|
||
1.60. Which RNG should I use?
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#how-to-use-the-mersenne-generator">
|
||
1.61. How to use the Mersenne generator
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#why-blocking">
|
||
1.62. Why blocking?
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#id16">
|
||
1.63. Why blocking?
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#code-to-demonstrate-the-calculation-of-the-autocorrelation-function">
|
||
1.64. Code to demonstrate the calculation of the autocorrelation function
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#what-is-blocking">
|
||
1.65. What is blocking?
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#id17">
|
||
1.66. What is blocking?
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#id18">
|
||
1.67. What is blocking?
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#implementation">
|
||
1.68. Implementation
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#actual-implementation-with-code-main-function">
|
||
1.69. Actual implementation with code, main function
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#the-bootstrap-method">
|
||
1.70. The Bootstrap method
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#bootstrapping">
|
||
1.71. Bootstrapping
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#bootstrapping-recipe">
|
||
1.72. Bootstrapping, recipe
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#bootstrapping-code">
|
||
1.73. Bootstrapping, code
|
||
</a>
|
||
</li>
|
||
<li class="toc-h2 nav-item toc-entry">
|
||
<a class="reference internal nav-link" href="#jackknife-code">
|
||
1.74. Jackknife, code
|
||
</a>
|
||
</li>
|
||
</ul>
|
||
|
||
</nav>
|
||
</div>
|
||
</div>
|
||
</div>
|
||
<div id="main-content" class="row">
|
||
<div class="col-12 col-md-9 pl-md-3 pr-md-0">
|
||
|
||
<div>
|
||
|
||
<div class="section" 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="Permalink to this headline">¶</a></h1>
|
||
<div class="section" id="domains-and-probabilities">
|
||
<h2><span class="section-number">1.1. </span>Domains and probabilities<a class="headerlink" href="#domains-and-probabilities" title="Permalink to this headline">¶</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>
|
||
</div>
|
||
<div class="section" id="tossing-the-dice">
|
||
<h2><span class="section-number">1.2. </span>Tossing the dice<a class="headerlink" href="#tossing-the-dice" title="Permalink to this headline">¶</a></h2>
|
||
<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>
|
||
</div>
|
||
<div class="section" id="stochastic-variables">
|
||
<h2><span class="section-number">1.3. </span>Stochastic variables<a class="headerlink" href="#stochastic-variables" title="Permalink to this headline">¶</a></h2>
|
||
<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>
|
||
</div>
|
||
<div class="section" id="stochastic-variables-and-the-main-concepts-the-discrete-case">
|
||
<h2><span class="section-number">1.4. </span>Stochastic variables and the main concepts, the discrete case<a class="headerlink" href="#stochastic-variables-and-the-main-concepts-the-discrete-case" title="Permalink to this headline">¶</a></h2>
|
||
<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>
|
||
</div>
|
||
<div class="section" id="stochastic-variables-and-the-main-concepts-the-continuous-case">
|
||
<h2><span class="section-number">1.5. </span>Stochastic variables and the main concepts, the continuous case<a class="headerlink" href="#stochastic-variables-and-the-main-concepts-the-continuous-case" title="Permalink to this headline">¶</a></h2>
|
||
<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>
|
||
</div>
|
||
<div class="section" id="the-cumulative-probability">
|
||
<h2><span class="section-number">1.6. </span>The cumulative probability<a class="headerlink" href="#the-cumulative-probability" title="Permalink to this headline">¶</a></h2>
|
||
<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>
|
||
</div>
|
||
<div class="section" id="properties-of-pdfs">
|
||
<h2><span class="section-number">1.7. </span>Properties of PDFs<a class="headerlink" href="#properties-of-pdfs" title="Permalink to this headline">¶</a></h2>
|
||
<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) & = 1,\\
|
||
\int_{x\in\mathbb D} p(x)\,dx & = 1.
|
||
\end{align*}
|
||
\end{split}\]</div>
|
||
</div>
|
||
<div class="section" id="important-distributions-the-uniform-distribution">
|
||
<h2><span class="section-number">1.8. </span>Important distributions, the uniform distribution<a class="headerlink" href="#important-distributions-the-uniform-distribution" title="Permalink to this headline">¶</a></h2>
|
||
<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 & \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>
|
||
</div>
|
||
<div class="section" id="gaussian-distribution">
|
||
<h2><span class="section-number">1.9. </span>Gaussian distribution<a class="headerlink" href="#gaussian-distribution" title="Permalink to this headline">¶</a></h2>
|
||
<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-ipython3 notranslate"><div class="highlight"><pre><span></span><span class="o">%</span><span class="k">matplotlib</span> inline
|
||
|
||
<span class="kn">import</span> <span class="nn">numpy</span> <span class="k">as</span> <span class="nn">np</span>
|
||
<span class="kn">from</span> <span class="nn">math</span> <span class="kn">import</span> <span class="n">acos</span><span class="p">,</span> <span class="n">exp</span><span class="p">,</span> <span class="n">sqrt</span>
|
||
<span class="kn">from</span> <span class="nn">matplotlib</span> <span class="kn">import</span> <span class="n">pyplot</span> <span class="k">as</span> <span class="n">plt</span>
|
||
<span class="kn">from</span> <span class="nn">matplotlib</span> <span class="kn">import</span> <span class="n">rc</span><span class="p">,</span> <span class="n">rcParams</span>
|
||
<span class="kn">import</span> <span class="nn">matplotlib.units</span> <span class="k">as</span> <span class="nn">units</span>
|
||
<span class="kn">import</span> <span class="nn">matplotlib.ticker</span> <span class="k">as</span> <span class="nn">ticker</span>
|
||
<span class="n">rc</span><span class="p">(</span><span class="s1">'text'</span><span class="p">,</span><span class="n">usetex</span><span class="o">=</span><span class="kc">True</span><span class="p">)</span>
|
||
<span class="n">rc</span><span class="p">(</span><span class="s1">'font'</span><span class="p">,</span><span class="o">**</span><span class="p">{</span><span class="s1">'family'</span><span class="p">:</span><span class="s1">'serif'</span><span class="p">,</span><span class="s1">'serif'</span><span class="p">:[</span><span class="s1">'Gaussian distribution'</span><span class="p">]})</span>
|
||
<span class="n">font</span> <span class="o">=</span> <span class="p">{</span><span class="s1">'family'</span> <span class="p">:</span> <span class="s1">'serif'</span><span class="p">,</span>
|
||
<span class="s1">'color'</span> <span class="p">:</span> <span class="s1">'darkred'</span><span class="p">,</span>
|
||
<span class="s1">'weight'</span> <span class="p">:</span> <span class="s1">'normal'</span><span class="p">,</span>
|
||
<span class="s1">'size'</span> <span class="p">:</span> <span class="mi">16</span><span class="p">,</span>
|
||
<span class="p">}</span>
|
||
<span class="n">pi</span> <span class="o">=</span> <span class="n">acos</span><span class="p">(</span><span class="o">-</span><span class="mf">1.0</span><span class="p">)</span>
|
||
<span class="n">mu0</span> <span class="o">=</span> <span class="mf">0.0</span>
|
||
<span class="n">sigma0</span> <span class="o">=</span> <span class="mf">1.0</span>
|
||
<span class="n">mu1</span><span class="o">=</span> <span class="mf">1.0</span>
|
||
<span class="n">sigma1</span> <span class="o">=</span> <span class="mf">2.0</span>
|
||
<span class="n">mu2</span> <span class="o">=</span> <span class="mf">2.0</span>
|
||
<span class="n">sigma2</span> <span class="o">=</span> <span class="mf">4.0</span>
|
||
|
||
<span class="n">x</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">linspace</span><span class="p">(</span><span class="o">-</span><span class="mf">20.0</span><span class="p">,</span> <span class="mf">20.0</span><span class="p">)</span>
|
||
<span class="n">v0</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">exp</span><span class="p">(</span><span class="o">-</span><span class="p">(</span><span class="n">x</span><span class="o">*</span><span class="n">x</span><span class="o">-</span><span class="mi">2</span><span class="o">*</span><span class="n">x</span><span class="o">*</span><span class="n">mu0</span><span class="o">+</span><span class="n">mu0</span><span class="o">*</span><span class="n">mu0</span><span class="p">)</span><span class="o">/</span><span class="p">(</span><span class="mi">2</span><span class="o">*</span><span class="n">sigma0</span><span class="o">*</span><span class="n">sigma0</span><span class="p">))</span><span class="o">/</span><span class="n">sqrt</span><span class="p">(</span><span class="mi">2</span><span class="o">*</span><span class="n">pi</span><span class="o">*</span><span class="n">sigma0</span><span class="o">*</span><span class="n">sigma0</span><span class="p">)</span>
|
||
<span class="n">v1</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">exp</span><span class="p">(</span><span class="o">-</span><span class="p">(</span><span class="n">x</span><span class="o">*</span><span class="n">x</span><span class="o">-</span><span class="mi">2</span><span class="o">*</span><span class="n">x</span><span class="o">*</span><span class="n">mu1</span><span class="o">+</span><span class="n">mu1</span><span class="o">*</span><span class="n">mu1</span><span class="p">)</span><span class="o">/</span><span class="p">(</span><span class="mi">2</span><span class="o">*</span><span class="n">sigma1</span><span class="o">*</span><span class="n">sigma1</span><span class="p">))</span><span class="o">/</span><span class="n">sqrt</span><span class="p">(</span><span class="mi">2</span><span class="o">*</span><span class="n">pi</span><span class="o">*</span><span class="n">sigma1</span><span class="o">*</span><span class="n">sigma1</span><span class="p">)</span>
|
||
<span class="n">v2</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">exp</span><span class="p">(</span><span class="o">-</span><span class="p">(</span><span class="n">x</span><span class="o">*</span><span class="n">x</span><span class="o">-</span><span class="mi">2</span><span class="o">*</span><span class="n">x</span><span class="o">*</span><span class="n">mu2</span><span class="o">+</span><span class="n">mu2</span><span class="o">*</span><span class="n">mu2</span><span class="p">)</span><span class="o">/</span><span class="p">(</span><span class="mi">2</span><span class="o">*</span><span class="n">sigma2</span><span class="o">*</span><span class="n">sigma2</span><span class="p">))</span><span class="o">/</span><span class="n">sqrt</span><span class="p">(</span><span class="mi">2</span><span class="o">*</span><span class="n">pi</span><span class="o">*</span><span class="n">sigma2</span><span class="o">*</span><span class="n">sigma2</span><span class="p">)</span>
|
||
<span class="n">plt</span><span class="o">.</span><span class="n">plot</span><span class="p">(</span><span class="n">x</span><span class="p">,</span> <span class="n">v0</span><span class="p">,</span> <span class="s1">'b-'</span><span class="p">,</span> <span class="n">x</span><span class="p">,</span> <span class="n">v1</span><span class="p">,</span> <span class="s1">'r-'</span><span class="p">,</span> <span class="n">x</span><span class="p">,</span> <span class="n">v2</span><span class="p">,</span> <span class="s1">'g-'</span><span class="p">)</span>
|
||
<span class="n">plt</span><span class="o">.</span><span class="n">title</span><span class="p">(</span><span class="sa">r</span><span class="s1">'{\bf Gaussian distributions}'</span><span class="p">,</span> <span class="n">fontsize</span><span class="o">=</span><span class="mi">20</span><span class="p">)</span>
|
||
<span class="n">plt</span><span class="o">.</span><span class="n">text</span><span class="p">(</span><span class="o">-</span><span class="mi">19</span><span class="p">,</span> <span class="mf">0.3</span><span class="p">,</span> <span class="sa">r</span><span class="s1">'Parameters: $\mu = 0$, $\sigma = 1$'</span><span class="p">,</span> <span class="n">fontdict</span><span class="o">=</span><span class="n">font</span><span class="p">)</span>
|
||
<span class="n">plt</span><span class="o">.</span><span class="n">text</span><span class="p">(</span><span class="o">-</span><span class="mi">19</span><span class="p">,</span> <span class="mf">0.18</span><span class="p">,</span> <span class="sa">r</span><span class="s1">'Parameters: $\mu = 1$, $\sigma = 2$'</span><span class="p">,</span> <span class="n">fontdict</span><span class="o">=</span><span class="n">font</span><span class="p">)</span>
|
||
<span class="n">plt</span><span class="o">.</span><span class="n">text</span><span class="p">(</span><span class="o">-</span><span class="mi">19</span><span class="p">,</span> <span class="mf">0.08</span><span class="p">,</span> <span class="sa">r</span><span class="s1">'Parameters: $\mu = 2$, $\sigma = 4$'</span><span class="p">,</span> <span class="n">fontdict</span><span class="o">=</span><span class="n">font</span><span class="p">)</span>
|
||
<span class="n">plt</span><span class="o">.</span><span class="n">xlabel</span><span class="p">(</span><span class="sa">r</span><span class="s1">'$x$'</span><span class="p">,</span><span class="n">fontsize</span><span class="o">=</span><span class="mi">20</span><span class="p">)</span>
|
||
<span class="n">plt</span><span class="o">.</span><span class="n">ylabel</span><span class="p">(</span><span class="sa">r</span><span class="s1">'$p(x)$ [MeV]'</span><span class="p">,</span><span class="n">fontsize</span><span class="o">=</span><span class="mi">20</span><span class="p">)</span>
|
||
|
||
<span class="c1"># Tweak spacing to prevent clipping of ylabel </span>
|
||
<span class="n">plt</span><span class="o">.</span><span class="n">subplots_adjust</span><span class="p">(</span><span class="n">left</span><span class="o">=</span><span class="mf">0.15</span><span class="p">)</span>
|
||
<span class="n">plt</span><span class="o">.</span><span class="n">savefig</span><span class="p">(</span><span class="s1">'gaussian.pdf'</span><span class="p">,</span> <span class="nb">format</span><span class="o">=</span><span class="s1">'pdf'</span><span class="p">)</span>
|
||
<span class="n">plt</span><span class="o">.</span><span class="n">show</span><span class="p">()</span>
|
||
</pre></div>
|
||
</div>
|
||
</div>
|
||
<div class="cell_output docutils container">
|
||
<img alt="_images/statistics_29_0.png" src="_images/statistics_29_0.png" />
|
||
</div>
|
||
</div>
|
||
</div>
|
||
<div class="section" id="exponential-distribution">
|
||
<h2><span class="section-number">1.10. </span>Exponential distribution<a class="headerlink" href="#exponential-distribution" title="Permalink to this headline">¶</a></h2>
|
||
<p>Another important distribution in science is the exponential distribution</p>
|
||
<div class="math notranslate nohighlight">
|
||
\[
|
||
p(x) = \alpha\exp{-(\alpha x)}.
|
||
\]</div>
|
||
</div>
|
||
<div class="section" id="expectation-values">
|
||
<h2><span class="section-number">1.11. </span>Expectation values<a class="headerlink" href="#expectation-values" title="Permalink to this headline">¶</a></h2>
|
||
<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>
|
||
</div>
|
||
<div class="section" id="stochastic-variables-and-the-main-concepts-mean-values">
|
||
<h2><span class="section-number">1.12. </span>Stochastic variables and the main concepts, mean values<a class="headerlink" href="#stochastic-variables-and-the-main-concepts-mean-values" title="Permalink to this headline">¶</a></h2>
|
||
<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>
|
||
</div>
|
||
<div class="section" id="stochastic-variables-and-the-main-concepts-central-moments-the-variance">
|
||
<h2><span class="section-number">1.13. </span>Stochastic variables and the main concepts, central moments, the variance<a class="headerlink" href="#stochastic-variables-and-the-main-concepts-central-moments-the-variance" title="Permalink to this headline">¶</a></h2>
|
||
<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 &=\mathrm{Var}(X) = \langle (x-\langle x\rangle)^2\rangle =
|
||
\int (x-\langle x\rangle)^2 p(x)dx\\
|
||
& = \int\left(x^2 - 2 x \langle x\rangle^{2} +\langle x\rangle^2\right)p(x)dx\\
|
||
& = \langle x^2\rangle - 2 \langle x\rangle\langle x\rangle + \langle x\rangle^2\\
|
||
& = \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>
|
||
</div>
|
||
<div class="section" id="probability-distribution-functions">
|
||
<h2><span class="section-number">1.14. </span>Probability Distribution Functions<a class="headerlink" href="#probability-distribution-functions" title="Permalink to this headline">¶</a></h2>
|
||
<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 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>
|
||
</div>
|
||
<div class="section" id="id1">
|
||
<h2><span class="section-number">1.15. </span>Probability Distribution Functions<a class="headerlink" href="#id1" title="Permalink to this headline">¶</a></h2>
|
||
<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>
|
||
</div>
|
||
<div class="section" id="the-three-famous-probability-distribution-functions">
|
||
<h2><span class="section-number">1.16. </span>The three famous Probability Distribution Functions<a class="headerlink" href="#the-three-famous-probability-distribution-functions" title="Permalink to this headline">¶</a></h2>
|
||
<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>
|
||
</div>
|
||
<div class="section" id="probability-distribution-functions-the-normal-distribution">
|
||
<h2><span class="section-number">1.17. </span>Probability Distribution Functions, the normal distribution<a class="headerlink" href="#probability-distribution-functions-the-normal-distribution" title="Permalink to this headline">¶</a></h2>
|
||
<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>
|
||
</div>
|
||
<div class="section" id="id2">
|
||
<h2><span class="section-number">1.18. </span>Probability Distribution Functions, the normal distribution<a class="headerlink" href="#id2" title="Permalink to this headline">¶</a></h2>
|
||
<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>
|
||
</div>
|
||
<div class="section" id="probability-distribution-functions-the-cumulative-distribution">
|
||
<h2><span class="section-number">1.19. </span>Probability Distribution Functions, the cumulative distribution<a class="headerlink" href="#probability-distribution-functions-the-cumulative-distribution" title="Permalink to this headline">¶</a></h2>
|
||
<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>
|
||
</div>
|
||
<div class="section" id="probability-distribution-functions-other-important-distribution">
|
||
<h2><span class="section-number">1.20. </span>Probability Distribution Functions, other important distribution<a class="headerlink" href="#probability-distribution-functions-other-important-distribution" title="Permalink to this headline">¶</a></h2>
|
||
<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>
|
||
</div>
|
||
<div class="section" id="probability-distribution-functions-the-binomial-distribution">
|
||
<h2><span class="section-number">1.21. </span>Probability Distribution Functions, the binomial distribution<a class="headerlink" href="#probability-distribution-functions-the-binomial-distribution" title="Permalink to this headline">¶</a></h2>
|
||
<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>
|
||
</div>
|
||
<div class="section" id="probability-distribution-functions-poisson-s-distribution">
|
||
<h2><span class="section-number">1.22. </span>Probability Distribution Functions, Poisson’s distribution<a class="headerlink" href="#probability-distribution-functions-poisson-s-distribution" title="Permalink to this headline">¶</a></h2>
|
||
<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 > 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>
|
||
</div>
|
||
<div class="section" id="id3">
|
||
<h2><span class="section-number">1.23. </span>Probability Distribution Functions, Poisson’s distribution<a class="headerlink" href="#id3" title="Permalink to this headline">¶</a></h2>
|
||
<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>
|
||
</div>
|
||
<div class="section" id="meet-the-covariance">
|
||
<h2><span class="section-number">1.24. </span>Meet the covariance!<a class="headerlink" href="#meet-the-covariance" title="Permalink to this headline">¶</a></h2>
|
||
<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>
|
||
</div>
|
||
<div class="section" id="meet-the-covariance-in-matrix-disguise">
|
||
<h2><span class="section-number">1.25. </span>Meet the covariance in matrix disguise<a class="headerlink" href="#meet-the-covariance-in-matrix-disguise" title="Permalink to this headline">¶</a></h2>
|
||
<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>
|
||
<div class="section" id="covariance">
|
||
<h2><span class="section-number">1.26. </span>Covariance<a class="headerlink" href="#covariance" title="Permalink to this headline">¶</a></h2>
|
||
<div class="cell docutils container">
|
||
<div class="cell_input docutils container">
|
||
<div class="highlight-ipython3 notranslate"><div class="highlight"><pre><span></span><span class="c1"># Importing various packages</span>
|
||
<span class="kn">from</span> <span class="nn">math</span> <span class="kn">import</span> <span class="n">exp</span><span class="p">,</span> <span class="n">sqrt</span>
|
||
<span class="kn">from</span> <span class="nn">random</span> <span class="kn">import</span> <span class="n">random</span><span class="p">,</span> <span class="n">seed</span>
|
||
<span class="kn">import</span> <span class="nn">numpy</span> <span class="k">as</span> <span class="nn">np</span>
|
||
<span class="kn">import</span> <span class="nn">matplotlib.pyplot</span> <span class="k">as</span> <span class="nn">plt</span>
|
||
|
||
<span class="k">def</span> <span class="nf">covariance</span><span class="p">(</span><span class="n">x</span><span class="p">,</span> <span class="n">y</span><span class="p">,</span> <span class="n">n</span><span class="p">):</span>
|
||
<span class="nb">sum</span> <span class="o">=</span> <span class="mf">0.0</span>
|
||
<span class="n">mean_x</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">mean</span><span class="p">(</span><span class="n">x</span><span class="p">)</span>
|
||
<span class="n">mean_y</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">mean</span><span class="p">(</span><span class="n">y</span><span class="p">)</span>
|
||
<span class="k">for</span> <span class="n">i</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="mi">0</span><span class="p">,</span> <span class="n">n</span><span class="p">):</span>
|
||
<span class="nb">sum</span> <span class="o">+=</span> <span class="p">(</span><span class="n">x</span><span class="p">[(</span><span class="n">i</span><span class="p">)]</span><span class="o">-</span><span class="n">mean_x</span><span class="p">)</span><span class="o">*</span><span class="p">(</span><span class="n">y</span><span class="p">[</span><span class="n">i</span><span class="p">]</span><span class="o">-</span><span class="n">mean_y</span><span class="p">)</span>
|
||
<span class="k">return</span> <span class="nb">sum</span><span class="o">/</span><span class="n">n</span>
|
||
|
||
<span class="n">n</span> <span class="o">=</span> <span class="mi">10</span>
|
||
|
||
<span class="n">x</span><span class="o">=</span><span class="n">np</span><span class="o">.</span><span class="n">random</span><span class="o">.</span><span class="n">normal</span><span class="p">(</span><span class="n">size</span><span class="o">=</span><span class="n">n</span><span class="p">)</span>
|
||
<span class="n">y</span> <span class="o">=</span> <span class="mi">4</span><span class="o">+</span><span class="mi">3</span><span class="o">*</span><span class="n">x</span><span class="o">+</span><span class="n">np</span><span class="o">.</span><span class="n">random</span><span class="o">.</span><span class="n">normal</span><span class="p">(</span><span class="n">size</span><span class="o">=</span><span class="n">n</span><span class="p">)</span>
|
||
<span class="n">covxy</span> <span class="o">=</span> <span class="n">covariance</span><span class="p">(</span><span class="n">x</span><span class="p">,</span><span class="n">y</span><span class="p">,</span><span class="n">n</span><span class="p">)</span>
|
||
<span class="nb">print</span><span class="p">(</span><span class="n">covxy</span><span class="p">)</span>
|
||
<span class="n">z</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">vstack</span><span class="p">((</span><span class="n">x</span><span class="p">,</span> <span class="n">y</span><span class="p">))</span>
|
||
<span class="n">c</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">cov</span><span class="p">(</span><span class="n">z</span><span class="o">.</span><span class="n">T</span><span class="p">)</span>
|
||
|
||
<span class="nb">print</span><span class="p">(</span><span class="n">c</span><span class="p">)</span>
|
||
</pre></div>
|
||
</div>
|
||
</div>
|
||
<div class="cell_output docutils container">
|
||
<div class="output stream highlight-myst-ansi notranslate"><div class="highlight"><pre><span></span>2.4890704720929593
|
||
[[11.60990408 7.82544198 14.11622885 1.78979284 6.76376813 14.09060397
|
||
14.04624055 3.77253612 5.11082097 14.99821162]
|
||
[ 7.82544198 5.27459502 9.51478402 1.20637689 4.55899331 9.49751204
|
||
9.4676097 2.54280848 3.44485473 10.10926826]
|
||
[14.11622885 9.51478402 17.16361442 2.17617003 8.22391797 17.13245767
|
||
17.07851716 4.58694429 6.21413561 18.23599801]
|
||
[ 1.78979284 1.20637689 2.17617003 0.27591601 1.04270833 2.17221967
|
||
2.16538057 0.58157742 0.78788857 2.31213725]
|
||
[ 6.76376813 4.55899331 8.22391797 1.04270833 3.94047694 8.20898927
|
||
8.18314376 2.19782691 2.97749299 8.73774883]
|
||
[14.09060397 9.49751204 17.13245767 2.17221967 8.20898927 17.10135749
|
||
17.04751489 4.57861771 6.20285522 18.20289459]
|
||
[14.04624055 9.4676097 17.07851716 2.16538057 8.18314376 17.04751489
|
||
16.99384182 4.56420221 6.18332591 18.14558387]
|
||
[ 3.77253612 2.54280848 4.58694429 0.58157742 2.19782691 4.57861771
|
||
4.56420221 1.2258524 1.66071628 4.87353683]
|
||
[ 5.11082097 3.44485473 6.21413561 0.78788857 2.97749299 6.20285522
|
||
6.18332591 1.66071628 2.24984554 6.60239515]
|
||
[14.99821162 10.10926826 18.23599801 2.31213725 8.73774883 18.20289459
|
||
18.14558387 4.87353683 6.60239515 19.37538418]]
|
||
</pre></div>
|
||
</div>
|
||
</div>
|
||
</div>
|
||
</div>
|
||
<div class="section" id="meet-the-covariance-uncorrelated-events">
|
||
<h2><span class="section-number">1.27. </span>Meet the covariance, uncorrelated events<a class="headerlink" href="#meet-the-covariance-uncorrelated-events" title="Permalink to this headline">¶</a></h2>
|
||
<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) &= \langle (x_i-\langle x_i\rangle)(x_j-\langle x_j\rangle)\rangle\\
|
||
&=\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\\
|
||
&=\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 \\
|
||
&=\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 \\
|
||
&=\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>
|
||
</div>
|
||
<div class="section" id="numerical-experiments-and-the-covariance">
|
||
<h2><span class="section-number">1.28. </span>Numerical experiments and the covariance<a class="headerlink" href="#numerical-experiments-and-the-covariance" title="Permalink to this headline">¶</a></h2>
|
||
<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>
|
||
</div>
|
||
<div class="section" id="id4">
|
||
<h2><span class="section-number">1.29. </span>Numerical experiments and the covariance<a class="headerlink" href="#id4" title="Permalink to this headline">¶</a></h2>
|
||
<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>
|
||
</div>
|
||
<div class="section" id="numerical-experiments-and-the-covariance-actual-situations">
|
||
<h2><span class="section-number">1.30. </span>Numerical experiments and the covariance, actual situations<a class="headerlink" href="#numerical-experiments-and-the-covariance-actual-situations" title="Permalink to this headline">¶</a></h2>
|
||
<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>
|
||
</div>
|
||
<div class="section" id="numerical-experiments-and-the-covariance-our-observables">
|
||
<h2><span class="section-number">1.31. </span>Numerical experiments and the covariance, our observables<a class="headerlink" href="#numerical-experiments-and-the-covariance-our-observables" title="Permalink to this headline">¶</a></h2>
|
||
<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>
|
||
</div>
|
||
<div class="section" id="numerical-experiments-and-the-covariance-the-sample-variance">
|
||
<h2><span class="section-number">1.32. </span>Numerical experiments and the covariance, the sample variance<a class="headerlink" href="#numerical-experiments-and-the-covariance-the-sample-variance" title="Permalink to this headline">¶</a></h2>
|
||
<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>
|
||
</div>
|
||
<div class="section" id="numerical-experiments-and-the-covariance-central-limit-theorem">
|
||
<h2><span class="section-number">1.33. </span>Numerical experiments and the covariance, central limit theorem<a class="headerlink" href="#numerical-experiments-and-the-covariance-central-limit-theorem" title="Permalink to this headline">¶</a></h2>
|
||
<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>
|
||
</div>
|
||
<div class="section" id="definition-of-correlation-functions-and-standard-deviation">
|
||
<h2><span class="section-number">1.34. </span>Definition of Correlation Functions and Standard Deviation<a class="headerlink" href="#definition-of-correlation-functions-and-standard-deviation" title="Permalink to this headline">¶</a></h2>
|
||
<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 external" href="#eq:exptvariance">7</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<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>
|
||
</div>
|
||
<div class="section" id="id5">
|
||
<h2><span class="section-number">1.35. </span>Definition of Correlation Functions and Standard Deviation<a class="headerlink" href="#id5" title="Permalink to this headline">¶</a></h2>
|
||
<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>
|
||
</div>
|
||
<div class="section" id="id6">
|
||
<h2><span class="section-number">1.36. </span>Definition of Correlation Functions and Standard Deviation<a class="headerlink" href="#id6" title="Permalink to this headline">¶</a></h2>
|
||
<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<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>
|
||
</div>
|
||
<div class="section" id="id7">
|
||
<h2><span class="section-number">1.37. </span>Definition of Correlation Functions and Standard Deviation<a class="headerlink" href="#id7" title="Permalink to this headline">¶</a></h2>
|
||
<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>
|
||
</div>
|
||
<div class="section" id="definition-of-correlation-functions-and-standard-deviation-sample-variance">
|
||
<h2><span class="section-number">1.38. </span>Definition of Correlation Functions and Standard Deviation, sample variance<a class="headerlink" href="#definition-of-correlation-functions-and-standard-deviation-sample-variance" title="Permalink to this headline">¶</a></h2>
|
||
<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>
|
||
</div>
|
||
<div class="section" id="id8">
|
||
<h2><span class="section-number">1.39. </span>Definition of Correlation Functions and Standard Deviation<a class="headerlink" href="#id8" title="Permalink to this headline">¶</a></h2>
|
||
<p>From the point of view of
|
||
Eq. (<a class="reference external" href="#eq:error_estimate_corr_time">10</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>
|
||
<div class="section" id="code-to-compute-the-covariance-matrix-and-the-covariance">
|
||
<h2><span class="section-number">1.40. </span>Code to compute the Covariance matrix and the Covariance<a class="headerlink" href="#code-to-compute-the-covariance-matrix-and-the-covariance" title="Permalink to this headline">¶</a></h2>
|
||
<div class="cell docutils container">
|
||
<div class="cell_input docutils container">
|
||
<div class="highlight-ipython3 notranslate"><div class="highlight"><pre><span></span><span class="c1"># Importing various packages</span>
|
||
<span class="kn">from</span> <span class="nn">math</span> <span class="kn">import</span> <span class="n">exp</span><span class="p">,</span> <span class="n">sqrt</span>
|
||
<span class="kn">from</span> <span class="nn">random</span> <span class="kn">import</span> <span class="n">random</span><span class="p">,</span> <span class="n">seed</span>
|
||
<span class="kn">import</span> <span class="nn">numpy</span> <span class="k">as</span> <span class="nn">np</span>
|
||
<span class="kn">import</span> <span class="nn">matplotlib.pyplot</span> <span class="k">as</span> <span class="nn">plt</span>
|
||
|
||
<span class="c1"># Sample covariance, note the factor 1/(n-1)</span>
|
||
<span class="k">def</span> <span class="nf">covariance</span><span class="p">(</span><span class="n">x</span><span class="p">,</span> <span class="n">y</span><span class="p">,</span> <span class="n">n</span><span class="p">):</span>
|
||
<span class="nb">sum</span> <span class="o">=</span> <span class="mf">0.0</span>
|
||
<span class="n">mean_x</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">mean</span><span class="p">(</span><span class="n">x</span><span class="p">)</span>
|
||
<span class="n">mean_y</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">mean</span><span class="p">(</span><span class="n">y</span><span class="p">)</span>
|
||
<span class="k">for</span> <span class="n">i</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="mi">0</span><span class="p">,</span> <span class="n">n</span><span class="p">):</span>
|
||
<span class="nb">sum</span> <span class="o">+=</span> <span class="p">(</span><span class="n">x</span><span class="p">[(</span><span class="n">i</span><span class="p">)]</span><span class="o">-</span><span class="n">mean_x</span><span class="p">)</span><span class="o">*</span><span class="p">(</span><span class="n">y</span><span class="p">[</span><span class="n">i</span><span class="p">]</span><span class="o">-</span><span class="n">mean_y</span><span class="p">)</span>
|
||
<span class="k">return</span> <span class="nb">sum</span><span class="o">/</span><span class="p">(</span><span class="n">n</span><span class="o">-</span><span class="mf">1.</span><span class="p">)</span>
|
||
|
||
<span class="n">n</span> <span class="o">=</span> <span class="mi">100</span>
|
||
<span class="n">x</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">random</span><span class="o">.</span><span class="n">normal</span><span class="p">(</span><span class="n">size</span><span class="o">=</span><span class="n">n</span><span class="p">)</span>
|
||
<span class="nb">print</span><span class="p">(</span><span class="n">np</span><span class="o">.</span><span class="n">mean</span><span class="p">(</span><span class="n">x</span><span class="p">))</span>
|
||
<span class="n">y</span> <span class="o">=</span> <span class="mi">4</span><span class="o">+</span><span class="mi">3</span><span class="o">*</span><span class="n">x</span><span class="o">+</span><span class="n">np</span><span class="o">.</span><span class="n">random</span><span class="o">.</span><span class="n">normal</span><span class="p">(</span><span class="n">size</span><span class="o">=</span><span class="n">n</span><span class="p">)</span>
|
||
<span class="nb">print</span><span class="p">(</span><span class="n">np</span><span class="o">.</span><span class="n">mean</span><span class="p">(</span><span class="n">y</span><span class="p">))</span>
|
||
<span class="n">z</span> <span class="o">=</span> <span class="n">x</span><span class="o">**</span><span class="mi">3</span><span class="o">+</span><span class="n">np</span><span class="o">.</span><span class="n">random</span><span class="o">.</span><span class="n">normal</span><span class="p">(</span><span class="n">size</span><span class="o">=</span><span class="n">n</span><span class="p">)</span>
|
||
<span class="nb">print</span><span class="p">(</span><span class="n">np</span><span class="o">.</span><span class="n">mean</span><span class="p">(</span><span class="n">z</span><span class="p">))</span>
|
||
<span class="n">covxx</span> <span class="o">=</span> <span class="n">covariance</span><span class="p">(</span><span class="n">x</span><span class="p">,</span><span class="n">x</span><span class="p">,</span><span class="n">n</span><span class="p">)</span>
|
||
<span class="n">covyy</span> <span class="o">=</span> <span class="n">covariance</span><span class="p">(</span><span class="n">y</span><span class="p">,</span><span class="n">y</span><span class="p">,</span><span class="n">n</span><span class="p">)</span>
|
||
<span class="n">covzz</span> <span class="o">=</span> <span class="n">covariance</span><span class="p">(</span><span class="n">z</span><span class="p">,</span><span class="n">z</span><span class="p">,</span><span class="n">n</span><span class="p">)</span>
|
||
<span class="n">covxy</span> <span class="o">=</span> <span class="n">covariance</span><span class="p">(</span><span class="n">x</span><span class="p">,</span><span class="n">y</span><span class="p">,</span><span class="n">n</span><span class="p">)</span>
|
||
<span class="n">covxz</span> <span class="o">=</span> <span class="n">covariance</span><span class="p">(</span><span class="n">x</span><span class="p">,</span><span class="n">z</span><span class="p">,</span><span class="n">n</span><span class="p">)</span>
|
||
<span class="n">covyz</span> <span class="o">=</span> <span class="n">covariance</span><span class="p">(</span><span class="n">y</span><span class="p">,</span><span class="n">z</span><span class="p">,</span><span class="n">n</span><span class="p">)</span>
|
||
<span class="nb">print</span><span class="p">(</span><span class="n">covxx</span><span class="p">,</span><span class="n">covyy</span><span class="p">,</span> <span class="n">covzz</span><span class="p">)</span>
|
||
<span class="nb">print</span><span class="p">(</span><span class="n">covxy</span><span class="p">,</span><span class="n">covxz</span><span class="p">,</span> <span class="n">covyz</span><span class="p">)</span>
|
||
<span class="n">w</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">vstack</span><span class="p">((</span><span class="n">x</span><span class="p">,</span> <span class="n">y</span><span class="p">,</span> <span class="n">z</span><span class="p">))</span>
|
||
<span class="c1">#print(w)</span>
|
||
<span class="n">c</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">cov</span><span class="p">(</span><span class="n">w</span><span class="p">)</span>
|
||
<span class="nb">print</span><span class="p">(</span><span class="n">c</span><span class="p">)</span>
|
||
<span class="c1">#eigen = np.zeros(n)</span>
|
||
<span class="n">Eigvals</span><span class="p">,</span> <span class="n">Eigvecs</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">linalg</span><span class="o">.</span><span class="n">eig</span><span class="p">(</span><span class="n">c</span><span class="p">)</span>
|
||
<span class="nb">print</span><span class="p">(</span><span class="n">Eigvals</span><span class="p">)</span>
|
||
</pre></div>
|
||
</div>
|
||
</div>
|
||
<div class="cell_output docutils container">
|
||
<div class="output stream highlight-myst-ansi notranslate"><div class="highlight"><pre><span></span>-0.000171990965637209
|
||
3.967581388294481
|
||
-0.025521898567525465
|
||
0.853884005749657 8.208733310166176 9.244798847617261
|
||
2.497355838775297 2.088848443911507 6.174325810031482
|
||
[[0.85388401 2.49735584 2.08884844]
|
||
[2.49735584 8.20873331 6.17432581]
|
||
[2.08884844 6.17432581 9.24479885]]
|
||
[15.62869983 0.07851446 2.60020187]
|
||
</pre></div>
|
||
</div>
|
||
</div>
|
||
</div>
|
||
</div>
|
||
<div class="section" id="random-numbers">
|
||
<h2><span class="section-number">1.41. </span>Random Numbers<a class="headerlink" href="#random-numbers" title="Permalink to this headline">¶</a></h2>
|
||
<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>
|
||
</div>
|
||
<div class="section" id="random-numbers-better-name-pseudo-random-numbers">
|
||
<h2><span class="section-number">1.42. </span>Random Numbers, better name: pseudo random numbers<a class="headerlink" href="#random-numbers-better-name-pseudo-random-numbers" title="Permalink to this headline">¶</a></h2>
|
||
<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>
|
||
</div>
|
||
<div class="section" id="random-number-generator-rng">
|
||
<h2><span class="section-number">1.43. </span>Random number generator RNG<a class="headerlink" href="#random-number-generator-rng" title="Permalink to this headline">¶</a></h2>
|
||
<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>
|
||
</div>
|
||
<div class="section" id="random-number-generator-rng-and-periodic-outputs">
|
||
<h2><span class="section-number">1.44. </span>Random number generator RNG and periodic outputs<a class="headerlink" href="#random-number-generator-rng-and-periodic-outputs" title="Permalink to this headline">¶</a></h2>
|
||
<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>
|
||
</div>
|
||
<div class="section" id="random-number-generator-rng-and-its-period">
|
||
<h2><span class="section-number">1.45. </span>Random number generator RNG and its period<a class="headerlink" href="#random-number-generator-rng-and-its-period" title="Permalink to this headline">¶</a></h2>
|
||
<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>
|
||
</div>
|
||
<div class="section" id="random-number-generator-rng-other-examples">
|
||
<h2><span class="section-number">1.46. </span>Random number generator RNG, other examples<a class="headerlink" href="#random-number-generator-rng-other-examples" title="Permalink to this headline">¶</a></h2>
|
||
<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>
|
||
</div>
|
||
<div class="section" id="id9">
|
||
<h2><span class="section-number">1.47. </span>Random number generator RNG, other examples<a class="headerlink" href="#id9" title="Permalink to this headline">¶</a></h2>
|
||
<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>
|
||
</div>
|
||
<div class="section" id="random-number-generator-rng-ran0">
|
||
<h2><span class="section-number">1.48. </span>Random number generator RNG, RAN0<a class="headerlink" href="#random-number-generator-rng-ran0" title="Permalink to this headline">¶</a></h2>
|
||
<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 < q\)</span>.</p>
|
||
</div>
|
||
<div class="section" id="id10">
|
||
<h2><span class="section-number">1.49. </span>Random number generator RNG, RAN0<a class="headerlink" href="#id10" title="Permalink to this headline">¶</a></h2>
|
||
<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>
|
||
</div>
|
||
<div class="section" id="id11">
|
||
<h2><span class="section-number">1.50. </span>Random number generator RNG, RAN0<a class="headerlink" href="#id11" title="Permalink to this headline">¶</a></h2>
|
||
<p>We can now rewrite Eq. (<a class="reference external" href="#eq:rntrick1">14</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
|
||
in</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>
|
||
</div>
|
||
<div class="section" id="id12">
|
||
<h2><span class="section-number">1.51. </span>Random number generator RNG, RAN0<a class="headerlink" href="#id12" title="Permalink to this headline">¶</a></h2>
|
||
<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 < 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))< 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>
|
||
</div>
|
||
<div class="section" id="random-number-generator-rng-ran0-code">
|
||
<h2><span class="section-number">1.52. </span>Random number generator RNG, RAN0 code<a class="headerlink" href="#random-number-generator-rng-ran0-code" title="Permalink to this headline">¶</a></h2>
|
||
<div class="highlight-none notranslate"><div class="highlight"><pre><span></span> /*
|
||
** The function
|
||
** ran0()
|
||
** is an "Minimal" random number generator of Park and Miller
|
||
** Set or reset the input value
|
||
** idum to any integer value (except the unlikely value MASK)
|
||
** to initialize the sequence; idum must not be altered between
|
||
** calls for sucessive deviates in a sequence.
|
||
** The function returns a uniform deviate between 0.0 and 1.0.
|
||
*/
|
||
double ran0(long &idum)
|
||
{
|
||
const int a = 16807, m = 2147483647, q = 127773;
|
||
const int r = 2836, MASK = 123459876;
|
||
const double am = 1./m;
|
||
long k;
|
||
double ans;
|
||
idum ^= MASK;
|
||
k = (*idum)/q;
|
||
idum = a*(idum - k*q) - r*k;
|
||
// add m if negative difference
|
||
if(idum < 0) idum += m;
|
||
ans=am*(idum);
|
||
idum ^= MASK;
|
||
return ans;
|
||
} // End: function ran0()
|
||
</pre></div>
|
||
</div>
|
||
</div>
|
||
<div class="section" id="properties-of-selected-random-number-generators">
|
||
<h2><span class="section-number">1.53. </span>Properties of Selected Random Number Generators<a class="headerlink" href="#properties-of-selected-random-number-generators" title="Permalink to this headline">¶</a></h2>
|
||
<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>
|
||
</div>
|
||
<div class="section" id="id13">
|
||
<h2><span class="section-number">1.54. </span>Properties of Selected Random Number Generators<a class="headerlink" href="#id13" title="Permalink to this headline">¶</a></h2>
|
||
<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>
|
||
</div>
|
||
<div class="section" id="id14">
|
||
<h2><span class="section-number">1.55. </span>Properties of Selected Random Number Generators<a class="headerlink" href="#id14" title="Permalink to this headline">¶</a></h2>
|
||
<p>The various random number generators produce results which agree rather well with
|
||
these limiting values.</p>
|
||
<table 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>
|
||
</div>
|
||
<div class="section" id="simple-demonstration-of-rngs-using-python">
|
||
<h2><span class="section-number">1.56. </span>Simple demonstration of RNGs using python<a class="headerlink" href="#simple-demonstration-of-rngs-using-python" title="Permalink to this headline">¶</a></h2>
|
||
<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-ipython3 notranslate"><div class="highlight"><pre><span></span><span class="ch">#!/usr/bin/env python</span>
|
||
<span class="kn">import</span> <span class="nn">numpy</span> <span class="k">as</span> <span class="nn">np</span>
|
||
<span class="kn">import</span> <span class="nn">matplotlib.mlab</span> <span class="k">as</span> <span class="nn">mlab</span>
|
||
<span class="kn">import</span> <span class="nn">matplotlib.pyplot</span> <span class="k">as</span> <span class="nn">plt</span>
|
||
<span class="kn">import</span> <span class="nn">random</span>
|
||
|
||
<span class="c1"># initialize the rng with a seed</span>
|
||
<span class="n">random</span><span class="o">.</span><span class="n">seed</span><span class="p">()</span>
|
||
<span class="n">counts</span> <span class="o">=</span> <span class="mi">10000</span>
|
||
<span class="n">values</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">zeros</span><span class="p">(</span><span class="n">counts</span><span class="p">)</span>
|
||
<span class="k">for</span> <span class="n">i</span> <span class="ow">in</span> <span class="nb">range</span> <span class="p">(</span><span class="mi">1</span><span class="p">,</span> <span class="n">counts</span><span class="p">,</span> <span class="mi">1</span><span class="p">):</span>
|
||
<span class="n">values</span><span class="p">[</span><span class="n">i</span><span class="p">]</span> <span class="o">=</span> <span class="n">random</span><span class="o">.</span><span class="n">random</span><span class="p">()</span>
|
||
|
||
<span class="c1"># the histogram of the data</span>
|
||
<span class="n">n</span><span class="p">,</span> <span class="n">bins</span><span class="p">,</span> <span class="n">patches</span> <span class="o">=</span> <span class="n">plt</span><span class="o">.</span><span class="n">hist</span><span class="p">(</span><span class="n">values</span><span class="p">,</span> <span class="mi">10</span><span class="p">,</span> <span class="n">facecolor</span><span class="o">=</span><span class="s1">'green'</span><span class="p">)</span>
|
||
|
||
<span class="n">plt</span><span class="o">.</span><span class="n">xlabel</span><span class="p">(</span><span class="s1">'$x$'</span><span class="p">)</span>
|
||
<span class="n">plt</span><span class="o">.</span><span class="n">ylabel</span><span class="p">(</span><span class="s1">'Number of counts'</span><span class="p">)</span>
|
||
<span class="n">plt</span><span class="o">.</span><span class="n">title</span><span class="p">(</span><span class="sa">r</span><span class="s1">'Test of uniform distribution'</span><span class="p">)</span>
|
||
<span class="n">plt</span><span class="o">.</span><span class="n">axis</span><span class="p">([</span><span class="mi">0</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="mi">1100</span><span class="p">])</span>
|
||
<span class="n">plt</span><span class="o">.</span><span class="n">grid</span><span class="p">(</span><span class="kc">True</span><span class="p">)</span>
|
||
<span class="n">plt</span><span class="o">.</span><span class="n">show</span><span class="p">()</span>
|
||
</pre></div>
|
||
</div>
|
||
</div>
|
||
<div class="cell_output docutils container">
|
||
<img alt="_images/statistics_178_0.png" src="_images/statistics_178_0.png" />
|
||
</div>
|
||
</div>
|
||
</div>
|
||
<div class="section" id="id15">
|
||
<h2><span class="section-number">1.57. </span>Properties of Selected Random Number Generators<a class="headerlink" href="#id15" title="Permalink to this headline">¶</a></h2>
|
||
<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 external" href="#eq:autocorrelformal">9</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>
|
||
</div>
|
||
<div class="section" id="autocorrelation-function">
|
||
<h2><span class="section-number">1.58. </span>Autocorrelation function<a class="headerlink" href="#autocorrelation-function" title="Permalink to this headline">¶</a></h2>
|
||
<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-ipython3 notranslate"><div class="highlight"><pre><span></span><span class="c1"># Importing various packages</span>
|
||
<span class="kn">from</span> <span class="nn">math</span> <span class="kn">import</span> <span class="n">exp</span><span class="p">,</span> <span class="n">sqrt</span>
|
||
<span class="kn">from</span> <span class="nn">random</span> <span class="kn">import</span> <span class="n">random</span><span class="p">,</span> <span class="n">seed</span>
|
||
<span class="kn">import</span> <span class="nn">numpy</span> <span class="k">as</span> <span class="nn">np</span>
|
||
<span class="kn">import</span> <span class="nn">matplotlib.pyplot</span> <span class="k">as</span> <span class="nn">plt</span>
|
||
|
||
<span class="k">def</span> <span class="nf">autocovariance</span><span class="p">(</span><span class="n">x</span><span class="p">,</span> <span class="n">n</span><span class="p">,</span> <span class="n">k</span><span class="p">,</span> <span class="n">mean_x</span><span class="p">):</span>
|
||
<span class="nb">sum</span> <span class="o">=</span> <span class="mf">0.0</span>
|
||
<span class="k">for</span> <span class="n">i</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="mi">0</span><span class="p">,</span> <span class="n">n</span><span class="o">-</span><span class="n">k</span><span class="p">):</span>
|
||
<span class="nb">sum</span> <span class="o">+=</span> <span class="p">(</span><span class="n">x</span><span class="p">[(</span><span class="n">i</span><span class="o">+</span><span class="n">k</span><span class="p">)]</span><span class="o">-</span><span class="n">mean_x</span><span class="p">)</span><span class="o">*</span><span class="p">(</span><span class="n">x</span><span class="p">[</span><span class="n">i</span><span class="p">]</span><span class="o">-</span><span class="n">mean_x</span><span class="p">)</span>
|
||
<span class="k">return</span> <span class="nb">sum</span><span class="o">/</span><span class="n">n</span>
|
||
|
||
<span class="n">n</span> <span class="o">=</span> <span class="mi">1000</span>
|
||
<span class="n">x</span><span class="o">=</span><span class="n">np</span><span class="o">.</span><span class="n">random</span><span class="o">.</span><span class="n">normal</span><span class="p">(</span><span class="n">size</span><span class="o">=</span><span class="n">n</span><span class="p">)</span>
|
||
<span class="n">autocor</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">zeros</span><span class="p">(</span><span class="n">n</span><span class="p">)</span>
|
||
<span class="n">figaxis</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">zeros</span><span class="p">(</span><span class="n">n</span><span class="p">)</span>
|
||
<span class="n">mean_x</span><span class="o">=</span><span class="n">np</span><span class="o">.</span><span class="n">mean</span><span class="p">(</span><span class="n">x</span><span class="p">)</span>
|
||
<span class="n">var_x</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">var</span><span class="p">(</span><span class="n">x</span><span class="p">)</span>
|
||
<span class="nb">print</span><span class="p">(</span><span class="n">mean_x</span><span class="p">,</span> <span class="n">var_x</span><span class="p">)</span>
|
||
<span class="k">for</span> <span class="n">i</span> <span class="ow">in</span> <span class="nb">range</span> <span class="p">(</span><span class="mi">0</span><span class="p">,</span> <span class="n">n</span><span class="p">):</span>
|
||
<span class="n">figaxis</span><span class="p">[</span><span class="n">i</span><span class="p">]</span> <span class="o">=</span> <span class="n">i</span>
|
||
<span class="n">autocor</span><span class="p">[</span><span class="n">i</span><span class="p">]</span><span class="o">=</span><span class="p">(</span><span class="n">autocovariance</span><span class="p">(</span><span class="n">x</span><span class="p">,</span> <span class="n">n</span><span class="p">,</span> <span class="n">i</span><span class="p">,</span> <span class="n">mean_x</span><span class="p">))</span><span class="o">/</span><span class="n">var_x</span>
|
||
|
||
<span class="n">plt</span><span class="o">.</span><span class="n">plot</span><span class="p">(</span><span class="n">figaxis</span><span class="p">,</span> <span class="n">autocor</span><span class="p">,</span> <span class="s2">"r-"</span><span class="p">)</span>
|
||
<span class="n">plt</span><span class="o">.</span><span class="n">axis</span><span class="p">([</span><span class="mi">0</span><span class="p">,</span><span class="n">n</span><span class="p">,</span><span class="o">-</span><span class="mf">0.1</span><span class="p">,</span> <span class="mf">1.0</span><span class="p">])</span>
|
||
<span class="n">plt</span><span class="o">.</span><span class="n">xlabel</span><span class="p">(</span><span class="sa">r</span><span class="s1">'$i$'</span><span class="p">)</span>
|
||
<span class="n">plt</span><span class="o">.</span><span class="n">ylabel</span><span class="p">(</span><span class="sa">r</span><span class="s1">'$\gamma_i$'</span><span class="p">)</span>
|
||
<span class="n">plt</span><span class="o">.</span><span class="n">title</span><span class="p">(</span><span class="sa">r</span><span class="s1">'Autocorrelation function'</span><span class="p">)</span>
|
||
<span class="n">plt</span><span class="o">.</span><span class="n">show</span><span class="p">()</span>
|
||
</pre></div>
|
||
</div>
|
||
</div>
|
||
<div class="cell_output docutils container">
|
||
<div class="output stream highlight-myst-ansi notranslate"><div class="highlight"><pre><span></span>0.005481494389472717 1.029117514654528
|
||
</pre></div>
|
||
</div>
|
||
<img alt="_images/statistics_184_1.png" src="_images/statistics_184_1.png" />
|
||
</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>
|
||
</div>
|
||
<div class="section" id="correlation-function-and-which-random-number-generators-should-i-use">
|
||
<h2><span class="section-number">1.59. </span>Correlation function and which random number generators should I use<a class="headerlink" href="#correlation-function-and-which-random-number-generators-should-i-use" title="Permalink to this headline">¶</a></h2>
|
||
<p>The program here computes the correlation function for one of the standard functions included with the c++ compiler.</p>
|
||
<div class="highlight-none notranslate"><div class="highlight"><pre><span></span> // This function computes the autocorrelation function for
|
||
// the standard c++ random number generator
|
||
|
||
#include <fstream>
|
||
#include <iomanip>
|
||
#include <iostream>
|
||
#include <cmath>
|
||
using namespace std;
|
||
// output file as global variable
|
||
ofstream ofile;
|
||
|
||
// Main function begins here
|
||
int main(int argc, char* argv[])
|
||
{
|
||
int n;
|
||
char *outfilename;
|
||
|
||
cin >> n;
|
||
double MCint = 0.; double MCintsqr2=0.;
|
||
double invers_period = 1./RAND_MAX; // initialise the random number generator
|
||
srand(time(NULL)); // This produces the so-called seed in MC jargon
|
||
// Compute the variance and the mean value of the uniform distribution
|
||
// Compute also the specific values x for each cycle in order to be able to
|
||
// the covariance and the correlation function
|
||
// Read in output file, abort if there are too few command-line arguments
|
||
if( argc <= 2 ){
|
||
cout << "Bad Usage: " << argv[0] <<
|
||
" read also output file and number of cycles on same line" << endl;
|
||
exit(1);
|
||
}
|
||
else{
|
||
outfilename=argv[1];
|
||
}
|
||
ofile.open(outfilename);
|
||
// Get the number of Monte-Carlo samples
|
||
n = atoi(argv[2]);
|
||
double *X;
|
||
X = new double[n];
|
||
for (int i = 0; i < n; i++){
|
||
double x = double(rand())*invers_period;
|
||
X[i] = x;
|
||
MCint += x;
|
||
MCintsqr2 += x*x;
|
||
}
|
||
double Mean = MCint/((double) n );
|
||
MCintsqr2 = MCintsqr2/((double) n );
|
||
double STDev = sqrt(MCintsqr2-Mean*Mean);
|
||
double Variance = MCintsqr2-Mean*Mean;
|
||
// Write mean value and standard deviation
|
||
cout << " Standard deviation= " << STDev << " Integral = " << Mean << endl;
|
||
|
||
// Now we compute the autocorrelation function
|
||
double *autocor; autocor = new double[n];
|
||
for (int j = 0; j < n; j++){
|
||
double sum = 0.0;
|
||
for (int k = 0; k < (n-j); k++){
|
||
sum += (X[k]-Mean)*(X[k+j]-Mean);
|
||
}
|
||
autocor[j] = sum/Variance/((double) n );
|
||
ofile << setiosflags(ios::showpoint | ios::uppercase);
|
||
ofile << setw(15) << setprecision(8) << j;
|
||
ofile << setw(15) << setprecision(8) << autocor[j] << endl;
|
||
}
|
||
ofile.close(); // close output file
|
||
return 0;
|
||
} // end of main program
|
||
</pre></div>
|
||
</div>
|
||
</div>
|
||
<div class="section" id="which-rng-should-i-use">
|
||
<h2><span class="section-number">1.60. </span>Which RNG should I use?<a class="headerlink" href="#which-rng-should-i-use" title="Permalink to this headline">¶</a></h2>
|
||
<ul class="simple">
|
||
<li><p>C++ has a class called <strong>random</strong>. The <a class="reference external" href="http://www.cplusplus.com/reference/random/">random class</a> contains a large selection of RNGs and is highly recommended. Some of these RNGs have very large periods making it thereby very safe to use these RNGs in case one is performing large calculations. In particular, the <a class="reference external" href="http://www.cplusplus.com/reference/random/mersenne_twister_engine/">Mersenne twister random number engine</a> has a period of <span class="math notranslate nohighlight">\(2^{19937}\)</span>.</p></li>
|
||
<li><p>Add RNGs in Python</p></li>
|
||
</ul>
|
||
</div>
|
||
<div class="section" id="how-to-use-the-mersenne-generator">
|
||
<h2><span class="section-number">1.61. </span>How to use the Mersenne generator<a class="headerlink" href="#how-to-use-the-mersenne-generator" title="Permalink to this headline">¶</a></h2>
|
||
<p>The following part of a c++ code (from project 4) sets up the uniform distribution for <span class="math notranslate nohighlight">\(x\in [0,1]\)</span>.</p>
|
||
<div class="highlight-none notranslate"><div class="highlight"><pre><span></span> /*
|
||
|
||
// You need this
|
||
#include <random>
|
||
|
||
// Initialize the seed and call the Mersienne algo
|
||
std::random_device rd;
|
||
std::mt19937_64 gen(rd());
|
||
// Set up the uniform distribution for x \in [[0, 1]
|
||
std::uniform_real_distribution<double> RandomNumberGenerator(0.0,1.0);
|
||
|
||
// Now use the RNG
|
||
int ix = (int) (RandomNumberGenerator(gen)*NSpins);
|
||
</pre></div>
|
||
</div>
|
||
</div>
|
||
<div class="section" id="why-blocking">
|
||
<h2><span class="section-number">1.62. </span>Why blocking?<a class="headerlink" href="#why-blocking" title="Permalink to this headline">¶</a></h2>
|
||
<p><strong>Statistical analysis.</strong></p>
|
||
<div class="highlight-none notranslate"><div class="highlight"><pre><span></span>* Monte Carlo simulations can be treated as *computer experiments*
|
||
|
||
* The results can be analysed with the same statistical tools as we would use analysing experimental data.
|
||
|
||
* As in all experiments, we are looking for expectation values and an estimate of how accurate they are, i.e., possible sources for errors.
|
||
</pre></div>
|
||
</div>
|
||
<p>A very good article which explains blocking is H. Flyvbjerg and H. G. Petersen, <em>Error estimates on averages of correlated data</em>, <a class="reference external" href="http://scitation.aip.org/content/aip/journal/jcp/91/1/10.1063/1.457480">Journal of Chemical Physics 91, 461-466 (1989)</a>.</p>
|
||
</div>
|
||
<div class="section" id="id16">
|
||
<h2><span class="section-number">1.63. </span>Why blocking?<a class="headerlink" href="#id16" title="Permalink to this headline">¶</a></h2>
|
||
<p><strong>Statistical analysis.</strong></p>
|
||
<div class="highlight-none notranslate"><div class="highlight"><pre><span></span>* As in other experiments, Monte Carlo experiments have two classes of errors:
|
||
|
||
* Statistical errors
|
||
|
||
* Systematical errors
|
||
|
||
|
||
* Statistical errors can be estimated using standard tools from statistics
|
||
|
||
* Systematical errors are method specific and must be treated differently from case to case. (In VMC a common source is the step length or time step in importance sampling)
|
||
</pre></div>
|
||
</div>
|
||
</div>
|
||
<div class="section" id="code-to-demonstrate-the-calculation-of-the-autocorrelation-function">
|
||
<h2><span class="section-number">1.64. </span>Code to demonstrate the calculation of the autocorrelation function<a class="headerlink" href="#code-to-demonstrate-the-calculation-of-the-autocorrelation-function" title="Permalink to this headline">¶</a></h2>
|
||
<p>The following code computes the autocorrelation function, the covariance and the standard deviation
|
||
for standard RNG.
|
||
The <a class="reference external" href="https://github.com/CompPhysics/ComputationalPhysics2/tree/gh-pages/doc/Programs/LecturePrograms/programs/Blocking/autocorrelation.cpp">following file</a> gives the code.</p>
|
||
<div class="highlight-none notranslate"><div class="highlight"><pre><span></span> // This function computes the autocorrelation function for
|
||
// the Mersenne random number generator with a uniform distribution
|
||
#include <iostream>
|
||
#include <fstream>
|
||
#include <iomanip>
|
||
#include <cstdlib>
|
||
#include <random>
|
||
#include <armadillo>
|
||
#include <string>
|
||
#include <cmath>
|
||
using namespace std;
|
||
using namespace arma;
|
||
// output file
|
||
ofstream ofile;
|
||
|
||
// Main function begins here
|
||
int main(int argc, char* argv[])
|
||
{
|
||
int MonteCarloCycles;
|
||
string filename;
|
||
if (argc > 1) {
|
||
filename=argv[1];
|
||
MonteCarloCycles = atoi(argv[2]);
|
||
string fileout = filename;
|
||
string argument = to_string(MonteCarloCycles);
|
||
fileout.append(argument);
|
||
ofile.open(fileout);
|
||
}
|
||
|
||
// Compute the variance and the mean value of the uniform distribution
|
||
// Compute also the specific values x for each cycle in order to be able to
|
||
// compute the covariance and the correlation function
|
||
|
||
vec X = zeros<vec>(MonteCarloCycles);
|
||
double MCint = 0.; double MCintsqr2=0.;
|
||
std::random_device rd;
|
||
std::mt19937_64 gen(rd());
|
||
// Set up the uniform distribution for x \in [[0, 1]
|
||
std::uniform_real_distribution<double> RandomNumberGenerator(0.0,1.0);
|
||
for (int i = 0; i < MonteCarloCycles; i++){
|
||
double x = RandomNumberGenerator(gen);
|
||
X(i) = x;
|
||
MCint += x;
|
||
MCintsqr2 += x*x;
|
||
}
|
||
double Mean = MCint/((double) MonteCarloCycles );
|
||
MCintsqr2 = MCintsqr2/((double) MonteCarloCycles );
|
||
double STDev = sqrt(MCintsqr2-Mean*Mean);
|
||
double Variance = MCintsqr2-Mean*Mean;
|
||
// Write mean value and variance
|
||
cout << " Sample variance= " << Variance << " Mean value = " << Mean << endl;
|
||
// Now we compute the autocorrelation function
|
||
vec autocorrelation = zeros<vec>(MonteCarloCycles);
|
||
for (int j = 0; j < MonteCarloCycles; j++){
|
||
double sum = 0.0;
|
||
for (int k = 0; k < (MonteCarloCycles-j); k++){
|
||
sum += (X(k)-Mean)*(X(k+j)-Mean);
|
||
}
|
||
autocorrelation(j) = sum/Variance/((double) MonteCarloCycles );
|
||
ofile << setiosflags(ios::showpoint | ios::uppercase);
|
||
ofile << setw(15) << setprecision(8) << j;
|
||
ofile << setw(15) << setprecision(8) << autocorrelation(j) << endl;
|
||
}
|
||
// Now compute the exact covariance using the autocorrelation function
|
||
double Covariance = 0.0;
|
||
for (int j = 0; j < MonteCarloCycles; j++){
|
||
Covariance += autocorrelation(j);
|
||
}
|
||
Covariance *= 2.0/((double) MonteCarloCycles);
|
||
// Compute now the total variance, including the covariance, and obtain the standard deviation
|
||
double TotalVariance = (Variance/((double) MonteCarloCycles ))+Covariance;
|
||
cout << "Covariance =" << Covariance << "Totalvariance= " << TotalVariance << "Sample Variance/n= " << (Variance/((double) MonteCarloCycles )) << endl;
|
||
cout << " STD from sample variance= " << sqrt(Variance/((double) MonteCarloCycles )) << " STD with covariance = " << sqrt(TotalVariance) << endl;
|
||
|
||
ofile.close(); // close output file
|
||
return 0;
|
||
} // end of main program
|
||
</pre></div>
|
||
</div>
|
||
</div>
|
||
<div class="section" id="what-is-blocking">
|
||
<h2><span class="section-number">1.65. </span>What is blocking?<a class="headerlink" href="#what-is-blocking" title="Permalink to this headline">¶</a></h2>
|
||
<p><strong>Blocking.</strong></p>
|
||
<div class="highlight-none notranslate"><div class="highlight"><pre><span></span>* Say that we have a set of samples from a Monte Carlo experiment
|
||
|
||
* Assuming (wrongly) that our samples are uncorrelated our best estimate of the standard deviation of the mean $\langle \mathbf{M}\rangle$ is given by
|
||
</pre></div>
|
||
</div>
|
||
<div class="math notranslate nohighlight">
|
||
\[
|
||
\sigma=\sqrt{\frac{1}{n}\left(\langle \mathbf{M}^2\rangle-\langle \mathbf{M}\rangle^2\right)}
|
||
\]</div>
|
||
<ul class="simple">
|
||
<li><p>If the samples are correlated we can rewrite our results to show that</p></li>
|
||
</ul>
|
||
<div class="math notranslate nohighlight">
|
||
\[
|
||
\sigma=\sqrt{\frac{1+2\tau/\Delta t}{n}\left(\langle \mathbf{M}^2\rangle-\langle \mathbf{M}\rangle^2\right)}
|
||
\]</div>
|
||
<p>where <span class="math notranslate nohighlight">\(\tau\)</span> is the correlation time (the time between a sample and the next uncorrelated sample) and <span class="math notranslate nohighlight">\(\Delta t\)</span> is time between each sample</p>
|
||
</div>
|
||
<div class="section" id="id17">
|
||
<h2><span class="section-number">1.66. </span>What is blocking?<a class="headerlink" href="#id17" title="Permalink to this headline">¶</a></h2>
|
||
<p><strong>Blocking.</strong></p>
|
||
<div class="highlight-none notranslate"><div class="highlight"><pre><span></span>* If $\Delta t\gg\tau$ our first estimate of $\sigma$ still holds
|
||
|
||
* Much more common that $\Delta t<\tau$
|
||
|
||
* In the method of data blocking we divide the sequence of samples into blocks
|
||
|
||
* We then take the mean $\langle \mathbf{M}_i\rangle$ of block $i=1\ldots n_{blocks}$ to calculate the total mean and variance
|
||
|
||
* The size of each block must be so large that sample $j$ of block $i$ is not correlated with sample $j$ of block $i+1$
|
||
|
||
* The correlation time $\tau$ would be a good choice
|
||
</pre></div>
|
||
</div>
|
||
</div>
|
||
<div class="section" id="id18">
|
||
<h2><span class="section-number">1.67. </span>What is blocking?<a class="headerlink" href="#id18" title="Permalink to this headline">¶</a></h2>
|
||
<p><strong>Blocking.</strong></p>
|
||
<div class="highlight-none notranslate"><div class="highlight"><pre><span></span>* Problem: We don't know $\tau$ or it is too expensive to compute
|
||
|
||
* Solution: Make a plot of std. dev. as a function of blocksize
|
||
|
||
* The estimate of std. dev. of correlated data is too low $\to$ the error will increase with increasing block size until the blocks are uncorrelated, where we reach a plateau
|
||
|
||
* When the std. dev. stops increasing the blocks are uncorrelated
|
||
</pre></div>
|
||
</div>
|
||
</div>
|
||
<div class="section" id="implementation">
|
||
<h2><span class="section-number">1.68. </span>Implementation<a class="headerlink" href="#implementation" title="Permalink to this headline">¶</a></h2>
|
||
<div class="highlight-none notranslate"><div class="highlight"><pre><span></span>* Do a Monte Carlo simulation, storing all samples to file
|
||
|
||
* Do the statistical analysis on this file, independently of your Monte Carlo program
|
||
|
||
* Read the file into an array
|
||
|
||
* Loop over various block sizes
|
||
|
||
* For each block size $n_b$, loop over the array in steps of $n_b$ taking the mean of elements $i n_b,\ldots,(i+1) n_b$
|
||
|
||
* Take the mean and variance of the resulting array
|
||
|
||
* Write the results for each block size to file for later
|
||
analysis
|
||
</pre></div>
|
||
</div>
|
||
</div>
|
||
<div class="section" id="actual-implementation-with-code-main-function">
|
||
<h2><span class="section-number">1.69. </span>Actual implementation with code, main function<a class="headerlink" href="#actual-implementation-with-code-main-function" title="Permalink to this headline">¶</a></h2>
|
||
<p>When the file gets large, it can be useful to write your data in binary mode instead of ascii characters.
|
||
The <a class="reference external" href="https://github.com/CompPhysics/MachineLearning/blob/master/doc/Programs/Sampling/analysis.py">following python file</a> reads data from file with the output from every Monte Carlo cycle.</p>
|
||
<div class="cell docutils container">
|
||
<div class="cell_input docutils container">
|
||
<div class="highlight-ipython3 notranslate"><div class="highlight"><pre><span></span><span class="c1"># Blocking</span>
|
||
<span class="nd">@timeFunction</span>
|
||
<span class="k">def</span> <span class="nf">blocking</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span> <span class="n">blockSizeMax</span> <span class="o">=</span> <span class="mi">500</span><span class="p">):</span>
|
||
<span class="n">blockSizeMin</span> <span class="o">=</span> <span class="mi">1</span>
|
||
|
||
<span class="bp">self</span><span class="o">.</span><span class="n">blockSizes</span> <span class="o">=</span> <span class="p">[]</span>
|
||
<span class="bp">self</span><span class="o">.</span><span class="n">meanVec</span> <span class="o">=</span> <span class="p">[]</span>
|
||
<span class="bp">self</span><span class="o">.</span><span class="n">varVec</span> <span class="o">=</span> <span class="p">[]</span>
|
||
|
||
<span class="k">for</span> <span class="n">i</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="n">blockSizeMin</span><span class="p">,</span> <span class="n">blockSizeMax</span><span class="p">):</span>
|
||
<span class="k">if</span><span class="p">(</span><span class="nb">len</span><span class="p">(</span><span class="bp">self</span><span class="o">.</span><span class="n">data</span><span class="p">)</span> <span class="o">%</span> <span class="n">i</span> <span class="o">!=</span> <span class="mi">0</span><span class="p">):</span>
|
||
<span class="k">pass</span><span class="c1">#continue</span>
|
||
<span class="n">blockSize</span> <span class="o">=</span> <span class="n">i</span>
|
||
<span class="n">meanTempVec</span> <span class="o">=</span> <span class="p">[]</span>
|
||
<span class="n">varTempVec</span> <span class="o">=</span> <span class="p">[]</span>
|
||
<span class="n">startPoint</span> <span class="o">=</span> <span class="mi">0</span>
|
||
<span class="n">endPoint</span> <span class="o">=</span> <span class="n">blockSize</span>
|
||
|
||
<span class="k">while</span> <span class="n">endPoint</span> <span class="o"><=</span> <span class="nb">len</span><span class="p">(</span><span class="bp">self</span><span class="o">.</span><span class="n">data</span><span class="p">):</span>
|
||
<span class="n">meanTempVec</span><span class="o">.</span><span class="n">append</span><span class="p">(</span><span class="n">np</span><span class="o">.</span><span class="n">average</span><span class="p">(</span><span class="bp">self</span><span class="o">.</span><span class="n">data</span><span class="p">[</span><span class="n">startPoint</span><span class="p">:</span><span class="n">endPoint</span><span class="p">]))</span>
|
||
<span class="n">startPoint</span> <span class="o">=</span> <span class="n">endPoint</span>
|
||
<span class="n">endPoint</span> <span class="o">+=</span> <span class="n">blockSize</span>
|
||
<span class="n">mean</span><span class="p">,</span> <span class="n">var</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">average</span><span class="p">(</span><span class="n">meanTempVec</span><span class="p">),</span> <span class="n">np</span><span class="o">.</span><span class="n">var</span><span class="p">(</span><span class="n">meanTempVec</span><span class="p">)</span><span class="o">/</span><span class="nb">len</span><span class="p">(</span><span class="n">meanTempVec</span><span class="p">)</span>
|
||
<span class="bp">self</span><span class="o">.</span><span class="n">meanVec</span><span class="o">.</span><span class="n">append</span><span class="p">(</span><span class="n">mean</span><span class="p">)</span>
|
||
<span class="bp">self</span><span class="o">.</span><span class="n">varVec</span><span class="o">.</span><span class="n">append</span><span class="p">(</span><span class="n">var</span><span class="p">)</span>
|
||
<span class="bp">self</span><span class="o">.</span><span class="n">blockSizes</span><span class="o">.</span><span class="n">append</span><span class="p">(</span><span class="n">blockSize</span><span class="p">)</span>
|
||
|
||
<span class="bp">self</span><span class="o">.</span><span class="n">blockingAvg</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">average</span><span class="p">(</span><span class="bp">self</span><span class="o">.</span><span class="n">meanVec</span><span class="p">[</span><span class="o">-</span><span class="mi">200</span><span class="p">:])</span>
|
||
<span class="bp">self</span><span class="o">.</span><span class="n">blockingVar</span> <span class="o">=</span> <span class="p">(</span><span class="n">np</span><span class="o">.</span><span class="n">average</span><span class="p">(</span><span class="bp">self</span><span class="o">.</span><span class="n">varVec</span><span class="p">[</span><span class="o">-</span><span class="mi">200</span><span class="p">:]))</span>
|
||
<span class="bp">self</span><span class="o">.</span><span class="n">blockingStd</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">sqrt</span><span class="p">(</span><span class="bp">self</span><span class="o">.</span><span class="n">blockingVar</span><span class="p">)</span>
|
||
</pre></div>
|
||
</div>
|
||
</div>
|
||
<div class="cell_output docutils container">
|
||
<div class="output traceback highlight-ipythontb notranslate"><div class="highlight"><pre><span></span><span class="gt"> File</span><span class="nn"> "<ipython-input-6-2ff97f4bf03b>"</span><span class="gt">, line </span><span class="mi">2</span>
|
||
<span class="nd">@timeFunction</span>
|
||
<span class="o">^</span>
|
||
<span class="ne">IndentationError</span>: unexpected indent
|
||
</pre></div>
|
||
</div>
|
||
</div>
|
||
</div>
|
||
</div>
|
||
<div class="section" id="the-bootstrap-method">
|
||
<h2><span class="section-number">1.70. </span>The Bootstrap method<a class="headerlink" href="#the-bootstrap-method" title="Permalink to this headline">¶</a></h2>
|
||
<p>The Bootstrap resampling method is also very popular. It is very simple:</p>
|
||
<ol class="simple">
|
||
<li><p>Start with your sample of measurements and compute the sample variance and the mean values</p></li>
|
||
<li><p>Then start again but pick in a random way the numbers in the sample and recalculate the mean and the sample variance.</p></li>
|
||
<li><p>Repeat this <span class="math notranslate nohighlight">\(K\)</span> times.</p></li>
|
||
</ol>
|
||
<p>It can be shown, see the article by <a class="reference external" href="https://projecteuclid.org/download/pdf_1/euclid.aos/1176344552">Efron</a>
|
||
that it produces the correct standard deviation.</p>
|
||
<p>This method is very useful for small ensembles of data points.</p>
|
||
</div>
|
||
<div class="section" id="bootstrapping">
|
||
<h2><span class="section-number">1.71. </span>Bootstrapping<a class="headerlink" href="#bootstrapping" title="Permalink to this headline">¶</a></h2>
|
||
<p>Given a set of <span class="math notranslate nohighlight">\(N\)</span> data, assume that we are interested in some
|
||
observable <span class="math notranslate nohighlight">\(\theta\)</span> which may be estimated from that set. This observable can also be for example the result of a fit based on all <span class="math notranslate nohighlight">\(N\)</span> raw data.
|
||
Let us call the value of the observable obtained from the original
|
||
data set <span class="math notranslate nohighlight">\(\hat{\theta}\)</span>. One recreates from the sample repeatedly
|
||
other samples by choosing randomly <span class="math notranslate nohighlight">\(N\)</span> data out of the original set.
|
||
This costs essentially nothing, since we just recycle the original data set for the building of new sets.</p>
|
||
</div>
|
||
<div class="section" id="bootstrapping-recipe">
|
||
<h2><span class="section-number">1.72. </span>Bootstrapping, recipe<a class="headerlink" href="#bootstrapping-recipe" title="Permalink to this headline">¶</a></h2>
|
||
<p>Let us assume we have done this <span class="math notranslate nohighlight">\(K\)</span> times and thus have <span class="math notranslate nohighlight">\(K\)</span> sets of <span class="math notranslate nohighlight">\(N\)</span>
|
||
data values each.
|
||
Of course some values will enter more than once in the new sets. For each of these sets one computes the observable <span class="math notranslate nohighlight">\(\theta\)</span> resulting in values <span class="math notranslate nohighlight">\(\theta_k\)</span> with <span class="math notranslate nohighlight">\(k = 1,...,K\)</span>. Then one determines</p>
|
||
<div class="math notranslate nohighlight">
|
||
\[
|
||
\tilde{\theta} = \frac{1}{K} \sum_{k=1}^K \theta_k,
|
||
\]</div>
|
||
<p>and</p>
|
||
<div class="math notranslate nohighlight">
|
||
\[
|
||
sigma^2_{\tilde{\theta}} = \frac{1}{K} \sum_{k=1}^K \left(\theta_k-\tilde{\theta}\right)^2.
|
||
\]</div>
|
||
<p>These are estimators for <span class="math notranslate nohighlight">\(\angle\theta\rangle\)</span> and its variance. They are not unbiased and therefore
|
||
<span class="math notranslate nohighlight">\(\tilde{\theta}\neq\hat{\theta}\)</span> for finite K.</p>
|
||
<p>The difference is called bias and gives an idea on how far away the result may be from
|
||
the true <span class="math notranslate nohighlight">\(\angle\theta\rangle\)</span>. As final result for the observable one quotes <span class="math notranslate nohighlight">\(\angle\theta\rangle = \tilde{\theta} \pm \sigma_{\tilde{\theta}}\)</span> .</p>
|
||
</div>
|
||
<div class="section" id="bootstrapping-code">
|
||
<h2><span class="section-number">1.73. </span>Bootstrapping, <a class="reference external" href="https://github.com/CompPhysics/MachineLearning/blob/master/doc/Programs/Sampling/analysis.py">code</a><a class="headerlink" href="#bootstrapping-code" title="Permalink to this headline">¶</a></h2>
|
||
<div class="highlight-none notranslate"><div class="highlight"><pre><span></span> # Bootstrap
|
||
@timeFunction
|
||
def bootstrap(self, nBoots = 1000):
|
||
bootVec = np.zeros(nBoots)
|
||
for k in range(0,nBoots):
|
||
bootVec[k] = np.average(np.random.choice(self.data, len(self.data)))
|
||
self.bootAvg = np.average(bootVec)
|
||
self.bootVar = np.var(bootVec)
|
||
self.bootStd = np.std(bootVec)
|
||
</pre></div>
|
||
</div>
|
||
</div>
|
||
<div class="section" id="jackknife-code">
|
||
<h2><span class="section-number">1.74. </span>Jackknife, <a class="reference external" href="https://github.com/CompPhysics/MachineLearning/blob/master/doc/Programs/Sampling/analysis.py">code</a><a class="headerlink" href="#jackknife-code" title="Permalink to this headline">¶</a></h2>
|
||
<div class="highlight-none notranslate"><div class="highlight"><pre><span></span> # Jackknife
|
||
@timeFunction
|
||
def jackknife(self):
|
||
jackknVec = np.zeros(len(self.data))
|
||
for k in range(0,len(self.data)):
|
||
jackknVec[k] = np.average(np.delete(self.data, k))
|
||
self.jackknAvg = self.avg - (len(self.data) - 1) * (np.average(jackknVec) - self.avg)
|
||
self.jackknVar = float(len(self.data) - 1) * np.var(jackknVec)
|
||
self.jackknStd = np.sqrt(self.jackknVar)
|
||
</pre></div>
|
||
</div>
|
||
</div>
|
||
</div>
|
||
|
||
<script type="text/x-thebe-config">
|
||
{
|
||
requestKernel: true,
|
||
binderOptions: {
|
||
repo: "binder-examples/jupyter-stacks-datascience",
|
||
ref: "master",
|
||
},
|
||
codeMirrorConfig: {
|
||
theme: "abcdef",
|
||
mode: "python"
|
||
},
|
||
kernelOptions: {
|
||
kernelName: "python3",
|
||
path: "./."
|
||
},
|
||
predefinedOutput: true
|
||
}
|
||
</script>
|
||
<script>kernelName = 'python3'</script>
|
||
|
||
</div>
|
||
|
||
|
||
<div class='prev-next-bottom'>
|
||
|
||
<a class='left-prev' id="prev-link" href="textbooks.html" title="previous page">Textbooks</a>
|
||
<a class='right-next' id="next-link" href="linalg.html" title="next page"><span class="section-number">2. </span>Linear Algebra, Handling of Arrays and more Python Features</a>
|
||
|
||
</div>
|
||
|
||
</div>
|
||
</div>
|
||
<footer class="footer mt-5 mt-md-0">
|
||
<div class="container">
|
||
<p>
|
||
|
||
By Morten Hjorth-Jensen<br/>
|
||
|
||
© Copyright 2021.<br/>
|
||
</p>
|
||
</div>
|
||
</footer>
|
||
</main>
|
||
|
||
|
||
</div>
|
||
</div>
|
||
|
||
<script src="_static/js/index.1c5a1a01449ed65a7b51.js"></script>
|
||
|
||
|
||
</body>
|
||
</html> |