Contents

Chapter 1
Logistics

1.1 Prerequisites

Computer Science 106(Introd. to Computer Programming) or equivalent
Stat 116 (Probability) or equivalent
Math 113 (Matrix Algebra) or equivalent

1.2 Text

Computational Statistics, Geof Givens and Jennifer Hoeting. Wiley 2005

1.3 Assignments

There will be four homeworks/labs to hand in 40 % and one big project ( final) 60 %

1.4 Projects

Projects

1. Pick a paper from the projects below or on the website on the course webpage at
   http://www-stat.stanford.edu/~susan/courses/s227/.
2. Write a scientific introduction to the paper you have chosen: context and state of the art at the time of the paper (eventually if the paper is old, what has changed in the area since the paper was written).
3. Write an algorithm explaining what the method in the paper does.
4. Simulate data from a known distribution appropriate for the analysis in the paper.
5. Replicate the method in the paper by writing an R program that implements one of the methods implemented in the paper and show it does what the authors say it does on your simulated data.
6. Analyze some real data with your R program.
7. Provide a complete bibliography.
   Bonus: Make a R package with the simulated data, example data and documented programs.

   Hint: the more original the work the higher the grade.

Papers for projects

Gibbs Sampler

Estimating the number of essential genes in a genome by random transposon mutagenesis Natalie J. Blades and Karl W. Broman, technical report, Johns Hopkins.

Coupling from the Past

Exact Sampling with Coupled Markov Chains and Applications to Statistical Mechanics, (1996) Jim Propp and David Wilson
http://www.math.wisc.edu/~propp/sample-short.ps.gz

Sequential Inportance Sampling for Tables

Chen, Y., Diaconis, P., Holmes, S. and Liu, J.S. (2005). Sequential Monte Carlo Methods for Statistical Analysis of Tables. Journal of the American Statistical Association, 100, 109-120.
http://www-stat.stanford.edu/~susan/papers/chenJASA05.pdf
Chen, Y., Dinwoodie, I. H., and Sullivant, S. (2004).
http://www.stat.duke.edu/~yuguo/paper/sis-4.pdf

Mixtures

Bayesian analysis of Poisson mixtures (Peter Green and Valerie Viallefont and Sylvia Richardson)
http://www.stats.bris.ac.uk/~peter/papers/poisjns2.ps
Bayesian Modelling and Inference on Mixtures of Distributions Jean-Michel Marin, Kerrie Mengersen and Christian P. Robert Handbook of Statistics 25, Elsevier
http://eprints.qut.edu.au/archive/00000901/02/Bayesian_Modelling.pdf

Reversible jump Monte Carlo

Reversible jump Markov chain Monte Carlo computation and Bayesian model determination, Biometrika, 82, 711-732 (1995).
http://www.stats.bris.ac.uk/~peter/papers/RJMCMCBka.pdf

Exhaustive Methods

http://www-stat.stanford.edu/~susan/courses/s208/node12.html
P. Diaconis and S. Holmes. Gray codes for randomization procedures. Statistics and Computing, 4:287-302, 1994.

Particle Filtering

Doucet, A., de Freitas, J.F.G. and Gordon, N.J. (2001). Sequential Monte Carlo Methods in Practice. New York: Springer-Verlag.
Doucet, A., Godsill, S.J. and Andrieu, C. (2000). On sequential Monte Carlo sampling methods for Bayesian filtering. Statist. Comput. 10 197-208.
Kitagawa, G. (1996). Monte Carlo filter and smoother for non-Gaussian nonlinear state space models. J. Comput. Graph. Statist. 5 1-25.
Liu, J.S. and Chen, R. (1998). Sequential Monte Carlo methods for dynamic systems. J. Amer. Statist. Assoc. 93 1032-1044.
Liu, J.S. (2001). Monte Carlo Strategies in Scientific Computing. New York: Springer-Verlag.
Pitt, M.K. and Shephard, N. (1999). Filtering via simulation: auxiliary particle filter. J. Amer. Statist. Assoc. 94 590-599.

Data Repository

http://www.ics.uci.edu/~mlearn/MLRepository.html

1.5 Office Hours

Professor Susan Holmes: Wednesday at 1.45, in Sequoia 102.
By email appointment to susan at stat.stanford.edu.
TA: Patrick Perry (patperry at stanford.edu), room 108.
Office hours: Monday 2.30-4.00pm.

1.6 Course Web Site

This will contain a bulletin board, homeworks, course summary, project description list, reading list, links to useful sites with in particular tutorials, software information, etc...

Address:http://www-stat.edu/~susan/courses/s227/

Weekly consultation of the web site will be necessary and expected of all students. Homework and lab. assignments and solutions will be on the course's coursework site.

1.7 Startup computer labs Sequoia 201

There will be 3 sets of labs with exercices on accessing the web and usin R.

Times: TBA

1.8 Outline

The course will contain 3 parts, here is a list containing both subjects covered in class and potential project matter:

1.9 Statistical Methods

Maximum Likelihood, EM algorithm, Missing Data,
Bootstrap, Cross-Validation
Step by step multivariate regression, robust regression,
Non-parametric Regression, Alternate Conditional Expectation, Projection Pursuit,
(Neural Nets)
Principal Components, Correspondence Analysis,
Classification and Regression Trees, Clustering,
Nonparametric Confidence Regions, (convex hulls),
Approximate Counting, Exhaustive Enumeration, Non parametric Tests.

1.10 Algorithms

Matrix manipulations, decompositions, (cholesky, QR,...),
Singular Value Decompositions, Eigenanalysis,
Iterative methods, Gauss-Seidel, Conjugate Gradient,
Downhill Simplex minimization,
Smoothing,
Random Number Generation ,
(Uniform, Normal, Beta, Inverse CDF, Acceptance Rejection, Metropolis),
Markov Chain approach to Combinatorial Optimization (Simulated Annealing), Updating,
Markov Chain Monte Carlo, Gray Codes
Sorting.

1.11 Software

Most methods exist in high level programming environments such as R.

We will learn about several additional R functions/packages such as utils (Sweave), mva,ade4, boot, MCMCpack, Rggobi, rpart, rrcov, modreg,....

