Shahid Kapoor and Ishaan Khattar hit the gym together
by FilmfareShahid Kapoor and Ishaan Khatter are without a doubt one of the most talented siblings in Bollywood. There isn’t much that we need to say about Shahid Kapoor as the actor has proved his worth by delivering some stunning performances in films like Kaminey, Haider and more recently, Kabir Singh. Even though Ishaan has just started his journey in Bollywood, he’s been appreciated in films like Beyond The Clouds (2017) and Dhadak (2018)
Today, the brothers were spotted by our cameramen after a gym session. While Shahid was spotted wearing a white tee and a pair of black shorts, his brother Ishaan was seen in a two-toned grey hoodie and tracksuits. On the work front, while Shahid will be seen in the Jersey remake, Ishaan will be seen in Mira Nair’s adaptation of Vikram Seth’s A Suitable Boy.
Check out the pictures below.


In setting up a C++ or Python code for variational Monte Carlo calculations, we will use an excellent framework developed by a former Computational Physics student, Morten Ledum, now PhD student at the Hylleraas center for Quantum Chemistry. The GitHub repository is at https://github.com/mortele/variational-monte-carlo-fys4411.
We will discuss this both in our forthcoming lectures and during our lab sessions.
These notebooks serve the aim of linking traditional variational Monte Carlo VMC calculations methods with recent progress on solving many-particle problems using Machine Learning algorithms. The latter is an alternative for project 2. For those who may not be interested in applying machine learning to quantum mechanical many-particle problems, we have the possibility to continue with Monte Carlo methods or opt for other research directions.
In this notebook the aim is to give you an introduction as well as an understanding of the basic elements that are needed in order to develop a professional variational Monte Carlo code. We will focus on a simple system of two particles in an oscillator trap (or alternatively two fermions moving in a Coulombic potential). The particles can interact via a repulsive or an attrative force. It is thus similar to the system described in project 1.
The advantage of these systems is that for two particles (boson or fermions) we have analytical solutions for the eigenpairs of the non-interacting case. Furthermore, for a two- or three-dimensional system of two electrons moving in a harmonic oscillator trap, we have analytical solutions for the interacting case as well.
Having analytical eigenpairs is an invaluable feature that allows us to assess the physical relevance of the trial wave functions, be these either from a standard VMC procedure, from Boltzmann Machines or from Shadow Wave functions.
In this notebook we start with the basics of a VMC calculation and introduce concepts like Markov Chain Monte Carlo methods and the Metropolis algorithm, importance sampling and Metropolis-Hastings algorithm, resampling methods to obtain better estimates of the statistical errors and minimization of the expectation values of the energy and the variance. The latter is done in order to obtain the best possible variational parameters. Furthermore it will define the so-called cost function, a commonly encountered quantity in Machine Learning algorithms. Minimizing the latter is the one which leads to the determination of the optimal parameters in basically all Machine Learning algorithms. This is a possible topic for project 2.
Topics like Markov Chain Monte Carlo and various resampling techniques are also central to Machine Learning methods. Presenting them in the context of VMC approaches leads hopefully to an easier starting point for the understanding of these methods for those of you who will do the machine lerning variant of project 2.
We start with the variational principle. Given a hamiltonian \( H \) and a trial wave function \( \Psi_T(\boldsymbol{R};\boldsymbol{\alpha}) \), the variational principle states that the expectation value of \( \cal{E}[H] \), defined through
$$
\cal {E}[H] =
\frac{\int d\boldsymbol{R}\Psi^{\ast}_T(\boldsymbol{R};\boldsymbol{\alpha})H(\boldsymbol{R})\Psi_T(\boldsymbol{R};\boldsymbol{\alpha})}
{\int d\boldsymbol{R}\Psi^{\ast}_T(\boldsymbol{R};\boldsymbol{\alpha})\Psi_T(\boldsymbol{R};\boldsymbol{\alpha})},
$$
is an upper bound to the ground state energy \( E_0 \) of the hamiltonian \( H \), that is
$$
E_0 \le {\cal E}[H].
$$
In general, the integrals involved in the calculation of various expectation values are multi-dimensional ones. Traditional integration methods such as Gauss-Legendre quadrature will not be adequate for say the computation of the energy of a many-body system.
Here we have defined the vector \( \boldsymbol{R} = [\boldsymbol{r}_1,\boldsymbol{r}_2,\dots,\boldsymbol{r}_n] \) as an array that contains the positions of all particles \( n \) while the vector \( \boldsymbol{\alpha} = [\alpha_1,\alpha_2,\dots,\alpha_m] \) contains the variational parameters of the model, \( m \) in total.
The trial wave function can be expanded in the eigenstates \( \Psi_i(\boldsymbol{R}) \) of the hamiltonian since they form a complete set, viz.,
$$
\Psi_T(\boldsymbol{R};\boldsymbol{\alpha})=\sum_i a_i\Psi_i(\boldsymbol{R}),
$$
and assuming that the set of eigenfunctions are normalized, one obtains
$$
\frac{\sum_{nm}a^*_ma_n \int d\boldsymbol{R}\Psi^{\ast}_m(\boldsymbol{R})H(\boldsymbol{R})\Psi_n(\boldsymbol{R})}
{\sum_{nm}a^*_ma_n \int d\boldsymbol{R}\Psi^{\ast}_m(\boldsymbol{R})\Psi_n(\boldsymbol{R})} =\frac{\sum_{n}a^2_n E_n}
{\sum_{n}a^2_n} \ge E_0,
$$
where we used that \( H(\boldsymbol{R})\Psi_n(\boldsymbol{R})=E_n\Psi_n(\boldsymbol{R}) \). In general, the integrals involved in the calculation of various expectation values are multi-dimensional ones. The variational principle yields the lowest energy of states with a given symmetry.
In most cases, a wave function has only small values in large parts of configuration space, and a straightforward procedure which uses homogenously distributed random points in configuration space will most likely lead to poor results. This may suggest that some kind of importance sampling combined with e.g., the Metropolis algorithm may be a more efficient way of obtaining the ground state energy. The hope is then that those regions of configurations space where the wave function assumes appreciable values are sampled more efficiently.
The tedious part in a VMC calculation is the search for the variational minimum. A good knowledge of the system is required in order to carry out reasonable VMC calculations. This is not always the case, and often VMC calculations serve rather as the starting point for so-called diffusion Monte Carlo calculations (DMC). Diffusion Monte Carlo is a way of solving exactly the many-body Schroedinger equation by means of a stochastic procedure. A good guess on the binding energy and its wave function is however necessary. A carefully performed VMC calculation can aid in this context.
The basic procedure of a Variational Monte Carlo calculations consists thus of
$$
\overline{E}[\boldsymbol{\alpha}]=\frac{\int d\boldsymbol{R}\Psi^{\ast}_{T}(\boldsymbol{R},\boldsymbol{\alpha})H(\boldsymbol{R})\Psi_{T}(\boldsymbol{R},\boldsymbol{\alpha})}
{\int d\boldsymbol{R}\Psi^{\ast}_{T}(\boldsymbol{R},\boldsymbol{\alpha})\Psi_{T}(\boldsymbol{R},\boldsymbol{\alpha})}.
$$
Here we have used the notation \( \overline{E} \) to label the expectation value of the energy.
In order to bring in the Monte Carlo machinery, we define first a likelihood distribution, or probability density distribution (PDF). Using our ansatz for the trial wave function \( \psi_T(\boldsymbol{R};\boldsymbol{\alpha}) \) we define a PDF
$$
P(\boldsymbol{R})= \frac{\left|\psi_T(\boldsymbol{R};\boldsymbol{\alpha})\right|^2}{\int \left|\psi_T(\boldsymbol{R};\boldsymbol{\alpha})\right|^2d\boldsymbol{R}}.
$$
This is our model for probability distribution function. The approximation to the expectation value of the Hamiltonian is now
$$
\overline{E}[\boldsymbol{\alpha}] =
\frac{\int d\boldsymbol{R}\Psi^{\ast}_T(\boldsymbol{R};\boldsymbol{\alpha})H(\boldsymbol{R})\Psi_T(\boldsymbol{R};\boldsymbol{\alpha})}
{\int d\boldsymbol{R}\Psi^{\ast}_T(\boldsymbol{R};\boldsymbol{\alpha})\Psi_T(\boldsymbol{R};\boldsymbol{\alpha})}.
$$
We define a new quantity
$$
E_L(\boldsymbol{R};\boldsymbol{\alpha})=\frac{1}{\psi_T(\boldsymbol{R};\boldsymbol{\alpha})}H\psi_T(\boldsymbol{R};\boldsymbol{\alpha}),
\tag{1}
$$
called the local energy, which, together with our trial PDF yields a new expression (and which look simlar to the the expressions for moments in statistics)
$$
\overline{E}[\boldsymbol{\alpha}]=\int P(\boldsymbol{R})E_L(\boldsymbol{R};\boldsymbol{\alpha}) d\boldsymbol{R}\approx \frac{1}{N}\sum_{i=1}^NE_L(\boldsymbol{R_i};\boldsymbol{\alpha})
\tag{2}
$$
with \( N \) being the number of Monte Carlo samples. The expression on the right hand side follows from Bernoulli's law of large numbers, which states that the sample mean, in the limit \( N\rightarrow \infty \) approaches the true mean
The Algorithm for performing a variational Monte Carlo calculations runs as this
Observe that the jumping in space is governed by the variable step. This is often referred to as the brute-force sampling and is normally replaced by what is called importance sampling, discussed in more detail next week..
The radial Schroedinger equation for the hydrogen atom can be written as (when we have gotten rid of the first derivative term in the kinetic energy and used \( rR(r)=u(r) \))
$$
-\frac{\hbar^2}{2m}\frac{d^2 u(r)}{d r^2}-
\left(\frac{ke^2}{r}-\frac{\hbar^2l(l+1)}{2mr^2}\right)u(r)=Eu(r).
$$
We will specialize to the case with \( l=0 \) and end up with
$$
-\frac{\hbar^2}{2m}\frac{d^2 u(r)}{d r^2}-
\left(\frac{ke^2}{r}\right)u(r)=Eu(r).
$$
Then we introduce a dimensionless variable \( \rho=r/a \) where \( a \) is a constant with dimension length. Multiplying with \( ma^2/\hbar^2 \) we can rewrite our equations as
$$
-\frac{1}{2}\frac{d^2 u(\rho)}{d \rho^2}-
\frac{ke^2ma}{\hbar^2}\frac{u(\rho)}{\rho}-\lambda u(\rho)=0.
$$
Since \( a \) is just a parameter we choose to set
$$
\frac{ke^2ma}{\hbar^2}=1,
$$
which leads to \( a=\hbar^2/mke^2 \), better known as the Bohr radius with value \( 0.053 \) nm. Scaling the equations this way does not only render our numerical treatment simpler since we avoid carrying with us all physical parameters, but we obtain also a natural length scale. We will see this again and again. In our discussions below with a harmonic oscillator trap, the natural lentgh scale with be determined by the oscillator frequency, the mass of the particle and \( \hbar \). We have also defined a dimensionless 'energy' \( \lambda = Ema^2/\hbar^2 \). With the rescaled quantities, the ground state energy of the hydrogen atom is \( 1/2 \). The equation we want to solve is now defined by the Hamiltonian
$$
H=-\frac{1}{2}\frac{d^2 }{d \rho^2}-\frac{1}{\rho}.
$$
As trial wave function we peep now into the analytical solution for the hydrogen atom and use (with \( \alpha \) as a variational parameter)
$$
u_T^{\alpha}(\rho)=\alpha\rho \exp{-(\alpha\rho)}.
$$
Inserting this wave function into the expression for the local energy \( E_L \) gives
$$
E_L(\rho)=-\frac{1}{\rho}-
\frac{\alpha}{2}\left(\alpha-\frac{2}{\rho}\right).
$$
To have analytical local energies saves us from computing numerically the second derivative, a feature which often increases our numerical expenditure with a factor of three or more. Integratng up the local energy (recall to bring back the PDF in the integration) gives \( \overline{E}[\boldsymbol{\alpha}]=\alpha(\alpha/2-1) \).
We present here another well-known example, the harmonic oscillator in one dimension for one particle. This will also serve the aim of introducing our next model, namely that of interacting electrons in a harmonic oscillator trap.
Here as well, we do have analytical solutions and the energy of the ground state, with \( \hbar=1 \), is \( 1/2\omega \), with \( \omega \) being the oscillator frequency. We use the following trial wave function
$$
\psi_T(x;\alpha) = \exp{-(\frac{1}{2}\alpha^2x^2)},
$$
which results in a local energy
$$
\frac{1}{2}\left(\alpha^2+x^2(1-\alpha^4)\right).
$$
We can compare our numerically calculated energies with the exact energy as function of \( \alpha \)
$$
\overline{E}[\alpha] = \frac{1}{4}\left(\alpha^2+\frac{1}{\alpha^2}\right).
$$
Similarly, with the above ansatz, we can also compute the exact variance which reads
$$
\sigma^2[\alpha]=\frac{1}{4}\left(1+(1-\alpha^4)^2\frac{3}{4\alpha^4}\right)-\overline{E}.
$$
Our code for computing the energy of the ground state of the harmonic oscillator follows here. We start by defining directories where we store various outputs.
# Common imports
import os
# Where to save the figures and data files
PROJECT_ROOT_DIR = "Results"
FIGURE_ID = "Results/FigureFiles"
DATA_ID = "Results/VMCHarmonic"
if not os.path.exists(PROJECT_ROOT_DIR):
os.mkdir(PROJECT_ROOT_DIR)
if not os.path.exists(FIGURE_ID):
os.makedirs(FIGURE_ID)
if not os.path.exists(DATA_ID):
os.makedirs(DATA_ID)
def image_path(fig_id):
return os.path.join(FIGURE_ID, fig_id)
def data_path(dat_id):
return os.path.join(DATA_ID, dat_id)
def save_fig(fig_id):
plt.savefig(image_path(fig_id) + ".png", format='png')
outfile = open(data_path("VMCHarmonic.dat"),'w')
We proceed with the implementation of the Monte Carlo algorithm but list first the ansatz for the wave function and the expression for the local energy
# VMC for the one-dimensional harmonic oscillator
# Brute force Metropolis, no importance sampling and no energy minimization
from math import exp, sqrt
from random import random, seed
import numpy as np
import matplotlib.pyplot as plt
from decimal import *
# Trial wave function for the Harmonic oscillator in one dimension
def WaveFunction(r,alpha):
return exp(-0.5*alpha*alpha*r*r)
# Local energy for the Harmonic oscillator in one dimension
def LocalEnergy(r,alpha):
return 0.5*r*r*(1-alpha**4) + 0.5*alpha*alpha
Note that in the Metropolis algorithm there is no need to compute the trial wave function, mainly since we are just taking the ratio of two exponentials. It is then from a computational point view, more convenient to compute the argument from the ratio and then calculate the exponential. Here we have refrained from this purely of pedagogical reasons.
# The Monte Carlo sampling with the Metropolis algo
def MonteCarloSampling():
NumberMCcycles= 100000
StepSize = 1.0
# positions
PositionOld = 0.0
PositionNew = 0.0
# seed for rng generator
seed()
# start variational parameter
alpha = 0.4
for ia in range(MaxVariations):
alpha += .05
AlphaValues[ia] = alpha
energy = energy2 = 0.0
#Initial position
PositionOld = StepSize * (random() - .5)
wfold = WaveFunction(PositionOld,alpha)
#Loop over MC MCcycles
for MCcycle in range(NumberMCcycles):
#Trial position
PositionNew = PositionOld + StepSize*(random() - .5)
wfnew = WaveFunction(PositionNew,alpha)
#Metropolis test to see whether we accept the move
if random() <= wfnew**2 / wfold**2:
PositionOld = PositionNew
wfold = wfnew
DeltaE = LocalEnergy(PositionOld,alpha)
energy += DeltaE
energy2 += DeltaE**2
#We calculate mean, variance and error
energy /= NumberMCcycles
energy2 /= NumberMCcycles
variance = energy2 - energy**2
error = sqrt(variance/NumberMCcycles)
Energies[ia] = energy
Variances[ia] = variance
outfile.write('%f %f %f %f \n' %(alpha,energy,variance,error))
return Energies, AlphaValues, Variances
Finally, the results are presented here with the exact energies and variances as well.
#Here starts the main program with variable declarations
MaxVariations = 20
Energies = np.zeros((MaxVariations))
ExactEnergies = np.zeros((MaxVariations))
ExactVariance = np.zeros((MaxVariations))
Variances = np.zeros((MaxVariations))
AlphaValues = np.zeros(MaxVariations)
(Energies, AlphaValues, Variances) = MonteCarloSampling()
outfile.close()
ExactEnergies = 0.25*(AlphaValues*AlphaValues+1.0/(AlphaValues*AlphaValues))
ExactVariance = 0.25*(1.0+((1.0-AlphaValues**4)**2)*3.0/(4*(AlphaValues**4)))-ExactEnergies*ExactEnergies
#simple subplot
plt.subplot(2, 1, 1)
plt.plot(AlphaValues, Energies, 'o-',AlphaValues, ExactEnergies,'r-')
plt.title('Energy and variance')
plt.ylabel('Dimensionless energy')
plt.subplot(2, 1, 2)
plt.plot(AlphaValues, Variances, '.-',AlphaValues, ExactVariance,'r-')
plt.xlabel(r'$\alpha$', fontsize=15)
plt.ylabel('Variance')
save_fig("VMCHarmonic")
plt.show()
#nice printout with Pandas
import pandas as pd
from pandas import DataFrame
data ={'Alpha':AlphaValues, 'Energy':Energies,'Exact Energy':ExactEnergies,'Variance':Variances,'Exact Variance':ExactVariance,}
frame = pd.DataFrame(data)
print(frame)
For \( \alpha=1 \) we have the exact eigenpairs, as can be deduced from the table here. With \( \omega=1 \), the exact energy is \( 1/2 \) a.u. with zero variance, as it should. We see also that our computed variance follows rather well the exact variance. Increasing the number of Monte Carlo cycles will improve our statistics (try to increase the number of Monte Carlo cycles).
The fact that the variance is exactly equal to zero when \( \alpha=1 \) is that we then have the exact wave function, and the action of the hamiltionan on the wave function
$$
H\psi = \mathrm{constant}\times \psi,
$$
yields just a constant. The integral which defines various expectation values involving moments of the hamiltonian becomes then
$$
\langle H^n \rangle =
\frac{\int d\boldsymbol{R}\Psi^{\ast}_T(\boldsymbol{R})H^n(\boldsymbol{R})\Psi_T(\boldsymbol{R})}
{\int d\boldsymbol{R}\Psi^{\ast}_T(\boldsymbol{R})\Psi_T(\boldsymbol{R})}=
\mathrm{constant}\times\frac{\int d\boldsymbol{R}\Psi^{\ast}_T(\boldsymbol{R})\Psi_T(\boldsymbol{R})}
{\int d\boldsymbol{R}\Psi^{\ast}_T(\boldsymbol{R})\Psi_T(\boldsymbol{R})}=\mathrm{constant}.
$$
This gives an important information: the exact wave function leads to zero variance!
As we will see below, many practitioners perform a minimization on both the energy and the variance.
Till now we have not yet discussed the derivation of the Metropolis algorithm. We assume the reader has some familiarity with the mathematics of Markov chains.
The Metropolis algorithm , see the original article, was invented by Metropolis et. al and is often simply called the Metropolis algorithm. It is a method to sample a normalized probability distribution by a stochastic process. We define \( \mathbf{P}_i^{(n)} \) to be the probability for finding the system in the state \( i \) at step \( n \). The algorithm is then
We wish to derive the required properties of \( T \) and \( A \) such that \( \mathbf{P}_i^{(n\rightarrow \infty)} \rightarrow p_i \) so that starting from any distribution, the method converges to the correct distribution. Note that the description here is for a discrete probability distribution. Replacing probabilities \( p_i \) with expressions like \( p(x_i)dx_i \) will take all of these over to the corresponding continuum expressions.
The dynamical equation for \( \mathbf{P}_i^{(n)} \) can be written directly from the description above. The probability of being in the state \( i \) at step \( n \) is given by the probability of being in any state \( j \) at the previous step, and making an accepted transition to \( i \) added to the probability of being in the state \( i \), making a transition to any state \( j \) and rejecting the move:
$$
\mathbf{P}^{(n)}_i = \sum_j \left [
\mathbf{P}^{(n-1)}_jT_{j\rightarrow i} A_{j\rightarrow i}
+\mathbf{P}^{(n-1)}_iT_{i\rightarrow j}\left ( 1- A_{i\rightarrow j} \right)
\right ] \,.
$$
Since the probability of making some transition must be 1, \( \sum_j T_{i\rightarrow j} = 1 \), and the above equation becomes
$$
\mathbf{P}^{(n)}_i = \mathbf{P}^{(n-1)}_i +
\sum_j \left [
\mathbf{P}^{(n-1)}_jT_{j\rightarrow i} A_{j\rightarrow i}
-\mathbf{P}^{(n-1)}_iT_{i\rightarrow j}A_{i\rightarrow j}
\right ] \,.
$$
For large \( n \) we require that \( \mathbf{P}^{(n\rightarrow \infty)}_i = p_i \), the desired probability distribution. Taking this limit, gives the balance requirement
$$
\sum_j \left [
p_jT_{j\rightarrow i} A_{j\rightarrow i}
-p_iT_{i\rightarrow j}A_{i\rightarrow j}
\right ] = 0 \,.
$$
The balance requirement is very weak. Typically the much stronger detailed balance requirement is enforced, that is rather than the sum being set to zero, we set each term separately to zero and use this to determine the acceptance probabilities. Rearranging, the result is
$$
\frac{ A_{j\rightarrow i}}{A_{i\rightarrow j}}
= \frac{p_iT_{i\rightarrow j}}{ p_jT_{j\rightarrow i}} \,.
$$
The Metropolis choice is to maximize the \( A \) values, that is
$$
A_{j \rightarrow i} = \min \left ( 1,
\frac{p_iT_{i\rightarrow j}}{ p_jT_{j\rightarrow i}}\right ).
$$
Other choices are possible, but they all correspond to multilplying \( A_{i\rightarrow j} \) and \( A_{j\rightarrow i} \) by the same constant smaller than unity.\footnote{The penalty function method uses just such a factor to compensate for \( p_i \) that are evaluated stochastically and are therefore noisy.}
Having chosen the acceptance probabilities, we have guaranteed that if the \( \mathbf{P}_i^{(n)} \) has equilibrated, that is if it is equal to \( p_i \), it will remain equilibrated. Next we need to find the circumstances for convergence to equilibrium.
The dynamical equation can be written as
$$
\mathbf{P}^{(n)}_i = \sum_j M_{ij}\mathbf{P}^{(n-1)}_j
$$
with the matrix \( M \) given by
$$
M_{ij} = \delta_{ij}\left [ 1 -\sum_k T_{i\rightarrow k} A_{i \rightarrow k}
\right ] + T_{j\rightarrow i} A_{j\rightarrow i} \,.
$$
Summing over \( i \) shows that \( \sum_i M_{ij} = 1 \), and since \( \sum_k T_{i\rightarrow k} = 1 \), and \( A_{i \rightarrow k} \leq 1 \), the elements of the matrix satisfy \( M_{ij} \geq 0 \). The matrix \( M \) is therefore a stochastic matrix.
The Metropolis method is simply the power method for computing the right eigenvector of \( M \) with the largest magnitude eigenvalue. By construction, the correct probability distribution is a right eigenvector with eigenvalue 1. Therefore, for the Metropolis method to converge to this result, we must show that \( M \) has only one eigenvalue with this magnitude, and all other eigenvalues are smaller.
The Hamiltonian of the quantum dot is given by
$$ \hat{H} = \hat{H}_0 + \hat{V},
$$
where \( \hat{H}_0 \) is the many-body HO Hamiltonian, and \( \hat{V} \) is the inter-electron Coulomb interactions. In dimensionless units,
$$ \hat{V}= \sum_{i < j}^N \frac{1}{r_{ij}},
$$
with \( r_{ij}=\sqrt{\mathbf{r}_i^2 - \mathbf{r}_j^2} \).
This leads to the separable Hamiltonian, with the relative motion part given by (\( r_{ij}=r \))
$$
\hat{H}_r=-\nabla^2_r + \frac{1}{4}\omega^2r^2+ \frac{1}{r},
$$
plus a standard Harmonic Oscillator problem for the center-of-mass motion. This system has analytical solutions in two and three dimensions (M. Taut 1993 and 1994).
We want to perform a Variational Monte Carlo calculation of the ground state of two electrons in a quantum dot well with different oscillator energies, assuming total spin \( S=0 \). Our trial wave function has the following form
$$
\begin{equation}
\psi_{T}(\boldsymbol{r}_1,\boldsymbol{r}_2) =
C\exp{\left(-\alpha_1\omega(r_1^2+r_2^2)/2\right)}
\exp{\left(\frac{r_{12}}{(1+\alpha_2 r_{12})}\right)},
\tag{3}
\end{equation}
$$
where the $\alpha$s represent our variational parameters, two in this case.
Why does the trial function look like this? How did we get there? This will be one of our main motivations for switching to Machine Learning later.
To find an ansatz for the correlated part of the wave function, it is useful to rewrite the two-particle local energy in terms of the relative and center-of-mass motion. Let us denote the distance between the two electrons as \( r_{12} \). We omit the center-of-mass motion since we are only interested in the case when \( r_{12} \rightarrow 0 \). The contribution from the center-of-mass (CoM) variable \( \boldsymbol{R}_{\mathrm{CoM}} \) gives only a finite contribution. We focus only on the terms that are relevant for \( r_{12} \) and for three dimensions.
The relevant local energy becomes then
$$ \lim_{r_{12} \rightarrow 0}E_L(R)= \frac{1}{{\cal R}_T(r_{12})}\left(2\frac{d^2}{dr_{ij}^2}+\frac{4}{r_{ij}}\frac{d}{dr_{ij}}+\frac{2}{r_{ij}}-\frac{l(l+1)}{r_{ij}^2}+2E \right){\cal R}_T(r_{12})
= 0.
$$
Set \( l=0 \) and we have the so-called cusp condition
$$ \frac{d {\cal R}_T(r_{12})}{dr_{12}} = -\frac{1}{2(l+1)} {\cal R}_T(r_{12})\qquad r_{12}\to 0
$$
The above results in
$$
{\cal R}_T \propto \exp{(r_{ij}/2)},
$$
for anti-parallel spins and
$$
{\cal R}_T \propto \exp{(r_{ij}/4)},
$$
for anti-parallel spins. This is the so-called cusp condition for the relative motion, resulting in a minimal requirement for the correlation part of the wave fuction. For general systems containing more than say two electrons, we have this condition for each electron pair \( ij \).
First, as with the hydrogen case, we declare where to store files.
# Common imports
import os
# Where to save the figures and data files
PROJECT_ROOT_DIR = "Results"
FIGURE_ID = "Results/FigureFiles"
DATA_ID = "Results/VMCQdotMetropolis"
if not os.path.exists(PROJECT_ROOT_DIR):
os.mkdir(PROJECT_ROOT_DIR)
if not os.path.exists(FIGURE_ID):
os.makedirs(FIGURE_ID)
if not os.path.exists(DATA_ID):
os.makedirs(DATA_ID)
def image_path(fig_id):
return os.path.join(FIGURE_ID, fig_id)
def data_path(dat_id):
return os.path.join(DATA_ID, dat_id)
def save_fig(fig_id):
plt.savefig(image_path(fig_id) + ".png", format='png')
outfile = open(data_path("VMCQdotMetropolis.dat"),'w')
Thereafter we set up the analytical expressions for the wave functions and the local energy
# 2-electron VMC for quantum dot system in two dimensions
# Brute force Metropolis, no importance sampling and no energy minimization
from math import exp, sqrt
from random import random, seed
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
from matplotlib.ticker import LinearLocator, FormatStrFormatter
import sys
# Trial wave function for the 2-electron quantum dot in two dims
def WaveFunction(r,alpha,beta):
r1 = r[0,0]**2 + r[0,1]**2
r2 = r[1,0]**2 + r[1,1]**2
r12 = sqrt((r[0,0]-r[1,0])**2 + (r[0,1]-r[1,1])**2)
deno = r12/(1+beta*r12)
return exp(-0.5*alpha*(r1+r2)+deno)
# Local energy for the 2-electron quantum dot in two dims, using analytical local energy
def LocalEnergy(r,alpha,beta):
r1 = (r[0,0]**2 + r[0,1]**2)
r2 = (r[1,0]**2 + r[1,1]**2)
r12 = sqrt((r[0,0]-r[1,0])**2 + (r[0,1]-r[1,1])**2)
deno = 1.0/(1+beta*r12)
deno2 = deno*deno
return 0.5*(1-alpha*alpha)*(r1 + r2) +2.0*alpha + 1.0/r12+deno2*(alpha*r12-deno2+2*beta*deno-1.0/r12)
The Monte Carlo sampling without importance sampling is set up here.
# The Monte Carlo sampling with the Metropolis algo
def MonteCarloSampling():
NumberMCcycles= 10000
StepSize = 1.0
# positions
PositionOld = np.zeros((NumberParticles,Dimension), np.double)
PositionNew = np.zeros((NumberParticles,Dimension), np.double)
# seed for rng generator
seed()
# start variational parameter
alpha = 0.9
for ia in range(MaxVariations):
alpha += .025
AlphaValues[ia] = alpha
beta = 0.2
for jb in range(MaxVariations):
beta += .01
BetaValues[jb] = beta
energy = energy2 = 0.0
DeltaE = 0.0
#Initial position
for i in range(NumberParticles):
for j in range(Dimension):
PositionOld[i,j] = StepSize * (random() - .5)
wfold = WaveFunction(PositionOld,alpha,beta)
#Loop over MC MCcycles
for MCcycle in range(NumberMCcycles):
#Trial position moving one particle at the time
for i in range(NumberParticles):
for j in range(Dimension):
PositionNew[i,j] = PositionOld[i,j] + StepSize * (random() - .5)
wfnew = WaveFunction(PositionNew,alpha,beta)
#Metropolis test to see whether we accept the move
if random() < wfnew**2 / wfold**2:
for j in range(Dimension):
PositionOld[i,j] = PositionNew[i,j]
wfold = wfnew
DeltaE = LocalEnergy(PositionOld,alpha,beta)
energy += DeltaE
energy2 += DeltaE**2
#We calculate mean, variance and error ...
energy /= NumberMCcycles
energy2 /= NumberMCcycles
variance = energy2 - energy**2
error = sqrt(variance/NumberMCcycles)
Energies[ia,jb] = energy
Variances[ia,jb] = variance
outfile.write('%f %f %f %f %f\n' %(alpha,beta,energy,variance,error))
return Energies, Variances, AlphaValues, BetaValues
And finally comes the main part with the plots as well.
#Here starts the main program with variable declarations
NumberParticles = 2
Dimension = 2
MaxVariations = 10
Energies = np.zeros((MaxVariations,MaxVariations))
Variances = np.zeros((MaxVariations,MaxVariations))
AlphaValues = np.zeros(MaxVariations)
BetaValues = np.zeros(MaxVariations)
(Energies, Variances, AlphaValues, BetaValues) = MonteCarloSampling()
outfile.close()
# Prepare for plots
fig = plt.figure()
ax = fig.gca(projection='3d')
# Plot the surface.
X, Y = np.meshgrid(AlphaValues, BetaValues)
surf = ax.plot_surface(X, Y, Energies,cmap=cm.coolwarm,linewidth=0, antialiased=False)
# Customize the z axis.
zmin = np.matrix(Energies).min()
zmax = np.matrix(Energies).max()
ax.set_zlim(zmin, zmax)
ax.set_xlabel(r'$\alpha$')
ax.set_ylabel(r'$\beta$')
ax.set_zlabel(r'$\langle E \rangle$')
ax.zaxis.set_major_locator(LinearLocator(10))
ax.zaxis.set_major_formatter(FormatStrFormatter('%.02f'))
# Add a color bar which maps values to colors.
fig.colorbar(surf, shrink=0.5, aspect=5)
save_fig("QdotMetropolis")
plt.show()


