Using results from Monte Carlo models for machine learning
-
-
Introduction
-
-
-The aim of this project is to use an already developed Monte Carlo program (either the ising Model or a variational Monte Carlo code) to produce, in case of the Ising model, the energy as function of temperature. For the variational Monte Carlo calculation of interacting electrons in an oscilaltor trap, the data are represented by the ground state energies as functions of the variational parameters.
-
-
-In its simplest form
-the energy of the Ising model is expressed as, without an externally applied magnetic field,
-$$
-E=-J\sum_{< kl >}^{N}s_ks_l
-$$
-
-with
-\( s_k=\pm 1 \). The quantity \( N \) represents the total number of spins and \( J \) is a coupling
-constant expressing the strength of the interaction between
-neighboring spins. The symbol \( < kl> \) indicates that we sum over
-nearest neighbors only. We will assume that we have a ferromagnetic
-ordering, viz \( J> 0 \). We will use periodic boundary conditions and
-the Metropolis algorithm only. Alternatively we can use the supplied variational Monte Carlo program which solves Schroedinger's equation for two interacting electrons in a harmonic oscillator trap. Both codes can be found at the webpage of the course under programs.
-
-
Part a): Producing the data
-
-
-If we opt for the Ising model code, we need to generate for every span over the lattice and output of the total energy and the magnetic moment (if we wish to study this quantity as well), that is compute and write to file as function of temperature
-\( \langle E\rangle \) and \( \langle \vert M\vert \rangle \).
-We will use a fixed lattice siste of \( L\times L = 40 \times 40 \). Make sure the calculations have been equilibrated and compute these expectation values for
-for \( T\in [1.0,3.0] \) with a step in
-temperature \( \Delta T=0.1 \) or smaller. This are the data you will use to estimate the standard deviation in the next part of the project.
-You should keep a file for the different temperatures as you will need to compute the standard deviation for each temperature.
-
-
-Alternatively you can run the variational Monte Carlo program for two interacting electrons confined to move in a harmonic oscillator trap. Here you will produce the variational expectaction value of the energy as a function of the variational parameters \( \alpha \) and \( \beta \). You could single out only a set of such parameters after you have found the minimum energy.
-
-
-These two cases form then our training data which we will use for our estimates of the standard deviation.
-
-
Part b): Estimating the standard deviation
-
-
-This part deals with widely used resampling methods to find the standard deviation on our data sets. The methods we can use are the Blocking method, the Bootstrap method and the Jackknife method. The latter two are examples of so-called resampling methods. The functions provided under the program folder encode all three methods. The typical situation is to generate data using a C++ or Fortan code, as in part a) above. These data are then used by a post-analysis program in order to perform a statistical analysis.
-
-
-In case of the Ising model, the final product will be a table of the expectation values of the energy (or the magnetic moment) as function of temperature and with a proper standard deviation. These data will then enter our analysis in parts c-e). We will call these data the training data.
-
-
-If you opt for the variational Monte Carlo program, you need to make a table of the ground state energy as function of various variational parameters (you can limit yourself to \( \alpha \) only) with proper standard deviation estimations.
-
-
-Your task here is thus to use the data from part a) and generate the standard deviation using one (or more) of the methods discussed above. In general, the Bootstrap method is the most widely used one.
-
-
-These data will then be used in our next step, where we will try to obtain a functional approximation to our data sets.
-
-
Part c): Fitting the data using regression analysis and other methods
-
-
-With the data from part b), with a proper determination of the standard deviation, the task here is to find a fit to the data using regression methods. Here we will focus on the least square method with a straight line fit first and thereafter using singular value decomposition to fit a function that reproduces the data.
-
-
-Since obtaining these data points may not be trivial, we want to use these data to fit a function which can allow us to make predictions for values of \( y \) which are not in the present set. The perhaps simplest approach is to assume we can parametrize our function in terms of a polynomial of degree \( n-1 \) with \( n \) points, that is
-$$
-y=y(x) \rightarrow y(x_i)=\tilde{y}_i+\epsilon_i=\sum_{j=0}^{n-1} \beta_i x_i^j+\epsilon_i,
-$$
-
-where \( \epsilon_i \) is the error in our approximation.
-
-
-In order to find the optimal parameters \( \beta_i \) instead of solving the above linear algebra problem, we define a function which gives a measure of the spread between the values \( y_i \) (which represent hopefully the exact values) and the parametrized values \( \tilde{y}_i \), namely
-$$
-Q(\hat{\beta})=\sum_{i=0}^{n-1}\left(y_i-\tilde{y}_i\right)^2=\left(\hat{y}-\hat{\tilde{y}}\right)^T\left(\hat{y}-\hat{\tilde{y}}\right),
-$$
-
-or using the matrix \( \hat{X} \) as
-$$
-Q(\hat{\beta})=\left(\hat{y}-\hat{X}\hat{\beta}\right)^T\left(\hat{y}-\hat{X}\hat{\beta}\right).
-$$
-
-
-The function
-$$
-Q(\hat{\beta})=\left(\hat{y}-\hat{X}\hat{\beta}\right)^T\left(\hat{y}-\hat{X}\hat{\beta}\right),
-$$
-
-can be linked to the variance of the quantity \( y_i \) if we interpret the latter as the mean value of for example a numerical experiment. When linking below with the maximum likelihood approach below, we will indeed interpret \( y_i \) as a mean value
-$$
-y_{i}=\langle y_i \rangle = \beta_0x_{i,0}+\beta_1x_{i,1}+\beta_2x_{i,2}+\dots+\beta_{n-1}x_{i,n-1}+\epsilon_i,
-$$
-
-where \( \langle y_i \rangle \) is the mean value. Keep in mind also that till now we have treated \( y_i \) as the exact value. Normally, the response (dependent or outcome) variable \( y_i \) the outcome of a numerical experiment or another type of experiment and is thus only an approximation to the true value. It is then always accompanied by an error estimate, often limited to a statistical error estimate given by the standard deviation discussed earlier. In the discussion here we will treat \( y_i \) as our exact value for the response variable.
-
-
-Introducing the standard deviation \( \sigma_i \) for each measurement \( y_i \), we define now the \( \chi^2 \) function as
-$$
-\chi^2(\hat{\beta})=\sum_{i=0}^{n-1}\frac{\left(y_i-\tilde{y}_i\right)^2}{\sigma_i^2}=\left(\hat{y}-\hat{\tilde{y}}\right)^T\frac{1}{\hat{\Sigma^2}}\left(\hat{y}-\hat{\tilde{y}}\right),
-$$
-
-where the matrix \( \hat{\Sigma} \) is a diagonal matrix with \( \sigma_i \) as matrix elements.
-
-
-In order to find the parameters \( \beta_i \) we will then minimize the spread of \( \chi^2(\hat{\beta}) \) by requiring
-$$
-\frac{\partial \chi^2(\hat{\beta})}{\partial \beta_j} = \frac{\partial }{\partial \beta_j}\left[ \sum_{i=0}^{n-1}\left(\frac{y_i-\beta_0x_{i,0}-\beta_1x_{i,1}-\beta_2x_{i,2}-\dots-\beta_{n-1}x_{i,n-1}}{\sigma_i}\right)^2\right]=0,
-$$
-
-which results in
-$$
-\frac{\partial \chi^2(\hat{\beta})}{\partial \beta_j} = -2\left[ \sum_{i=0}^{n-1}\frac{x_{ij}}{\sigma_i}\left(\frac{y_i-\beta_0x_{i,0}-\beta_1x_{i,1}-\beta_2x_{i,2}-\dots-\beta_{n-1}x_{i,n-1}}{\sigma_i}\right)\right]=0,
-$$
-
-or in a matrix-vector form as
-$$
-\frac{\partial \chi^2(\hat{\beta})}{\partial \hat{\beta}} = 0 = \hat{A}^T\left( \hat{b}-\hat{A}\hat{\beta}\right).
-$$
-
-where we have defined the matrix \( \hat{A} =\hat{X}/\hat{\Sigma} \) with matrix elements \( a_{ij} = x_{ij}/\sigma_i \) and the vector \( \hat{b} \) with elements \( b_i = y_i/\sigma_i \).
-
-
-We can rewrite
-$$
-\frac{\partial \chi^2(\hat{\beta})}{\partial \hat{\beta}} = 0 = \hat{A}^T\left( \hat{b}-\hat{A}\hat{\beta}\right),
-$$
-
-as
-$$
-\hat{A}^T\hat{b} = \hat{A}^T\hat{A}\hat{\beta},
-$$
-
-and if the matrix \( \hat{A}^T\hat{A} \) is invertible we have the solution
-$$
-\hat{\beta} =\left(\hat{A}^T\hat{A}\right)^{-1}\hat{A}^T\hat{b}.
-$$
-
-If we then introduce the matrix
-$$
-\hat{H} = \hat{A}^T\hat{A},
-$$
-
-we have then the following expression for the parameters \( \beta_j \) (the matrix elements of \( \hat{H} \) are \( h_{ij} \))
-$$
-\beta_j = \sum_{k=0}^{p-1}h_{jk}\sum_{i=0}^{n-1}\frac{y_i}{\sigma_i}\frac{x_{ik}}{\sigma_i} = \sum_{k=0}^{p-1}h_{jk}\sum_{i=0}^{n-1}b_ia_{ik}
-$$
-
-We state without proof the expression for the uncertainty in the parameters \( \beta_j \) as
-$$
-\sigma^2(\beta_j) = \sum_{i=0}^{n-1}\sigma_i^2\left( \frac{\partial \beta_j}{\partial y_i}\right)^2,
-$$
-
-resulting in
-$$
-\sigma^2(\beta_j) = \left(\sum_{k=0}^{p-1}h_{jk}\sum_{i=0}^{n-1}a_{ik}\right)\left(\sum_{l=0}^{p-1}h_{jl}\sum_{m=0}^{n-1}a_{ml}\right) = h_{jj}!
-$$
-
-
-The first step here is to approximate the function \( y \) with a first-order polynomial, that is we write
-$$
-y=y(x) \rightarrow y(x_i) \approx \beta_0+\beta_1 x_i.
-$$
-
-By computing the derivatives of \( \chi^2 \) with respect to \( \beta_0 \) and \( \beta_1 \) show that these are given by
-$$
-\frac{\partial \chi^2(\hat{\beta})}{\partial \beta_0} = -2\left[ \sum_{i=0}^{1}\left(\frac{y_i-\beta_0-\beta_1x_{i}}{\sigma_i^2}\right)\right]=0,
-$$
-
-and
-$$
-\frac{\partial \chi^2(\hat{\beta})}{\partial \beta_0} = -2\left[ \sum_{i=0}^{1}x_i\left(\frac{y_i-\beta_0-\beta_1x_{i}}{\sigma_i^2}\right)\right]=0.
-$$
-
-Define then
-$$
-\gamma = \sum_{i=0}^{1}\frac{1}{\sigma_i^2},
-$$
-
-
-$$
-\gamma_x = \sum_{i=0}^{1}\frac{x_{i}}{\sigma_i^2},
-$$
-
-$$
-\gamma_y = \sum_{i=0}^{1}\left(\frac{y_i}{\sigma_i^2}\right),
-$$
-
-$$
-\gamma_{xx} = \sum_{i=0}^{1}\frac{x_ix_{i}}{\sigma_i^2},
-$$
-
-$$
-\gamma_{xy} = \sum_{i=0}^{1}\frac{y_ix_{i}}{\sigma_i^2},
-$$
-
-and show that
-$$
-\beta_0 = \frac{\gamma_{xx}\gamma_y-\gamma_x\gamma_y}{\gamma\gamma_{xx}-\gamma_x^2},
-$$
-
-$$
-\beta_1 = \frac{\gamma_{xy}\gamma-\gamma_x\gamma_y}{\gamma\gamma_{xx}-\gamma_x^2}.
-$$
-
-Find these parameters for your data set and compare the fitted function \( y=y(x) \rightarrow y(x_i) \approx \beta_0+\beta_1 x_i \) and find the error in the parameters \( \beta_0 \) and \( \beta_1 \) as well.
-How does your parametrization compare with the input data \( y_i \)? Give a critical discussion of the method.
-
-
-Develop now a program where you extract the parameters \( \beta_i \) and their corresponding errors using the general Least Square method (LSM), with polynomials of degree \( p=2 \) and \( p=3 \). Which polyomials gives the best fit to your data? Here you should prepare an input file with the number of data points in your training data set, the corresponding standard deviation \( \sigma_i \) as well as the degree \( p \) of the polynomial you want to fit.
-
-
-The LSM suffers often from both being underdetermined and overdetermined in the unknown coefficients \( \beta_i \). A better approach is to use the Singular Value Decomposition (SVD) method discussed in the lecture notes on regression analysis. Here you should use the SVD programs and use the same input file you used for the general LSM. In the main functions you should specifiy which function you wish to approximate the date with. The example programs discussed above show you examples on how to use these codes as well as how to write the functional form of the polynomial or fitting function.
-
-
Part d): Introducing Bayesian statistics
-
-More text to come
-
-
Part e): Studying the Ising model or the VMC results with Neural networks
-
-More text to come
-
-
Background literature
-
-
-If you wish to read more about the Ising model and statistical physics here are three suggestions.
-
-
-Here follows a brief recipe and recommendation on how to write a report for each
-project.
-
-
-
Give a short description of the nature of the problem and the eventual numerical methods you have used.
-
Describe the algorithm you have used and/or developed. Here you may find it convenient to use pseudocoding. In many cases you can describe the algorithm in the program itself.
-
Include the source code of your program. Comment your program properly.
-
If possible, try to find analytic solutions, or known limits in order to test your program when developing the code.
-
Include your results either in figure form or in a table. Remember to label your results. All tables and figures should have relevant captions and labels on the axes.
-
Try to evaluate the reliabilty and numerical stability/precision of your results. If possible, include a qualitative and/or quantitative discussion of the numerical stability, eventual loss of precision etc.
-
Try to give an interpretation of you results in your answers to the problems.
-
Critique: if possible include your comments and reflections about the exercise, whether you felt you learnt something, ideas for improvements and other thoughts you've made when solving the exercise. We wish to keep this course at the interactive level and your comments can help us improve it.
-
Try to establish a practice where you log your work at the computerlab. You may find such a logbook very handy at later stages in your work, especially when you don't properly remember what a previous test version of your program did. Here you could also record the time spent on solving the exercise, various algorithms you may have tested or other topics which you feel worthy of mentioning.
Using results from Monte Carlo models for machine learning
-
-
Introduction
-
-
-The aim of this project is to use an already developed Monte Carlo program (either the ising Model or a variational Monte Carlo code) to produce, in case of the Ising model, the energy as function of temperature. For the variational Monte Carlo calculation of interacting electrons in an oscilaltor trap, the data are represented by the ground state energies as functions of the variational parameters.
-
-
-In its simplest form
-the energy of the Ising model is expressed as, without an externally applied magnetic field,
-$$
-E=-J\sum_{< kl >}^{N}s_ks_l
-$$
-
-with
-\( s_k=\pm 1 \). The quantity \( N \) represents the total number of spins and \( J \) is a coupling
-constant expressing the strength of the interaction between
-neighboring spins. The symbol \( < kl> \) indicates that we sum over
-nearest neighbors only. We will assume that we have a ferromagnetic
-ordering, viz \( J> 0 \). We will use periodic boundary conditions and
-the Metropolis algorithm only. Alternatively we can use the supplied variational Monte Carlo program which solves Schroedinger's equation for two interacting electrons in a harmonic oscillator trap. Both codes can be found at the webpage of the course under programs.
-
-
Part a): Producing the data
-
-
-If we opt for the Ising model code, we need to generate for every span over the lattice and output of the total energy and the magnetic moment (if we wish to study this quantity as well), that is compute and write to file as function of temperature
-\( \langle E\rangle \) and \( \langle \vert M\vert \rangle \).
-We will use a fixed lattice siste of \( L\times L = 40 \times 40 \). Make sure the calculations have been equilibrated and compute these expectation values for
-for \( T\in [1.0,3.0] \) with a step in
-temperature \( \Delta T=0.1 \) or smaller. This are the data you will use to estimate the standard deviation in the next part of the project.
-You should keep a file for the different temperatures as you will need to compute the standard deviation for each temperature.
-
-
-Alternatively you can run the variational Monte Carlo program for two interacting electrons confined to move in a harmonic oscillator trap. Here you will produce the variational expectaction value of the energy as a function of the variational parameters \( \alpha \) and \( \beta \). You could single out only a set of such parameters after you have found the minimum energy.
-
-
-These two cases form then our training data which we will use for our estimates of the standard deviation.
-
-
Part b): Estimating the standard deviation
-
-
-This part deals with widely used resampling methods to find the standard deviation on our data sets. The methods we can use are the Blocking method, the Bootstrap method and the Jackknife method. The latter two are examples of so-called resampling methods. The functions provided under the program folder encode all three methods. The typical situation is to generate data using a C++ or Fortan code, as in part a) above. These data are then used by a post-analysis program in order to perform a statistical analysis.
-
-
-In case of the Ising model, the final product will be a table of the expectation values of the energy (or the magnetic moment) as function of temperature and with a proper standard deviation. These data will then enter our analysis in parts c-e). We will call these data the training data.
-
-
-If you opt for the variational Monte Carlo program, you need to make a table of the ground state energy as function of various variational parameters (you can limit yourself to \( \alpha \) only) with proper standard deviation estimations.
-
-
-Your task here is thus to use the data from part a) and generate the standard deviation using one (or more) of the methods discussed above. In general, the Bootstrap method is the most widely used one.
-
-
-These data will then be used in our next step, where we will try to obtain a functional approximation to our data sets.
-
-
Part c): Fitting the data using regression analysis and other methods
-
-
-With the data from part b), with a proper determination of the standard deviation, the task here is to find a fit to the data using regression methods. Here we will focus on the least square method with a straight line fit first and thereafter using singular value decomposition to fit a function that reproduces the data.
-
-
-Since obtaining these data points may not be trivial, we want to use these data to fit a function which can allow us to make predictions for values of \( y \) which are not in the present set. The perhaps simplest approach is to assume we can parametrize our function in terms of a polynomial of degree \( n-1 \) with \( n \) points, that is
-$$
-y=y(x) \rightarrow y(x_i)=\tilde{y}_i+\epsilon_i=\sum_{j=0}^{n-1} \beta_i x_i^j+\epsilon_i,
-$$
-
-where \( \epsilon_i \) is the error in our approximation.
-
-
-In order to find the optimal parameters \( \beta_i \) instead of solving the above linear algebra problem, we define a function which gives a measure of the spread between the values \( y_i \) (which represent hopefully the exact values) and the parametrized values \( \tilde{y}_i \), namely
-$$
-Q(\hat{\beta})=\sum_{i=0}^{n-1}\left(y_i-\tilde{y}_i\right)^2=\left(\hat{y}-\hat{\tilde{y}}\right)^T\left(\hat{y}-\hat{\tilde{y}}\right),
-$$
-
-or using the matrix \( \hat{X} \) as
-$$
-Q(\hat{\beta})=\left(\hat{y}-\hat{X}\hat{\beta}\right)^T\left(\hat{y}-\hat{X}\hat{\beta}\right).
-$$
-
-
-The function
-$$
-Q(\hat{\beta})=\left(\hat{y}-\hat{X}\hat{\beta}\right)^T\left(\hat{y}-\hat{X}\hat{\beta}\right),
-$$
-
-can be linked to the variance of the quantity \( y_i \) if we interpret the latter as the mean value of for example a numerical experiment. When linking below with the maximum likelihood approach below, we will indeed interpret \( y_i \) as a mean value
-$$
-y_{i}=\langle y_i \rangle = \beta_0x_{i,0}+\beta_1x_{i,1}+\beta_2x_{i,2}+\dots+\beta_{n-1}x_{i,n-1}+\epsilon_i,
-$$
-
-where \( \langle y_i \rangle \) is the mean value. Keep in mind also that till now we have treated \( y_i \) as the exact value. Normally, the response (dependent or outcome) variable \( y_i \) the outcome of a numerical experiment or another type of experiment and is thus only an approximation to the true value. It is then always accompanied by an error estimate, often limited to a statistical error estimate given by the standard deviation discussed earlier. In the discussion here we will treat \( y_i \) as our exact value for the response variable.
-
-
-Introducing the standard deviation \( \sigma_i \) for each measurement \( y_i \), we define now the \( \chi^2 \) function as
-$$
-\chi^2(\hat{\beta})=\sum_{i=0}^{n-1}\frac{\left(y_i-\tilde{y}_i\right)^2}{\sigma_i^2}=\left(\hat{y}-\hat{\tilde{y}}\right)^T\frac{1}{\hat{\Sigma^2}}\left(\hat{y}-\hat{\tilde{y}}\right),
-$$
-
-where the matrix \( \hat{\Sigma} \) is a diagonal matrix with \( \sigma_i \) as matrix elements.
-
-
-In order to find the parameters \( \beta_i \) we will then minimize the spread of \( \chi^2(\hat{\beta}) \) by requiring
-$$
-\frac{\partial \chi^2(\hat{\beta})}{\partial \beta_j} = \frac{\partial }{\partial \beta_j}\left[ \sum_{i=0}^{n-1}\left(\frac{y_i-\beta_0x_{i,0}-\beta_1x_{i,1}-\beta_2x_{i,2}-\dots-\beta_{n-1}x_{i,n-1}}{\sigma_i}\right)^2\right]=0,
-$$
-
-which results in
-$$
-\frac{\partial \chi^2(\hat{\beta})}{\partial \beta_j} = -2\left[ \sum_{i=0}^{n-1}\frac{x_{ij}}{\sigma_i}\left(\frac{y_i-\beta_0x_{i,0}-\beta_1x_{i,1}-\beta_2x_{i,2}-\dots-\beta_{n-1}x_{i,n-1}}{\sigma_i}\right)\right]=0,
-$$
-
-or in a matrix-vector form as
-$$
-\frac{\partial \chi^2(\hat{\beta})}{\partial \hat{\beta}} = 0 = \hat{A}^T\left( \hat{b}-\hat{A}\hat{\beta}\right).
-$$
-
-where we have defined the matrix \( \hat{A} =\hat{X}/\hat{\Sigma} \) with matrix elements \( a_{ij} = x_{ij}/\sigma_i \) and the vector \( \hat{b} \) with elements \( b_i = y_i/\sigma_i \).
-
-
-We can rewrite
-$$
-\frac{\partial \chi^2(\hat{\beta})}{\partial \hat{\beta}} = 0 = \hat{A}^T\left( \hat{b}-\hat{A}\hat{\beta}\right),
-$$
-
-as
-$$
-\hat{A}^T\hat{b} = \hat{A}^T\hat{A}\hat{\beta},
-$$
-
-and if the matrix \( \hat{A}^T\hat{A} \) is invertible we have the solution
-$$
-\hat{\beta} =\left(\hat{A}^T\hat{A}\right)^{-1}\hat{A}^T\hat{b}.
-$$
-
-If we then introduce the matrix
-$$
-\hat{H} = \hat{A}^T\hat{A},
-$$
-
-we have then the following expression for the parameters \( \beta_j \) (the matrix elements of \( \hat{H} \) are \( h_{ij} \))
-$$
-\beta_j = \sum_{k=0}^{p-1}h_{jk}\sum_{i=0}^{n-1}\frac{y_i}{\sigma_i}\frac{x_{ik}}{\sigma_i} = \sum_{k=0}^{p-1}h_{jk}\sum_{i=0}^{n-1}b_ia_{ik}
-$$
-
-We state without proof the expression for the uncertainty in the parameters \( \beta_j \) as
-$$
-\sigma^2(\beta_j) = \sum_{i=0}^{n-1}\sigma_i^2\left( \frac{\partial \beta_j}{\partial y_i}\right)^2,
-$$
-
-resulting in
-$$
-\sigma^2(\beta_j) = \left(\sum_{k=0}^{p-1}h_{jk}\sum_{i=0}^{n-1}a_{ik}\right)\left(\sum_{l=0}^{p-1}h_{jl}\sum_{m=0}^{n-1}a_{ml}\right) = h_{jj}!
-$$
-
-
-The first step here is to approximate the function \( y \) with a first-order polynomial, that is we write
-$$
-y=y(x) \rightarrow y(x_i) \approx \beta_0+\beta_1 x_i.
-$$
-
-By computing the derivatives of \( \chi^2 \) with respect to \( \beta_0 \) and \( \beta_1 \) show that these are given by
-$$
-\frac{\partial \chi^2(\hat{\beta})}{\partial \beta_0} = -2\left[ \sum_{i=0}^{1}\left(\frac{y_i-\beta_0-\beta_1x_{i}}{\sigma_i^2}\right)\right]=0,
-$$
-
-and
-$$
-\frac{\partial \chi^2(\hat{\beta})}{\partial \beta_0} = -2\left[ \sum_{i=0}^{1}x_i\left(\frac{y_i-\beta_0-\beta_1x_{i}}{\sigma_i^2}\right)\right]=0.
-$$
-
-Define then
-$$
-\gamma = \sum_{i=0}^{1}\frac{1}{\sigma_i^2},
-$$
-
-
-$$
-\gamma_x = \sum_{i=0}^{1}\frac{x_{i}}{\sigma_i^2},
-$$
-
-$$
-\gamma_y = \sum_{i=0}^{1}\left(\frac{y_i}{\sigma_i^2}\right),
-$$
-
-$$
-\gamma_{xx} = \sum_{i=0}^{1}\frac{x_ix_{i}}{\sigma_i^2},
-$$
-
-$$
-\gamma_{xy} = \sum_{i=0}^{1}\frac{y_ix_{i}}{\sigma_i^2},
-$$
-
-and show that
-$$
-\beta_0 = \frac{\gamma_{xx}\gamma_y-\gamma_x\gamma_y}{\gamma\gamma_{xx}-\gamma_x^2},
-$$
-
-$$
-\beta_1 = \frac{\gamma_{xy}\gamma-\gamma_x\gamma_y}{\gamma\gamma_{xx}-\gamma_x^2}.
-$$
-
-Find these parameters for your data set and compare the fitted function \( y=y(x) \rightarrow y(x_i) \approx \beta_0+\beta_1 x_i \) and find the error in the parameters \( \beta_0 \) and \( \beta_1 \) as well.
-How does your parametrization compare with the input data \( y_i \)? Give a critical discussion of the method.
-
-
-Develop now a program where you extract the parameters \( \beta_i \) and their corresponding errors using the general Least Square method (LSM), with polynomials of degree \( p=2 \) and \( p=3 \). Which polyomials gives the best fit to your data? Here you should prepare an input file with the number of data points in your training data set, the corresponding standard deviation \( \sigma_i \) as well as the degree \( p \) of the polynomial you want to fit.
-
-
-The LSM suffers often from both being underdetermined and overdetermined in the unknown coefficients \( \beta_i \). A better approach is to use the Singular Value Decomposition (SVD) method discussed in the lecture notes on regression analysis. Here you should use the SVD programs and use the same input file you used for the general LSM. In the main functions you should specifiy which function you wish to approximate the date with. The example programs discussed above show you examples on how to use these codes as well as how to write the functional form of the polynomial or fitting function.
-
-
Part d): Introducing Bayesian statistics
-
-More text to come
-
-
Part e): Studying the Ising model or the VMC results with Neural networks
-
-More text to come
-
-
Background literature
-
-
-If you wish to read more about the Ising model and statistical physics here are three suggestions.
-
-
-Here follows a brief recipe and recommendation on how to write a report for each
-project.
-
-
-
Give a short description of the nature of the problem and the eventual numerical methods you have used.
-
Describe the algorithm you have used and/or developed. Here you may find it convenient to use pseudocoding. In many cases you can describe the algorithm in the program itself.
-
Include the source code of your program. Comment your program properly.
-
If possible, try to find analytic solutions, or known limits in order to test your program when developing the code.
-
Include your results either in figure form or in a table. Remember to label your results. All tables and figures should have relevant captions and labels on the axes.
-
Try to evaluate the reliabilty and numerical stability/precision of your results. If possible, include a qualitative and/or quantitative discussion of the numerical stability, eventual loss of precision etc.
-
Try to give an interpretation of you results in your answers to the problems.
-
Critique: if possible include your comments and reflections about the exercise, whether you felt you learnt something, ideas for improvements and other thoughts you've made when solving the exercise. We wish to keep this course at the interactive level and your comments can help us improve it.
-
Try to establish a practice where you log your work at the computerlab. You may find such a logbook very handy at later stages in your work, especially when you don't properly remember what a previous test version of your program did. Here you could also record the time spent on solving the exercise, various algorithms you may have tested or other topics which you feel worthy of mentioning.
Using results from Monte Carlo models for machine learning
-
-
Introduction
-
-
-The aim of this project is to use an already developed Monte Carlo program (either the ising Model or a variational Monte Carlo code) to produce, in case of the Ising model, the energy as function of temperature. For the variational Monte Carlo calculation of interacting electrons in an oscilaltor trap, the data are represented by the ground state energies as functions of the variational parameters.
-
-
-In its simplest form
-the energy of the Ising model is expressed as, without an externally applied magnetic field,
-$$
-E=-J\sum_{< kl >}^{N}s_ks_l
-$$
-
-with
-\( s_k=\pm 1 \). The quantity \( N \) represents the total number of spins and \( J \) is a coupling
-constant expressing the strength of the interaction between
-neighboring spins. The symbol \( < kl> \) indicates that we sum over
-nearest neighbors only. We will assume that we have a ferromagnetic
-ordering, viz \( J> 0 \). We will use periodic boundary conditions and
-the Metropolis algorithm only. Alternatively we can use the supplied variational Monte Carlo program which solves Schroedinger's equation for two interacting electrons in a harmonic oscillator trap. Both codes can be found at the webpage of the course under programs.
-
-
Part a): Producing the data
-
-
-If we opt for the Ising model code, we need to generate for every span over the lattice and output of the total energy and the magnetic moment (if we wish to study this quantity as well), that is compute and write to file as function of temperature
-\( \langle E\rangle \) and \( \langle \vert M\vert \rangle \).
-We will use a fixed lattice siste of \( L\times L = 40 \times 40 \). Make sure the calculations have been equilibrated and compute these expectation values for
-for \( T\in [1.0,3.0] \) with a step in
-temperature \( \Delta T=0.1 \) or smaller. This are the data you will use to estimate the standard deviation in the next part of the project.
-You should keep a file for the different temperatures as you will need to compute the standard deviation for each temperature.
-
-
-Alternatively you can run the variational Monte Carlo program for two interacting electrons confined to move in a harmonic oscillator trap. Here you will produce the variational expectaction value of the energy as a function of the variational parameters \( \alpha \) and \( \beta \). You could single out only a set of such parameters after you have found the minimum energy.
-
-
-These two cases form then our training data which we will use for our estimates of the standard deviation.
-
-
Part b): Estimating the standard deviation
-
-
-This part deals with widely used resampling methods to find the standard deviation on our data sets. The methods we can use are the Blocking method, the Bootstrap method and the Jackknife method. The latter two are examples of so-called resampling methods. The functions provided under the program folder encode all three methods. The typical situation is to generate data using a C++ or Fortan code, as in part a) above. These data are then used by a post-analysis program in order to perform a statistical analysis.
-
-
-In case of the Ising model, the final product will be a table of the expectation values of the energy (or the magnetic moment) as function of temperature and with a proper standard deviation. These data will then enter our analysis in parts c-e). We will call these data the training data.
-
-
-If you opt for the variational Monte Carlo program, you need to make a table of the ground state energy as function of various variational parameters (you can limit yourself to \( \alpha \) only) with proper standard deviation estimations.
-
-
-Your task here is thus to use the data from part a) and generate the standard deviation using one (or more) of the methods discussed above. In general, the Bootstrap method is the most widely used one.
-
-
-These data will then be used in our next step, where we will try to obtain a functional approximation to our data sets.
-
-
Part c): Fitting the data using regression analysis and other methods
-
-
-With the data from part b), with a proper determination of the standard deviation, the task here is to find a fit to the data using regression methods. Here we will focus on the least square method with a straight line fit first and thereafter using singular value decomposition to fit a function that reproduces the data.
-
-
-Since obtaining these data points may not be trivial, we want to use these data to fit a function which can allow us to make predictions for values of \( y \) which are not in the present set. The perhaps simplest approach is to assume we can parametrize our function in terms of a polynomial of degree \( n-1 \) with \( n \) points, that is
-$$
-y=y(x) \rightarrow y(x_i)=\tilde{y}_i+\epsilon_i=\sum_{j=0}^{n-1} \beta_i x_i^j+\epsilon_i,
-$$
-
-where \( \epsilon_i \) is the error in our approximation.
-
-
-In order to find the optimal parameters \( \beta_i \) instead of solving the above linear algebra problem, we define a function which gives a measure of the spread between the values \( y_i \) (which represent hopefully the exact values) and the parametrized values \( \tilde{y}_i \), namely
-$$
-Q(\hat{\beta})=\sum_{i=0}^{n-1}\left(y_i-\tilde{y}_i\right)^2=\left(\hat{y}-\hat{\tilde{y}}\right)^T\left(\hat{y}-\hat{\tilde{y}}\right),
-$$
-
-or using the matrix \( \hat{X} \) as
-$$
-Q(\hat{\beta})=\left(\hat{y}-\hat{X}\hat{\beta}\right)^T\left(\hat{y}-\hat{X}\hat{\beta}\right).
-$$
-
-
-The function
-$$
-Q(\hat{\beta})=\left(\hat{y}-\hat{X}\hat{\beta}\right)^T\left(\hat{y}-\hat{X}\hat{\beta}\right),
-$$
-
-can be linked to the variance of the quantity \( y_i \) if we interpret the latter as the mean value of for example a numerical experiment. When linking below with the maximum likelihood approach below, we will indeed interpret \( y_i \) as a mean value
-$$
-y_{i}=\langle y_i \rangle = \beta_0x_{i,0}+\beta_1x_{i,1}+\beta_2x_{i,2}+\dots+\beta_{n-1}x_{i,n-1}+\epsilon_i,
-$$
-
-where \( \langle y_i \rangle \) is the mean value. Keep in mind also that till now we have treated \( y_i \) as the exact value. Normally, the response (dependent or outcome) variable \( y_i \) the outcome of a numerical experiment or another type of experiment and is thus only an approximation to the true value. It is then always accompanied by an error estimate, often limited to a statistical error estimate given by the standard deviation discussed earlier. In the discussion here we will treat \( y_i \) as our exact value for the response variable.
-
-
-Introducing the standard deviation \( \sigma_i \) for each measurement \( y_i \), we define now the \( \chi^2 \) function as
-$$
-\chi^2(\hat{\beta})=\sum_{i=0}^{n-1}\frac{\left(y_i-\tilde{y}_i\right)^2}{\sigma_i^2}=\left(\hat{y}-\hat{\tilde{y}}\right)^T\frac{1}{\hat{\Sigma^2}}\left(\hat{y}-\hat{\tilde{y}}\right),
-$$
-
-where the matrix \( \hat{\Sigma} \) is a diagonal matrix with \( \sigma_i \) as matrix elements.
-
-
-In order to find the parameters \( \beta_i \) we will then minimize the spread of \( \chi^2(\hat{\beta}) \) by requiring
-$$
-\frac{\partial \chi^2(\hat{\beta})}{\partial \beta_j} = \frac{\partial }{\partial \beta_j}\left[ \sum_{i=0}^{n-1}\left(\frac{y_i-\beta_0x_{i,0}-\beta_1x_{i,1}-\beta_2x_{i,2}-\dots-\beta_{n-1}x_{i,n-1}}{\sigma_i}\right)^2\right]=0,
-$$
-
-which results in
-$$
-\frac{\partial \chi^2(\hat{\beta})}{\partial \beta_j} = -2\left[ \sum_{i=0}^{n-1}\frac{x_{ij}}{\sigma_i}\left(\frac{y_i-\beta_0x_{i,0}-\beta_1x_{i,1}-\beta_2x_{i,2}-\dots-\beta_{n-1}x_{i,n-1}}{\sigma_i}\right)\right]=0,
-$$
-
-or in a matrix-vector form as
-$$
-\frac{\partial \chi^2(\hat{\beta})}{\partial \hat{\beta}} = 0 = \hat{A}^T\left( \hat{b}-\hat{A}\hat{\beta}\right).
-$$
-
-where we have defined the matrix \( \hat{A} =\hat{X}/\hat{\Sigma} \) with matrix elements \( a_{ij} = x_{ij}/\sigma_i \) and the vector \( \hat{b} \) with elements \( b_i = y_i/\sigma_i \).
-
-
-We can rewrite
-$$
-\frac{\partial \chi^2(\hat{\beta})}{\partial \hat{\beta}} = 0 = \hat{A}^T\left( \hat{b}-\hat{A}\hat{\beta}\right),
-$$
-
-as
-$$
-\hat{A}^T\hat{b} = \hat{A}^T\hat{A}\hat{\beta},
-$$
-
-and if the matrix \( \hat{A}^T\hat{A} \) is invertible we have the solution
-$$
-\hat{\beta} =\left(\hat{A}^T\hat{A}\right)^{-1}\hat{A}^T\hat{b}.
-$$
-
-If we then introduce the matrix
-$$
-\hat{H} = \hat{A}^T\hat{A},
-$$
-
-we have then the following expression for the parameters \( \beta_j \) (the matrix elements of \( \hat{H} \) are \( h_{ij} \))
-$$
-\beta_j = \sum_{k=0}^{p-1}h_{jk}\sum_{i=0}^{n-1}\frac{y_i}{\sigma_i}\frac{x_{ik}}{\sigma_i} = \sum_{k=0}^{p-1}h_{jk}\sum_{i=0}^{n-1}b_ia_{ik}
-$$
-
-We state without proof the expression for the uncertainty in the parameters \( \beta_j \) as
-$$
-\sigma^2(\beta_j) = \sum_{i=0}^{n-1}\sigma_i^2\left( \frac{\partial \beta_j}{\partial y_i}\right)^2,
-$$
-
-resulting in
-$$
-\sigma^2(\beta_j) = \left(\sum_{k=0}^{p-1}h_{jk}\sum_{i=0}^{n-1}a_{ik}\right)\left(\sum_{l=0}^{p-1}h_{jl}\sum_{m=0}^{n-1}a_{ml}\right) = h_{jj}!
-$$
-
-
-The first step here is to approximate the function \( y \) with a first-order polynomial, that is we write
-$$
-y=y(x) \rightarrow y(x_i) \approx \beta_0+\beta_1 x_i.
-$$
-
-By computing the derivatives of \( \chi^2 \) with respect to \( \beta_0 \) and \( \beta_1 \) show that these are given by
-$$
-\frac{\partial \chi^2(\hat{\beta})}{\partial \beta_0} = -2\left[ \sum_{i=0}^{1}\left(\frac{y_i-\beta_0-\beta_1x_{i}}{\sigma_i^2}\right)\right]=0,
-$$
-
-and
-$$
-\frac{\partial \chi^2(\hat{\beta})}{\partial \beta_0} = -2\left[ \sum_{i=0}^{1}x_i\left(\frac{y_i-\beta_0-\beta_1x_{i}}{\sigma_i^2}\right)\right]=0.
-$$
-
-Define then
-$$
-\gamma = \sum_{i=0}^{1}\frac{1}{\sigma_i^2},
-$$
-
-
-$$
-\gamma_x = \sum_{i=0}^{1}\frac{x_{i}}{\sigma_i^2},
-$$
-
-$$
-\gamma_y = \sum_{i=0}^{1}\left(\frac{y_i}{\sigma_i^2}\right),
-$$
-
-$$
-\gamma_{xx} = \sum_{i=0}^{1}\frac{x_ix_{i}}{\sigma_i^2},
-$$
-
-$$
-\gamma_{xy} = \sum_{i=0}^{1}\frac{y_ix_{i}}{\sigma_i^2},
-$$
-
-and show that
-$$
-\beta_0 = \frac{\gamma_{xx}\gamma_y-\gamma_x\gamma_y}{\gamma\gamma_{xx}-\gamma_x^2},
-$$
-
-$$
-\beta_1 = \frac{\gamma_{xy}\gamma-\gamma_x\gamma_y}{\gamma\gamma_{xx}-\gamma_x^2}.
-$$
-
-Find these parameters for your data set and compare the fitted function \( y=y(x) \rightarrow y(x_i) \approx \beta_0+\beta_1 x_i \) and find the error in the parameters \( \beta_0 \) and \( \beta_1 \) as well.
-How does your parametrization compare with the input data \( y_i \)? Give a critical discussion of the method.
-
-
-Develop now a program where you extract the parameters \( \beta_i \) and their corresponding errors using the general Least Square method (LSM), with polynomials of degree \( p=2 \) and \( p=3 \). Which polyomials gives the best fit to your data? Here you should prepare an input file with the number of data points in your training data set, the corresponding standard deviation \( \sigma_i \) as well as the degree \( p \) of the polynomial you want to fit.
-
-
-The LSM suffers often from both being underdetermined and overdetermined in the unknown coefficients \( \beta_i \). A better approach is to use the Singular Value Decomposition (SVD) method discussed in the lecture notes on regression analysis. Here you should use the SVD programs and use the same input file you used for the general LSM. In the main functions you should specifiy which function you wish to approximate the date with. The example programs discussed above show you examples on how to use these codes as well as how to write the functional form of the polynomial or fitting function.
-
-
Part d): Introducing Bayesian statistics
-
-More text to come
-
-
Part e): Studying the Ising model or the VMC results with Neural networks
-
-More text to come
-
-
Background literature
-
-
-If you wish to read more about the Ising model and statistical physics here are three suggestions.
-
-
-Here follows a brief recipe and recommendation on how to write a report for each
-project.
-
-
-
Give a short description of the nature of the problem and the eventual numerical methods you have used.
-
Describe the algorithm you have used and/or developed. Here you may find it convenient to use pseudocoding. In many cases you can describe the algorithm in the program itself.
-
Include the source code of your program. Comment your program properly.
-
If possible, try to find analytic solutions, or known limits in order to test your program when developing the code.
-
Include your results either in figure form or in a table. Remember to label your results. All tables and figures should have relevant captions and labels on the axes.
-
Try to evaluate the reliabilty and numerical stability/precision of your results. If possible, include a qualitative and/or quantitative discussion of the numerical stability, eventual loss of precision etc.
-
Try to give an interpretation of you results in your answers to the problems.
-
Critique: if possible include your comments and reflections about the exercise, whether you felt you learnt something, ideas for improvements and other thoughts you've made when solving the exercise. We wish to keep this course at the interactive level and your comments can help us improve it.
-
Try to establish a practice where you log your work at the computerlab. You may find such a logbook very handy at later stages in your work, especially when you don't properly remember what a previous test version of your program did. Here you could also record the time spent on solving the exercise, various algorithms you may have tested or other topics which you feel worthy of mentioning.
-
-
-
-
-
-
diff --git a/doc/Projects/2017/Project/ipynb b/doc/Projects/2017/Project/ipynb
deleted file mode 100644
index cd0e9b8c3..000000000
Binary files a/doc/Projects/2017/Project/ipynb and /dev/null differ
diff --git a/doc/Projects/2017/Project/pdf/Project.p.tex b/doc/Projects/2017/Project/pdf/Project.p.tex
deleted file mode 100644
index 9154847f6..000000000
--- a/doc/Projects/2017/Project/pdf/Project.p.tex
+++ /dev/null
@@ -1,379 +0,0 @@
-%%
-%% Automatically generated file from DocOnce source
-%% (https://github.com/hplgit/doconce/)
-%%
-%%
-% #ifdef PTEX2TEX_EXPLANATION
-%%
-%% The file follows the ptex2tex extended LaTeX format, see
-%% ptex2tex: http://code.google.com/p/ptex2tex/
-%%
-%% Run
-%% ptex2tex myfile
-%% or
-%% doconce ptex2tex myfile
-%%
-%% to turn myfile.p.tex into an ordinary LaTeX file myfile.tex.
-%% (The ptex2tex program: http://code.google.com/p/ptex2tex)
-%% Many preprocess options can be added to ptex2tex or doconce ptex2tex
-%%
-%% ptex2tex -DMINTED myfile
-%% doconce ptex2tex myfile envir=minted
-%%
-%% ptex2tex will typeset code environments according to a global or local
-%% .ptex2tex.cfg configure file. doconce ptex2tex will typeset code
-%% according to options on the command line (just type doconce ptex2tex to
-%% see examples). If doconce ptex2tex has envir=minted, it enables the
-%% minted style without needing -DMINTED.
-% #endif
-
-% #define PREAMBLE
-
-% #ifdef PREAMBLE
-%-------------------- begin preamble ----------------------
-
-\documentclass[%
-oneside, % oneside: electronic viewing, twoside: printing
-final, % draft: marks overfull hboxes, figures with paths
-10pt]{article}
-
-\listfiles % print all files needed to compile this document
-
-\usepackage{relsize,makeidx,color,setspace,amsmath,amsfonts,amssymb}
-\usepackage[table]{xcolor}
-\usepackage{bm,ltablex,microtype}
-
-\usepackage[pdftex]{graphicx}
-
-\usepackage[T1]{fontenc}
-%\usepackage[latin1]{inputenc}
-\usepackage{ucs}
-\usepackage[utf8x]{inputenc}
-
-\usepackage{lmodern} % Latin Modern fonts derived from Computer Modern
-
-% Hyperlinks in PDF:
-\definecolor{linkcolor}{rgb}{0,0,0.4}
-\usepackage{hyperref}
-\hypersetup{
- breaklinks=true,
- colorlinks=true,
- linkcolor=linkcolor,
- urlcolor=linkcolor,
- citecolor=black,
- filecolor=black,
- %filecolor=blue,
- pdfmenubar=true,
- pdftoolbar=true,
- bookmarksdepth=3 % Uncomment (and tweak) for PDF bookmarks with more levels than the TOC
- }
-%\hyperbaseurl{} % hyperlinks are relative to this root
-
-\setcounter{tocdepth}{2} % levels in table of contents
-
-% --- fancyhdr package for fancy headers ---
-\usepackage{fancyhdr}
-\fancyhf{} % sets both header and footer to nothing
-\renewcommand{\headrulewidth}{0pt}
-\fancyfoot[LE,RO]{\thepage}
-% Ensure copyright on titlepage (article style) and chapter pages (book style)
-\fancypagestyle{plain}{
- \fancyhf{}
- \fancyfoot[C]{{\footnotesize \copyright\ 1999-2017, "Data Analysis and Machine Learning FYS-MAT3155/FYS4155":"http://www.uio.no/studier/emner/matnat/fys/FYS3155/index-eng.html". Released under CC Attribution-NonCommercial 4.0 license}}
-% \renewcommand{\footrulewidth}{0mm}
- \renewcommand{\headrulewidth}{0mm}
-}
-% Ensure copyright on titlepages with \thispagestyle{empty}
-\fancypagestyle{empty}{
- \fancyhf{}
- \fancyfoot[C]{{\footnotesize \copyright\ 1999-2017, "Data Analysis and Machine Learning FYS-MAT3155/FYS4155":"http://www.uio.no/studier/emner/matnat/fys/FYS3155/index-eng.html". Released under CC Attribution-NonCommercial 4.0 license}}
- \renewcommand{\footrulewidth}{0mm}
- \renewcommand{\headrulewidth}{0mm}
-}
-
-\pagestyle{fancy}
-
-
-% prevent orhpans and widows
-\clubpenalty = 10000
-\widowpenalty = 10000
-
-% --- end of standard preamble for documents ---
-
-
-% insert custom LaTeX commands...
-
-\raggedbottom
-\makeindex
-\usepackage[totoc]{idxlayout} % for index in the toc
-\usepackage[nottoc]{tocbibind} % for references/bibliography in the toc
-
-%-------------------- end preamble ----------------------
-
-\begin{document}
-
-% matching end for #ifdef PREAMBLE
-% #endif
-
-\newcommand{\exercisesection}[1]{\subsection*{#1}}
-
-
-% ------------------- main content ----------------------
-
-
-
-% ----------------- title -------------------------
-
-\thispagestyle{empty}
-
-\begin{center}
-{\LARGE\bf
-\begin{spacing}{1.25}
-Project on Machine Learning
-\end{spacing}
-}
-\end{center}
-
-% ----------------- author(s) -------------------------
-
-\begin{center}
-{\bf \href{{http://www.uio.no/studier/emner/matnat/fys/FYS3155/index-eng.html}}{Data Analysis and Machine Learning FYS-MAT3155/FYS4155}}
-\end{center}
-
- \begin{center}
-% List of all institutions:
-\centerline{{\small Department of Physics, University of Oslo, Norway}}
-\end{center}
-
-% ----------------- end author(s) -------------------------
-
-% --- begin date ---
-\begin{center}
-Fall semester 2017
-\end{center}
-% --- end date ---
-
-\vspace{1cm}
-
-
-\subsection{Using results from Monte Carlo models for machine learning}
-
-\paragraph{Introduction.}
-The aim of this project is to use an already developed Monte Carlo program (either the ising Model or a variational Monte Carlo code) to produce, in case of the \href{{https://github.com/CompPhysics/MachineLearning/tree/master/doc/Programs/IsingModel}}{Ising model}, the energy as function of temperature. For the \href{{https://github.com/CompPhysics/MachineLearning/tree/master/doc/Programs/VMC2Electrons/}}{variational Monte Carlo calculation} of interacting electrons in an oscilaltor trap, the data are represented by the ground state energies as functions of the variational parameters.
-
-In its simplest form
-the energy of the Ising model is expressed as, without an externally applied magnetic field,
-\[
-E=-J\sum_{< kl >}^{N}s_ks_l
-\]
-with
-$s_k=\pm 1$. The quantity $N$ represents the total number of spins and $J$ is a coupling
-constant expressing the strength of the interaction between
-neighboring spins. The symbol $$ indicates that we sum over
-nearest neighbors only. We will assume that we have a ferromagnetic
-ordering, viz $J> 0$. We will use periodic boundary conditions and
-the Metropolis algorithm only. Alternatively we can use the supplied variational Monte Carlo program which solves Schroedinger's equation for two interacting electrons in a harmonic oscillator trap. Both codes can be found at the webpage of the \href{{https://github.com/CompPhysics/MachineLearning/tree/master/doc/Programs/}}{course under programs}.
-
-\paragraph{Part a): Producing the data.}
-If we opt for the Ising model code, we need to generate for every span over the lattice and output of the total energy and the magnetic moment (if we wish to study this quantity as well), that is compute and write to file as function of temperature
-$\langle E\rangle$ and $\langle \vert M\vert \rangle$.
-We will use a fixed lattice siste of $L\times L = 40 \times 40$. Make sure the calculations have been equilibrated and compute these expectation values for
-for $T\in [1.0,3.0]$ with a step in
-temperature $\Delta T=0.1$ or smaller. This are the data you will use to estimate the standard deviation in the next part of the project.
-You should keep a file for the different temperatures as you will need to compute the standard deviation for each temperature.
-
-
-Alternatively you can run the variational Monte Carlo program for two interacting electrons confined to move in a harmonic oscillator trap. Here you will produce the variational expectaction value of the energy as a function of the variational parameters $\alpha$ and $\beta$. You could single out only a set of such parameters after you have found the minimum energy.
-
-These two cases form then our training data which we will use for our estimates of the standard deviation.
-
-
-\paragraph{Part b): Estimating the standard deviation.}
-This part deals with widely used resampling methods to find the standard deviation on our data sets. The methods we can use are the Blocking method, the Bootstrap method and the Jackknife method. The latter two are examples of so-called resampling methods. \href{{https://github.com/CompPhysics/MachineLearning/tree/master/doc/Programs/Sampling}}{The functions provided under the program folder} encode all three methods. The typical situation is to generate data using a C++ or Fortan code, as in part a) above. These data are then used by a post-analysis program in order to perform a statistical analysis.
-
-In case of the Ising model, the final product will be a table of the expectation values of the energy (or the magnetic moment) as function of temperature and with a proper standard deviation. These data will then enter our analysis in parts c-e). We will call these data the training data.
-
-If you opt for the variational Monte Carlo program, you need to make a table of the ground state energy as function of various variational parameters (you can limit yourself to $\alpha$ only) with proper standard deviation estimations.
-
-Your task here is thus to use the data from part a) and generate the standard deviation using one (or more) of the methods discussed above. In general, the Bootstrap method is the most widely used one.
-
-These data will then be used in our next step, where we will try to obtain a functional approximation to our data sets.
-
-\paragraph{Part c): Fitting the data using regression analysis and other methods.}
-With the data from part b), with a proper determination of the standard deviation, the task here is to find a fit to the data using regression methods. Here we will focus on the least square method with a straight line fit first and thereafter using singular value decomposition to fit a function that reproduces the data.
-
-Since obtaining these data points may not be trivial, we want to use these data to fit a function which can allow us to make predictions for values of $y$ which are not in the present set. The perhaps simplest approach is to assume we can parametrize our function in terms of a polynomial of degree $n-1$ with $n$ points, that is
-\[
-y=y(x) \rightarrow y(x_i)=\tilde{y}_i+\epsilon_i=\sum_{j=0}^{n-1} \beta_i x_i^j+\epsilon_i,
-\]
-where $\epsilon_i$ is the error in our approximation.
-
-In order to find the optimal parameters $\beta_i$ instead of solving the above linear algebra problem, we define a function which gives a measure of the spread between the values $y_i$ (which represent hopefully the exact values) and the parametrized values $\tilde{y}_i$, namely
-\[
-Q(\hat{\beta})=\sum_{i=0}^{n-1}\left(y_i-\tilde{y}_i\right)^2=\left(\hat{y}-\hat{\tilde{y}}\right)^T\left(\hat{y}-\hat{\tilde{y}}\right),
-\]
-or using the matrix $\hat{X}$ as
-\[
-Q(\hat{\beta})=\left(\hat{y}-\hat{X}\hat{\beta}\right)^T\left(\hat{y}-\hat{X}\hat{\beta}\right).
-\]
-
-The function
-\[
-Q(\hat{\beta})=\left(\hat{y}-\hat{X}\hat{\beta}\right)^T\left(\hat{y}-\hat{X}\hat{\beta}\right),
-\]
-can be linked to the variance of the quantity $y_i$ if we interpret the latter as the mean value of for example a numerical experiment. When linking below with the maximum likelihood approach below, we will indeed interpret $y_i$ as a mean value
-\[
-y_{i}=\langle y_i \rangle = \beta_0x_{i,0}+\beta_1x_{i,1}+\beta_2x_{i,2}+\dots+\beta_{n-1}x_{i,n-1}+\epsilon_i,
-\]
-where $\langle y_i \rangle$ is the mean value. Keep in mind also that till now we have treated $y_i$ as the exact value. Normally, the response (dependent or outcome) variable $y_i$ the outcome of a numerical experiment or another type of experiment and is thus only an approximation to the true value. It is then always accompanied by an error estimate, often limited to a statistical error estimate given by the standard deviation discussed earlier. In the discussion here we will treat $y_i$ as our exact value for the response variable.
-
-Introducing the standard deviation $\sigma_i$ for each measurement $y_i$, we define now the $\chi^2$ function as
-\[
-\chi^2(\hat{\beta})=\sum_{i=0}^{n-1}\frac{\left(y_i-\tilde{y}_i\right)^2}{\sigma_i^2}=\left(\hat{y}-\hat{\tilde{y}}\right)^T\frac{1}{\hat{\Sigma^2}}\left(\hat{y}-\hat{\tilde{y}}\right),
-\]
-where the matrix $\hat{\Sigma}$ is a diagonal matrix with $\sigma_i$ as matrix elements.
-
-
-
-
-In order to find the parameters $\beta_i$ we will then minimize the spread of $\chi^2(\hat{\beta})$ by requiring
-\[
-\frac{\partial \chi^2(\hat{\beta})}{\partial \beta_j} = \frac{\partial }{\partial \beta_j}\left[ \sum_{i=0}^{n-1}\left(\frac{y_i-\beta_0x_{i,0}-\beta_1x_{i,1}-\beta_2x_{i,2}-\dots-\beta_{n-1}x_{i,n-1}}{\sigma_i}\right)^2\right]=0,
-\]
-which results in
-\[
-\frac{\partial \chi^2(\hat{\beta})}{\partial \beta_j} = -2\left[ \sum_{i=0}^{n-1}\frac{x_{ij}}{\sigma_i}\left(\frac{y_i-\beta_0x_{i,0}-\beta_1x_{i,1}-\beta_2x_{i,2}-\dots-\beta_{n-1}x_{i,n-1}}{\sigma_i}\right)\right]=0,
-\]
-or in a matrix-vector form as
-\[
-\frac{\partial \chi^2(\hat{\beta})}{\partial \hat{\beta}} = 0 = \hat{A}^T\left( \hat{b}-\hat{A}\hat{\beta}\right).
-\]
-where we have defined the matrix $\hat{A} =\hat{X}/\hat{\Sigma}$ with matrix elements $a_{ij} = x_{ij}/\sigma_i$ and the vector $\hat{b}$ with elements $b_i = y_i/\sigma_i$.
-
-
-We can rewrite
-\[
-\frac{\partial \chi^2(\hat{\beta})}{\partial \hat{\beta}} = 0 = \hat{A}^T\left( \hat{b}-\hat{A}\hat{\beta}\right),
-\]
-as
-\[
-\hat{A}^T\hat{b} = \hat{A}^T\hat{A}\hat{\beta},
-\]
-and if the matrix $\hat{A}^T\hat{A}$ is invertible we have the solution
-\[
-\hat{\beta} =\left(\hat{A}^T\hat{A}\right)^{-1}\hat{A}^T\hat{b}.
-\]
-If we then introduce the matrix
-\[
-\hat{H} = \hat{A}^T\hat{A},
-\]
-we have then the following expression for the parameters $\beta_j$ (the matrix elements of $\hat{H}$ are $h_{ij}$)
-\[
-\beta_j = \sum_{k=0}^{p-1}h_{jk}\sum_{i=0}^{n-1}\frac{y_i}{\sigma_i}\frac{x_{ik}}{\sigma_i} = \sum_{k=0}^{p-1}h_{jk}\sum_{i=0}^{n-1}b_ia_{ik}
-\]
-We state without proof the expression for the uncertainty in the parameters $\beta_j$ as
-\[
-\sigma^2(\beta_j) = \sum_{i=0}^{n-1}\sigma_i^2\left( \frac{\partial \beta_j}{\partial y_i}\right)^2,
-\]
-resulting in
-\[
-\sigma^2(\beta_j) = \left(\sum_{k=0}^{p-1}h_{jk}\sum_{i=0}^{n-1}a_{ik}\right)\left(\sum_{l=0}^{p-1}h_{jl}\sum_{m=0}^{n-1}a_{ml}\right) = h_{jj}!
-\]
-
-
-The first step here is to approximate the function $y$ with a first-order polynomial, that is we write
-\[
-y=y(x) \rightarrow y(x_i) \approx \beta_0+\beta_1 x_i.
-\]
-By computing the derivatives of $\chi^2$ with respect to $\beta_0$ and $\beta_1$ show that these are given by
-\[
-\frac{\partial \chi^2(\hat{\beta})}{\partial \beta_0} = -2\left[ \sum_{i=0}^{1}\left(\frac{y_i-\beta_0-\beta_1x_{i}}{\sigma_i^2}\right)\right]=0,
-\]
-and
-\[
-\frac{\partial \chi^2(\hat{\beta})}{\partial \beta_0} = -2\left[ \sum_{i=0}^{1}x_i\left(\frac{y_i-\beta_0-\beta_1x_{i}}{\sigma_i^2}\right)\right]=0.
-\]
-Define then
-\[
-\gamma = \sum_{i=0}^{1}\frac{1}{\sigma_i^2},
-\]
-
-\[
-\gamma_x = \sum_{i=0}^{1}\frac{x_{i}}{\sigma_i^2},
-\]
-\[
-\gamma_y = \sum_{i=0}^{1}\left(\frac{y_i}{\sigma_i^2}\right),
-\]
-\[
-\gamma_{xx} = \sum_{i=0}^{1}\frac{x_ix_{i}}{\sigma_i^2},
-\]
-\[
-\gamma_{xy} = \sum_{i=0}^{1}\frac{y_ix_{i}}{\sigma_i^2},
-\]
-and show that
-\[
-\beta_0 = \frac{\gamma_{xx}\gamma_y-\gamma_x\gamma_y}{\gamma\gamma_{xx}-\gamma_x^2},
-\]
-\[
-\beta_1 = \frac{\gamma_{xy}\gamma-\gamma_x\gamma_y}{\gamma\gamma_{xx}-\gamma_x^2}.
-\]
-Find these parameters for your data set and compare the fitted function $y=y(x) \rightarrow y(x_i) \approx \beta_0+\beta_1 x_i$ and find the error in the parameters $\beta_0$ and $\beta_1$ as well.
-How does your parametrization compare with the input data $y_i$? Give a critical discussion of the method.
-
-Develop now a program where you extract the parameters $\beta_i$ and their corresponding errors using the general Least Square method (LSM), with polynomials of degree $p=2$ and $p=3$. Which polyomials gives the best fit to your data? Here you should prepare an input file with the number of data points in your training data set, the corresponding standard deviation $\sigma_i$ as well as the degree $p$ of the polynomial you want to fit.
-
-The LSM suffers often from both being underdetermined and overdetermined in the unknown coefficients $\beta_i$. A better approach is to use the Singular Value Decomposition (SVD) method discussed in the lecture notes on regression analysis. Here you should use the \href{{https://github.com/CompPhysics/MachineLearning/tree/master/doc/Programs/SVD}}{SVD programs} and use the same input file you used for the general LSM. In the main functions you should specifiy which function you wish to approximate the date with. The example programs discussed above show you examples on how to use these codes as well as how to write the functional form of the polynomial or fitting function.
-
-
-\paragraph{Part d): Introducing Bayesian statistics.}
-More text to come
-
-\paragraph{Part e): Studying the Ising model or the VMC results with Neural networks.}
-More text to come
-
-\subsection{Background literature}
-
-If you wish to read more about the Ising model and statistical physics here are three suggestions.
-
-\begin{itemize}
- \item \href{{http://www.worldscientific.com/worldscibooks/10.1142/5660}}{M. Plischke and B. Bergersen}, \emph{Equilibrium Statistical Physics}, World Scientific, see chapters 5 and 6.
-
- \item \href{{http://www.cambridge.org/no/academic/subjects/physics/computational-science-and-modelling/guide-monte-carlo-simulations-statistical-physics-4th-edition?format=HB}}{D. P. Landau and K. Binder}, \emph{A Guide to Monte Carlo Simulations in Statistical Physics}, Cambridge, see chapters 2,3 and 4.
-
- \item \href{{https://global.oup.com/academic/product/monte-carlo-methods-in-statistical-physics-9780198517979?cc=no&lang=en&}}{M. E. J. Newman and T. Barkema}, \emph{Monte Carlo Methods in Statistical Physics}, Oxford, see chapters 3 and 4.
-\end{itemize}
-
-\noindent
-\subsection{Introduction to numerical projects}
-
-Here follows a brief recipe and recommendation on how to write a report for each
-project.
-
-\begin{itemize}
- \item Give a short description of the nature of the problem and the eventual numerical methods you have used.
-
- \item Describe the algorithm you have used and/or developed. Here you may find it convenient to use pseudocoding. In many cases you can describe the algorithm in the program itself.
-
- \item Include the source code of your program. Comment your program properly.
-
- \item If possible, try to find analytic solutions, or known limits in order to test your program when developing the code.
-
- \item Include your results either in figure form or in a table. Remember to label your results. All tables and figures should have relevant captions and labels on the axes.
-
- \item Try to evaluate the reliabilty and numerical stability/precision of your results. If possible, include a qualitative and/or quantitative discussion of the numerical stability, eventual loss of precision etc.
-
- \item Try to give an interpretation of you results in your answers to the problems.
-
- \item Critique: if possible include your comments and reflections about the exercise, whether you felt you learnt something, ideas for improvements and other thoughts you've made when solving the exercise. We wish to keep this course at the interactive level and your comments can help us improve it.
-
- \item Try to establish a practice where you log your work at the computerlab. You may find such a logbook very handy at later stages in your work, especially when you don't properly remember what a previous test version of your program did. Here you could also record the time spent on solving the exercise, various algorithms you may have tested or other topics which you feel worthy of mentioning.
-\end{itemize}
-
-\noindent
-
-% ------------------- end of main content ---------------
-
-% #ifdef PREAMBLE
-\end{document}
-% #endif
-
diff --git a/doc/Projects/2017/Project/pdf/Project.pdf b/doc/Projects/2017/Project/pdf/Project.pdf
deleted file mode 100644
index 41570e4bb..000000000
Binary files a/doc/Projects/2017/Project/pdf/Project.pdf and /dev/null differ
diff --git a/doc/Projects/2017/Project/pdf/Project.tex b/doc/Projects/2017/Project/pdf/Project.tex
deleted file mode 100644
index 00b1c57da..000000000
--- a/doc/Projects/2017/Project/pdf/Project.tex
+++ /dev/null
@@ -1,351 +0,0 @@
-%%
-%% Automatically generated file from DocOnce source
-%% (https://github.com/hplgit/doconce/)
-%%
-%%
-
-
-%-------------------- begin preamble ----------------------
-
-\documentclass[%
-oneside, % oneside: electronic viewing, twoside: printing
-final, % draft: marks overfull hboxes, figures with paths
-10pt]{article}
-
-\listfiles % print all files needed to compile this document
-
-\usepackage{relsize,makeidx,color,setspace,amsmath,amsfonts,amssymb}
-\usepackage[table]{xcolor}
-\usepackage{bm,ltablex,microtype}
-
-\usepackage[pdftex]{graphicx}
-
-\usepackage[T1]{fontenc}
-%\usepackage[latin1]{inputenc}
-\usepackage{ucs}
-\usepackage[utf8x]{inputenc}
-
-\usepackage{lmodern} % Latin Modern fonts derived from Computer Modern
-
-% Hyperlinks in PDF:
-\definecolor{linkcolor}{rgb}{0,0,0.4}
-\usepackage{hyperref}
-\hypersetup{
- breaklinks=true,
- colorlinks=true,
- linkcolor=linkcolor,
- urlcolor=linkcolor,
- citecolor=black,
- filecolor=black,
- %filecolor=blue,
- pdfmenubar=true,
- pdftoolbar=true,
- bookmarksdepth=3 % Uncomment (and tweak) for PDF bookmarks with more levels than the TOC
- }
-%\hyperbaseurl{} % hyperlinks are relative to this root
-
-\setcounter{tocdepth}{2} % levels in table of contents
-
-% --- fancyhdr package for fancy headers ---
-\usepackage{fancyhdr}
-\fancyhf{} % sets both header and footer to nothing
-\renewcommand{\headrulewidth}{0pt}
-\fancyfoot[LE,RO]{\thepage}
-% Ensure copyright on titlepage (article style) and chapter pages (book style)
-\fancypagestyle{plain}{
- \fancyhf{}
- \fancyfoot[C]{{\footnotesize \copyright\ 1999-2017, "Data Analysis and Machine Learning FYS-MAT3155/FYS4155":"http://www.uio.no/studier/emner/matnat/fys/FYS3155/index-eng.html". Released under CC Attribution-NonCommercial 4.0 license}}
-% \renewcommand{\footrulewidth}{0mm}
- \renewcommand{\headrulewidth}{0mm}
-}
-% Ensure copyright on titlepages with \thispagestyle{empty}
-\fancypagestyle{empty}{
- \fancyhf{}
- \fancyfoot[C]{{\footnotesize \copyright\ 1999-2017, "Data Analysis and Machine Learning FYS-MAT3155/FYS4155":"http://www.uio.no/studier/emner/matnat/fys/FYS3155/index-eng.html". Released under CC Attribution-NonCommercial 4.0 license}}
- \renewcommand{\footrulewidth}{0mm}
- \renewcommand{\headrulewidth}{0mm}
-}
-
-\pagestyle{fancy}
-
-
-% prevent orhpans and widows
-\clubpenalty = 10000
-\widowpenalty = 10000
-
-% --- end of standard preamble for documents ---
-
-
-% insert custom LaTeX commands...
-
-\raggedbottom
-\makeindex
-\usepackage[totoc]{idxlayout} % for index in the toc
-\usepackage[nottoc]{tocbibind} % for references/bibliography in the toc
-
-%-------------------- end preamble ----------------------
-
-\begin{document}
-
-% matching end for #ifdef PREAMBLE
-
-\newcommand{\exercisesection}[1]{\subsection*{#1}}
-
-
-% ------------------- main content ----------------------
-
-
-
-% ----------------- title -------------------------
-
-\thispagestyle{empty}
-
-\begin{center}
-{\LARGE\bf
-\begin{spacing}{1.25}
-Project on Machine Learning
-\end{spacing}
-}
-\end{center}
-
-% ----------------- author(s) -------------------------
-
-\begin{center}
-{\bf \href{{http://www.uio.no/studier/emner/matnat/fys/FYS3155/index-eng.html}}{Data Analysis and Machine Learning FYS-MAT3155/FYS4155}}
-\end{center}
-
- \begin{center}
-% List of all institutions:
-\centerline{{\small Department of Physics, University of Oslo, Norway}}
-\end{center}
-
-% ----------------- end author(s) -------------------------
-
-% --- begin date ---
-\begin{center}
-Fall semester 2017
-\end{center}
-% --- end date ---
-
-\vspace{1cm}
-
-
-\subsection*{Using results from Monte Carlo models for machine learning}
-
-\paragraph{Introduction.}
-The aim of this project is to use an already developed Monte Carlo program (either the ising Model or a variational Monte Carlo code) to produce, in case of the \href{{https://github.com/CompPhysics/MachineLearning/tree/master/doc/Programs/IsingModel}}{Ising model}, the energy as function of temperature. For the \href{{https://github.com/CompPhysics/MachineLearning/tree/master/doc/Programs/VMC2Electrons/}}{variational Monte Carlo calculation} of interacting electrons in an oscilaltor trap, the data are represented by the ground state energies as functions of the variational parameters.
-
-In its simplest form
-the energy of the Ising model is expressed as, without an externally applied magnetic field,
-\[
-E=-J\sum_{< kl >}^{N}s_ks_l
-\]
-with
-$s_k=\pm 1$. The quantity $N$ represents the total number of spins and $J$ is a coupling
-constant expressing the strength of the interaction between
-neighboring spins. The symbol $$ indicates that we sum over
-nearest neighbors only. We will assume that we have a ferromagnetic
-ordering, viz $J> 0$. We will use periodic boundary conditions and
-the Metropolis algorithm only. Alternatively we can use the supplied variational Monte Carlo program which solves Schroedinger's equation for two interacting electrons in a harmonic oscillator trap. Both codes can be found at the webpage of the \href{{https://github.com/CompPhysics/MachineLearning/tree/master/doc/Programs/}}{course under programs}.
-
-\paragraph{Part a): Producing the data.}
-If we opt for the Ising model code, we need to generate for every span over the lattice and output of the total energy and the magnetic moment (if we wish to study this quantity as well), that is compute and write to file as function of temperature
-$\langle E\rangle$ and $\langle \vert M\vert \rangle$.
-We will use a fixed lattice siste of $L\times L = 40 \times 40$. Make sure the calculations have been equilibrated and compute these expectation values for
-for $T\in [1.0,3.0]$ with a step in
-temperature $\Delta T=0.1$ or smaller. This are the data you will use to estimate the standard deviation in the next part of the project.
-You should keep a file for the different temperatures as you will need to compute the standard deviation for each temperature.
-
-
-Alternatively you can run the variational Monte Carlo program for two interacting electrons confined to move in a harmonic oscillator trap. Here you will produce the variational expectaction value of the energy as a function of the variational parameters $\alpha$ and $\beta$. You could single out only a set of such parameters after you have found the minimum energy.
-
-These two cases form then our training data which we will use for our estimates of the standard deviation.
-
-
-\paragraph{Part b): Estimating the standard deviation.}
-This part deals with widely used resampling methods to find the standard deviation on our data sets. The methods we can use are the Blocking method, the Bootstrap method and the Jackknife method. The latter two are examples of so-called resampling methods. \href{{https://github.com/CompPhysics/MachineLearning/tree/master/doc/Programs/Sampling}}{The functions provided under the program folder} encode all three methods. The typical situation is to generate data using a C++ or Fortan code, as in part a) above. These data are then used by a post-analysis program in order to perform a statistical analysis.
-
-In case of the Ising model, the final product will be a table of the expectation values of the energy (or the magnetic moment) as function of temperature and with a proper standard deviation. These data will then enter our analysis in parts c-e). We will call these data the training data.
-
-If you opt for the variational Monte Carlo program, you need to make a table of the ground state energy as function of various variational parameters (you can limit yourself to $\alpha$ only) with proper standard deviation estimations.
-
-Your task here is thus to use the data from part a) and generate the standard deviation using one (or more) of the methods discussed above. In general, the Bootstrap method is the most widely used one.
-
-These data will then be used in our next step, where we will try to obtain a functional approximation to our data sets.
-
-\paragraph{Part c): Fitting the data using regression analysis and other methods.}
-With the data from part b), with a proper determination of the standard deviation, the task here is to find a fit to the data using regression methods. Here we will focus on the least square method with a straight line fit first and thereafter using singular value decomposition to fit a function that reproduces the data.
-
-Since obtaining these data points may not be trivial, we want to use these data to fit a function which can allow us to make predictions for values of $y$ which are not in the present set. The perhaps simplest approach is to assume we can parametrize our function in terms of a polynomial of degree $n-1$ with $n$ points, that is
-\[
-y=y(x) \rightarrow y(x_i)=\tilde{y}_i+\epsilon_i=\sum_{j=0}^{n-1} \beta_i x_i^j+\epsilon_i,
-\]
-where $\epsilon_i$ is the error in our approximation.
-
-In order to find the optimal parameters $\beta_i$ instead of solving the above linear algebra problem, we define a function which gives a measure of the spread between the values $y_i$ (which represent hopefully the exact values) and the parametrized values $\tilde{y}_i$, namely
-\[
-Q(\hat{\beta})=\sum_{i=0}^{n-1}\left(y_i-\tilde{y}_i\right)^2=\left(\hat{y}-\hat{\tilde{y}}\right)^T\left(\hat{y}-\hat{\tilde{y}}\right),
-\]
-or using the matrix $\hat{X}$ as
-\[
-Q(\hat{\beta})=\left(\hat{y}-\hat{X}\hat{\beta}\right)^T\left(\hat{y}-\hat{X}\hat{\beta}\right).
-\]
-
-The function
-\[
-Q(\hat{\beta})=\left(\hat{y}-\hat{X}\hat{\beta}\right)^T\left(\hat{y}-\hat{X}\hat{\beta}\right),
-\]
-can be linked to the variance of the quantity $y_i$ if we interpret the latter as the mean value of for example a numerical experiment. When linking below with the maximum likelihood approach below, we will indeed interpret $y_i$ as a mean value
-\[
-y_{i}=\langle y_i \rangle = \beta_0x_{i,0}+\beta_1x_{i,1}+\beta_2x_{i,2}+\dots+\beta_{n-1}x_{i,n-1}+\epsilon_i,
-\]
-where $\langle y_i \rangle$ is the mean value. Keep in mind also that till now we have treated $y_i$ as the exact value. Normally, the response (dependent or outcome) variable $y_i$ the outcome of a numerical experiment or another type of experiment and is thus only an approximation to the true value. It is then always accompanied by an error estimate, often limited to a statistical error estimate given by the standard deviation discussed earlier. In the discussion here we will treat $y_i$ as our exact value for the response variable.
-
-Introducing the standard deviation $\sigma_i$ for each measurement $y_i$, we define now the $\chi^2$ function as
-\[
-\chi^2(\hat{\beta})=\sum_{i=0}^{n-1}\frac{\left(y_i-\tilde{y}_i\right)^2}{\sigma_i^2}=\left(\hat{y}-\hat{\tilde{y}}\right)^T\frac{1}{\hat{\Sigma^2}}\left(\hat{y}-\hat{\tilde{y}}\right),
-\]
-where the matrix $\hat{\Sigma}$ is a diagonal matrix with $\sigma_i$ as matrix elements.
-
-
-
-
-In order to find the parameters $\beta_i$ we will then minimize the spread of $\chi^2(\hat{\beta})$ by requiring
-\[
-\frac{\partial \chi^2(\hat{\beta})}{\partial \beta_j} = \frac{\partial }{\partial \beta_j}\left[ \sum_{i=0}^{n-1}\left(\frac{y_i-\beta_0x_{i,0}-\beta_1x_{i,1}-\beta_2x_{i,2}-\dots-\beta_{n-1}x_{i,n-1}}{\sigma_i}\right)^2\right]=0,
-\]
-which results in
-\[
-\frac{\partial \chi^2(\hat{\beta})}{\partial \beta_j} = -2\left[ \sum_{i=0}^{n-1}\frac{x_{ij}}{\sigma_i}\left(\frac{y_i-\beta_0x_{i,0}-\beta_1x_{i,1}-\beta_2x_{i,2}-\dots-\beta_{n-1}x_{i,n-1}}{\sigma_i}\right)\right]=0,
-\]
-or in a matrix-vector form as
-\[
-\frac{\partial \chi^2(\hat{\beta})}{\partial \hat{\beta}} = 0 = \hat{A}^T\left( \hat{b}-\hat{A}\hat{\beta}\right).
-\]
-where we have defined the matrix $\hat{A} =\hat{X}/\hat{\Sigma}$ with matrix elements $a_{ij} = x_{ij}/\sigma_i$ and the vector $\hat{b}$ with elements $b_i = y_i/\sigma_i$.
-
-
-We can rewrite
-\[
-\frac{\partial \chi^2(\hat{\beta})}{\partial \hat{\beta}} = 0 = \hat{A}^T\left( \hat{b}-\hat{A}\hat{\beta}\right),
-\]
-as
-\[
-\hat{A}^T\hat{b} = \hat{A}^T\hat{A}\hat{\beta},
-\]
-and if the matrix $\hat{A}^T\hat{A}$ is invertible we have the solution
-\[
-\hat{\beta} =\left(\hat{A}^T\hat{A}\right)^{-1}\hat{A}^T\hat{b}.
-\]
-If we then introduce the matrix
-\[
-\hat{H} = \hat{A}^T\hat{A},
-\]
-we have then the following expression for the parameters $\beta_j$ (the matrix elements of $\hat{H}$ are $h_{ij}$)
-\[
-\beta_j = \sum_{k=0}^{p-1}h_{jk}\sum_{i=0}^{n-1}\frac{y_i}{\sigma_i}\frac{x_{ik}}{\sigma_i} = \sum_{k=0}^{p-1}h_{jk}\sum_{i=0}^{n-1}b_ia_{ik}
-\]
-We state without proof the expression for the uncertainty in the parameters $\beta_j$ as
-\[
-\sigma^2(\beta_j) = \sum_{i=0}^{n-1}\sigma_i^2\left( \frac{\partial \beta_j}{\partial y_i}\right)^2,
-\]
-resulting in
-\[
-\sigma^2(\beta_j) = \left(\sum_{k=0}^{p-1}h_{jk}\sum_{i=0}^{n-1}a_{ik}\right)\left(\sum_{l=0}^{p-1}h_{jl}\sum_{m=0}^{n-1}a_{ml}\right) = h_{jj}!
-\]
-
-
-The first step here is to approximate the function $y$ with a first-order polynomial, that is we write
-\[
-y=y(x) \rightarrow y(x_i) \approx \beta_0+\beta_1 x_i.
-\]
-By computing the derivatives of $\chi^2$ with respect to $\beta_0$ and $\beta_1$ show that these are given by
-\[
-\frac{\partial \chi^2(\hat{\beta})}{\partial \beta_0} = -2\left[ \sum_{i=0}^{1}\left(\frac{y_i-\beta_0-\beta_1x_{i}}{\sigma_i^2}\right)\right]=0,
-\]
-and
-\[
-\frac{\partial \chi^2(\hat{\beta})}{\partial \beta_0} = -2\left[ \sum_{i=0}^{1}x_i\left(\frac{y_i-\beta_0-\beta_1x_{i}}{\sigma_i^2}\right)\right]=0.
-\]
-Define then
-\[
-\gamma = \sum_{i=0}^{1}\frac{1}{\sigma_i^2},
-\]
-
-\[
-\gamma_x = \sum_{i=0}^{1}\frac{x_{i}}{\sigma_i^2},
-\]
-\[
-\gamma_y = \sum_{i=0}^{1}\left(\frac{y_i}{\sigma_i^2}\right),
-\]
-\[
-\gamma_{xx} = \sum_{i=0}^{1}\frac{x_ix_{i}}{\sigma_i^2},
-\]
-\[
-\gamma_{xy} = \sum_{i=0}^{1}\frac{y_ix_{i}}{\sigma_i^2},
-\]
-and show that
-\[
-\beta_0 = \frac{\gamma_{xx}\gamma_y-\gamma_x\gamma_y}{\gamma\gamma_{xx}-\gamma_x^2},
-\]
-\[
-\beta_1 = \frac{\gamma_{xy}\gamma-\gamma_x\gamma_y}{\gamma\gamma_{xx}-\gamma_x^2}.
-\]
-Find these parameters for your data set and compare the fitted function $y=y(x) \rightarrow y(x_i) \approx \beta_0+\beta_1 x_i$ and find the error in the parameters $\beta_0$ and $\beta_1$ as well.
-How does your parametrization compare with the input data $y_i$? Give a critical discussion of the method.
-
-Develop now a program where you extract the parameters $\beta_i$ and their corresponding errors using the general Least Square method (LSM), with polynomials of degree $p=2$ and $p=3$. Which polyomials gives the best fit to your data? Here you should prepare an input file with the number of data points in your training data set, the corresponding standard deviation $\sigma_i$ as well as the degree $p$ of the polynomial you want to fit.
-
-The LSM suffers often from both being underdetermined and overdetermined in the unknown coefficients $\beta_i$. A better approach is to use the Singular Value Decomposition (SVD) method discussed in the lecture notes on regression analysis. Here you should use the \href{{https://github.com/CompPhysics/MachineLearning/tree/master/doc/Programs/SVD}}{SVD programs} and use the same input file you used for the general LSM. In the main functions you should specifiy which function you wish to approximate the date with. The example programs discussed above show you examples on how to use these codes as well as how to write the functional form of the polynomial or fitting function.
-
-
-\paragraph{Part d): Introducing Bayesian statistics.}
-More text to come
-
-\paragraph{Part e): Studying the Ising model or the VMC results with Neural networks.}
-More text to come
-
-\subsection*{Background literature}
-
-If you wish to read more about the Ising model and statistical physics here are three suggestions.
-
-\begin{itemize}
- \item \href{{http://www.worldscientific.com/worldscibooks/10.1142/5660}}{M. Plischke and B. Bergersen}, \emph{Equilibrium Statistical Physics}, World Scientific, see chapters 5 and 6.
-
- \item \href{{http://www.cambridge.org/no/academic/subjects/physics/computational-science-and-modelling/guide-monte-carlo-simulations-statistical-physics-4th-edition?format=HB}}{D. P. Landau and K. Binder}, \emph{A Guide to Monte Carlo Simulations in Statistical Physics}, Cambridge, see chapters 2,3 and 4.
-
- \item \href{{https://global.oup.com/academic/product/monte-carlo-methods-in-statistical-physics-9780198517979?cc=no&lang=en&}}{M. E. J. Newman and T. Barkema}, \emph{Monte Carlo Methods in Statistical Physics}, Oxford, see chapters 3 and 4.
-\end{itemize}
-
-\noindent
-\subsection*{Introduction to numerical projects}
-
-Here follows a brief recipe and recommendation on how to write a report for each
-project.
-
-\begin{itemize}
- \item Give a short description of the nature of the problem and the eventual numerical methods you have used.
-
- \item Describe the algorithm you have used and/or developed. Here you may find it convenient to use pseudocoding. In many cases you can describe the algorithm in the program itself.
-
- \item Include the source code of your program. Comment your program properly.
-
- \item If possible, try to find analytic solutions, or known limits in order to test your program when developing the code.
-
- \item Include your results either in figure form or in a table. Remember to label your results. All tables and figures should have relevant captions and labels on the axes.
-
- \item Try to evaluate the reliabilty and numerical stability/precision of your results. If possible, include a qualitative and/or quantitative discussion of the numerical stability, eventual loss of precision etc.
-
- \item Try to give an interpretation of you results in your answers to the problems.
-
- \item Critique: if possible include your comments and reflections about the exercise, whether you felt you learnt something, ideas for improvements and other thoughts you've made when solving the exercise. We wish to keep this course at the interactive level and your comments can help us improve it.
-
- \item Try to establish a practice where you log your work at the computerlab. You may find such a logbook very handy at later stages in your work, especially when you don't properly remember what a previous test version of your program did. Here you could also record the time spent on solving the exercise, various algorithms you may have tested or other topics which you feel worthy of mentioning.
-\end{itemize}
-
-\noindent
-
-% ------------------- end of main content ---------------
-
-\end{document}
-
diff --git a/doc/pub/Statistics/html/._Statistics-bs000.html b/doc/pub/Statistics/html/._Statistics-bs000.html
deleted file mode 100644
index 967e7f14f..000000000
--- a/doc/pub/Statistics/html/._Statistics-bs000.html
+++ /dev/null
@@ -1,451 +0,0 @@
-
-
-
-
-
-
-
-Data Analysis and Machine Learning: Elements of Probability Theory
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-Consider the following simple example, namely the tossing of a dice, resulting in the following possible values
-$$
-\begin{equation*}
-\{2,3,4,5,6,7,8,9,10,11,12\}.
-\end{equation*}
-$$
-
-These values are called the domain.
-To this domain we have the corresponding probabilities
-$$
-\begin{equation*}
-\{1/36,2/36/3/36,4/36,5/36,6/36,5/36,4/36,3/36,2/36,1/36\}.
-\end{equation*}
-$$
-
-The numbers in the domain are the outcomes of the physical process tossing the dice.
-We cannot tell beforehand whether the outcome is 3 or 5 or any other number in this domain.
-This defines the randomness of the outcome, or unexpectedness or any other synonimous word which
-encompasses the uncertitude of the final outcome.
-
-
-The only thing we can tell beforehand
-is that say the outcome 2 has a certain probability.
-If our favorite hobby is to spend an hour every evening throwing dice and
-registering the sequence of outcomes, we will note that the numbers in the above domain
-$$
-\begin{equation*}
-\{2,3,4,5,6,7,8,9,10,11,12\},
-\end{equation*}
-$$
-
-appear in a random order. After 11 throws the results may look like
-
-$$
-\begin{equation*}
-\{10,8,6,3,6,9,11,8,12,4,5\}.
-\end{equation*}
-$$
-
-Random variables are characterized by a domain which contains all possible values that the random value may take. This domain has a corresponding PDF.
-
Stochastic variables and the main concepts, the discrete case
-
-
-
-There are two main concepts associated with a stochastic variable. The
-domain is the set \( \mathbb D = \{x\} \) of all accessible values
-the variable can assume, so that \( X \in \mathbb D \). An example of a
-discrete domain is the set of six different numbers that we may get by
-throwing of a dice, \( x\in\{1,\,2,\,3,\,4,\,5,\,6\} \).
-
-
-The probability distribution function (PDF) is a function
-\( p(x) \) on the domain which, in the discrete case, gives us the
-probability or relative frequency with which these values of \( X \)
-occur
-$$
-\begin{equation*}
-p(x) = \mathrm{Prob}(X=x).
-\end{equation*}
-$$
-
Stochastic variables and the main concepts, the continuous case
-
-
-
-In the continuous case, the PDF does not directly depict the
-actual probability. Instead we define the probability for the
-stochastic variable to assume any value on an infinitesimal interval
-around \( x \) to be \( p(x)dx \). The continuous function \( p(x) \) then gives us
-the density of the probability rather than the probability
-itself. The probability for a stochastic variable to assume any value
-on a non-infinitesimal interval \( [a,\,b] \) is then just the integral
-
-$$
-\begin{equation*}
-\mathrm{Prob}(a\leq X\leq b) = \int_a^b p(x)dx.
-\end{equation*}
-$$
-
-Qualitatively speaking, a stochastic variable represents the values of
-numbers chosen as if by chance from some specified PDF so that the
-selection of a large set of these numbers reproduces this PDF.
-
-Of interest to us is the cumulative probability
-distribution function (CDF), \( P(x) \), which is just the probability
-for a stochastic variable \( X \) to assume any value less than \( x \)
-$$
-\begin{equation*}
-P(x)=\mathrm{Prob(}X\leq x\mathrm{)} =
-\int_{-\infty}^x p(x^{\prime})dx^{\prime}.
-\end{equation*}
-$$
-
-The relation between a CDF and its corresponding PDF is then
-
-$$
-\begin{equation*}
-p(x) = \frac{d}{dx}P(x).
-\end{equation*}
-$$
-
-There are two properties that all PDFs must satisfy. The first one is
-positivity (assuming that the PDF is normalized)
-
-$$
-\begin{equation*}
-0 \leq p(x) \leq 1.
-\end{equation*}
-$$
-
-Naturally, it would be nonsensical for any of the values of the domain
-to occur with a probability greater than \( 1 \) or less than \( 0 \). Also,
-the PDF must be normalized. That is, all the probabilities must add up
-to unity. The probability of "anything" to happen is always unity. For
-both discrete and continuous PDFs, this condition is
-$$
-\begin{align*}
-\sum_{x_i\in\mathbb D} p(x_i) & = 1,\\
-\int_{x\in\mathbb D} p(x)\,dx & = 1.
-\end{align*}
-$$
-
-The first one
-is the most basic PDF; namely the uniform distribution
-$$
-\begin{equation}
-p(x) = \frac{1}{b-a}\theta(x-a)\theta(b-x),
-\tag{1}
-\end{equation}
-$$
-
-with
-$$
-\begin{equation*}
-\begin{array}{ll}
-\theta(x)=0 & x < 0 \\
-\theta(x)=\frac{1}{b-a} & \in [a,b].
-\end{array}
-\end{equation*}
-$$
-
-The normal distribution with \( b=1 \) and \( a=0 \) is used to generate random numbers.
-
-The second one is the Gaussian Distribution
-$$
-\begin{equation*}
-p(x) = \frac{1}{\sigma\sqrt{2\pi}} \exp{(-\frac{(x-\mu)^2}{2\sigma^2})},
-\end{equation*}
-$$
-
-with mean value \( \mu \) and standard deviation \( \sigma \). If \( \mu=0 \) and \( \sigma=1 \), it is normally called the standard normal distribution
-$$
-\begin{equation*}
-p(x) = \frac{1}{\sqrt{2\pi}} \exp{(-\frac{x^2}{2})},
-\end{equation*}
-$$
-
-
-The following simple Python code plots the above distribution for different values of \( \mu \) and \( \sigma \).
-
-Another important distribution in science is the exponential distribution
-$$
-\begin{equation*}
-p(x) = \alpha\exp{-(\alpha x)}.
-\end{equation*}
-$$
-
-Let \( h(x) \) be an arbitrary continuous function on the domain of the stochastic
-variable \( X \) whose PDF is \( p(x) \). We define the expectation value
-of \( h \) with respect to \( p \) as follows
-
-$$
-\begin{equation}
-\langle h \rangle_X \equiv \int\! h(x)p(x)\,dx
-\tag{2}
-\end{equation}
-$$
-
-Whenever the PDF is known implicitly, like in this case, we will drop
-the index \( X \) for clarity.
-A particularly useful class of special expectation values are the
-moments. The \( n \)-th moment of the PDF \( p \) is defined as
-follows
-$$
-\begin{equation*}
-\langle x^n \rangle \equiv \int\! x^n p(x)\,dx
-\end{equation*}
-$$
-
Stochastic variables and the main concepts, mean values
-
-
-
-The zero-th moment \( \langle 1\rangle \) is just the normalization condition of
-\( p \). The first moment, \( \langle x\rangle \), is called the mean of \( p \)
-and often denoted by the letter \( \mu \)
-$$
-\begin{equation*}
-\langle x\rangle = \mu \equiv \int x p(x)dx,
-\end{equation*}
-$$
-
-for a continuous distribution and
-$$
-\begin{equation*}
-\langle x\rangle = \mu \equiv \frac{1}{N}\sum_{i=1}^N x_i p(x_i),
-\end{equation*}
-$$
-
-for a discrete distribution.
-Qualitatively it represents the centroid or the average value of the
-PDF and is therefore simply called the expectation value of \( p(x) \).
-
Stochastic variables and the main concepts, central moments, the variance
-
-
-
-
-
-A special version of the moments is the set of central moments, the n-th central moment defined as
-$$
-\begin{equation*}
-\langle (x-\langle x\rangle )^n\rangle \equiv \int\! (x-\langle x\rangle)^n p(x)\,dx
-\end{equation*}
-$$
-
-The zero-th and first central moments are both trivial, equal \( 1 \) and
-\( 0 \), respectively. But the second central moment, known as the
-variance of \( p \), is of particular interest. For the stochastic
-variable \( X \), the variance is denoted as \( \sigma^2_X \) or \( \mathrm{Var}(X) \)
-$$
-\begin{align*}
-\sigma^2_X &=\mathrm{Var}(X) = \langle (x-\langle x\rangle)^2\rangle =
-\int (x-\langle x\rangle)^2 p(x)dx\\
-& = \int\left(x^2 - 2 x \langle x\rangle^{2} +\langle x\rangle^2\right)p(x)dx\\
-& = \langle x^2\rangle\rangle - 2 \langle x\rangle\langle x\rangle + \langle x\rangle^2\\
-& = \langle x^2 \rangle - \langle x\rangle^2
-\end{align*}
-$$
-
-The square root of the variance, \( \sigma =\sqrt{\langle (x-\langle x\rangle)^2\rangle} \) is called the
-standard deviation of \( p \). It is the RMS (root-mean-square)
-value of the deviation of the PDF from its mean value, interpreted
-qualitatively as the "spread" of \( p \) around its mean.
-
-The following table collects properties of probability distribution functions.
-In our notation we reserve the label \( p(x) \) for the probability of a certain event,
-while \( P(x) \) is the cumulative probability.
-
-
-With a PDF we can compute expectation values of selected quantities such as
-
-$$
-\begin{equation*}
- \langle x^k\rangle=\frac{1}{N}\sum_{i=1}^{N}x_i^kp(x_i),
-\end{equation*}
-$$
-
-if we have a discrete PDF or
-
-$$
-\begin{equation*}
- \langle x^k\rangle=\int_a^b x^kp(x)dx,
-\end{equation*}
-$$
-
-in the case of a continuous PDF. We have already defined the mean value \( \mu \)
-and the variance \( \sigma^2 \).
-
The three famous Probability Distribution Functions
-
-
-
-
-
-There are at least three PDFs which one may encounter. These are the
-
-
-Uniform distribution
-$$
-\begin{equation*}
-p(x)=\frac{1}{b-a}\Theta(x-a)\Theta(b-x),
-\end{equation*}
-$$
-
-yielding probabilities different from zero in the interval \( [a,b] \).
-
-
-The exponential distribution
-$$
-\begin{equation*}
-p(x)=\alpha \exp{(-\alpha x)},
-\end{equation*}
-$$
-
-yielding probabilities different from zero in the interval \( [0,\infty) \) and with mean value
-$$
-\begin{equation*}
-\mu = \int_0^{\infty}xp(x)dx=\int_0^{\infty}x\alpha \exp{(-\alpha x)}dx=\frac{1}{\alpha},
-\end{equation*}
-$$
-
Probability Distribution Functions, the normal distribution
-
-
-
-Finally, we have the so-called univariate normal distribution, or just the normal distribution
-$$
-\begin{equation*}
-p(x)=\frac{1}{b\sqrt{2\pi}}\exp{\left(-\frac{(x-a)^2}{2b^2}\right)}
-\end{equation*}
-$$
-
-with probabilities different from zero in the interval \( (-\infty,\infty) \).
-The integral \( \int_{-\infty}^{\infty}\exp{\left(-(x^2\right)}dx \) appears in many calculations, its value
-is \( \sqrt{\pi} \), a result we will need when we compute the mean value and the variance.
-The mean value is
-$$
-\begin{equation*}
- \mu = \int_0^{\infty}xp(x)dx=\frac{1}{b\sqrt{2\pi}}\int_{-\infty}^{\infty}x \exp{\left(-\frac{(x-a)^2}{2b^2}\right)}dx,
-\end{equation*}
-$$
-
-which becomes with a suitable change of variables
-$$
-\begin{equation*}
- \mu =\frac{1}{b\sqrt{2\pi}}\int_{-\infty}^{\infty}b\sqrt{2}(a+b\sqrt{2}y)\exp{-y^2}dy=a.
-\end{equation*}
-$$
-
Probability Distribution Functions, the normal distribution
-
-
-
-Similarly, the variance becomes
-$$
-\begin{equation*}
- \sigma^2 = \frac{1}{b\sqrt{2\pi}}\int_{-\infty}^{\infty}(x-\mu)^2 \exp{\left(-\frac{(x-a)^2}{2b^2}\right)}dx,
-\end{equation*}
-$$
-
-and inserting the mean value and performing a variable change we obtain
-
-$$
-\begin{equation*}
- \sigma^2 = \frac{1}{b\sqrt{2\pi}}\int_{-\infty}^{\infty}b\sqrt{2}(b\sqrt{2}y)^2\exp{\left(-y^2\right)}dy=
-\frac{2b^2}{\sqrt{\pi}}\int_{-\infty}^{\infty}y^2\exp{\left(-y^2\right)}dy,
-\end{equation*}
-$$
-
-and performing a final integration by parts we obtain the well-known result \( \sigma^2=b^2 \).
-It is useful to introduce the standard normal distribution as well, defined by \( \mu=a=0 \), viz. a distribution
-centered around zero and with a variance \( \sigma^2=1 \), leading to
-
-$$
-\begin{equation}
- p(x)=\frac{1}{\sqrt{2\pi}}\exp{\left(-\frac{x^2}{2}\right)}.
-\tag{3}
-\end{equation}
-$$
-
Probability Distribution Functions, the cumulative distribution
-
-
-
-
-
-The exponential and uniform distributions have simple cumulative functions,
-whereas the normal distribution does not, being proportional to the so-called
-error function \( erf(x) \), given by
-
-$$
-\begin{equation*}
-P(x) = \frac{1}{\sqrt{2\pi}}\int_{-\infty}^x\exp{\left(-\frac{t^2}{2}\right)}dt,
-\end{equation*}
-$$
-
-which is difficult to evaluate in a quick way.
-
Probability Distribution Functions, other important distribution
-
-
-
-
-
-Some other PDFs which one encounters often in the natural sciences are the binomial distribution
-$$
-\begin{equation*}
- p(x) = \left(\begin{array}{c} n \\ x\end{array}\right)y^x(1-y)^{n-x} \hspace{0.5cm}x=0,1,\dots,n,
-\end{equation*}
-$$
-
-where \( y \) is the probability for a specific event, such as the tossing of a coin or moving left or right
-in case of a random walker. Note that \( x \) is a discrete stochastic variable.
-
-
-The sequence of binomial trials is characterized by the following definitions
-
-
-
Every experiment is thought to consist of \( N \) independent trials.
-
In every independent trial one registers if a specific situation happens or not, such as the jump to the left or right of a random walker.
-
The probability for every outcome in a single trial has the same value, for example the outcome of tossing (either heads or tails) a coin is always \( 1/2 \).
Probability Distribution Functions, the binomial distribution
-
-
-
-
-
-In order to compute the mean and variance we need to recall Newton's binomial
-formula
-$$
-\begin{equation*}
- (a+b)^m=\sum_{n=0}^m \left(\begin{array}{c} m \\ n\end{array}\right)a^nb^{m-n},
-\end{equation*}
-$$
-
-which can be used to show that
-
-$$
-\begin{equation*}
-\sum_{x=0}^n\left(\begin{array}{c} n \\ x\end{array}\right)y^x(1-y)^{n-x} = (y+1-y)^n = 1,
-\end{equation*}
-$$
-
-the PDF is normalized to one.
-The mean value is
-$$
-\begin{equation*}
-\mu = \sum_{x=0}^n x\left(\begin{array}{c} n \\ x\end{array}\right)y^x(1-y)^{n-x} =
-\sum_{x=0}^n x\frac{n!}{x!(n-x)!}y^x(1-y)^{n-x},
-\end{equation*}
-$$
-
-resulting in
-$$
-\begin{equation*}
-\mu =
-\sum_{x=0}^n x\frac{(n-1)!}{(x-1)!(n-1-(x-1))!}y^{x-1}(1-y)^{n-1-(x-1)},
-\end{equation*}
-$$
-
-which we rewrite as
-
-$$
-\begin{equation*}
-\mu=ny\sum_{\nu=0}^n\left(\begin{array}{c} n-1 \\ \nu\end{array}\right)y^{\nu}(1-y)^{n-1-\nu} =ny(y+1-y)^{n-1}=ny.
-\end{equation*}
-$$
-
-
-
-The variance is slightly trickier to get. It reads \( \sigma^2=ny(1-y) \).
-
-
Probability Distribution Functions, Poisson's distribution
-
-
-
-
-
-Another important distribution with discrete stochastic variables \( x \) is
-the Poisson model, which resembles the exponential distribution and reads
-$$
-\begin{equation*}
- p(x) = \frac{\lambda^x}{x!} e^{-\lambda} \hspace{0.5cm}x=0,1,\dots,;\lambda > 0.
-\end{equation*}
-$$
-
-In this case both the mean value and the variance are easier to calculate,
-
-$$
-\begin{equation*}
-\mu = \sum_{x=0}^{\infty} x \frac{\lambda^x}{x!} e^{-\lambda} = \lambda e^{-\lambda}\sum_{x=1}^{\infty}
-\frac{\lambda^{x-1}}{(x-1)!}=\lambda,
-\end{equation*}
-$$
-
-and the variance is \( \sigma^2=\lambda \).
-
Probability Distribution Functions, Poisson's distribution
-
-
-
-An example of applications of the Poisson distribution could be the counting
-of the number of \( \alpha \)-particles emitted from a radioactive source in a given time interval.
-In the limit of \( n\rightarrow \infty \) and for small probabilities \( y \), the binomial distribution
-approaches the Poisson distribution. Setting \( \lambda = ny \), with \( y \) the probability for an event in
-the binomial distribution we can show that
-
-$$
-\begin{equation*}
-\lim_{n\rightarrow \infty}\left(\begin{array}{c} n \\ x\end{array}\right)y^x(1-y)^{n-x} e^{-\lambda}=\sum_{x=1}^{\infty}\frac{\lambda^x}{x!} e^{-\lambda}.
-\end{equation*}
-$$
-
-An important quantity in a statistical analysis is the so-called covariance.
-
-
-Consider the set \( \{X_i\} \) of \( n \)
-stochastic variables (not necessarily uncorrelated) with the
-multivariate PDF \( P(x_1,\dots,x_n) \). The covariance of two
-of the stochastic variables, \( X_i \) and \( X_j \), is defined as follows
-
-$$
-\begin{align}
-\mathrm{Cov}(X_i,\,X_j) & = \langle (x_i-\langle x_i\rangle)(x_j-\langle x_j\rangle)\rangle
-\tag{4}\\
-&=\int\cdots\int (x_i-\langle x_i\rangle)(x_j-\langle x_j\rangle)P(x_1,\dots,x_n)\,dx_1\dots dx_n,
-\tag{5}
-\end{align}
-$$
-
-with
-$$
-\begin{equation*}
-\langle x_i\rangle =
-\int\cdots\int x_i P(x_1,\dots,x_n)\,dx_1\dots dx_n.
-\end{equation*}
-$$
-
-If we consider the above covariance as a matrix
-$$
-C_{ij} =\mathrm{Cov}(X_i,\,X_j),
-$$
-
-then the diagonal elements are just the familiar
-variances, \( C_{ii} = \mathrm{Cov}(X_i,\,X_i) = \mathrm{Var}(X_i) \). It turns out that
-all the off-diagonal elements are zero if the stochastic variables are
-uncorrelated.
-
-Now that we have constructed an idealized mathematical framework, let
-us try to apply it to empirical observations. Examples of relevant
-physical phenomena may be spontaneous decays of nuclei, or a purely
-mathematical set of numbers produced by some deterministic
-mechanism. It is the latter we will deal with, using so-called pseudo-random
-number generators. In general our observations will contain only a limited set of
-observables. We remind the reader that
-a stochastic process is a process that produces sequentially a
-chain of values
-$$
-\begin{equation*}
-\{x_1, x_2,\dots\,x_k,\dots\}.
-\end{equation*}
-$$
-
-We will call these
-values our measurements and the entire set as our measured
-sample. The action of measuring all the elements of a sample
-we will call a stochastic experiment (since, operationally,
-they are often associated with results of empirical observation of
-some physical or mathematical phenomena; precisely an experiment). We
-assume that these values are distributed according to some
-PDF \( p_X^{\phantom X}(x) \), where \( X \) is just the formal symbol for the
-stochastic variable whose PDF is \( p_X^{\phantom X}(x) \). Instead of
-trying to determine the full distribution \( p \) we are often only
-interested in finding the few lowest moments, like the mean
-\( \mu_X^{\phantom X} \) and the variance \( \sigma_X^{\phantom X} \).
-
Numerical experiments and the covariance, actual situations
-
-
-
-In practical situations however, a sample is always of finite size. Let that
-size be \( n \). The expectation value of a sample \( \alpha \), the sample mean, is then defined as follows
-$$
-\begin{equation*}
-\langle x_{\alpha} \rangle \equiv \frac{1}{n}\sum_{k=1}^n x_{\alpha,k}.
-\end{equation*}
-$$
-
-The sample variance is:
-$$
-\begin{equation*}
-\mathrm{Var}(x) \equiv \frac{1}{n}\sum_{k=1}^n (x_{\alpha,k} - \langle x_{\alpha} \rangle)^2,
-\end{equation*}
-$$
-
-with its square root being the standard deviation of the sample.
-
Numerical experiments and the covariance, our observables
-
-
-
-You can think of the above observables as a set of quantities which define
-a given experiment. This experiment is then repeated several times, say \( m \) times.
-The total average is then
-$$
-\begin{equation}
-\langle X_m \rangle= \frac{1}{m}\sum_{\alpha=1}^mx_{\alpha}=\frac{1}{mn}\sum_{\alpha, k} x_{\alpha,k},
-\tag{6}
-\end{equation}
-$$
-
-where the last sums end at \( m \) and \( n \).
-The total variance is
-$$
-\begin{equation*}
-\sigma^2_m= \frac{1}{mn^2}\sum_{\alpha=1}^m(\langle x_{\alpha} \rangle-\langle X_m \rangle)^2,
-\end{equation*}
-$$
-
-which we rewrite as
-$$
-\begin{equation}
-\sigma^2_m=\frac{1}{m}\sum_{\alpha=1}^m\sum_{kl=1}^n (x_{\alpha,k}-\langle X_m \rangle)(x_{\alpha,l}-\langle X_m \rangle).
-\tag{7}
-\end{equation}
-$$
-
Numerical experiments and the covariance, the sample variance
-
-
-
-
-
-We define also the sample variance \( \sigma^2 \) of all \( mn \) individual experiments as
-$$
-\begin{equation}
-\sigma^2=\frac{1}{mn}\sum_{\alpha=1}^m\sum_{k=1}^n (x_{\alpha,k}-\langle X_m \rangle)^2.
-\tag{8}
-\end{equation}
-$$
-
-
-These quantities, being known experimental values or the results from our calculations,
-may differ, in some cases
-significantly, from the similarly named
-exact values for the mean value \( \mu_X \), the variance \( \mathrm{Var}(X) \)
-and the covariance \( \mathrm{Cov}(X,Y) \).
-
Numerical experiments and the covariance, central limit theorem
-
-
-
-
-
-The central limit theorem states that the PDF \( \tilde{p}(z) \) of
-the average of \( m \) random values corresponding to a PDF \( p(x) \)
-is a normal distribution whose mean is the
-mean value of the PDF \( p(x) \) and whose variance is the variance
-of the PDF \( p(x) \) divided by \( m \), the number of values used to compute \( z \).
-
-
-The central limit theorem leads then to the well-known expression for the
-standard deviation, given by
-$$
-\begin{equation*}
- \sigma_m=
-\frac{\sigma}{\sqrt{m}}.
-\end{equation*}
-$$
-
-
-In many cases the above estimate for the standard deviation, in particular if correlations are strong, may be too simplistic. We need therefore a more precise defintion of the error and the variance in our results.
-
Definition of Correlation Functions and Standard Deviation
-
-
-
-Our estimate of the true average \( \mu_{X} \) is the sample mean \( \langle X_m \rangle \)
-
-$$
-\begin{equation*}
-\mu_{X}^{\phantom X} \approx X_m=\frac{1}{mn}\sum_{\alpha=1}^m\sum_{k=1}^n x_{\alpha,k}.
-\end{equation*}
-$$
-
-
-We can then use Eq. (7)
-$$
-\begin{equation*}
-\sigma^2_m=\frac{1}{mn^2}\sum_{\alpha=1}^m\sum_{kl=1}^n (x_{\alpha,k}-\langle X_m \rangle)(x_{\alpha,l}-\langle X_m \rangle),
-\end{equation*}
-$$
-
-and rewrite it as
-$$
-\begin{equation*}
-\sigma^2_m=\frac{\sigma^2}{n}+\frac{2}{mn^2}\sum_{\alpha=1}^m\sum_{k < l}^n (x_{\alpha,k}-\langle X_m \rangle)(x_{\alpha,l}-\langle X_m \rangle),
-\end{equation*}
-$$
-
-where the first term is the sample variance of all \( mn \) experiments divided by \( n \)
-and the last term is nothing but the covariance which arises when \( k\ne l \).
-
Definition of Correlation Functions and Standard Deviation
-
-
-
-Our estimate of the true average \( \mu_{X} \) is the sample mean \( \langle X_m \rangle \)
-
-
-If the
-observables are uncorrelated, then the covariance is zero and we obtain a total variance
-which agrees with the central limit theorem. Correlations may often be present in our data set, resulting in a non-zero covariance. The first term is normally called the uncorrelated
-contribution.
-Computationally the uncorrelated first term is much easier to treat
-efficiently than the second.
-We just accumulate separately the values \( x^2 \) and \( x \) for every
-measurement \( x \) we receive. The correlation term, though, has to be
-calculated at the end of the experiment since we need all the
-measurements to calculate the cross terms. Therefore, all measurements
-have to be stored throughout the experiment.
-
Definition of Correlation Functions and Standard Deviation
-
-
-
-
-
-Let us analyze the problem by splitting up the correlation term into
-partial sums of the form
-
-$$
-\begin{equation*}
-f_d = \frac{1}{nm}\sum_{\alpha=1}^m\sum_{k=1}^{n-d}(x_{\alpha,k}-\langle X_m \rangle)(x_{\alpha,k+d}-\langle X_m \rangle),
-\end{equation*}
-$$
-
-The correlation term of the total variance can now be rewritten in terms of
-\( f_d \)
-
-$$
-\begin{equation*}
-\frac{2}{mn^2}\sum_{\alpha=1}^m\sum_{k < l}^n (x_{\alpha,k}-\langle X_m \rangle)(x_{\alpha,l}-\langle X_m \rangle)=
-\frac{2}{n}\sum_{d=1}^{n-1} f_d
-\end{equation*}
-$$
-
Definition of Correlation Functions and Standard Deviation
-
-
-
-The value of \( f_d \) reflects the correlation between measurements
-separated by the distance \( d \) in the samples. Notice that for
-\( d=0 \), \( f \) is just the sample variance, \( \sigma^2 \). If we divide \( f_d \)
-by \( \sigma^2 \), we arrive at the so called autocorrelation function
-
-$$
-\begin{equation}
-\kappa_d = \frac{f_d}{\sigma^2}
-\tag{9}
-\end{equation}
-$$
-
-which gives us a useful measure of the correlation pair correlation
-starting always at \( 1 \) for \( d=0 \).
-
Definition of Correlation Functions and Standard Deviation, sample variance
-
-
-
-
-
-The sample variance of the \( mn \) experiments can now be
-written in terms of the autocorrelation function
-
-$$
-\begin{equation}
-\sigma_m^2=\frac{\sigma^2}{n}+\frac{2}{n}\cdot\sigma^2\sum_{d=1}^{n-1}
-\frac{f_d}{\sigma^2}=\left(1+2\sum_{d=1}^{n-1}\kappa_d\right)\frac{1}{n}\sigma^2=\frac{\tau}{n}\cdot\sigma^2
-\tag{10}
-\end{equation}
-$$
-
-and we see that \( \sigma_m \) can be expressed in terms of the
-uncorrelated sample variance times a correction factor \( \tau \) which
-accounts for the correlation between measurements. We call this
-correction factor the autocorrelation time
-
-$$
-\begin{equation}
-\tau = 1+2\sum_{d=1}^{n-1}\kappa_d
-\tag{11}
-\end{equation}
-$$
-
-
-
-For a correlation free experiment, \( \tau \)
-equals 1.
-
Definition of Correlation Functions and Standard Deviation
-
-
-
-From the point of view of
-Eq. (10) we can interpret a sequential
-correlation as an effective reduction of the number of measurements by
-a factor \( \tau \). The effective number of measurements becomes
-$$
-\begin{equation*}
-n_\mathrm{eff} = \frac{n}{\tau}
-\end{equation*}
-$$
-
-To neglect the autocorrelation time \( \tau \) will always cause our
-simple uncorrelated estimate of \( \sigma_m^2\approx \sigma^2/n \) to
-be less than the true sample error. The estimate of the error will be
-too "good". On the other hand, the calculation of the full
-autocorrelation time poses an efficiency problem if the set of
-measurements is very large. The solution to this problem is given by
-more practically oriented methods like the blocking technique.
-
-
-Uniform deviates are just random numbers that lie within a specified range
-(typically 0 to 1), with any one number in the range just as likely as any other. They
-are, in other words, what you probably think random numbers are. However,
-we want to distinguish uniform deviates from other sorts of random numbers, for
-example numbers drawn from a normal (Gaussian) distribution of specified mean
-and standard deviation. These other sorts of deviates are almost always generated by
-performing appropriate operations on one or more uniform deviates, as we will see
-in subsequent sections. So, a reliable source of random uniform deviates, the subject
-of this section, is an essential building block for any sort of stochastic modeling
-or Monte Carlo computer work.
-
Random Numbers, better name: pseudo random numbers
-
-
-
-
-
-A disclaimer is however appropriate. It should be fairly obvious that
-something as deterministic as a computer cannot generate purely random numbers.
-
-
-Numbers generated by any of the standard algorithms are in reality pseudo random
-numbers, hopefully abiding to the following criteria:
-
-
-
they produce a uniform distribution in the interval [0,1].
-
correlations between random numbers are negligible
-
the period before the same sequence of random numbers is repeated is as large as possible and finally
- The most common random number generators are based on so-called
-Linear congruential relations of the type
-
-$$
-\begin{equation*}
- N_i=(aN_{i-1}+c) \mathrm{MOD} (M),
-\end{equation*}
-$$
-
-which yield a number in the interval [0,1] through
-
-$$
-\begin{equation*}
- x_i=N_i/M
-\end{equation*}
-$$
-
-
-The number
-\( M \) is called the period and it should be as large as possible
- and
-\( N_0 \) is the starting value, or seed. The function \( \mathrm{MOD} \) means the remainder,
-that is if we were to evaluate \( (13)\mathrm{MOD}(9) \), the outcome is the remainder
-of the division \( 13/9 \), namely \( 4 \).
-
-The problem with such generators is that their outputs are periodic;
-they
-will start to repeat themselves with a period that is at most \( M \). If however
-the parameters \( a \) and \( c \) are badly chosen, the period may be even shorter.
-
-
-Consider the following example
-
-$$
-\begin{equation*}
- N_i=(6N_{i-1}+7) \mathrm{MOD} (5),
-\end{equation*}
-$$
-
-with a seed \( N_0=2 \). This generator produces the sequence
-\( 4,1,3,0,2,4,1,3,0,2,...\dots \), i.e., a sequence with period \( 5 \).
-However, increasing \( M \) may not guarantee a larger period as the following
-example shows
-
-$$
-\begin{equation*}
- N_i=(27N_{i-1}+11) \mathrm{MOD} (54),
-\end{equation*}
-$$
-
-which still, with \( N_0=2 \), results in \( 11,38,11,38,11,38,\dots \), a period of
-just \( 2 \).
-
-Typical periods for the random generators provided in the program library
-are of the order of \( \sim 10^9 \) or larger. Other random number generators which have
-become increasingly popular are so-called shift-register generators.
-In these generators each successive number depends on many preceding
-values (rather than the last values as in the linear congruential
-generator).
-For example, you could make a shift register generator whose $l$th
-number is the sum of the $l-i$th and $l-j$th values with modulo \( M \),
-$$
-\begin{equation*}
- N_l=(aN_{l-i}+cN_{l-j})\mathrm{MOD}(M).
-\end{equation*}
-$$
-
-Such a generator again produces a sequence of pseudorandom numbers
-but this time with a period much larger than \( M \).
-It is also possible to construct more elaborate algorithms by including
-more than two past terms in the sum of each iteration.
-One example is the generator of Marsaglia and Zaman
-which consists of two congruential relations
-
-$$
-\begin{equation}
- N_l=(N_{l-3}-N_{l-1})\mathrm{MOD}(2^{31}-69),
-\tag{12}
-\end{equation}
-$$
-
-followed by
-$$
-\begin{equation}
- N_l=(69069N_{l-1}+1013904243)\mathrm{MOD}(2^{32}),
-\tag{13}
-\end{equation}
-$$
-
-which according to the authors has a period larger than \( 2^{94} \).
-
-Instead of using modular addition, we could use the bitwise
-exclusive-OR (\( \oplus \)) operation so that
-
-$$
-\begin{equation*}
- N_l=(N_{l-i})\oplus (N_{l-j})
-\end{equation*}
-$$
-
-where the bitwise action of \( \oplus \) means that if \( N_{l-i}=N_{l-j} \) the result is
-\( 0 \) whereas if \( N_{l-i}\ne N_{l-j} \) the result is
-\( 1 \). As an example, consider the case where \( N_{l-i}=6 \) and \( N_{l-j}=11 \). The first
-one has a bit representation (using 4 bits only) which reads \( 0110 \) whereas the
-second number is \( 1011 \). Employing the \( \oplus \) operator yields
-\( 1101 \), or \( 2^3+2^2+2^0=13 \).
-
-
-In Fortran90, the bitwise \( \oplus \) operation is coded through the intrinsic
-function \( \mathrm{IEOR}(m,n) \) where \( m \) and \( n \) are the input numbers, while in \( C \)
-it is given by \( m\wedge n \).
-
-We show here how the linear congruential algorithm can be implemented, namely
-$$
-\begin{equation*}
- N_i=(aN_{i-1}) \mathrm{MOD} (M).
-\end{equation*}
-$$
-
-However, since \( a \) and \( N_{i-1} \) are integers and their multiplication
-could become greater than the standard 32 bit integer, there is a trick via
-Schrage's algorithm which approximates the multiplication
-of large integers through the factorization
-$$
-\begin{equation*}
- M=aq+r,
-\end{equation*}
-$$
-
-where we have defined
-
-$$
-\begin{equation*}
- q=[M/a],
-\end{equation*}
-$$
-
-and
-$$
-\begin{equation*}
- r = M\hspace{0.1cm}\mathrm{MOD} \hspace{0.1cm}a.
-\end{equation*}
-$$
-
-where the brackets denote integer division. In the code below the numbers
-\( q \) and \( r \) are chosen so that \( r < q \).
-
-To see how this works we note first that
-$$
-\begin{equation}
-(aN_{i-1}) \mathrm{MOD} (M)= (aN_{i-1}-[N_{i-1}/q]M)\mathrm{MOD} (M),
-\tag{14}
-\end{equation}
-$$
-
-since we can add or subtract any integer multiple of \( M \) from \( aN_{i-1} \).
-The last term \( [N_{i-1}/q]M\mathrm{MOD}(M) \) is zero since the integer division
-\( [N_{i-1}/q] \) just yields a constant which is multiplied with \( M \).
-
-The term \( [N_{i-1}/q]r \) is always smaller or equal \( N_{i-1}(r/q) \) and with \( r < q \) we obtain always a
-number smaller than \( N_{i-1} \), which is smaller than \( M \).
-And since the number \( N_{i-1}\mathrm{MOD} (q) \) is between zero and \( q-1 \) then
-\( a(N_{i-1}\mathrm{MOD} (q)) < aq \). Combined with our definition of \( q=[M/a] \) ensures that
-this term is also smaller than \( M \) meaning that both terms fit into a
-32-bit signed integer. None of these two terms can be negative, but their difference could.
-The algorithm below adds \( M \) if their difference is negative.
-Note that the program uses the bitwise \( \oplus \) operator to generate
-the starting point for each generation of a random number. The period
-of \( ran0 \) is \( \sim 2.1\times 10^{9} \). A special feature of this
-algorithm is that is should never be called with the initial seed
-set to \( 0 \).
-
/*
- ** The function
- ** ran0()
- ** is an "Minimal" random number generator of Park and Miller
- ** Set or reset the input value
- ** idum to any integer value (except the unlikely value MASK)
- ** to initialize the sequence; idum must not be altered between
- ** calls for sucessive deviates in a sequence.
- ** The function returns a uniform deviate between 0.0 and 1.0.
- */
-doubleran0(long&idum)
-{
- constint a =16807, m =2147483647, q =127773;
- constint r =2836, MASK =123459876;
- constdouble am =1./m;
- long k;
- double ans;
- idum ^= MASK;
- k = (*idum)/q;
- idum = a*(idum - k*q) - r*k;
- // add m if negative difference
- if(idum <0) idum += m;
- ans=am*(idum);
- idum ^= MASK;
- return ans;
-} // End: function ran0()
-
-As mentioned previously, the underlying PDF for the generation of
-random numbers is the uniform distribution, meaning that the
-probability for finding a number \( x \) in the interval [0,1] is \( p(x)=1 \).
-
-
-A random number generator should produce numbers which are uniformly distributed
-in this interval. The table shows the distribution of \( N=10000 \) random
-numbers generated by the functions in the program library.
-We note in this table that the number of points in the various
-intervals \( 0.0-0.1 \), \( 0.1-0.2 \) etc are fairly close to \( 1000 \), with some minor
-deviations.
-
-
-Two additional measures are the standard deviation \( \sigma \) and the mean
-\( \mu=\langle x\rangle \).
-
-For the uniform distribution, the mean value \( \mu \) is then
-
-$$
-\begin{equation*}
- \mu=\langle x\rangle=\frac{1}{2}
-\end{equation*}
-$$
-
-while the standard deviation is
-
-$$
-\begin{equation*}
- \sigma=\sqrt{\langle x^2\rangle-\mu^2}=\frac{1}{\sqrt{12}}=0.2886.
-\end{equation*}
-$$
-
-The following simple Python code plots the distribution of the produced random numbers using the linear congruential RNG employed by Python. The trend displayed in the previous table is seen rather clearly.
-
-Since our random numbers, which are typically generated via a linear congruential algorithm,
-are never fully independent, we can then define
-an important test which measures the degree of correlation, namely the so-called
-auto-correlation function defined previously, see again Eq. (9).
-We rewrite it here as
-$$
-\begin{equation*}
- C_k=\frac{f_d}
- {\sigma^2},
-\end{equation*}
-$$
-
-with \( C_0=1 \). Recall that
-\( \sigma^2=\langle x_i^2\rangle-\langle x_i\rangle^2 \) and that
-$$
-\begin{equation*}
-f_d = \frac{1}{nm}\sum_{\alpha=1}^m\sum_{k=1}^{n-d}(x_{\alpha,k}-\langle X_m \rangle)(x_{\alpha,k+d}-\langle X_m \rangle),
-\end{equation*}
-$$
-
-
-The non-vanishing of \( C_k \) for \( k\ne 0 \) means that the random
-numbers are not independent. The independence of the random numbers is crucial
-in the evaluation of other expectation values. If they are not independent, our
-assumption for approximating \( \sigma_N \) is no longer valid.
-
-
Correlation function and which random number generators should I use
-
-
-
-The program here computes the correlation function for one of the standard functions included with the c++ compiler.
-
-
-
-
// This function computes the autocorrelation function for
-// the standard c++ random number generator
-
-#include<fstream>
-#include<iomanip>
-#include<iostream>
-#include<cmath>
-usingnamespace std;
-// output file as global variable
-ofstream ofile;
-
-// Main function begins here
-intmain(int argc, char* argv[])
-{
- int n;
- char*outfilename;
-
- cin >> n;
- double MCint =0.; double MCintsqr2=0.;
- double invers_period =1./RAND_MAX; // initialise the random number generator
- srand(time(NULL)); // This produces the so-called seed in MC jargon
- // Compute the variance and the mean value of the uniform distribution
- // Compute also the specific values x for each cycle in order to be able to
- // the covariance and the correlation function
- // Read in output file, abort if there are too few command-line arguments
- if( argc <=2 ){
- cout <<"Bad Usage: "<< argv[0] <<
- " read also output file and number of cycles on same line"<< endl;
- exit(1);
- }
- else{
- outfilename=argv[1];
- }
- ofile.open(outfilename);
- // Get the number of Monte-Carlo samples
- n = atoi(argv[2]);
- double*X;
- X =newdouble[n];
- for (int i =0; i < n; i++){
- double x =double(rand())*invers_period;
- X[i] = x;
- MCint += x;
- MCintsqr2 += x*x;
- }
- double Mean = MCint/((double) n );
- MCintsqr2 = MCintsqr2/((double) n );
- double STDev = sqrt(MCintsqr2-Mean*Mean);
- double Variance = MCintsqr2-Mean*Mean;
-// Write mean value and standard deviation
- cout <<" Standard deviation= "<< STDev <<" Integral = "<< Mean << endl;
-
- // Now we compute the autocorrelation function
- double*autocor; autocor =newdouble[n];
- for (int j =0; j < n; j++){
- double sum =0.0;
- for (int k =0; k < (n-j); k++){
- sum += (X[k]-Mean)*(X[k+j]-Mean);
- }
- autocor[j] = sum/Variance/((double) n );
- ofile << setiosflags(ios::showpoint | ios::uppercase);
- ofile << setw(15) << setprecision(8) << j;
- ofile << setw(15) << setprecision(8) << autocor[j] << endl;
- }
- ofile.close(); // close output file
- return0;
-} // end of main program
-
In the library files lib.cpp and lib.h we have included four popular RNGs taken from the widely used textbook Numerical Recipes. These are called ran0, ran1, ran2 and ran3.
-
C++ has a class called random. The random class contains a large selection of RNGs and is highly recommended. Some of these RNGs have very large periods making it thereby very safe to use these RNGs in case one is performing large calculations. In particular, the Mersenne twister random number engine has a period of \( 2^{19937} \).
-The following part of a c++ code (from project 4) sets up the uniform distribution for \( x\in [0,1] \).
-
-
-
-
/*
-
-// You need this
-#include<random>
-
-// Initialize the seed and call the Mersienne algo
-std::random_device rd;
-std::mt19937_64 gen(rd());
-// Set up the uniform distribution for x \in [[0, 1]
-std::uniform_real_distribution<double> RandomNumberGenerator(0.0,1.0);
-
-// Now use the RNG
-int ix = (int) (RandomNumberGenerator(gen)*NSpins);
-
As in other experiments, Monte Carlo experiments have two classes of errors:
-
-
-
Statistical errors
-
Systematical errors
-
-
-
Statistical errors can be estimated using standard tools from statistics
-
Systematical errors are method specific and must be treated differently from case to case. (In VMC a common source is the step length or time step in importance sampling)
Code to demonstrate the calculation of the autocorrelation function
-The following code computes the autocorrelation function, the covariance and the standard deviation
-for standard RNG.
-The following file gives the code.
-
-
-
-
// This function computes the autocorrelation function for
-// the Mersenne random number generator with a uniform distribution
-#include<iostream>
-#include<fstream>
-#include<iomanip>
-#include<cstdlib>
-#include<random>
-#include<armadillo>
-#include<string>
-#include<cmath>
-usingnamespace std;
-usingnamespace arma;
-// output file
-ofstream ofile;
-
-// Main function begins here
-intmain(int argc, char* argv[])
-{
- int MonteCarloCycles;
- string filename;
- if (argc >1) {
- filename=argv[1];
- MonteCarloCycles = atoi(argv[2]);
- string fileout = filename;
- string argument = to_string(MonteCarloCycles);
- fileout.append(argument);
- ofile.open(fileout);
- }
-
- // Compute the variance and the mean value of the uniform distribution
- // Compute also the specific values x for each cycle in order to be able to
- // compute the covariance and the correlation function
-
- vec X = zeros<vec>(MonteCarloCycles);
- double MCint =0.; double MCintsqr2=0.;
- std::random_device rd;
- std::mt19937_64 gen(rd());
- // Set up the uniform distribution for x \in [[0, 1]
- std::uniform_real_distribution<double> RandomNumberGenerator(0.0,1.0);
- for (int i =0; i < MonteCarloCycles; i++){
- double x = RandomNumberGenerator(gen);
- X(i) = x;
- MCint += x;
- MCintsqr2 += x*x;
- }
- double Mean = MCint/((double) MonteCarloCycles );
- MCintsqr2 = MCintsqr2/((double) MonteCarloCycles );
- double STDev = sqrt(MCintsqr2-Mean*Mean);
- double Variance = MCintsqr2-Mean*Mean;
- // Write mean value and variance
- cout <<" Sample variance= "<< Variance <<" Mean value = "<< Mean << endl;
- // Now we compute the autocorrelation function
- vec autocorrelation = zeros<vec>(MonteCarloCycles);
- for (int j =0; j < MonteCarloCycles; j++){
- double sum =0.0;
- for (int k =0; k < (MonteCarloCycles-j); k++){
- sum += (X(k)-Mean)*(X(k+j)-Mean);
- }
- autocorrelation(j) = sum/Variance/((double) MonteCarloCycles );
- ofile << setiosflags(ios::showpoint | ios::uppercase);
- ofile << setw(15) << setprecision(8) << j;
- ofile << setw(15) << setprecision(8) << autocorrelation(j) << endl;
- }
- // Now compute the exact covariance using the autocorrelation function
- double Covariance =0.0;
- for (int j =0; j < MonteCarloCycles; j++){
- Covariance += autocorrelation(j);
- }
- Covariance *=2.0/((double) MonteCarloCycles);
- // Compute now the total variance, including the covariance, and obtain the standard deviation
- double TotalVariance = (Variance/((double) MonteCarloCycles ))+Covariance;
- cout <<"Covariance ="<< Covariance <<"Totalvariance= "<< TotalVariance <<"Sample Variance/n= "<< (Variance/((double) MonteCarloCycles )) << endl;
- cout <<" STD from sample variance= "<< sqrt(Variance/((double) MonteCarloCycles )) <<" STD with covariance = "<< sqrt(TotalVariance) << endl;
-
- ofile.close(); // close output file
- return0;
-} // end of main program
-
Say that we have a set of samples from a Monte Carlo experiment
-
Assuming (wrongly) that our samples are uncorrelated our best estimate of the standard deviation of the mean \( \langle \mathbf{M}\rangle \) is given by
If the samples are correlated we can rewrite our results to show that
-
-
-$$
-\sigma=\sqrt{\frac{1+2\tau/\Delta t}{n}\left(\langle \mathbf{M}^2\rangle-\langle \mathbf{M}\rangle^2\right)}
-$$
-
- where \( \tau \) is the correlation time (the time between a sample and the next uncorrelated sample) and \( \Delta t \) is time between each sample
-
Problem: We don't know \( \tau \) or it is too expensive to compute
-
Solution: Make a plot of std. dev. as a function of blocksize
-
The estimate of std. dev. of correlated data is too low \( \to \) the error will increase with increasing block size until the blocks are uncorrelated, where we reach a plateau
-
When the std. dev. stops increasing the blocks are uncorrelated
-When the file gets large, it can be useful to write your data in binary mode instead of ascii characters.
-The following python file reads data from file with the output from every Monte Carlo cycle.
-
-Given a set of \( N \) data, assume that we are interested in some
-observable \( \theta \) which may be estimated from that set. This observable can also be for example the result of a fit based on all \( N \) raw data.
-Let us call the value of the observable obtained from the original
-data set \( \hat{\theta} \). One recreates from the sample repeatedly
-other samples by choosing randomly \( N \) data out of the original set.
-This costs essentially nothing, since we just recycle the original data set for the building of new sets.
-
-
-Let us assume we have done this \( K \) times and thus have \( K \) sets of \( N \)
-data values each.
-Of course some values will enter more than once in the new sets. For each of these sets one computes the observable \( \theta \) resulting in values \( \theta_k \) with \( k = 1,...,K \). Then one determines
-$$
-\tilde{\theta} = \frac{1}{K} \sum_{k=1}^K \theta_k,
-$$
-
-and
-$$
-sigma^2_{\tilde{\theta}} = \frac{1}{K} \sum_{k=1}^K \left(\theta_k-\tilde{\theta}\right)^2.
-$$
-
-
-These are estimators for \( \angle\theta\rangle \) and its variance. They are not unbiased and therefore
-\( \tilde{\theta}\neq\hat{\theta} \) for finite K.
-
-
-The difference is called bias and gives an idea on how far away the result may be from
-the true \( \angle\theta\rangle \). As final result for the observable one quotes \( \angle\theta\rangle = \tilde{\theta} \pm \sigma_{\tilde{\theta}} \) .
-
-
-The numbers in the domain are the outcomes of the physical process tossing the dice.
-We cannot tell beforehand whether the outcome is 3 or 5 or any other number in this domain.
-This defines the randomness of the outcome, or unexpectedness or any other synonimous word which
-encompasses the uncertitude of the final outcome.
-
-
-The only thing we can tell beforehand
-is that say the outcome 2 has a certain probability.
-If our favorite hobby is to spend an hour every evening throwing dice and
-registering the sequence of outcomes, we will note that the numbers in the above domain
-
-Random variables are characterized by a domain which contains all possible values that the random value may take. This domain has a corresponding PDF.
-
-
-
-
-
-
Stochastic variables and the main concepts, the discrete case
-
-
-
-There are two main concepts associated with a stochastic variable. The
-domain is the set \( \mathbb D = \{x\} \) of all accessible values
-the variable can assume, so that \( X \in \mathbb D \). An example of a
-discrete domain is the set of six different numbers that we may get by
-throwing of a dice, \( x\in\{1,\,2,\,3,\,4,\,5,\,6\} \).
-
-
-The probability distribution function (PDF) is a function
-\( p(x) \) on the domain which, in the discrete case, gives us the
-probability or relative frequency with which these values of \( X \)
-occur
-
Stochastic variables and the main concepts, the continuous case
-
-
-
-In the continuous case, the PDF does not directly depict the
-actual probability. Instead we define the probability for the
-stochastic variable to assume any value on an infinitesimal interval
-around \( x \) to be \( p(x)dx \). The continuous function \( p(x) \) then gives us
-the density of the probability rather than the probability
-itself. The probability for a stochastic variable to assume any value
-on a non-infinitesimal interval \( [a,\,b] \) is then just the integral
-
-
-
-Qualitatively speaking, a stochastic variable represents the values of
-numbers chosen as if by chance from some specified PDF so that the
-selection of a large set of these numbers reproduces this PDF.
-
-
-
-
-
-
The cumulative probability
-
-
-
-Of interest to us is the cumulative probability
-distribution function (CDF), \( P(x) \), which is just the probability
-for a stochastic variable \( X \) to assume any value less than \( x \)
-
-
-Naturally, it would be nonsensical for any of the values of the domain
-to occur with a probability greater than \( 1 \) or less than \( 0 \). Also,
-the PDF must be normalized. That is, all the probabilities must add up
-to unity. The probability of "anything" to happen is always unity. For
-both discrete and continuous PDFs, this condition is
-
-
-with mean value \( \mu \) and standard deviation \( \sigma \). If \( \mu=0 \) and \( \sigma=1 \), it is normally called the standard normal distribution
-
-Let \( h(x) \) be an arbitrary continuous function on the domain of the stochastic
-variable \( X \) whose PDF is \( p(x) \). We define the expectation value
-of \( h \) with respect to \( p \) as follows
-
-
-
-Whenever the PDF is known implicitly, like in this case, we will drop
-the index \( X \) for clarity.
-A particularly useful class of special expectation values are the
-moments. The \( n \)-th moment of the PDF \( p \) is defined as
-follows
-
Stochastic variables and the main concepts, mean values
-
-
-
-The zero-th moment \( \langle 1\rangle \) is just the normalization condition of
-\( p \). The first moment, \( \langle x\rangle \), is called the mean of \( p \)
-and often denoted by the letter \( \mu \)
-
-
-for a discrete distribution.
-Qualitatively it represents the centroid or the average value of the
-PDF and is therefore simply called the expectation value of \( p(x) \).
-
-
-
-
-
-
Stochastic variables and the main concepts, central moments, the variance
-
-
-
-A special version of the moments is the set of central moments, the n-th central moment defined as
-
-
-The zero-th and first central moments are both trivial, equal \( 1 \) and
-\( 0 \), respectively. But the second central moment, known as the
-variance of \( p \), is of particular interest. For the stochastic
-variable \( X \), the variance is denoted as \( \sigma^2_X \) or \( \mathrm{Var}(X) \)
-
-
-The square root of the variance, \( \sigma =\sqrt{\langle (x-\langle x\rangle)^2\rangle} \) is called the
-standard deviation of \( p \). It is the RMS (root-mean-square)
-value of the deviation of the PDF from its mean value, interpreted
-qualitatively as the "spread" of \( p \) around its mean.
-
-
-
-
-
-
Probability Distribution Functions
-
-
-
-The following table collects properties of probability distribution functions.
-In our notation we reserve the label \( p(x) \) for the probability of a certain event,
-while \( P(x) \) is the cumulative probability.
-
-
-
-
-
Discrete PDF
Continuous PDF
-
-
-
Domain
\( \left\{x_1, x_2, x_3, \dots, x_N\right\} \)
\( [a,b] \)
-
Probability
\( p(x_i) \)
\( p(x)dx \)
-
Cumulative
\( P_i=\sum_{l=1}^ip(x_l) \)
\( P(x)=\int_a^xp(t)dt \)
-
Positivity
$ 0\le p(x_i)\le 1$
$ p(x) \ge 0$
-
Positivity
$ 0\le P_i\le 1$
$ 0\le P(x)\le 1$
-
Monotonic
\( P_i\ge P_j \) if \( x_i\ge x_j \)
\( P(x_i)\ge P(x_j) \) if \( x_i\ge x_j \)
-
Normalization
\( P_N=1 \)
\( P(b)=1 \)
-
-
-
-
-
-
-
-
-
Probability Distribution Functions
-
-
-
-With a PDF we can compute expectation values of selected quantities such as
-
-
-
-with probabilities different from zero in the interval \( (-\infty,\infty) \).
-The integral \( \int_{-\infty}^{\infty}\exp{\left(-(x^2\right)}dx \) appears in many calculations, its value
-is \( \sqrt{\pi} \), a result we will need when we compute the mean value and the variance.
-The mean value is
-
-
-and performing a final integration by parts we obtain the well-known result \( \sigma^2=b^2 \).
-It is useful to introduce the standard normal distribution as well, defined by \( \mu=a=0 \), viz. a distribution
-centered around zero and with a variance \( \sigma^2=1 \), leading to
-
-
Probability Distribution Functions, the cumulative distribution
-
-
-
-The exponential and uniform distributions have simple cumulative functions,
-whereas the normal distribution does not, being proportional to the so-called
-error function \( erf(x) \), given by
-
-
-
-where \( y \) is the probability for a specific event, such as the tossing of a coin or moving left or right
-in case of a random walker. Note that \( x \) is a discrete stochastic variable.
-
-
-The sequence of binomial trials is characterized by the following definitions
-
-
-
-
Every experiment is thought to consist of \( N \) independent trials.
-
-
In every independent trial one registers if a specific situation happens or not, such as the jump to the left or right of a random walker.
-
-
The probability for every outcome in a single trial has the same value, for example the outcome of tossing (either heads or tails) a coin is always \( 1/2 \).
-
-
-
-
-
-
-
Probability Distribution Functions, the binomial distribution
-
-
-
-In order to compute the mean and variance we need to recall Newton's binomial
-formula
-
-
-The variance is slightly trickier to get. It reads \( \sigma^2=ny(1-y) \).
-
-
-
-
-
Probability Distribution Functions, Poisson's distribution
-
-
-
-Another important distribution with discrete stochastic variables \( x \) is
-the Poisson model, which resembles the exponential distribution and reads
-
Probability Distribution Functions, Poisson's distribution
-
-
-
-An example of applications of the Poisson distribution could be the counting
-of the number of \( \alpha \)-particles emitted from a radioactive source in a given time interval.
-In the limit of \( n\rightarrow \infty \) and for small probabilities \( y \), the binomial distribution
-approaches the Poisson distribution. Setting \( \lambda = ny \), with \( y \) the probability for an event in
-the binomial distribution we can show that
-
-
-An important quantity in a statistical analysis is the so-called covariance.
-
-
-Consider the set \( \{X_i\} \) of \( n \)
-stochastic variables (not necessarily uncorrelated) with the
-multivariate PDF \( P(x_1,\dots,x_n) \). The covariance of two
-of the stochastic variables, \( X_i \) and \( X_j \), is defined as follows
-
-
-If we consider the above covariance as a matrix
-
-$$
-C_{ij} =\mathrm{Cov}(X_i,\,X_j),
-$$
-
-
-then the diagonal elements are just the familiar
-variances, \( C_{ii} = \mathrm{Cov}(X_i,\,X_i) = \mathrm{Var}(X_i) \). It turns out that
-all the off-diagonal elements are zero if the stochastic variables are
-uncorrelated.
-
-
-
-
-
-
Meet the covariance, uncorrelated events
-
-
-
-This is easy to show, keeping in mind the linearity of
-the expectation value. Consider the stochastic variables \( X_i \) and
-\( X_j \), (\( i\neq j \))
-
-Now that we have constructed an idealized mathematical framework, let
-us try to apply it to empirical observations. Examples of relevant
-physical phenomena may be spontaneous decays of nuclei, or a purely
-mathematical set of numbers produced by some deterministic
-mechanism. It is the latter we will deal with, using so-called pseudo-random
-number generators. In general our observations will contain only a limited set of
-observables. We remind the reader that
-a stochastic process is a process that produces sequentially a
-chain of values
-
-We will call these
-values our measurements and the entire set as our measured
-sample. The action of measuring all the elements of a sample
-we will call a stochastic experiment (since, operationally,
-they are often associated with results of empirical observation of
-some physical or mathematical phenomena; precisely an experiment). We
-assume that these values are distributed according to some
-PDF \( p_X^{\phantom X}(x) \), where \( X \) is just the formal symbol for the
-stochastic variable whose PDF is \( p_X^{\phantom X}(x) \). Instead of
-trying to determine the full distribution \( p \) we are often only
-interested in finding the few lowest moments, like the mean
-\( \mu_X^{\phantom X} \) and the variance \( \sigma_X^{\phantom X} \).
-
-
-
-
-
-
Numerical experiments and the covariance, actual situations
-
-
-
-In practical situations however, a sample is always of finite size. Let that
-size be \( n \). The expectation value of a sample \( \alpha \), the sample mean, is then defined as follows
-
-
-with its square root being the standard deviation of the sample.
-
-
-
-
-
-
Numerical experiments and the covariance, our observables
-
-
-
-You can think of the above observables as a set of quantities which define
-a given experiment. This experiment is then repeated several times, say \( m \) times.
-The total average is then
-
-These quantities, being known experimental values or the results from our calculations,
-may differ, in some cases
-significantly, from the similarly named
-exact values for the mean value \( \mu_X \), the variance \( \mathrm{Var}(X) \)
-and the covariance \( \mathrm{Cov}(X,Y) \).
-
-
-
-
-
-
Numerical experiments and the covariance, central limit theorem
-
-
-
-The central limit theorem states that the PDF \( \tilde{p}(z) \) of
-the average of \( m \) random values corresponding to a PDF \( p(x) \)
-is a normal distribution whose mean is the
-mean value of the PDF \( p(x) \) and whose variance is the variance
-of the PDF \( p(x) \) divided by \( m \), the number of values used to compute \( z \).
-
-
-The central limit theorem leads then to the well-known expression for the
-standard deviation, given by
-
-In many cases the above estimate for the standard deviation, in particular if correlations are strong, may be too simplistic. We need therefore a more precise defintion of the error and the variance in our results.
-
-
-
-
-
-
Definition of Correlation Functions and Standard Deviation
-
-
-
-Our estimate of the true average \( \mu_{X} \) is the sample mean \( \langle X_m \rangle \)
-
-
-
-where the first term is the sample variance of all \( mn \) experiments divided by \( n \)
-and the last term is nothing but the covariance which arises when \( k\ne l \).
-
-
-
-
-
-
Definition of Correlation Functions and Standard Deviation
-
-
-
-Our estimate of the true average \( \mu_{X} \) is the sample mean \( \langle X_m \rangle \)
-
-
-If the
-observables are uncorrelated, then the covariance is zero and we obtain a total variance
-which agrees with the central limit theorem. Correlations may often be present in our data set, resulting in a non-zero covariance. The first term is normally called the uncorrelated
-contribution.
-Computationally the uncorrelated first term is much easier to treat
-efficiently than the second.
-We just accumulate separately the values \( x^2 \) and \( x \) for every
-measurement \( x \) we receive. The correlation term, though, has to be
-calculated at the end of the experiment since we need all the
-measurements to calculate the cross terms. Therefore, all measurements
-have to be stored throughout the experiment.
-
-
-
-
-
-
Definition of Correlation Functions and Standard Deviation
-
-
-
-Let us analyze the problem by splitting up the correlation term into
-partial sums of the form
-
-
Definition of Correlation Functions and Standard Deviation
-
-
-
-The value of \( f_d \) reflects the correlation between measurements
-separated by the distance \( d \) in the samples. Notice that for
-\( d=0 \), \( f \) is just the sample variance, \( \sigma^2 \). If we divide \( f_d \)
-by \( \sigma^2 \), we arrive at the so called autocorrelation function
-
-
-
-and we see that \( \sigma_m \) can be expressed in terms of the
-uncorrelated sample variance times a correction factor \( \tau \) which
-accounts for the correlation between measurements. We call this
-correction factor the autocorrelation time
-
-
Definition of Correlation Functions and Standard Deviation
-
-
-
-From the point of view of
-Eq. (10) we can interpret a sequential
-correlation as an effective reduction of the number of measurements by
-a factor \( \tau \). The effective number of measurements becomes
-
-
-To neglect the autocorrelation time \( \tau \) will always cause our
-simple uncorrelated estimate of \( \sigma_m^2\approx \sigma^2/n \) to
-be less than the true sample error. The estimate of the error will be
-too "good". On the other hand, the calculation of the full
-autocorrelation time poses an efficiency problem if the set of
-measurements is very large. The solution to this problem is given by
-more practically oriented methods like the blocking technique.
-
-
-
-
-
-
-
Random Numbers
-
-
-
-Uniform deviates are just random numbers that lie within a specified range
-(typically 0 to 1), with any one number in the range just as likely as any other. They
-are, in other words, what you probably think random numbers are. However,
-we want to distinguish uniform deviates from other sorts of random numbers, for
-example numbers drawn from a normal (Gaussian) distribution of specified mean
-and standard deviation. These other sorts of deviates are almost always generated by
-performing appropriate operations on one or more uniform deviates, as we will see
-in subsequent sections. So, a reliable source of random uniform deviates, the subject
-of this section, is an essential building block for any sort of stochastic modeling
-or Monte Carlo computer work.
-
-
-
-
-
-
Random Numbers, better name: pseudo random numbers
-
-
-
-A disclaimer is however appropriate. It should be fairly obvious that
-something as deterministic as a computer cannot generate purely random numbers.
-
-
-Numbers generated by any of the standard algorithms are in reality pseudo random
-numbers, hopefully abiding to the following criteria:
-
-
-
-
they produce a uniform distribution in the interval [0,1].
-
-
correlations between random numbers are negligible
-
-
the period before the same sequence of random numbers is repeated is as large as possible and finally
-
-
the algorithm should be fast.
-
-
-
-
-
-
-
Random number generator RNG
-
-
-
- The most common random number generators are based on so-called
-Linear congruential relations of the type
-
-
-The number
-\( M \) is called the period and it should be as large as possible
- and
-\( N_0 \) is the starting value, or seed. The function \( \mathrm{MOD} \) means the remainder,
-that is if we were to evaluate \( (13)\mathrm{MOD}(9) \), the outcome is the remainder
-of the division \( 13/9 \), namely \( 4 \).
-
-
-
-
-
-
Random number generator RNG and periodic outputs
-
-
-
-The problem with such generators is that their outputs are periodic;
-they
-will start to repeat themselves with a period that is at most \( M \). If however
-the parameters \( a \) and \( c \) are badly chosen, the period may be even shorter.
-
-
-
-with a seed \( N_0=2 \). This generator produces the sequence
-\( 4,1,3,0,2,4,1,3,0,2,...\dots \), i.e., a sequence with period \( 5 \).
-However, increasing \( M \) may not guarantee a larger period as the following
-example shows
-
-
-
-which still, with \( N_0=2 \), results in \( 11,38,11,38,11,38,\dots \), a period of
-just \( 2 \).
-
-
-
-
-
-
Random number generator RNG and its period
-
-
-
-Typical periods for the random generators provided in the program library
-are of the order of \( \sim 10^9 \) or larger. Other random number generators which have
-become increasingly popular are so-called shift-register generators.
-In these generators each successive number depends on many preceding
-values (rather than the last values as in the linear congruential
-generator).
-For example, you could make a shift register generator whose $l$th
-number is the sum of the $l-i$th and $l-j$th values with modulo \( M \),
-
-Such a generator again produces a sequence of pseudorandom numbers
-but this time with a period much larger than \( M \).
-It is also possible to construct more elaborate algorithms by including
-more than two past terms in the sum of each iteration.
-One example is the generator of Marsaglia and Zaman
-which consists of two congruential relations
-
-
-
-where the bitwise action of \( \oplus \) means that if \( N_{l-i}=N_{l-j} \) the result is
-\( 0 \) whereas if \( N_{l-i}\ne N_{l-j} \) the result is
-\( 1 \). As an example, consider the case where \( N_{l-i}=6 \) and \( N_{l-j}=11 \). The first
-one has a bit representation (using 4 bits only) which reads \( 0110 \) whereas the
-second number is \( 1011 \). Employing the \( \oplus \) operator yields
-\( 1101 \), or \( 2^3+2^2+2^0=13 \).
-
-
-In Fortran90, the bitwise \( \oplus \) operation is coded through the intrinsic
-function \( \mathrm{IEOR}(m,n) \) where \( m \) and \( n \) are the input numbers, while in \( C \)
-it is given by \( m\wedge n \).
-
-
-
-
-
-
Random number generator RNG, RAN0
-
-
-
-We show here how the linear congruential algorithm can be implemented, namely
-
-
-However, since \( a \) and \( N_{i-1} \) are integers and their multiplication
-could become greater than the standard 32 bit integer, there is a trick via
-Schrage's algorithm which approximates the multiplication
-of large integers through the factorization
-
-
-since we can add or subtract any integer multiple of \( M \) from \( aN_{i-1} \).
-The last term \( [N_{i-1}/q]M\mathrm{MOD}(M) \) is zero since the integer division
-\( [N_{i-1}/q] \) just yields a constant which is multiplied with \( M \).
-
-The term \( [N_{i-1}/q]r \) is always smaller or equal \( N_{i-1}(r/q) \) and with \( r < q \) we obtain always a
-number smaller than \( N_{i-1} \), which is smaller than \( M \).
-And since the number \( N_{i-1}\mathrm{MOD} (q) \) is between zero and \( q-1 \) then
-\( a(N_{i-1}\mathrm{MOD} (q)) < aq \). Combined with our definition of \( q=[M/a] \) ensures that
-this term is also smaller than \( M \) meaning that both terms fit into a
-32-bit signed integer. None of these two terms can be negative, but their difference could.
-The algorithm below adds \( M \) if their difference is negative.
-Note that the program uses the bitwise \( \oplus \) operator to generate
-the starting point for each generation of a random number. The period
-of \( ran0 \) is \( \sim 2.1\times 10^{9} \). A special feature of this
-algorithm is that is should never be called with the initial seed
-set to \( 0 \).
-
-
-
-
-
-
Random number generator RNG, RAN0 code
-
-
-
-
-
-
/*
- ** The function
- ** ran0()
- ** is an "Minimal" random number generator of Park and Miller
- ** Set or reset the input value
- ** idum to any integer value (except the unlikely value MASK)
- ** to initialize the sequence; idum must not be altered between
- ** calls for sucessive deviates in a sequence.
- ** The function returns a uniform deviate between 0.0 and 1.0.
- */
-doubleran0(long &idum)
-{
- constint a = 16807, m = 2147483647, q = 127773;
- constint r = 2836, MASK = 123459876;
- constdouble am = 1./m;
- long k;
- double ans;
- idum ^= MASK;
- k = (*idum)/q;
- idum = a*(idum - k*q) - r*k;
- // add m if negative difference
- if(idum < 0) idum += m;
- ans=am*(idum);
- idum ^= MASK;
- return ans;
-} // End: function ran0()
-
-
-
-
-
-
-
-
Properties of Selected Random Number Generators
-
-
-
-As mentioned previously, the underlying PDF for the generation of
-random numbers is the uniform distribution, meaning that the
-probability for finding a number \( x \) in the interval [0,1] is \( p(x)=1 \).
-
-
-A random number generator should produce numbers which are uniformly distributed
-in this interval. The table shows the distribution of \( N=10000 \) random
-numbers generated by the functions in the program library.
-We note in this table that the number of points in the various
-intervals \( 0.0-0.1 \), \( 0.1-0.2 \) etc are fairly close to \( 1000 \), with some minor
-deviations.
-
-
-Two additional measures are the standard deviation \( \sigma \) and the mean
-\( \mu=\langle x\rangle \).
-
-
-
-
-
-
Properties of Selected Random Number Generators
-
-
-
-For the uniform distribution, the mean value \( \mu \) is then
-
-
-The various random number generators produce results which agree rather well with
-these limiting values.
-
-
-
-
-
\( x \)-bin
ran0
ran1
ran2
ran3
-
-
-
0.0-0.1
1013
991
938
1047
-
0.1-0.2
1002
1009
1040
1030
-
0.2-0.3
989
999
1030
993
-
0.3-0.4
939
960
1023
937
-
0.4-0.5
1038
1001
1002
992
-
0.5-0.6
1037
1047
1009
1009
-
0.6-0.7
1005
989
1003
989
-
0.7-0.8
986
962
985
954
-
0.8-0.9
1000
1027
1009
1023
-
0.9-1.0
991
1015
961
1026
-
\( \mu \)
0.4997
0.5018
0.4992
0.4990
-
\( \sigma \)
0.2882
0.2892
0.2861
0.2915
-
-
-
-
-
-
-
-
-
Simple demonstration of RNGs using python
-
-
-
-The following simple Python code plots the distribution of the produced random numbers using the linear congruential RNG employed by Python. The trend displayed in the previous table is seen rather clearly.
-
-
-
-
-
-
-
-
-
-
-
Properties of Selected Random Number Generators
-
-
-
-Since our random numbers, which are typically generated via a linear congruential algorithm,
-are never fully independent, we can then define
-an important test which measures the degree of correlation, namely the so-called
-auto-correlation function defined previously, see again Eq. (9).
-We rewrite it here as
-
-The non-vanishing of \( C_k \) for \( k\ne 0 \) means that the random
-numbers are not independent. The independence of the random numbers is crucial
-in the evaluation of other expectation values. If they are not independent, our
-assumption for approximating \( \sigma_N \) is no longer valid.
-
-
-
-
-
-
-
-
Correlation function and which random number generators should I use
-
-
-
-The program here computes the correlation function for one of the standard functions included with the c++ compiler.
-
-
-
-
// This function computes the autocorrelation function for
-// the standard c++ random number generator
-
-#include<fstream>
-#include<iomanip>
-#include<iostream>
-#include<cmath>
-usingnamespace std;
-// output file as global variable
-ofstream ofile;
-
-// Main function begins here
-intmain(int argc, char* argv[])
-{
- int n;
- char *outfilename;
-
- cin >> n;
- double MCint = 0.; double MCintsqr2=0.;
- double invers_period = 1./RAND_MAX; // initialise the random number generator
- srand(time(NULL)); // This produces the so-called seed in MC jargon
- // Compute the variance and the mean value of the uniform distribution
- // Compute also the specific values x for each cycle in order to be able to
- // the covariance and the correlation function
- // Read in output file, abort if there are too few command-line arguments
- if( argc <= 2 ){
- cout << "Bad Usage: " << argv[0] <<
- " read also output file and number of cycles on same line" << endl;
- exit(1);
- }
- else{
- outfilename=argv[1];
- }
- ofile.open(outfilename);
- // Get the number of Monte-Carlo samples
- n = atoi(argv[2]);
- double *X;
- X = newdouble[n];
- for (int i = 0; i < n; i++){
- double x = double(rand())*invers_period;
- X[i] = x;
- MCint += x;
- MCintsqr2 += x*x;
- }
- double Mean = MCint/((double) n );
- MCintsqr2 = MCintsqr2/((double) n );
- double STDev = sqrt(MCintsqr2-Mean*Mean);
- double Variance = MCintsqr2-Mean*Mean;
-// Write mean value and standard deviation
- cout << " Standard deviation= " << STDev << " Integral = " << Mean << endl;
-
- // Now we compute the autocorrelation function
- double *autocor; autocor = newdouble[n];
- for (int j = 0; j < n; j++){
- double sum = 0.0;
- for (int k = 0; k < (n-j); k++){
- sum += (X[k]-Mean)*(X[k+j]-Mean);
- }
- autocor[j] = sum/Variance/((double) n );
- ofile << setiosflags(ios::showpoint | ios::uppercase);
- ofile << setw(15) << setprecision(8) << j;
- ofile << setw(15) << setprecision(8) << autocor[j] << endl;
- }
- ofile.close(); // close output file
- return0;
-} // end of main program
-
-
-
-
-
-
-
-
Correlation function and which random number generators should I use
-
-
-
-The following Python code plots the results for the correlation function from the above program.
-
-
-
-
-
-
-
-
-
-
-
Which RNG should I use?
-
-
-
-
In the library files lib.cpp and lib.h we have included four popular RNGs taken from the widely used textbook Numerical Recipes. These are called ran0, ran1, ran2 and ran3.
-
C++ has a class called random. The random class contains a large selection of RNGs and is highly recommended. Some of these RNGs have very large periods making it thereby very safe to use these RNGs in case one is performing large calculations. In particular, the Mersenne twister random number engine has a period of \( 2^{19937} \).
-
-
-
-
-
-
-
How to use the Mersenne generator
-
-
-
-The following part of a c++ code (from project 4) sets up the uniform distribution for \( x\in [0,1] \).
-
-
-
-
/*
-
-// You need this
-#include<random>
-
-// Initialize the seed and call the Mersienne algo
-std::random_device rd;
-std::mt19937_64 gen(rd());
-// Set up the uniform distribution for x \in [[0, 1]
-std::uniform_real_distribution<double> RandomNumberGenerator(0.0,1.0);
-
-// Now use the RNG
-int ix = (int) (RandomNumberGenerator(gen)*NSpins);
-
-
-
-
-
-
-
-
Why blocking?
-
-Statistical analysis.
-
-
-
Monte Carlo simulations can be treated as computer experiments
-
-
The results can be analysed with the same statistical tools as we would use analysing experimental data.
-
-
As in all experiments, we are looking for expectation values and an estimate of how accurate they are, i.e., possible sources for errors.
-
-
-
-A very good article which explains blocking is H. Flyvbjerg and H. G. Petersen, Error estimates on averages of correlated data, Journal of Chemical Physics 91, 461-466 (1989).
-
-
-
-
-
-
-
-
Why blocking?
-
-Statistical analysis.
-
-
-
As in other experiments, Monte Carlo experiments have two classes of errors:
-
-
-
-
Statistical errors
-
-
Systematical errors
-
-
Statistical errors can be estimated using standard tools from statistics
-
-
Systematical errors are method specific and must be treated differently from case to case. (In VMC a common source is the step length or time step in importance sampling)
-
-
-
-
-
-
-
Code to demonstrate the calculation of the autocorrelation function
-The following code computes the autocorrelation function, the covariance and the standard deviation
-for standard RNG.
-The following file gives the code.
-
-
-
-
// This function computes the autocorrelation function for
-// the Mersenne random number generator with a uniform distribution
-#include<iostream>
-#include<fstream>
-#include<iomanip>
-#include<cstdlib>
-#include<random>
-#include<armadillo>
-#include<string>
-#include<cmath>
-usingnamespace std;
-usingnamespace arma;
-// output file
-ofstream ofile;
-
-// Main function begins here
-intmain(int argc, char* argv[])
-{
- int MonteCarloCycles;
- string filename;
- if (argc > 1) {
- filename=argv[1];
- MonteCarloCycles = atoi(argv[2]);
- string fileout = filename;
- string argument = to_string(MonteCarloCycles);
- fileout.append(argument);
- ofile.open(fileout);
- }
-
- // Compute the variance and the mean value of the uniform distribution
- // Compute also the specific values x for each cycle in order to be able to
- // compute the covariance and the correlation function
-
- vec X = zeros<vec>(MonteCarloCycles);
- double MCint = 0.; double MCintsqr2=0.;
- std::random_device rd;
- std::mt19937_64 gen(rd());
- // Set up the uniform distribution for x \in [[0, 1]
- std::uniform_real_distribution<double> RandomNumberGenerator(0.0,1.0);
- for (int i = 0; i < MonteCarloCycles; i++){
- double x = RandomNumberGenerator(gen);
- X(i) = x;
- MCint += x;
- MCintsqr2 += x*x;
- }
- double Mean = MCint/((double) MonteCarloCycles );
- MCintsqr2 = MCintsqr2/((double) MonteCarloCycles );
- double STDev = sqrt(MCintsqr2-Mean*Mean);
- double Variance = MCintsqr2-Mean*Mean;
- // Write mean value and variance
- cout << " Sample variance= " << Variance << " Mean value = " << Mean << endl;
- // Now we compute the autocorrelation function
- vec autocorrelation = zeros<vec>(MonteCarloCycles);
- for (int j = 0; j < MonteCarloCycles; j++){
- double sum = 0.0;
- for (int k = 0; k < (MonteCarloCycles-j); k++){
- sum += (X(k)-Mean)*(X(k+j)-Mean);
- }
- autocorrelation(j) = sum/Variance/((double) MonteCarloCycles );
- ofile << setiosflags(ios::showpoint | ios::uppercase);
- ofile << setw(15) << setprecision(8) << j;
- ofile << setw(15) << setprecision(8) << autocorrelation(j) << endl;
- }
- // Now compute the exact covariance using the autocorrelation function
- double Covariance = 0.0;
- for (int j = 0; j < MonteCarloCycles; j++){
- Covariance += autocorrelation(j);
- }
- Covariance *= 2.0/((double) MonteCarloCycles);
- // Compute now the total variance, including the covariance, and obtain the standard deviation
- double TotalVariance = (Variance/((double) MonteCarloCycles ))+Covariance;
- cout << "Covariance =" << Covariance << "Totalvariance= " << TotalVariance << "Sample Variance/n= " << (Variance/((double) MonteCarloCycles )) << endl;
- cout << " STD from sample variance= " << sqrt(Variance/((double) MonteCarloCycles )) << " STD with covariance = " << sqrt(TotalVariance) << endl;
-
- ofile.close(); // close output file
- return0;
-} // end of main program
-
-
-
-
-
-
What is blocking?
-
-Blocking.
-
-
-
Say that we have a set of samples from a Monte Carlo experiment
-
-
Assuming (wrongly) that our samples are uncorrelated our best estimate of the standard deviation of the mean \( \langle \mathbf{M}\rangle \) is given by
-
- where \( \tau \) is the correlation time (the time between a sample and the next uncorrelated sample) and \( \Delta t \) is time between each sample
-
-
-
-
-
-
What is blocking?
-
-Blocking.
-
-
-
If \( \Delta t\gg\tau \) our first estimate of \( \sigma \) still holds
-
-
Much more common that \( \Delta t < \tau \)
-
-
In the method of data blocking we divide the sequence of samples into blocks
-
-
We then take the mean \( \langle \mathbf{M}_i\rangle \) of block \( i=1\ldots n_{blocks} \) to calculate the total mean and variance
-
-
The size of each block must be so large that sample \( j \) of block \( i \) is not correlated with sample \( j \) of block \( i+1 \)
-
-
The correlation time \( \tau \) would be a good choice
-
-
-
-
-
-
-
What is blocking?
-
-Blocking.
-
-
-
Problem: We don't know \( \tau \) or it is too expensive to compute
-
-
Solution: Make a plot of std. dev. as a function of blocksize
-
-
The estimate of std. dev. of correlated data is too low \( \to \) the error will increase with increasing block size until the blocks are uncorrelated, where we reach a plateau
-
-
When the std. dev. stops increasing the blocks are uncorrelated
-
-
-
-
-
-
-
Implementation
-
-
-
-
-
Do a Monte Carlo simulation, storing all samples to file
-
-
Do the statistical analysis on this file, independently of your Monte Carlo program
-
-
Read the file into an array
-
-
Loop over various block sizes
-
-
For each block size \( n_b \), loop over the array in steps of \( n_b \) taking the mean of elements \( i n_b,\ldots,(i+1) n_b \)
-
-
Take the mean and variance of the resulting array
-
-
Write the results for each block size to file for later
- analysis
-
-
-
-
-
-
-
Actual implementation with code, main function
-When the file gets large, it can be useful to write your data in binary mode instead of ascii characters.
-The following python file reads data from file with the output from every Monte Carlo cycle.
-
-The Bootstrap resampling method is also very popular. It is very simple:
-
-
-
Start with your sample of measurements and compute the sample variance and the mean values
-
Then start again but pick in a random way the numbers in the sample and recalculate the mean and the sample variance.
-
Repeat this \( K \) times.
-
-
-
-It can be shown, see the article by Efron
-that it produces the correct standard deviation.
-
-
-This method is very useful for small ensembles of data points.
-
-
-
-
-
Bootstrapping
-Given a set of \( N \) data, assume that we are interested in some
-observable \( \theta \) which may be estimated from that set. This observable can also be for example the result of a fit based on all \( N \) raw data.
-Let us call the value of the observable obtained from the original
-data set \( \hat{\theta} \). One recreates from the sample repeatedly
-other samples by choosing randomly \( N \) data out of the original set.
-This costs essentially nothing, since we just recycle the original data set for the building of new sets.
-
-
-
-
-
Bootstrapping, recipe
-Let us assume we have done this \( K \) times and thus have \( K \) sets of \( N \)
-data values each.
-Of course some values will enter more than once in the new sets. For each of these sets one computes the observable \( \theta \) resulting in values \( \theta_k \) with \( k = 1,...,K \). Then one determines
-
-These are estimators for \( \angle\theta\rangle \) and its variance. They are not unbiased and therefore
-\( \tilde{\theta}\neq\hat{\theta} \) for finite K.
-
-
-The difference is called bias and gives an idea on how far away the result may be from
-the true \( \angle\theta\rangle \). As final result for the observable one quotes \( \angle\theta\rangle = \tilde{\theta} \pm \sigma_{\tilde{\theta}} \) .
-