###################################################
# Computer Lab 2 – F71SM
###################################################
## In the tasks below, some (but not all) R code is given
## for guidance. In some places the code is also incomplete,
## so you will need to figure out how to fill in the details.
## Make sure you use the help menus in R (e.g.”?mean” will
## open up a help window for “mean()”).
# Tutorial 1, Question 2:
# Repeat the calculations
# The data can be found in file T1Q2.RData
load(“T1Q2.RData”)
sum(x); sum(x^2)
q = quantile(x,probs=c(0.25,0.75)); q
iqr = q[2] – q[1]; iqr
# Note that R by default uses a different definition of quantiles
# The definition in lectures notes can be obtained as
q.notes = quantile(x,probs=c(0.25,0.75), type=6); q.notes
iqr.notes = q.notes[2] – q.notes[1]; iqr.notes
###################################################
# Tutorial 1, Question 7:
# Confirm the answers by computing the mean, sd and skewness
# of the two samples
a = rep(c(1:4),c(80,60,40,20))
b = rep(c(1:4),c(110,50,30,10))
mean(a); mean(b)
sd(a); sd(b)
# For skewness, use function from lab 1
beta1.fn(a)
beta1.fn(b)
###################################################
# Worked example 2.4:
# Confirm the answers for all values of k
p = numeric(5)
# Note that p[0] is not defined for a vector in R
for(k in 0:4){
p[k+1] = choose(4,k)*choose(7,(6-k))/choose(11,6)
###################################################
# Simulate 10,000 values throws of a fair six-sided die
# Let X be the r.v. giving the score of the die
# Compute (approximately) the mean and variance of X
# (this is also worked example 4.1).
# A single value of X can be simulated as (check):
trunc(6*runif(1, 0, 1)) + 1
# So 10^4 values simulated as
x = trunc(6*runif(10000, 0, 1)) + 1
mean(x); var(x)
# Compute (approximately) P(X<=3)
length(x[x<=3]) / length(x)
# You can repeat these calculations to see how the anwsers
# change slightly
###################################################
# Worked example 3.3:
# Simulate 10,000 values of r.v. Y, where Y is
# number of throws until we get a "6"
# Confirm the mean and variance of Y
# A sigle value of y can be simulated as follows:
# Generate throws of the die as in Task 3
# and stop when a 6 is thrown
while(x!=6){
x = trunc(6*runif(1, 0, 1)) + 1
# Use for loop to simulate 10^4 values of Y
y = numeric(10^4)
for (i in 1:10^4){
while(x!=6){
x = trunc(6*runif(1, 0, 1)) + 1
y[i] = y[i]+1
# Again, repeat these calculations to see how the anwsers
# change slightly
###################################################
# Worked examples 3.5, 3.6:
# Plot the pdf and cdf of r.v. X
x1 = seq(1.001,4,by=0.01)
f.x1 = 3/x1^4
F.x1 = 1-x1^(-3)
# Change the plotting window layout
par(mfrow=c(1,2))
plot(x1,f.x1,type="l",xlab="x", ylab="f(x)", main="PDF")
plot(x1,F.x1,type="l",xlab="x", ylab="F(x)", main="CDF")
# Can also change the limits of the x axis:
plot(x1,f.x1,type="l",xlim=c(0,3.3), xlab="x", ylab="f(x)", main="PDF")
plot(x1,F.x1,type="l",xlim=c(-1,3.3),xlab="x", ylab="F(x)", main="CDF")
###################################################
# Worked example 4.2
# X ~ Binomial(10,0.2)
# Calculate pmf f(3), cdf F(3), i.e. P(X=3), P(X<=3)
dbinom(3,10,0.2)
pbinom(3,10,0.2)
# Simulate 10^3 values of X and confirm E(X), Var(X)
# E(X)=... ; Var(X)=...
x = rbinom(10^3, size=10,prob=0.2)
mean(x); var(x)
###################################################
# Worked example 4.3
# X ~ Poisson(7)
# Calculate pmf f(9), P(X>=8) = …
dpois(9,7)
1-ppois(7,7)
# Simulate N=10^4 values of X and confirm E(X), Var(X)
# E(X) = V(X) = …
# Use rpois() …
x = rpois(10^4,7)
mean(x); var(x)
# Repeat above with N=10^2. Comment.
###################################################
# Worked example 4.4
# X ~ Binomial(200,0.02)
# Calcualate P(X = 5 or X = 6) using the binomial above and a Poisson approximation
# P(X=5 or X=6) = …
pbinom(6,200,0.02) – pbinom(4,200,0.02)
ppois(6,200*0.02) – ppois(4,200*0.02)
###################################################
## Task 10
# Worked example 4.5
# X ~ Geo(0.3)
# Calcualate P(X = 1, 2 or 3)
pgeom(3,0.3) – pgeom(0,0.3)
###################################################
## Task 11
# Worked example 4.6
# X ~ Exp(0.2)
# Calcualate P(6 < X < 10), P(6 <= X < 10)
pexp(10,0.2) - pexp(6,0.2)
# Simulate N=10^4 values of X and confirm E(X), Var(X), P(6 < X < 10)
# Use rexp() ...
x = rexp(10^4, rate=0.2)
mean(x); var(x)
( length(x[x<10]) - length(x[x<6]) ) / length(x)
# Also check to compute this using the CDF:
exp(-1.2) - exp(-2)
# Plot the pdf and cdf of r.v. X
x1 = seq(0.001,30,by=0.01)
f.x1 = 0.2*exp(-0.2*x1)
# or f.x1 = dexp(x1,0.2)
F.x1 = 1-exp(-0.2*x1)
# or F.x1 = pexp(x1,0.2)
# Change the plotting window layout
par(mfrow=c(1,2))
plot(x1,f.x1,type="l",xlab="x", ylab="f(x)", main="PDF")
plot(x1,F.x1,type="l",xlab="x", ylab="F(x)", main="CDF")
###################################################
## Task 12
# Worked example 4.7
# X ~ N(500,10)
# Calculate P(493 < X < 505), 99th quantile
pnorm(505,500,sqrt(10)) - pnorm(493,500,sqrt(10))
qnorm(0.99,500,sqrt(10))
###################################################
## Task 13
# X ~ Gamma(4,3)
# Calculate 0.95 quantile
qgamma(0.95, shape=4, rate=3)
# Y = 6*X ~ Chi-square(8)
# Calculate 0.95 quantile of X and confirm that it is the same as above
qchisq(0.95,8)/6