STAT 2450 Assignment 1 (35 points)
Banner ID, First name, Last name
Intro
Problem 1: Using variables [2 pts] Problem 2: Compute mathematical formulas [5 pts] Problem 3: Matrix operations [10 pts] Problem 4: dataframes and basic plots [10 pts] Problem 5: Loops [8 pts]
Upload two files (one by one) 1. this Rmd file (completed) 2. a pdf print of the knitted report of this file
Do not forget to replace your ID in the header above: author: “Banner ID B009556129 First name kenina Last name wang”
Problem 1 [2 pts]
In this problem, you will use R to calculate the mass of the proton in kilograms (kg). The goal is to express all values of physical quantities in the International System of Units, internationally known by the abbreviation SI.
https://melusine.eu.org/syracuse/immae/mp/physique-chimie/constantes_physiques_fondamentales.pdf
According to Einstein, the following relation holds between the energy of a particle (E), the mass of the particle (m) and the speed of light (c):
\[E=mc^2\]
In nuclear or atomic physics, energy is often expressed in electronvolts. 1 Electronvolt is the energy carried by a charge equal to the unit of charge (defined as the absolute value of the charge of the electron) when this charge passes through a voltage of V=1 volt.
More generally, the energy E of a charge q passing in a voltage of V Volt is:
\[E=qV\]
where q is expressed in Coulomb (symbol: C), which is the unit of charge in the standard system SI of unit (SI). This energy is expressed in Joules (symbol: J), which is the unit of energy in SI.
One electronvolt is therefore equal to E=qx1 in Joule (J) (SI).
You need to write a small R code to compute the mass in kilogram (kg in SI) that corresponds to the energy of the proton. The energy of the proton is 938 MeV, where 1 Mev = 1 million of electronvolts. You simply need to apply the equation of Einstein and solve for the mass of the proton.
Your code should use the following variables: m (mass of the proton), q (charge of the proton), c (speed of light).
Note that the charge of the proton is positive and equal to the unit of charge (which itself is the absolute value of the charge of the electron). Again the charge is expressed in Coulomb (SI symbol: C).
I give you the following numerical values:
\[q_e=1.602176634 \times 10^{-19} \quad C \\ c=299792458 \quad m/s\]
# leave this code unchanged here
m=NA # charge of proton in Coulomb (C)
q=1.602176634e-19
# speed of light in meters per second
c=299792458
# energy of proton in Joule (J)
E=q*938*1e6
# mass of proton in kg
m=E/(c^2)
# print the mass of the proton
print(m)## [1] 1.672137e-27The mass of the proton in kg is 1.6721369^{-27}
Problem 2 [5 pts]
Evaluate a mathematical formula in R.
Consider the following mathematical function:
\[f(x)=\frac{x^{\alpha-1}*e^{-x/\sigma}}{(\alpha-1)! *\sigma^\alpha}\]
Use R to calculate this function for alpha=5 and sigma=2 at the value x=1.3 Use the variable names alpha,sigma,x,fx to store the values of alpha,sigma,x and the value fx of the function at x. After calculating fx, print the value of fx. Hint: you can use the built-in functions factorial and exp (for exponential) plus usual arithmetic operations.
# leave this code unchanged here
fx=NA# define the shape parameter : alpha=
alpha=5
# define the scale parameter: sigma=
sigma=2
# define the value of x
x=1.3
# compute fx
fx=((x^(alpha-1))*exp(-x/sigma))/(factorial(alpha-1)*sigma^alpha)
# print the value of fx
print(fx)## [1] 0.001941426The value of fx is 0.0019414
Note: you can check your result using the built-in R function “dgamma”.
Run the next cell of code to see if this gives the same result that you computed above forfx.
dgamma(x,shape=alpha,scale=sigma)# if you want to see the documentation on function dgamma:
?dgammaProblem 3 [10 pts]
Matrix operations. In this exercise, you need to fill in R code to answer simple questions on matrix algebra in R.
We want to create a matrix named M1 that is equal to:
0 2 1 6
0 0 1 2
0 0 0 3
0 0 0 0# use c() to define a vector V of data row by row
v=c(0,2,1,6,0,0,1,2,0,0,0,3,0,0,0,0)
# use matrix() to define the matrix M1. enter the data by row:
M1=matrix(v,nrow = 4,ncol=4,byrow = TRUE)
# print the matrix M1
M1 # 1 pt## [,1] [,2] [,3] [,4]
## [1,] 0 2 1 6
## [2,] 0 0 1 2
## [3,] 0 0 0 3
## [4,] 0 0 0 0Use the names (M2, Mprod etc.) that are indicated.
# Compute M2 = the transpose of matrix M1
M2=t(M1)
# print matrix M2
print(M2) # 1 pt ## [,1] [,2] [,3] [,4]
## [1,] 0 0 0 0
## [2,] 2 0 0 0
## [3,] 1 1 0 0
## [4,] 6 2 3 0# compute the product of matrices in the sense of linear algebra
# Mprod=M1 . M2 . M1 . M1 do not use * for the multiplication!
Mprod=M1%*%M2%*%M1%*%M1
# print the matrix Mprod
print(Mprod) # 2 pts## [,1] [,2] [,3] [,4]
## [1,] 0 0 82 326
## [2,] 0 0 26 106
## [3,] 0 0 36 144
## [4,] 0 0 0 0# compute the determinant (detMprod) of matrix Mprod
detMprod=det(Mprod)
# print the determinant of Mprod
print(detMprod) # 2 pts## [1] 0# compute the dimension (Mdim) of Mprod
Mdim=dim(Mprod)
# print the dimension Mdim of Mprod
print(Mdim) # 2 pts## [1] 4 4# compute the sum (Vs) of the first row of M1 and the fourth column of M2:
Vs=sum(M1[1,]+M2[,4])
# print the sum Vs
print(Vs) # 2 pts ## [1] 9Problem 4 [10 pts]
On dataframes
Q1 [3 pt]
- Use the read.csv command to read the data from the following URL into a dataframe named d.
- print the names of the columns of dataframe object d
- print the structure of d
# a. read the data filem use the read.csv function [1 pt]
d=read.csv("https://gist.githubusercontent.com/seankross/a412dfbd88b3db70b74b/raw/5f23f993cd87c283ce766e7ac6b329ee7cc2e1d1/mtcars.csv")
# b. print the names of the columns of d [1 pt]
print(names(d))## [1] "model" "mpg" "cyl" "disp" "hp" "drat" "wt" "qsec" "vs"
## [10] "am" "gear" "carb"# c. print the structure of d [1 pt]
str(d)## 'data.frame': 32 obs. of 12 variables:
## $ model: chr "Mazda RX4" "Mazda RX4 Wag" "Datsun 710" "Hornet 4 Drive" ...
## $ mpg : num 21 21 22.8 21.4 18.7 18.1 14.3 24.4 22.8 19.2 ...
## $ cyl : int 6 6 4 6 8 6 8 4 4 6 ...
## $ disp : num 160 160 108 258 360 ...
## $ hp : int 110 110 93 110 175 105 245 62 95 123 ...
## $ drat : num 3.9 3.9 3.85 3.08 3.15 2.76 3.21 3.69 3.92 3.92 ...
## $ wt : num 2.62 2.88 2.32 3.21 3.44 ...
## $ qsec : num 16.5 17 18.6 19.4 17 ...
## $ vs : int 0 0 1 1 0 1 0 1 1 1 ...
## $ am : int 1 1 1 0 0 0 0 0 0 0 ...
## $ gear : int 4 4 4 3 3 3 3 4 4 4 ...
## $ carb : int 4 4 1 1 2 1 4 2 2 4 ...Q2 [2 pt]
Write one line of R code to query the R documentation about the dataframe mtcars. MAKE SURE TO SET eval to F to prevent this cell of code from being knitted!
?mtcarsWhat does the documentation tell you about the format of column number 6 (wt)?
# to answer this question, just paste here the line that you read in the documentation (lower right panel) for column 6 in the Format section in the R documentation.
wt Weight (1000 lbs)Q3 [1 pt]
Use the summary command to get summary statistics for all of the variables in the data set.
summary(mtcars)## mpg cyl disp hp
## Min. :10.40 Min. :4.000 Min. : 71.1 Min. : 52.0
## 1st Qu.:15.43 1st Qu.:4.000 1st Qu.:120.8 1st Qu.: 96.5
## Median :19.20 Median :6.000 Median :196.3 Median :123.0
## Mean :20.09 Mean :6.188 Mean :230.7 Mean :146.7
## 3rd Qu.:22.80 3rd Qu.:8.000 3rd Qu.:326.0 3rd Qu.:180.0
## Max. :33.90 Max. :8.000 Max. :472.0 Max. :335.0
## drat wt qsec vs
## Min. :2.760 Min. :1.513 Min. :14.50 Min. :0.0000
## 1st Qu.:3.080 1st Qu.:2.581 1st Qu.:16.89 1st Qu.:0.0000
## Median :3.695 Median :3.325 Median :17.71 Median :0.0000
## Mean :3.597 Mean :3.217 Mean :17.85 Mean :0.4375
## 3rd Qu.:3.920 3rd Qu.:3.610 3rd Qu.:18.90 3rd Qu.:1.0000
## Max. :4.930 Max. :5.424 Max. :22.90 Max. :1.0000
## am gear carb
## Min. :0.0000 Min. :3.000 Min. :1.000
## 1st Qu.:0.0000 1st Qu.:3.000 1st Qu.:2.000
## Median :0.0000 Median :4.000 Median :2.000
## Mean :0.4062 Mean :3.688 Mean :2.812
## 3rd Qu.:1.0000 3rd Qu.:4.000 3rd Qu.:4.000
## Max. :1.0000 Max. :5.000 Max. :8.000 # (1 pt)Q4 [2 pt]
Make a boxplot
Make a figure showing 3 side-by-side boxplots of the weight of cars (column wt) by groups of cylinders numbers. Make sure to add appropriate labels to the X and Y axes, and add a title to the plot.
In the function boxplot, you can use the argument xlab=““ to add a label for the x axis, ylab=““ to add a label to the y axis, and *main=“” to add a title to your plot.
Use the following values:
Title: Weight vs #cylinders Label of the x axis: Number of cylinders Label of the y axis: Weight (1000 lbs)
boxplot(wt ~ cyl,data=mtcars,
xlab="Number of cylinders",
ylab="Weight (1000 lbs)",
main="Weight vs #cylinders"
)
Q7 [2 pts]
Make a scatterplot of weights of cars (column wt) versus horsepower (column hp).
Title: Cars: weight vs horsepower Label the y-axis: Weight (1000 lbs) Label the x-axis: Horsepower
plot(wt ~ hp,data=mtcars,
xlab="Horsepower",
ylab="Weight (1000 lbs)",
main="weight vs horsepower"
)
Problem 5 [8 pts]
On loops and control flow
Q5.1 [4 pts]
On control flow
Write a small R code that loops over all integer numbers (i) from 1 to 100 and prints the loop index (i) if it is divisible by 12 Your code should have 3 lines.
for(i in (1:100)){
if(i%%12==0){
print(i)
}
}## [1] 12
## [1] 24
## [1] 36
## [1] 48
## [1] 60
## [1] 72
## [1] 84
## [1] 96Q5.1 [4 pts]
On loops.
We define the following suite of numbers:
\[S_n=\sum_{k=0}^n \frac{1}{16^k}( \frac{4}{8k+1} -\frac{2}{8k+4}-\frac{1}{8k+5}-\frac{1}{8k+6} )\] The goal of this problem is to compute \(S_8\).
HINT: you can use the following structure for your code:
Your code should have 8 lines (not including the comment lines that start with a # sign).
# give an initial value to your sum : S
S=0
# initialize an empty vector of errors: e
e=vector("numeric",8)
# start your loop
for(i in (0:7)){
uk=(1/(16^i))*(4/(8*i+1)-2/(8*i+4)-1/(8*i+5)-1/(8*i+6))
# <- add code here
# compute your sum by adding the k-th term of the series to rsum : S=S+uk
S=S+uk
# compute the error at rank n: difference between S and pi
en=abs(S-pi)
# append en to the vector of errors e
e[i+1]=en
# finish your loop
} # <- add code here
# print the error vectors as a column vector. You can use cbind().
print(cbind(e))## e
## [1,] 8.259320e-03
## [2,] 1.701872e-04
## [3,] 5.263243e-06
## [4,] 1.960224e-07
## [5,] 8.129457e-09
## [6,] 3.617053e-10
## [7,] 1.691225e-11
## [8,] 8.202328e-13# print the value of the sum S
print(S)## [1] 3.141593# print the value of pi
pi## [1] 3.141593