Symbolic software will be introduced, if not used much. (maple)

Chapter 2
Homeworks

Hw 1 to hand in before Tuesday, April 19th at 4.30pm Rlab1.pdf
See coursework for solutions.

Hw 2 to hand in before Thursday, April 28th at 4.15pm hw2.pdf

Hw 3 to hand in before Monday, May 9th at 4.15pm Chapter 6, problems 6.4, 6.5 and 6.6.

Hint: If you are having trouble with the problem 6.4, there is a dierent approach using Gibbs sampling that is given in detail here:
http://psblade.ucdavis.edu/BMSBSA/chapter09.ps, this might help you somewhat with a simplification of your calculations.

Hw 4 to hand in before Monday, May 23th at 4.15pm

1. Pick a paper from the projects below or on the website on the cours e webpage at
   http://www-stat.stanford.edu/~susan/courses/s227/.
2. Write a summary of the paper you have chosen in 2 paragraphs.
3. Explain what you plan to do in your final project.
4. Chapter 9, do problem 9.6.
5. Chapter 10, number 10.1 except for (e).
   Hint: Look at the packages quantreg, locfit and lokern for adaptive bandwidth methods.
   Very nice set of notes that will help you (in particular page 11):
   http://www.stat.berkeley.edu/~stark/Teach/S240-s03/Notes/ch8.pdf

Chapter 3
Numerical Analysis for Statisticians

3.1 If you missed some lectures:

Don’t make a habit of it :75 minutes gained will probably cost you the triple to catch up!

I won't be explicitly writing all my lectures as notes on the web, but for the benefit of those who missed a lecture, here is a little of what I said, so you get an idea about at what level the course will be:

3.2 Real numbers: Floating Point Arithmetic

These are normalized so that the first digit of the mantissa is non zero, in the binary mode, ( = 2), this means we don't have to code the first one, sometimes it is left out (phantom of the digit).

If x < ^Emin-1 (smallest number representable in machine)
we have underflow,
if x > ^Emax(1 - ^-t) overflow.

Otherwise x is in the representable range, maybe we won't be able to represent it exactly. Macine Accuracy:
Smallest fp number added to 1 that gives an answer dierent from 1. (eps in matlab)

Now if we can't represent x without representation error, we'll call x*

Taylor series show how much we are going to be o.

When error growth is linear, the sensitivity is small: f is said to be well-conditionned. relative to the data x and the problem S: is the constant that bounds the relative error.

Relative to the norm, the problem and that particular data x.

Worse still, we usually don't have f, but an approximation of it. this induces a second type of error: generated error, also called truncation error because we stop computing before we have the exact right answer(for an integral computed by trapezoids, or any iterative method for computing functions by series, ...) for example.

An unstable algorithm is one where the roundo error gets mixed into the calculation at an early stage and is successively magnified until it swamps the true answer.

Example: Want to compute the powers of the golden mean:

Error analysis:

• Forward: take the errors and keep track of all the fp operations and how they make the errors evolve.
• Backward analysis: suppose you wanted to compute the root of an equation: F(x) = 0, you only get x*, then see how much you are o by computing F(x*) if it's large, your method is not so good.
• Try to compare on a simpler problem where you know the answer.

Moment computations, (see Pages 612-613 of Numerical Recipes)

A little example of how mistakes tote up on a total.

Vectors These are the basic building blocks of the new generation of languages, unlike C or fortran whose are fp numbers, the gains are important in writing all computations vectorially.

3.3 Matrices, and their decompositions

Matrices Notations, operations, partitionning, careful with the defaults from one system to another.

The Botany of Matrices Diagonal, Symmetric, Banded, Hessenberg, Circulant, Hankel, Triangular, Sparse, Orthogonal, semi-positive definite,

My notations for: transpose, inner product,norm, matrix norms,

3.4 The linear equation solving problem

Suppose we have a linear system of equtions, typified by the matrix of coefficients A, unknowns, x and the right hand sides, the vector b.

3.4.1 Square systems

Singularity, Degenerate Systems

If the rows are linear combinations of one another, or the columns linear combinations or one another then even if there seems to be the same number of unknowns as equations, the system is in fact equivalent to a singular set of equations.

Close singularity is just as bad because the roundo could take us there, we will have to measure this.

Single precision is usually enough for 30-50 equations.

Double precision makes it possible to solve for a few hundred.

This does NOT hold for pathological systems we will look into these and how to recognize them when we talk about conditionning.

A full square system

:

Suppose the matrix of coefficients is a square matrix, (of full rank also, as we saw, non-degenerate) then the equation looks like this:

Formally you would say : invert it, in fact a computational mistake, except if you have multiple right hand sides.

Doing what one would first think of doing if you had to solve a system by hand.

The easiest system: triangular.

Bringing to triangular shape.

Pivotting. This is NOT to solve ill-conditionning, but unlucky diagonal elements that are very small make the algorithm unstable.

LU decomposition, with a permutation matrix (prepivotting) is what is best to do.

Flops: Of an order (the same as the cost of the multiplication of two matrices, not so bad).

Condition number for the inversion problem

Try a small perturbation of b to b^* = b + b, call the solution of the perturbed system: x^* = x + x,

The condition number for the data A and the linear solving problem is thus: (depends on the norm, not on the algorithm). For the two norm we have a better explanation of how to compute this, we need the following to do it:

Singular Value Decomposition

This is the most important decomposition in statistics.

We will study it several times, here is my first introduction: a little made-up example with matlab to start with:

>> u=[3 1 -1 2]'   u =        3        1       -1        2   >> v=(1:4)'   v =        1        2        3        4   >> X=u*v'   X =        3     6     9    12        1     2     3     4       -1    -2    -3    -4        2     4     6     8   >>svd(X)   ans =      21.2132       0.0000            0            0     >> E=10^(-3)*randn(4);   >> XE=X+E   XE =         3.0012    5.9993    9.0003   12.0012       1.0006    2.0017    3.0009    3.9994      -0.9999   -1.9999   -3.0014   -3.9994       2.0004    4.0018    5.9993    7.9996     >> svd(XE)   ans =      21.2143       0.0028       0.0019       0.0005   >> cond(X)   Condition is infinite   ans =      Inf   >> cond(XE)   ans =        4.2567e+04   %An example you can't see with your bare eyes:   >> X2=u'*v   X2 =       0.2468    0.2499    0.0166    0.2487       0.3225    0.3266    0.0218    0.3250       0.1495    0.1514    0.0101    0.1506       0.4367    0.4423    0.0295    0.4401     >> svd(X2)   ans =         1.0730       0.0000            0            0   >> [U,S,V]=svd(X2);   >> U*S*V'   ans =       0.2468    0.2499    0.0166    0.2487       0.3225    0.3266    0.0218    0.3250       0.1495    0.1514    0.0101    0.1506       0.4367    0.4423    0.0295    0.4401   >> 10000*(X2-U*S*V')   ans =       1.0e-11 *      -0.0278    0.0833   -0.0035    0.0555      -0.0555    0.0555         0    0.0555            0    0.0555    0.0017    0.0278      -0.0555    0.1665         0    0.1110

>> A=rand(4)   A =       0.5045    0.4940    0.0737    0.9138       0.5163    0.2661    0.5007    0.5297       0.3190    0.0907    0.3841    0.4644       0.9866    0.9478    0.2771    0.9410   >> flops(0)   >> [L,U,P]=lu(A)   L =        1.0000         0         0         0       0.5233    1.0000         0         0       0.5114   -0.0406    1.0000         0       0.3234    0.9388    0.7363    1.0000   U =        0.9866    0.9478    0.2771    0.9410            0   -0.2298    0.3557    0.0373            0         0   -0.0535    0.4342            0         0         0   -0.1945   >> flops   ans =       34   >> P*A   ans =       0.9866    0.9478    0.2771    0.9410       0.5163    0.2661    0.5007    0.5297       0.5045    0.4940    0.0737    0.9138       0.3190    0.0907    0.3841    0.4644   >> L*U   ans =        0.9866    0.9478    0.2771    0.9410       0.5163    0.2661    0.5007    0.5297       0.5045    0.4940    0.0737    0.9138       0.3190    0.0907    0.3841    0.4644   >> P*A-L*U   ans =    1.0e-15 *            0         0         0         0      -0.1110         0         0         0            0         0         0         0      -0.0555   -0.0139         0         0

Here is what we need to remember:

Actually the singular values are the square roots of the eigenvalues of X'X.

The condition number in ² norm has the simple interpretation as a function of the singular values, ration of the largest to the smallest:

Read 2.6 of NR (pages 61-71).

Tridiagonal

Of an order of N steps are sufficient, see NR 51.

Iterative Improvement

This enables, a the price of storing A, to improve the solution by using the residual. Do the following in double precision:

Cholesky Decomposition

For a symmetric positive definite matrix A the LU decomposition can be rewritten: with D a diagonal matrix with only positive entries.

3.4.2 Rectangular system: QR decomposition and least squares

When there are more equations than variables, one can't hope for an exact solution, the system is overdetermined.

We try to minimize the error which is usually measured as the norm

Important remark: For any orthogonal Q, this norm is unchanged by trnasformation by Q, so

3.4.3 Examples

Least Squares in Two dierent ways:

3.4.4 Givens Rotations

What are Given's rotations good for?

You can use them to zero out individual isolated elements in any matrix, without changing any of the norms of the vectors, these transformations are orthogonal.

They are based on simple 2 × 2 ones of the form:

You can choose c and s so that a certain vector is orthogonally transformed (actually, rotated) into a vector y = (y₁,y₂) whose second coordinate is zero.

This is of course easy by solving the equations

The examples that follow show two ways of doing this:

Now in order to extend this to a general context you just notice that if you fill out the rest of the matrix with 1's on the diagonal, it will not change the other coordinates thus:

3.4.5 R examples handout

3.4.6 Backsubstitution

3.5 Householder Functions

3.5.1 Householder Multiplication-Row

3.6 QR decomposition in R

3.7 Singular value decompostion and pc regression

3.7.1 Ridge Regression

Ridge regression in R

(see Faraways chapter 9 )

Example of a matlab ridge regression function:

3.7.2 Principal Components Regression

If we decompose the centered and rescaled matrix

3.7.3 Principal Components

A quick summary, all is detailed below:

• X centered 1_n X = 0,
• SVD best approximation to x,
• Cols of V new variables,
• Principal Components C = US C^'C = S² C = XV ,
• Principal axes Z = U^'X = SV ^',
• Supplementary points x_k^.V = x_k^.Z^'S^-1,
• Distance between two points,
• Transition Formulae Z = S^-1C^'X C = XZ^'S^-1,
• Scale problem solution: Q changes the metric,
• Generalized PCA 2nd generalized svd,
• B = D X Q , U = D^- E, V = D^- A,
• Special case: correspondence analysis,
• N V = , r = N ^. 1_p, c = N^'. 1 _n, X = rc^',
• (-) ,
• Number of Components. PRESS(M) = ^{n p2}.

Maximization of u^'Au with u^'Mu = 1

• for u eigenvector of M^-1A,
• F.

1) X centered, all points (observations) same weight.

The columns of V are the directions along which the variances are maximal.

Definition: Principal components are the coordinates of the observations on the basis of the new variables (namely the columns of V ) and they are the rows of G = XV = US. The components are orthogonal and their lengths are the singular values C^'C = SU^'US = S².

In the same way the principal axes are defined as the rows of the matrix Z = U^'X = SV ^'. Coordinates of the variables on the basis of columns of U.

The new model will be:

Example with R
First rescale to have unit variance:

Then do the singular value decompositions:

3.8 More than I said on the Singular Value Decomposition

Singular Value Decomposition

Of numerical and statistical interest is the decomposition of X known as the s.v.d.

Up to now we have shown that X could be written X = UR, U orthogonal n × n R = p × p upper trianglar. The transformation of y by U^' is c = U^'y + , and the solution to the U's problem is provided by s such that

Columns of U are left-singular vectors, columns of V right-singular vectors.

Remark: When X is symmetric positive semi-definite V = V and the singular values of X are its eigenvalues. On the other hand, X^'X = V S²V ^'. The eigenvalues of X^'X are the square of the singular values of X (and XX^' U).