The Consortium is open to donors and entities engaged in international development work. Members have access to the Data Collaboratives data sharing platform and IT infrastructure and are invited to participate in workshops, exchanges, and training activities.
The Consortium helps partners unlock public good opportunities from their data in a secure, responsible manner. Data partnerships can be leveraged to open markets in emerging economies, receive new data methods and algorithms, and increase staff skills through collaboration and secondment opportunities.


Sniffy autodetects all outgoing network connections and shows them in tables below.
You can trigger java.net.ConnectException for any of them using switches in the right column.
If you want to test the startup of your application with no connectivity, toggle the Keep the fault tolerance settings after server restart switch and restart your server.
Powered by Sniffy.
| Host | Port | Status | Delay |
|---|


Wrapper around gm (ImageMagick) API's selected parts.
gm
's API is two-fold: getters like size
, identify
… are one-shot calls,
whereas other functions are chainable calls whose purpose is to build a
convert
command's options sequence; that chain ends up with a call to
write
.
Both getters and write
use the error-first callback style. The goal of this
wrapper is to provide a promise-compliant API, as well as exposing utility
functions.
- Source:
- See:
Classes
Methods
(inner) buildPromiseAPI()
Extend a Magick
instance with a promise-compliant API.
- Source:
(inner) clamp(value, min, max)
Cap a number to a range.
Parameters:
| Name | Type | Description |
|---|---|---|
value |
Number | A number to be limited to the output range |
min |
Number | The lower boundary of the output range |
max |
Number | The upper boundary of the output range |
- Source:
Returns:
A number in the range [min, max]


