1 Instructions

1.  Reports should be handed in as a single PDF file using Blackboard, by noon on the due date.

RMarkdown, image, word, zip files, for example, will not be marked.

2. You can work alone or in a group with up to two other people.

3.  One person per group should hand in a report. The names and student numbers of all group members must be on the first page.

4. Your answers should combine text and code snippets in R. It is recommended that you use RMarkdown to prepare your reports, since this is typically easier for students, but this is not mandatory.

5. You must explain what you are doing clearly to obtain full marks on each question.  You can use comments (which start with #) to annotate your code. Mathematical derivations can be written using LaTeX commands in RMarkdown or on paper, with a photo then appended to the end of the PDF being submitted.

6. This practical counts for 10% of your assessment for Statistics 2.

2 Populations of world cities

The data is population sizes of the world’s cities and towns, according to the World Cities Database. To load the raw data and remove missing (NA) entries, we run the following commands.

raw.data  <- read.csv ( "worldcities.csv")

raw.populations  <- raw.data$population

populations  <- raw.populations[ ! is.na(raw.populations)]

sum (populations)

##  [1]  5076372338

There are around 5.076 billion people counted, which is off by about 3 billion since there are about 8 billion people in the world. In any case we can try to understand the distribution of population sizes for the data we do have. Now we remove populations under 10000.

populations  <- populations[populations>=10000]

length(populations)

##  [1]  38615

sum (populations)

##  [1]  5030436026

We can see that removing the smaller towns still leaves us considering around 5.03 billion people, and we have observations of population sizes from 38615 cities.

To finish cleaning the data, we now divide populations by 10000 and subtract 1 so that the values start from 0 and we essentially counting population sizes in units of 10000.

populations  <- populations/10000-1

3    Visualizing the data

The data have a large range, and so visualizing a density estimate on the original scale is not very informative. Instead, and purely for the purpose of visualization, we can plot the density of the logarithms of the  observations.

plot (density (log(populations)))

density(x = log(populations))

−10                       −5                         0                          5

N = 38615   Bandwidth = 0.1845

4 Model

We will model populations as realizations of independent and identically distributed random variables from a member of the family of distributions with density

f (x; θ) = I(x ≥ 0),                                                            (1)

with two-dimensional statistical parameter θ = (a, b) ∈ (0, ∞ )2. This is a regular statistical model. This PDF is non-negative and integrates to 1 because for any b > 0 and t ∈ R, we have

l0 ∞ dx = {

1

(t−1)bt − 1

∞

t > 1

t ≤ 1. (2)

We shall write X ∼ Dist(a, b) to mean that X is a random variable with density (1). To simulate a vector of realizations of independent Dist(a, b) random variables, we can use the following code. If you are interested, this is an application of inverse transform sampling.

rdist  <- function (n,  a,  b)  {

u  <- runif (n)

(uˆ(-1/a)-1)*b

}

To evaluate the PDF at each point in a vector, we can use the following code. If the log argument is true then the log density is returned rather than the density, as is usual in R code.

ddist  <- function (xs,a,b,log=FALSE)  {

if (log)  {

log(a) + a*log(b) - (a+1)*log(xs+b)

} else {

a*bˆa/ (xs+b)ˆ(a+1)

}

}

5 Estimation

Using (2) it is quite easy to show that the first and second moments of X ∼ Dist( · ; θ) are

E(X;θ) =

and

E(X2 ; θ) = − b(−2)

a > 1,

a ≤ 1.

a > 2,

a ≤ 2.

Question 1.  [2 marks] Derive the method of moments estimators of a and b.  Report the estimates for this dataset.

I suggest that you write a function mom.estimate that takes as input some data and returns a vector containing the estimated values of a and b. For example, something like this:

mom.estimate  <- function (xs)  {

## your  code  here

a.hat  <- 0 # obviously wrong

b.hat  <- 0 #  obviously  wrong

c (a.hat,  b.hat)

}

You can check if your estimate makes sense by considering a large dataset with known parameters.

fake.data  <- rdist (100000 ,  10 ,  20)

fake.theta.hat  <- mom.estimate (fake.data)

fake.theta.hat

##  [1]  0  0

Having computed the method of moments estimate for the populations data, we can visualize the fitted density on the log scale. Even if you have computed the estimates properly, the fit should not look great. As you go through the rest of the practical, you may be able to guess why this is.

On the logarithmic scale, we can plot the fitted density by using the transformed PDF suggested by Exercise 2.1 of the lecture notes, i.e. if Y = log(X) and fX  is the PDF for X then the PDF for Y is

fY (y) = fX (exp(y)) exp(y).

theta.mom  <- mom.estimate (populations)

plot (density (log(populations)),ylim=c (0 ,0.3))

zs  <- seq (-10 ,10 ,0.01)

lines (zs,ddist (exp (zs),theta.mom[1],theta.mom[2])*exp(zs),col="red")

density(x = log(populations))

−10

N = 38615   Bandwidth = 0.1845

Question 2.  [2 marks] Show that for fixed b, the log-likelihood is maximized by taking

a =

The log-likelihood cannot be maximized analytically, so we will use a numerical optimization procedure. In particular, we will make use of the fact that we can maximize a function of b only since we know that for any b the maximizing a is given by (3).

ml.estimate  <- function (xs)  {

n  <- length(xs)

get.a  <- function (b)  {

n/ (sum (log(xs+b))-n*log(b))

}

ell  <- function (b)  {

lb  <- log(b)

a  <- get.a (b)

la  <- log(a)

n*la+a*n*lb- (a+1)*sum (log(xs+b))

}

out  <- optim (1 ,  ell,  control  = list (fnscale=-1), method= "Brent" ,  lower=0 ,  upper=100) c (get.a (out$par),out$par)

}

Now we can visualize the fitted density using the ML estimate. It looks pretty good!

theta.ml  <- ml.estimate (populations)

plot (density (log(populations)))

zs  <- seq (-10 ,10 ,0.01)

lines (zs,ddist (exp (zs),theta.ml[1],theta.ml[2])*exp(zs),col="red")

density(x = log(populations))

−10

N = 38615   Bandwidth = 0.1845

6 Confidence Intervals

In lectures we have seen that for regular statistical models with a one-dimensional parameter θ, the ML estimator θ(ˆ)n  is asymptotically normal with

*nI(θ)(θ(ˆ)n − θ) →D(·;θ)  Z ∼ N (0, 1).

This convergence in distribution justifies the construction of Wald confidence intervals for θ .

assumptions, the ML estimator θ(ˆ)n  is asymptotically (multivariate) normal in the sense that

√n(θ(ˆ)n − θ) →D(·;θ)  W ∼ N (0, I(θ) — 1 ),

where I(θ) is the Fisher information matrix

I(θ) = −E[∇2ℓ(θ; X1 ); θ],

and the expectation is taken element-wise. That is, the ijth entry of I(θ) is the negative expectation of the corresponding ijth second order partial derivative of the log-likelihood associated with 1 observation. Notice that W is a two-dimensional normal random vector with mean 0 and variance-covariance matrix I(θ) — 1 .

One can deduce from this multivariate asymptotic normality that for j ∈ {1, . . . , 2},

√n(θ(ˆ)n,j − θj ) →D(·;θ)  W1 ∼ N (0, (I(θ) — 1 )jj ),

which we can rewrite as

(θ(ˆ)n,j − θj ) →D(·;θ)  Z ∼ N (0, 1),

where θ(ˆ)n,j  denotes the jth component ofθ(ˆ)n  and θj  denotes the jth component of θ .

Notice that (I(θ) — 1 )jj  is the jth diagonal entry of the inverse of the Fisher information matrix, and is not in general equal to (I(θ)jj ) — 1 , the inverse of the jth diagonal entry of the Fisher information matrix1 .  In R you can compute numerically the inverse of a matrix using the solve command.

As in the one-dimensional parameter case, one can replace I(θ) with I(θ(ˆ)n ) if θ(ˆ)n  is a consistent sequence of

estimators of θ to obtain

( (I(θ(ˆ)n— 1 )jj (θ(ˆ)n,j − θj ) →D(·;θ)  Z ∼ N (0, 1),

Question 3.  [2 marks] Show that the Fisher information for one observation is

I(θ) = ( − b 2b()   ) ,

and compute the Fisher information for θ =θ(ˆ)ML .

I recommend that you write a function to compute the Fisher information for any value of θ, structured like this:

info  <- function (theta)  {

a  <- theta[1]

b  <- theta[2]

# your code

i11  <- 1 # obviously wrong

i12  <- 0 # obviously wrong

i21  <- 0 # obviously wrong

i22  <- 1 # obviously wrong

mtx  <- matrix (c (i11,i21,i12,i22),2 ,2)

return (mtx)

}

Question 4.  [2 marks] Compute and report the endpoints of the observed 98% Wald confidence intervals for a and b individually.

I recommend that you write a function that computes the observed confidence intervals’ end points, which should call the info function.

compute.CI.endpoints  <- function (xs,alpha)  {

# your code

L.a  <- 0 # obviously wrong

U.a  <- 1 # obviously wrong

L.b  <- 0 # obviously wrong

U.b  <- 1 # obviously wrong

list (a=c (L.a,U.a),b=c(L.b,U.b))

}

You may find it helpful to know that you can use the solve function to compute the inverse of a matrix.