Introduction to R - Exercise Solutions

Emmanuelle Comets

1 Course 1
1.1 First steps
Simple calculations, use of help:

10*9*8*7*6*5*4*3*2*1 log?
log(2) log10(2)
log(2)/log(10) log(2, base=10)

1.2 Vectors
Creation of a vector (1)

vec<-c(10,3,4,5,6,10,100,100,10,20, vec<-c(10,3:6,10,rep(100,2),
30,40) seq(10,40,by=10)

Creation of a vector (2)

x<-c(1,4,5) y[-2]
x xy <- y[c(1,4,5)]
y <- seq(1, 9, 2) # or alternatively : or since x=c(1,4,5)
y = -2 * (1 to 5) - 1 xy <- y[x]
y xy
y[2] y[2:4]

Operation on vectors

y+1 (y+1)[1:4]
y[1:4]+1 x<-2:4
or an equivalent writing y multiplied by x

Vector manipulation

vec<-rnorm(10) yvec<-log(vec)
vec vec2 <- yvec[]
length(vec[vec>0]) length(vec2)
1
1.3 Tables
Creation of a matrix

mat<-matrix(1:15,ncol=5,byrow=T) mat[mat[:,1]<3,1]
mat mat[mat[,1]<3,]
mat[2:3, c(2,4)]

2 Course 2
Warm-up (1):

rep(c(0,6),3)
c(1:4)*3-2 1+(1:10)/10
rep(1:3, 4) c(0:9)*2+1
rep(1:3,1:3) rep(-2:2,2)
rep(1:3,3:1) rep(-2:2, each=2)
1 + 0 : 2 * 4.5 (1:10)*10
rep(1:3,each=4)

Warm-up (2):

d<-rep(-2:2,each=2) sum([Link](d))
d[d<0]<-NA d<-rep(-2:2, each=2)
length(d[[Link](d)]) d[d<0] <-- 10
or

Warm-up (3):

x<-matrix(1:120,ncol=12) x[apply(x,1,mean)<60,apply(x,2,sum)<500]
x[2*(1:5),] x[apply(x, 1, mean) < 60 & x[, 1] != 3,
x[apply(x,1,mean)<60,] apply(x,2,sum)<500
x[,apply(x,2,sum)<500]

In the last writing, we use the fact thatx[apply(x,1,mean)<60,apply(x,2,sum)<500]

itself a table, so we can reference its elements... A less concise and more
Claire would be to make 2 lines using an intermediate table:

x1 <- x[apply(x, 1, mean) < 60, apply(x, 2, sum) < 500]

x1[-3,]

2
2.1 Dataframe
Data extraction

air quality
air1 <- airquality[airquality$Temp > 92, ]
air1 <- transform(air1, logTemp = log(Temp))
air1 <- transform(air1, ftemp = ifelse(Temp > 94, 1, 0))
air2<-air1[ & air1$Temp<=94,]
air2

air3 <- airquality[,]

air3[1:10,] # to see the first 10 lines of air3
head(air3) # to see the first 6 lines of air3

air3 <- transform(air3, monindic = ifelse(Month < 4 & Temp > 80, 1, 0))

Exercise on match

dates<-c("1971-01-20","1971-01-28","1971-02-03","1971-02-11","1971-02-18",
"1973-01-17","1973-01-25","1973-01-31","1973-02-17","1974-01-07","1974-01-10",
"1974-01-15","1974-01-22","1974-01-29","1974-02-05","1974-02-12","1974-02-19")
mesure<-c(64,69,71,71,71,32,42,28,32,18,25,29,34,36,42,50,61)

dates[match(unique(measure),measure)]
mat<-cbind(measure=measure,dates=dates)
mat<-mat[order(mat[,1]),]

2.2 Descriptive statistics

Quantiles :

x<-matrix(1:100,ncol=4)
apply(x, 2, quantile, c(0.1, 0.9))
apply(x, 1, quantile, c(0.1, 0.9))

Average, variance:

x<-matrix(1:100,ncol=4) apply(x, 2, var)

mean(x) apply(x[1:3,], 1, mean)
apply(x, 2, mean) apply(x[1:3,],1,var)

3
Correlation:

ToothGrowth xmat<-matrix(ToothGrowth[,1],ncol=6)
head(ToothGrowth) cor(xmat)
str(ToothGrowth)

2.3 Probability laws and simulation

Sample simulation:

Drawing 6 values from N(2,5), 4 from N(-1,4)

vec <- c(rnorm(6, 2, 5), rnorm(4, -1, 4))
Drawing from the mixture using a variable in N(0,1) (symmetric distribution)
v1 <- c()
for(i in 1:10) {
i1<-rnorm(1)
v1 <- c(v1, ifelse(i1 > 0, rnorm(1, 2, 5), rnorm(1, -1, 4)))
}
m1 <- mean(v1)
Note: we can do it without the loop
i1<-rnorm(10)
v1 <- ifelse(i1 > 0, rnorm(10, 2, 5), rnorm(10, -1, 4))

Calculation of the average m1 of v1, and repetition (a new loop) 10 times

m1<-c()
for(j in 1:10) {
v1 <- c()
for(i in 1:10) {
i1<-rnorm(1)
v1 <- c(v1, ifelse(i1 > 0, rnorm(1, 2, 5), rnorm(1, -1, 4)))
}
m1<-c(m1,mean(v1))
}
Drawing in a mixture of probabilities 10-90%
m2 <- c()
for(j in 1:10) {
v2<-c()
for(i in 1:10) {
i1<-rnorm(1)

4
 v2 <- c(v2, ifelse(i1 > qnorm(0.9), rnorm(1, 2, 5), rnorm(1, -1, 4)))
}
m2<-c(m2,mean(v2))
}
Averages of m1 and m2
mean(m1)
mean(m2)

Probability that the average is less than 130, n=10 subjects:

Let's note X the variable 'creatinine clearance'.
X¯- µ 130 - µ 130−µ
¯
P(X≤130) =P( =P(Z≤ (1)
σ/n√ ≤ s/n√ √ )
sigma/n
where Z is the standardized variable associated with X, which follows a lawN(0,1).
So we are looking for the probability:

¯ 130−120
P(X≤130) = P(Z≤ P(Z≤0.79)
√
40/10
moy<-120 #By relating it to N(0,1)
ect<-40 z <- (a1 - moy) / (ect / sqrt(nsuj))
a1<-130 pnorm(z)
nsuj<-10 #Using pnorm
pnorm(a1, mean, std/sqrt(sample size))

Probability that the average is between 120 and 130:

¯ 130 - µ
120 - µ X - µ 130−µ 120 - µ
¯
P(120≤X≤130) =P( =P(Z≤ )−P(Z≤ ) (2)
σ/n√ ≤ σ/n√ ≤ σ/n√ σ/n√ √
sigma/n
And we deduce:

b1<-120 pnorm(z1)-pnorm(z2)
nsuj<-10 Directly
z1 <- (a1 - mean) / (ect / sqrt(nsuj)) pnorm(a1, mean, ect/sqrt(nsuj))-
z2 <- (b1 - moy) / (ect / sqrt(nsuj)) pnorm(b1, mean, standard_deviation/sqrt(number_of_subjects))

Number of subjects needed for the probability (1) to be at least 95%:

nsn <- (qnorm(0.95) * ect / (a1 - moy)) ** 2

nsn<-ceiling(nsn)

z <- (a1 - moy) / (ect / sqrt(nsn))

pnorm(z)

5
3 Course 3
Warm-up (1):

vec1<-rnorm(20,70,sqrt(10)) age=vec2,douleur=vec3)
vec2 <- rnorm(20, 25, sqrt(4)) mean(essai$weight)
vec3 <- trunc(runif(20, 0, 5)) var(trial$weight)
test <- [Link](weight = vec1,

3.1 Statistical tests

Tests of mean and variance

With the above base:

[Link](essay[,1]~[Link](essay[,3]>=2))
[Link](essay[,1]~[Link](essay[,3]>=2))
#ou
[Link](experiment[experiment[,3]>=2,1],experiment[experiment[,3]<2,1])
[Link](essai[,1]~[Link](essai[,3]>=2),exact=T)

With a Wilcoxon test

A<-c(0,1,2) A<-c(0,1,2,3,4)
B<-c(100,150,5000) [Link](A,B)
[Link](A,B) B<-c(100,150,5000,6000,500)
B<-c(100,150,5000,6000) A<-c(0,1,2,3,4,5)
[Link](A, B) [Link](A,B)

Despite the difference in averages, at least 4 elements from A and B are needed to detect.
a difference.
The t-test does even worse...

A <- c(0,1,2) [Link](A,B)

B <- c(100, 150, 5000) B<-c(100,150,5000,6000,500)
[Link](A,B) A<-c(0,1,2,3,4,5)
B<-c(100,150,5000,6000) [Link](A,B)
A <- c(0, 1, 2, 3, 4)

6
ANOVA

airquality[1:10,]
[Link](airquality$Ozone[airquality$Month==5],
airquality$Ozone[airquality$Month==8]
#ou
ozconc <- airquality$Ozone
ozmonth<-airquality$Month
[Link](ozconc[ozmonth==5],ozconc[ozmonth==8])
[Link](ozconc[ozmonth==5],ozconc[ozmonth==8])

anova(lm(Ozone~[Link](Month),data=airquality))
[Link](Ozone~[Link](Month),data=airquality)

Distribution tests, χ2
Exercise onair quality:

ozconc <- airquality$Ozone

oztemp <- airquality$Temp
v1 <- ozconc[ & ]
oztemp <- oztemp[ & ]
ozconc<-v1
x1<-matrix(c(length(ozconc[ozconc>75 & oztemp>85]),
length(ozconc[ozconc>75 & oztemp<=85]),length(ozconc[ozconc<=75 & oztemp>85]),
length(ozconc[ozconc<=75 & oztemp<=85])),ncol=2,byrow=T)

[Link](x1)
[Link](airquality$Ozone,"pnorm")

3.2 Graphs
Graph of ozone concentrations for the first 20 days of May

Ozone observed in the first 20 days of May

ozc <- airquality$Ozone
ozm <- airquality$Month
ozd<-airquality$Day

plot(ozd[ozm==5],ozc[ozm==5],type="b",xlab="days",
xlim=c(0,20),ylim=c(0,50),pch=3,ylab="ozone",main=tit)

7
Graph of the relationship between ozone and wind, according to temperature

tit<-"Relation entre l’ozone et le vent, selon la température"

plot(ozw[ & ozt<85],ozc[ & ozt<85],xlab="Vent (mph)",
ylab="Ozone (ppb)",pch=2,xlim=c(2,15),ylim=c(35,125),main=tit)
points(ozw[ & ozt>=85],ozc[ & ozt>=85],pch=6)
legend(12,120,c("Temp>=85˚F","Temp<85˚F"),pch=c(2,6))

Histograms of the swiss database

idens<-c(-1,-1,25,-1)
icol<-c("orange","white","red","peachpuff")
["gray0","red","red","red"]
par(mfrow=c(2,2))
for(i in 2:5) {
hist(swiss[,i],xlab=names(swiss)[i],main="",breaks=20,border=ibord[i-1],
col=icol[i-1],density=idens[i-1]
}

Swiss base mustache boxes

binedu <- [Link](swiss$Education > 10)

binmor<-[Link](swiss$[Link]>20)
par(mfrow=c(1,2))
boxplot(swiss$Fertility~binedu,xlab="Education",ylab="Fertility")
boxwex=0.2,col="orange"
legend(0.5,45,c("0: Edu<10","1: Edu>10"))
boxplot(swiss$Fertility~binmor,xlab="Education",ylab="Fertility",
boxwex=0.2, col="red"
legend(1.2,45,c("0: [Link]<20","1: [Link]>20"))

Exercise
Density of 4 laws:

par(mfrow=c(2,2))
#N(0,1)
xpl<-c(-50:50)/10
ypl<-dnorm(xpl)
plot(xpl,ypl,xlab="X",ylab="Density of normal law",type="l")

log-normal law

8
xpl<-c(0:50)/10
ypl<-dlnorm(xpl)
plot(xpl, ypl, xlab='X', ylab='Density of the log-normal law', type='l')

Poisson's law
lambda <- 2
xpl<-c(0:40)/2
ypl<-dpois(xpl,lambda)
plot(xpl, ypl, xlab="X", ylab=paste("Poisson distribution density (lambda=", lambda, ")"))
sep=""),type="l")

#Gamma law
gam<-2
bet<-1
xpl<-c(0:50)/5
ypl <- dgamma(xpl, gam, bet)
plot(xpl, ypl, xlab="X", ylab=paste("Density of Gamma distribution (gamma=", gam, ", beta=")
bet, ")", sep="" ), type="l"

4 Course 4
Leg warm-up (1):

women?
data<-women
dat[,1]<-dat[,1]*2.5
dat[,2]<-dat[,2]/2
plot(dat[,1],dat[,2],xlab="Size (cm)",ylab="Weight (kg)",type="b"
main="Taille et poids moyen chez des femmes américaines de 30 à 39 ans")

Warm-up (2):

?CO2
data<-CO2
par(mfrow=c(1,2))
hist(dat[dat[,3]=="chilled",5],xlab="CO2 uptake",main="Chilled")
hist(dat[dat[,3]=="nonchilled",5],xlab="CO2 uptake",main="Non-chilled")
[Link](dat[,5]~dat[,3])

9
4.1 Statistical analyses
Linear regression:

library(MASS)
pairs(cats)
[Link] <- lm(Hwt~Bwt*Sex,data=cats)
summary([Link])

par(mfrow=c(1,1))
qqnorm(studres([Link]))
qqline(studres([Link]))

plot(predict([Link]), studres([Link]), xlab="Predictions",

ylabel="Standardized residuals"
abline(h=0)
plot([Link], which=4)

attributes(summary([Link]))
summary([Link])$[Link]
summary([Link])$df
summary([Link])$[Link]

cats.lm1 <- update([Link], Hwt ~ Bwt + Sex)

summary(cats.lm1)
anova(cats.lm1, [Link])

Logistic regression:

library(MASS)
?birthwt
lbwt <- [Link](low=factor(birthwt$low),age=birthwt$age,ptl=birthwt$ptl,
smoke=(birthwt$smoke>0),ht=(birthwt$ht>0))
[Link] <- glm(low~age+ptl+smoke+ht,data=birthwt,family=binomial)
summary([Link])
drop1(glm(low~age+ptl+smoke+ht,data=birthwt,family=binomial),test="Chisq")

lbwt.glm1 <- update([Link], . ~ . - smoke)

anova(lbwt.glm1, [Link], test="Chisq")
drop1(lbwt.glm1, test="Chisq")

10
lbwt.glm2 <- update(lbwt.glm1, . ~ . - age)
anova(lbwt.glm2, lbwt.glm1, test="Chisq")
drop1(lbwt.glm2, test="Chisq")
lbwt.glm3<-update(lbwt.glm2,.~ptl*ht)
summary(lbwt.glm3)
we stay with lbwt.glm2
summary(lbwt.glm2)

plot(lbwt.glm2)

p1<-predict(lbwt.glm2,type="response")
p1 <- [Link](p1 >= 0.5)
table(p1, birthwt$low)

With step
[Link]<-paste(names(birthwt)[2:10],collapse="+")
[Link] <- paste(names(birthwt)[1], [Link], sep="~")
[Link] <- glm([Link]([Link]), family=binomial, data=birthwt)
step([Link])

Calculation of the necessary number of subjects:

[Link](power=0.95, p1=0.4, p2=0.6)

xpmin<-0.42
delmin <- xpmin - 0.4
for(xp2 in seq(0.6,xpmin,by=-0.02)) {
y<-[Link](n=1000,p1=0.4,p2=xp2)
cat("delta=", xp2 - 0.4, " power=", y$power, "\n")
if(y$power>=0.95) delmin<-xp2-0.4
}
y <- [Link](n = 1000, p1 = 0.4, p2 = 0.4 + delmin)
cat("With a power of ",y$power," we detect an effect of ",delmin,"\n")

4.2 Loops
Arrange the students:

Take a student at random and stand them up

11
For each student still seated:
to get up
compare themselves to the student standing farthest to the left

If the student is taller, move forward one step.

repeat until the student has reached the end of the line
Stop the algorithm when all students are standing.
Tests :
{
x<-scan("",nlines=1)
if(x>0) cat("This number is positive\n") else cat("This number is negative\n")
}
Braces allow you to execute a block of commands at once. Compare with and without
the braces.
The following script must be executed in two stages, first the first 2 lines then the following ones:
itir<-sample(0:100,1)
x <- scan("", nlines=1)
if(x==itir) cat("It's won!\n") else {
if(x>itir) cat("Too high\n") else cat("Too low\n")
cat("The number was:", itir, "\n")
}
A nicer way to do this is to define a function (see lesson 5):
guess<-function() {
itir<-sample(0:100,1)
x<-scan("",nlines=1)
if(x==itir) cat("It's won!\n") else {
if(x>itir) cat("Too high\n") else cat("Too low\n")
cat("The number was :",itir,"\n")
}
}
then execute it:
guess()
1:45
Read 1 item
Too high
The number was: 35

12
Additional exercise: modify the function to play (requires loops in the rest of the course
4) to guess the correct number by successive attempts.

Boucles

x<-matrix(1:120,ncol=12)
Do not forget to define and initialize lmoy before the loop
lmoy<-c()
for(i in 1:dim(x)[1]) lmoy<-c(lmoy,mean(x[i,]))
another possibility
lmoy<-NULL
# or lmoy<-0
for(i in 1:dim(x)[1]) lmoy[i]<-mean(x[i,])

cmoy<-c()
for(i in 1:dim(x)[2]) cmoy<-c(cmoy,mean(x[,i]))
Note: we can also use apply
apply(x, 1, mean)
apply(x, 2, mean)

Exercise:

{ x<-matrix(1:120,ncol=12)
for(i in 1:3) cat(LETTERS[i]) row sums of x
apply(x, 1, sum)
} colSums(x)
{ apply(x, 2, sum)
for(i in c(5,24,5,18,3,9,3,5)) rowMeans(x)
cat(letters[i]) apply(x, 1, mean)
colMeans(x)
} apply(x, 2, mean)

4.3 TCL
Illustration exercise of TLC, draw with nobs=3 observations:

x<-rexp(500,5)
hist(x)
nobs<-3

13
x <- rexp(nobs, 5)
mean(x)
my<-c()
for(i in 1:500) mean_values <- c(mean_values, mean(rexp(nobs, 5))

hist(moy)

Loop to vary nobs:

nobs<-c(3,5,10,30) vec<-rexp(i,5)
par(mfrow=c(2,2)) moy <- c(moy, mean(vec))
for (i in nobs) { }
moy<-c() hist(moy, breaks=20)
for (j in 1:50) { }

This exercise is an illustration of the central limit theorem: when nobs is small, the distri-
The distribution of the average of n observations is close to the distribution of the sampled variable.
On the other hand, when nobs increases, the distribution of the means approaches that of a normal distribution.
male. It is noted that for the exponential law, one must wait at least nobs=30 to start
to have something satisfying.

Same exercise for a uniform law:

par(mfrow=c(2,2)) moy<-c(moy,mean(vec))
for (i in nobs) { }
moy<-c() hist(moy, breaks=20)
for (j in 1:100) { }
vec<-runif(i)

Here, the histogram begins to resemble that of a normal distribution almost immediately.

5 Course 5
Warm-up:

library("ISwR")
plot(bp~obese,pch = ifelse(sex==1, 1,2), data = [Link],xlab="Obesity ratio"
ylab="Blood pressure (mm Hg)"
legend(2,200,c("Women","Men"),pch=1:2)

y<-lm(bp~obese,data = [Link])

14
summary(y)
y2 <- lm(bp ~ obese + sex, data = [Link])
summary(y2)
anova(y2, y, test="Chisq")
co<-coef(y2)
plot(bp~obese,pch = ifelse(sex==1, 1,2), data = [Link],xlab="Obesity ratio"
ylab="Blood pressure (mm Hg)"
abline(y)
abline(a=co[1]+co[3],b=co[2],col="red",lty=2)
abline(a=co[1],b=co[2],col="blue",lty=2)
legend(1.6,200,c("Both (model 1)","Women (model 2)","Men (model 2)"),lty=c(1,2,2),
col=c("black","red","blue"))

5.1 Strings
grep("er$",[Link]) # returns the months ending in er (in English)
grep("^[A-J]", [Link]) # returns the months starting with a letter
between A and J (in English)
grep('^[^J]', [Link]) # returns the months not starting with J

chain<-"chain"
gsub("i","e",chaine)
gsub("[a-e]","x",chaine)

Exercises ongrep:

This is a string
gsub("i","i",chaine)
gsub("[a-e]", "x", string)

String manipulation:

["Anne Dubois","Julie Bertrand","Emmanuelle Comets","Caroline Bazzoli"]

Hervé Le Nagard

Here I define 2 functions to extract the 2 parts first name/last name (see second part of the course)
5). I suppose not to complicate matters that if the first name is composed, it is in the form
Jean-Pierre (and therefore the first box of the extracted vector contains the full name).

extract.first_name<-function(vec) {
return(vec[1])}

15
[Link] <- function(vec) {
return(paste(vec[2:length(vec)],collapse=" "))}

[Link]<-strsplit(student," ")
student.first_name <- unlist(lapply([Link], extract.first_name))
[Link] <- unlist(lapply([Link], [Link]))
print(student.first_name)
print([Link])

[Link]<-paste(student.first_name,student.last_name,sep="-")
print([Link])

student.first_name[grep("i",student.first_name) | ]grep("I",student.first_name)

student.first_name[grep("^[A-M]", student.first_name)]

Lecture/writing:

[Link]<-c("20-2-1978","12-10-1978","25-05-1971","8-8-1979","1-6-1971")
anniv<-[Link](prenom=[Link],nom=[Link],naiss=[Link])
[Link](anniv,"[Link]",[Link]=F)

anniv2 <- [Link]("[Link]", header = T)

[Link](anniv,anniv2)

vec <- strsplit([Link], "-")

[Link]<-unlist(lapply(vec,paste,collapse="/"))
anniv2[,3] <- [Link]
anniv[order(anniv[,2]),]

demog <- [Link](first_name = student.first_name[order(student.last_name)],

nom=sort([Link]),taille=c(1.7,1.85,1.75,1.65))
[Link](demog,"[Link]",[Link]=F)

#File containing names and sizes

demog <- [Link](first_name = student.first_name[order(student.last_name)],
nom=sort([Link]),taille=c(1.7,1.85,1.75,1.65))
[Link](demog,"[Link]",[Link]=F)

Number of students born this month

16
[Link] <- matrix([Link](unlist(strsplit([Link], "-"))), ncol=3, byrow=T)
[Link]<-rep(1,12)
for (i in 1:12) {
nba <- length([Link][[Link][,2] == i, 1])
if(nba>0) cat("Number of birthdays in", [Link][i],":", nba, "\n")
[Link][i]<-nba
}
plot(c(1:12),[Link],type="h")

if(length(unique([Link]))<length([Link]))
At least 2 students were born on the same day

Birth Year Histogram

hist([Link][,3])

Age Calculation
now <- date()
da1<-unlist(strsplit(now," "))
da1<-da1[da1!=""]
mon <- pmatch(da1[2], [Link])
day<-[Link](da1[3])
year<-[Link](da1[5])
age <- year - [Link][,3]
for(i in 1:length(age)) {
if([Link][i,2]>mon | ([Link][i,2]==mon & [Link][i,1]>day))
age[i]<-age[i]-1
}
Next and last birthday of the year
now <- date()
da1<-unlist(strsplit(now," "))
da1<-da1[da1!=""]
[Link] <- [Link](da1[5])
today<-[Link](pmatch(da1[2],[Link]),[Link](da1[3]),[Link])

[Link]<-[Link]([Link][,2],[Link][,1],[Link])
vec<-[Link]
We check that there has been at least one birthday this year.

17
#else we consider the birthdays of the previous year
if(length(vec[vec>=0])<=0)
[Link] <- [Link]([Link][,2], [Link][,1], [Link] - 1)
vec <- today - [Link]

indx<-1:length(vec)
ind1 <- indx[vec == min(vec[vec >= 0])]
The last birthday took place on

Is there at least one birthday that has passed this year?

#Otherwise we consider next year's birthdays
if(length(vec[vec<0])<=0)
[Link] <- [Link]([Link][,2], [Link][,1], [Link] + 1)
vec <- today - [Link]

The next birthday will be in

5.2 Functions
Biased variance

[Link] <- function(x) { sum((x-mean(x))**2)/length(x) else

sum((x-mean(x))**2)/length(x) sum((x-mean(x))**2)/(length(x)-1)
} }
x<-1:100 x <- 1:100
[Link](x) [Link](x)
var(x) var(x)
[Link]<-function(x,biased=F) { [Link](x,biased=True)
if(biased)

Density

phi <- function(x) { phi(0)

exp(-x^2/2) / sqrt(2 * pi) dnorm(0)
}

Function to correct. To debug, run the function while carefully reading the error messages.
(syntax, object not found, warnings), and gradually modify the function.

Original function

18
toto<-myFunction(x,y=true) {
if(y=T) mean(x*a) else print(mean(x*a,[Link]=T)
return(res=sum(x),)
}
#Corrected function
toto <- function(x, a, y = TRUE) {
if(y) print(mean(x*a)) else print(mean(x*a,[Link]=T))
return(res=sum(x))
}

This is the final grade

Original function
[Link] <- function(grades, p) {
netud <- nrow(notes)
neval <- ncol(notes)
final <- (1:netud) * 0
for(i in 1:netud) {
for(j in 1:neval) {
final[i] <- final[i] + notes[i, j] * p[j]
}
}
final
}
Optimized function
final_notes2 <- function(notes, p)
apply(t(notes)*p, 2, sum)
Simulation of a grades and weights table to test the function
ns<-20;np<-5
notes <- matrix(trunc(runif(ns * np, 3, 21)), ncol = np)
pond <- trunc(runif(4, 1, 4)) / 20
pond <- c(pond, 1 - sum(pond))
pond
We check that our two functions return the same result...
[Link](notes,weight)
final_notes2(notes,score)

19
The return of the TCL:

Definition of the function }

extcl<-function(nobs, ntir=1) { Execution
moy<-c() nobs<-c(3,5,10,30)
for(i in 1:ntir) { ntir<-50
vec <- runif(nobs) par(mfrow=c(2,2))
moy<-c(moy,mean(vec)) for (i in nobs)
} hist(extcl(nobs,ntir),breaks=20)
return(average)

We can, but it is much trickier to include the type of distribution and its parameters in the
definition of the function. This calls for a new function,[Link], which accepts as an argument
the name of a function and a list with its arguments. Here, we will give the name of the function as
a string, and use the "..." argument to allow specifying arguments to be given
to this function. WhenRreceives arguments in place of '...', it passes them to the function where
these "..." appear here in [Link].

extcl2<-function(nobs,ntir=1,fonc="rnorm",...) {
moy<-c()
for(i in 1:ntir) {
vec <- [Link](function, list(nobs, ...))
moy<-c(moy,mean(vec))
}
return(moy)
}

We can then use this function with all the known sampling functions in R:

nobs<-c(3,5,10,30)
ntir<-50
par(mfrow=c(2,2))
for (i in nobs)
hist(extcl2(i, ntir, 'runif', 0, 1), xlab=paste(i, 'draws'))
main="Uniform distribution",breaks=20)
par(mfrow=c(2,2))
for (i in nobs)
hist(extcl2(i, ntir, "rexp", 5), xlab=paste(i, "draws"),
main="Exponential distribution",breaks=20)

20
5.3 Advanced Concepts
rpareto <- function(n, alpha, lambda)
lambda * (runif(n)^(-1/alpha) - 1)
malist<-vector(length=5,mode="list")
i1<-1
for(i in c(100,150,200,250,300)) {
malist[[i1]]<-rpareto(i,2,5000)
i1<-i1+1
}
names(malist) <- paste("sample", 1:5, sep="")
vec <- lapply(malist, mean)
pareto <- function(x, alpha, lambda) {
1-exp(alpha*log(lambda/(x+lambda)))
}
dpareto<-function(x,alpha,lambda) {
alpha*(lambda^alpha)/exp((alpha+1)*log(x+lambda))
}
vec1<-sort(ppareto(unlist(malist),2,5000))
hist(malist[[5]])
hist(malist$sample5)
vec2<-unlist(lapply(malist, function(x) sort(ppareto(x, 2,5000))))
vec3 <- lapply(lapply(malist, sort), ppareto, alpha = 2, lambda = 5000)

5.4 Libraries
Combined bookstore:

library(combinat)
x<-c("rouge","orange","jaune","vert","bleu","indigo","violet")
tab<-permn(x)
length(tab)
#note: this number is obviously equal to
factorial(length(x))

tab<-combn(x,4)
dim(tab)
choose(7,4)

21
Library date:

now <- date()

naiss<-[Link](5,25,1971)
today <- [Link](pmatch(da1[2], [Link]), [Link](da1[3]), [Link](da1[5]))
I was born

mon <- pmatch(da1[2], [Link]) + 2

yea<-[Link](da1[5])
while(mon>=13) {
mon<-mon-12
yea <- yea + 1 }
unmois<-paste([Link](da1[3]),mon,yea,sep="/")
cat("In two months, it will be the", unmois, "\n")

da1<-unlist(strsplit(now," "))
da1 <- da1[da1 != ""]
today<-paste(da1[3],da1[2],da1[5])
cat("Today, we are on", today, "\n")
cat(" and it is
heur <- [Link](unlist(strsplit(da1[4], ":")))
heur[1]<-heur[1]+1
if(heur[1]==24) heur[1]<-0
You never know, you might be studying at midnight.
heur2 <- paste(heur, collapse = ":")
cat(" in an hour it will be

22