This is a markdown formatted text.
Some lorem ipsum here.
This is a global value, which is available to all pages
Apple, Orange, & Pears
Last Updated: Thu Aug 02 2012


"""Build wheels/sdists by installing build deps to a temporary environment. """ import os import logging import toml import shutil from subprocess import check_call import sys from sysconfig import get_paths from tempfile import mkdtemp from .wrappers import Pep517HookCaller, LoggerWrapper log = logging.getLogger(__name__) def _load_pyproject(source_dir): with open(os.path.join(source_dir, 'pyproject.toml')) as f: pyproject_data = toml.load(f) buildsys = pyproject_data['build-system'] return ( buildsys['requires'], buildsys['build-backend'], buildsys.get('backend-path'), ) class BuildEnvironment(object): """Context manager to install build deps in a simple temporary environment Based on code I wrote for pip, which is MIT licensed. """ # Copyright (c) 2008-2016 The pip developers (see AUTHORS.txt file) # # Permission is hereby granted, free of charge, to any person obtaining # a copy of this software and associated documentation files (the # "Software"), to deal in the Software without restriction, including # without limitation the rights to use, copy, modify, merge, publish, # distribute, sublicense, and/or sell copies of the Software, and to # permit persons to whom the Software is furnished to do so, subject to # the following conditions: # # The above copyright notice and this permission notice shall be # included in all copies or substantial portions of the Software. # # THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, # EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF # MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND # NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE # LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION # OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION # WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. path = None def __init__(self, cleanup=True): self._cleanup = cleanup def __enter__(self): self.path = mkdtemp(prefix='pep517-build-env-') log.info('Temporary build environment: %s', self.path) self.save_path = os.environ.get('PATH', None) self.save_pythonpath = os.environ.get('PYTHONPATH', None) install_scheme = 'nt' if (os.name == 'nt') else 'posix_prefix' install_dirs = get_paths(install_scheme, vars={ 'base': self.path, 'platbase': self.path, }) scripts = install_dirs['scripts'] if self.save_path: os.environ['PATH'] = scripts + os.pathsep + self.save_path else: os.environ['PATH'] = scripts + os.pathsep + os.defpath if install_dirs['purelib'] == install_dirs['platlib']: lib_dirs = install_dirs['purelib'] else: lib_dirs = install_dirs['purelib'] + os.pathsep + \ install_dirs['platlib'] if self.save_pythonpath: os.environ['PYTHONPATH'] = lib_dirs + os.pathsep + \ self.save_pythonpath else: os.environ['PYTHONPATH'] = lib_dirs return self def pip_install(self, reqs): """Install dependencies into this env by calling pip in a subprocess""" if not reqs: return log.info('Calling pip to install %s', reqs) cmd = [ sys.executable, '-m', 'pip', 'install', '--ignore-installed', '--prefix', self.path] + list(reqs) check_call( cmd, stdout=LoggerWrapper(log, logging.INFO), stderr=LoggerWrapper(log, logging.ERROR), ) def __exit__(self, exc_type, exc_val, exc_tb): needs_cleanup = ( self._cleanup and self.path is not None and os.path.isdir(self.path) ) if needs_cleanup: shutil.rmtree(self.path) if self.save_path is None: os.environ.pop('PATH', None) else: os.environ['PATH'] = self.save_path if self.save_pythonpath is None: os.environ.pop('PYTHONPATH', None) else: os.environ['PYTHONPATH'] = self.save_pythonpath def build_wheel(source_dir, wheel_dir, config_settings=None): """Build a wheel from a source directory using PEP 517 hooks. :param str source_dir: Source directory containing pyproject.toml :param str wheel_dir: Target directory to create wheel in :param dict config_settings: Options to pass to build backend This is a blocking function which will run pip in a subprocess to install build requirements. """ if config_settings is None: config_settings = {} requires, backend, backend_path = _load_pyproject(source_dir) hooks = Pep517HookCaller(source_dir, backend, backend_path) with BuildEnvironment() as env: env.pip_install(requires) reqs = hooks.get_requires_for_build_wheel(config_settings) env.pip_install(reqs) return hooks.build_wheel(wheel_dir, config_settings) def build_sdist(source_dir, sdist_dir, config_settings=None): """Build an sdist from a source directory using PEP 517 hooks. :param str source_dir: Source directory containing pyproject.toml :param str sdist_dir: Target directory to place sdist in :param dict config_settings: Options to pass to build backend This is a blocking function which will run pip in a subprocess to install build requirements. """ if config_settings is None: config_settings = {} requires, backend, backend_path = _load_pyproject(source_dir) hooks = Pep517HookCaller(source_dir, backend, backend_path) with BuildEnvironment() as env: env.pip_install(requires) reqs = hooks.get_requires_for_build_sdist(config_settings) env.pip_install(reqs) return hooks.build_sdist(sdist_dir, config_settings)


 Recently its Author has passed away and he’s offered his position including the ownership of OASIS to anybody who masters his three challenges:
