Mathematical tools in R

rm(list=ls())

Arithmetic operators

2+1
## [1] 3
4-2
## [1] 2
2*3
## [1] 6
1/5
## [1] 0.2
2^3
## [1] 8

Equality testing “==”

4-2==2
## [1] TRUE
1==.9999999999999
## [1] FALSE

Objects (or variables)assignment “=”

a=2
b=5
a
## [1] 2
b
## [1] 5
c=3;c
## [1] 3
a+b
## [1] 7
a/b
## [1] 0.4

Relational

a>b
## [1] FALSE
a<b
## [1] TRUE

Assignment “<-” also works

c<-4
a*c<b*c
## [1] TRUE
c=-4
a*c<b*c
## [1] FALSE

Number of digits on an object

round(1/3,4) # 4 digits
## [1] 0.3333
round(2*4.7536,2) # 2 digits
## [1] 9.51

A number raised to the power of zero is equal to one

2^0
## [1] 1
223^0
## [1] 1

Infinity

2/0 # Positive inf
## [1] Inf
-10/0 # Negative inf
## [1] -Inf

x^Inf when 0 < x < 1 is equal to zero

0.8^Inf
## [1] 0
2^Inf
## [1] Inf

Scientific notation

a=510000 ;a
## [1] 510000
b=0.00000034 ;b
## [1] 3.4e-07
a*b
## [1] 0.1734
a/b
## [1] 1.5e+12
options(scipen=15)
a/b
## [1] 1500000000000
b/a
## [1] 0.0000000000006666667
options(scipen=F)
b/a
## [1] 6.666667e-13

Logarithms

ls()
## [1] "a" "b" "c"
rm(list=ls())

x=c() # An empty object
x[1] = 1          # x(1) is 1 
for (i in 1:7) {  # Loop i from 1 to 7
  x[i]=10^(i-1)   # x(i) 
}
options(scipen=999)
x
## [1]       1      10     100    1000   10000  100000 1000000
cbind(x,log(x))
##            x          
## [1,]       1  0.000000
## [2,]      10  2.302585
## [3,]     100  4.605170
## [4,]    1000  6.907755
## [5,]   10000  9.210340
## [6,]  100000 11.512925
## [7,] 1000000 13.815511
table=cbind(x,log(x));table
##            x          
## [1,]       1  0.000000
## [2,]      10  2.302585
## [3,]     100  4.605170
## [4,]    1000  6.907755
## [5,]   10000  9.210340
## [6,]  100000 11.512925
## [7,] 1000000 13.815511
colnames(table)[2]="log(x)"
table
##            x    log(x)
## [1,]       1  0.000000
## [2,]      10  2.302585
## [3,]     100  4.605170
## [4,]    1000  6.907755
## [5,]   10000  9.210340
## [6,]  100000 11.512925
## [7,] 1000000 13.815511
log(1000*10000)
## [1] 16.1181
log(1000)+log(10000)
## [1] 16.1181
log(1000*10000)==log(1000)+log(10000)
## [1] TRUE
a=log(1000*10000)
a
## [1] 16.1181
exp(a)
## [1] 10000000
curve(log(x),from=0,to=100, main="Natural logarithm")

curve(log(x),from=0,to=1000000)

Defining a function

f <- function(x) { log(x) }
f
## function(x) { log(x) }
f(100)
## [1] 4.60517
log(100)
## [1] 4.60517

Decimal and percentages

y=c(3,3.02)
y
## [1] 3.00 3.02
diff(y)
## [1] 0.02
diff(y)/y[1]
## [1] 0.006666667
100*diff(y)/y[1]
## [1] 0.6666667
lny=log(y)
100*diff(lny)
## [1] 0.6644543
y=seq(1,1.25,0.05)
y=c(1,y[1]+0.01,y[-1])
y
## [1] 1.00 1.01 1.05 1.10 1.15 1.20 1.25
percent=100*(y-y[1])/y[1]
logdiff=100*(log(y)-log(y[1]))
approx.error=100*((percent-logdiff)/logdiff)

table=cbind(y,percent,logdiff,approx.error)
table[-1,]
##         y percent    logdiff approx.error
## [1,] 1.01       1  0.9950331    0.4991708
## [2,] 1.05       5  4.8790164    2.4796716
## [3,] 1.10      10  9.5310180    4.9205869
## [4,] 1.15      15 13.9761942    7.3253544
## [5,] 1.20      20 18.2321557    9.6962990
## [6,] 1.25      25 22.3143551   12.0355029

A linear relationship

rm(list=ls())
f <- function(x) { 1+1*x }

Intercept

f(0)
## [1] 1
par(pty="s") # set the aspect ratio in the plot to be square
curve(f(x), -3,10, main="A linear relationship")
abline(h=0,v=0, col="red")

x1=2
x1;f(x1)
## [1] 2
## [1] 3
points(x1,f(x1),pch=20, col="blue")

x2=4
x2;f(x2)
## [1] 4
## [1] 5
points(x2,f(x2),pch=20, col="green")

Slope

m=(f(x2)-f(x1))/(x2-x1)
m
## [1] 1
rise=f(x2)-f(x1);rise
## [1] 2
run=x2-x1;run
## [1] 2
rise/run
## [1] 1

