tag:blogger.com,1999:blog-75669454907707988402024-03-23T03:16:04.388-07:00Get Your Data OnData Science, Signal Processing, Technology, LifeChris Carbonehttp://www.blogger.com/profile/10458001230470047172noreply@blogger.comBlogger30125tag:blogger.com,1999:blog-7566945490770798840.post-29558331849430865252023-11-18T20:14:00.000-08:002023-11-18T20:14:06.233-08:00November Fog<div class="separator" style="clear: both;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhdiK7_B2uR-6q6sJHKWqYF6OiK92WnDA0RNe7yTtYWVn1ZGTE9YKUxI9eHSGcenoob-hymBRBq5vDAtFUgIpanZQQCwWclMuPOurw27-UN_vXfskkqiwXvaebncoPJVb0AfX3xpAa0S51hMdg6wBvOltLtDrq706S4S-ORSF4pXqA8XKoVldBLJJZNDN08/s7702/NovemberFog.jpeg" style="display: block; padding: 1em 0; text-align: center; "><img alt="" border="0" width="400" data-original-height="4546" data-original-width="7702" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhdiK7_B2uR-6q6sJHKWqYF6OiK92WnDA0RNe7yTtYWVn1ZGTE9YKUxI9eHSGcenoob-hymBRBq5vDAtFUgIpanZQQCwWclMuPOurw27-UN_vXfskkqiwXvaebncoPJVb0AfX3xpAa0S51hMdg6wBvOltLtDrq706S4S-ORSF4pXqA8XKoVldBLJJZNDN08/s400/NovemberFog.jpeg"/></a></div>Chris Carbonehttp://www.blogger.com/profile/10458001230470047172noreply@blogger.com0tag:blogger.com,1999:blog-7566945490770798840.post-46104494410950030972023-10-29T20:00:00.001-07:002023-10-29T20:00:04.326-07:00Pumpkins! <div class="separator" style="clear: both;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEiYVWNd2imsMNwt1sA3tBxdpLzjuQT4TpterAaOVHcpG5L-SEWKCcZ_-wR5RwEC3uD9Q5elHeQZUsz3hN7Ck3sIBA6pNQ77_X8CNuRn61JTnNdESxD7ZAiyuwVKIGXMBxCHNzf03nPmWZUds0ZJ4pw_08z96Gc0ipfeioGbuKphIG-qlGHJxX-3aDImpqMY/s1200/Pumpkin%20Night.jpg" style="display: block; padding: 1em 0; text-align: center; "><img alt="" border="0" width="400" data-original-height="867" data-original-width="1200" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEiYVWNd2imsMNwt1sA3tBxdpLzjuQT4TpterAaOVHcpG5L-SEWKCcZ_-wR5RwEC3uD9Q5elHeQZUsz3hN7Ck3sIBA6pNQ77_X8CNuRn61JTnNdESxD7ZAiyuwVKIGXMBxCHNzf03nPmWZUds0ZJ4pw_08z96Gc0ipfeioGbuKphIG-qlGHJxX-3aDImpqMY/s400/Pumpkin%20Night.jpg"/></a></div>Chris Carbonehttp://www.blogger.com/profile/10458001230470047172noreply@blogger.com0tag:blogger.com,1999:blog-7566945490770798840.post-45952141091928456452023-07-26T18:48:00.002-07:002023-07-26T19:13:33.579-07:00Wednesday Linkorama!<!DOCTYPE html>
<html xmlns="http://www.w3.org/1999/xhtml" lang="en" xml:lang="en">
<head>
<meta charset="utf-8">
<meta name="generator" content="quarto-1.1.189">
<meta name="viewport" content="width=device-width, initial-scale=1.0, user-scalable=yes">
<style>
code{white-space: pre-wrap;}
span.smallcaps{font-variant: small-caps;}
div.columns{display: flex; gap: min(4vw, 1.5em);}
div.column{flex: auto; overflow-x: auto;}
div.hanging-indent{margin-left: 1.5em; text-indent: -1.5em;}
ul.task-list{list-style: none;}
ul.task-list li input[type="checkbox"] {
width: 0.8em;
margin: 0 0.8em 0.2em -1.6em;
vertical-align: middle;
}
.display.math{display: block; text-align: center; margin: 0.5rem auto;}
</style>
<link href="linkarama_files/libs/quarto-html/quarto-html.min.css" rel="stylesheet" data-mode="light">
<link href="linkarama_files/libs/quarto-html/quarto-syntax-highlighting.css" rel="stylesheet" id="quarto-text-highlighting-styles">
</head>
<body>
<header id="title-block-header">
</header>
<h4 id="round-up-of-things-i-found-interesting.">Round up of things I found interesting.</h4>
<!-- The link was from the story page. The w_1315 was changed to w_400 so the figure would fit nicly -->
<p><img src="https://i.kinja-img.com/gawker-media/image/upload/c_fit,f_auto,g_center,q_60,w_400/eac93f71fadfad3b7570b788c64f9c8e.jpg"></p>
<p>The image above came from a <a href="https://gizmodo.com/measuring-helium-in-distant-galaxies-may-give-physicist-1850679794">Gizmodo</a> article explaining how measuring helium from distant galaxies offers clues to the matter/antimatter imbalance. </p>
<p>The four-color problem asks, "Are four colors sufficient to color the countries of a map so that no two adjacent countries share the same color?" A few minutes of trying this with some paper and markers would lead you to think that four is enough. But it was hard to rigorously prove that 4 is the minimum. In 1976, with the help of a computer program, the question was answered. This article talks about the <a href="https://www.scientificamerican.com/article/how-a-doodlers-problem-sparked-a-controversy-in-math/"> problem</a> and the controversy surrounding the use of computer programs in math proofs.</p>
<p>The news about the room temperature <a href="https://www.newscientist.com/article/2384782-room-temperature-superconductor-breakthrough-met-with-scepticism/"> superconductor</a> would be the biggest science and technology news of the week, but there is quite a bit of skepticism about the result. If true, it would be a huge deal. It could speed up computers, improve batteries, bigger electro magnetic (floating trains!) and so on.
</body></html>Chris Carbonehttp://www.blogger.com/profile/10458001230470047172noreply@blogger.com0tag:blogger.com,1999:blog-7566945490770798840.post-6273540240295817702023-07-25T17:33:00.005-07:002023-07-27T16:23:57.892-07:00Autoregressive Moving Average Models and Power Spectral Densities<!DOCTYPE html>
<html xmlns="http://www.w3.org/1999/xhtml" lang="en" xml:lang="en"><head>
<meta charset="utf-8">
<meta name="generator" content="quarto-1.2.280">
<meta name="viewport" content="width=device-width, initial-scale=1.0, user-scalable=yes">
<style>
code{white-space: pre-wrap;}
span.smallcaps{font-variant: small-caps;}
div.columns{display: flex; gap: min(4vw, 1.5em);}
div.column{flex: auto; overflow-x: auto;}
div.hanging-indent{margin-left: 1.5em; text-indent: -1.5em;}
ul.task-list{list-style: none;}
ul.task-list li input[type="checkbox"] {
width: 0.8em;
margin: 0 0.8em 0.2em -1.6em;
vertical-align: middle;
}
pre > code.sourceCode { white-space: pre; position: relative; }
pre > code.sourceCode > span { display: inline-block; line-height: 1.25; }
pre > code.sourceCode > span:empty { height: 1.2em; }
.sourceCode { overflow: visible; }
code.sourceCode > span { color: inherit; text-decoration: inherit; }
div.sourceCode { margin: 1em 0; }
pre.sourceCode { margin: 0; }
@media screen {
div.sourceCode { overflow: auto; }
}
@media print {
pre > code.sourceCode { white-space: pre-wrap; }
pre > code.sourceCode > span { text-indent: -5em; padding-left: 5em; }
}
pre.numberSource code
{ counter-reset: source-line 0; }
pre.numberSource code > span
{ position: relative; left: -4em; counter-increment: source-line; }
pre.numberSource code > span > a:first-child::before
{ content: counter(source-line);
position: relative; left: -1em; text-align: right; vertical-align: baseline;
border: none; display: inline-block;
-webkit-touch-callout: none; -webkit-user-select: none;
-khtml-user-select: none; -moz-user-select: none;
-ms-user-select: none; user-select: none;
padding: 0 4px; width: 4em;
color: #aaaaaa;
}
pre.numberSource { margin-left: 3em; border-left: 1px solid #aaaaaa; padding-left: 4px; }
div.sourceCode
{ }
@media screen {
pre > code.sourceCode > span > a:first-child::before { text-decoration: underline; }
}
code span.al { color: #ff0000; font-weight: bold; } /* Alert */
code span.an { color: #60a0b0; font-weight: bold; font-style: italic; } /* Annotation */
code span.at { color: #7d9029; } /* Attribute */
code span.bn { color: #40a070; } /* BaseN */
code span.bu { color: #008000; } /* BuiltIn */
code span.cf { color: #007020; font-weight: bold; } /* ControlFlow */
code span.ch { color: #4070a0; } /* Char */
code span.cn { color: #880000; } /* Constant */
code span.co { color: #60a0b0; font-style: italic; } /* Comment */
code span.cv { color: #60a0b0; font-weight: bold; font-style: italic; } /* CommentVar */
code span.do { color: #ba2121; font-style: italic; } /* Documentation */
code span.dt { color: #902000; } /* DataType */
code span.dv { color: #40a070; } /* DecVal */
code span.er { color: #ff0000; font-weight: bold; } /* Error */
code span.ex { } /* Extension */
code span.fl { color: #40a070; } /* Float */
code span.fu { color: #06287e; } /* Function */
code span.im { color: #008000; font-weight: bold; } /* Import */
code span.in { color: #60a0b0; font-weight: bold; font-style: italic; } /* Information */
code span.kw { color: #007020; font-weight: bold; } /* Keyword */
code span.op { color: #666666; } /* Operator */
code span.ot { color: #007020; } /* Other */
code span.pp { color: #bc7a00; } /* Preprocessor */
code span.sc { color: #4070a0; } /* SpecialChar */
code span.ss { color: #bb6688; } /* SpecialString */
code span.st { color: #4070a0; } /* String */
code span.va { color: #19177c; } /* Variable */
code span.vs { color: #4070a0; } /* VerbatimString */
code span.wa { color: #60a0b0; font-weight: bold; font-style: italic; } /* Warning */
div.csl-bib-body { }
div.csl-entry {
clear: both;
}
.hanging div.csl-entry {
margin-left:2em;
text-indent:-2em;
}
div.csl-left-margin {
min-width:2em;
float:left;
}
div.csl-right-inline {
margin-left:2em;
padding-left:1em;
}
div.csl-indent {
margin-left: 2em;
}
</style>
<script src="ArmaEstimation_files/libs/clipboard/clipboard.min.js"></script>
<script src="ArmaEstimation_files/libs/quarto-html/quarto.js"></script>
<script src="ArmaEstimation_files/libs/quarto-html/popper.min.js"></script>
<script src="ArmaEstimation_files/libs/quarto-html/tippy.umd.min.js"></script>
<script src="ArmaEstimation_files/libs/quarto-html/anchor.min.js"></script>
<link href="ArmaEstimation_files/libs/quarto-html/tippy.css" rel="stylesheet">
<link href="ArmaEstimation_files/libs/quarto-html/quarto-syntax-highlighting.css" rel="stylesheet" id="quarto-text-highlighting-styles">
<script src="ArmaEstimation_files/libs/bootstrap/bootstrap.min.js"></script>
<link href="ArmaEstimation_files/libs/bootstrap/bootstrap-icons.css" rel="stylesheet">
<link href="ArmaEstimation_files/libs/bootstrap/bootstrap.min.css" rel="stylesheet" id="quarto-bootstrap" data-mode="light">
<script src="https://cdn.jsdelivr.net/npm/mathjax@3/es5/tex-chtml-full.js" type="text/javascript"></script>
</head>
<body class="fullcontent">
<div id="quarto-content" class="page-columns page-rows-contents page-layout-article">
<main class="content" id="quarto-document-content">
<header id="title-block-header" class="quarto-title-block default">
<div class="quarto-title">
</div>
<div class="quarto-title-meta">
</div>
</header>
<section id="sec-introduction" class="level2" data-number="1">
<h2 data-number="1" class="anchored" data-anchor-id="sec-introduction"><span class="header-section-number">1</span> Introduction</h2>
<p>The simplest way to estimate a power spectral density (PSD) is to use a periodogram, defined as</p>
<div class="cell">
<div class="sourceCode cell-code" id="cb1"><pre class="sourceCode r code-with-copy"><code class="sourceCode r"><span id="cb1-1"><a href="#cb1-1" aria-hidden="true" tabindex="-1"></a>periodogram <span class="ot"><-</span> <span class="cf">function</span>(x) (<span class="dv">1</span><span class="sc">/</span><span class="fu">length</span>(x))<span class="sc">*</span><span class="fu">abs</span>(<span class="fu">fft</span>(x))<span class="sc">^</span><span class="dv">2</span></span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
</div>
<p>where <code>x</code> is the data samples. The periodogram estimates the power level at <span class="math inline">\(N\)</span> frequency locations where <span class="math inline">\(N\)</span> is the number of samples in x. The periodogram uses <span class="math inline">\(N\)</span> data samples to estimate <span class="math inline">\(N\)</span> power levels. This leads to high variance estimates. There are ways to reduce the variance, but they also reduce the resolution.</p>
<p>A different approach is to estimate the PSD using a model. A good model can capture the important features of the PSD and be controlled by a few parameters. The autoregressive moving average (ARMA) model and the closely related autoregressive (AR) and moving average (MA) models can model a wide variety of PSDs and can be controlled by a few parameters.</p>
<p>This article defines the ARMA, AR, and MA models and then shows how to calculate their PSDs. We then cover how they are used to estimate PSDs, how to select between the three models, and how to choose the number of parameters (the model order). By the end of this article, you will understand how ARMA/AR/MA models are used for PSD estimation.</p>
<p>Much of this article comes from Chapter 8 of <a href="https://www.amazon.com/Power-Spectral-Density-Estimation-Example/dp/B0BHWDC4VW">Power Spectral Density Estimation by Example</a> <span class="citation" data-cites="Carbone2022">Carbone (<a href="#ref-Carbone2022" role="doc-biblioref">2022</a>)</span>.</p>
</section>
<section id="sec-armaModel" class="level2" data-number="2">
<h2 data-number="2" class="anchored" data-anchor-id="sec-armaModel"><span class="header-section-number">2</span> ARMA Model</h2>
<p>An ARMA model of order <span class="math inline">\(p\)</span>, <span class="math inline">\(q\)</span> is defined as</p>
<p><span id="eq-armaDef"><span class="math display">\[
a[0]y[n] = -\sum_{k=1}^p a[k]y[n-k] +\sum_{k=0}^q b[k]x[n-k],
\tag{1}\]</span></span></p>
<p>where <span class="math inline">\(y[n]\)</span> is the current output, <span class="math inline">\(a[k]\)</span> and <span class="math inline">\(b[k]\)</span> are the coefficients of the model, <span class="math inline">\(p\)</span> is the number of <span class="math inline">\(a\)</span> coefficients, <span class="math inline">\(q\)</span> is number of <span class="math inline">\(b\)</span> coefficients, <span class="math inline">\(y[n − k]\)</span> are the past outputs, and <span class="math inline">\(x[n − k]\)</span> are past (and present) inputs. The coefficient <span class="math inline">\(a[0]\)</span> is always assumed to be 1, so from here on; we replace <span class="math inline">\(a[0]\)</span> with 1. If <span class="math inline">\(a[k]\)</span> <span class="math inline">\(\neq\)</span> 1 we could divide <a href="#eq-armaDef">Equation 1</a> by <span class="math inline">\(a[0]\)</span>.</p>
<p>The input <span class="math inline">\(x\)</span> is a sequence of samples from a white noise process with variance <span class="math inline">\(\sigma^2\)</span>. The input process could be any white noise process, but for the sake of specificity, we assume the process is white Gaussian noise (WGN). The ARMA model of order <span class="math inline">\(p\)</span>, <span class="math inline">\(q\)</span> is written as ARMA[<span class="math inline">\(p,q\)</span>]. The parameters of the ARMA model are the autoregressive coefficients <span class="math inline">\(a[k]\)</span>, the moving average coefficients <span class="math inline">\(b[k]\)</span>, and the variance of the input <span class="math inline">\(\sigma^2\)</span>.</p>
<p>Below, we use the ARMA[1,1] model to produce simulated data. In this case, we let <span class="math inline">\(a[1] = 0.9\)</span>, <span class="math inline">\(b[1] = −0.9\)</span>, and the input <span class="math inline">\(x[n]\)</span> is WGN with variance 1. The model specified above is written as <span id="eq-arma11Example"><span class="math display">\[
y[n] = −0.9y[n − 1] − 0.9x[n − 1] + x[n].
\tag{2}\]</span></span></p>
<p>Note the <span class="math inline">\(a[1]\)</span> coefficient is 0.9. The negative sign in front of 0.9 is part of <a href="#eq-armaDef">Equation 1</a>. Below we use <code>rnorm(N)</code> to produce <code>N</code> samples of WGN with variance 1.</p>
<div class="cell">
<div class="sourceCode cell-code" id="cb2"><pre class="sourceCode r code-with-copy"><code class="sourceCode r"><span id="cb2-1"><a href="#cb2-1" aria-hidden="true" tabindex="-1"></a><span class="fu">set.seed</span>(<span class="dv">1</span>) </span>
<span id="cb2-2"><a href="#cb2-2" aria-hidden="true" tabindex="-1"></a>N <span class="ot"><-</span> <span class="dv">512</span> </span>
<span id="cb2-3"><a href="#cb2-3" aria-hidden="true" tabindex="-1"></a>x <span class="ot"><-</span> <span class="fu">rnorm</span>(N)</span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
</div>
<p><a href="#fig-inputToArmaProcess">Figure 1</a> is a plot of <code>x</code>, the input into the ARMA[1,1] process defined in <a href="#eq-arma11Example">Equation 2</a>.</p>
<!--# fig-inputToArmaProcess XXX -->
<div class="cell">
<div class="cell-output-display">
<div id="fig-inputToArmaProcess" class="quarto-figure quarto-figure-center anchored">
<figure class="figure">
<p><img src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhWKV7eUUhUlES5L1Sry9XBBr7Uhg5NJNf6s8FZcL3jm4IQzG8XdD3rKzcDU1h71OWGxC5NWayZ3yW-3xMRff6e5TDB0mt6z8qnYGGydWN8Onk2OE9_1XBW23SJMUmn1_Vm0erz84mL4F8-weWYUxqphF8ZGomFgqHg5n2NTZ3ysM7s36baiGj4E79f9Jof/s1600/fig-inputToArmaProcess-1.png" class="img-fluid figure-img" width="768"></p>
<p></p><figcaption class="figure-caption">Figure 1: The figure is a plot of the <strong>input</strong> x to the ARMA[1,1] process defined in <a href="#eq-arma11Example">Equation 2</a></figcaption><p></p>
</figure>
</div>
</div>
</div>
<!--# fig-inputToArmaProcess -->
<p>Below is the ARMA[1,1] model from <a href="#eq-arma11Example">Equation 2</a> implemented in R.</p>
<div class="cell">
<div class="sourceCode cell-code" id="cb3"><pre class="sourceCode r code-with-copy"><code class="sourceCode r"><span id="cb3-1"><a href="#cb3-1" aria-hidden="true" tabindex="-1"></a>y_arma11 <span class="ot"><-</span> <span class="fu">rep</span>(<span class="dv">0</span>,N)</span>
<span id="cb3-2"><a href="#cb3-2" aria-hidden="true" tabindex="-1"></a>a <span class="ot"><-</span> <span class="fl">0.9</span> </span>
<span id="cb3-3"><a href="#cb3-3" aria-hidden="true" tabindex="-1"></a>b <span class="ot"><-</span> <span class="sc">-</span><span class="fl">0.9</span> </span>
<span id="cb3-4"><a href="#cb3-4" aria-hidden="true" tabindex="-1"></a></span>
<span id="cb3-5"><a href="#cb3-5" aria-hidden="true" tabindex="-1"></a><span class="cf">for</span> (n <span class="cf">in</span> <span class="dv">2</span><span class="sc">:</span>N) </span>
<span id="cb3-6"><a href="#cb3-6" aria-hidden="true" tabindex="-1"></a> y_arma11[n] <span class="ot"><-</span> <span class="sc">-</span>a[<span class="dv">1</span>]<span class="sc">*</span>y_arma11[n<span class="dv">-1</span>]<span class="sc">+</span> b[<span class="dv">1</span>]<span class="sc">*</span>x[n<span class="dv">-1</span>] <span class="sc">+</span> x[n] </span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
</div>
<p><a href="#fig-samplesArma11">Figure 2</a> is a plot of <code>y_arma11</code>, the model’s output.</p>
<div class="cell">
<div class="cell-output-display">
<div id="fig-samplesArma11" class="quarto-figure quarto-figure-center anchored">
<figure class="figure">
<p><img src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjNRfXeEqPk1dPcKa674IdakU8SDZW9a16QVA0tQa6pEn5tSsozMFCTqNKzGO8ZnbmtqKblZbry4ORymngimGc01tCdS2DDMFTK8MlVyBGQQMZ3_-WWr6vVGdarKWvO6nA7SBwMH4rf1iXao-5yz-BXcLOMuxwPq0vTVGVrcNHRodelrsKqsFH-xrSY_Uvj/s1600/fig-samplesArma11-1.png" class="img-fluid figure-img" width="768"></p>
<p></p><figcaption class="figure-caption">Figure 2: The figure is a plot of the <strong>output</strong> produced by the ARMA[1,1] process in <a href="#eq-arma11Example">Equation 2</a>.</figcaption><p></p>
</figure>
</div>
</div>
</div>
<p>The PSD <span class="math inline">\(P(f)\)</span> of an ARMA model is <span id="eq-armaPsdDef"><span class="math display">\[
P_{\text{ARMA}}(f) = \sigma^2 \left |\frac{ 1 + \sum_{k=1}^q b[k]\exp(-j2\pi fk) }
{ 1 + \sum_{k=1}^p a[k]\exp(-j2\pi fk) } \right |^2,
\tag{3}\]</span></span> where the <span class="math inline">\(a\)</span>’s and <span class="math inline">\(b\)</span>’s are the coefficients of the model and <span class="math inline">\(\sigma^2\)</span> is the variance of the WGN input. Notice the numerator and denominator in <a href="#eq-armaPsdDef">Equation 3</a> are polynomials. The roots (zeros) of these polynomials control the shape of the PSD. The denominator’s roots are called poles, and they control the peaks (high parts). The numerator’s roots are called zeros and control the valleys (low parts). Note that the peaks are more important than the valleys in most (not all!) applications.</p>
<p>Because ARMA models have poles and zeros, they are good at modeling PSDs with peaks and valleys. When used for PSD estimation, ARMA estimators are called “high resolution” estimators because they are better than the periodogram at resolving closely spaced peaks.</p>
<p>If we let <span class="math display">\[B(f) = 1 + \sum_{k=1}^q b[k]\exp(-j2\pi fk)\]</span> and <span class="math display">\[A(f) = 1 + \sum_{k=1}^p a[k]\exp(-j2\pi fk),\]</span> then we can write <a href="#eq-armaPsdDef">Equation 3</a> as <span id="eq-ArmaPsdCompact"><span class="math display">\[
P_{\text{ARMA}}(f) = \sigma^2 \left |\frac{B(f) }
{ A(f) } \right |^2 = \sigma^2 |B(f)|^2 \frac{1}{|A(f)|^2}.
\tag{4}\]</span></span> Notice <span class="math inline">\(A(f)\)</span> is the Fourier transform of the <span class="math inline">\(a\)</span>’s, and <span class="math inline">\(B(f)\)</span> is the Fourier transform of the <span class="math inline">\(b\)</span>’s. <a href="#eq-ArmaPsdCompact">Equation 4</a> makes clear that the ARMA PSD is the product of the input variance <span class="math inline">\(\sigma^2\)</span>, <span class="math inline">\(|B(f)|^2\)</span>, and <span class="math inline">\(\frac{1}{|A(f)|^2}\)</span>. Note the PSD of white noise with variance <span class="math inline">\(\sigma^2\)</span> is <span class="math inline">\(P_{\sigma^2}(f) = \sigma^2\)</span>, which is a constant.</p>
<p>The function <code>psdArma</code> inputs the ARMA parameters and output the PSD calculated at evenly spaced point along the frequency axis. The function inputs: the AR coefficients <code>arCoefficients</code>, the MA coefficients <code>maCoefficients</code>, the variance of the input process <code>inputVar</code>, and the number of points in the <code>fft</code> used to calculate the PSD. The function uses the fft function to find <span class="math inline">\(A(f)\)</span> and <span class="math inline">\(B(f)\)</span> in <a href="#eq-ArmaPsdCompact">Equation 4</a>. The function uses <code>zeroPad</code> (see below) to zero pad the AR and MA coefficients to produce <code>Nfs</code> frequency samples. The parameters must be known or estimated before using <code>psdArma</code>. If the parameters are known, the function outputs <em>the</em> PSD. If the parameters are estimated, the function outputs an <em>estimated</em> PSD.</p>
<div class="cell">
<div class="sourceCode cell-code" id="cb4"><pre class="sourceCode r code-with-copy"><code class="sourceCode r"><span id="cb4-1"><a href="#cb4-1" aria-hidden="true" tabindex="-1"></a>zeroPad <span class="ot"><-</span> <span class="cf">function</span>(x,N) {</span>
<span id="cb4-2"><a href="#cb4-2" aria-hidden="true" tabindex="-1"></a> <span class="fu">c</span>(x,<span class="fu">rep</span>(<span class="dv">0</span>,<span class="fu">max</span>(N<span class="sc">-</span><span class="fu">length</span>(x),<span class="dv">0</span>)))</span>
<span id="cb4-3"><a href="#cb4-3" aria-hidden="true" tabindex="-1"></a>}</span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
</div>
<div class="cell">
<div class="sourceCode cell-code" id="cb5"><pre class="sourceCode r code-with-copy"><code class="sourceCode r"><span id="cb5-1"><a href="#cb5-1" aria-hidden="true" tabindex="-1"></a>psdArma <span class="ot"><-</span> <span class="cf">function</span>(arCoefficients,maCoefficients,inputVar,<span class="at">Nfs=</span><span class="dv">256</span>){</span>
<span id="cb5-2"><a href="#cb5-2" aria-hidden="true" tabindex="-1"></a> A <span class="ot"><-</span> <span class="fu">fft</span>(<span class="fu">zeroPad</span>(arCoefficients,Nfs))</span>
<span id="cb5-3"><a href="#cb5-3" aria-hidden="true" tabindex="-1"></a> B <span class="ot"><-</span> <span class="fu">fft</span>(<span class="fu">zeroPad</span>(maCoefficients,Nfs))</span>
<span id="cb5-4"><a href="#cb5-4" aria-hidden="true" tabindex="-1"></a> inputVar<span class="sc">*</span><span class="fu">abs</span>(B)<span class="sc">^</span><span class="dv">2</span><span class="sc">/</span><span class="fu">abs</span>(A)<span class="sc">^</span><span class="dv">2</span></span>
<span id="cb5-5"><a href="#cb5-5" aria-hidden="true" tabindex="-1"></a>}</span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
</div>
<p>Let’s use <code>psdArma</code> to calculate the PSD of the ARMA[1,1] process in <a href="#eq-arma11Example">Equation 2</a>.</p>
<div class="cell">
<div class="sourceCode cell-code" id="cb6"><pre class="sourceCode r code-with-copy"><code class="sourceCode r"><span id="cb6-1"><a href="#cb6-1" aria-hidden="true" tabindex="-1"></a>psdArma11 <span class="ot"><-</span> <span class="dv">10</span><span class="sc">*</span><span class="fu">log10</span>(<span class="fu">psdArma</span>(<span class="fu">c</span>(<span class="dv">1</span>,a),<span class="fu">c</span>(<span class="dv">1</span>,b),<span class="at">inputVar=</span><span class="dv">1</span>,<span class="at">Nfs=</span><span class="fu">length</span>(x)))</span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
</div>
<p><a href="#fig-psdArma11">Figure 3</a> contains two plots. The blue plot is the PSD of the ARMA[1,1] process that produced the data in <a href="#fig-samplesArma11">Figure 2</a>. The black plot is the PSD of the WGN input.</p>
<div class="cell">
<div class="cell-output-display">
<div id="fig-psdArma11" class="quarto-figure quarto-figure-center anchored">
<figure class="figure">
<p><img src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjTjq97ay-x3qiAEx-XJ42e8HQ-JWuV_6Xz4u1tUtGK1nMZXcMAjzRNBdNsu6Ex71jjENoxFG7xnZEa2Z9XxpeQWYuCyTd4XvcTO2tgOGmG-glUZEJnl4Q7Za9KLiYUX4FXRTzyPCzyeWJQe-xltR4cOlpDk4p54vqdNtbuLQGUFt4hjUGLPDFPZNDEeL9P/s1600/fig-psdArma11-1.png" class="img-fluid figure-img" width="768"></p>
<p></p><figcaption class="figure-caption">Figure 3: This figure contains two PSDs. The black line at 0 dB is the PSD of the WGN input, and the blue line is the PSD of the ARMA[1,1] process from <a href="#eq-arma11Example">Equation 2</a>. The ARMA[1,1] PSD has a peak at 0.5 and a valley at 0.0 (or 1.0). Remember, the PSD is periodic.</figcaption><p></p>
</figure>
</div>
</div>
</div>
<p>There are three things to notice about the PSD plots. First, the frequency units on the x-axis are a fraction of the sampling frequency (cycles/sample). If you know the sampling frequency (samples/second), you could recover frequency in Hz (cycles/second) by multiplying frequency by the sampling frequency, as <span class="math display">\[
(\text{cycles/sample})(\text{samples/second})=(\text{cycles/second}).\]</span></p>
<p>Second, the y-axis is power in decibels (dB), calculated as <span class="math display">\[\text{power dB} = 10\log_{10}(\text{power}).\]</span></p>
<p>Third, PSDs are periodic with the sampling frequency. So, they can be plotted along any region that shows one period. Below is an example of the PSD from above plotted from -0.5 to 0.5.</p>
<div class="cell">
<div class="cell-output-display">
<div id="fig-psdArma11Alt" class="quarto-figure quarto-figure-center anchored">
<figure class="figure">
<p><img src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjNBuHGtn14dF9D7s6KD3i1Nv1SM81a0WT-h24GYxq7uPcLcwJUtws-Ntut43cj7HLEgKaNs9J3uMZNq_nW09G3h0ZWgJvXXGsnnSN0fkttGaHi5Dat7K6Vx87LwNOpbd7ZfQ_I5hMVNOxvRpws3NdVTfkey7KlZs8G1DqjTEMIXWnr_8sqkxlRMmWg4-nM/s1600/fig-psdArma11Alt-1.png" class="img-fluid figure-img" width="768"></p>
<p></p><figcaption class="figure-caption">Figure 4: This figure is the same as the one above, except it is plotted -0.5 to 0.5.</figcaption><p></p>
</figure>
</div>
</div>
</div>
</section>
<section id="autoregresive-model" class="level2" data-number="3">
<h2 data-number="3" class="anchored" data-anchor-id="autoregresive-model"><span class="header-section-number">3</span> Autoregresive Model</h2>
<p>ARMA[<span class="math inline">\(p\)</span>,0] models are known as autoregressive (AR) models of order <span class="math inline">\(p\)</span>. The AR[<span class="math inline">\(p\)</span>] model is written as <span id="eq-ARmodelEq"><span class="math display">\[
y[n] = -\sum_{k=1}^{p} a[k]y[n-k] + x[n],
\tag{5}\]</span></span> where the <span class="math inline">\(a[k]\)</span>’s are the AR coefficients, <span class="math inline">\(x\)</span> is the WGN input, and <span class="math inline">\(p\)</span> is the process’s order.</p>
<p>We use the AR[1] model below to produce simulated data. We let <span class="math inline">\(a[1]=0.9\)</span>, and the input <span class="math inline">\(x[n]\)</span> is WGN with variance 1. The model is written as <span id="eq-AR1example"><span class="math display">\[
y[n] = -0.9y[n-1] + x[n],
\tag{6}\]</span></span></p>
<p>We use a <code>for</code> loop to implement the ARMA model. We could do the same thing for the AR model, but using the <code>filter</code> command from the <code>gsignal</code> is faster and makes the code more readable. The function inputs the AR coefficients <code>a</code>, the MA coefficients <code>b</code>, and the samples to be filtered <code>x</code>. The function outputs the data filtered by the ARMA coefficients.</p>
<p><a href="#fig-samplesAr1">Figure 5</a> is a plot of the samples produced by the code below.</p>
<div class="cell">
<div class="sourceCode cell-code" id="cb7"><pre class="sourceCode r code-with-copy"><code class="sourceCode r"><span id="cb7-1"><a href="#cb7-1" aria-hidden="true" tabindex="-1"></a><span class="fu">set.seed</span>(<span class="dv">1</span>)</span>
<span id="cb7-2"><a href="#cb7-2" aria-hidden="true" tabindex="-1"></a>N <span class="ot"><-</span> <span class="dv">512</span></span>
<span id="cb7-3"><a href="#cb7-3" aria-hidden="true" tabindex="-1"></a>x <span class="ot"><-</span> <span class="fu">rnorm</span>(N)</span>
<span id="cb7-4"><a href="#cb7-4" aria-hidden="true" tabindex="-1"></a>a <span class="ot"><-</span> <span class="fu">c</span>(<span class="fl">1.0</span>, <span class="fl">0.9</span>)</span>
<span id="cb7-5"><a href="#cb7-5" aria-hidden="true" tabindex="-1"></a>b <span class="ot"><-</span> <span class="dv">1</span></span>
<span id="cb7-6"><a href="#cb7-6" aria-hidden="true" tabindex="-1"></a>y_ar1 <span class="ot"><-</span> gsignal<span class="sc">::</span><span class="fu">filter</span>(b,a,x)</span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
</div>
<div class="cell">
<div class="cell-output-display">
<div id="fig-samplesAr1" class="quarto-figure quarto-figure-center anchored">
<figure class="figure">
<p><img src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgObsQ8AlGkUeEG6b_0KP1Cy0LNF0pREarrAsrNkGjtIlGmKXDyynKGQ06YtUVMO0yBdKC8uctYGb16yNpN7VlKmTEdW0D05zoikDXdZZIXsHlLl_PwGgFbZfSaa25HfiFzPjcfB7tmyBTvRUl5c1mDgSuMMPOq6tOwaSKBuYIBeXN0JsgP31xdcWfKXcav/s1600/fig-samplesAr1-1.png" class="img-fluid figure-img" width="768"></p>
<p></p><figcaption class="figure-caption">Figure 5: This is a plot of the samples produced by the AR[1] model in <a href="#eq-AR1example">Equation 6</a> produced by the code above.</figcaption><p></p>
</figure>
</div>
</div>
</div>
<p>The PSD of an AR process is <span id="eq-ArPsd1def"><span class="math display">\[
P_{\text{AR}}(f) = \frac{\sigma^2 }
{|1 + \sum_{k=1}^p a[k]\exp(-j2\pi fk)|^2 }.
\tag{7}\]</span></span> Since <a href="#eq-ArPsd1def">Equation 7</a> has no zeros in the numerator; the AR PSD only has poles. The poles make AR models good at modeling PSDs with peaks. AR PSD estimators are high resolution estimators because they are good at resolving closely spaced peaks. AR models are good at modeling peaks and are simpler than ARMA models, so they are used more often than the other two.</p>
<p>Below we find the PSD of the AR[1] model above.</p>
<div class="cell">
<div class="sourceCode cell-code" id="cb8"><pre class="sourceCode r code-with-copy"><code class="sourceCode r"><span id="cb8-1"><a href="#cb8-1" aria-hidden="true" tabindex="-1"></a>psdAr1 <span class="ot"><-</span> <span class="dv">10</span><span class="sc">*</span><span class="fu">log10</span>(<span class="fu">psdArma</span>(a,b,<span class="at">inputVar=</span><span class="dv">1</span>,<span class="at">Nfs=</span><span class="fu">length</span>(x)))</span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
</div>
<p><a href="#fig-samplesAr1">Figure 5</a> contains two plots. The blue plot is the PSD of the AR[1] process that produced the data in <a href="#fig-samplesAr1">Figure 5</a>. The black plot is the PSD of the WGN input.</p>
<div class="cell">
<div class="sourceCode cell-code" id="cb9"><pre class="sourceCode r code-with-copy"><code class="sourceCode r"><span id="cb9-1"><a href="#cb9-1" aria-hidden="true" tabindex="-1"></a>psdAr1 <span class="ot"><-</span> <span class="dv">10</span><span class="sc">*</span><span class="fu">log10</span>(<span class="fu">psdArma</span>(a,b,<span class="at">inputVar=</span><span class="dv">1</span>,<span class="at">Nfs=</span><span class="fu">length</span>(x)))</span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
</div>
<div class="cell">
<div class="cell-output-display">
<div id="fig-psdAr1" class="quarto-figure quarto-figure-center anchored">
<figure class="figure">
<p><img src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEg22w9mObiMJ1zSffPUwTwRz9RlPVwgS3hZAlHTC30Mq1dfiSfGUn9LWrf9ahvTjb0U5fsCkwnhdVz-jsCusg_WA0keBedMHynS1i5rE803xFv0B8UIlatnPe3HJfuh4e67z9EOS07Uy5GXjNI1ApXpqgexF91lUR8oxv9k2gQLeXzUtunNDncxQcg_8oKJ/s1600/fig-psdAr1-1.png" class="img-fluid figure-img" width="768"></p>
<p></p><figcaption class="figure-caption">Figure 6: This figure contains two PSDs. The black line at 0 dB is the PSD of the WGN input, and the blue line is the PSD of the AR[1] process from <a href="#eq-ArPsd1def">Equation 7</a>. This PSD has one pole at 0.5 and no zeros.</figcaption><p></p>
</figure>
</div>
</div>
</div>
</section>
<section id="moving-average-model" class="level2" data-number="4">
<h2 data-number="4" class="anchored" data-anchor-id="moving-average-model"><span class="header-section-number">4</span> Moving Average Model</h2>
<p>ARMA[0,<span class="math inline">\(q\)</span>] models are known as moving average (MA) model of order <span class="math inline">\(q\)</span>. An MA[<span class="math inline">\(q\)</span>] is written as <span id="eq-MAmodle"><span class="math display">\[
y[n] = \sum_{k=0}^q b[k]x[n-k],
\tag{8}\]</span></span> where the <span class="math inline">\(b[k]\)</span>’s are the MA coefficients, <span class="math inline">\(x[k]\)</span> is WGN input, and the <span class="math inline">\(q\)</span> is the process’s order. The simple averaging filter is a special case of the MA model, where the MA parameters are <span class="math inline">\(\frac{1}{q+1}\)</span>. The simple averaging filter is written as <span class="math display">\[
y[n] = \frac{1}{q+1}\sum_{k=0}^q x[n-k].
\]</span> Below we use the MA[1] model to produce simulated data. We let <span class="math inline">\(b[1]=-0.9\)</span>, and the input <span class="math inline">\(x[n]\)</span> is WGN with variance 1. The model is written as <span id="eq-MA1example"><span class="math display">\[
y[n] = -0.9x[n-1] + x[n],
\tag{9}\]</span></span> Now we filter the complex WGN samples.</p>
<div class="cell">
<div class="sourceCode cell-code" id="cb10"><pre class="sourceCode r code-with-copy"><code class="sourceCode r"><span id="cb10-1"><a href="#cb10-1" aria-hidden="true" tabindex="-1"></a><span class="fu">set.seed</span>(<span class="dv">1</span>)</span>
<span id="cb10-2"><a href="#cb10-2" aria-hidden="true" tabindex="-1"></a>N <span class="ot"><-</span> <span class="dv">512</span></span>
<span id="cb10-3"><a href="#cb10-3" aria-hidden="true" tabindex="-1"></a>x <span class="ot"><-</span> <span class="fu">rnorm</span>(N)</span>
<span id="cb10-4"><a href="#cb10-4" aria-hidden="true" tabindex="-1"></a>a <span class="ot"><-</span> <span class="dv">1</span></span>
<span id="cb10-5"><a href="#cb10-5" aria-hidden="true" tabindex="-1"></a>b <span class="ot"><-</span> <span class="fu">c</span>(<span class="dv">1</span>,<span class="sc">-</span><span class="fl">0.9</span>)</span>
<span id="cb10-6"><a href="#cb10-6" aria-hidden="true" tabindex="-1"></a>y_ma1 <span class="ot"><-</span> gsignal<span class="sc">::</span><span class="fu">filter</span>(b,a,x)</span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
</div>
<p><a href="#fig-samplesMA1">Figure 7</a> is a plot of the samples produced by the code above.</p>
<div class="cell">
<div class="cell-output-display">
<div id="fig-samplesMA1" class="quarto-figure quarto-figure-center anchored">
<figure class="figure">
<p><img src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgToEx87wEJ0910A_8clvrUhJR4QvRABLhiqGBtZu2RWs2mGkqEKGtS6IUd2As9nZMBwwo0Os2PxDsGFfvK7K61zZ99ERqCAU-l1FiE_MTxRzHPEcjJ_lT-EyNvOnZw-xC_bPZkwa7KOjJpAzeuJbkrry2v06UMMSf3bov2JTujmjQTLfq7TZOHf7un10Ky/s1600/fig-samplesMA1-1.png" class="img-fluid figure-img" width="768"></p>
<p></p><figcaption class="figure-caption">Figure 7: This is a plot of the samples produced by the Ma[1] model in <a href="#eq-MA1example">Equation 9</a>.</figcaption><p></p>
</figure>
</div>
</div>
</div>
<p>The PSD of an MA process is <span id="eq-MaPsd"><span class="math display">\[
P_{\text{MA}}(f) = \sigma^2 |1 + \sum_{k=1}^q b[k]\exp(-j2\pi fk)|^2.
\tag{10}\]</span></span> The MA PSD has zeros, but no poles. So, the MA model is good at modeling PSDs with deep valleys. The MA PSD estimator is not able to resolve closely spaced peaks.</p>
<p>Below we find the PSD of the MA[1] model above.</p>
<div class="cell">
<div class="sourceCode cell-code" id="cb11"><pre class="sourceCode r code-with-copy"><code class="sourceCode r"><span id="cb11-1"><a href="#cb11-1" aria-hidden="true" tabindex="-1"></a>psdMa1 <span class="ot"><-</span> <span class="dv">10</span><span class="sc">*</span><span class="fu">log10</span>(<span class="fu">psdArma</span>(a,b,<span class="at">inputVar=</span><span class="dv">1</span>,<span class="at">Nfs=</span><span class="fu">length</span>(x)))</span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
</div>
<p><a href="#fig-psdMa1">Figure 8</a> contains two plots. The blue plot is the PSD of the MA[1] process that produced the data in <a href="#fig-samplesMA1">Figure 7</a>. The black plot is the PSD of the WGN input.</p>
<div class="cell">
<div class="cell-output-display">
<div id="fig-psdMa1" class="quarto-figure quarto-figure-center anchored">
<figure class="figure">
<p><img src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjM-IbRRxY7ZwXqR8nMw7SANOzs7OdbdPu-lodnEIR-unbJ1yObdZtGWTV87pETzWkRbPBV7rDiIibHxSY2vRUv9LK3IbsD8rtN8vvQztDZbrR4qesZf9Fu7CpyEyUlrBlfXGAXM76xQiwi6Zc2GMklATODJA3yFv1fogIzA0My9BWIjMeY9mU8Bx2M1qTh/s1600/fig-psdMa1-1.png" class="img-fluid figure-img" width="768"></p>
<p></p><figcaption class="figure-caption">Figure 8: This figure contains two PSDs. The black line at 0 dB is the PSD of the WGN input, and the blue line is the PSD of the MA[1] process from <a href="#eq-MA1example">Equation 9</a>. You can see the zero at 0.0 (or 1.0).</figcaption><p></p>
</figure>
</div>
</div>
</div>
</section>
<section id="estimating-psds-using-models" class="level2" data-number="5">
<h2 data-number="5" class="anchored" data-anchor-id="estimating-psds-using-models"><span class="header-section-number">5</span> Estimating PSDs using Models</h2>
<p>ARMA/AR/MA PSD estimation is performed using the following steps. First, choose the model you wish to use. Second, use the data to estimate the model parameters. Parameter estimation is beyond the scope of this article. You can find information in various places, including <span class="citation" data-cites="Carbone2022">Carbone (<a href="#ref-Carbone2022" role="doc-biblioref">2022</a>)</span>, <span class="citation" data-cites="Kay88">Kay (<a href="#ref-Kay88" role="doc-biblioref">1988</a>)</span>, or <span class="citation" data-cites="Marple87">Marple (<a href="#ref-Marple87" role="doc-biblioref">1987</a>)</span>. Finally, calculate the estimated PSD using the estimated parameters in the <code>psdArma</code> function.</p>
</section>
<section id="sec-choosingAModel" class="level2" data-number="6">
<h2 data-number="6" class="anchored" data-anchor-id="sec-choosingAModel"><span class="header-section-number">6</span> Choosing a Model</h2>
<p>If you are only interested in peaks, use the AR model. If only the valleys are important, use the MA model. If both the peaks and valleys are important, use the ARMA model. As mentioned above, peaks are usually the more important feature, so the AR model is often chosen.</p>
</section>
<section id="model-order-selection" class="level2" data-number="7">
<h2 data-number="7" class="anchored" data-anchor-id="model-order-selection"><span class="header-section-number">7</span> Model Order Selection</h2>
<p>An ARMA model can model at most <span class="math inline">\(p\)</span> peaks and <span class="math inline">\(q\)</span> nulls. So, if you know the number of peaks and valleys in the PSD, choosing the model order is straightforward. For instance, if the PSD had one peak and two nulls, you would use an ARMA[1,2] model. Unfortunately, we don’t always know the number of peaks and nulls, but we can use a model order estimator to help us choose the model order.</p>
<p>Model order estimators are mathematical tools that let you estimate the order of the process that produced the data. The Akaike Information Criterion (AIC) is a model order estimator built into the <code>arima</code> function. The AIC calculates how well the model fits the data (lower is better).</p>
<p>The <code>arima</code> function produces several outputs, but we will use the AIC. So we use the <code>arima</code> function with <code>$aic</code>, which will only output the AIC value. The <code>arima$aic</code> function inputs the data and the model order and then outputs the AIC associated with that model order. Below is an example of using the function to find the model order.</p>
<p>Below we use the AR[1] data shown in <a href="#fig-samplesAr1">Figure 5</a> and calculate the AIC for AR model orders 1 to 10.</p>
<div class="cell">
<div class="sourceCode cell-code" id="cb12"><pre class="sourceCode r code-with-copy"><code class="sourceCode r"><span id="cb12-1"><a href="#cb12-1" aria-hidden="true" tabindex="-1"></a>aicAR <span class="ot"><-</span> <span class="fu">rep</span>(<span class="dv">0</span>,<span class="dv">10</span>)</span>
<span id="cb12-2"><a href="#cb12-2" aria-hidden="true" tabindex="-1"></a><span class="cf">for</span> (modelOrder <span class="cf">in</span> <span class="dv">1</span><span class="sc">:</span><span class="dv">10</span>) {</span>
<span id="cb12-3"><a href="#cb12-3" aria-hidden="true" tabindex="-1"></a> aicAR[modelOrder] <span class="ot"><-</span> <span class="fu">arima</span>(y_ar1, <span class="at">order =</span> <span class="fu">c</span>(<span class="at">AR=</span>modelOrder, <span class="dv">0</span>, <span class="at">MA=</span><span class="dv">0</span>))<span class="sc">$</span>aic</span>
<span id="cb12-4"><a href="#cb12-4" aria-hidden="true" tabindex="-1"></a>}</span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
</div>
<p><a href="#fig-aicAR1">Figure 9</a> is a plot of the AIC for AR models 1 to 10. Note the minimum AIC is at AR[1], which is the correct model order.</p>
<div class="cell">
<div class="cell-output-display">
<div id="fig-aicAR1" class="quarto-figure quarto-figure-center anchored">
<figure class="figure">
<p><img src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEi3aNoWb4zNqHyXX4kTvRPjN5t6O1AogHOX9zMxyNdufH3DTPnhyDnttNqlZ320LzvRRHqjKMiiyQ64IHAYc8kyr4-kot0sK_b7sEr4I9zqVyoQu3bnIYPbmgcW6BcVan91OUgAZipOBMGmnOkeXvn2pkjsEffFLZDoNy-ojBHp_eAMIEwDziH_u5kkUOZa/s1600/fig-aicAR1-1.png" class="img-fluid figure-img" width="768"></p>
<p></p><figcaption class="figure-caption">Figure 9: This figure contains a plot of the AIC for model orders 1 to 10. The lowest AIC value is at AR[1], which is the true model order.</figcaption><p></p>
</figure>
</div>
</div>
</div>
<p>However, we choose the model order; we don’t want to overfit the model. We want to choose the minimum model order that can model the important features. Overfitting the model requires more parameters and can add bias and variance to our estimate.</p>
<div id="quarto-appendix" class="default"><section class="quarto-appendix-contents" role="doc-bibliography"><h2 class="anchored quarto-appendix-heading">References</h2><div id="refs" class="references csl-bib-body hanging-indent" role="doc-bibliography">
<div id="ref-Carbone2022" class="csl-entry" role="doc-biblioentry">
Carbone, C. 2022. <em>Power Spectral Density Estimation by Example</em>. Amazon.
</div>
<div id="ref-Kay88" class="csl-entry" role="doc-biblioentry">
Kay, Steven M. 1988. <em>Modern Spectral Estimation</em>. Englewood Cliffs, NJ: Prentice-Hall.
</div>
<div id="ref-Marple87" class="csl-entry" role="doc-biblioentry">
Marple, S. L. 1987. <em>Digital Spectral Analysis with Applications</em>. Englewood Cliffs, NJ: Prentice-Hall.
</div>
</div></section></div></main>
<!-- /main column -->
<script id="quarto-html-after-body" type="application/javascript">
window.document.addEventListener("DOMContentLoaded", function (event) {
const toggleBodyColorMode = (bsSheetEl) => {
const mode = bsSheetEl.getAttribute("data-mode");
const bodyEl = window.document.querySelector("body");
if (mode === "dark") {
bodyEl.classList.add("quarto-dark");
bodyEl.classList.remove("quarto-light");
} else {
bodyEl.classList.add("quarto-light");
bodyEl.classList.remove("quarto-dark");
}
}
const toggleBodyColorPrimary = () => {
const bsSheetEl = window.document.querySelector("link#quarto-bootstrap");
if (bsSheetEl) {
toggleBodyColorMode(bsSheetEl);
}
}
toggleBodyColorPrimary();
const icon = "";
const anchorJS = new window.AnchorJS();
anchorJS.options = {
placement: 'right',
icon: icon
};
anchorJS.add('.anchored');
const clipboard = new window.ClipboardJS('.code-copy-button', {
target: function(trigger) {
return trigger.previousElementSibling;
}
});
clipboard.on('success', function(e) {
// button target
const button = e.trigger;
// don't keep focus
button.blur();
// flash "checked"
button.classList.add('code-copy-button-checked');
var currentTitle = button.getAttribute("title");
button.setAttribute("title", "Copied!");
let tooltip;
if (window.bootstrap) {
button.setAttribute("data-bs-toggle", "tooltip");
button.setAttribute("data-bs-placement", "left");
button.setAttribute("data-bs-title", "Copied!");
tooltip = new bootstrap.Tooltip(button,
{ trigger: "manual",
customClass: "code-copy-button-tooltip",
offset: [0, -8]});
tooltip.show();
}
setTimeout(function() {
if (tooltip) {
tooltip.hide();
button.removeAttribute("data-bs-title");
button.removeAttribute("data-bs-toggle");
button.removeAttribute("data-bs-placement");
}
button.setAttribute("title", currentTitle);
button.classList.remove('code-copy-button-checked');
}, 1000);
// clear code selection
e.clearSelection();
});
function tippyHover(el, contentFn) {
const config = {
allowHTML: true,
content: contentFn,
maxWidth: 500,
delay: 100,
arrow: false,
appendTo: function(el) {
return el.parentElement;
},
interactive: true,
interactiveBorder: 10,
theme: 'quarto',
placement: 'bottom-start'
};
window.tippy(el, config);
}
const noterefs = window.document.querySelectorAll('a[role="doc-noteref"]');
for (var i=0; i<noterefs.length; i++) {
const ref = noterefs[i];
tippyHover(ref, function() {
// use id or data attribute instead here
let href = ref.getAttribute('data-footnote-href') || ref.getAttribute('href');
try { href = new URL(href).hash; } catch {}
const id = href.replace(/^#\/?/, "");
const note = window.document.getElementById(id);
return note.innerHTML;
});
}
const findCites = (el) => {
const parentEl = el.parentElement;
if (parentEl) {
const cites = parentEl.dataset.cites;
if (cites) {
return {
el,
cites: cites.split(' ')
};
} else {
return findCites(el.parentElement)
}
} else {
return undefined;
}
};
var bibliorefs = window.document.querySelectorAll('a[role="doc-biblioref"]');
for (var i=0; i<bibliorefs.length; i++) {
const ref = bibliorefs[i];
const citeInfo = findCites(ref);
if (citeInfo) {
tippyHover(citeInfo.el, function() {
var popup = window.document.createElement('div');
citeInfo.cites.forEach(function(cite) {
var citeDiv = window.document.createElement('div');
citeDiv.classList.add('hanging-indent');
citeDiv.classList.add('csl-entry');
var biblioDiv = window.document.getElementById('ref-' + cite);
if (biblioDiv) {
citeDiv.innerHTML = biblioDiv.innerHTML;
}
popup.appendChild(citeDiv);
});
return popup.innerHTML;
});
}
}
});
</script>
</div> <!-- /content -->
</body></html>Chris Carbonehttp://www.blogger.com/profile/10458001230470047172noreply@blogger.com0tag:blogger.com,1999:blog-7566945490770798840.post-10883600873868181092023-07-18T19:04:00.011-07:002023-07-19T05:11:24.513-07:00Wednesday Linkorama!<!DOCTYPE html>
<html xmlns="http://www.w3.org/1999/xhtml" lang="en" xml:lang="en">
<head>
<meta charset="utf-8">
<meta name="generator" content="quarto-1.1.189">
<meta name="viewport" content="width=device-width, initial-scale=1.0, user-scalable=yes">
<style>
code{white-space: pre-wrap;}
span.smallcaps{font-variant: small-caps;}
div.columns{display: flex; gap: min(4vw, 1.5em);}
div.column{flex: auto; overflow-x: auto;}
div.hanging-indent{margin-left: 1.5em; text-indent: -1.5em;}
ul.task-list{list-style: none;}
ul.task-list li input[type="checkbox"] {
width: 0.8em;
margin: 0 0.8em 0.2em -1.6em;
vertical-align: middle;
}
.display.math{display: block; text-align: center; margin: 0.5rem auto;}
</style>
<link href="linkarama_files/libs/quarto-html/quarto-html.min.css" rel="stylesheet" data-mode="light">
<link href="linkarama_files/libs/quarto-html/quarto-syntax-highlighting.css" rel="stylesheet" id="quarto-text-highlighting-styles">
</head>
<body>
<header id="title-block-header">
</header>
<h4 id="round-up-of-things-i-found-interesting.">Round up of things I found interesting.</h4>
<div class="separator" style="clear: both;"><a href="https://www.sciencenews.org/wp-content/uploads/2023/07/070623_no_anomalocaris_feat.jpg" style="display: block; padding: 1em 0; text-align: center; "><img alt="" border="0" width="400" data-original-height="450" data-original-width="800" src="https://www.sciencenews.org/wp-content/uploads/2023/07/070623_no_anomalocaris_feat.jpg"/></a></div>
<p>The image above came from a <a href="https://www.sciencenews.org/article/ancient-apex-predator-prey">Sciencenews</a> article and it is possibly the first predator. </p>
<p>As a big fan of diet soda, I'm not happy to hear the bad news on <a href="https://www.scientificamerican.com/article/aspartame-declared-possible-carcinogen-heres-what-that-really-means/">aspartame</a>. On the other hand "The ruling does not mean you need to stop consuming all aspartame-containing products, nor does it change the acceptable daily intake put forth by the World Health Organization (WHO). In a separate ruling released at the same time, the Joint Food and Agriculture Organization/WHO Expert Committee on Food Additives (JECFA), which evaluates the levels of exposure to food additives that pose a risk, found there was no convincing evidence for harm with aspartame consumption below the current acceptable daily intake limit of 40 milligrams per kilogram of body weight. That limit was established by JECFA in 1981. For a 70-kilogram (155-pound) person, it’s equivalent to about 14.5 cans of Diet Coke." </p>
<p>This <a href="https://www.youtube.com/watch?v=_EZQx87DyzM"> video</a> shows 40 years of Boston Dynamics progress. It shows that 40 years of incremental progress really adds up.</p>
<p>I liked season 1 of Foundation, and I'm looking forward to <a href="https://www.digitaltrends.com/movies/foundation-interview-season-2-david-goyer-cast/"> Foundation season 2</a>.</p>
</body></html>Chris Carbonehttp://www.blogger.com/profile/10458001230470047172noreply@blogger.com0tag:blogger.com,1999:blog-7566945490770798840.post-60965345817435445112023-07-16T20:19:00.003-07:002023-07-18T15:58:51.225-07:00Autoregressive Estimation and the AIC<!DOCTYPE html>
<html xmlns="http://www.w3.org/1999/xhtml" lang="en" xml:lang="en"><head>
<meta charset="utf-8">
<meta name="generator" content="quarto-1.2.280">
<meta name="viewport" content="width=device-width, initial-scale=1.0, user-scalable=yes">
<!--<title>Autoregressive Estimation and the AIC</title>-->
<style>
code{white-space: pre-wrap;}
span.smallcaps{font-variant: small-caps;}
div.columns{display: flex; gap: min(4vw, 1.5em);}
div.column{flex: auto; overflow-x: auto;}
div.hanging-indent{margin-left: 1.5em; text-indent: -1.5em;}
ul.task-list{list-style: none;}
ul.task-list li input[type="checkbox"] {
width: 0.8em;
margin: 0 0.8em 0.2em -1.6em;
vertical-align: middle;
}
pre > code.sourceCode { white-space: pre; position: relative; }
pre > code.sourceCode > span { display: inline-block; line-height: 1.25; }
pre > code.sourceCode > span:empty { height: 1.2em; }
.sourceCode { overflow: visible; }
code.sourceCode > span { color: inherit; text-decoration: inherit; }
div.sourceCode { margin: 1em 0; }
pre.sourceCode { margin: 0; }
@media screen {
div.sourceCode { overflow: auto; }
}
@media print {
pre > code.sourceCode { white-space: pre-wrap; }
pre > code.sourceCode > span { text-indent: -5em; padding-left: 5em; }
}
pre.numberSource code
{ counter-reset: source-line 0; }
pre.numberSource code > span
{ position: relative; left: -4em; counter-increment: source-line; }
pre.numberSource code > span > a:first-child::before
{ content: counter(source-line);
position: relative; left: -1em; text-align: right; vertical-align: baseline;
border: none; display: inline-block;
-webkit-touch-callout: none; -webkit-user-select: none;
-khtml-user-select: none; -moz-user-select: none;
-ms-user-select: none; user-select: none;
padding: 0 4px; width: 4em;
color: #aaaaaa;
}
pre.numberSource { margin-left: 3em; border-left: 1px solid #aaaaaa; padding-left: 4px; }
div.sourceCode
{ }
@media screen {
pre > code.sourceCode > span > a:first-child::before { text-decoration: underline; }
}
code span.al { color: #ff0000; font-weight: bold; } /* Alert */
code span.an { color: #60a0b0; font-weight: bold; font-style: italic; } /* Annotation */
code span.at { color: #7d9029; } /* Attribute */
code span.bn { color: #40a070; } /* BaseN */
code span.bu { color: #008000; } /* BuiltIn */
code span.cf { color: #007020; font-weight: bold; } /* ControlFlow */
code span.ch { color: #4070a0; } /* Char */
code span.cn { color: #880000; } /* Constant */
code span.co { color: #60a0b0; font-style: italic; } /* Comment */
code span.cv { color: #60a0b0; font-weight: bold; font-style: italic; } /* CommentVar */
code span.do { color: #ba2121; font-style: italic; } /* Documentation */
code span.dt { color: #902000; } /* DataType */
code span.dv { color: #40a070; } /* DecVal */
code span.er { color: #ff0000; font-weight: bold; } /* Error */
code span.ex { } /* Extension */
code span.fl { color: #40a070; } /* Float */
code span.fu { color: #06287e; } /* Function */
code span.im { color: #008000; font-weight: bold; } /* Import */
code span.in { color: #60a0b0; font-weight: bold; font-style: italic; } /* Information */
code span.kw { color: #007020; font-weight: bold; } /* Keyword */
code span.op { color: #666666; } /* Operator */
code span.ot { color: #007020; } /* Other */
code span.pp { color: #bc7a00; } /* Preprocessor */
code span.sc { color: #4070a0; } /* SpecialChar */
code span.ss { color: #bb6688; } /* SpecialString */
code span.st { color: #4070a0; } /* String */
code span.va { color: #19177c; } /* Variable */
code span.vs { color: #4070a0; } /* VerbatimString */
code span.wa { color: #60a0b0; font-weight: bold; font-style: italic; } /* Warning */
</style>
<script src="Model Order Selection_files/libs/clipboard/clipboard.min.js"></script>
<script src="Model Order Selection_files/libs/quarto-html/quarto.js"></script>
<script src="Model Order Selection_files/libs/quarto-html/popper.min.js"></script>
<script src="Model Order Selection_files/libs/quarto-html/tippy.umd.min.js"></script>
<script src="Model Order Selection_files/libs/quarto-html/anchor.min.js"></script>
<link href="Model Order Selection_files/libs/quarto-html/tippy.css" rel="stylesheet">
<link href="Model Order Selection_files/libs/quarto-html/quarto-syntax-highlighting.css" rel="stylesheet" id="quarto-text-highlighting-styles">
<script src="Model Order Selection_files/libs/bootstrap/bootstrap.min.js"></script>
<link href="Model Order Selection_files/libs/bootstrap/bootstrap-icons.css" rel="stylesheet">
<link href="Model Order Selection_files/libs/bootstrap/bootstrap.min.css" rel="stylesheet" id="quarto-bootstrap" data-mode="light">
</head>
<body class="fullcontent">
<div id="quarto-content" class="page-columns page-rows-contents page-layout-article">
<main class="content" id="quarto-document-content">
<!--<header id="title-block-header" class="quarto-title-block default">
<div class="quarto-title">
<h1 class="title">Autoregressive Estimation and the AIC</h1>
</div>-->
<div class="quarto-title-meta">
</div>
</header>
<section id="model-order-selection" class="level2" data-number="1">
<h2 data-number="1" class="anchored" data-anchor-id="model-order-selection"><span class="header-section-number">1</span> Model Order Selection</h2>
<p> This section provides a brief introduction to autoregressive (AR) model order estimation. An AR model of order p (AR[p]) can model up to p peaks in a PSD. So, if we wanted to estimate a PSD with 4 peaks we would use an AR[4] model. Unfortunately, we don’t always know the number of peaks. Note, however we choose the model order we don’t want to over fit the model. Meaning we want to choose the minimum model order that can model the important features. Over fitting the model requires more parameters and adds bias and variance to our estimate.</p>
<p>Model order estimators are mathematical tools that let you estimate the order of the process that produced the data. The Akaike Information Criterion (AIC) is a model order estimator built into the <code>arima</code> function. The AIC calculates how well model fits the data (lower is better).</p>
<p>The <code>arima</code> function inputs the model order and the data. The function produces a number of outputs, but we only use the AIC. Using the <code>arima</code> function with <code>$aic</code> outputs the AIC associated with the input model order. Below is an example of using the <code>arima</code> function to find the model order.</p>
</section>
<section id="example" class="level2" data-number="2">
<h2 data-number="2" class="anchored" data-anchor-id="example"><span class="header-section-number">2</span> Example</h2>
<p>In this example we use data from an AR model and then use the AIC to estimate the order of the AR model that produced the data. Note the example is implemented in the R programing language. First, we create data using an AR[1] model.</p>
<div class="cell">
<div class="sourceCode cell-code" id="cb1"><pre class="sourceCode r code-with-copy"><code class="sourceCode r"><span id="cb1-1"><a href="#cb1-1" aria-hidden="true" tabindex="-1"></a><span class="fu">set.seed</span>(<span class="dv">1</span>) </span>
<span id="cb1-2"><a href="#cb1-2" aria-hidden="true" tabindex="-1"></a>N <span class="ot"><-</span> <span class="dv">512</span></span>
<span id="cb1-3"><a href="#cb1-3" aria-hidden="true" tabindex="-1"></a>y_ar1 <span class="ot"><-</span> <span class="fu">rep</span>(<span class="dv">0</span>,N)</span>
<span id="cb1-4"><a href="#cb1-4" aria-hidden="true" tabindex="-1"></a>x <span class="ot"><-</span> <span class="fu">rnorm</span>(N)</span>
<span id="cb1-5"><a href="#cb1-5" aria-hidden="true" tabindex="-1"></a>a <span class="ot"><-</span> <span class="fl">0.9</span> </span>
<span id="cb1-6"><a href="#cb1-6" aria-hidden="true" tabindex="-1"></a></span>
<span id="cb1-7"><a href="#cb1-7" aria-hidden="true" tabindex="-1"></a><span class="cf">for</span> (n <span class="cf">in</span> <span class="dv">2</span><span class="sc">:</span>N) </span>
<span id="cb1-8"><a href="#cb1-8" aria-hidden="true" tabindex="-1"></a> y_ar1[n] <span class="ot"><-</span> <span class="sc">-</span>a<span class="sc">*</span>y_ar1[n<span class="dv">-1</span>] <span class="sc">+</span> x[n] </span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
</div>
<p>Next, we use the data to estimate the model order. We calculate the AIC for AR models AR[1] to AR[10]. The model with the lowest AIC is the model that fits the data. If two models have similar AICs we should use the lower order model.</p>
<div class="cell">
<div class="sourceCode cell-code" id="cb2"><pre class="sourceCode r code-with-copy"><code class="sourceCode r"><span id="cb2-1"><a href="#cb2-1" aria-hidden="true" tabindex="-1"></a>aicAR <span class="ot"><-</span> <span class="fu">rep</span>(<span class="dv">0</span>,<span class="dv">10</span>)</span>
<span id="cb2-2"><a href="#cb2-2" aria-hidden="true" tabindex="-1"></a><span class="cf">for</span> (modelOrder <span class="cf">in</span> <span class="dv">1</span><span class="sc">:</span><span class="dv">10</span>) {</span>
<span id="cb2-3"><a href="#cb2-3" aria-hidden="true" tabindex="-1"></a> aicAR[modelOrder] <span class="ot"><-</span> <span class="fu">arima</span>(y_ar1, <span class="at">order =</span> <span class="fu">c</span>(<span class="at">AR=</span>modelOrder, <span class="dv">0</span>, <span class="at">MA=</span><span class="dv">0</span>))<span class="sc">$</span>aic</span>
<span id="cb2-4"><a href="#cb2-4" aria-hidden="true" tabindex="-1"></a>}</span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
</div>
<p><a href="#fig-aicAR1">Figure 1</a> show the AIC value for AR model orders from 1 to 10.</p>
<div class="cell">
<div class="cell-output-display">
<div id="fig-aicAR1" class="quarto-figure quarto-figure-center anchored">
<figure class="figure"><p><img src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEiRwwJOMomofLBGQCmp8NU84cMxJ8rFbWwLi4w8I7M5lg75DjEo9DL4XhRliU55alHqP3efmU5G1Y4e1cROFdIIGO-H7bAjTN9DOkXaLoWq9uNsRBWMqlhnLnj2IicOwFjwb0fWmKdMWr7K60jrdSn-37O_OphCNo8VTpCup-yhrZmlTM_x5NLEKuroywSG/s1600/fig-aicAR1-1.png" class="img-fluid figure-img" width="768"></p>
<p></p><figcaption class="figure-caption">Figure 1: This figure contains a plot of the AIC for model orders 1 to 10. The lowest AIC value is from the AR[1] model, which is the true model order. The AIC is not perfect and sometimes will not get the exact right answer.</figcaption><p></p>
</figure>
</div>
</div>
</div>
<p>Note a similar method could be use to find the model order of moving average (MA) and autoregressive moving avarage (ARMA) models.</p>
</section>
</main>
<!-- /main column -->
<script id="quarto-html-after-body" type="application/javascript">
window.document.addEventListener("DOMContentLoaded", function (event) {
const toggleBodyColorMode = (bsSheetEl) => {
const mode = bsSheetEl.getAttribute("data-mode");
const bodyEl = window.document.querySelector("body");
if (mode === "dark") {
bodyEl.classList.add("quarto-dark");
bodyEl.classList.remove("quarto-light");
} else {
bodyEl.classList.add("quarto-light");
bodyEl.classList.remove("quarto-dark");
}
}
const toggleBodyColorPrimary = () => {
const bsSheetEl = window.document.querySelector("link#quarto-bootstrap");
if (bsSheetEl) {
toggleBodyColorMode(bsSheetEl);
}
}
toggleBodyColorPrimary();
const icon = "";
const anchorJS = new window.AnchorJS();
anchorJS.options = {
placement: 'right',
icon: icon
};
anchorJS.add('.anchored');
const clipboard = new window.ClipboardJS('.code-copy-button', {
target: function(trigger) {
return trigger.previousElementSibling;
}
});
clipboard.on('success', function(e) {
// button target
const button = e.trigger;
// don't keep focus
button.blur();
// flash "checked"
button.classList.add('code-copy-button-checked');
var currentTitle = button.getAttribute("title");
button.setAttribute("title", "Copied!");
let tooltip;
if (window.bootstrap) {
button.setAttribute("data-bs-toggle", "tooltip");
button.setAttribute("data-bs-placement", "left");
button.setAttribute("data-bs-title", "Copied!");
tooltip = new bootstrap.Tooltip(button,
{ trigger: "manual",
customClass: "code-copy-button-tooltip",
offset: [0, -8]});
tooltip.show();
}
setTimeout(function() {
if (tooltip) {
tooltip.hide();
button.removeAttribute("data-bs-title");
button.removeAttribute("data-bs-toggle");
button.removeAttribute("data-bs-placement");
}
button.setAttribute("title", currentTitle);
button.classList.remove('code-copy-button-checked');
}, 1000);
// clear code selection
e.clearSelection();
});
function tippyHover(el, contentFn) {
const config = {
allowHTML: true,
content: contentFn,
maxWidth: 500,
delay: 100,
arrow: false,
appendTo: function(el) {
return el.parentElement;
},
interactive: true,
interactiveBorder: 10,
theme: 'quarto',
placement: 'bottom-start'
};
window.tippy(el, config);
}
const noterefs = window.document.querySelectorAll('a[role="doc-noteref"]');
for (var i=0; i<noterefs.length; i++) {
const ref = noterefs[i];
tippyHover(ref, function() {
// use id or data attribute instead here
let href = ref.getAttribute('data-footnote-href') || ref.getAttribute('href');
try { href = new URL(href).hash; } catch {}
const id = href.replace(/^#\/?/, "");
const note = window.document.getElementById(id);
return note.innerHTML;
});
}
const findCites = (el) => {
const parentEl = el.parentElement;
if (parentEl) {
const cites = parentEl.dataset.cites;
if (cites) {
return {
el,
cites: cites.split(' ')
};
} else {
return findCites(el.parentElement)
}
} else {
return undefined;
}
};
var bibliorefs = window.document.querySelectorAll('a[role="doc-biblioref"]');
for (var i=0; i<bibliorefs.length; i++) {
const ref = bibliorefs[i];
const citeInfo = findCites(ref);
if (citeInfo) {
tippyHover(citeInfo.el, function() {
var popup = window.document.createElement('div');
citeInfo.cites.forEach(function(cite) {
var citeDiv = window.document.createElement('div');
citeDiv.classList.add('hanging-indent');
citeDiv.classList.add('csl-entry');
var biblioDiv = window.document.getElementById('ref-' + cite);
if (biblioDiv) {
citeDiv.innerHTML = biblioDiv.innerHTML;
}
popup.appendChild(citeDiv);
});
return popup.innerHTML;
});
}
}
});
</script>
</div> <!-- /content -->
</body></html>Chris Carbonehttp://www.blogger.com/profile/10458001230470047172noreply@blogger.com0tag:blogger.com,1999:blog-7566945490770798840.post-86121627188869574892022-12-07T17:54:00.006-08:002023-07-17T04:15:30.477-07:00Wednesday Linkorama!<!DOCTYPE html>
<html xmlns="http://www.w3.org/1999/xhtml" lang="en" xml:lang="en">
<head>
<meta charset="utf-8">
<meta name="generator" content="quarto-1.1.189">
<meta name="viewport" content="width=device-width, initial-scale=1.0, user-scalable=yes">
<style>
code{white-space: pre-wrap;}
span.smallcaps{font-variant: small-caps;}
div.columns{display: flex; gap: min(4vw, 1.5em);}
div.column{flex: auto; overflow-x: auto;}
div.hanging-indent{margin-left: 1.5em; text-indent: -1.5em;}
ul.task-list{list-style: none;}
ul.task-list li input[type="checkbox"] {
width: 0.8em;
margin: 0 0.8em 0.2em -1.6em;
vertical-align: middle;
}
.display.math{display: block; text-align: center; margin: 0.5rem auto;}
</style>
<link href="linkarama_files/libs/quarto-html/quarto-html.min.css" rel="stylesheet" data-mode="light">
<link href="linkarama_files/libs/quarto-html/quarto-syntax-highlighting.css" rel="stylesheet" id="quarto-text-highlighting-styles">
</head>
<body>
<header id="title-block-header">
</header>
<h4 id="round-up-of-things-i-found-interesting.">Round up of things I found interesting.</h4>
<p><img src="https://i.kinja-img.com/gawker-media/image/upload/c_fit,f_auto,g_center,pg_1,q_60,w_400/daf21914a855f26a9c7d03b840fca927.jpg"></p>
<p>The <a href="https://gizmodo.com/orion-artemis-1-moon-flyby-photos-nasa-1849863344">photo</a> above came from a Gizomodo article about NASA’s Latest Artemis 1 Moon Images. Highly recomed you take a look!</p>
<p><a href="https://www.techdirt.com/2022/12/07/san-francisco-legislators-greenlight-killing-of-residents-by-police-robots/">San Francisco</a> narrowly avoids killer robots! Look like the city was thinking about deploying police robots that carried explosives.</p>
<p>How <a href="https://www.scientificamerican.com/article/the-linguistics-of-swearing-explain-why-we-substitute-darn-for-damn/">damn became darn</a>. The explanation of how and why use substitute swears.</p>
</body></html>Chris Carbonehttp://www.blogger.com/profile/10458001230470047172noreply@blogger.com0tag:blogger.com,1999:blog-7566945490770798840.post-17675571590295963592022-12-04T08:57:00.002-08:002022-12-04T11:01:31.814-08:00The Periodogram or Something Else?<p> </p>
<p>The periodogram is the simplest and most common spectral estimator and, in R is implemented as </p>
<code>findPer <- function(x){(1/length(x))*abs(fft(x))^2}, </code>
<p>but there are may othere spectral estimators<sup>*</sup>.</p>
<p>So, which spectral estimator should you pick? The flow chart below presents a reasonable way to pick between spectral estimators. The first question you should ask is, do you want to assume a model for the data? If you don’t know anything about the data or don’t trust what you know, then move to the left side of the flow chart. If the Fourier methods have enough resolution, use the periodogram or the Blackman-Tukey spectral estimator (BTSE). If they don’t, you could try the minimum variance spectral estimator (MVSE). If the MVSE does not have enough resolution, try a high order autoregressive (AR) estimate. The AR method does assume a model for the data, but a high enough AR model can estimate any PSD.</p>
<p>If you want to assume a model, we move to the right side of the flow chart. If we know the data consists of sinusoids in noise, then the multiple signal clarifier (MUSIC) estimator is a good choice. If the important features are peaks, but the data is not made up of sinusoids in noise, the AR estimator is a good choice. If nulls are the important features, then use a moving average (MA) estimator. If the important features are nulls and peaks, use the autoregressive moving average (ARMA) estimator. Finally, if you have no idea about the PSD shape and want to use a model, you could try a high order AR.</p>
<p> </p>
<div class="separator" style="clear: both;"><a href="https://github.com/carbone1853/psdFlowChart/blob/main/FlowChart.png?raw=true" style="display: block; padding: 1em 0px; text-align: center;"><img alt="" border="0" data-original-height="725" data-original-width="800" src="https://github.com/carbone1853/psdFlowChart/blob/main/FlowChart.png?raw=true" width="400" /></a></div>
<p> * For more info on spectral estimation, see my new book! -> <a href="https://www.amazon.com/Power-Spectral-Density-Estimation-Example/dp/B0BHWDC4VW">Spectral Density Estimation by Example</a></p>Chris Carbonehttp://www.blogger.com/profile/10458001230470047172noreply@blogger.com0tag:blogger.com,1999:blog-7566945490770798840.post-31071086759526717012022-11-30T17:35:00.001-08:002022-11-30T17:35:31.716-08:00Wednesday Linkorama<h4>Round up of things I found interesting. </h4>
<p> <a href="https://www.scientificamerican.com/article/space-elevators-are-less-sci-fi-than-you-think/" target="_blank">Space elevators</a> have been a part of SiFi for a long time, this article talks about possibly building them for real! Space elevators are cables that go from the ground into space that a vehicle can use to climb into space. This article describes the basic ideas about how they work and why you may want to use one instead of rockets.</p>
<p> <a href="https://www.wired.com/story/us-congress-ftx-cryptocurrency-regulation/" target="_blank">Crypto</a> is starting to get the attention of regulators. </p>
<p> <a href="https://arstechnica.com/information-technology/2022/11/meta-researchers-create-ai-that-masters-diplomacy-tricking-human-players/" target="_blank">Diplomacy</a> is a classic social board game focused on trickery and negotiations between players. The first games computers got really good at (like chess) required a lot of calculation. This article talks about Cicero, an AI that (allegedly) has human level performance in the game Diplomacy. So, computers may be getting better at other types of games. Chris Carbonehttp://www.blogger.com/profile/10458001230470047172noreply@blogger.com0tag:blogger.com,1999:blog-7566945490770798840.post-85017273217229303552022-11-02T18:57:00.002-07:002022-11-02T18:57:08.645-07:00Wednesday Linkorama
<h4><span style="font-family: arial;">Round up of things I found interesting. </span></h4>
<p> <span style="font-family: arial;"><a href="https://www.techdirt.com/2022/10/27/anti-cheat-software-continues-to-be-the-new-drm-in-pissing-off-legit-customers/"></a></span></p>
<p><span style="font-family: arial;">A Reddit user bought an <a href="https://arstechnica.com/information-technology/2022/10/redditor-acquires-decommissioned-netflix-cache-server-with-262tb-of-storage/"><b>old Netflix server</b></a> and reconfigured it for personal storage.</span></p><p><span style="font-family: arial;"><a href="https://arstechnica.com/information-technology/2022/11/metas-ai-powered-audio-codec-promises-10x-compression-over-mp3/">New compression</a> claims 10x over MP3. MP3 is already compressed so ... If this is 10x and sounds as good then they really have something.</span></p><p><a href="https://techcrunch.com/2022/11/02/aurora-says-it-has-enough-cash-to-commercialize-autonomous-trucks-by-2024/">Normally</a><span style="font-family: arial;"> autonomous trucks are 5 years out, but Aurora says they are coming in 2024. So, maybe they are on to something?</span><br /></p><p><span style="font-family: arial;"><br /></span></p>
<p> </p>
Chris Carbonehttp://www.blogger.com/profile/10458001230470047172noreply@blogger.com0tag:blogger.com,1999:blog-7566945490770798840.post-70992677498936360822022-10-26T17:46:00.010-07:002022-10-26T22:50:10.319-07:00Wednesday Linkorama
<h4>Round up of things I found interesting. </h4>
<p> <a href="https://cran.r-project.org/web/packages/gsignal/vignettes/gsignal.html"
target="_blank">gsignal</a
> is an R package that is an improvement over the signal package. The signal package is a group of
signal processing functions like you might find in Matlab. It was helpful but
was missing some common functions, and some functions did not play well with complex numbers. The gsingal packaged added some key functions and works with complex numbers. It has the poly function but seems to be missing the roots function. </p>
<p> <a href="https://arstechnica.com/gadgets/2022/10/macos-13-ventura-the-ars-technica-review/" target="_blank">Ventura</a>, the new mac operating system, is available. I'm going for it. Wish me luck! </p>
<p><a href="https://github.com/features/copilot" target="_blank">Copilot</a> is an app that helps write code. It uses what you wrote to suggest code snippets as you write the code. It looks interesting; I'm going to give it a try. </p>
Chris Carbonehttp://www.blogger.com/profile/10458001230470047172noreply@blogger.com0tag:blogger.com,1999:blog-7566945490770798840.post-58349718645281145722022-10-25T10:08:00.010-07:002023-11-20T06:29:50.397-08:00My new book is published!
<a href="https://www.amazon.com/Power-Spectral-Density-Estimation-Example/dp/B0BHWDC4VW"><div class="separator" style="clear: both;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEggOyYMQ98YkVN3OTOqAtpIVT6Pjp02XeDE0kQZ_5asDO8ftnQMJ3yRj92w8gUztVHbRIj-3DzOAtdxmI9EhcXZiFr8BzmOPo6-5ctk-7-AEG2Y8CgyaeQLSJi6JJetHZKDCB4xF0Ri9324VahlOgWzFDIxFv32H2Mpvyzrtd1mT-daJQF13y5739l97sGm/s1600/book%20cover%20link.jpeg" style="display: block; padding: 1em 0; text-align: center; "><img alt="" border="0" height="320" data-original-height="3399" data-original-width="2473" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEggOyYMQ98YkVN3OTOqAtpIVT6Pjp02XeDE0kQZ_5asDO8ftnQMJ3yRj92w8gUztVHbRIj-3DzOAtdxmI9EhcXZiFr8BzmOPo6-5ctk-7-AEG2Y8CgyaeQLSJi6JJetHZKDCB4xF0Ri9324VahlOgWzFDIxFv32H2Mpvyzrtd1mT-daJQF13y5739l97sGm/s1600/book%20cover%20link.jpeg"/></a></div></a>
<p>This book provides the tools to use, understand, and apply power spectral density (PSD) estimation. By the end of this book, you will be able to select an appropriate PSD estimation method based on your data, perform the estimation, and understand the results. To accomplish this, we cover the methods and just enough theory so you can use the estimation techniques correctly and interpret the results meaningfully.</p>
Get it here -> <a href="https://www.amazon.com/Power-Spectral-Density-Estimation-Example/dp/B0BHWDC4VW">Spectral Density Estimation by Example</a>
Chris Carbonehttp://www.blogger.com/profile/10458001230470047172noreply@blogger.com0tag:blogger.com,1999:blog-7566945490770798840.post-7209759979733920762020-08-27T19:59:00.009-07:002020-08-27T20:04:55.181-07:00Filtering Complex Data with R<html>
<head>
<style type="text/css">
.inline {
background-color: #f7f7f7;
border:solid 1px #B0B0B0;
}
.error {
font-weight: bold;
color: #FF0000;
}
.warning {
font-weight: bold;
}
.message {
font-style: italic;
}
.source, .output, .warning, .error, .message {
padding: 0 1em;
border:solid 1px #F7F7F7;
}
.source {
background-color: #f5f5f5;
}
.left {
text-align: left;
}
.right {
text-align: right;
}
.center {
text-align: center;
}
.hl.num {
color: #AF0F91;
}
.hl.str {
color: #317ECC;
}
.hl.com {
color: #AD95AF;
font-style: italic;
}
.hl.opt {
color: #000000;
}
.hl.std {
color: #585858;
}
.hl.kwa {
color: #295F94;
font-weight: bold;
}
.hl.kwb {
color: #B05A65;
}
.hl.kwc {
color: #55aa55;
}
.hl.kwd {
color: #BC5A65;
font-weight: bold;
}
</style>
<title> </title>
</head>
<body>
<h4>Filtering</h4>
<p>Base R has a function filter that can perform moving average (MA) or filters (AR). But it cannot do both at the same time; it cannot do an autoregresive ARMA filter. It also only works on real data. The filter function in the signal package can perform ARMA filters, but is is built on the filter function in the base package. So, it also cannot filter complex data. Below is an example of filtering some data with the filter from the signal package. We are using a 6 sample MA filter as specified by the b=rep(1/6,6).</p>
<div class="chunk" id="unnamed-chunk-1"><div class="rcode"><div class="source"><pre class="knitr r"><span class="hl std">N</span> <span class="hl kwb"><-</span> <span class="hl num">64</span>
<span class="hl std">x</span> <span class="hl kwb"><-</span> <span class="hl kwd">rnorm</span><span class="hl std">(N)</span>
<span class="hl std">b</span> <span class="hl kwb"><-</span> <span class="hl kwd">rep</span><span class="hl std">(</span><span class="hl num">1</span><span class="hl opt">/</span><span class="hl num">6</span><span class="hl std">,</span><span class="hl num">6</span><span class="hl std">)</span>
<span class="hl std">a</span> <span class="hl kwb"><-</span> <span class="hl num">1</span>
<span class="hl std">y</span> <span class="hl kwb"><-</span> <span class="hl std">signal</span><span class="hl opt">::</span><span class="hl kwd">filter</span><span class="hl std">(b,a,x)</span>
</pre></div>
</div></div>
<p>The figure below has two plots. The blue circles are the Gaussian samples. The green circles are the filtered Gaussian samples. The filtered data has less variance.<p>
<!-- ******** PLOT filtered y ********** -->
<div class="separator" style="clear: both;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjARzO43-Oe_u2-N0B6fj9LhhMN8YMNd-lJ-CWLA5L6YjjkGlnOGbubq5UpcV4ZDBWE3k95woYXrPlfvtMoQT30oMIaTFwf-qjo14-i7s3AY-2SNslhINu_GrLfCWS2y9CiOMBAu-le_NuF/s800/rNorm.png" style="display: block; padding: 1em 0; text-align: none;"><img alt="" border="0" width="400" data-original-height="600" data-original-width="800" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjARzO43-Oe_u2-N0B6fj9LhhMN8YMNd-lJ-CWLA5L6YjjkGlnOGbubq5UpcV4ZDBWE3k95woYXrPlfvtMoQT30oMIaTFwf-qjo14-i7s3AY-2SNslhINu_GrLfCWS2y9CiOMBAu-le_NuF/s400/rNorm.png"/></a></div>
<h4>Complex Data</h4>
<p>In engineering the most common kind of complex data is the circularly symmetric gaussian noise. The circularly symmetric part means that the real and complex part are independent and identically distributed. I described the makeCG function in another post, but as you can see the real part and the imaginary part are made using separate calls to the rnorm function.</p>
<div class="chunk" id="unnamed-chunk-3"><div class="rcode"><div class="source"><pre class="knitr r"><span class="hl std">makeCG</span> <span class="hl kwb"><-</span> <span class="hl kwa">function</span><span class="hl std">(</span><span class="hl kwc">N</span><span class="hl std">,</span><span class="hl kwc">v</span><span class="hl std">=</span><span class="hl num">1</span><span class="hl std">) {(</span><span class="hl kwd">sqrt</span><span class="hl std">(v</span><span class="hl opt">/</span><span class="hl num">2</span><span class="hl std">))</span><span class="hl opt">*</span><span class="hl kwd">rnorm</span><span class="hl std">(N)</span> <span class="hl opt">+</span> <span class="hl std">(</span><span class="hl kwd">sqrt</span><span class="hl std">(v</span><span class="hl opt">/</span><span class="hl num">2</span><span class="hl std">))</span><span class="hl opt">*</span><span class="hl num">1i</span><span class="hl opt">*</span><span class="hl kwd">rnorm</span><span class="hl std">(N)}</span>
</pre></div>
</div></div>
<p>The function filterComplex filters real or complex data. The function inputs the MA parameters b and the AR parameters a, and the data x. If the data is complex, it makes two filter calls: one with the real part and one with the imaginary part. If the data is real it make one call.</p>
<div class="chunk" id="unnamed-chunk-4"><div class="rcode"><div class="source"><pre class="knitr r"><span class="hl std">filterComplex</span> <span class="hl kwb"><-</span> <span class="hl kwa">function</span><span class="hl std">(</span><span class="hl kwc">b</span><span class="hl std">,</span> <span class="hl kwc">a</span><span class="hl std">,</span> <span class="hl kwc">x</span><span class="hl std">) {</span>
<span class="hl kwa">if</span> <span class="hl std">(</span><span class="hl kwd">is.complex</span><span class="hl std">(x)) {</span>
<span class="hl std">y</span> <span class="hl kwb"><-</span> <span class="hl std">signal</span><span class="hl opt">::</span><span class="hl kwd">filter</span><span class="hl std">(b,a,</span><span class="hl kwd">Re</span><span class="hl std">(x))</span> <span class="hl opt">+</span> <span class="hl num">1i</span><span class="hl opt">*</span><span class="hl std">signal</span><span class="hl opt">::</span><span class="hl kwd">filter</span><span class="hl std">(b,a,</span><span class="hl kwd">Im</span><span class="hl std">(x))</span>
<span class="hl std">}</span> <span class="hl kwa">else</span> <span class="hl std">y</span> <span class="hl kwb"><-</span> <span class="hl std">signal</span><span class="hl opt">::</span><span class="hl kwd">filter</span><span class="hl std">(b,a,x)</span>
<span class="hl kwd">return</span><span class="hl std">(y)</span>
<span class="hl std">}</span>
</pre></div>
</div></div>
<p>Here we simulate complex Gaussian data and use the function filterComplex to filter it.</p>
<div class="chunk" id="unnamed-chunk-5"><div class="rcode"><div class="source"><pre class="knitr r"><span class="hl std">N</span> <span class="hl kwb"><-</span> <span class="hl num">64</span>
<span class="hl std">x</span> <span class="hl kwb"><-</span> <span class="hl kwd">makeCG</span><span class="hl std">(N)</span>
<span class="hl std">b</span> <span class="hl kwb"><-</span> <span class="hl kwd">rep</span><span class="hl std">(</span><span class="hl num">1</span><span class="hl opt">/</span><span class="hl num">6</span><span class="hl std">,</span><span class="hl num">6</span><span class="hl std">)</span>
<span class="hl std">a</span> <span class="hl kwb"><-</span> <span class="hl num">1</span>
<span class="hl std">y</span> <span class="hl kwb"><-</span> <span class="hl kwd">filterComplex</span><span class="hl std">(b,a,x)</span>
</pre></div>
</div></div>
<p>The figure below has two plots. The blue circles are the real part of the complex Gaussian samples. The green circles are the real part of filtered complex Gaussian samples. As expected, the filtered data has less variance.<p>
<!-- ******** PLOT complex filtered data. ********** -->
<div class="separator" style="clear: both;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjEfkkxp6G3Z0EqWVsG9BN1MDPcA5iQzUqCKokSqDoZW7HVwjRoJig6FTZ_YUXRjTqHhQhOHTFIESwJVf-nVgNyoLByzXcIbniJBuI8pWNxkZg7b-ECiOFTlCFonVYmkaUPXiq0KkwBbMlp/s800/rComplex.png" style="display: block; padding: 1em 0; text-align: none;"><img alt="" border="0" width="400" data-original-height="600" data-original-width="800" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjEfkkxp6G3Z0EqWVsG9BN1MDPcA5iQzUqCKokSqDoZW7HVwjRoJig6FTZ_YUXRjTqHhQhOHTFIESwJVf-nVgNyoLByzXcIbniJBuI8pWNxkZg7b-ECiOFTlCFonVYmkaUPXiq0KkwBbMlp/s400/rComplex.png"/></a></div>
<p>Now the imaginary part. The blue circles are the imaginary part of the complex Gaussian samples. The green circles are the imaginary part of filtered complex Gaussian samples.<p>
<div class="separator" style="clear: both;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEg6zwG3R3SOs9BxWvzz7anr8mgalJJvvTqwhScg32rDiXu0Xf07odKOfUU2WosmEv3SdyigsT12Tx3zAvB0MAkdEBVaS8siFyYkGpbMAMnikxVig9N2TXecqgsnB9FVDXKCdil1il9iTgMd/s800/rComplexIm.png" style="display: block; padding: 1em 0; text-align: none;"><img alt="" border="0" width="400" data-original-height="600" data-original-width="800" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEg6zwG3R3SOs9BxWvzz7anr8mgalJJvvTqwhScg32rDiXu0Xf07odKOfUU2WosmEv3SdyigsT12Tx3zAvB0MAkdEBVaS8siFyYkGpbMAMnikxVig9N2TXecqgsnB9FVDXKCdil1il9iTgMd/s400/rComplexIm.png"/></a></div>
</body>
</html>
Chris Carbonehttp://www.blogger.com/profile/10458001230470047172noreply@blogger.com0tag:blogger.com,1999:blog-7566945490770798840.post-14890787446988564442020-06-24T11:17:00.002-07:002020-06-24T11:17:34.477-07:00Dungeons and Dragons: Advantage<html>
<head>
<style type="text/css">
.inline {
background-color: #f7f7f7;
border:solid 1px #B0B0B0;
}
.error {
font-weight: bold;
color: #FF0000;
}
.warning {
font-weight: bold;
}
.message {
font-style: italic;
}
.source, .output, .warning, .error, .message {
padding: 0 1em;
border:solid 1px #F7F7F7;
}
.source {
background-color: #f5f5f5;
}
.left {
text-align: left;
}
.right {
text-align: right;
}
.center {
text-align: center;
}
.hl.num {
color: #AF0F91;
}
.hl.str {
color: #317ECC;
}
.hl.com {
color: #AD95AF;
font-style: italic;
}
.hl.opt {
color: #000000;
}
.hl.std {
color: #585858;
}
.hl.kwa {
color: #295F94;
font-weight: bold;
}
.hl.kwb {
color: #B05A65;
}
.hl.kwc {
color: #55aa55;
}
.hl.kwd {
color: #BC5A65;
font-weight: bold;
}
</style>
<title> </title>
</head>
<body>
<h4>D20 and Random Events</h4>
<p>
In the game Dungeons and Dragons, the success or failure of an event is determined by rolling a 20 sided die (D20): higher is better. If you need to roll 11 or higher you have a 50% chance of success. If another event requires 10 or better you now have a 55% chance for success. Each point the roll goes up or down is worth 5%. Often times, rolling a 20 is a
critical success and a 1 is a critical failure; critical means it's extra good/bad.
</p>
<h4>Advantage</h4>
<p>
Sometimes things are really going your way and you roll with advantage.
Sometimes things are not in your favor and you roll with disadvantage. When rolling with advantage, roll two dice and pick the higher. When rolling with disadvantage, roll two dice and pick the lower. How do advantage and disadvantage affect the chance of getting a 20 or 1?
</p>
<p>
When rolling without advantage or disadvantage, the probability of getting
a 20 or a 1 is 1/20 = 0.05 or 5%. The probability of not getting a twenty 1 -
1/20 = 19/20.
</p>
<div class="chunk" id="unnamed-chunk-1">
<div class="rcode">
<div class="source">
<pre class="knitr r"><span class="hl std">P20</span> <span class="hl kwb"><-</span> <span class="hl num">0.05</span>
<span class="hl std">Pnot20</span> <span class="hl kwb"><-</span> <span class="hl num">1</span><span class="hl opt">-</span><span class="hl std">P20</span>
</pre>
</div>
</div>
</div>
<p>
When you have advantage, to not get a 20, you have to not roll a 20 twice.
The probability of getting a 20 is (1 - the proability of not getting a 20 twice), and as you can see below,
is almost 10%. Advantage about doubles you chance of getting a 20. But you
probability guessed that since you're rolling twice. :) The chance of getting a 1 is 0.05^2=0.0025, almost 0!
</p>
<div class="chunk" id="unnamed-chunk-2">
<div class="rcode">
<div class="source">
<pre class="knitr r"><span class="hl std">Pnot20</span><span class="hl opt">*</span><span class="hl std">Pnot20</span>
</pre>
</div>
<div class="output">
<pre class="knitr r">
## [1] 0.9025
</pre>
</div>
<div class="source">
<pre class="knitr r"><span class="hl num">1</span> <span class="hl opt">-</span> <span class="hl std">Pnot20</span><span class="hl opt">*</span><span class="hl std">Pnot20</span>
</pre>
</div>
<div class="output">
<pre class="knitr r">
## [1] 0.0975
</pre>
</div>
</div>
</div>
<p>
How will advantage and disadvantage affect the average roll? This time
let's estimate the answer using a simulation. A 1000 trials will give a
good estimate.
</p>
<div class="chunk" id="unnamed-chunk-3">
<div class="rcode">
<div class="source">
<pre class="knitr r"><span class="hl std">Trials</span> <span class="hl kwb"><-</span><span class="hl num">10000</span>
</pre>
</div>
</div>
</div>
<p>Below we simulate rolling two D20s</p>
<div class="chunk" id="unnamed-chunk-4">
<div class="rcode">
<div class="source">
<pre class="knitr r"><span class="hl kwd">set.seed</span><span class="hl std">(</span><span class="hl num">1</span><span class="hl std">)</span>
<span class="hl std">x</span> <span class="hl kwb"><-</span> <span class="hl kwd">sample</span><span class="hl std">(</span><span class="hl num">1</span><span class="hl opt">:</span><span class="hl num">20</span><span class="hl std">,Trials,</span><span class="hl kwc">replace</span> <span class="hl std">=</span> <span class="hl num">TRUE</span><span class="hl std">)</span>
<span class="hl std">y</span> <span class="hl kwb"><-</span> <span class="hl kwd">sample</span><span class="hl std">(</span><span class="hl num">1</span><span class="hl opt">:</span><span class="hl num">20</span><span class="hl std">,Trials,</span><span class="hl kwc">replace</span> <span class="hl std">=</span> <span class="hl num">TRUE</span><span class="hl std">)</span>
</pre>
</div>
</div>
</div>
<p>
When we have advantage, roll two D20 and pick the max. With disadvantage
roll two D20 and pick the min.
</p>
<div class="chunk" id="unnamed-chunk-5">
<div class="rcode">
<div class="source">
<pre class="knitr r"><span class="hl std">RollsWithAdvantage</span> <span class="hl kwb"><-</span> <span class="hl kwd">apply</span><span class="hl std">(</span><span class="hl kwd">cbind</span><span class="hl std">(x,y),</span> <span class="hl num">1</span><span class="hl std">,max)</span>
<span class="hl std">Advantage</span> <span class="hl kwb"><-</span><span class="hl kwd">mean</span><span class="hl std">(RollsWithAdvantage)</span>
<span class="hl std">RollsWithDisadvantage</span> <span class="hl kwb"><-</span> <span class="hl kwd">apply</span><span class="hl std">(</span><span class="hl kwd">cbind</span><span class="hl std">(x,y),</span> <span class="hl num">1</span><span class="hl std">,min)</span>
<span class="hl std">Disadvantage</span> <span class="hl kwb"><-</span><span class="hl kwd">mean</span><span class="hl std">(RollsWithDisadvantage)</span>
</pre>
</div>
</div>
</div>
<p>
Below is the simulated advantage mean, the calculated mean of a regular
D20 roll, and the simulated disadvantage mean. Advantage adds about 3.4
points and disadvantage subtracts 3.4 or about +/- 17%.
</p>
<div class="chunk" id="unnamed-chunk-6">
<div class="rcode">
<div class="source">
<pre class="knitr r"><span class="hl std">Advantage</span>
</pre>
</div>
<div class="output">
<pre class="knitr r">
## [1] 13.8863
</pre>
</div>
<div class="source">
<pre class="knitr r"><span class="hl kwd">mean</span><span class="hl std">(</span><span class="hl num">1</span><span class="hl opt">:</span><span class="hl num">20</span><span class="hl std">)</span>
</pre>
</div>
<div class="output">
<pre class="knitr r">
## [1] 10.5
</pre>
</div>
<div class="source">
<pre class="knitr r"><span class="hl std">Disadvantage</span>
</pre>
</div>
<div class="output">
<pre class="knitr r">
## [1] 7.1287
</pre>
</div>
</div>
</div>
<p>
Below are 3 figures with histograms. The first histogram is made with 10,000 simulated D20 advantage rolls. Advantage moves a lot of probability to the right. The second histogram is made using 10,000 simulated D20 regular rolls. The histogram is approximately flat and each number is at about 0.05% and that's what we calculated. The last histogram is made with 10,000 simulated D20 disadvantage rolls. Disadvantage moves an equal amount to the left.
</p>
<!-- ******** PLOT Estimated Density of ADV ********** -->
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgbYsBIZC02EgAO3ff_-PjM9JfoZIX2QNRBhpXh5dB8eGt3bzpeUp0bO2MM8Qco4vTQO-kMOEqa8cngf2TSwN3HkWUuAZQV9XCmTFAbh9ves5N-F8RZqAtcHi9P9tDoDMuWvvmcmW6ePgiZ/s1600/HistOfAdvantage.png" imageanchor="1" ><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgbYsBIZC02EgAO3ff_-PjM9JfoZIX2QNRBhpXh5dB8eGt3bzpeUp0bO2MM8Qco4vTQO-kMOEqa8cngf2TSwN3HkWUuAZQV9XCmTFAbh9ves5N-F8RZqAtcHi9P9tDoDMuWvvmcmW6ePgiZ/s400/HistOfAdvantage.png" width="400" height="300" data-original-width="800" data-original-height="600" /></a>
<!-- ******** PLOT Estimated Density of x ********** -->
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEh1Mkl7t_gkokp_5Es0GZ8KuiKtzZsB8KF5zIhrDxxEBqIOtlcLWBiiIHrY6vOU-11UPAZhlk5ovcSiaY3IlOh0eOeNhLinR0ZdV2cQTuN5UpI9DqrhRZtQTl3Cn8RjnLe2ZVLofxLsAjvb/s1600/HistOfx.png" imageanchor="1" ><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEh1Mkl7t_gkokp_5Es0GZ8KuiKtzZsB8KF5zIhrDxxEBqIOtlcLWBiiIHrY6vOU-11UPAZhlk5ovcSiaY3IlOh0eOeNhLinR0ZdV2cQTuN5UpI9DqrhRZtQTl3Cn8RjnLe2ZVLofxLsAjvb/s400/HistOfx.png" width="400" height="300" data-original-width="800" data-original-height="600" /></a>
<!-- ******** PLOT Estimated Density of DIS_ADV ********** -->
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEijXPeXsztbucgNWh6I5cWGm3nai_7TxJMIUoLvBDeVuAtMiA7bX49h9nxC-vimB0gcC3P_reOLeHYwkSRwA8M65B9WXYguUERbKeQUnfLeZKah_MmXa9qJ3d84ryrsjDwYpfDfvbX0FySw/s1600/HistOfDisadvantage.png" imageanchor="1" ><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEijXPeXsztbucgNWh6I5cWGm3nai_7TxJMIUoLvBDeVuAtMiA7bX49h9nxC-vimB0gcC3P_reOLeHYwkSRwA8M65B9WXYguUERbKeQUnfLeZKah_MmXa9qJ3d84ryrsjDwYpfDfvbX0FySw/s400/HistOfDisadvantage.png" width="400" height="300" data-original-width="800" data-original-height="600" /></a>
<!-- ******** ********** -->
<p>
The probability of getting 11 or better with advantage is 1- probability
of getting 10 or less twice. Rolling with advantage moves a 50% chance
to a 75% chance!
</p>
<div class="chunk" id="unnamed-chunk-10">
<div class="rcode">
<div class="source">
<pre class="knitr r"><span class="hl std">p10orLess</span> <span class="hl kwb"><-</span> <span class="hl num">0.5</span>
<span class="hl num">1</span> <span class="hl opt">-</span> <span class="hl std">p10orLess</span><span class="hl opt">^</span><span class="hl num">2</span>
</pre>
</div>
<div class="output">
<pre class="knitr r">
## [1] 0.75
</pre>
</div>
</div>
</div>
</body>
</html>
Chris Carbonehttp://www.blogger.com/profile/10458001230470047172noreply@blogger.com1tag:blogger.com,1999:blog-7566945490770798840.post-50940818102560336552020-06-20T08:03:00.001-07:002020-06-20T08:03:41.167-07:00Complex Normal Samples In R<html>
<head>
<style type="text/css">
.inline {
background-color: #f7f7f7;
border:solid 1px #B0B0B0;
}
.error {
font-weight: bold;
color: #FF0000;
}
.warning {
font-weight: bold;
}
.message {
font-style: italic;
}
.source, .output, .warning, .error, .message {
padding: 0 1em;
border:solid 1px #F7F7F7;
}
.source {
background-color: #f5f5f5;
}
.left {
text-align: left;
}
.right {
text-align: right;
}
.center {
text-align: center;
}
.hl.num {
color: #AF0F91;
}
.hl.str {
color: #317ECC;
}
.hl.com {
color: #AD95AF;
font-style: italic;
}
.hl.opt {
color: #000000;
}
.hl.std {
color: #585858;
}
.hl.kwa {
color: #295F94;
font-weight: bold;
}
.hl.kwb {
color: #B05A65;
}
.hl.kwc {
color: #55aa55;
}
.hl.kwd {
color: #BC5A65;
font-weight: bold;
}
</style>
</head>
<body>
<h4>Normal Samples</h4>
<p>If we want 10 samples from a Gaussian or normal random process with variance 4 can use rnorm(10,sd=2). Remember the standard deviation (sd) is the square root of the variance.</p>
<div class="chunk" id="unnamed-chunk-1"><div class="rcode"><div class="source"><pre class="knitr r"><span class="hl kwd">set.seed</span><span class="hl std">(</span><span class="hl num">2</span><span class="hl std">)</span>
<span class="hl std">x</span> <span class="hl kwb"><-</span> <span class="hl kwd">rnorm</span><span class="hl std">(</span><span class="hl num">10</span><span class="hl std">,</span><span class="hl kwc">sd</span><span class="hl std">=</span><span class="hl num">2</span><span class="hl std">)</span>
<span class="hl std">x</span>
</pre></div>
<div class="output"><pre class="knitr r">## [1] -1.7938291 0.3696984 3.1756907 -2.2607513 -0.1605035 0.2648406
## [7] 1.4159095 -0.4793960 3.9689479 -0.2775740
</pre></div>
<div class="source"><pre class="knitr r"><span class="hl kwd">var</span><span class="hl std">(x)</span>
</pre></div>
<div class="output"><pre class="knitr r">## [1] 3.880803
</pre></div>
</div></div>
<p>The var() function produces an estimate of the variance, if we want a better estimate we need more samples.</p>
<div class="chunk" id="unnamed-chunk-2"><div class="rcode"><div class="source"><pre class="knitr r"><span class="hl kwd">var</span><span class="hl std">(</span><span class="hl kwd">rnorm</span><span class="hl std">(</span><span class="hl num">1000</span><span class="hl std">,</span><span class="hl kwc">sd</span><span class="hl std">=</span><span class="hl num">2</span><span class="hl std">))</span>
</pre></div>
<div class="output"><pre class="knitr r">## [1] 4.105966
</pre></div>
</div></div>
<h4>Complex Normal Samples</h4>
<p> If we are using base R and want complex normal (CN) samples, we need to write our own function. When the signal processing literature refers to CN they are usually referring to circularly-symmetric CN. Circularly-symmetric means the samples are independent and their mean is 0.</p>
<p>The function produces N CN samples with variance v. The real and imaginary parts are independent, because they are produced by different calls to rnorm(). Let x,y be independent. The var(ax) = a^2var(x) and var(x+y)=var(x)+var(y). So, if we want a variance of 1 would have to start a variance of sqrt(1/2).</p>
<div class="chunk" id="unnamed-chunk-3"><div class="rcode"><div class="source"><pre class="knitr r"><span class="hl std">makeCN</span> <span class="hl kwb"><-</span> <span class="hl kwa">function</span><span class="hl std">(</span><span class="hl kwc">N</span><span class="hl std">,</span><span class="hl kwc">v</span><span class="hl std">=</span><span class="hl num">1</span><span class="hl std">) {(</span><span class="hl kwd">sqrt</span><span class="hl std">(v</span><span class="hl opt">/</span><span class="hl num">2</span><span class="hl std">))</span><span class="hl opt">*</span><span class="hl kwd">rnorm</span><span class="hl std">(N)</span> <span class="hl opt">+</span> <span class="hl std">(</span><span class="hl kwd">sqrt</span><span class="hl std">(v</span><span class="hl opt">/</span><span class="hl num">2</span><span class="hl std">))</span><span class="hl opt">*</span><span class="hl num">1i</span><span class="hl opt">*</span><span class="hl kwd">rnorm</span><span class="hl std">(N)}</span>
<span class="hl kwd">makeCN</span><span class="hl std">(</span><span class="hl num">10</span><span class="hl std">)</span>
</pre></div>
<div class="output"><pre class="knitr r">## [1] 0.0023376-0.2079938i 0.7613032+0.6053620i 0.3946671-0.4049715i
## [4] 0.4892950-0.1207824i 0.4651165-0.2871364i -0.2312504+0.9408834i
## [7] -0.2153405-0.9648887i -1.0994866+1.0119199i 1.0396552+0.7824796i
## [10] 0.1147878+0.9059002i
</pre></div>
</div></div>
<p>If we want to check the variance, we can't use var() directly.</p>
<div class="chunk" id="unnamed-chunk-4"><div class="rcode"><div class="source"><pre class="knitr r"><span class="hl kwd">var</span><span class="hl std">(</span><span class="hl kwd">makeCN</span><span class="hl std">(</span><span class="hl num">10</span><span class="hl std">))</span>
</pre></div>
<div class="warning"><pre class="knitr r">## Warning in var(makeCN(10)): imaginary parts discarded in coercion
</pre></div>
<div class="output"><pre class="knitr r">## [1] 0.6039204
</pre></div>
</div></div>
<p>But the real and imaginary parts are independent, so we can calculate the variance separately.</p>
<div class="chunk" id="unnamed-chunk-5"><div class="rcode"><div class="source"><pre class="knitr r"><span class="hl std">z</span> <span class="hl kwb"><-</span> <span class="hl kwd">makeCN</span><span class="hl std">(</span><span class="hl num">10</span><span class="hl std">)</span>
<span class="hl kwd">var</span><span class="hl std">(</span><span class="hl kwd">Re</span><span class="hl std">(z))</span> <span class="hl opt">+</span> <span class="hl kwd">var</span><span class="hl std">(</span><span class="hl kwd">Im</span><span class="hl std">(z))</span>
</pre></div>
<div class="output"><pre class="knitr r">## [1] 0.5623849
</pre></div>
</div></div>
<p>To make this easier, we can create a function to find the variance.</p>
<div class="chunk" id="unnamed-chunk-6"><div class="rcode"><div class="source"><pre class="knitr r"><span class="hl std">varComplex</span> <span class="hl kwb"><-</span> <span class="hl kwa">function</span><span class="hl std">(</span><span class="hl kwc">z</span><span class="hl std">)</span> <span class="hl kwd">var</span><span class="hl std">(</span><span class="hl kwd">Re</span><span class="hl std">(z))</span> <span class="hl opt">+</span> <span class="hl kwd">var</span><span class="hl std">(</span><span class="hl kwd">Im</span><span class="hl std">(z))</span>
</pre></div>
</div></div>
<p>To get a good estimate we need a-lot of samples.</p>
<div class="chunk" id="unnamed-chunk-7"><div class="rcode"><div class="source"><pre class="knitr r"><span class="hl kwd">varComplex</span><span class="hl std">(</span><span class="hl kwd">makeCN</span><span class="hl std">(</span><span class="hl num">1000</span><span class="hl std">))</span>
</pre></div>
<div class="output"><pre class="knitr r">## [1] 1.016615
</pre></div>
</div></div>
<p>Let's set the variance to 2 and then estimate the variance of the samples.</p>
<div class="chunk" id="unnamed-chunk-8"><div class="rcode"><div class="source"><pre class="knitr r"><span class="hl kwd">varComplex</span><span class="hl std">(</span><span class="hl kwd">makeCN</span><span class="hl std">(</span><span class="hl num">1000</span><span class="hl std">,</span><span class="hl kwc">v</span><span class="hl std">=</span><span class="hl num">2</span><span class="hl std">))</span>
</pre></div>
<div class="output"><pre class="knitr r">## [1] 1.925119
</pre></div>
</div></div>
<p>Success!</p>
</body>
</html>
Chris Carbonehttp://www.blogger.com/profile/10458001230470047172noreply@blogger.com0tag:blogger.com,1999:blog-7566945490770798840.post-78551428361413792642020-05-12T18:47:00.000-07:002020-05-12T18:47:10.994-07:00Book Review: Introductory Time Series with RI'm a big fan of R and time series analysis, so I was excited to read the book "Introductory Time Series with R. I've been using the book for about 9 years, so I thought it was about time for a review! In this review, I'm going to cover the following topics: the amount of R content, the subject content, who is the book for, and my overall recommendation.
<br />
<br />
R Content<br />
The R content is high. All the ideas in the book are heavily illustrated with R code. At the beginning of the book, the authors point out that they use Sweave to embed the code and plots. They also make use of online data sets, so you can type in examples from the book and reproduce the calculations and figures.
<br />
<br />
Time Series Content<br />
The book covers all the time series topics you'd want in an introduction, plus a few specialty topics like multivariate models. Each chapter is a solid introduction to a topic in time series analysis.
<br />
<br />
Who is the Book For<br />
1) People who what to learn time series analysis. It covers the theory, application, and has plenty of opportunities for hands-on learning. 2) People who want to learn to do time series analysis using R. That's why I bought the book! I'm familiar with time series analysis and I bought the book to learn how to do time series stuff with R. The book definitely delivered on that account. 3) People who want to teach a course on time series analysis. The book has plenty of examples and exercises (not sure if there is a solution manual).
<br />
<br />
Overall Recommendation<br />
The book is well written. The R code is clear and well presented. The figures are numerous and informative. The book is wonderful and I highly recommend it!
<br />
<br />
<a href="https://www.amazon.com/gp/product/0387886974/ref=ppx_yo_dt_b_search_asin_title?ie=UTF8&psc=1">https://www.amazon.com/gp/product/0387886974/ref=ppx_yo_dt_b_search_asin_title?ie=UTF8&psc=1</a>Chris Carbonehttp://www.blogger.com/profile/10458001230470047172noreply@blogger.com1tag:blogger.com,1999:blog-7566945490770798840.post-36516431148799865612019-07-29T18:16:00.000-07:002020-05-12T17:07:21.147-07:00Book Review: Advanced RWhen I first started using R, my code looked just like MATLAB code. I read an intro to R programming, did some data analysis using R, and I was using the arrow assignment operator! But otherwise, my R code looked very much like MATLAB. Then I read Advanced R by Hadley Wickham. I found out, while I wasn't doing anything wrong, I wasn't getting the most out of R.<br />
<br />
Advanced R is very readable if you've used R and maybe some other programming experience. If you're self-taught in R and have no other programming experience, you may find the book a bit tough. The book does a great job describing the important features and how to make the most of them. I particularly liked the chapter on functions and the section on functional programming.<br />
<br />
I highly recommend it!<br />
<br />
<a href="https://www.amazon.com/dp/1466586966/ref=cm_sw_su_dp?tag=devtools-20">https://www.amazon.com/dp/1466586966/ref=cm_sw_su_dp?tag=devtools-20</a><br />
<br />
<a href="https://www.amazon.com/dp/1466586966/ref=cm_sw_su_dp?tag=devtools-20"></a> Update: A commenters pointed out there is a new version of the book out. The link is below.<br />
<br />
<br />
<a href="https://www.amazon.com/Advanced-Second-Chapman-Hall-CRC/dp/0815384572/ref=pd_lpo_sbs_14_t_0?_encoding=UTF8&psc=1&refRID=11FT84CEP3CE0YQ6XYBJ">https://www.amazon.com/Advanced-Second-Chapman-Hall-CRC/dp/0815384572/ref=pd_lpo_sbs_14_t_0?_encoding=UTF8&psc=1&refRID=11FT84CEP3CE0YQ6XYBJ</a>Chris Carbonehttp://www.blogger.com/profile/10458001230470047172noreply@blogger.com5tag:blogger.com,1999:blog-7566945490770798840.post-42436006863931035792019-06-26T16:58:00.001-07:002019-06-27T22:20:06.269-07:00Blackman-Tukey Spectral Estimator in R
<html>
<head>
<style type="text/css">
.inline {
background-color: #f7f7f7;
border:solid 1px #B0B0B0;
}
.error {
font-weight: bold;
color: #FF0000;
}
.warning {
font-weight: bold;
}
.message {
font-style: italic;
}
.source, .output, .warning, .error, .message {
padding: 0 1em;
border:solid 1px #F7F7F7;
}
.source {
background-color: #f5f5f5;
}
.left {
text-align: left;
}
.right {
text-align: right;
}
.center {
text-align: center;
}
.hl.num {
color: #AF0F91;
}
.hl.str {
color: #317ECC;
}
.hl.com {
color: #AD95AF;
font-style: italic;
}
.hl.opt {
color: #000000;
}
.hl.std {
color: #585858;
}
.hl.kwa {
color: #295F94;
font-weight: bold;
}
.hl.kwb {
color: #B05A65;
}
.hl.kwc {
color: #55aa55;
}
.hl.kwd {
color: #BC5A65;
font-weight: bold;
}
</style>
<title>Blackman-Tukey Spectral Estimator in R!</title>
</head>
<body>
<p> </p>
<p>There are two definitions of the power spectral density (PSD). Both definitions are mathematically nearly identical and define a function that describes the distribution of power over the frequency components in our data set. The periodogram PSD estimator is based on the first definition of the PSD (see periodogram post). The Blackman-Tukey spectral estimator (BTSE) is based on the second definition. The second definition says, find the PSD by calculating the Fourier transform of the autocorrelation sequence (ACS). In R this definition is written as </p>
<div class="chunk" id="unnamed-chunk-1"><div class="rcode"><div class="source"><pre class="knitr r"><span class="hl std">PSD</span> <span class="hl kwb"><-</span> <span class="hl kwa">function</span><span class="hl std">(</span><span class="hl kwc">rxx</span><span class="hl std">) {</span>
<span class="hl kwd">fft</span><span class="hl std">(rxx)</span>
<span class="hl std">}</span>
</pre></div>
</div></div>
<p>where fft is the R implementation of the fast Fourier transform, rxx is the autocorrelation sequence (ACS), the k'th element of the ACS rxx[k] = E[x[0]x[k]], k -infinity to +infinity, and E is the expectation operator. The xx in rxx[k] is a reminder r is a correlation between x and itself. The rxx[k]s are sometimes called lags. The ACS has the propriety that rxx[-k]=rxx[k]*, where * is the complex conjugate. In the post, we will only use real numbers and I'll drop the * from here forward.</p>
<p>So, to find the PSD we just calculate rxx and take its fft! Unfortunately, in practice, we cannot do this. Calculating the expected value requires the probability density function (PDF) of x, which we don't know and we need an infinite amount of data, which we don't have. So, we can't calculate the PSD: we're doomed!</p>
<p>No, we are not doomed. We can't calculate the PSD, but we estimate it! We can derive an estimator for the PSD from the definition of the PSD. First, we replace rxx with an estimate of rxx. We replace the expected value, which gives the exact rxx, with an average, which gives us an estimate of rxx. The E[x[0]x[k]] is replaced with (1/N)(x[1]x[1+k]+x[2]x[2+k]+...+x[N-1-k]x[N-1]+x[N-k]x[N]), where N is the number of data samples. For example; if k=0, then rxx[k]=(1/N)*sum(x*x). In R code the estimate is written as</p>
<div class="chunk" id="unnamed-chunk-2"><div class="rcode"><div class="source"><pre class="knitr r"><span class="hl std">lagEstimate</span> <span class="hl kwb"><-</span><span class="hl kwa">function</span><span class="hl std">(</span><span class="hl kwc">x</span><span class="hl std">,</span><span class="hl kwc">k</span><span class="hl std">,</span><span class="hl kwc">N</span><span class="hl std">=</span><span class="hl kwd">length</span><span class="hl std">(x)){</span>
<span class="hl std">(</span><span class="hl num">1</span><span class="hl opt">/</span><span class="hl std">N)</span><span class="hl opt">*</span><span class="hl kwd">sum</span><span class="hl std">(x[</span><span class="hl num">1</span><span class="hl opt">:</span><span class="hl std">(N</span><span class="hl opt">-</span><span class="hl std">k)]</span><span class="hl opt">*</span><span class="hl std">x[(k</span><span class="hl opt">+</span><span class="hl num">1</span><span class="hl std">)</span><span class="hl opt">:</span><span class="hl std">N])</span>
<span class="hl std">}</span>
</pre></div>
</div></div>
<p>If we had an infinite amount of data, N=infinity, we could use lagEstimate to estimate the entire infinite ACS. Unfortunately we don't have an infinite amount of data, even if we did it wouldn't fit into a computer. So, we can only estimate a finite amount of ASC elements. The function below calculates lags 0 to kMax.</p>
<div class="chunk" id="unnamed-chunk-3"><div class="rcode"><div class="source"><pre class="knitr r"><span class="hl std">Lags</span> <span class="hl kwb"><-</span><span class="hl kwa">function</span><span class="hl std">(</span><span class="hl kwc">x</span><span class="hl std">,</span><span class="hl kwc">kMax</span><span class="hl std">) {</span>
<span class="hl kwd">sapply</span><span class="hl std">(</span><span class="hl num">0</span><span class="hl opt">:</span><span class="hl std">kMax,lagEstimate,</span><span class="hl kwc">x</span><span class="hl std">=x)</span>
<span class="hl std">}</span>
</pre></div>
</div></div>
<p>Before we can try these functions out we need data. In this case the data came from a random process with the PSD plotted in the figure below. The x-axis is normalized frequency(frequency divided by the sampling rate). So, if the sampling rate was 1000 Hz, you could multiply the normalized frequency by 1000 Hz and then the frequency axis would read 0 Hz to 1000 Hz. The y-axis in in dB (10log10(amplitude)). You can see six large sharp peaks in the plot and a gradual dip towards 0 Hz and then back up. Some of the peaks are close together and will be hard to resolve.</p>
<!-- ******** PLOT PsdPlotNumerical ********** -->
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhgZ-GsG_RzL2CnZcP56C7MCw7MclQnAJi1WBH2HMaoSAFxhXb1vEE3ahzUoZysVzRVj13jJLuBonpIVYKr8h6anLaTRXB8s5EaZgwtLS257Wcj_FBKyjSSd7XL8zaH8Ii35T_a8we52Lej/s1600/PsdPlotNumerical.png" imageanchor="1" ><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhgZ-GsG_RzL2CnZcP56C7MCw7MclQnAJi1WBH2HMaoSAFxhXb1vEE3ahzUoZysVzRVj13jJLuBonpIVYKr8h6anLaTRXB8s5EaZgwtLS257Wcj_FBKyjSSd7XL8zaH8Ii35T_a8we52Lej/s400/PsdPlotNumerical.png" width="400" height="300" data-original-width="800" data-original-height="600" /></a>
<p>The data produced by the random process is plotted below. This is the data we will use through this post. The data is posted here: <a href="https://github.com/carbone1853/GetYourDataOn">x.RData</a></p>
<!-- ******** PLOT PlotOfData ********** -->
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEiTzna3BXYYhe5Ng70AuuyldIxZE4LDa2QUSFiqWvvWYkU_0TTpXyGSGXU4AYXbS07w-3MXtzop6bHYVIJKZF9zhkjL1OYV7X0TBzNK3iyYZQvECyvdF3uQGiVrjZQFLx4vrcw7OUJFiBXS/s1600/PlotOfData.png" imageanchor="1" ><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEiTzna3BXYYhe5Ng70AuuyldIxZE4LDa2QUSFiqWvvWYkU_0TTpXyGSGXU4AYXbS07w-3MXtzop6bHYVIJKZF9zhkjL1OYV7X0TBzNK3iyYZQvECyvdF3uQGiVrjZQFLx4vrcw7OUJFiBXS/s400/PlotOfData.png" width="400" height="300" data-original-width="800" data-original-height="600" /></a>
<p>Let's calculate the the ACS up to the 5th lag using the data.</p>
<div class="chunk" id="unnamed-chunk-6"><div class="rcode"><div class="source"><pre class="knitr r"><span class="hl kwd">Lags</span><span class="hl std">(x,</span><span class="hl kwc">kMax</span><span class="hl std">=</span><span class="hl num">5</span><span class="hl std">)</span>
</pre></div>
<div class="output"><pre class="knitr r">## [1] 6.095786 -1.368336 3.341608 1.738122 -1.737459 3.651765
</pre></div>
</div></div>
<p>A kMax of 5 gives us 6 lags: {r[0], r[1], r[2], r[3], r[4], r[5]}. These 6 lags are not an ACS, but are part of an ACS.</p>
<p>We used Lags to estimate the positive lags up to kMax, but the ACS is even sequence, r[-k]=r[k] for all k. So, let's write a function to make a sequence consisting of lags from r[-kMax] to r[kMax]. This is a windowed ACS, values outside of the +/- kMax are replaced with 0. Where it won't cause confusion, I'll refer to the windowed ACS, as the ACS.</p>
<div class="chunk" id="unnamed-chunk-7"><div class="rcode"><div class="source"><pre class="knitr r"><span class="hl std">acsWindowed</span> <span class="hl kwb"><-</span> <span class="hl kwa">function</span><span class="hl std">(</span><span class="hl kwc">x</span><span class="hl std">,</span><span class="hl kwc">kMax</span><span class="hl std">,</span><span class="hl kwc">Nzero</span><span class="hl std">=</span><span class="hl num">0</span><span class="hl std">){</span>
<span class="hl std">rHalf</span> <span class="hl kwb"><-</span> <span class="hl kwd">c</span><span class="hl std">(</span><span class="hl kwd">Lags</span><span class="hl std">(x,kMax),</span><span class="hl kwd">rep</span><span class="hl std">(</span><span class="hl num">0</span><span class="hl std">,Nzero))</span>
<span class="hl kwd">c</span><span class="hl std">(</span><span class="hl kwd">rev</span><span class="hl std">(rHalf[</span><span class="hl num">2</span><span class="hl opt">:</span><span class="hl kwd">length</span><span class="hl std">(rHalf)]),rHalf)</span>
<span class="hl std">}</span>
</pre></div>
</div></div>
<p>Let's try this function out.</p>
<div class="chunk" id="unnamed-chunk-8"><div class="rcode"><div class="source"><pre class="knitr r"><span class="hl std">acsW</span> <span class="hl kwb"><-</span> <span class="hl kwd">acsWindowed</span><span class="hl std">(x,</span><span class="hl num">9</span><span class="hl std">)</span>
</pre></div>
</div></div>
<p>In the figure below you can see the r[0] lag, the maximum, is plotted in the middle of the plot.</p>
<!-- ******** PLOT acsW ********** -->
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhijAF6AGNiN9vYUZsXlMpeFz908HirxGBtr0s0nZnHTkVA9Gg8Svoik7QrKDDActOxq4pdgb0tAA7JpZnOIkSuQyHexeFvDr6FTgfOLUFyV5Q-qXoLKH8BTDn3hVZ1k6dFDnONf_Humm7K/s1600/acsW.png" imageanchor="1" ><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhijAF6AGNiN9vYUZsXlMpeFz908HirxGBtr0s0nZnHTkVA9Gg8Svoik7QrKDDActOxq4pdgb0tAA7JpZnOIkSuQyHexeFvDr6FTgfOLUFyV5Q-qXoLKH8BTDn3hVZ1k6dFDnONf_Humm7K/s400/acsW.png" width="400" height="300" data-original-width="800" data-original-height="600" /></a>
<p>The ACS in the figure above is how the ACS is usually plotted in textbooks. In textbooks the sum in the Fourier transform ranges from -N/2 to (N-1)/2. So, the r[0] lag should be in the center of the plot. In R the sum in the Fourier transform ranges from 1 to N, so the 0'th lag has to be first. We could just make the sequence in R form, but it is often handy to start in textbook from and switch to R form. We can write a function to make switching from textbook to R easy.</p>
<div class="chunk" id="unnamed-chunk-10"><div class="rcode"><div class="source"><pre class="knitr r"><span class="hl std">Textbook2R</span> <span class="hl kwb"><-</span> <span class="hl kwa">function</span><span class="hl std">(</span><span class="hl kwc">x</span><span class="hl std">,</span><span class="hl kwc">N</span><span class="hl std">=</span><span class="hl kwd">length</span><span class="hl std">(x),</span><span class="hl kwc">foldN</span><span class="hl std">=</span><span class="hl kwd">ceiling</span><span class="hl std">(N</span><span class="hl opt">/</span><span class="hl num">2</span><span class="hl std">)) {</span>
<span class="hl kwd">c</span><span class="hl std">(x[foldN</span><span class="hl opt">:</span><span class="hl std">N],x[</span><span class="hl num">1</span><span class="hl opt">:</span><span class="hl std">(foldN</span><span class="hl opt">-</span><span class="hl num">1</span><span class="hl std">)])</span>
<span class="hl std">}</span>
</pre></div>
</div></div>
<p>Notice in the figure below the maximum lag r[0], is plotted at the beginning.</p>
<!-- ******** PLOT acsWtextbook ********** -->
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgcQTCMKb9lFQZu5giUs8VZC2gduJd_bwYqbJQ77GzJx8q6oo-XFmTpWo2i3GAa04VJvHyP07Sj1bpM3e5Mkk0r2baaao72l0GsB5s32yghDGZE-2kRrDPpzHtqXa6JpmjUAmgIehs-iE4h/s1600/acsWtextbook.png" imageanchor="1" ><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgcQTCMKb9lFQZu5giUs8VZC2gduJd_bwYqbJQ77GzJx8q6oo-XFmTpWo2i3GAa04VJvHyP07Sj1bpM3e5Mkk0r2baaao72l0GsB5s32yghDGZE-2kRrDPpzHtqXa6JpmjUAmgIehs-iE4h/s400/acsWtextbook.png" width="400" height="300" data-original-width="800" data-original-height="600" /></a>
<p>Let's imagine we have an infinite amount of data and used it to estimated and infinite number of ACS lags. Let's call that sequence rAll. We make a windowed ACS by setting rW=rAll*W, where W=1 for our 9 lags and 0 everywhere else. W is called the rectangular window, because, as you can see in the plot below, it's plot looks like a rectangle. By default when we estimate a finite number of lags we are using a rectangular window.</p>
<div class="chunk" id="unnamed-chunk-12"><div class="rcode"><div class="source"><pre class="knitr r"><span class="hl std">W</span> <span class="hl kwb"><-</span> <span class="hl kwd">c</span><span class="hl std">(</span><span class="hl kwd">rep</span><span class="hl std">(</span><span class="hl num">0</span><span class="hl std">,</span><span class="hl num">9</span><span class="hl std">),</span><span class="hl kwd">rep</span><span class="hl std">(</span><span class="hl num">1</span><span class="hl std">,</span><span class="hl num">9</span><span class="hl std">),</span><span class="hl kwd">rep</span><span class="hl std">(</span><span class="hl num">0</span><span class="hl std">,</span><span class="hl num">9</span><span class="hl std">))</span>
</pre></div>
</div></div>
<!-- ******** PLOT W ********** -->
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEiO0GLgaHh5L0nk6w6ZAbi9ghvTFXkZid55UT6Ontgb95iVs9J60TgHnA7Aa8UGXCBMB3r1QlqOThM_Nc2zWAOy2rRmADHcY3VqWEy13XvMap3d8rGFE4cIVqhPAAeOtiQxWbydBshKCwGS/s1600/W.png" imageanchor="1" ><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEiO0GLgaHh5L0nk6w6ZAbi9ghvTFXkZid55UT6Ontgb95iVs9J60TgHnA7Aa8UGXCBMB3r1QlqOThM_Nc2zWAOy2rRmADHcY3VqWEy13XvMap3d8rGFE4cIVqhPAAeOtiQxWbydBshKCwGS/s400/W.png" width="400" height="300" data-original-width="800" data-original-height="600" /></a>
<p>The reason we can not use a rectangular window is its Fourier Transform is not always positive. As you can see in the plot below there are several values below zero, indicated with dotted line. Re() functions removes some small imaginary numbers due to numerical error, some imaginary dust we have to sweep up. </p>
<div class="chunk" id="unnamed-chunk-14"><div class="rcode"><div class="source"><pre class="knitr r"><span class="hl std">FFT_W</span> <span class="hl kwb"><-</span> <span class="hl kwd">Re</span><span class="hl std">(</span><span class="hl kwd">fft</span><span class="hl std">(</span><span class="hl kwd">Textbook2R</span><span class="hl std">(W)))</span>
</pre></div>
</div></div>
<!-- ******** PLOT FFT_W ********** -->
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgpAGSedpQ6FWfYhRXsfCguN_iJe1DFFIAlO7lyuzvD7d-AJZzHMuyBxJv2J9hF627i6BFBG4KaYJ1CmZTFOpXH4Bl2r6Bu8Nte1VdX0kq2DjnRvPO9WyuT3d6UvzMq-HASUrO0X3jK_-wm/s1600/FFT_W.png" imageanchor="1" ><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgpAGSedpQ6FWfYhRXsfCguN_iJe1DFFIAlO7lyuzvD7d-AJZzHMuyBxJv2J9hF627i6BFBG4KaYJ1CmZTFOpXH4Bl2r6Bu8Nte1VdX0kq2DjnRvPO9WyuT3d6UvzMq-HASUrO0X3jK_-wm/s400/FFT_W.png" width="400" height="300" data-original-width="800" data-original-height="600" /></a>
<p>Even though the fft of the ACS rAll is positive , the produce rAll and a rectangular window might not be positive! The Bartlett window is a simple window whos fft is positive. </p>
<div class="chunk" id="unnamed-chunk-16"><div class="rcode"><div class="source"><pre class="knitr r"><span class="hl std">BartlettWindow</span> <span class="hl kwb"><-</span> <span class="hl kwa">function</span><span class="hl std">(</span><span class="hl kwc">N</span><span class="hl std">,</span><span class="hl kwc">n</span><span class="hl std">=</span><span class="hl kwd">seq</span><span class="hl std">(</span><span class="hl num">0</span><span class="hl std">, (N</span><span class="hl opt">-</span><span class="hl num">1</span><span class="hl std">))) {</span>
<span class="hl num">1</span> <span class="hl opt">-</span> <span class="hl kwd">abs</span><span class="hl std">( (n</span><span class="hl opt">-</span><span class="hl std">(N</span><span class="hl opt">-</span><span class="hl num">1</span><span class="hl std">)</span><span class="hl opt">/</span><span class="hl num">2</span><span class="hl std">)</span> <span class="hl opt">/</span> <span class="hl std">( (N</span><span class="hl opt">-</span><span class="hl num">1</span><span class="hl std">)</span><span class="hl opt">/</span><span class="hl num">2</span><span class="hl std">) )</span>
<span class="hl std">}</span>
<span class="hl std">Wb</span> <span class="hl kwb"><-</span> <span class="hl kwd">BartlettWindow</span><span class="hl std">(</span><span class="hl num">19</span><span class="hl std">)</span>
</pre></div>
</div></div>
<p>As you can see in the plot below the Fourier transform of the Bartlett window is positive.</p>
<div class="chunk" id="unnamed-chunk-17"><div class="rcode"><div class="source"><pre class="knitr r"><span class="hl std">WbFft</span> <span class="hl kwb"><-</span> <span class="hl kwd">Re</span><span class="hl std">(</span><span class="hl kwd">fft</span><span class="hl std">(</span><span class="hl kwd">Textbook2R</span><span class="hl std">(Wb)))</span>
</pre></div>
</div></div>
<!-- ******** PLOT WbFft ********** -->
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhvYgnuuuLgZlS3DOBDUSwRYEVcyeQB1V9EV-NYjX3cG_H2KLs4faZCM-RLAih2_1CjtxkQ8JFxYd912w5i26kbac1L_eU1xOTy6C8i_9GTlbJSWVznd3LjQrqMMfOBM1iEyo4ZIhege0U2/s1600/WbFft.png" imageanchor="1" ><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhvYgnuuuLgZlS3DOBDUSwRYEVcyeQB1V9EV-NYjX3cG_H2KLs4faZCM-RLAih2_1CjtxkQ8JFxYd912w5i26kbac1L_eU1xOTy6C8i_9GTlbJSWVznd3LjQrqMMfOBM1iEyo4ZIhege0U2/s400/WbFft.png" width="400" height="300" data-original-width="800" data-original-height="600" /></a>
<h2>Calculating the BTSE with R </h2>
<p>Now that we can estimate the ACS and window our estimate, we are ready to estimate the PSD of our data. The BTSE is written as</p>
<div class="chunk" id="unnamed-chunk-19"><div class="rcode"><div class="source"><pre class="knitr r"><span class="hl std">Btse</span> <span class="hl kwb"><-</span> <span class="hl kwa">function</span><span class="hl std">(</span><span class="hl kwc">rHat</span><span class="hl std">,</span><span class="hl kwc">Wb</span><span class="hl std">) {</span>
<span class="hl kwd">Re</span><span class="hl std">(</span><span class="hl kwd">fft</span><span class="hl std">(rHat</span><span class="hl opt">*</span><span class="hl std">Wb))</span>
<span class="hl std">}</span>
</pre></div>
</div></div>
<p>Note the Re() is correcting for numerical error.</p>
<p>In the first example we use a 19 point ACS lag sequence. </p>
<div class="chunk" id="unnamed-chunk-20"><div class="rcode"><div class="source"><pre class="knitr r"><span class="hl std">rHat</span> <span class="hl kwb"><-</span> <span class="hl kwd">Textbook2R</span><span class="hl std">(</span><span class="hl kwd">acsWindowed</span><span class="hl std">(x,</span><span class="hl kwc">kMax</span><span class="hl std">=</span><span class="hl num">9</span><span class="hl std">))</span>
<span class="hl std">Wb</span> <span class="hl kwb"><-</span> <span class="hl kwd">Textbook2R</span><span class="hl std">(</span><span class="hl kwd">BartlettWindow</span><span class="hl std">(</span><span class="hl kwd">length</span><span class="hl std">(rHat)))</span>
<span class="hl std">Pbtse9</span> <span class="hl kwb"><-</span> <span class="hl kwd">Btse</span><span class="hl std">(rHat,Wb)</span>
</pre></div>
</div></div>
<p>In the figure below is the BTSE calculated with a maximum lag of 9. The dotted lines indicate the locations of the peaks in the PSD we are trying to estimate. The estimate with a maximum lage of only 9 produces a poor estimate.</p>
<!-- ******** PLOT Pbtse9 ********** -->
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjiK2Vsj62s87ma1qu6ZmRa10vsfE-QTyF2bxmPZ7Wa_GxsPZWMSh6K5VXBEBO8yhVodNgI_MGgLGiLy9TQKlNYMAzQaLj71h5J1pfj8ypL2wukICf00tVHGcWlnULXc6_eRZPvppRVZbj6/s1600/Pbtse9.png" imageanchor="1" ><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjiK2Vsj62s87ma1qu6ZmRa10vsfE-QTyF2bxmPZ7Wa_GxsPZWMSh6K5VXBEBO8yhVodNgI_MGgLGiLy9TQKlNYMAzQaLj71h5J1pfj8ypL2wukICf00tVHGcWlnULXc6_eRZPvppRVZbj6/s400/Pbtse9.png" width="400" height="300" data-original-width="800" data-original-height="600" /></a>
<p>We calculate a new estimate with a maximum lag of 18.</p>
<div class="chunk" id="unnamed-chunk-22"><div class="rcode"><div class="source"><pre class="knitr r"><span class="hl std">rHat</span> <span class="hl kwb"><-</span> <span class="hl kwd">Textbook2R</span><span class="hl std">(</span><span class="hl kwd">acsWindowed</span><span class="hl std">(x,</span><span class="hl kwc">kMax</span><span class="hl std">=</span><span class="hl num">18</span><span class="hl std">))</span>
<span class="hl std">Wb</span> <span class="hl kwb"><-</span> <span class="hl kwd">Textbook2R</span><span class="hl std">(</span><span class="hl kwd">BartlettWindow</span><span class="hl std">(</span><span class="hl kwd">length</span><span class="hl std">(rHat)))</span>
<span class="hl std">Pbtse18</span> <span class="hl kwb"><-</span> <span class="hl kwd">Btse</span><span class="hl std">(rHat,Wb)</span>
</pre></div>
</div></div>
<p>The next estimate is made with a maximum lag of 18. This estimate is better, the peaks around 0.4 and 0.6 are not resolved. We still need to increase the maximum lag.</p>
<!-- ******** PLOT Pbtse18 ********** -->
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhHhdIOC8sAOZX4aoPX8oqJ4A_i-Arnl7SWY6JjT1oSsyyIe_ekBee5B1Xai8JC_iVUmaZbY50JQe0zaJ3pwZFVoz4v0YDmczpojwqWcXEsPhZqBNggn0-urfJxNFTYit0_JZVVONBJrYXU/s1600/Pbtse18.png" imageanchor="1" ><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhHhdIOC8sAOZX4aoPX8oqJ4A_i-Arnl7SWY6JjT1oSsyyIe_ekBee5B1Xai8JC_iVUmaZbY50JQe0zaJ3pwZFVoz4v0YDmczpojwqWcXEsPhZqBNggn0-urfJxNFTYit0_JZVVONBJrYXU/s400/Pbtse18.png" width="400" height="300" data-original-width="800" data-original-height="600" /></a>
<p>Finally we increase the maximum lag to 65 and recalculate the estimate.</p>
<div class="chunk" id="unnamed-chunk-24"><div class="rcode"><div class="source"><pre class="knitr r"><span class="hl std">rHat</span> <span class="hl kwb"><-</span> <span class="hl kwd">Textbook2R</span><span class="hl std">(</span><span class="hl kwd">acsWindowed</span><span class="hl std">(x,</span><span class="hl kwc">kMax</span><span class="hl std">=</span><span class="hl num">65</span><span class="hl std">))</span>
<span class="hl std">Wb</span> <span class="hl kwb"><-</span> <span class="hl kwd">Textbook2R</span><span class="hl std">(</span><span class="hl kwd">BartlettWindow</span><span class="hl std">(</span><span class="hl kwd">length</span><span class="hl std">(rHat)))</span>
<span class="hl std">Pbtse65</span> <span class="hl kwb"><-</span> <span class="hl kwd">Btse</span><span class="hl std">(rHat,Wb)</span>
</pre></div>
</div></div>
<p>This finial estimate is very good. All six peaks are resolved and the location of our estimated peaks are very close to the true peak locations locations. </p>
<!-- ******** PLOT Pbtse65 ********** -->
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEh_t1U8f9Di1vSuwZk0cBm1nkleypJ5gLOqJzeB72zWB4DZNCoZDzDD7tyIKWSWoNlD3o_9lgmSW4hqVyklYQ5d48WH_MCqMwdhAwCRP6AOoZz8EPRoj-YRLCopXsPe7uu-IxhZmbAwq-eb/s1600/Pbtse65.png" imageanchor="1" ><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEh_t1U8f9Di1vSuwZk0cBm1nkleypJ5gLOqJzeB72zWB4DZNCoZDzDD7tyIKWSWoNlD3o_9lgmSW4hqVyklYQ5d48WH_MCqMwdhAwCRP6AOoZz8EPRoj-YRLCopXsPe7uu-IxhZmbAwq-eb/s400/Pbtse65.png" width="400" height="300" data-original-width="800" data-original-height="600" /></a>
<h2>Final Thoughts</h2>
<p>Could we use 500 lags in the BTSE? In this case we could, since we have a lot of data, but the higher lags get estimated with less data and therefore have more variance. Using the high variance lags will produce a higher variance estimate.</p>
<p>Are there other ways to improve the BTSE other than using more lags? Yes! There are a few other ways. For instance, we could zero pad the lags. Basically add zeros to the end of our lag sequence. This will make the fft, in the BTSE estimator, evaluate the estimate at more frequencies and we will be able to see more details in the estimated PSD.</p>
<p>Also keep in mind there are other PSD estimation methods that do better on other PSD features. For instance, if you we more interested finding deep nulls rather than peaks, a moving average PSD estimation would be better.</p>
</body>
</html>
Chris Carbonehttp://www.blogger.com/profile/10458001230470047172noreply@blogger.com3tag:blogger.com,1999:blog-7566945490770798840.post-47867937507423919492019-06-13T18:05:00.000-07:002019-10-25T18:25:05.241-07:00Periodogram with R<html>
<head>
<style type="text/css">
.inline {
background-color: #f7f7f7;
border:solid 1px #B0B0B0;
}
.error {
font-weight: bold;
color: #FF0000;
}
.warning {
font-weight: bold;
}
.message {
font-style: italic;
}
.source, .output, .warning, .error, .message {
padding: 0 1em;
border:solid 1px #F7F7F7;
}
.source {
background-color: #f5f5f5;
}
.left {
text-align: left;
}
.right {
text-align: right;
}
.center {
text-align: center;
}
.hl.num {
color: #AF0F91;
}
.hl.str {
color: #317ECC;
}
.hl.com {
color: #AD95AF;
font-style: italic;
}
.hl.opt {
color: #000000;
}
.hl.std {
color: #585858;
}
.hl.kwa {
color: #295F94;
font-weight: bold;
}
.hl.kwb {
color: #B05A65;
}
.hl.kwc {
color: #55aa55;
}
.hl.kwd {
color: #BC5A65;
font-weight: bold;
}
</style>
<title>Periodogram with R</title>
</head>
<body>
<br />
The power spectral density (PSD) is a function that describes the distribution of power over the frequency components composing our data set. If we knew the process that generated the data, we could just calculate the PSD; we would not have to estimate it. Unfortunately, in practice we won't have access to the random process, only the samples (data) produced by the process. So, we can't get the true PSD, we can only get an estimate of the true PSD. In our example below, we made the data, so we know what the true PSD should look like.<br />
The PSD is written as<br />
<!-- Define PSD EQ -->
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEiYMuih6IsqmzwXzob3B1Mf_-0SCAUnGvGT-0OrjQnOqLYHjATMDmz4y6-6uoPRq10p8XsOSOsAtcz6q5LfROkkGuZq20v5u5BeHUdAwF9dXhR-TdLFLgoD_kBQAXoZvub8NFyw7wDENVM5/s1600/PSD.png" imageanchor="1"><img border="0" data-original-height="217" data-original-width="800" height="109" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEiYMuih6IsqmzwXzob3B1Mf_-0SCAUnGvGT-0OrjQnOqLYHjATMDmz4y6-6uoPRq10p8XsOSOsAtcz6q5LfROkkGuZq20v5u5BeHUdAwF9dXhR-TdLFLgoD_kBQAXoZvub8NFyw7wDENVM5/s400/PSD.png" width="400" /></a>
where -1/2 <= f < 1/2 x is the data, e is the complex exponential, the sum of the data x multiplied by e is the Fourier transform. In English, calculate the PSD by finding the expected value of the magnitude squared of the Fourier transform as the amount of data increases to infinity. Well, that's a mouthful! The frequency f is expressed in normalized frequency. Normalized frequency is the frequency divided by the sampling rate. Normalized frequency is often used for simplicity. To convert to frequency in Hz multiply by the sampling rate.<br />
The periodogram estimator is based on the definition in equation above. The definition has a limit and an expected value. To use the definition directly we would need an infinite amount of data, which we don't have, and we would need to know the probability density function (PDF) of the data which we don't know. To use this definition to get an estimator of the PSD, we must drop the limit and we must drop the expectation operator ($E$). We also can only evaluate a finite number of frequencies. After dropping the limit and the expectation we are left with<br />
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEg4xJMG3ZY9L4P2hPzTOFef7RGJF_uftLard0VxYAkS9MOp3vtXbHMsnyVo5GM1l5m9-ULx_imKyM-34GmkLYZbG-w2NxoseKlVub5ffsI39rAK7aP3KwgStdRGcUDwqv6ynWzNTlqllo4A/s1600/Period.png" imageanchor="1"><img border="0" data-original-height="420" data-original-width="1600" height="105" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEg4xJMG3ZY9L4P2hPzTOFef7RGJF_uftLard0VxYAkS9MOp3vtXbHMsnyVo5GM1l5m9-ULx_imKyM-34GmkLYZbG-w2NxoseKlVub5ffsI39rAK7aP3KwgStdRGcUDwqv6ynWzNTlqllo4A/s400/Period.png" width="400" /></a>
<!-- Define Periodogram EQ -->
where 0 <= f < 1 and N is the number of data samples. Notice, the range of f changes from (-1/2,1/2) in the first equation to (0,1) in the second equation. Textbooks usually use the former and R uses the latter. <br />
The periodogram is very easy to implement in R, but before we do we need to simulate some data. The code below first uses the set.seed command so R will produce the same "random" numbers each time. Then it creates a 32 normally distributed numbers and 32 points of a sine wave with a normalized frequency of 0.4 and a amplitude of 2. The signal is made up of a sine wave and the random points added together.<br />
The figure below is a plot of the data generated above. To me it looks random and it is not obvious there is a sine wave in there.<br />
<div class="chunk" id="unnamed-chunk-1">
<div class="rcode">
<div class="source">
<pre class="knitr r"><span class="hl kwd">set.seed</span><span class="hl std">(</span><span class="hl num">0</span><span class="hl std">)</span>
<span class="hl std">N</span> <span class="hl kwb"><-</span><span class="hl num">32</span>
<span class="hl std">n</span><span class="hl kwb"><-</span> <span class="hl num">0</span><span class="hl opt">:</span><span class="hl std">(N</span><span class="hl opt">-</span><span class="hl num">1</span><span class="hl std">)</span>
<span class="hl std">w</span> <span class="hl kwb"><-</span> <span class="hl kwd">rnorm</span><span class="hl std">(</span><span class="hl num">1</span><span class="hl opt">:</span><span class="hl std">N)</span>
<span class="hl std">f1</span> <span class="hl kwb"><-</span> <span class="hl num">0.4</span>
<span class="hl std">A1</span> <span class="hl kwb"><-</span> <span class="hl num">2</span>
<span class="hl std">s</span> <span class="hl kwb"><-</span> <span class="hl std">A1</span><span class="hl opt">*</span><span class="hl kwd">sin</span><span class="hl std">(</span><span class="hl num">2</span><span class="hl opt">*</span><span class="hl std">pi</span><span class="hl opt">*</span><span class="hl std">f1</span><span class="hl opt">*</span><span class="hl std">n)</span>
<span class="hl std">x</span> <span class="hl kwb"><-</span> <span class="hl std">s</span> <span class="hl opt">+</span> <span class="hl std">w</span>
</pre>
</div>
</div>
</div>
Now we are going to implement the periodogram in R. It is pretty much a direct translation from the math. We start with the 1/N and multiply that by the absolute value of the square of the Fourier transform of the data. Writing that in English was more work than implementing the periodogram in R!<br />
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEiuXefGJa896qu0gVsrRkbWddnqIESPN7SiLIOpeQ0xM6xqwMRnHg-MtO366VJpyT_zODugDkJIiYcp7NuzcCMWg0K7ztVVRSC8_QzLtrhlFG57c3MwGNVyb9EHL-V-umGYZ4l6E4BvLPiI/s1600/x.png" imageanchor="1"><img border="0" data-original-height="600" data-original-width="800" height="300" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEiuXefGJa896qu0gVsrRkbWddnqIESPN7SiLIOpeQ0xM6xqwMRnHg-MtO366VJpyT_zODugDkJIiYcp7NuzcCMWg0K7ztVVRSC8_QzLtrhlFG57c3MwGNVyb9EHL-V-umGYZ4l6E4BvLPiI/s400/x.png" width="400" /></a>
<div class="chunk" id="unnamed-chunk-2">
<div class="rcode">
<div class="source">
<pre class="knitr r"><span class="hl std">xPer</span> <span class="hl kwb"><-</span> <span class="hl std">(</span><span class="hl num">1</span><span class="hl opt">/</span><span class="hl std">N)</span><span class="hl opt">*</span><span class="hl kwd">abs</span><span class="hl std">(</span><span class="hl kwd">fft</span><span class="hl std">(x)</span><span class="hl opt">^</span><span class="hl num">2</span><span class="hl std">)</span>
<span class="hl std">f</span> <span class="hl kwb"><-</span> <span class="hl kwd">seq</span><span class="hl std">(</span><span class="hl num">0</span><span class="hl std">,</span><span class="hl num">1.0</span><span class="hl opt">-</span><span class="hl num">1</span><span class="hl opt">/</span><span class="hl std">N,</span><span class="hl kwc">by</span><span class="hl std">=</span><span class="hl num">1</span><span class="hl opt">/</span><span class="hl std">N)</span>
</pre>
</div>
</div>
</div>
The figure below is a plot of the periodogram of the data. The dotted line marks the location of the frequency f1, the frequency of the sine wave of the data. Now the sine wave component really stands out!. Notice two things. First, the peaks of the periodogram seems to be a bit off the true values. Second, the plot looks jagged. Both of these things are caused by the same thing: the periodogram is only evaluated at 32 frequency bins. The frequencies we evaluate the periodogram on are called bins. We can fix both problems by evaluating the periodogram at more bins. One way to evaluate the periodogram at more points is to get more data! That would certainly fix both problems. Unfortunately, we often are struck with the data we have. There is another way.<br />
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgrjsd6LO9cxTJ7s5GZhn6eprHf9KP26Bpf-AJ67Mr_J2TIVprW6duI4s-goOm4Z41wS7eCJL3xlikNr_AgMHBj1PEb4OTDkOJl6NAiXqDmSI6hRIi7uRA4CLbnETFSgTosLaYoxqkW0U4k/s1600/xPer.png" imageanchor="1"><img border="0" data-original-height="600" data-original-width="800" height="300" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgrjsd6LO9cxTJ7s5GZhn6eprHf9KP26Bpf-AJ67Mr_J2TIVprW6duI4s-goOm4Z41wS7eCJL3xlikNr_AgMHBj1PEb4OTDkOJl6NAiXqDmSI6hRIi7uRA4CLbnETFSgTosLaYoxqkW0U4k/s400/xPer.png" width="400" /></a>
Since f is a continuous variable we can evaluate the periodogram at as many points as we want. The way we do this is to pretend we have more data by sticking zeros at the end of our data, we zero pad it. The R code below adds 968 zeros to the end of x zero-padding it to a total of 1000 "data" points. As you can see in the figure below, we have fixed both problems. In the 32 point periodogram missed the peak, because the periodogram was not evaluated at enough points. The zero-padding does not guarantee we hit the peak, but it will get closer. Also, we fixed the jaggedness by evaluating the periodogram at many more bins.<br />
<div class="chunk" id="unnamed-chunk-3">
<div class="rcode">
<div class="source">
<pre class="knitr r"><span class="hl std">xzp</span> <span class="hl kwb"><-</span> <span class="hl kwd">c</span><span class="hl std">(x,</span><span class="hl kwd">rep</span><span class="hl std">(</span><span class="hl num">0</span><span class="hl std">,</span><span class="hl num">1000</span><span class="hl opt">-</span><span class="hl std">N))</span>
<span class="hl std">Nzp</span> <span class="hl kwb"><-</span> <span class="hl kwd">length</span><span class="hl std">(xzp)</span>
</pre>
</div>
</div>
</div>
<div class="chunk" id="unnamed-chunk-4">
<div class="rcode">
<div class="source">
<pre class="knitr r"><span class="hl std">xPerZp</span> <span class="hl kwb"><-</span> <span class="hl std">(</span><span class="hl num">1</span><span class="hl opt">/</span><span class="hl std">N)</span><span class="hl opt">*</span><span class="hl kwd">abs</span><span class="hl std">(</span><span class="hl kwd">fft</span><span class="hl std">(xzp)</span><span class="hl opt">^</span><span class="hl num">2</span><span class="hl std">)</span>
<span class="hl std">fzp</span> <span class="hl kwb"><-</span> <span class="hl kwd">seq</span><span class="hl std">(</span><span class="hl num">0</span><span class="hl std">,</span><span class="hl num">1.0</span><span class="hl opt">-</span><span class="hl num">1</span><span class="hl opt">/</span><span class="hl std">Nzp,</span><span class="hl kwc">by</span><span class="hl std">=</span><span class="hl num">1</span><span class="hl opt">/</span><span class="hl std">Nzp)</span>
</pre>
</div>
</div>
</div>
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEiHau748a6oiXfs8uu5QHz8WWGcUmef70AB-udc3CZfNTE3gOJ_2XbeW-en_a80x4zkRhvch8U1-jTXSdTHXdI2MWELu625YxBd_NDOGDABxJOq1iWtzqYKHY94KyTuQguzncPXK9kXCreZ/s1600/xPerZp.png" imageanchor="1"><img border="0" data-original-height="600" data-original-width="800" height="300" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEiHau748a6oiXfs8uu5QHz8WWGcUmef70AB-udc3CZfNTE3gOJ_2XbeW-en_a80x4zkRhvch8U1-jTXSdTHXdI2MWELu625YxBd_NDOGDABxJOq1iWtzqYKHY94KyTuQguzncPXK9kXCreZ/s400/xPerZp.png" width="400" /></a>
</body>
</html>
Chris Carbonehttp://www.blogger.com/profile/10458001230470047172noreply@blogger.com1tag:blogger.com,1999:blog-7566945490770798840.post-4483834904121133342019-06-03T23:42:00.001-07:002019-06-03T23:44:48.008-07:00JavaScript is Fast?<html>
<head>
<title>JavaScript is Fast?</title>
</head>
<body>
<p>OK, I admit it. I'm always looking for the one true programming language. Mostly I use programming languages for signal processing, data science, that kind of thing. So, for me I want a language to be good for prototyping, flexible, fast, have lots of math/statistics libraries. The first language I use for signal processing was C, not much to say about that. I then became a heavy MATLAB user, great for a lot of things. I used MathCad. It was great for making readable code! I even used Mathematica. I never got the hang of it's syntax, it always felt award. Lately, I've been doing most of my stuff in R.</p>
<p>So, I find myself at the Julia site, at the time is was the next new mathematical programming language. They have a benchmark <a href="http://julialang.org/benchmarks/">page</a> showing how fast it is and it is fast! On that page they show JavaScript as faster than a bunch of standard math mathematically oriented languages. I would not have guessed that. Maybe, JavaScript should be the shiny "new" programming language for math! </p>
</body>
</html>
Chris Carbonehttp://www.blogger.com/profile/10458001230470047172noreply@blogger.com1tag:blogger.com,1999:blog-7566945490770798840.post-62509581172714769252019-06-02T17:15:00.001-07:002019-06-02T21:21:18.467-07:00What is a Statistic?<html>
<head>
<style type="text/css">
.inline {
background-color: #f7f7f7;
border:solid 1px #B0B0B0;
}
.error {
font-weight: bold;
color: #FF0000;
}
.warning {
font-weight: bold;
}
.message {
font-style: italic;
}
.source, .output, .warning, .error, .message {
padding: 0 1em;
border:solid 1px #F7F7F7;
}
.source {
background-color: #f5f5f5;
}
.left {
text-align: left;
}
.right {
text-align: right;
}
.center {
text-align: center;
}
.hl.num {
color: #AF0F91;
}
.hl.str {
color: #317ECC;
}
.hl.com {
color: #AD95AF;
font-style: italic;
}
.hl.opt {
color: #000000;
}
.hl.std {
color: #585858;
}
.hl.kwa {
color: #295F94;
font-weight: bold;
}
.hl.kwb {
color: #B05A65;
}
.hl.kwc {
color: #55aa55;
}
.hl.kwd {
color: #BC5A65;
font-weight: bold;
}
</style>
<title>What is a Statistic?</title>
</head>
<body>
<p>A statistic is a mathematical operation on a data set, performed to get information from the data.</p>
<p>Below is the R code that generates 20 random samples. The samples are uniformly distributed between 0 and 1. Uniformly means all the data samples are equally likely, like when you flip a coin heads and tails are equally likely.</p>
<div class="chunk" id="unnamed-chunk-2"><div class="rcode"><div class="source"><pre class="knitr r"><span class="hl std">x</span> <span class="hl kwb"><-</span> <span class="hl kwd">round</span><span class="hl std">(</span><span class="hl kwd">runif</span><span class="hl std">(</span><span class="hl num">20</span><span class="hl std">),</span><span class="hl num">2</span><span class="hl std">)</span>
<span class="hl std">x</span>
</pre></div>
<div class="output"><pre class="knitr r">## [1] 0.90 0.27 0.37 0.57 0.91 0.20 0.90 0.94 0.66 0.63 0.06 0.21 0.18
## [14] 0.69 0.38 0.77 0.50 0.72 0.99 0.38
</pre></div>
</div></div>
<p>One thing we might want to know about this data set it the expected value (EV). The EV is the typical value of the data. Since we know the data is uniformly distributed between 0 and 1, the EV is just the middle of the range, 0.5. In practice, we have the data, but don't know the distribution the data came from. So, we can't calculate the EV, but we could use the average as an estimate of the EV. The R code below is the average.</p>
<div class="chunk" id="unnamed-chunk-3"><div class="rcode"><div class="source"><pre class="knitr r"><span class="hl kwd">sum</span><span class="hl std">(x)</span><span class="hl opt">/</span><span class="hl num">20</span>
</pre></div>
<div class="output"><pre class="knitr r">## [1] 0.5615
</pre></div>
</div></div>
<p>The average of the data is about 0.56. The average is not the only way to estimate the EV, there are many ways! We could just use the first element of the data set 0.90 as the estimate. Another way is to find the maximum of the data and divide it by 2. The code below finds the max of the data and divides it by 2. Notice the max/2 produces an estimate that is closer to the true value of 0.5 than the average.</p>
<div class="chunk" id="unnamed-chunk-4"><div class="rcode"><div class="source"><pre class="knitr r"><span class="hl kwd">max</span><span class="hl std">(x)</span><span class="hl opt">/</span><span class="hl num">2</span>
</pre></div>
<div class="output"><pre class="knitr r">## [1] 0.495
</pre></div>
</div></div>
<p>There are many way to estimate the EV, but some are better than others. Estimators of random data are also random. So, the only way to compare estimators is statistically. The case above where the max/2 produced a better estimate than the average could have been luck. The code below generates 20 data samples from the same uniform distribution above 20 times. We then find the estimate the EV in three ways: the average, the first sample, and the max(x)/2.</p>
<div class="chunk" id="unnamed-chunk-5"><div class="rcode"><div class="source"><pre class="knitr r"><span class="hl std">trials</span> <span class="hl kwb"><-</span> <span class="hl num">100</span>
<span class="hl std">avg</span> <span class="hl kwb"><-</span> <span class="hl kwd">matrix</span><span class="hl std">(</span><span class="hl num">0</span><span class="hl std">,</span><span class="hl num">1</span><span class="hl std">,trials)</span>
<span class="hl std">first</span> <span class="hl kwb"><-</span> <span class="hl kwd">matrix</span><span class="hl std">(</span><span class="hl num">0</span><span class="hl std">,</span><span class="hl num">1</span><span class="hl std">,trials)</span>
<span class="hl std">halfMax</span> <span class="hl kwb"><-</span> <span class="hl kwd">matrix</span><span class="hl std">(</span><span class="hl num">0</span><span class="hl std">,</span><span class="hl num">1</span><span class="hl std">,trials)</span>
<span class="hl kwa">for</span> <span class="hl std">(indx</span> <span class="hl kwa">in</span> <span class="hl num">1</span><span class="hl opt">:</span><span class="hl std">trials) {</span>
<span class="hl std">x</span> <span class="hl kwb"><-</span> <span class="hl kwd">runif</span><span class="hl std">(</span><span class="hl num">20</span><span class="hl std">)</span>
<span class="hl std">avg[indx]</span> <span class="hl kwb"><-</span> <span class="hl kwd">sum</span><span class="hl std">(x)</span><span class="hl opt">/</span><span class="hl num">20</span>
<span class="hl std">first[indx]</span> <span class="hl kwb"><-</span> <span class="hl std">x[</span><span class="hl num">1</span><span class="hl std">]</span>
<span class="hl std">halfMax[indx]</span> <span class="hl kwb"><-</span> <span class="hl kwd">max</span><span class="hl std">(x)</span><span class="hl opt">/</span><span class="hl num">2</span>
<span class="hl std">}</span>
</pre></div>
</div></div>
<p>The figure below contains 100 averages plotted with circles and true EV plotted in a line at 0.5. Most of the averages are between 0.4 and 0.6 and many are closer. The average seems to be a good estimate of the EV. Which shouldn't be surprising since the average is the most common statistic.</p>
<!-- *************************** -->
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEja9Z-JIw4HbdEqZeirKLkWx_MSEg5wGrApebBTsSVX8AY7Agi72mDYqGUVl_tdk4KCv_iLlPslqHCjt6dkZKQatoRSEcKs2R_1gtiLTx1QYuO9DsL24JDtWKYuq2iUfEqxeGnjAfRZ8tg7/s1600/Avarage.png" imageanchor="1" ><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEja9Z-JIw4HbdEqZeirKLkWx_MSEg5wGrApebBTsSVX8AY7Agi72mDYqGUVl_tdk4KCv_iLlPslqHCjt6dkZKQatoRSEcKs2R_1gtiLTx1QYuO9DsL24JDtWKYuq2iUfEqxeGnjAfRZ8tg7/s400/Avarage.png" width="400" height="300" data-original-width="800" data-original-height="600" /></a>
<p>This figure is the estimate based only on the first element of the data set. The estimates are between just over 0 and almost 1. A few estimates are close to 0.5, but most are not. This is not a good estimator for EV. This is probably not surprising since this estimate is based on only one sample and the average is based on 20.</p>
<!-- *************************** -->
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjiphxATjNiHK45SL-rrRTEfD074bHj1N2QxpbtYuvqQZQipeYyf7jhORL8wlLm2JrILxUZleL6ccgCj9Eft1QOh1FxT-MRwWwqy7vgvgSnoZR7dFrwJIvftrsT_OStj4Mepo0YXqeHRhHR/s1600/First.png" imageanchor="1" ><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjiphxATjNiHK45SL-rrRTEfD074bHj1N2QxpbtYuvqQZQipeYyf7jhORL8wlLm2JrILxUZleL6ccgCj9Eft1QOh1FxT-MRwWwqy7vgvgSnoZR7dFrwJIvftrsT_OStj4Mepo0YXqeHRhHR/s400/First.png" width="400" height="300" data-original-width="800" data-original-height="600" /></a>
<p>This figure is the plot of the estimate based on the max/2. Most of the estimates are very close to the true value. It looks even better than the average! This may seem surprising, since it looks like it is based on only one sample, the maximum sample, but it is actually makes use of all the samples. This is because you need to look at all the samples to know which one is the maximum.</p>
<!-- *************************** -->
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjrMRMhYaYbcH5jz8kHtOSo3X411-IN1VXfAZyngdzK6ugw_ZQ2xd-NtZjRtEUdsi3Crrh0GLW6Zo2LTz_69vC2KUTrGY9yxqQiakYUZPcx4GLMkB-U8L_QmT1gU9wh3NR0hy0-CJCsIWjV/s1600/Max.png" imageanchor="1" ><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjrMRMhYaYbcH5jz8kHtOSo3X411-IN1VXfAZyngdzK6ugw_ZQ2xd-NtZjRtEUdsi3Crrh0GLW6Zo2LTz_69vC2KUTrGY9yxqQiakYUZPcx4GLMkB-U8L_QmT1gU9wh3NR0hy0-CJCsIWjV/s400/Max.png" width="400" height="300" data-original-width="800" data-original-height="600" /></a>
<p>The average is a good choice (and often the best choice) for estimating the EV, but it is not always the best estimator for EV. In the case of this particular distribution, the max/2 is better estimator of EV than the average.</p>
<p>So, we now know what a statistic is and we worked with a few. As a field of study, statistics uses visualization, probability and other math tools to find the best ways to get information from data.</p>
</body>
</html>
Chris Carbonehttp://www.blogger.com/profile/10458001230470047172noreply@blogger.com0tag:blogger.com,1999:blog-7566945490770798840.post-77045140147915617972019-05-27T19:34:00.001-07:002019-06-27T18:44:47.487-07:00Book Review: Data Analysis with Open Source Tools<br />
<span style="font-size: large;"></span><br />
<br />
<div style="margin-bottom: .0001pt; margin: 0in;">
<span style="color: black; font-family: "Times","serif"; font-size: 13.5pt;">This book gives a very good
overview of the kinds of things data scientists do and the tools they
use. All the concepts are illustrated with hands on programing examples
you can follow along with. <o:p></o:p></span></div>
<div style="-webkit-text-stroke-width: 0px; font-variant-caps: normal; font-variant-ligatures: normal; margin-bottom: .0001pt; margin: 0in; orphans: 2; text-align: start; text-decoration-color: initial; text-decoration-style: initial; widows: 2; word-spacing: 0px;">
<br /></div>
<div style="-webkit-text-stroke-width: 0px; font-variant-caps: normal; font-variant-ligatures: normal; margin-bottom: .0001pt; margin: 0in; orphans: 2; text-align: start; text-decoration-color: initial; text-decoration-style: initial; widows: 2; word-spacing: 0px;">
<span style="color: black; font-family: "Times","serif"; font-size: 13.5pt;">I
think this book would be very good for two groups of people. First, someone working in a technical field who wants to know what data science is all about. This would give them a good idea of the skills they
will need and the type of work they may do as a data scientist.<o:p></o:p></span></div>
<div style="-webkit-text-stroke-width: 0px; font-variant-caps: normal; font-variant-ligatures: normal; margin-bottom: .0001pt; margin: 0in; orphans: 2; text-align: start; text-decoration-color: initial; text-decoration-style: initial; widows: 2; word-spacing: 0px;">
<br /></div>
<div style="-webkit-text-stroke-width: 0px; font-variant-caps: normal; font-variant-ligatures: normal; margin-bottom: .0001pt; margin: 0in; orphans: 2; text-align: start; text-decoration-color: initial; text-decoration-style: initial; widows: 2; word-spacing: 0px;">
<span style="color: black; font-family: "Times","serif"; font-size: 13.5pt;">The
second group that would benefit from this book are people engaged in some field
of data science and who are looking for a broad overview of the field. Since
many people came to data science from other fields and learned hands on they
may have trouble seeing the forest for the trees. <o:p></o:p></span></div>
<div style="-webkit-text-stroke-width: 0px; font-variant-caps: normal; font-variant-ligatures: normal; margin-bottom: .0001pt; margin: 0in; orphans: 2; text-align: start; text-decoration-color: initial; text-decoration-style: initial; widows: 2; word-spacing: 0px;">
<br /></div>
<!--[if gte mso 9]><xml>
<o:OfficeDocumentSettings>
<o:AllowPNG/>
</o:OfficeDocumentSettings>
</xml><![endif]-->
<!--[if gte mso 9]><xml>
<w:WordDocument>
<w:View>Normal</w:View>
<w:Zoom>0</w:Zoom>
<w:TrackMoves/>
<w:TrackFormatting/>
<w:PunctuationKerning/>
<w:ValidateAgainstSchemas/>
<w:SaveIfXMLInvalid>false</w:SaveIfXMLInvalid>
<w:IgnoreMixedContent>false</w:IgnoreMixedContent>
<w:AlwaysShowPlaceholderText>false</w:AlwaysShowPlaceholderText>
<w:DoNotPromoteQF/>
<w:LidThemeOther>EN-US</w:LidThemeOther>
<w:LidThemeAsian>JA</w:LidThemeAsian>
<w:LidThemeComplexScript>X-NONE</w:LidThemeComplexScript>
<w:Compatibility>
<w:BreakWrappedTables/>
<w:SnapToGridInCell/>
<w:WrapTextWithPunct/>
<w:UseAsianBreakRules/>
<w:DontGrowAutofit/>
<w:SplitPgBreakAndParaMark/>
<w:EnableOpenTypeKerning/>
<w:DontFlipMirrorIndents/>
<w:OverrideTableStyleHps/>
<w:UseFELayout/>
</w:Compatibility>
<m:mathPr>
<m:mathFont m:val="Cambria Math"/>
<m:brkBin m:val="before"/>
<m:brkBinSub m:val="--"/>
<m:smallFrac m:val="off"/>
<m:dispDef/>
<m:lMargin m:val="0"/>
<m:rMargin m:val="0"/>
<m:defJc m:val="centerGroup"/>
<m:wrapIndent m:val="1440"/>
<m:intLim m:val="subSup"/>
<m:naryLim m:val="undOvr"/>
</m:mathPr></w:WordDocument>
</xml><![endif]--><!--[if gte mso 9]><xml>
<w:LatentStyles DefLockedState="false" DefUnhideWhenUsed="true"
DefSemiHidden="true" DefQFormat="false" DefPriority="99"
LatentStyleCount="276">
<w:LsdException Locked="false" Priority="0" SemiHidden="false"
UnhideWhenUsed="false" QFormat="true" Name="Normal"/>
<w:LsdException Locked="false" Priority="9" SemiHidden="false"
UnhideWhenUsed="false" QFormat="true" Name="heading 1"/>
<w:LsdException Locked="false" Priority="9" QFormat="true" Name="heading 2"/>
<w:LsdException Locked="false" Priority="9" QFormat="true" Name="heading 3"/>
<w:LsdException Locked="false" Priority="9" QFormat="true" Name="heading 4"/>
<w:LsdException Locked="false" Priority="9" QFormat="true" Name="heading 5"/>
<w:LsdException Locked="false" Priority="9" QFormat="true" Name="heading 6"/>
<w:LsdException Locked="false" Priority="9" QFormat="true" Name="heading 7"/>
<w:LsdException Locked="false" Priority="9" QFormat="true" Name="heading 8"/>
<w:LsdException Locked="false" Priority="9" QFormat="true" Name="heading 9"/>
<w:LsdException Locked="false" Priority="39" Name="toc 1"/>
<w:LsdException Locked="false" Priority="39" Name="toc 2"/>
<w:LsdException Locked="false" Priority="39" Name="toc 3"/>
<w:LsdException Locked="false" Priority="39" Name="toc 4"/>
<w:LsdException Locked="false" Priority="39" Name="toc 5"/>
<w:LsdException Locked="false" Priority="39" Name="toc 6"/>
<w:LsdException Locked="false" Priority="39" Name="toc 7"/>
<w:LsdException Locked="false" Priority="39" Name="toc 8"/>
<w:LsdException Locked="false" Priority="39" Name="toc 9"/>
<w:LsdException Locked="false" Priority="35" QFormat="true" Name="caption"/>
<w:LsdException Locked="false" Priority="10" SemiHidden="false"
UnhideWhenUsed="false" QFormat="true" Name="Title"/>
<w:LsdException Locked="false" Priority="1" Name="Default Paragraph Font"/>
<w:LsdException Locked="false" Priority="11" SemiHidden="false"
UnhideWhenUsed="false" QFormat="true" Name="Subtitle"/>
<w:LsdException Locked="false" Priority="22" SemiHidden="false"
UnhideWhenUsed="false" QFormat="true" Name="Strong"/>
<w:LsdException Locked="false" Priority="20" SemiHidden="false"
UnhideWhenUsed="false" QFormat="true" Name="Emphasis"/>
<w:LsdException Locked="false" Priority="59" SemiHidden="false"
UnhideWhenUsed="false" Name="Table Grid"/>
<w:LsdException Locked="false" UnhideWhenUsed="false" Name="Placeholder Text"/>
<w:LsdException Locked="false" Priority="1" SemiHidden="false"
UnhideWhenUsed="false" QFormat="true" Name="No Spacing"/>
<w:LsdException Locked="false" Priority="60" SemiHidden="false"
UnhideWhenUsed="false" Name="Light Shading"/>
<w:LsdException Locked="false" Priority="61" SemiHidden="false"
UnhideWhenUsed="false" Name="Light List"/>
<w:LsdException Locked="false" Priority="62" SemiHidden="false"
UnhideWhenUsed="false" Name="Light Grid"/>
<w:LsdException Locked="false" Priority="63" SemiHidden="false"
UnhideWhenUsed="false" Name="Medium Shading 1"/>
<w:LsdException Locked="false" Priority="64" SemiHidden="false"
UnhideWhenUsed="false" Name="Medium Shading 2"/>
<w:LsdException Locked="false" Priority="65" SemiHidden="false"
UnhideWhenUsed="false" Name="Medium List 1"/>
<w:LsdException Locked="false" Priority="66" SemiHidden="false"
UnhideWhenUsed="false" Name="Medium List 2"/>
<w:LsdException Locked="false" Priority="67" SemiHidden="false"
UnhideWhenUsed="false" Name="Medium Grid 1"/>
<w:LsdException Locked="false" Priority="68" SemiHidden="false"
UnhideWhenUsed="false" Name="Medium Grid 2"/>
<w:LsdException Locked="false" Priority="69" SemiHidden="false"
UnhideWhenUsed="false" Name="Medium Grid 3"/>
<w:LsdException Locked="false" Priority="70" SemiHidden="false"
UnhideWhenUsed="false" Name="Dark List"/>
<w:LsdException Locked="false" Priority="71" SemiHidden="false"
UnhideWhenUsed="false" Name="Colorful Shading"/>
<w:LsdException Locked="false" Priority="72" SemiHidden="false"
UnhideWhenUsed="false" Name="Colorful List"/>
<w:LsdException Locked="false" Priority="73" SemiHidden="false"
UnhideWhenUsed="false" Name="Colorful Grid"/>
<w:LsdException Locked="false" Priority="60" SemiHidden="false"
UnhideWhenUsed="false" Name="Light Shading Accent 1"/>
<w:LsdException Locked="false" Priority="61" SemiHidden="false"
UnhideWhenUsed="false" Name="Light List Accent 1"/>
<w:LsdException Locked="false" Priority="62" SemiHidden="false"
UnhideWhenUsed="false" Name="Light Grid Accent 1"/>
<w:LsdException Locked="false" Priority="63" SemiHidden="false"
UnhideWhenUsed="false" Name="Medium Shading 1 Accent 1"/>
<w:LsdException Locked="false" Priority="64" SemiHidden="false"
UnhideWhenUsed="false" Name="Medium Shading 2 Accent 1"/>
<w:LsdException Locked="false" Priority="65" SemiHidden="false"
UnhideWhenUsed="false" Name="Medium List 1 Accent 1"/>
<w:LsdException Locked="false" UnhideWhenUsed="false" Name="Revision"/>
<w:LsdException Locked="false" Priority="34" SemiHidden="false"
UnhideWhenUsed="false" QFormat="true" Name="List Paragraph"/>
<w:LsdException Locked="false" Priority="29" SemiHidden="false"
UnhideWhenUsed="false" QFormat="true" Name="Quote"/>
<w:LsdException Locked="false" Priority="30" SemiHidden="false"
UnhideWhenUsed="false" QFormat="true" Name="Intense Quote"/>
<w:LsdException Locked="false" Priority="66" SemiHidden="false"
UnhideWhenUsed="false" Name="Medium List 2 Accent 1"/>
<w:LsdException Locked="false" Priority="67" SemiHidden="false"
UnhideWhenUsed="false" Name="Medium Grid 1 Accent 1"/>
<w:LsdException Locked="false" Priority="68" SemiHidden="false"
UnhideWhenUsed="false" Name="Medium Grid 2 Accent 1"/>
<w:LsdException Locked="false" Priority="69" SemiHidden="false"
UnhideWhenUsed="false" Name="Medium Grid 3 Accent 1"/>
<w:LsdException Locked="false" Priority="70" SemiHidden="false"
UnhideWhenUsed="false" Name="Dark List Accent 1"/>
<w:LsdException Locked="false" Priority="71" SemiHidden="false"
UnhideWhenUsed="false" Name="Colorful Shading Accent 1"/>
<w:LsdException Locked="false" Priority="72" SemiHidden="false"
UnhideWhenUsed="false" Name="Colorful List Accent 1"/>
<w:LsdException Locked="false" Priority="73" SemiHidden="false"
UnhideWhenUsed="false" Name="Colorful Grid Accent 1"/>
<w:LsdException Locked="false" Priority="60" SemiHidden="false"
UnhideWhenUsed="false" Name="Light Shading Accent 2"/>
<w:LsdException Locked="false" Priority="61" SemiHidden="false"
UnhideWhenUsed="false" Name="Light List Accent 2"/>
<w:LsdException Locked="false" Priority="62" SemiHidden="false"
UnhideWhenUsed="false" Name="Light Grid Accent 2"/>
<w:LsdException Locked="false" Priority="63" SemiHidden="false"
UnhideWhenUsed="false" Name="Medium Shading 1 Accent 2"/>
<w:LsdException Locked="false" Priority="64" SemiHidden="false"
UnhideWhenUsed="false" Name="Medium Shading 2 Accent 2"/>
<w:LsdException Locked="false" Priority="65" SemiHidden="false"
UnhideWhenUsed="false" Name="Medium List 1 Accent 2"/>
<w:LsdException Locked="false" Priority="66" SemiHidden="false"
UnhideWhenUsed="false" Name="Medium List 2 Accent 2"/>
<w:LsdException Locked="false" Priority="67" SemiHidden="false"
UnhideWhenUsed="false" Name="Medium Grid 1 Accent 2"/>
<w:LsdException Locked="false" Priority="68" SemiHidden="false"
UnhideWhenUsed="false" Name="Medium Grid 2 Accent 2"/>
<w:LsdException Locked="false" Priority="69" SemiHidden="false"
UnhideWhenUsed="false" Name="Medium Grid 3 Accent 2"/>
<w:LsdException Locked="false" Priority="70" SemiHidden="false"
UnhideWhenUsed="false" Name="Dark List Accent 2"/>
<w:LsdException Locked="false" Priority="71" SemiHidden="false"
UnhideWhenUsed="false" Name="Colorful Shading Accent 2"/>
<w:LsdException Locked="false" Priority="72" SemiHidden="false"
UnhideWhenUsed="false" Name="Colorful List Accent 2"/>
<w:LsdException Locked="false" Priority="73" SemiHidden="false"
UnhideWhenUsed="false" Name="Colorful Grid Accent 2"/>
<w:LsdException Locked="false" Priority="60" SemiHidden="false"
UnhideWhenUsed="false" Name="Light Shading Accent 3"/>
<w:LsdException Locked="false" Priority="61" SemiHidden="false"
UnhideWhenUsed="false" Name="Light List Accent 3"/>
<w:LsdException Locked="false" Priority="62" SemiHidden="false"
UnhideWhenUsed="false" Name="Light Grid Accent 3"/>
<w:LsdException Locked="false" Priority="63" SemiHidden="false"
UnhideWhenUsed="false" Name="Medium Shading 1 Accent 3"/>
<w:LsdException Locked="false" Priority="64" SemiHidden="false"
UnhideWhenUsed="false" Name="Medium Shading 2 Accent 3"/>
<w:LsdException Locked="false" Priority="65" SemiHidden="false"
UnhideWhenUsed="false" Name="Medium List 1 Accent 3"/>
<w:LsdException Locked="false" Priority="66" SemiHidden="false"
UnhideWhenUsed="false" Name="Medium List 2 Accent 3"/>
<w:LsdException Locked="false" Priority="67" SemiHidden="false"
UnhideWhenUsed="false" Name="Medium Grid 1 Accent 3"/>
<w:LsdException Locked="false" Priority="68" SemiHidden="false"
UnhideWhenUsed="false" Name="Medium Grid 2 Accent 3"/>
<w:LsdException Locked="false" Priority="69" SemiHidden="false"
UnhideWhenUsed="false" Name="Medium Grid 3 Accent 3"/>
<w:LsdException Locked="false" Priority="70" SemiHidden="false"
UnhideWhenUsed="false" Name="Dark List Accent 3"/>
<w:LsdException Locked="false" Priority="71" SemiHidden="false"
UnhideWhenUsed="false" Name="Colorful Shading Accent 3"/>
<w:LsdException Locked="false" Priority="72" SemiHidden="false"
UnhideWhenUsed="false" Name="Colorful List Accent 3"/>
<w:LsdException Locked="false" Priority="73" SemiHidden="false"
UnhideWhenUsed="false" Name="Colorful Grid Accent 3"/>
<w:LsdException Locked="false" Priority="60" SemiHidden="false"
UnhideWhenUsed="false" Name="Light Shading Accent 4"/>
<w:LsdException Locked="false" Priority="61" SemiHidden="false"
UnhideWhenUsed="false" Name="Light List Accent 4"/>
<w:LsdException Locked="false" Priority="62" SemiHidden="false"
UnhideWhenUsed="false" Name="Light Grid Accent 4"/>
<w:LsdException Locked="false" Priority="63" SemiHidden="false"
UnhideWhenUsed="false" Name="Medium Shading 1 Accent 4"/>
<w:LsdException Locked="false" Priority="64" SemiHidden="false"
UnhideWhenUsed="false" Name="Medium Shading 2 Accent 4"/>
<w:LsdException Locked="false" Priority="65" SemiHidden="false"
UnhideWhenUsed="false" Name="Medium List 1 Accent 4"/>
<w:LsdException Locked="false" Priority="66" SemiHidden="false"
UnhideWhenUsed="false" Name="Medium List 2 Accent 4"/>
<w:LsdException Locked="false" Priority="67" SemiHidden="false"
UnhideWhenUsed="false" Name="Medium Grid 1 Accent 4"/>
<w:LsdException Locked="false" Priority="68" SemiHidden="false"
UnhideWhenUsed="false" Name="Medium Grid 2 Accent 4"/>
<w:LsdException Locked="false" Priority="69" SemiHidden="false"
UnhideWhenUsed="false" Name="Medium Grid 3 Accent 4"/>
<w:LsdException Locked="false" Priority="70" SemiHidden="false"
UnhideWhenUsed="false" Name="Dark List Accent 4"/>
<w:LsdException Locked="false" Priority="71" SemiHidden="false"
UnhideWhenUsed="false" Name="Colorful Shading Accent 4"/>
<w:LsdException Locked="false" Priority="72" SemiHidden="false"
UnhideWhenUsed="false" Name="Colorful List Accent 4"/>
<w:LsdException Locked="false" Priority="73" SemiHidden="false"
UnhideWhenUsed="false" Name="Colorful Grid Accent 4"/>
<w:LsdException Locked="false" Priority="60" SemiHidden="false"
UnhideWhenUsed="false" Name="Light Shading Accent 5"/>
<w:LsdException Locked="false" Priority="61" SemiHidden="false"
UnhideWhenUsed="false" Name="Light List Accent 5"/>
<w:LsdException Locked="false" Priority="62" SemiHidden="false"
UnhideWhenUsed="false" Name="Light Grid Accent 5"/>
<w:LsdException Locked="false" Priority="63" SemiHidden="false"
UnhideWhenUsed="false" Name="Medium Shading 1 Accent 5"/>
<w:LsdException Locked="false" Priority="64" SemiHidden="false"
UnhideWhenUsed="false" Name="Medium Shading 2 Accent 5"/>
<w:LsdException Locked="false" Priority="65" SemiHidden="false"
UnhideWhenUsed="false" Name="Medium List 1 Accent 5"/>
<w:LsdException Locked="false" Priority="66" SemiHidden="false"
UnhideWhenUsed="false" Name="Medium List 2 Accent 5"/>
<w:LsdException Locked="false" Priority="67" SemiHidden="false"
UnhideWhenUsed="false" Name="Medium Grid 1 Accent 5"/>
<w:LsdException Locked="false" Priority="68" SemiHidden="false"
UnhideWhenUsed="false" Name="Medium Grid 2 Accent 5"/>
<w:LsdException Locked="false" Priority="69" SemiHidden="false"
UnhideWhenUsed="false" Name="Medium Grid 3 Accent 5"/>
<w:LsdException Locked="false" Priority="70" SemiHidden="false"
UnhideWhenUsed="false" Name="Dark List Accent 5"/>
<w:LsdException Locked="false" Priority="71" SemiHidden="false"
UnhideWhenUsed="false" Name="Colorful Shading Accent 5"/>
<w:LsdException Locked="false" Priority="72" SemiHidden="false"
UnhideWhenUsed="false" Name="Colorful List Accent 5"/>
<w:LsdException Locked="false" Priority="73" SemiHidden="false"
UnhideWhenUsed="false" Name="Colorful Grid Accent 5"/>
<w:LsdException Locked="false" Priority="60" SemiHidden="false"
UnhideWhenUsed="false" Name="Light Shading Accent 6"/>
<w:LsdException Locked="false" Priority="61" SemiHidden="false"
UnhideWhenUsed="false" Name="Light List Accent 6"/>
<w:LsdException Locked="false" Priority="62" SemiHidden="false"
UnhideWhenUsed="false" Name="Light Grid Accent 6"/>
<w:LsdException Locked="false" Priority="63" SemiHidden="false"
UnhideWhenUsed="false" Name="Medium Shading 1 Accent 6"/>
<w:LsdException Locked="false" Priority="64" SemiHidden="false"
UnhideWhenUsed="false" Name="Medium Shading 2 Accent 6"/>
<w:LsdException Locked="false" Priority="65" SemiHidden="false"
UnhideWhenUsed="false" Name="Medium List 1 Accent 6"/>
<w:LsdException Locked="false" Priority="66" SemiHidden="false"
UnhideWhenUsed="false" Name="Medium List 2 Accent 6"/>
<w:LsdException Locked="false" Priority="67" SemiHidden="false"
UnhideWhenUsed="false" Name="Medium Grid 1 Accent 6"/>
<w:LsdException Locked="false" Priority="68" SemiHidden="false"
UnhideWhenUsed="false" Name="Medium Grid 2 Accent 6"/>
<w:LsdException Locked="false" Priority="69" SemiHidden="false"
UnhideWhenUsed="false" Name="Medium Grid 3 Accent 6"/>
<w:LsdException Locked="false" Priority="70" SemiHidden="false"
UnhideWhenUsed="false" Name="Dark List Accent 6"/>
<w:LsdException Locked="false" Priority="71" SemiHidden="false"
UnhideWhenUsed="false" Name="Colorful Shading Accent 6"/>
<w:LsdException Locked="false" Priority="72" SemiHidden="false"
UnhideWhenUsed="false" Name="Colorful List Accent 6"/>
<w:LsdException Locked="false" Priority="73" SemiHidden="false"
UnhideWhenUsed="false" Name="Colorful Grid Accent 6"/>
<w:LsdException Locked="false" Priority="19" SemiHidden="false"
UnhideWhenUsed="false" QFormat="true" Name="Subtle Emphasis"/>
<w:LsdException Locked="false" Priority="21" SemiHidden="false"
UnhideWhenUsed="false" QFormat="true" Name="Intense Emphasis"/>
<w:LsdException Locked="false" Priority="31" SemiHidden="false"
UnhideWhenUsed="false" QFormat="true" Name="Subtle Reference"/>
<w:LsdException Locked="false" Priority="32" SemiHidden="false"
UnhideWhenUsed="false" QFormat="true" Name="Intense Reference"/>
<w:LsdException Locked="false" Priority="33" SemiHidden="false"
UnhideWhenUsed="false" QFormat="true" Name="Book Title"/>
<w:LsdException Locked="false" Priority="37" Name="Bibliography"/>
<w:LsdException Locked="false" Priority="39" QFormat="true" Name="TOC Heading"/>
</w:LatentStyles>
</xml><![endif]-->
<style>
<!--
/* Font Definitions */
@font-face
{font-family:Times;
panose-1:2 0 5 0 0 0 0 0 0 0;
mso-font-alt:"Times Roman";
mso-font-charset:77;
mso-generic-font-family:roman;
mso-font-format:other;
mso-font-pitch:variable;
mso-font-signature:3 0 0 0 1 0;}
@font-face
{font-family:"MS 明朝";
mso-font-charset:78;
mso-generic-font-family:auto;
mso-font-pitch:variable;
mso-font-signature:1 134676480 16 0 131072 0;}
@font-face
{font-family:"MS 明朝";
mso-font-charset:78;
mso-generic-font-family:auto;
mso-font-pitch:variable;
mso-font-signature:1 134676480 16 0 131072 0;}
@font-face
{font-family:Cambria;
panose-1:2 4 5 3 5 4 6 3 2 4;
mso-font-charset:0;
mso-generic-font-family:auto;
mso-font-pitch:variable;
mso-font-signature:-536870145 1073743103 0 0 415 0;}
/* Style Definitions */
p.MsoNormal, li.MsoNormal, div.MsoNormal
{mso-style-unhide:no;
mso-style-qformat:yes;
mso-style-parent:"";
margin:0in;
margin-bottom:.0001pt;
mso-pagination:widow-orphan;
font-size:12.0pt;
font-family:Cambria;
mso-ascii-font-family:Cambria;
mso-ascii-theme-font:minor-latin;
mso-fareast-font-family:"MS 明朝";
mso-fareast-theme-font:minor-fareast;
mso-hansi-font-family:Cambria;
mso-hansi-theme-font:minor-latin;
mso-bidi-font-family:"Times New Roman";
mso-bidi-theme-font:minor-bidi;}
p
{mso-style-noshow:yes;
mso-style-priority:99;
mso-margin-top-alt:auto;
margin-right:0in;
mso-margin-bottom-alt:auto;
margin-left:0in;
mso-pagination:widow-orphan;
font-size:10.0pt;
font-family:"Times New Roman";
mso-fareast-font-family:"MS 明朝";
mso-fareast-theme-font:minor-fareast;}
.MsoChpDefault
{mso-style-type:export-only;
mso-default-props:yes;
font-family:Cambria;
mso-ascii-font-family:Cambria;
mso-ascii-theme-font:minor-latin;
mso-fareast-font-family:"MS 明朝";
mso-fareast-theme-font:minor-fareast;
mso-hansi-font-family:Cambria;
mso-hansi-theme-font:minor-latin;
mso-bidi-font-family:"Times New Roman";
mso-bidi-theme-font:minor-bidi;}
@page WordSection1
{size:8.5in 11.0in;
margin:1.0in 1.25in 1.0in 1.25in;
mso-header-margin:.5in;
mso-footer-margin:.5in;
mso-paper-source:0;}
div.WordSection1
{page:WordSection1;}
-->
</style>
<!--[if gte mso 10]>
<style>
/* Style Definitions */
table.MsoNormalTable
{mso-style-name:"Table Normal";
mso-tstyle-rowband-size:0;
mso-tstyle-colband-size:0;
mso-style-noshow:yes;
mso-style-priority:99;
mso-style-parent:"";
mso-padding-alt:0in 5.4pt 0in 5.4pt;
mso-para-margin:0in;
mso-para-margin-bottom:.0001pt;
mso-pagination:widow-orphan;
font-size:12.0pt;
font-family:Cambria;
mso-ascii-font-family:Cambria;
mso-ascii-theme-font:minor-latin;
mso-hansi-font-family:Cambria;
mso-hansi-theme-font:minor-latin;}
</style>
<![endif]-->
<!--StartFragment-->
<!--EndFragment--><br />
<div style="-webkit-text-stroke-width: 0px; font-variant-caps: normal; font-variant-ligatures: normal; margin-bottom: .0001pt; margin: 0in; orphans: 2; text-align: start; text-decoration-color: initial; text-decoration-style: initial; widows: 2; word-spacing: 0px;">
<span style="color: black; font-family: "Times","serif"; font-size: 13.5pt;">This
is not a pop science book. If you looking for an armchairs introduction to the wonderful
world of data science, this is not the book for you. If you want to learn by
doing and know a bit of math and have some familiarity with programing, then
this <i>is</i> the book for you. <o:p></o:p></span></div>
<br />
<a href="https://www.amazon.com/Data-Analysis-Source-Tools-Hands/dp/0596802358">https://www.amazon.com/Data-Analysis-Source-Tools-Hands/dp/0596802358</a>Chris Carbonehttp://www.blogger.com/profile/10458001230470047172noreply@blogger.com0tag:blogger.com,1999:blog-7566945490770798840.post-41579089990242152442019-05-26T19:17:00.001-07:002019-05-26T19:17:54.400-07:00Mirror image data with MATLABA sequence is even if x[midpoint-n] = x[midpoint+n] for some midpoint and for all n. There are many cases where it is handy to have a even sequence. For Instance, the Fourier transform of an even sequence is real and it is handy to have a real function in the transform domain.<br />
<br />
Below is the some MATLAB code that will make a sequence even by appending a mirror image section. The function will mirror image one, two and three dimensional data.<br />
<br />
<style type="text/css">
p.p1 {margin: 0.0px 0.0px 0.0px 0.0px; font: 10.0px Courier}
p.p2 {margin: 0.0px 0.0px 0.0px 0.0px; font: 10.0px Courier; color: #228b22}
p.p3 {margin: 0.0px 0.0px 0.0px 0.0px; font: 10.0px Courier; color: #228b22; min-height: 12.0px}
p.p4 {margin: 0.0px 0.0px 0.0px 0.0px; font: 10.0px Courier; min-height: 12.0px}
p.p5 {margin: 0.0px 0.0px 0.0px 0.0px; font: 10.0px Courier; color: #0000ff}
p.p6 {margin: 0.0px 0.0px 0.0px 0.0px; font: 10.0px Courier; color: #0000ff; min-height: 12.0px}
span.s1 {color: #0000ff}
</style>
<br />
<div class="p1">
<span class="s1">function</span> DataMirrorImage = MirrorImageData(Data)</div>
<div class="p2">
% This function mirror images one, two and three dimensional data in all</div>
<div class="p2">
% in one, two or three dimensions.</div>
<div class="p2">
%</div>
<div class="p2">
% INPUT:<span class="Apple-converted-space"> </span></div>
<div class="p2">
% Data3D................The data to be mirror imaged.</div>
<div class="p2">
%<span class="Apple-converted-space"> </span></div>
<div class="p2">
% OUTPUT:</div>
<div class="p2">
% DataMirrorImage.....The mirror image data.</div>
<div class="p3">
<br /></div>
<div class="p2">
% Get the dimensions of the data.</div>
<div class="p1">
[xLen, yLen, zLen] = size(Data);</div>
<div class="p4">
<br /></div>
<div class="p2">
% Mirror image the data in 1D or 2D.</div>
<div class="p1">
<span class="s1">for</span> zNdx = 1:zLen</div>
<div class="p1">
<span class="Apple-converted-space"> </span>Data2D = squeeze(Data(:,:,zNdx));</div>
<div class="p1">
<span class="Apple-converted-space"> </span>Data2D = [ Data2D fliplr(Data2D(:,2:yLen))];</div>
<div class="p1">
<span class="Apple-converted-space"> </span>Data2DMirrorImage = [ Data2D; flipud(Data2D(2:xLen,:))];</div>
<div class="p1">
<span class="Apple-converted-space"> </span>DataMirrorImage(:,:,zNdx) = Data2DMirrorImage;</div>
<div class="p5">
end</div>
<div class="p6">
<br /></div>
<div class="p2">
% If data is 3D , mirror image the data in 3D.</div>
<div class="p1">
<span class="s1">if</span> zLen > 1</div>
<div class="p1">
<span class="Apple-converted-space"> </span>Data3DMirrorImage(:,:,1:zLen) = DataMirrorImage;</div>
<div class="p1">
<span class="Apple-converted-space"> </span>Data3DMirrorImage(:,:,(zLen+1):(2*zLen-1)) = flipdim(DataMirrorImage(:,:,2:zLen),3);</div>
<div class="p1">
<span class="Apple-converted-space"> </span>DataMirrorImage = Data3DMirrorImage;</div>
<div class="p5">
end</div>
<br /><br />
Example:<br />
<br />
>> x = randn(1,9)<br />
x =<br />
1.4172 0.6715 -1.2075 0.7172 1.6302 0.4889 1.0347 0.7269 -0.3034<br />
>> xMi = MirrorImageData(x)<br />
xMi =<br />
Columns 1 through 10<br />
1.4172 0.6715 -1.2075 0.7172 1.6302 0.4889 1.0347 0.7269 -0.3034 -0.3034<br />
Columns 11 through 17<br />
0.7269 1.0347 0.4889 1.6302 0.7172 -1.2075 0.6715<br />
<br />
>> x = randn(9);<br />
>> imagesc(x)<br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEh7pMrJKPPq7X4kH2Hq5nj9kF6v88Zly_UrmI78-L_NQFYkrn2Ii9buCFNoKFI56mZozgsx_HPq51OMRciW5kfuV2dG5vPOfGfPnQ_nI_4nTOuPlPGn2-ggBme73Bt7Xo5Yij-oVtHaU3HR/s1600/randImage.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" data-original-height="517" data-original-width="673" height="245" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEh7pMrJKPPq7X4kH2Hq5nj9kF6v88Zly_UrmI78-L_NQFYkrn2Ii9buCFNoKFI56mZozgsx_HPq51OMRciW5kfuV2dG5vPOfGfPnQ_nI_4nTOuPlPGn2-ggBme73Bt7Xo5Yij-oVtHaU3HR/s320/randImage.png" width="320" /></a></div>
<br />
<br />
>> imagesc(MirrorImageData(x))<br />
<br />
<br />
<br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjjBv01Ql_MKufsCun2bLVYYZqo2a74QsMokWgzzRMSdz03T94XXF3-gWLVnH5AQM6d3fQKP1hGR0cqvMbEi-h85BoIz3sZ5y1rncKO2XAj5GUwVWb_57vMXpnuUMLJEGjuzyELmr8zKODc/s1600/randImageMi.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" data-original-height="517" data-original-width="673" height="305" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjjBv01Ql_MKufsCun2bLVYYZqo2a74QsMokWgzzRMSdz03T94XXF3-gWLVnH5AQM6d3fQKP1hGR0cqvMbEi-h85BoIz3sZ5y1rncKO2XAj5GUwVWb_57vMXpnuUMLJEGjuzyELmr8zKODc/s400/randImageMi.png" width="400" /></a></div>
<br />
<br />
<br />Chris Carbonehttp://www.blogger.com/profile/10458001230470047172noreply@blogger.com0tag:blogger.com,1999:blog-7566945490770798840.post-54138724472573120112019-05-25T13:21:00.002-07:002019-05-25T13:21:51.947-07:00'Semi-infinite' Sounds like a lot!Looks like we'll be making computers for a while.
https://www.independent.co.uk/environment/rare-earth-metals-japan-semi-infinite-ocean-mobile-phones-electric-cars-a8301966.htmlChris Carbonehttp://www.blogger.com/profile/10458001230470047172noreply@blogger.com0tag:blogger.com,1999:blog-7566945490770798840.post-81100318172765672042019-05-24T14:50:00.000-07:002019-05-24T14:50:30.415-07:00Random Autocorrelation Sequences R version<!DOCTYPE html>
<html>
<head>
<title>Random Autocorrelation Sequences R version </title>
</head>
<body>
<h3>What is an autocorrelation sequence?</h3>
<p>Autocorrelation sequences (ACSs) are super common when doing anything in probability and statistics. Autocorrelation is a sequence of measurements of how similar a sequence is to it self. In math the autocorrelation sequence r[k] is </p>
<p> r[k] = E[x[n]x[n+k]] for k={0,1,...N-1},</p>
<p>where N is the number of data samples, E is the expected value, x[n] is a data sample and k is the lag. The lag is the separation in samples.</p>
<h3>Why make a random autocorrelation sequence?</h3>
<p>When testing an algorithm or conducting simulations it is often useful to use a random ACS. Generating random a random ACS can be difficult because they have a lot of special properties and if you select a sequence at random, the chance it is a valid ACS is small. </p>
<h3>Trick to making a random autocorrelation sequence</h3>
<p>We can use the following property of ACSs to make generating random ACSs easy. The ACS and the power spectral density (PSD) are Fourier transform (FT) pairs. For our purpose here, a PSD is just a function that is positive everywhere. "FT pair" means the FT of an ACS is a PSD and the inverse FT of a PSD is an ACS. </p>
<p>So we can generate a random ACS using the following steps. First, generate a random sequence. Second, square each element, so the sequence is positive. Finally, find the inverse FT of the squared sequence.</p>
<h3>The R code that produces a random ACS</h3>
<p>The ACS could be any size, but in this case we want a 9 element sequence.</p>
<p>
N <- 9 <br />
PSD <- rnorm(N)^2 <br />
ACS <- fft(PSD,inverse = TRUE) <br />
</p>
<p>The line below outputs the ACS and as you can see it is a complex sequence.</p>
<p>
ACS </p>
<p>
[1] 0.6183715+0.0000000i -0.1375219+0.1960568i -0.1672163-0.2084656i 0.2199730-0.0977208i <br />
[5] -0.0281983+0.2615475i -0.0281983-0.2615475i 0.2199730+0.0977208i -0.1672163+0.2084656i <br />
[9] -0.1375219-0.1960568i <br />
</p>
<h3>What if I want a real ACS</h3>
<p>If you want a real ACS then the PSD has to be even. So, let's make the sequence even!</p>
<p>
PSDeven <- c(PSD,PSD[N:2]) <br />
PSDeven
</p>
[1] 0.39244438 0.03372487 0.69827518 2.54492084 0.10857537 0.67316837 0.23758708 0.54512337<br />
[9] 0.33152416 0.33152416 0.54512337 0.23758708 0.67316837 0.10857537 2.54492084 0.69827518<br />
[17] 0.03372487<br />
<p>
</p>
<p>Notice the ACS is still complex. Numerical error causes some imaginary dust we need to clean up.</p>
<p>
ACS <- fft(PSDeven,inverse = TRUE)/N<br/>
ACS</p>
<p>
[1] 1.1931381+0i 0.1714080+0i -0.3200109+0i -0.5007558+0i -0.2372697+0i 0.2647028+0i<br />
[7] 0.4700409+0i 0.3009945+0i -0.3750372+0i -0.3750372+0i 0.3009945+0i 0.4700409+0i<br />
[13] 0.2647028+0i -0.2372697+0i -0.5007558+0i -0.3200109+0i 0.1714080+0i<br />
</p>
<p>Clean up the small imaginary part with Re() and now we are ready to plot.</p>
<p>
ACS <- Re(ACS)
</p>
<h3> Plot of ACS</h3>
<p>The figure below is a plot of the ACS from lag k = 0 to 16. In textbooks the ACS would have been plotted from k=-8 to 8, with r[0] in the center. </p>
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhHlmRjAg9czdb2tKqK-hsMvwmyeBnIrnHAZZtKDeUpUn47SQmyJDtsJXsNondzNAPDL3Avjelkm_tQNWZHi-gNL63kp_M5mabNbaWOoouQ-LmzMxgOSzs10Vxi500mUl67SQ22SjDWxN6g/s1600/ACSnum.png" imageanchor="1"><img border="0" data-original-height="600" data-original-width="800" height="300" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhHlmRjAg9czdb2tKqK-hsMvwmyeBnIrnHAZZtKDeUpUn47SQmyJDtsJXsNondzNAPDL3Avjelkm_tQNWZHi-gNL63kp_M5mabNbaWOoouQ-LmzMxgOSzs10Vxi500mUl67SQ22SjDWxN6g/s400/ACSnum.png" width="400" /></a>
<p>This is the plot of the ACS in the textbook style. Notice, the lag at 0 r[0] is positive and larger than the other lags, a standard property of ACSs. All is well!</p>
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgNRNzNIDIwsAAnBdyyri1tXGNyboNtGnKUATUXGQz7TYc40RHdZFagA3-I3lalHjg7VRCZOYnVBns3US1v8LFfwBZsDATlf93-eivvb-nernCTDIpy3KE1Pxka3h6-IuF_5-Kvggk384WS/s1600/ACStextbook.png" imageanchor="1"><img border="0" data-original-height="600" data-original-width="800" height="300" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgNRNzNIDIwsAAnBdyyri1tXGNyboNtGnKUATUXGQz7TYc40RHdZFagA3-I3lalHjg7VRCZOYnVBns3US1v8LFfwBZsDATlf93-eivvb-nernCTDIpy3KE1Pxka3h6-IuF_5-Kvggk384WS/s400/ACStextbook.png" width="400" /></a>
</-></body>
</html>
Chris Carbonehttp://www.blogger.com/profile/10458001230470047172noreply@blogger.com0