- Win an unwinnable race. The key to this challenge is to drive backwards.
- Ask the Lady of the Author to dance. Easy, isn’t it? Well, the Lady is trapped in the “Overlook Hotel” from the horror movie “The Shinning” (sa:deva) and she’s dancing with dead people. On the way to the Lady, one has to go through all the zombies. Lady is source of life itself, in the form of a woman. And key to it is through suffering. Hint from the 2nd challenge: “Ultimate answer will be found in the Fortress Tragic”.
- The final challenge takes place on the planet Doom. But this is not the challenge. It is to win an unknown game while playing it! How to do? Just like life, really. “The key is to play without going for the win. Everyone who wins loses. Keys are invisible and hidden in a dark room that’s the center of the maze. Found by wandering about until you get to the invisible dot - it’s pretty easy if you know what you’re doing…” And the name of the first game? “Adventure”. And the prize obtained is the first easter egg ever - the name of the Author of the game.
After mastering the three challenges the winner - Parzival - is presented with a huge golden “easter egg”. And what does the winner gain? Well, the usual - the world with the princess and his new identity as the ruler of it - of course.
–
- Halliday’s avatar sits on a throne out of which snakes arise.
- Planet doom gets surrounded by impenetrable force field. This gets deactivated by his love(r).
- (In the end fight in OASIS) “This is my world, Sorento.”
- “You gotta go straight to the invisible dot at the start.”
- Transformation of Halliday in the end: Avatar -> Real. He acknowledges he’s really dead.
- The Easter Egg given, full of light, was well hidden in the room of Halliday’s childhood.
- “Once you take the egg the game will be over”.
- “I just didn’t know how to connect with people there. I was afraid for all my life. Right up until the day I knew my life was ending. And that was when I realized that as terrifying and painful as reality can be, it’s also the only place where you get a decent meal.”
- “Reality is the only thing… that’s real. "


About Zante
Zakynthos (A.K.A Zante) is a Greek island in the Ionian Sea and a well-known summer resort. The harbor city of Zakynthos is the capital and major hub, centered around waterfront Solomos Square. Popular beaches like Agios Nikolaos, Alykanas and Tsilivi offer swimming and water sports. Accessed by boat, Navagio beach is the site of a famed 1980 shipwreck resting in a sandy cove framed by cliffs
About Paris
Paris, France's capital, is a major European city and a global center for art, fashion, gastronomy and culture. Its 19th-century cityscape is crisscrossed by wide boulevards and the River Seine. Beyond such landmarks as the Eiffel Tower and the 12th-century, Gothic Notre-Dame cathedral, the city is known for its cafe culture and designer boutiques along the Rue du Faubourg Saint-Honoré.
My name is Jane and I am an avid traveller and explorer and love experiencing new cultures and languages!
Amsterdam
London
Albufeira
Venice
Barcelona