Statistically: E(Y ) = X = USV ^' = US = XV . is the vector of regression coefficients on a new collection of p orthogonal regressors and these coefficients are rotated versions of the original p regression coefficients.

When is S_i small?

Courant-Fischer Theorem:
X_n×p, n > p and singular values s₁ > s₂> s_p > 0.

s_k is then the largest possible value for the length of X.

V _i is the vector for which the corresponding linear combination of the columns of X has greatest sum of squares. (d₁ is its square root).

If the column means have been removed from X then this corresponds to the linear combination that has the largest variance.

= (₁,,_p) parameter of interest
log likelihood () achieves its maximum value for () = = 0_p, where

The Fisher information or information matrix for is defined to be

e = Y - X r = rank(X) k₂(X) = d₁/d_r.

TSS = SSM + RSS
Total = Model + Residual
R² = SSM/TSS = 1 - RSS/TSS.

Theorem: = max < .

estimate for the perturbed problem of reg. coe.

ê estimate of residually perturbed problem of reg. coe.

For R² < 1

3.8.1 Generalized Inverses

What if X is not full rank?
This will often happen when coding for analysis of variance.

X^'X singular, the vector of coefficients is not estimable although some linear combinations of 's are estimable. (Those for which I(c^' : Y ) = c^'X^'Xc are nonzero.

c^' estimable if there exists a vector such that

(b) equivalent to c^'(X^'X)^-X^' = c^'.

+1

X⁺ is the Moore-Penrose pseudo-inverse of X and if X _n×p then X_p×n⁺ is defined as the unique matrix for which:

• (XX⁺) symmetric.
• (X⁺X) symmetric.
• XX⁺X = X.
• X⁺XX⁺ = X⁺.

The solution = X⁺Y is equivalent to = (X^'X)⁺X^'Y .

3.9 Principal Components

3.9.1 Eigenvalue Analysis

Summary:

• Motivation: PC are new variables, uncorrelated built from the old ones.
• Coefficients are obtained through eigenvalues of variance-covariance or correlation matrix.
• Eigenvalues represent variance explained.
• Also useful discriminant analysis, canonical variate analysis.
• Generalized eigenvalue problem det(A - B) AB^-1 B^- A B^- .
• Statistics symmetrical eigenvalue problems.
• A and Z^'AZ same eigenvalues Z tranforms A's structure.
• A = LR RL - L^'AL similarity.
• Power method: ₁ simple v₁,x₀ non v₂ X_i = A.
• x_n = _nA_?ⁿx ₀ = w_iiⁿ, || x_n || ₁.
• Advantageous for A, easy to compute, only provides v₁,₁.
• X₀ orthogonal columns, re-orthogonalizing at each step.
• T = AX_i-1 T = X_iR_i X_i [v₁v_k],R_i [_1k].
• Rate _k+1/_k.
• QR-algorithm A_i+1 = Q_i^'A _iQ_i where A_i = Q_iR_i.
• Improvement 1: start with A₀ = U^'AU tridiagonal (Householder).
• Tridiagonal form is not lost at each step, subdiagonal elmination through Givens to give the QR decomposition.
• Improvement 2: convergence depends upon _p+1/_p. A - I and A same eigenvectors and eigenvalues dier by .
• Choose so that (_k+1/_k) small ( close to _k+1).
• _k+1 not known, Q_pp good or eigenvalue of lower block of A.
• A_i+1 = R_iQ_i + I.
• Cubic convergence.
• Implicit lemma: A sym. non-singular B = Q^'AQ, Q, B positive o-diagonal elements then Q and B are determined as soon as the first col. of Q is.
• A_i+1 = R_iQ_i + _iI = Q^'(A ₂ - _iI)Q_i + _iI = Q_i^'A _iQ_i.
• R_k = Q_k^'(A _k - mu_kI).
• Q_k = J₁J_p,J_k = G_k+1k.
• First row of Q_k^' is that of J ₁^' because only one to touch first row.

Example of the Power method with Matlab

Chapter 4
EM algorithm

Motivation:

A method for finding maximum likelihood estimates,
-either in presence of missing data.
-or when the model can be simplified by adding `latent parameters'.

Some references:

Hartley, 1958, Biometrics, 174-194, is a good starting place because he does simple examples.

Apparently the method was discovered by Oscar Rothaus while he was at IDA, with Baum, who then published several papers of which an example is Baum, Petrie, Soules and Weiss, (1970), Ann. Math. Stat., 164-171 who use the idea for speech recognition.

The main reference that people cite is Dempster, Laird et al. (1977), JRSS B - they coined the name EM.

It is hard to motivate with a `trivial example' because the the likelihood is easy to dierentiate.

Most of the elementary examples come from the multinomial likelihoods.

A very recent reference is the book called: The EM algorithm and Extensions, by G.J. McLachlan and T. Krishnan,(1997), Wiley.

Here is a first motivating example

Example 1:
Multinomial counts of animals into 4 categories

Then we can consider that this is in fact a five-category problem with first two categories being split.

Call = (1 - )², then the density is proportional to

The solutions are:

Now, suppose we didn't have the equation solver available and that we want to nd the solution using the EM algorithm.

Split the rst cell into two cells with respective counts y₁₁ + y₁₂ = y₁, then the new likelihood is proportional to:

The maximisation step doesn't need any more computations, as now we have values for the unobsevrables that we just plug into

4.0.2 Implementation

Here is the little R function that I used to check the computations:

Results:

4.0.3 Estimating Mixture Proportions

Suppose the density is of the form:

,

,

If we observed the z's, the mle of _i would be _i = .

Take the new, complete data vector to be x = (y,z).

The E-step:

4.0.4 EM for exponential families

Suppose the complete density is of the form:

As an application, you can consider any multinomial as an exponential family, (see Lehmann, Theory of Point Estimation, page 28), so that they all enter this case.

A particular example was the one I did the rst day where the new parameter space had a natural interpretation, that of the gene counting example, where the rst mutinomial is the observed one of the frequencies of the observed phenotypes. y = (n_A,n_B,n_AB,n_O), the parameter of interest is the probability of each of the genotypes A,B,O = (p,q,r), we can usethe augmented model: x = (n_AA,n_BA,n_BB,n_OO)

Here is a link to some very serious papers about tomography: http://www.stats.bris.ac.uk/pub/reports/TOMO.

And here is a link to a nice application of the EM to Bayesian self-organizing map simulations: http://www.aist.go.jp/NIBH/ b0616/Lab/BSOM1/

Chapter 5
Monte Carlo Simulations

5.1 Monte Carlo methods

There is not necessarily a random component in the original problem that one wants to solve,usually a problem for which there is no analytical solution. An unknown parameter (deterministic) is expressed as a parameter of some random distribution, that is then simulated. The oldest well-known example is that of the stimation of by Buon, in his needle on the oorboards expeiment, where supposing a needle of the same length as the width between cracks we have:

1. The crude, or mean-value Monte Carlo method thus proposes to generate B numbers uniformly from (0, 1) and take their average: to estimate ,
2. The hit-or-miss Monte Carlo method generates random points in a bounded rectangle and counts the number of 'hits' or points that are in the region whose area we want to evaluate.

Which estimate is better?
This is similar to comparing statistical estimators in general. There are certain desirable properties that we want estimators to have, consistency which ensures that as the sample size increases the estimates converges to the right answer is ensured here by the properties of Riemann sums. Other properties of interest are:

• Unbiasedness.
• Minimal Variance.

Both the above methods are unbiased, that is when repeated many times their average values are centred in the actual value .

More precisely, the variance of crude Monte Carlo is

5.1.1 Antithetic Resampling

Suppose we have two random variables that provide estimators for , X and Y , that they have the same variance but that they are negatively correlated, then (X + Y ) will provide a better estimate for because it's variance will be smaller.

This the idea in antithetic resampling (see Hall, 1989). Suppose we are interested with a real-valued parameter, and that we have ordered our original samplex₁ < x₂x_n, for each resample X^* = {x _j1,x_j2,...,x_jn} and statistic s(X^*) we associate X^** by taking a special permutation of the j_i's that will make cov(X^**),s(X^*)) < 0, and as small as possible. If the statistic is a smooth function of mean for instance, then the 'reversal permutation' that maps 1 into n, 2 into n - 1, etc... is the best, the small sample values are transformed into the larger observations, and the average of these two estimates will give an estimate with smaller variance.

5.1.2 Importance Sampling

This is often used in simulation, and is a method to work around the small area problem. If we want to evaluate tail probabilities, or very small areas, we may have very few hits of our random number generator in that area. However we can modify the random number generator, make that area more likely as long as we take that into account we we do the summation. Importance sampling is based on the equalities:

5.2 Random Number Generation- Uniform[0, 1]

Try the following web site pages: http://www.maa.org/mathland/mathland_4_22.html http://stat.fsu.edu/~geo/diehard.html

And read chapter 7 of Numerical Recipes.

Motivation:
Random numbers are used:
1- To solve problems that are too complex (stochastic algorithms for NP complete problems)

2- For optimization (simulated annealing, genetic algorithms).

3- System Simulation (physics and biology, econometrics...)

4- Non-parametric tests, bootstrap, permutation tests.

5- Quadrature methods using Monte Carlo.

Solving a deterministic problem by rewriting it as equivalent to a problem with a stochastic component.

Little Example as a reminder of what Monte-Carlo is:

Most methods, even if they need non-uniform random variates, as in simulation, still use uniform random variables as their starting points.

So we will spend some time explaining the basics of random numbers.

The machine being deterministic, there is denitely a problem, we will never get purely random numbers, they will only look so.

(We already saw that there are some pretty bad ones, when we did the multidimensional visualization of the IBM randu generated numbers with xgobi).

Also beware of srand, the basic unix one is bad.

5.2.1 Congruential Methods

The way most of the algorithms work is they use an initial seed and a simple recursion to feed out the next in sequence, which will also serve as a seed .. etc.

As soon as one hits a number that has already been used we will be having a sequence that was already seen. -Not so good-

Period of the sequence= minimal length of a subsequence that by repetition will form the whole sequence.

Period of the generator will be the largest period of any seed for a given generator.

We need to choose the parameters so as to have a maximal period.

Proposition:
The congruential generator has a period of m i:

1. c is relatively prime to m.
2. mod 1(mod p) for any prime factor p of m
3. mod 1(mod 4) if 4 is a factor of m.

Actual choices:
They have to take into account both theoretical and practical aspects:

1. Repeatability
2. Portability
3. Feasability
4. Non-predictability

5.2.2 Other `Better?' Methods

Substract with Borrow

Instead of a single seed, the random number generator has a state: (35 words of memory), 32 fp numbers z between 0 and 1. The others make up an integer index : i, which varies between 1 and 32, and integer j and a 'borrow' ag.

j is used to initialize.

z_i = fp(z_i) xor j b is the borrow, small or zero. Complexity Justified Methods
Random number generation is akin to cryptography, knowing that factorising a number into its prime factors is a very HARD problem (computationally) allows one to think up various `multiplicative schemes', for instance:

and back to tables but on CDs now (new Numerical Recipes).

5.2.3 Random Number Generation- other than Uniform

Luc Devroye, Non-Uniform Random Number Generation, Springer Verlag.
A very useful home page is Luc's with links to places where you can download random number generator packgages etc... http://jeff.cs.mcgill.ca/~luc/ Web pages with good links to papers on random number generation, and encryption packages:
http://www.cs.berkeley.edu/~daw/netscape-randomness.html

5.2.4 Inversion Method

Based on the inverse CDF property :
X = F^-1(U) ~ F where U is uniform:

This means that when we have a closed form formula for the inverse of the cdf which is not too expensive, we can use a random uniform and do the transformation that makes it have the correct density.

Example 1: Exponential

Example 2: Correlated Random Variables
F-1(U) and G-1(U) have maximum positive correlation and have the distributions F and G, for negative correlation just take G-1(-U).

Example 3: Maxima
(Fⁿ)^-1(U)

If closed forms are not available, but a good numerical solver is, then you can use that.

5.2.5 Rejection Methods

Theorem: Suppose f is a non zero density on [a,b] and zero outside. (a,b may be innite). Suppose we know a function M(x) which majorizes f: M(x) > f(x) on [a,b] and take:

Algorithm:

1. Generate T with density m.
2. Generate U uniform on [0, 1], independently of T

This is why it works:

We need more than one uniform to generate one random variable with the right distribution.

How many are needed? This will depend on c = _a^bM(t)dt.

The number of necessary trials before accepting will be a geometric with p = .

Example 1:
Suppose f is a bounded density on a compact support, f(x) < M take c = M(b - a),
Repeat the folowing:

• Generate two independent uniforms, U and V on [0, 1].
• X u + (b - a)V

Until UM < f(X)
Return X.

This is usually very inecient, it is much better to take a bounding function that `ts'.

Example 2:
Simulating a normal using an exponential as a bounding function: we use the fact that the absolute value of a normal has density:

5.2.6 Table Lookup

If the probabilities are rational numbers and the variable discrete, or well approximated by a discrete distribution you can use table lookup, which is very expensive to set-up but can be useful if you a very a very large number of random variables to generate:

Algorithm:

• Given p₀,p₁,...,p_k = , generate a table of length M with n_k entries equal to k.
• Generate a uniform U.
• Return A[floor(M * U)]

Examples : histogram data.

5.3 Specialized Methods

5.3.1 Polar methods for the Normal

Box - Muller
Fact:
If X and Y are independent normals then their polar coordiantes R = , = arctan(Y/X) are independent, R² ~ expon(2), ~ U(0, 2).

Modified to short-cut cosines and sines:
Generate a point uniformly in the disk(1): (V ₁,V ₂).
Call S = V ₁² + V ₂²,

Return X = V ₁,Y = V ₂

5.4 Importance Sampling

See notes at:
http://www.ma.utexas.edu/documentation/nr/bookcpdf/c7-8.pdf (Chapter 7 of Numerical Recipes)
or a nice tutorial here: http://ib.berkeley.edu/labs/slatkin/eriq/classes/guest_lect/mc_lecture_notes.pdf

Chapter 6
The Bootstrap, Permutation Tests, Simulation

6.1 Motivation

We want to compute bias, variability, condence intervals for estimates for which the theory doesn't provide closed form solutions.

Examples :

1. The parameters in non parametric regressions (smooths, ppreg,...)
2. Correlation coecients for samples that don't necessarily come from a Normal distribution.
3. Testing when we don't know the distribution of the test statistic under the null hypothesis.

Suppose we are interested in the estimation of an unknown parameter (this doesn't mean that we are in a parametric context). The two main questions asked at that stage are :

• Question 1 What estimator should be used ?
• Question 2 Having chosen an estimator, how accurate is it ?

The second question has to be answered through information on the distribution or at least the variance of the estimator.

Of course there are answers in very simple contexts: for instance when the parameter of interest is the mean then the estimator has a known standard deviation : the estimated standard error noted sometimes

In maximum likelihood theory the question 1 is answered through using the mle and then the question 2 can be answered with an approximate standard error of

The bootstrap is a more general way to answer question 2 , with the following aspects:

• Less or no parametric modelling.
• More computation. (a factor 100 to 1000)
• Automatic, whatever the situation (can be complex).

If we had several samples from the unknown (true) distribution F then we could consider the variations of the estimator :

6.2 The bootstrap : the univariate context

6.3 In practice

From an original sample

If we are interested in the behaviour of a random variable = (X,F), then we can consider the sequence of B new values obtained through computation of B new bootstrap samples.

Practically speaking this will need generatation of an integer between 1 and n, each of these integers having the same probability through a function denoted by RANINT(1,n) hereof.

Here is an example of a line of matlab that does just that:

If we use S we won't need to generate the new observations one by one, the following command generates a n-vector with replacement in the vector of indices (1...n).

An approximation of the distribution of the estimate = (X_n,F) is provided by the distribution of

1. For b=1 to B do : /* B is the number of bootstrap samples */
   1. For i=1 to n do : /* Generation of Sample X_b^* */
      1. OBS=RANINT(1,n)
      2. Add observation number OBS to Sample X_b^*
   2. Compute _b^* = (X _b^*)

An example data set is LS (or laws) : the law school data consisting in 15 bivariate data points, each corresponding to an American law school as analysed in Efron 1982:

In order to repeat a large number of times the same instructions needed to construct B data sets one usually has to use a looping facility, actually matlab and S are easier to program in this circumstance because one can use ordinary loops and use owner-created functions, as well as trying to benet form the vectorized commands.

6.4 S examples for simple bootstraps

6.5 Parametric bootstrap

If one has a parametric distribution F which is well dened except for the pamameter , instead of drawing with replacement from the original sample, one can draw from the distribution F, where is the estimation obtained from the original sample.

6.6 Smoothed Bootstrap

An intermediate solution between parametric and nonparametric bootsrapping is the smoothed bootstrap.

Instead of resampling directly from the empirical distribution F_n, one smoothes it out rst then the smoothed version is used to generate the new samples.

If a kernel smoother is used, for instance a Gaussian kernel, then generation from the smoothed distribution _n is easy since it is sucient to resample with replacement and then perturb the sampled point by adding a small gaussian random variable (Silverman B.W. & Young G.A.)

6.7 Dierent Uses for bootstrapping

6.7.1 Quality of Estimates

Bias :

Variance :

Confidence Intervals : The percentile method As before we denote by G_n^*(t) = P _Fn the bootstrap distribution of ^*, approximated by

The percentile method consists in taking the 1 - 2 condence interval for as being

:

1. Input the level = 2 of the confidence interval
2. For b=1 to B do /* B is the number of bootstrap samples*/
3. For i=1 to n do /* Generation of Sample X_b */
   1. OBS=runif(1,n)
   2. Add observation number OBS to the Sample X_b
4. Compute _b^* = (X _b)
5. Combine the _b^* = (X _b) into a new data set THETA
6. Compute the and (1 - )’th quantiles for THETA.

6.8 Enhancements

As can be predicted, the bigger B, the number of bootstrap resamples, the better the approximations are, thus much recent work has been devoted to making the most of the computations :

• By choosing better quantities to bootstrap than the original estimate, through various transforms some of which require a double level bootstrap. (such as those in sections 2.1, 2.2 and 2.5 )
• By using techniques inspired from Monte Carlo methods : importance sampling, antithetic sampling. balanced sampling.
• By transforming the results obtained after an ordinary bootstrap process : bias correction Efron 1982, accelerated bias correction Efron 1987 for instance.

Preliminaries : Pivotal quantities
Construction of condence intervals is based (ideally) on a pivotal quantity R_n, that is a function of the sample and the parameter whose distribution is independant of the parameter , the sample, or any other unknown parameter. Thus for instance no knowledge is needed about the parameters and ² of a normal variable X ~ N(,²) to construct condence intervals for both of these parameters, because one has two pivotal quantities available :

Let c_n() denote the largest (1 - )th quantile of the distribution of the pivot R_n, and the parameter space, set of all possible values for then

(1-/2)

/2

1-/2

Such an ideal situation is of course rare and one has to do with approximate pivots, (sometimes called roots). Therefore errors will be made on the level (1 -) of the condence sets. Level error is sensitive to how far R_n is from being a pivot.

We will present several approaches for nding roots One approach is to studentize the root that is divide it by an estimate of its standard deviation. Another is to use a cdf transform constructed from a rst set of bootstrap simulations as explained below.

6.9 Bootstrap-t : Studentizing

Suppose that we are interested in a parameter , and that is an estimate based on the n-sample X_n.

We will suppose that the asymptotic variance of is ²/n with a corresponding estimate ²/n.

The studentized bootstrap or bootstrap-t method (Efron 1982, Hall 1988) is based on approximating the distribution function H_F (x) of n (- )/ by H_Fn(x), distribution of n (- )/^*.

Denote by c_n,^* the th quantile of the distribution of n (- ),

1. Compute the original estimate of : THETA(0)
2. For b=1 to B do :
   1. Generate a sample X_b^*
   2. Compute the estimate of based on X_b^*: THETA(b)
   3. Use X_b^* as original data and compute an estimate of the standard deviation of : SIGMA(b).
   4. Compute R(X_b^*)=(THETA(b)-THETA(0))/ SIGMA(b)
3. Compute CNA and CNCA the and (1 - )’th quantiles for the sample formed by all the R’s.
4. LO=THETA(0)-SQRT(N)*CNCA*SIGMA(0)
5. HI=THETA(0)-SQRT(N)*CNA*SIGMA(0)

6.10 R functions

6.11 Preliminary Transformations of the Estimate

One of the drawbacks of the above-dened method is that the variance estimate can tend to be quite erratic. In the particular case of bootstrapping a correlation coecient one sometimes ends up with condence intervals that contain [-1,1], which is rather disappointing.

However experiences have shown that preliminary transformation of such as Fishers z transform (=arctanh ), followed by the Bootstrap-t and then retransformation of the condence interval back to the original scale leads to much more dependable condence intervals. Unfortunately such a handy transform is not always available but R. Tibshirani, 1988 has suggested an algorithm for creating automatically dened variance stabilizing transformations.

6.12 Pre-pivoting or Double Bootstrap

Based on the inverse CDF property :

i.e. that if X is distributed according to F then F^-1(X) is uniform.

Denote by H_n = H_n(.,F) the c.d.f. of our root R_n() which is supposed to converge weakly to the c.d.f. H = H(.,F) as n increases.

As we have seen the idea behind the bootstrap is to replace the unknown distribution function F by its empirical counterpart F_n, thus using _n = _n(.,F_n) as an estimate for the unknown cumulative distribution. Once validation of the procedure has been insured through verifying that

:

1. For b=1 to B do
   1. Generate sample Y_b^* a bootstrap sample of size n from F _n
   2. Compute R(b)=R_n(Y_b^*, )
   3. Initialize Z(b)=0
   4. For k=1..K do
      1. Generate sample X_k^* a bootstrap sample of size n from F _n,b^*
      2. Compute R_n(X_k^*,)
      3. If R_n(X_k^*,) < R_n(Y_b^*, ) add one to Z(b)
   5. Divide Z(b) by K.
2. Compute CN1 :the (1 - )’th quantile of the Z’s
3. Compute the CN1’th quantile of the R’s.

6.13 Balanced Bootstrap : saving computations

Davison, Hinkley and Schechtmann, 1986 introduced the idea that one could reduce the amount of simulations (=B) necessary to attain a given precision by using each of the sample observations exactly equally often.

The algorithm proposed by the authors was however rather expensive timewise ( nB log nB). We have a more economical one implemented proposed by Gleason .

:

1. Form the list L =
2. For i = n * B to 2 do
   1. j = runif(1,i)
   2. Exchange L_j and L_i

6.14 Antithetic Sampling : more on saving computations

The idea comes from current monte-carlo techniques, and is detailed in unpublished work by Hall. Like other enhancements presented here, its objective is to reduce the estimate's variance.

Let ₁ and ₂ be two dierent estimates of an unknown parameter , if ₁ and ₂ can be chosen so as to be negatively correlated, then = ( ₁ + ₂) has half the variance of ₁ and ₂.

The bootstrap version of this is based on an antithetic permutation of the sample. B ordinary resampling operations then supplies the equivalent of 2B resampling operations.

6.15 Bootstrap and Randomization Tests

6.15.1 An example : Infants walking

(from Daniel W.W., p.216) :
Researchers wished to know if they could conclude that two populations of infants dier with respect to the mean age at which they walked alone. The following data (ages in months) were collected:
Sample from population A:
9.5, 10.5, 9.0, 9.75, 10.0, 13.0, 10.0, 13.5, 10.0, 9.5, 10.0, 9.75
Sample from population B:
12.5, 9.5, 13.5, 13.75, 12.5, 9.5, 12.0, 13.5, 12.0, 12.0
What should we conclude ? ( = 0.05)
The classical approach : a t-test

The software output from script in the appendix shows :
The mean for population A is : 10.375
The mean for population B is : 12.208
The dierence is : 1.833
As the variances can be considered as equal, a classical test under the normal assumption gives a t-value of -3.1561 indicating a signicant dierence because the P-value is 0.0046.

Assumptions :

1. Random sampling of observations from the population
2. Equal population standard deviations for males and females
3. Normal distributions for test scores within groups

A permutation test approach Based on the fact that if there is no dierence between the two populations then the result will be compatible to allocation at random of each observation to one of two groups (shuing). Any statistical measure of the dierence could be used here : mean dierence, studentised (t) dierence, or even a median dierence.

1. Compute the mean difference between samples A and B: MAB(0).
2. For p=1 to P do
   1. Randomly allocate 12 observations to sample A and 10 to sample B
   2. Compute the mean difference MAB(p)
3. Is MAB(0) a typical value from the randomization distribution MAB(p),p=1..P?
   (Compute #{MAB(p)<MAB(0)}/P and #{MAB(p)<MAB(0)}/P

A bootstrap approach Following the previous denitions of the bootstrap a possible bootstrap test could run as follows:

1. Compute MAB=mean difference between the two original samples
2. For b=1 to B=500 do :
   1. Draw a bootstrap sample from population A.
   2. Draw a bootstrap sample from population B.
   3. Compute a bootstrap estimate of the mean difference :MD(b).
3. Construct the bootstrap distribution of MD, eventually smoothed and recentred at 0.
   Situate MAB(0) with respect to this distribution.

6.15.2 More generally

Approximate randomisation tests are frequently used to test independance of two sets of variables (which is the null hypothesis). Estimation of the distribution under this hypothesis is obtained by shuing one variable (or set of variables) with respect to the other variable(s). (see Manly, 1990, for an easy overview).

However it is important to note that there is no assumption on the fact that the observations form a random sample from some population. There is no inference that can be carried automatically over to a hypothetic population.

On the contrary the assumption underlying the bootstrap is that the observations form a random sample whose distribution approximates that of the population. Thus the conclusions drawn from a bootstrap test can be used to infer on the population.

Noreen 1989 suggests using a combination of these two tests when inference is needed for testing whether two variables are stochastically independant in the population. However one must keep in mind that at the same time one tests whether the marginal distributions of the variables in the sample satisfactorially estimate those of the variables in the population.

6.16 Finding out how good a parametric bootstrap procedure is through a simulation experience :

Confidence Intervals for the expectation of an exponential distribution

1. Choose n
2. Choose an expectation for X :
3. For b=1 to B do :
   1. Create a sample X_b^* = X ₁^b,...,X _n^b having an exponential distribution
      (P = 1 - e^-d)
      (In fact we use U₁,U₂,...,U_n uniformly distributed on
      and then compute X_i = log (1 - U_i))
   2. Compute a bootstrap 0.95 confidence interval for
   3. Count c_b = 1 if is in the interval and c_b = 0 otherwise.
4. P = ₁^Sc _i

For B large enough, P provides the probability that is in the bootstrap condence interval.

In order for the experience to be really useful, one ought to retry with dierent sample sizes and dierent expectations, one could also retry with all sorts of dierent distributions.

Using the Bootstrap to build Condence Intervals

We call the unknown parameter 1 and our estimate . Suppose that we had an ideal (unrealistic) situation in which we knew the distribution of -1 , we will be interested especially in its quantiles : denote the quantile by and the 1 - quantile by . By denition we have: P( - 1 < ) = P( - 1 < ) = 1 -

It is by writing this that you can build the correct intervals.

Now I am going to compare an answer to an estimation probalem both through maximum likelihood and through bootstrapping. First let us remind ourselves of notations:

6.17 Maximum Likelihood Estimation

X₁,X₂,X₃,...X_n have joint density denoted

Given observed values X₁ = x₁,X₂ = x₂,...,X_n = x_n, the liklihood of is the function

If the distribution is discrete, f will be the frequency distribution function.

In words: lik()=probability of observing the given data as a function of . Denition:
The maximum likelihood estimate (mle) of is that value of that maximises lik(): it is the value that makes the observed data the \most probable".

If the X_i are iid, then the likelihood simplies to

Rather than maximising this product which can be quite tedious, we often use the fact that the logarithm is an increasing function so it will be equivalent to maximise the log likelihood:

6.17.1 Maximum Likelihood of Multinomial Cell Probabilities

X₁,X₂,...,X_m are counts in cells/ boxes 1 up to m, each box has a dierent probability (think of the boxes being bigger or smaller) and we x the number of balls that fall to be n:X₁ + X₂ + + X_m = n. The probability of each box is p_i, with also a constraint: p₁ + p₂ + + p_m = 1, this is a case in which the X_i's are NOT independent, the joint probability of a vector x₁,x₂,...,x_m is called the multinomial ,and has the form:

Each box taken separately against all the other boxes is a binomial, this is an extension thereof.

We study the log-likelihood of this :

We use

6.17.2 The Hardy Weinberg Example

This is a trinomial with three boxes: the probabilities are parametrized by:

Each of the counts is binomially distributed with probabilities as described above so that:

We can either use asymptotic maximum likelihood or the bootstrap, we are going to compare the two using Splus.

After some algebra we nd:

Of course as usual we don't know we have to plug in the estimate:

The bootstrap procedure we are going to follow generates multinomials with the estimated probabilites and then computes for these new counts the bootstrap estimate for which we then compute the standard deviation.

Here is the R code I used:

Chapter 7
Density Estimation

Chapter 8
References

CRAN
Numerical Recipes http://www.library.cornell.edu/nr/cbookcpdf.html
Anuj Srivastava's Notes: http://calais.stat.fsu.edu/anuj/PDF-les/STA5107-s-05/BOOK.pdf
Julian Faraways Anova and Regression with R: http://www.stat.lsa.umich.edu/~faraway/book/
http://www.quantlet.com/mdstat/scripts/csa/html/node45.html
http://www.neurosci.aist.go.jp/ akaho/MixtureEM.html

http://www.riskglossary.com/articles/monte_carlo_method.htm

MCMC notes:
The Gibbs Sampler Motif Page
Gibbs sampler for a Normal
Peter Green's Reversible Jump

Bootstrap Course:
Susan Holmes Stat 208 course

Density Estimation:
http://www.stat.berkeley.edu/~stark/Teach/S240-s03/Notes/ch8.pdf

Smoothing:
John Fox introduction to nonparametric regression

MAchine Learning Data Repository
http://www.ics.uci.edu/~mlearn/MLRepository.html