Slope as the derivative

D(expression(1+1*x), "x")
## [1] 1
curve(1*x/f(x),0,10,  main = "Elasticity of y=1+1x")
abline(v=2, col="red")
abline(h=2/f(2), col="blue")

curve(1*x/f(x),0,10,  main = "Elasticity of y=1+1x")
arrows(2,0,2,2/f(2), col="blue")
arrows(2,2/f(2),-0.3,2/f(2), col="green")

Nonlinear relationships

rm(list=ls())

Derivative rule 1

D(expression(c), "x")
## [1] 0

Derivative rule 2

D(expression(x^n), "x")
## x^(n - 1) * n

Derivative rule 4

D(expression(c*x^n), "x")
## c * (x^(n - 1) * n)

Derivative rule 7

D(expression(exp(x)), "x")
## exp(x)
D(expression(exp(a*x+b)), "x")
## exp(a * x + b) * a

Derivative rule 8

D(expression(log(x)), "x")
## 1/x
D(expression(log(a*x+b)), "x")
## a/(a * x + b)

Example A.1

D(expression(4*x+1), "x")
## [1] 4

Example A.2

curve(x^2-8*x+16, from=-1, to=10)

D(expression(x^2-8*x+16), "x")
## 2 * x - 8

The slope (derivative) of the y=log(x) function

D(expression(log(x)), "x")
## 1/x

The slope dy/dx as a function

d <- function(x) { 1/x }
x <- 1:50
d(x)
##  [1] 1.00000000 0.50000000 0.33333333 0.25000000 0.20000000 0.16666667
##  [7] 0.14285714 0.12500000 0.11111111 0.10000000 0.09090909 0.08333333
## [13] 0.07692308 0.07142857 0.06666667 0.06250000 0.05882353 0.05555556
## [19] 0.05263158 0.05000000 0.04761905 0.04545455 0.04347826 0.04166667
## [25] 0.04000000 0.03846154 0.03703704 0.03571429 0.03448276 0.03333333
## [31] 0.03225806 0.03125000 0.03030303 0.02941176 0.02857143 0.02777778
## [37] 0.02702703 0.02631579 0.02564103 0.02500000 0.02439024 0.02380952
## [43] 0.02325581 0.02272727 0.02222222 0.02173913 0.02127660 0.02083333
## [49] 0.02040816 0.02000000
plot(x,d(x),type="l")

Elasticity = slope * x/y of the log(x) function

e <- function(x) { (1/x)*(x/log(x)) }
plot(x,e(x),type="l")

e(1)
## [1] Inf
e(2)
## [1] 1.442695
curve((1/x)*(x/log(x)), from=0, to=5)

Example A.2

rm(list=ls())
curve(x^2-8*x+16, from=-1, to=10)

D(expression(x^2-8*x+16), "x")
## 2 * x - 8
f <- function(x) { x^2-8*x+16 }
dy <- function(x) { (2*x-8) }
e <- function(x) { (2*x-8)*(x/f(x)) }

x=seq(0,8,2)
f(x)
## [1] 16  4  0  4 16
dy(x)
## [1] -8 -4  0  4  8
e(x)
## [1]   0  -2 NaN   6   4
cbind(x,f(x),dy(x),e(x))
##      x          
## [1,] 0 16 -8   0
## [2,] 2  4 -4  -2
## [3,] 4  0  0 NaN
## [4,] 6  4  4   6
## [5,] 8 16  8   4
curve(e(x), 0,8)

Partial derivative

D(expression(a*x^2+b*x+c*z+d), "x")
## a * (2 * x) + b

Area under the curve

curve(2*x,0,1)

integrand <- function(x) { 2*x }
integrate(integrand, lower = 0, upper = 1)
## 1 with absolute error < 0.000000000000011
integrate(integrand, lower = 0.2, upper = 0.6)
## 0.32 with absolute error < 0.0000000000000036
0.6^2-0.2^2
## [1] 0.32

We would like to shade the region represented by P(0.2 < X < 0.6). The first vertex of our polygon is (0.2,0).

cord.x <- c(0.2)
cord.y <- c(0)

The second vertex is (0.2,f(0.2)), where f(0.2) is the value on the y axis of 2*x evaluated at 0.2.

cord.x <- c(cord.x,0.2) 
cord.y <- c(cord.y,integrand(0.2))

The third and fourth vertices is (0.6,f(0.6)) and (f(0.6),0).

cord.x <- c(cord.x,0.6,0.6)
cord.y <- c(cord.y,integrand(0.6),0)

curve(2*x,0,1)
polygon(cord.x,cord.y,col='aquamarine3')

num = integrate(integrand, lower = 0.2, upper = 0.6) ; num
## 0.32 with absolute error < 0.0000000000000036
str(num)
## List of 5
##  $ value       : num 0.32
##  $ abs.error   : num 0.00000000000000355
##  $ subdivisions: int 1
##  $ message     : chr "OK"
##  $ call        : language integrate(f = integrand, lower = 0.2, upper = 0.6)
##  - attr(*, "class")= chr "integrate"
text(0.4, 0.4, paste(num$value))