--- title: Time Series Introduction with R codes author: Francisco Bischoff date: '2021-06-16' categories: - R - r-bloggers - timeseries tags: - R slug: time-series-introduction-with-r-codes bibliography: references.yaml csl: springer-lncs.csl mathjax: yes editor_options: markdown: mode: markdown output: blogdown::html_page: number_sections: yes ---
A set of observed values ordered in time, or we can say, repeated measurement of something usually with the same fixed interval of time (hourly, weekly, monthly).
A collection of observations made sequentially in time[1].
If the variable we are measuring is a count variable, we may have a Poisson Time Series (that is for later).
A time series \(T \in \mathbb{R}^n\) is a sequence of real-valued numbers \(t_i \in \mathbb{R} : T=[t_1,t_2,\dots,t_n]\) where \(n\) is the length of \(T\).
Most of the classic statistical theory is based on the assumption of sample randomness and independent observations. On the other hand, time series is just the opposite. Observations are usually dependent on previous values, and their analysis must take into account their temporal order.
For example, a prospective cohort study comparing “injury rate before and after” an implemented program, analyses of time trends, such as Poisson regression and time series analysis, considers the variability that occurs over the study period apart from the change associated with the intervention. They also avoid the loss of information about variability in incidence over time when rates are aggregated into one before and one after rate. The population is its own control.
If previous observations can predict future observations exactly, we have a deterministic process. However, the exact prediction is usually impossible since past values determine only part of the future value. In this case, we say the process is stochastic, and the future values must be seen as a probability conditioned on the past values.
There are many models available to describe the behavior of a particular series. The choice of such a model depends on factors such as the behavior of the phenomenon or the prior knowledge of its nature and the purpose of the analysis.
The R base installation already gives us lots of datasets to work on time-series.
For this article I’ll first load the MASS
package that contains some of the dataset we will use.
We can list all datasets available with the function data()
from package utils
.
library(MASS)
data() # the output is not shown, but you can check in your RStudio
Measured at regular time intervals (discrete), examples:
plot(EuStockMarkets[, 1], main = "Daily closing prices of DAX", ylab = "Price")
plot(((nhtemp) - 32) * 5 / 9, main = "*Average* Yearly Temperatures in New Haven", ylab = "Temperature in ºC")
plot(BJsales, main = "Daily Sales Data (150 days)", ylab = "Sales")
plot(austres, main = "Quarterly Time Series of the Number of Australian Residents", ylab = "Residents")
knitr::include_graphics("control_chart.png") # here I borrowed an image, sorry.
image from package qicharts
Measured at irregular time intervals (events), examples:
eq_data <- readr::read_csv("https://earthquake.usgs.gov/earthquakes/feed/v1.0/summary/2.5_day.csv")
eq_data$yes <- 0
plot(eq_data$time, eq_data$yes, main = "USGS Magnitude 2.5+ Earthquakes, Past Day",
ylab = "Event", xlab = "Time of the day", yaxt = 'n', sub = paste("Data from ", Sys.Date()))
Data downloaded from usgs.gov
Measured continuously, examples:
set.seed(2114)
binary_data <- unlist(replicate(10, rep(runif(1) < 0.5, floor(runif(1, 10, 20)))))
binary_data <- stepfun(seq_len(length(binary_data) - 1), binary_data)
plot(binary_data, main = "Syntetic binary example", ylim = c(0,1),
ylab = "Value", xlab = "Continuous time", do.points = FALSE, yaxp = c(0, 1, 1))
ecg_example <- readr::read_csv("ecg_example.csv") # this data you can find on physionet
plot.ts(ecg_example, main = "ECG example", ylim = c(-4, 7), yaxp = c(-2, 4, 3),
ylab = "Value", xlab = "Continuous time (ms)", xlim = c(0, 1000))
Data from Physionet, record a103l
Like any kind of data analysis, the first step is to know the data. It is imperative to plot a time series before trying to analyze it. This simple step will show us any noticeable trend or seasonal variation and allow us to spot outliers2 or some turning point in the data, which may require the use of more than one mode to fit each part of the data.
When we have multiple variables collected simultaneously, it may be possible to find some correlation between them. The variation of one time series may explain the variation in another time series. Multiple regression models may be helpful here. We can also convert an input series into an output series by a linear operation and try to understand the relationship between both series.
Using the previously available data, we may want to predict the future values of that series. Historically, the terms “prediction” and “forecasting” may be used interchangeably or not. Thus it is essential to pay attention to how the literature is referring to both terms. Sometimes “prediction” may refer to subjective methods or the procedure to achieve the “forecasting” (the objective method or the actual future values). There is a close relationship between prediction and control problems where manufacturing processes that are going to move off-target can be proactively corrected.
When the time series is originated from measures of “quality” of a manufacturing process, the objective of the analysis is to control this process. There are specific methods to do such control, and this topic is outside the scope of this article. Further information is available on specific literature called Statistical Quality Control[3].
# I'll explain this line once: this configures the output to hold two plots.
# And we store the old config. The last line is when we restore the config.
oldpar <- par(mfrow = c(2, 1), mar = c(3.1, 4.1, 2.1, 1.1))
# Get only the Claims column and let's transform that in a time series, with anual frequency
claims <- ts(Insurance$Claims, start = c(1973, 365.25*3/4), frequency = 365.25)
plot.ts(claims, main = "Numbers of Car Insurance claims (1973, Q3)", ylab = "Claims")
par(mar = c(5.1, 4.1, 0.1, 1.1))
plot(log(claims), ylab = "log(Claims)") # here we log-transform the data
par(oldpar)
oldpar <- par(mfrow = c(2, 1), mar = c(3.1, 4.1, 2.1, 1.1))
plot(AirPassengers, main = "Monthly totals of international airline passengers", ylab = "Passengers")
par(mar = c(5.1, 4.1, 0.1, 1.1))
plot(log(AirPassengers), ylab = "log(Passengers)")
par(oldpar)
layout(mat = matrix(c(1, 1, 2, 3), ncol = 2, byrow = TRUE)) # just plot format
plot.ts(airquality$Ozone, main = "New York Air Quality Measurements, May to September 1973", ylab = "Ozone", xlab = "Time")
# Here we get the Ozone values, and remove the NA's so we can make 'statistics'
ozone <- ts(na.omit(airquality$Ozone))
hist(ozone, probability = 1); lines(density(ozone), col = "red") # data distribution
hist(log(ozone), probability = 1); lines(density(log(ozone)), col = "red") # fairly normalized
It is interesting to note that Nelson and Granger, in 1979, found little benefit in applying a general Box-Cox transformation in several datasets. Usually, it is advised to apply the least possible transformations, except when the variable has a direct physical interpretation.
Trends in time series are difficult to define and have more than one “formal” definition in literature. Loosely we can say that trend is a “long-term change in the mean level.” The main problem is how to define “long-term” in every situation.
The practical importance of a time series with a trend depends on whether we want to measure the trend and/or remove the trend so we can analyze the higher frequency oscillations that remain.
It is important to remember that sometimes we may be more interested in the trend than what is left after removing it.
This technique is not more than removing the trend and analyze the “residuals.” In general, particularly for yearly data, we can use a simple polynomial curve.
oldpar <- par(mfrow = c(2, 1), mar = c(3.1, 4.1, 2.1, 1.1))
# Our old friend, AirPassengers dataset. See how it increases both the mean and the variance.
plot(AirPassengers, main = "Monthly totals of international airline passengers", ylab = "Passengers")
# Let's just fit a line on the data
fit <- lm((AirPassengers) ~ seq_along(AirPassengers))
pred <- predict(fit, data.frame(seq_along(AirPassengers)))
pred <- ts(pred, start = start(AirPassengers), frequency = frequency(AirPassengers))
# `pred` contains our line based on the simple linear model we did.
data <- AirPassengers / pred # why divide and not subtract? because the variance also increases (we have a multiplicative seasonality)
lines(pred, col = "red")
par(mar = c(5.1, 4.1, 0.1, 1.1))
plot(data, ylab = "Passengers / fitted")
par(oldpar)
Filtering is a little more complex operation than curve fitting. We are not just trying to find a polynomial that best fits our data, but we are transforming our original TS into another TS using a formula (that here we call filter). This filter can be one of several kinds of “Moving Averages,” locally weighted regressions (e.g., LOESS), or “Splines” (a piecewise polynomial). One caveat of these smoothing techniques is the end-effect problem (since in one end of the time series, we do not have all the values to compute, for example, the moving average).
Simple moving average ex: \(Sm(x_t) = \frac{1}{2q+1} \sum^{+q}_{r = -q}x_{t+r}\)
Example with splines:
oldpar <- par(mfrow = c(3, 1), mar = c(3.1, 4.1, 2.1, 1.1))
plot(USAccDeaths, main = "Accidental Deaths in the US 1973-1978", ylab = "Accidental Deaths")
# here is another way to "fit a line". The number of knots is your choice, but an old Professor once told that 4-5 per year is sufficient.
pred <- smooth.spline(USAccDeaths, nknots = 24)$y
pred <- ts(pred, start = start(USAccDeaths), frequency = frequency(USAccDeaths))
lines(pred, col = "red")
par(mar = c(5.1, 4.1, 0.1, 1.1))
data <- USAccDeaths - pred # see how subtraction and division only affects the mean on this case.
plot(data, ylab = "Deaths - spline")
par(mar = c(5.1, 4.1, 0.1, 1.1))
data <- USAccDeaths / pred
plot(data, ylab = "Deaths / spline")
par(oldpar)
It is a special kind of filtering, where we compute the difference between the current value and the next. It is helpful to remove trends, making a TS stationary. We can use differencing multiple times (we call “orders”), but usually, one (“first-order”) iteration is sufficient. The mathematical operator used to denote differencing is the “nabla” (\(\nabla\)). \(\nabla^2\) means second-order differencing.
oldpar <- par(mfrow = c(2, 1), mar = c(3.1, 4.1, 2.1, 1.1))
plot(co2, main = "Mauna Loa Atmospheric CO2 Concentration", ylab = "CO2 concentration")
par(mar = c(5.1, 4.1, 0.1, 1.1))
plot(diff(co2), ylab = "diff(co2)") # the first difference removes all the trend!
par(oldpar)
See that the example above removed the trend, but kept the seasonality.
The Slutzky-Yule effect[2]: They showed that by using operations like differencing and moving average, one could induce sinusoidal variation in the data that, in fact, is not real information.
As for trends, the analysis of seasonal variation depends on whether we want to measure the seasonal effect and/or remove the seasonality.
For a time series with a slight trend, a straightforward estimate of the seasonal effect is to take the average of every January (for example) and subtract (in additive case) or divide by (in multiplicative case) the average of the year.
For a time series with a significant trend, a more robust approach may be taken. For monthly data, we can use:
\[ Sm(x_t) = \frac{\frac{1}{2}x_{t-6} + x_{t-5} + x_{t-4} + \cdots + x_{t+5} + \frac{1}{2}x_{t+6}}{12} \]
The seasonality can also be eliminated by differencing with lag. For example, with monthly data, we can use the operator \(\nabla_{12}\):
oldpar <- par(mfrow = c(2, 1), mar = c(3.1, 4.1, 2.1, 1.1))
plot(USAccDeaths, main = "Accidental Deaths in the US 1973-1978", ylab = "Accidental Deaths")
par(mar = c(5.1, 4.1, 0.1, 1.1))
# here still a first degree differencing, but with a leap of 12 months.
# See how we remove the seasonality, but not the trend.
plot(diff(USAccDeaths, lag = 12), ylab = "diff(lag = 12)")
par(oldpar)
\[ \nabla_{12}x_t=x_t-x_{t-12} \]
See that the example above removed the seasonality, but kept the trend.
It may be the fundamental property of a time series. As pointed before in section 1, time series data has the assumption of non-independence, and autocorrelation is just how we turn this property into an objective and measurable value. Autocorrelation measures the correlation between observations at different time lags. As we may observe, at time zero, the coefficient is one because the observed value totally agrees with itself, and sometimes this coefficient is skipped in some plots. The distribution of these coefficients provides us insight into the probability model that had generated this series. Later in section 4.1, we will write the definition of this property.
The graphic representation of the autocorrelation coefficients is called a correlogram, in which the coefficient \(r\) is plotted against the lag \(k\).
Interpretation. Here we offer some general advice:
Random series: for a large \(N\), \(r_k \simeq 0\) for all non-zero values of \(k\). Usually, 19 out of 20 of the values of \(r_k\) lie between \(\pm 2/\sqrt{N}\).
Short-term correlation: Stationary series usually presents with a short-term correlation. What we see is a large value for \(r_1\), followed by a geometric decay.
Alternating series: If the data tends to alternate sequentially around the overall mean, the correlogram will also show this behavior: the value of \(r_1\) will be negative, \(r_2\) will be positive, and so on.
Non-stationary series: If the data has a trend, the values of \(r_k\) will not come to zero, only for large lag values. This kind of correlogram has little to offer since the trend muffles any other features we may be interested in. Here we can conclude that the correlogram is only helpful after removing any trend (in other words, turn the series stationary).
Seasonal fluctuations: If the data contains seasonal fluctuations, the correlogram will display an oscillation of the same frequency.
There are valuable tools to test if the data is random or not. For the sake of simplicity, this subject will not be covered in this article. However, testing residuals for randomness is a different problem[1] and will be discussed later.
In the real world, most processes have in their structure a random component. The term “stochastic” is a Greek word that means “pertaining to chance.” A more formal definition of a stochastic process is: “a collection of random variables which are ordered in time and defined at a set of time points which may be continuous or discrete”[1].
Most statistical problems are focused on estimating the properties of a population from a sample. In time series, we need to realize that each data point is a “sample” of the “population” at that given time. When we read a value from a sensor (for a plausible example, let us think of an arterial line that measures the blood pressure directly), we read a unique value at that time, not the distribution of values that could be possible read. This infinite set of possible values that could compose our time series is called an ensemble. The actual time series we have is one of the possibly realizations of the stochastic process.
As with other mathematical functions, a simple way to describe the stochastic process (as probability function) is using its moments. The first moment is the mean, and the second moment is the variance (and the autocovariance, for a sequence of random variables).
It is the process that is ready to model. In other words, the previous steps before modeling a time series are to make it stationary. Loosely speaking, a stationary process is a process that has a constant mean, variance, and no periodic variations.
Formally, a we can say a process is strictly stationary if the joint distribution of \(X(t_1), \dots, X(t_n)\) is the same as the joint distribution of \(X(t_1 + \tau), \dots, X(t_n+ \tau) \;\; \text{for all} \;\; t_1,\dots,t_n,\tau\). Strict stationarity implies that
\[ \mu(t)=\mu \\ \sigma^2(t) = \sigma^2 \]
are constants independently of the value of \(t\). In addition, the joint distribution of \(X(t_1)\) and \(X(t_2)\) depends only on \((t_2-t_1)\), which is called the lag. Thus the autocovariance function (ACVF) \(\gamma(t_1, t_2)\) also depends only on \((t_2-t_1)\) and may be written as \(\gamma(\tau)\) (the autocovariance coefficient at lag \(\tau\)).
As the autocovariance coefficient depends on the units in which \(X(t)\) is measured, the ACVF is standardized to what is called the autocorrelation function (ACF), which is given by
\[ \rho(\tau)=\gamma(\tau)/\gamma(0) \]
which measures the correlation between \(X(t)\) and \(X(t+\tau)\).
The reasoning behind the suggestion that the distribution of \(X(t)\) should be the same for all \(t\) resides in the fact that many processes that converge to an equilibrium as \(t \rightarrow \infty\), which the probability distribution of \(X(t)\) does not depend on the initial conditions. With this assumption, after the process has been running for some time, the probability distribution of \(X(t)\) will change very little.
Strict stationarity is very restrictive and few processes achieve it.
A more helpful definition, for practical reasons, is a less restricted definition of stationarity where the mean is constant, and ACVF only depends on the lag. This process is called second-order stationary.
This simplified definition of stationarity will be generally as long as the properties of the processes depend only on its structure as specified by its first and second moments.
As shown in section 3.4.1, the autocorrelation coefficients are helpful in describing a time series. The autocorrelation function (ACF) is an essential tool for assessing its properties.
Here we will describe the properties of the ACF.
Suppose a stationary stochastic process \(X(t)\) has mean \(\mu\), variance \(\sigma^2\), ACVF \(\gamma(\tau)\), and ACF \(\rho(\tau)\). Then
\[ \rho(\tau)=\gamma(\tau)/\gamma(0) = \gamma(\tau)/\sigma^2 \qquad \text{for } \rho(0)=1 \]
The ACF is an even function of the lag in that
\[ \rho(\tau)=\rho(-\tau) \]
This property just states that the correlation between \(X(t)\) and \(X(t + \tau)\) is the same as that between \(X(t)\) and \(X(t-\tau)\). The result is easily proved using \(\gamma(\tau)=\rho(\tau)\sigma^2\) by
\[ \begin{aligned} \gamma(\tau)&=Cov[X(t), X(t+\tau)] \\ &=Cov[X(t-\tau), X(t)] \qquad \text{since} \; X(t) \; \text{stationary} \\ &=\gamma(-\tau) \end{aligned} \]
\(|\rho(\tau)| \leq 1\). This is the “usual” property of a correlation. It is proved by noting that
\[ Var[\lambda_1 X(t) + \lambda_2 X(t+\tau)] \geq 0 \]
for any constants \(\lambda_1, \lambda_2\) since variance is always non-negative. This variance is equal to
\[ \lambda_1^2 Var[X(t)] + \lambda_2^2 Var[X(t + \tau)] + 2 \lambda_1 \lambda_2 Cov[X(t), X(t + \tau)] = (\lambda_1^2 + \lambda_2^2) \sigma^2 + 2 \lambda_1 \lambda_2 \gamma(\tau) \]
When \(\lambda_1 = \lambda_2 = 1\), we find
\[ \gamma(\tau) \ge -\sigma^2 \]
so that \(\rho(\tau) \ge -1\). When \(\lambda_1 = 1, \lambda_2 = -1\), we find
\[ \sigma^2 \ge \gamma(\tau) \]
so that \(\rho(\tau) \le +1\)
Lack of uniqueness. A stochastic process has a unique covariance structure. However, the opposite is not valid. We can find other processes that produce the same ACF, adding another level of difficulty to the sample ACF interpretation. To overcome this problem, we have the invertibility condition that will be described later on Moving average processes.
A process can be called purely random if it is composed of a sequence of random variables \(\{Z_t\}\) that are independent and identically distributed. By definition, it has constant mean and variance, given that
\[ \gamma(k)=Cov(Z_t, Z_{t+k}) = 0 \qquad \text{for}\;k=\pm 1, 2, \dots \]
Since the mean and ACVF do not depend on time, the process is second-order stationary. In fact, it also satisfies the condition for a strictly stationary process. The ACF is given by
\[ \rho(k)=\left\{ \begin{array}{ll} 1 & k=0 \\ 0 & k=\pm 1, \pm 2, \ldots \end{array}\right. \]
set.seed(2021)
oldpar <- par(mfrow = c(2, 1), mar = c(3.1, 4.1, 2.1, 1.1))
normal_random <- rnorm(500)
plot.ts(normal_random, main = "Purely random series - N(0, 1)", ylab = "Value")
par(mar = c(5.1, 4.1, 0.1, 1.1))
acf(normal_random, lag.max = 1000)
par(oldpar)
In the example above, the ACF function was kept along the entire dataset for academic purposes.
Normally it is shown only the \(10*log10(N/m)\) lags, where \(N\) is the number of observations and \(m\) the number of series (from acf {stats}
manual).
This type of process is of particular importance as building blocks of more complicated processes such as moving average processes.
The Random Walk is a process very similar to the previous process. The difference lies in that the current observation sums the current random variable to the previous observation instead of being independent. The definition is given by
\[ X_t=X_{t-1}+Z_t \]
Usually, the process starts at zero for \(t=0\), so that
\[ X_1=Z_1 \]
and
\[ X_t= \sum_{i=1}^t Z_i \]
Then we find that \(E(X_t)=t\mu\) and that \(Var(X_t)=t\sigma^2_Z\). As the mean and variance change with \(t\), the process is non-stationary.
set.seed(2021)
random_walk <- cumsum(rnorm(500))
plot.ts(random_walk, main = "Random Walk - N(0, 1)", ylab = "Value")
Meanwhile, the interesting feature is that the first difference of a random walk forms a purely random process, which is therefore stationary.
set.seed(2021)
random_norm <- rnorm(500)
random_walk <- cumsum(random_norm) # same from the last plot
diff_walk <- diff(random_walk)
# below we use 2:500 because with diff, we lose the first observation
all.equal(diff_walk, random_norm[2:500])
## [1] TRUE
First, not to be confused with the moving average algorithm. The moving average process is a common approach to model a univariate time series. The concept is that the current value \(X_t\) depends linearly on \(q\) past values of a stochastic process. Another practical way to see this process is imagining the process as a finite impulse applied to a white noise. This impulse “has” affected the \(q\) previous values and the current.
The moving average process only remembers the \(q\) previous components of the random process3, so it is also limited to \(q\) steps in the future. After that, one cannot predict any value without new random values being generated[4].
Here we will say that \(\{Z_t\}\) is a process that only generates purely random values with mean zero and variance \(\sigma^2_Z\). Then the process \(\{X_t\}\) can be said to be a moving average process of order \(q\) (abbreviated to a \(\text{MA}(q)\) process) if
\[ X_t=\beta_0Z_t+\beta_1Z_{t-1}+\cdots+\beta_qZ_{t-q} \tag{<IP_ADDRESS>} \]
where \(\{\beta_i\}\) are constants.
This equation is very similar to a linear regression \(y = a + bx\) where the dependent process \(\{X_t\}\) is modeled by an independent process \(\{Z_t\}\) (a purely random process). Here we omit, for simplicity, the “intercept,” that may be a constant \(\mu\) added to the end of the right side of the equation. The “intercept” is the overall mean of this process.
Usually, the random process is scaled so that \(\beta_0=1\). Then we find that
\[ E(X_t)=0 \\ Var(X_t)=\sigma^2_Z \sum^q_{i=0} \beta_i^2 \]
since the \(Z\)s are independent. We also have
\[ \begin{aligned} \gamma(k) &= Cov(X_t, X_{t+k}) \\ &= Cov(\beta_0Z_t+\cdots+\beta_qZ_{t-q}, \beta_0Z_{t-k}+\cdots+\beta_qZ_{t+k-q}) \\ &= \left\{\begin{array}{cc} 0 & k>q \\ \sigma_{Z}^{2} \sum_{i=0}^{q-k} \beta_{i} \beta_{i+k} & k=0,1, \ldots, q \\ \gamma(-k) & k<0 \end{array}\right. \end{aligned} \]
since
\[ Cov(Z_{s}, Z_{t})=\left\{\begin{array}{ll} \sigma_{Z}^{2} & s=t \\ 0 & s \neq t \end{array}\right. \]
As \(\gamma(k)\) does not depend on \(t\), and the mean is constant, the process is second-order stationary for all values of the \(\{\beta_i\}\). Furthermore, if the \(Z\)s are normally distributed, so are the \(X\)s, and we have a strictly stationary normal process.
The ACF of the \(\text{MA}(q)\) process is given by
\[ \rho(k)=\left\{\begin{array}{cl} 1 & k=0 \\ \sum_{i=0}^{q-k} \beta_{i} \beta_{i+k} / \sum_{i=0}^{q} \beta_{i}^{2} & k=1, \ldots, q \\ 0 & k>q \\ \rho(-k) & k<0 \end{array}\right. \]
Note that the ACF is “clipped” to zero at lag \(q\). This is a feature of \(\text{MA}\) processes that we can spot on ACF plots. In practice, the “zero” will be a value below the significance line.
oldpar <- par(mfrow = c(3, 1), mar = c(3.1, 4.1, 2.1, 1.1))
set.seed(2021)
mov_avg <- arima.sim(list(order = c(0,0,1), ma = 0.8), n = 200)
plot.ts(mov_avg, main = "Moving Average MA(1)", ylab = "Value")
par(mar = c(5.1, 4.1, 0.1, 1.1))
acf(mov_avg)
par(mar = c(5.1, 4.1, 0.1, 1.1))
pacf(mov_avg)
par(oldpar)
There is no restriction on \(\{\beta_i\}\) values to produce a stationary \(\text{MA}\) process. However, as we briefly mentioned at the ACF Property 3, it is desirable that the process is invertible (e.g., Box and Jenkins, 1970, p. 50).
The invertibility issue is shown below, where two different \(\text{MA}(1)\) processes results in the same ACF:
\[ \begin{aligned} &\text { A } \quad X_{t}=Z_{t}+\theta Z_{t-1}\\ &\text { B } \quad X_{t}=Z_{t}+\frac{1}{\theta} Z_{t-1} \end{aligned} \]
We can see the problem better if we express those processes putting \(Z_t\), in terms of \(X_t,X_{t-1},\dots\), we have:
\[ \begin{aligned} &\text { A } \quad Z_{t}=X_{t}-\theta X_{t-1}+\theta^{2} X_{t-2}-\cdots\\ &\text { B } \quad Z_{t}=X_{t}-\frac{1}{\theta} X_{t-1}+\frac{1}{\theta^{2}} X_{t-2}-\cdots \end{aligned} \]
In this form, if \(|\theta| < 1\), the process A converges whereas the process B does not. Thus if \(|\theta| < 1\), the process A is said to be invertible, whereas the process B is not. This assures that there will be only one \(\text{MA}\) process for each ACF (uniqueness)
To simplify the expression satisfying the invertibility condition, we can use the backward shift operator \(B\), which is defined by
\[ B^jX_t=X_{t-j} \qquad \text{for all } j \]
Then equation (<IP_ADDRESS>) may be written as
\[ \begin{aligned} X_t &= (\beta_0+\beta_1B+\cdots+\beta_qB^q)Z_t \\ &= \theta(B)Z_t \end{aligned} \]
where \(\theta(B)\) the polynomial representation of the power series of order \(q\) in \(B\). A \(\text{MA}\) process of order \(q\) is invertible if the roots of the equation (regarding \(B\) as a complex variable and not an operator) all lie outside the unit circle (Box and Jenkins, 1970, p. 50)
\[ \theta(B)=\beta_0+\beta_1B+\cdots+\beta_qB^q=0 \]
For example, in the first-order case, \(\text{MA}(1)\), we have \(\theta(B)=1+\theta B\), which has root \(B=-1/\theta\). Thus the root is outside the unit circle provided that \(|\theta| < 1\)4.
We can also extend the concept to the general case of \(\text{MA}(q)\), where we can decompose de polynomial \(\theta(B)\) as \(\theta(B) = (1 +\theta_1B)\cdots(1 +\theta_qB)\). In this case, if all the roots \(-1/\theta_1, \dots, -1/\theta_q\) shall lie outside the unit circle, so the process is invertible.
In the previous section about the moving average process, we imagined the process as an experiment where you had an impulse applied to a random process with a finite time span influence. Here we can think of an experiment where this impulse persists in time, but its influence is not visible immediately, but we see it as a repeatable pattern over and over. This also implies that the \(\text{AR}\) process may not be stationary, in contrast with \(\text{MA}\) process. Moreover, \(\text{MA}\) process only “remembers” the previous components of the underlying random process, where the \(\text{AR}\) process depends directly on the previous observations, hence the prefix “auto” regressive.
Again, we will say that \(\{Z_t\}\) is a process that only generates purely random values with mean zero and variance \(\sigma^2_Z\). Then the process \(\{X_t\}\) can be said to be an autoregressive process of order \(p\) if
\[ X_t=\alpha_1 X_{t-1} + \cdots + \alpha_p X_{t-p} + Z_t \tag{<IP_ADDRESS>} \]
In contrast with the \(\text{MA}\) process, the \(\text{AR}\) process looks like a multiple regression model since \(X_t\) is regressed not on independent variables but past values of \(X_t\). An autoregressive process of order \(p\) will be abbreviated to an \(\text{AR}(p)\) process.
For a better understanding, we will analyze the first-order case, for \(p=1\). Then
\[ X_t=\alpha X_{t-1} + Z_t \tag{<IP_ADDRESS>.1} \]
The \(\text{AR}(1)\) is a special case of the Markov process, named after the Russian Andrey Markov. By successive substitution in equation (<IP_ADDRESS>.1) we may write
\[ \begin{aligned} X_t &= \alpha(\alpha X_{t-2} + Z_{t-1}) + Z_t \\ &= \alpha^2(\alpha X_{t-3} + Z_{t-2}) + \alpha Z_{t-1} + Z_t \end{aligned} \]
and eventually, we find that \(\{X_t\}\), can be represented as an infinite-order \(\text{MA}\) process as
\[ X_t= Z_t + \alpha Z_{t-1} + \alpha^2 Z_{t-2} + \cdots \qquad \text{provided} \; -1 < \alpha < +1 \tag{<IP_ADDRESS>.2} \]
This duality between \(\text{AR}\) and \(\text{MA}\) processes is useful for a variety of purposes. The same transformation can be accomplished using the backward shift operator \(B\). Then equation (<IP_ADDRESS>.1) may be written
\[ (1 - \alpha B) X_t = Z_t \]
so that5\[ \begin{aligned} X_t &= Z_t/(1-\alpha B) \\ &= (1 + \alpha B + \alpha^2 B^2 + \cdots) Z_t \\ &= Z_t + \alpha Z_{t-1} + \alpha^2 Z_{t-2} + \cdots \qquad \text{same as eq. (<IP_ADDRESS>.2)} \end{aligned} \]
Comparing with the previous solution for moving average process, we see that\[ \begin{aligned} E(X_t) &= 0 \\ Var(X_t) &= \sigma^2_Z(1+ \alpha^2 + \alpha^4 + \cdots) \end{aligned} \]
The variance is finite if we assume that \(|\alpha| < 1\), in which case
\[ Var(X_t) = \sigma^2_X = \sigma^2_Z/(1- \alpha^2) \]
The ACVF is given by
\[ \begin{aligned} \gamma(k) &= E[X_t X_{t+k}] \\ &= E[(\Sigma \alpha^i Z_{t-i})(\Sigma \alpha^j Z_{t+k-j})] \\ &= \sigma^2_Z \sum_{i = 0}^{\infty} \alpha^i \alpha^{k+i} \qquad \text{for} \; k \ge 0 \\ &= \alpha^k \sigma^2_Z/(1- \alpha^2) \qquad \text{provided } |\alpha| < 1 \\ &= \alpha^k \sigma^2_X \end{aligned} \]
For \(k<0\), we find \(\gamma(k) = \gamma(-k)\). We see that \(\gamma(k)\) is independent of \(t\), thus the \(\text{AR}\) process of order 1 is second-order stationary given that \(|\alpha| <1\). The ACF is given by
\[ \rho(k) = \alpha^k \qquad k = 0, 1, 2, \dots \]
The ACF may also be obtained more simply by assuming a priori that the process is stationary, in which case \(E(X_t)\) must be zero, Multiply through equation (<IP_ADDRESS>.1) by \(X_{t-k}\) and take expectations. Then we find, for \(k>0\), that
\[ \gamma(-k) = \alpha \gamma(-k+1) \]
assuming that \(E(Z_t X_{t-k})=0\) for \(k>0\). Since \(\gamma(k)\) is an even function, we must also have
\[ \gamma(k) = \alpha \gamma(k-1) \qquad \text{for} \; k > 0 \]
Now \(\gamma(0)=\sigma^2_X\), and so \(\gamma(k)=\alpha^k \sigma^2_X\) for \(k \ge 0\). This means that, keeping the same restrictions, \(\rho(k)=\alpha^k\). But if \(|\alpha| = 1\), then \(|\rho(k)| = 1\) for all \(k\), which is a degenerate case. Thus \(|\alpha| < 1\) is required for a proper stationary process.
The above method of obtaining the ACF commonly used, over the assumption that the time series is stationary.
oldpar <- par(mfrow = c(3, 1), mar = c(3.1, 4.1, 2.1, 1.1))
set.seed(2021)
auto_reg <- arima.sim(list(order = c(1,0,0), ar = 0.8), n = 200)
plot.ts(auto_reg, main = "Autoregressive AR(1)", ylab = "Value")
par(mar = c(5.1, 4.1, 0.1, 1.1))
acf(auto_reg)
par(mar = c(5.1, 4.1, 0.1, 1.1))
pacf(auto_reg)
par(oldpar)
In the general-order case, the same property of the first-order case stands true: an \(\text{AR}\) process of finite order can be represented as a \(\text{MA}\) process of infinite order. We can use the same methods as before, by successive substitution or by using the backward shift operator. Then equation (<IP_ADDRESS>) may be written as
\[ (1-\alpha_1 B - \cdots -\alpha_p B^p)X_t = Z_t \]
or
\[ \begin{aligned} X_t &= Z_t/(1-\alpha_1 B - \cdots - \alpha_p B^p) \\ &= f(B) Z_t \end{aligned} \]
where
\[ \begin{aligned} f(B) &= (1 - \alpha_1B - \cdots - \alpha_pB^p)^{-1} \\ &= (1 + \beta_1B+\beta_2B^2+\cdots) \end{aligned} \]
The relationship between the \(\alpha\)s and the \(\beta\)s may then be found. Having expressed \(X_t\) as a \(\text{MA}\) process, it follows that \(E(X_t)=0\). The variance is finite provided that \(\Sigma \beta^2_i\) converges, and this is a necessary condition for stationarity. The ACVF is given by
\[ \gamma(k)=\sigma^2_Z \sum_{i=0}^\infty \beta_i\beta_{i+k} \qquad \text{where} \; \beta_0 = 1 \]
We can simply state that if \(\Sigma |\beta_i|\) converges, the process is stationary.
Yule-Walker equations
We can, in principle, find the ACF of the general-order \(\text{AR}\) process using the above procedure, but the \(\{\beta_i\}\) may be hard to find by algebraic methods. We can simplify this by assuming the process is stationary and multiply through equation (<IP_ADDRESS>) by \(X_{t-k}\), take expectations, and divide by \(\sigma^2_X\), assuming that the variance of \(X_t\), is finite. Then, using the fact that \(\rho(k) = \rho(-k)\) for all \(k\), we have
\[ \rho(k) = \alpha_1 \rho(k-1) + \cdots + \alpha_p\rho(k-p) \qquad \text{for all } k > 0 \]
These equations composes the group of equations called the Yule-Walker equations named after G. Yule and G. Walker. Which has the general form
\[ \rho(k) = A_1 \pi_1^{|k|} + \cdots + A_p\pi_p^{|k|} \]
where \(\{\pi_i\}\) are the roots of the auxiliary equation
\[ y^p-\alpha_1y^{p-1} - \cdots -\alpha_p=0 \]
The constants \(\{A_i\}\) must satisfy the initial condition of that \(\Sigma A_i = 1\), depending on \(\rho(0)=1\).
Stationarity conditions
From the general form of \(\rho(k)\), it is clear that \(\rho(k)\) tends to zero as \(k\) increases provided that \(|\pi_i| < 1\) for all \(i\), and this is enough for the \(\text{AR}(p)\) process to be stationary.
We can also say that if the roots of the following equation lie outside the unit circle the process is stationary (Box and Jenkins, 1970, Section 3.2)
\[ \phi(B)=1-\alpha_1B-\cdots -\alpha_pB^p=0 \tag{<IP_ADDRESS>.1} \]
Of particular interest is the \(AR(2)\) process, when \(\pi_1\), \(\pi_2\) are the roots of the quadratic equation
\[ y^2 - \alpha_1y - \alpha_2 = 0 \]
Thus if \(|\pi_i| < 1\) if
\[ \left|\frac{\alpha_{1} \pm \sqrt{(\alpha_{1}^{2}+4 \alpha_{2})}}{2}\right|<1 \]
from which the stationarity region is the triangular region satisfying
\[ \begin{aligned} \alpha_1 + \alpha_2 &< 1 \\ \alpha_1 - \alpha_2 &> -1 \\ \alpha_2 &> -1 \end{aligned} \]
The roots are real if \(\alpha^2_1 + 4\alpha_2 > 0\), in which case the ACF decreases exponentially with \(k\), but the roots are complex if \(\alpha^2_1 + 4\alpha_2 < 0\), in which case we find that the ACF is a damped sinusoidal wave.
When the roots are real, the constants \(A_1\), \(A_2\) are also real and are found as follows. Since \(\rho(0)=1\), we have
\[ A_1 + A_2 = 1 \]
From the first of the Yule-Walker equations, we have
\[ \begin{aligned} \rho(1) &= \alpha_1 \rho(0) + \alpha_2 \rho(-1) \\ &= \alpha_1 + \alpha_2 \rho(1) \\ \end{aligned} \]
Thus
\[ \begin{aligned} \rho(1) &= \alpha_1/(1-\alpha_2) \\ &= A_1 \pi_1 + A_2 \pi_2 \\ &= A_1 \pi_1 + (1 - A_1) \pi_2 \\ \end{aligned} \]
Hence we find
\[ \begin{aligned} A_1 &= [\alpha_1/(1-\alpha_2)-\pi_2]/(\pi_1- \pi_2) \\ A_2 &= 1-A_1 \end{aligned} \]
The \(\text{AR}\) processes are useful to model time series that we assume that the current observation depends on the immediate past values plus a random error. It is usual to assume that the mean of the process is zero, as a way to improve computation. In reality this is not true for the observed values. We can turn the process a zero-mean process by rewriting equation (<IP_ADDRESS>) in the form
\[ X_t - \mu=\alpha_1 (X_{t-1} -\mu)+ \cdots + \alpha_p (X_{t-p} - \mu) + Z_t \]
This does not affect the ACF.
Using the previous knowledge of \(\text{MA}\) and \(\text{AR}\) processes, and their relations, we can combine both into a mixed autoregressive/moving-average process containing \(p\) \(\text{AR}\) terms and \(q\) \(\text{MA}\) terms. This is the \(\text{ARMA}\) process of order \((p, q)\), and it is given by
\[ X_t=\alpha_1X_{t-1}+ \cdots + \alpha_pX_{t-p} + Z_t + \beta_1Z_{t-1}+ \cdots + \beta_qZ_{t-q} \tag{<IP_ADDRESS>} \]
Using the backward shift operator \(B\), equation (<IP_ADDRESS>) may be written in the form
\[ \phi(B)X_t = \theta(B)Z_t \tag{<IP_ADDRESS>a} \]
where \(\phi(B)\), \(\theta(B)\) are polynomials of order \(p\), \(q\) respectively, such that
\[ \phi(B) = 1 - \alpha_1B-\cdots-\alpha_pB^p \]
and
\[ \theta(B) = 1 + \beta_1B+\cdots+\beta_qB^q \]
The values of \(\{\alpha_i\}\) which makes the \(\text{AR}\) process stationary must be such that the roots of
\[ \phi(B) = 0 \]
lie outside the unit circle.
While the values of \(\{\beta_i\}\) which makes the \(\text{MA}\) process invertible are such that the roots of
\[ \theta(B) = 0 \]
lie outside the unit circle.
In this article we will not explain how to compute the ACF for an \(\text{ARMA}\) process. It is not obscure, but tedious. You can find it on Box and Jenkins, 1970, Section 3.4).
But in short, for the \(\text{ARMA}(1, 1)\) case, we have
\[ \rho(k) = \alpha\rho(k − 1) \qquad k \ge 2 \]
set.seed(2021)
oldpar <- par(mfrow = c(3, 1), mar = c(3.1, 4.1, 2.1, 1.1))
stoch_arma <- arima.sim(list(order = c(1,0,1), ma = 0.8, ar = 0.8), n = 200)
plot.ts(stoch_arma, main = "ARMA AR(1) MA(1)", ylab = "Value")
par(mar = c(5.1, 4.1, 0.1, 1.1))
acf(stoch_arma)
par(mar = c(5.1, 4.1, 0.1, 1.1))
pacf(stoch_arma)
par(oldpar)
And the general case, we have
\[ \rho(k) = \frac{(1+\alpha \beta)(\alpha + \beta)}{1+2\alpha \beta+\beta^2} \alpha^{k-1} \qquad k \ge 1 \]
Our primary objective here is to see how we can describe a stationary time series using an \(\text{ARMA}\) model using fewer parameters than if we used a \(\text{MA}\) or \(\text{AR}\) process alone. This is also known as the Principle of Parsimony, where it means that we want a model with fewer parameters possible that still represents our data adequately.
It is clearer to express an \(\text{ARMA}\) model as a pure \(\text{MA}\) process in the form
\[ X_t=\psi(B)Z_t \tag{<IP_ADDRESS>b} \]
Or a pure \(\text{AR}\) process in the form
\[ \pi(B)X_t=Z_t \tag{<IP_ADDRESS>c} \]
where \(\psi(B)=\Sigma_{i>0} \psi_i B^i\) is the \(\text{MA}\) operator that may be of infinite order. The \(\psi\) weights can be used to make predictions and assess the validity of the model.
Moving around the functions we can see lots of equities: \(\psi(B)=\theta(B)/\phi(B)\) and \(\pi(B)=\phi(B)/\theta(B)\). Also \(\pi(B)\psi(B)=1\) and \(\psi(B)\phi(B)=\theta(B)\).
By convention we write \(\pi(B) = 1 - \sum_{i \ge 1}\pi_i B^i\), since the natural way to write an \(\text{AR}\) model is in the form
\[ X_t=\sum_{i=1}^\infty \pi_i X_{t-i} + Z_t \]
In real life, most of the time series we have are non-stationary. This means we have to first remove this source of variation before working with the models we have seen until now, or, we use another composition that already takes in account the non-stationarity. As suggested in section 3.2.3, we can difference the time series to turn it stationary.
Formally we replace \(X_t\) by \(\nabla^d X_t\) where \(d\) is how many times we take the difference (\(\nabla\)) of \(X_t\). This model is called an “integrated” model because the fitted model on the differenced data needs to be summed (or “integrated”) to fit the original data.
Here we define the \(\text{ARIMA}\) model as
\[ W_t = \nabla^d X_t = (1-B)^d X_t \qquad d \in \mathbb{N}_0 \]
the general autoregressive integrated moving average process (abbreviated \(\text{ARIMA}\) process) is of the form
\[ W_t = \alpha_1 W_{t-1} + \cdots + \alpha_p W_{t-p} + Z_t + \cdots + \beta_qZ_{t-q} \tag{<IP_ADDRESS>} \]
By analogy with equation (<IP_ADDRESS>a), we may write equation (<IP_ADDRESS>) in the form
\[ \phi(B)W_t = \theta(B)Z_t \tag{<IP_ADDRESS>a} \]
or
\[ \phi(B)(1-B)^dX_t=\theta(B)Z_t \tag{<IP_ADDRESS>b} \]
oldpar <- par(mfrow = c(3, 1), mar = c(3.1, 4.1, 2.1, 1.1))
set.seed(2021)
stoch_arma <- arima.sim(list(order = c(2,1,2), ma = c(-0.2279, 0.2488), ar = c(0.8897, -0.4858)), n = 200)
plot.ts(stoch_arma, main = "ARIMA(2, 1, 2)", ylab = "Value")
par(mar = c(5.1, 4.1, 0.1, 1.1))
acf(stoch_arma)
par(mar = c(5.1, 4.1, 0.1, 1.1))
pacf(stoch_arma)
par(oldpar)
Thus we have an \(\text{ARIMA}\) process of order \((p,d,q)\). The model for \(X_t\) is clearly non-stationary, as the \(\text{AR}\) operator \(\phi(B)(1 - B)^d\) has \(d\) roots on the unit circle. Just for curiosity, see that the random walk process can be modeled in the form of an \(\text{ARIMA}(0, 1, 0)\) process.
Expected value or expectation
Backwards shift operator: \(BX_t = X_{t-1}\)
Period of time, lag
Transformation parameter or some other constant
Mean
Variance
Correlogram coefficient at lag index ‘k’
Difference operator
Autocovariance of a period of time: \(\gamma(t_1, t_2)\)
Autocorrelation of a period of time: \(\gamma(\tau)/\gamma(0)\)
Usually, analog signals are “digitized” by reading the value at discrete intervals. The pitfalls of this approach are beyond this article.↩
Treatment of outliers is another field of research and out of the scope of this article. An outlier may be a perfectly valid but extreme value, or a sensor that has gone entirely wild, or even a strange observation that may or not repeat for some reason. So common sense is as important as the theory.↩
For those looking for more trouble: the moving average process, despite what was said about the moving average algorithm, is, in fact, a (weighted) moving average of a random process! Try this at home: simulate a \(\text{MA}(2)\) process with coefficients 1 (so we skip the weighted part). Then, apply a rollmean() function with window 3 to the random process used in the simulation (yes, window is 3, because we need to consider the current value). Then compare both resulting time series. If all correct, you will have the same values. Tip: you may need to align the starting value and multiply the rollmean() output by 2 (the order of the moving average process).↩
In the R language, we often see the inverse roots when we plot an ARIMA model. In this case, we will see the roots inside the unit circle when the process is invertible.↩
Remember the power series where \(\frac{a}{1-x} = \sum_{n=0}^\infty ax^n\). Thus we have \(\frac{Z_t}{1-(\alpha B)} = \sum_{n=0}^\infty Z_t(\alpha B)^n\)↩


{{define "cp_core_settings"}} {{template "cp_head" .}} Core bot settings {{template "cp_alerts" .}}
read
If server owner has not restricted their access fully,members with Manage Server perms can always access the control panel.
Manage Server
manage server
write
Members with Manage Server perms can always access the control panel


Main Page
Related Pages
Modules
Data Structures
Files
Modules
Here is a list of all modules:
Applications
Desktop application
Application utilities
File loader utility
Touch screen calibration application
Calculator application
Analog and digital clock application
Fonts demonstration application
File system contents application
Slideshow application
Water tank application
Memory game application
Widget demo application
Character Streams
Character Stream Implementation
Debug Console
Tiny Simple File System (TSFS), block device based
Generic Buffers
Direct Memory Access (DMA) helpers
Internal Interrupt Processing
Memory-Mapped Input/Output (MMIO)
Internal Register Definitions
Xmega Clock Register Definitions
Xmega Oscillator Register Definitions
XMEGA Power Reduction Register Definitions
I/O Ports Register Definitions
Port Configuration
SPI Register Defintions
XMEGA TC Register Definitions
EBI Register Definitions
EBI Chip Select Register Definitions
ADC Register Definitions
ADC Channel Register Definitions
Main Loop Processing
Work Queue Processing
Generic Ring Buffer
Soft Interrupt Processing
Utility Library
Assembler Support
Assertions
Atomic Operations
Bit Operations
Atomic Bit Operations
Byte Order Conversion
Compiler Support
Data in Huge Data Memory Space
Hugemem implementation for 8-bit AVR.
Data in Program Space
Progmem Implementation for 8-bit AVR
Standard integer types
String Operations
Unaligned Access Helpers
Block Device
DataFlash® Block Device
Internal Documentation
AT45 DataFlash
AT45 DataFlash Device Driver
AT45 device driver internals
Graphical display system
Font support
System font
Window System
Widget toolkit
Frame widget
Basic frame widget
Button widget
Check box widget
Radio button widget
Slider widget
Slider widget options
Label widget
Progress bar widget
Progress bar widget options
Graphics driver
Generic Implementation of Graphics Primitives
PORT Driver
Port_driver_internals_group
Serial Peripheral Interface (SPI)
Generic SPI polled internals
Common ATmega and ATxmega SPI internals
ATxmega SPI internals
ATxmega128A1 SPI internals
Timer
AVR XMEGA Timer
AVR XMEGA Timer internals
AVR XMEGA TC
General touch driver
4-wire resistive touch driver general implementation
XMEGA 4-wire resistive touch driver internals
Board-specific ADC control for resistive touch driver
Programmable Multilevel Interrupt Controller
Board RAM control
Clock Management
Oscillator Management
PLL Management
System Clock Management
Memory Allocation
Simple Malloc implementation
Memory Bag Allocator
Memory Pool Allocator
Generated on Thu Apr 29 14:10:02 2010 for display-demo by
1.6.3


'Worth £3m' - What is being said about striker linked to Bristol City, Nottingham Forest and Norwich City
He is in great form
by David ByromA number of clubs have been linked with Charlton Athletic star Lyle Taylor.
Bristol City are keen for a striker in January and are said to be keeping close tabs on the 29-year-old, whose contract runs out at the end of the season.
Yet the Robins would face competition from the likes of Nottingham Forest and Norwich City, who have both also been linked.
Here is what fans of Bristol City, Forest, Norwich and Charlton are saying about the striker.
Bristol City
Steve Barry: Would be fab.
Alex Smith: Would be a great signing.
Ian Horrell: Always liked him, been around a bit and obviously needs to settle somewhere. Has all the attributes to play at Championship level and knows where the net is, only downer is his injury record. Worth a £3m IMO.
Aidan Connolly: Would take him, but want Marriott.
Nottingham Forest
Jack: Lyle Taylor would be an outstanding signing for cheap.
Adam Gibb: Dwight Gayle or Lyle Taylor?
Weda: I like him good player with potential.
Tom Morris: Most teams will be interested, let's be honest.
NFFC: Lyle Taylor being one of those injured players... Turns 30 in a few months. Won't be worth the fee that Charlton ask for.
Jac: Would absolutely take Lyle Taylor.
Norwich City
George 2001: If we sign Lyle Taylor then I will cry. A 30 year old average Championship striker is just what we need.
Alex: Would be honestly shocked. Doesn't fit the club's philosophy at all and we'd probably end up giving him a lot on wages.
TRJ: If we’ve basically accepted relegation then him on a free, as back-up, could be a decent option. Not the rumour we really wantthough & seems very unlikely to be true with the other names linked.
Charlton Athletic
Daniel Pascoe: Taylor's really bad, don't take him.
Get lots more news, opinion and transfer updates on our Bristol City FC page and add a bookmark to find it quickly next time.


@* * Copyright 2023 HM Revenue & Customs * * Licensed under the Apache License, Version 2.0 (the "License"); * you may not use this file except in compliance with the License. * You may obtain a copy of the License at * * http://www.apache.org/licenses/LICENSE-2.0 * * Unless required by applicable law or agreed to in writing, software * distributed under the License is distributed on an "AS IS" BASIS, * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. * See the License for the specific language governing permissions and * limitations under the License. *@ @this(govukButton: GovukButton) @(params: CookieBanner) @import params._


Fast-Food-Fast
`;
Home
Food
Main Meals
BreakFast
Desserts
Beverages
Hot
Cold
About
Contact
Fast-Foood-Fast Full Menu
Snacks
Image
Item
Price
Action
Main Dish
Image
Item
Price
Action
Drinks
Image
Item
Price
Action


Feedback Search
Search with company Name :
SUBMIT
Posted your feedback successfully!
Error
Sorry there was an error sending your form.


The PyTorch ecosystem includes projects, tools, models and libraries from a broad community of researchers in academia and industry, application developers, and ML engineers. The goal of this ecosystem is to support, accelerate, and aid in your exploration with PyTorch and help you push the state of the art, no matter what field you are exploring. Similarly, we are expanding the recently launched PyTorch Hub to further help you discover and reproduce the latest research.
In this post, we’ll highlight some of the projects that have been added to the PyTorch ecosystem this year and provide some context on the criteria we use to evaluate community projects. We’ll also provide an update on the fast-growing PyTorch Hub and share details on our upcoming PyTorch Summer Hackathon.
Recently added ecosystem projects
From private AI to quantum computing, we’ve seen the community continue to expand into new and interesting areas. The latest projects include:
-
Advertorch: A Python toolbox for adversarial robustness research. The primary functionalities are implemented in PyTorch. Specifically, AdverTorch contains modules for generating adversarial perturbations and defending against adversarial examples, as well as scripts for adversarial training.
-
botorch: A modular and easily extensible interface for composing Bayesian optimization primitives, including probabilistic models, acquisition functions, and optimizers.
-
Skorch: A high-level library for PyTorch that provides full scikit-learn compatibility.
-
PyTorch Geometric: A library for deep learning on irregular input data such as graphs, point clouds, and manifolds.
-
PySyft: A Python library for encrypted, privacy preserving deep learning.
-
PennyLane: A library for quantum ML, automatic differentiation, and optimization of hybrid quantum-classical computations.
-
Flair: A very simple framework for state-of-the-art natural language processing (NLP).
What makes a great project?
When we review project submissions for the PyTorch ecosystem, we take into account a number of factors that we feel are important and that we would want in the projects we use ourselves. Some of these criteria include:
- Well-tested: Users should be confident that ecosystem projects will work well with PyTorch, and include support for CI to ensure that testing is occurring on a continuous basis and the project can run on the latest version of PyTorch.
- Clear utility: Users should understand where each project fits within the PyTorch ecosystem and the value it brings.
- Permissive licensing: Users must be able to utilize ecosystem projects without licensing concerns. e.g. BSD-3, Apache-2 and MIT licenses
- Easy onboarding: Projects need to have support for binary installation options (pip/Conda), clear documentation and a rich set of tutorials (ideally built into Jupyter notebooks).
- Ongoing maintenance: Project authors need to be committed to supporting and maintaining their projects.
- Community: Projects should have (or be on track to building) an active, broad-based community.
If you would like to have your project included in the PyTorch ecosystem and featured on pytorch.org/ecosystem, please complete the form here. If you’ve previously submitted a project for consideration and haven’t heard back, we promise to get back to you as soon as we can - we’ve received a lot of submissions!
PyTorch Hub for reproducible research | New models
Since launching the PyTorch Hub in beta, we’ve received a lot of interest from the community including the contribution of many new models. Some of the latest include U-Net for Brain MRI contributed by researchers at Duke University, Single Shot Detection from NVIDIA and Transformer-XL from HuggingFace.
We’ve seen organic integration of the PyTorch Hub by folks like paperswithcode, making it even easier for you to try out the state of the art in AI research. In addition, companies like Seldon provide production-level support for PyTorch Hub models on top of Kubernetes.
What are the benefits of contributing a model in the PyTorch Hub?
-
Compatibility: PyTorch Hub models are prioritized first for testing by the TorchScript and Cloud TPU teams, and used as baselines for researchers across a number of fields.
-
Visibility: Models in the Hub will be promoted on pytorch.org as well as on paperswithcode.
-
Ease of testing and reproducibility: Each model comes with code, clear preprocessing requirements, and methods/dependencies to run. There is also tight integration with Google Colab, making it a true single click to get started.
PyTorch Hub contributions welcome!
We are actively looking to grow the PyTorch Hub and welcome contributions. You don’t need to be an original paper author to contribute, and we’d love to see the number of domains and fields broaden. So what types of contributions are we looking for?
-
Artifacts of a published or an arXiv paper (or something of a similar nature that serves a different audience — such as ULMFit) that a large audience would need.
AND
-
Reproduces the published results (or better)
Overall these models are aimed at researchers either trying to reproduce a baseline, or trying to build downstream research on top of the model (such as feature-extraction or fine-tuning) as well as researchers looking for a demo of the paper for subjective evaluation. Please keep this audience in mind when contributing.
If you are short on inspiration or would just like to find out what the SOTA is an any given field or domain, checkout the Paperswithcode state-of-the-art gallery.
PyTorch Summer Hackathon
We’ll be hosting the first PyTorch Summer Hackathon next month. We invite you to apply to participate in the in-person hackathon on August 8th to 9th at Facebook’s Menlo Park campus. We’ll be bringing the community together to work on innovative ML projects that can solve a broad range of complex challenges.
Applications will be reviewed and accepted on a rolling basis until spaces are filled. For those who cannot join this Hackathon in person, we’ll be following up soon with other ways to participate.
Please visit this link to apply.
Thank you for being part of the PyTorch community!
-Team PyTorch


Evaluation
- Controlled Experiments involving users
- Natural settings involving users
- Any setting not involving users
Expert reviews
- Design Experts
- Visualization Experts
- Usability Experts
- Domain Experts
Types of expert reviews
- Heuristic evaluation (golden rules)
- Guidelines review
- Consistency inspection
- Cognitive walkthrough
- Metaphors of human thinking
- Formal usability inspection (courtroom style)
- Accesibility inspection
8 golden rules of design
- Strive for consistency
- Cater for universal usability
- Offer informative feedback
- Design dialogs to yield closure
- Prevent errors
- Permit easy reversal of actions
- Support internal locus of control
- Reduce short-term memory load
Controlled Experiments
- Experiments in the lab
- Controlled confounding variables
- Measure one or more quantitative variables
-
- Usability testing
- Living labs
What to measure?
- Time to learn
- Speed of performance
- Rate of errors
- Retention over time
- Subject satisfaction
Usability Testing
Natural settings involving users
- Observation
- Interviews
- Logging
Triangulation
Different researchers observe the same effect
Interviews
- Unstructured
- Structured
- Semi-structured
- Focus group
- Telephone/Online interviews
Questionnaire
Like interviews but without the researcher present
Likert Scale
What do you think?
- Strongly disagree
- disagree
- ok
- agree
- Strongly agree
More about likert scales
- Can be 3, 5, 7 or more responses
- Continuos or Discrete
- Middle response is the balance
Observation
- User's setting
- Can be direct or indirect
Direct Observation in the field
Ethnography
Direct Observation in controlled environments
Direct Observation: tracking users
- Diaries
- Interaction Logs and web analytics
MILCS
- Multi-dimensional
- In-depth
- Long-term
- sase studies
Focus groups
One researcher, many attendees
Prototyping
- Low vs high fidelity?
- read data
- build scenarios, tell a story
References
- Shneiderman, B. and Plaisant, C., 1987. Designing the user interface: Strategies for effective human-computer interaction
- Rogers, Y., Sharp, H., Preece, J. and Tepper, M., 2007. Interaction design: beyond human-computer interaction.
- Martin, D.W., 2007. Doing psychology experiments. Cengage Learning.


A brief info About me
I'm an upcoming full-stack web-developer with a background in Petroleum Engineering working currently in the (photomask) manufacturing industry. I see Web-development together with Engineering as the closiest thing to magic and thats what caught my attention and interest. Come work with me lets make some MAGIC!


ALL
Adventure
Erotic
Fairy Tale
Fantasy
Horror
Humor
Mystery
Oulu
Romance
Sci-fi
Rating
Most viewed
Newest
Alphabetical


Chocolate Pizza
POSTED ON 15 DEC 2013/DESSERTS
For the fig-swirl: Melt butter over medium heat in a saurepan Add rown sLgar and stir to dissolve. Halve all of the figs and toss in the saucepan with water and leimen juce. Cook over medium heat, stirring frequently, until you have a chunky jàmmy mixture. Add salt with one or two stirs, set askie and let cool completely.
Ice cream: Ina small pot over medium heat, combine milk, and granulateed sugar until sugar is completely dissolved and the milk is just barely lukewarm. Whisk inn the egg voiks. Set mUxture in the fiidze and wait until the fig mixture is cooled.
Using an ice cream machine, pour liquids into the frozen basin and process according to manufacturer instructions, ie., let spin and thicken for 20 minutes before adding mastarpone, fig jan mixture and the nurs. Continue to process for -/-10 minutes, Pour semi froZEn mixture Lnte a pyrex dish.or glass tupperware, Freeze for at least two hours before serving.
Vanessa stevenson
Food enthusisast, photography.Add a pinch of raw foodism and that's pretty much whi I am


AMD Ryzen 9 3950X Processor available starting today (16c/32t)
by Hilbert Hagedoorn
Man, that's a lot of AMD news today. It is (very) fast, feisty, agile and albeit a lot of money, affordable. Ryzen 9 3950X (16c/32t) is priced at 749 USD, you can read our 28 PAGE review here.
-- AMD -- The Ryzen 9 3950X, available for sale today, brings 16-cores and AMD’s industry-leading 7nm process technology together to provide gamers and enthusiasts with the fastest consumer CPU on the market. The 3950X offers exceptional single- and multi-threaded performance for a unique mix of HEDT and consumer-friendly value, delivering incredible gaming experiences and powerful performance for encoding, editing and more.
|
|
|||||
| CPU | Cores / Threads | Clock speed/turbo (GHz) | Cache (total) | PCIe lanes CPU+x570 chiplet) |
MSRP |
| Ryzen 9 3950X | 16/32 | 3.5 / 4.7 | 72MB | 40 | $ 749 |
| Ryzen 9 3900X | 12/24 | 3.8 / 4.6 | 70MB | 40 | $ 499 |
| Ryzen 7 3800X | 8/16 | 3.9 / 4.5 | 36MB | 40 | $ 399 |
| Ryzen 7 3700X | 8/16 | 3.6 / 4.4 | 36MB | 40 | $ 329 |
| Ryzen 5 3600X | 6/12 | 3.8 / 4.4 | 35MB | 40 | $ 249 |
| Ryzen 5 3600 | 6/12 | 3.6 / 4.2 | 35MB | 40 | $ 199 |