Some useful R functions
1. c()
2. print()
3. rep() and seq()
4. apply()
5. for ()
6. [Link](), runif() to generate random numbers from the uniform
distribution
7. length(), mean(), sd(), min(), max(), summary()
8. sum()
9. sample() and replicate()
10. paste()
To read a .csv file in RStudio in CBTF, in a data frame called ‘mydata’:
mydata=[Link](“~/Downloads/[Link]”, header=TRUE)
Numeric functions
1
Some statistical functions
______________________________________________________________________
Using R for one RV integration
Use R function integrate(), to integrate a one-dimensional integral over a finite or
infinite interval,
Example:
# define the integrated function
integrand <- function(x) {1/((x+1)*sqrt(x))}
# integrate the function from 0 to infinity
integrate(integrand, lower = 0, upper = Inf)
res <-integrate(integrand, lower = 0, upper = Inf)
print(res[[1]])
#R outputs res as a list of 5 items, the result is the first
#item. In order to retrieve it from a list you use double
#brackets
x is a vector. Integrand must accept a vector as input and produces a vector of function
evaluations at those points (which are in the range from lower to upper).
For example, if your function is defined as f(x)=2.3 in e.g., range [0,1], writing
function (x) {2.3} #won’t work
But
function(x) {2.3 * x^0} #will work
2
R functions for probability distributions
Recall the meaning of the first letter of the R functions for calculating probabilities:
d signifies probability distribution (gives the density)
p signifies the cumulative probability distribution (gives the distribution function) (for example,
pnorm(q) gives the cumulative normal probability for q, i.e., the area under the normal curve to
the left of q)
q signifies the inverse cumulative (i.e it will return the RV for certain cumulative probability. For
example qnorm(p) gives the normal quantile. i.e., value at the p percentile of the normal
distribution).
r signifies random (gives random numbers taken from the given distribution function)
[Link] = TRUE means P(X≤x) is returned consistent with how the book tables are
compiled.
(Numbers shown below are for example cases, replace as applicable to your problem. For
functions that not explicit arguments are given here, use the help ? utility in RStudio.)
A. CONTINUOUS RVs
Uniform pdf
dunif(x, min = a, max = b, log = FALSE)
punif(n1:n2, min = a, max = b, [Link] = TRUE)
Normal pdf
1 (x )2
f (x; , ) exp
2 2
2
dnorm(x, mean, sd) #gives height of the probability distribution
at each point for a given mean and standard deviation (i.e. PDF)
pnorm(x, mean, sd, [Link]=TRUE) #gives the probability of a
normally distributed random number to be less that the value of a
given number , i.e. the CDF
qnorm(p, mean, sd, [Link]=TRUE)# takes the probability value
and gives a number whose cumulative value matches the probability
value, i.e. the inverse CDF.
3
Student t pdf
dt(x, df, log = FALSE)
pt(q, df, [Link] = TRUE, log.p = FALSE)
qt(p, df, [Link] = TRUE, log.p = FALSE)
rt(n, df)
Lognormal pdf
If W is normally distributed with mean θ and standard deviation ω, and X = exp(W)
1
1 (ln( x ) )2
f (x) e 2 2
x 2
dlnorm(x, meanlog, sdlog, log = FALSE)
plnorm(q, meanlog, sdlog, [Link] = TRUE, log.p = FALSE)
qlnorm(p, meanlog, sdlog, [Link] = TRUE, log.p = FALSE)
Be careful because meanlog and sdlog refer to the parameters of the normal distribution. The
meaning of the first letter of each function is the same as in all other probability functions mentioned
in previous worksheets. Always check help in R by typing ?functionname in R prompt.
Gamma pdf
1
f (x;, ) x 1e x/ for x>0, and 0 elsewhere.
( )
=; 2 = 2
(n) = (n-1)! for integers; more complex for other real-valued numbers
shape = α rate = 1 / β
dgamma(x, shape, rate, log = FALSE)
pgamma(q, shape, rate, [Link] = TRUE, log.p = FALSE)
qgamma(p, shape, rate, [Link] = TRUE, log.p = FALSE)
rgamma(n, shape, rate)
Exponential pdf
The special gamma distribution for which α = 1 is the exponential distribution.
4
Chi-Squared pdf
The special gamma distribution for which α = ν/2 and β = 2 is the Chi-squared distribution,
where v degrees of freedom.
B. DISCRETE RVs
Uniform pmf
f(xi) = 1/n = 1/(b-a+1)
dunif(x, min = a, max = b, log = FALSE)
punif(q, min = a, max = b, [Link] = TRUE) #gives P(X≤q)
Binomial pmf
n x
b(x;n, p) p (1 p)nx np 2 np(1 p) npq
x
dbinom(x=3, size=20, prob=0.25) #P(X=3)
dbinom(x=0:5, size=20, prob=0.25) #will give probabilities for
#multiple values of x
pbinom(q=5, size =20, prob=0.25, [Link] = T) # P(X ≤ 5)
Hypergeometric pmf
k N k
k k k N n
h(x; N, n, k)
x n x n 2 n 1
N N N N N 1
n
dhyper(x,k,(N-k),n)
phyper(x,k,(N-k),n, [Link]=TRUE)
Negative Binomial pmf
k
x 1 k p
b * (x;k, p) p (1 p) xk
k 1 x ≥ k
k(1 p)
2
p2
dnbinom((x-k),k,p)
5
pnbinom((x-k),k,p)
Geometric pmf
1 (1 p)
g(x; p) p(1 p) x1 2
p p2
dgeom((x-1),p)
pgeom((x-1),p)
Poisson pmf
e t ( t) x t 2 t
p(x; t)
x!
(Remember what lambda represents in the R function for the Poisson
distribution)
dpois(x,λt)
ppois(x,λt,[Link]=TRUE)
______________________________________________________________________
Simple linear regression
Now use function lm to fit a linear regression model.
lmodel = lm(varY ~ varX) #notice the dependent variable goes first
A summary of this model:
print(summary (lmodel))
attributes(lmodel)
You obtain headings (descriptions) of the data stored in lmodel. You may pull out and see some of
these data by typing (for example)
lmodel$coefficients # (lmodel$coef does work too)
Alternatively, you may type
coef(lmodel)
and you can also get confidence intervals of these coefficients
confint(lmodel, level=0.95) #level is what confidence level you choose
If you want to add a line to your plot you can do the following
6
plot (varX, varY, main=’meaningful title’, xlab = ‘varX (units)’,
ylab=’varY (units)’)
abline(lmodel) #type help(?abline)in the command window to see
#options for changing line properties such as color
#and line width
Testing for correlation R is not part of the output of linear regression in R. You can use:
[Link]( ~ varY + vayX, method = "pearson", [Link] = 0.95)
#use the appropriate names of your variables.
Checking the assumptions of the model
par(mfrow=c(2,2))
plot(lmodel)