# Math Coding

## How to do Math Coding with Example

x<-seq(0,10,0.2)

> par(mfrow=c(1,2))

```> y <- exp(x)
> plot (x~y,type='l',main="Exponential")
> y <- log(x)
> plot(y~x,type="l",main="Logarithmic")
> plot (y~x,type='l',main="Exponential")```
```> y <- exp(x)
> plot (y~x,type='l',main="Exponential")```
```
> #Trigonometric functions
> windows(7,7)
> par(mfrow=c(2,2))
> x <- seq(0,2*pi,2*pi/100)
> y1 <- cos(x)
> y2 <- sin(x)
> y3 <- tan(x)
> plot(y1~x,type="l",main="cosine")```

```
> plot(y2~x,type="l",main="sine")
> plot(y1~x,type="l",main="tangent")
> plot(y3~x,type="l",ylim=c(-3,3),main="tangent")
> plot(y3~x,type="l",main="tangent")```
```
> #power laws
> x<-seq(0,10,.01)
> y
[1]     1.000000     1.221403     1.491825
[4]     1.822119     2.225541     2.718282
[7]     3.320117     4.055200     4.953032
[10]     6.049647     7.389056     9.025013
[13]    11.023176    13.463738    16.444647
[16]    20.085537    24.532530    29.964100
[19]    36.598234    44.701184    54.598150
[22]    66.686331    81.450869    99.484316
[25]   121.510418   148.413159   181.272242
[28]   221.406416   270.426407   330.299560
[31]   403.428793   492.749041   601.845038
[34]   735.095189   897.847292  1096.633158
[37]  1339.430764  1635.984430  1998.195895
[40]  2440.601978  2980.957987  3640.950307
[43]  4447.066748  5431.659591  6634.244006
[46]  8103.083928  9897.129059 12088.380730
[49] 14764.781566 18033.744928 22026.465795
> y=x^0.5
> plot(x,y,type="l",main="0<b<1")
> y <- x
> plot(x,y,type="l",main="b=1")
> y <- x^2
> plot(x,y,type="l",main="b>1")
> y <- 1/x
> plot(x,y,type="l",main="b<0")```

```# polynomial function > x<-seq(0,10,.01)
> y1 <- 2+5*x-0.2*x^2
> y2 <- 2+5*x-0.4*x^2
> y3 <- 2+4*x-0.6*x^2+0.04*x^3
> y4<- 2+4*x+2*x^2-0.6*x^3+0.04*x^4
> par(mfrow=c(2,2))
> plot(x,y1,type="l",ylab="y",main="decelerating")
> plot(x,y2,type="l",ylab="y",main="humped")
> plot(x,y3,type="l",ylab="y",main="inflection")
> plot(x,y4,type="l",ylab="y",main="local maximum")> y <- c(8,3,5,7,6,6,8,9,2,3,9,4,10,4,11)
> mean(y)
[1] 6.333333
> median(y)
[1] 6
> mode(y)
[1] "numeric"
> range(y)
[1]  2 11
> var(y)
[1] 7.809524
> quantile(y)
0%  25%  50%  75% 100%
2.0  4.0  6.0  8.5 11.0
> cumsum(y)
[1]  8 11 16 23 29 35 43 52 54 57 66 70 80 84 95
> colMeans(y)
Error in colMeans(y) : 'x' must be an array of at least two dimensions
> counts <- rnbinom(10000,mu=0.92,size=1.1)
> counts[1:30]
[1] 1 2 0 1 0 4 1 0 0 0 4 1 1 1 1 0 0 1 0 0 0 3 0
[24] 0 3 1 0 0 0 0
> X <- matrix(c(1,0,0,0,1,0,0,0,1),nrow=3)
> X
[,1] [,2] [,3]
[1,]    1    0    0
[2,]    0    1    0
[3,]    0    0    1
> class(X)
[1] "matrix" "array"
> attributes(X)
\$dim
[1] 3 3

> vector <- c(1,2,3,4,4,3,2,1)
> V <- matrix(vector,byrow=T,nrow=2)
> V
[,1] [,2] [,3] [,4]
[1,]    1    2    3    4
[2,]    4    3    2    1
> dim(vector) <- c(4,2)
> dim(vector)
[1] 4 2
> is.matrix(vector)
[1] TRUE
> vector
[,1] [,2]
[1,]    1    4
[2,]    2    3
[3,]    3    2
[4,]    4    1
> (vector <- t(vector))
[,1] [,2] [,3] [,4]
[1,]    1    2    3    4
[2,]    4    3    2    1Differential Equationdetach("package:Deriv", unload=TRUE)
> Lorenz <- function(t, state, parameters) {
+     with(as.list(c(state, parameters)), {
+         dX <-  a * X + Y * Z
+         dY <-  b * (Y - Z)
+         dZ <- -X * Y + c * Y - Z
+         list(c(dX, dY, dZ))
+     })
+ }
> parameters <- c(a = -8/3, b = -10, c = 28)
> state      <- c(X = 1, Y = 1, Z = 1)
> times      <- seq(0, 100, by = 0.01)
> out <- ode(y = state, times = times, func = Lorenz, parms = parameters)
> plot(out)#Maximum likelihood with the normal distribution
> par(mfrow=c(2,2))
> curve(dnorm,-3,3,xlab="z",ylab="Probability density",main="Density")
> curve(pnorm,-3,3,xlab="z",ylab="Probability",main="Probability")
> curve(qnorm,0,1,xlab="p",ylab="Quantile (z)",main="Quantiles")
> y <- rnorm(1000)
> hist(y,xlab="z",ylab="frequency",main="Random numbers")
> yvals <- rnorm(100,24,4)
> mean(yvals)
[1] 23.98589
> sd(yvals)
[1] 3.395763
> ydevs <- rnorm(100,0,1)
> ydevs <- (ydevs-mean(ydevs))/sd(ydevs)
> mean(ydevs)
[1] -1.587923e-17
> sd(ydevs)
[1] 1
> yvals <- 24 + ydevs*4
> mean(yvals)
[1] 24
> sd(yvals)
[1] 4
> windows(7,4)
> par(mfrow=c(1,2))
> x <- seq(0,30,.25)
> plot(x,pchisq(x,3,7.25),type="l",ylab="p(x)",xlab="x")
> plot(x,pchisq(x,5,10),type="l",ylab="p(x)",xlab="x")
> 8*10.2/qchisq(.975,8)
[1] 4.65367
> 8*10.2/qchisq(.025,8)
[1] 37.43582
> qf(.95,2,18)
[1] 3.554557
> x <- seq(0.05,4,0.05)
> plot(x,df(x,2,18),type="l",ylab="f(x)",xlab="x")
> plot(x,df(x,6,18),type="l",ylab="f(x)",xlab="x")
> windows(7,7)
> par(mfrow=c(1,1))
> df <- seq(1,30,.1)
> plot(df,qf(.95,df,30),type="l",ylab="Critical F")
> lines(df,qf(.95,df,10),lty=2)
> x <- seq(0.01,3,0.01)
> plot(x,df(x,1,10),type="l",ylim=c(0,1),ylab="f(x)")
> lines(x,df(x,2,10),lty=6,col="red")
> lines(x,df(x,5,10),lty=2,col="green")
> lines(x,df(x,30,10),lty=3,col="blue")
> legend(2,0.9,c("1","2","5","30"),col=(1:4),lty=c(1,6,2,3),
+        title="numerator d.f.")
> curve( (1+xË†2)Ë†(-0.5), -3, 3,ylab="t(x)",col="red")
Error: unexpected input in "curve( (1+xË†"
> curve( (1+x^2)^(-0.5), -3, 3,ylab="t(x)",col="red")
> plot(1:30,qt(0.975,1:30), ylim=c(0,12),type="l",
+      ylab="Students t value",xlab="d.f.",col="red")
> abline(h=2,lty=2,col="green")
> xvs <- seq(-4,4,0.01)
> plot(xvs,dnorm(xvs),type="l",lty=2,
+     ylab="Probability density",xlab="Deviates")
> lines(xvs,dt(xvs,df=5),col="red")
> qt(0.975,5)
[1] 2.570582
> #The gamma distribution
> x <- seq(0.01,4,.01)
> par(mfrow=c(2,2))
> y <- dgamma(x,.5,.5)
> plot(x,y,type="l",col="red",main="alpha = 0.5")
> y <- dgamma(x,.8,.8)
> plot(x,y,type="l",col="red", main="alpha = 0.8")
> x <- seq(0,1,0.01)
> fx <- dbeta(x,2,3)
> plot(x,fx,type="l",main="a=2 b=3",col="red")
> fx <- dbeta(x,0.5,2)
> plot(x,fx,type="l",main="a=0.5 b=2",col="red")
> fx <- dbeta(x,2,0.5)
> plot(x,fx,type="l",main="a=2 b=0.5",col="red")
> fx <- dbeta(x,0.5,0.5)
> plot(x,fx,type="l",main="a=0.5 b=0.5",col="red")
> rbeta(10,2,3)
[1] 0.10678527 0.46414831 0.09402682 0.74723789 0.43859319
[6] 0.49928575 0.83196791 0.13017362 0.17905408 0.61618885
> windows(7,4)
> par(mfrow=c(1,2))
> plot(-200:200,dcauchy(-200:200,0,10),type="l",ylab="p(x)",xlab="x",
+      col="red")
> plot(-200:200,dcauchy(-200:200,0,50),type="l",ylab="p(x)",xlab="x",
+      col="red")
> windows(7,7)
> plot(seq(0,10,0.05),dlnorm(seq(0,10,0.05)),
+      type="l",xlab="x",ylab="LogNormal f(x)",col="x")
Error in plot.xy(xy, type, ...) : invalid color name 'x'
> plot(seq(0,10,0.05),dlnorm(seq(0,10,0.05)),
+ type="l",xlab="x",ylab="LogNormal f(x)",col="x")
Error in plot.xy(xy, type, ...) : invalid color name 'x'
> windows(7,7)
> a <- 3
> l <- 1
> t <- seq(0,1.8,.05)
> ft <- a*l*tË†(a-1)*exp(-l*tË†a)
Error: unexpected input in "ft <- a*l*tË†"
> ft <- a*l*t^(a-1)*exp(-l*t^a)
> plot(t,ft,type="l",col="blue",ylab="f(t) ")
> a <- 1
> ft <- a*l*t^(a-1)*exp(-l*t^a)
> lines(t,ft,type="l",col="red")
> a <- 2
> ft <- a*l*t^(a-1)*exp(-l*t^a)
> lines(t,ft,type="l",col="green")
> legend(1.4,1.1,c("1","2","3"),title="alpha",lty=c(1,1,1),col=c(2,3,4))
> xy <- mvrnorm(1000,mu=c(50,60),matrix(c(4,3.7,3.7,9),2))
Error in mvrnorm(1000, mu = c(50, 60), matrix(c(4, 3.7, 3.7, 9), 2)) :
could not find function "mvrnorm"
> x <- ceiling(runif(10000)*6)
> table(x)
x
1    2    3    4    5    6
1702 1693 1639 1684 1591 1691
> hist(x,breaks=0.5:6.5,main="")
> # Matrix algebra
> a <- matrix(c(1,0,4,2,-1,1),nrow=3)
> a
[,1] [,2]
[1,]    1    2
[2,]    0   -1
[3,]    4    1
> b <- matrix(c(1,-1,2,1,1,0),nrow=2)
> b
[,1] [,2] [,3]
[1,]    1    2    1
[2,]   -1    1    0
> a[1,]
[1] 1 2
> a[1,]*b[,1]
[1]  1 -2
> sum(a[1,]*b[,1])
[1] -1
> a[1,]
[1] 1 2
> b[,2]
[1] 2 1
> a[1,]*b[,2]
[1] 2 2
> sum(a[1,]*b[,2])
[1] 4
> a[1,]*b[,3]
[1] 1 0
> sum(a[1,]*b[,3])
[1] 1
> a %*% b
[,1] [,2] [,3]
[1,]   -1    4    1
[2,]    1   -1    0
[3,]    3    9    4
> b %*% a
[,1] [,2]
[1,]    5    1
[2,]   -1   -3
> #Diagonals of matrices
> (ym <- diag(1,3,3))
[,1] [,2] [,3]
[1,]    1    0    0
[2,]    0    1    0
[3,]    0    0    1
> diag(ym) <- 1:3
> ym
[,1] [,2] [,3]
[1,]    1    0    0
[2,]    0    2    0
[3,]    0    0    3
> diag(ym)
[1] 1 2 3
> M <- cbind(X=1:5, Y=rnorm(5))
> M
X           Y
[1,] 1 -0.25561312
[2,] 2  0.48569734
[3,] 3 -0.92566107
[4,] 4 -2.07057266
[5,] 5 -0.04476449
> var(M)
X          Y
X  2.5000000 -0.5336432
Y -0.5336432  0.9667790
> diag(var(M))
X        Y
2.500000 0.966779
> #Determinant
> A <- matrix(c(1,2,4,2,1,1,3,1,2),nrow=3)
> A
[,1] [,2] [,3]
[1,]    1    2    3
[2,]    2    1    1
[3,]    4    1    2
> det(A)
[1] -5
> B <- A
> B[3,] <- 3*B[3,]
> B
[,1] [,2] [,3]
[1,]    1    2    3
[2,]    2    1    1
[3,]   12    3    6
> det(B)
[1] -15
> C <- A
> C[,2] <- 0
> C
[,1] [,2] [,3]
[1,]    1    0    3
[2,]    2    0    1
[3,]    4    0    2
> det(C)
[1] 0
> library(MASS)
> MASS
> ginv(A)
[,1] [,2] [,3]
[1,] -2.000000e-01  0.2  0.2
[2,] -5.828671e-16  2.0 -1.0
[3,]  4.000000e-01 -1.4  0.6
> ginv(ginv(A))
[,1] [,2] [,3]
[1,]    1    2    3
[2,]    2    1    1
[3,]    4    1    2
> 1/det(ginv(A))
[1] -5
> L <- c(0,0.7,0,0,6,0,0.5,0,3,0,0,0.3,1,0,0,0)
> L <- matrix(L,nrow=4)
> L
[,1] [,2] [,3] [,4]
[1,]  0.0  6.0  3.0    1
[2,]  0.7  0.0  0.0    0
[3,]  0.0  0.5  0.0    0
[4,]  0.0  0.0  0.3    0
> n <- c(45,20,17,3)
> n <- matrix(n,ncol=1)
> n
[,1]
[1,]   45
[2,]   20
[3,]   17
[4,]    3
> L %*% n
[,1]
[1,] 174.0
[2,]  31.5
[3,]  10.0
[4,]   5.1
> 45*0+20*6+17*3+3*1
[1] 174
> fun <- function(x) L %*% x
> n <- c(45,20,17,3)
> n
[1] 45 20 17  3
> n <- matrix(n,ncol=1)
> n
[,1]
[1,]   45
[2,]   20
[3,]   17
[4,]    3
> structure <- numeric(160)
> dim(structure) <- c(40,4)
> for (i in 1:40) {
+     n <- fun(n)
+     structure[i,] <- n
+ }
> matplot(1:40,log(structure),type="l")
> sum(structure[40,])/sum(structure[39,])
[1] 2.164035
> structure[40,]/sum(structure[40,])
[1] 0.709769309 0.230139847 0.052750539 0.007340305
> eigen(L)
eigen() decomposition
\$values
[1]  2.1694041+0.0000000i -1.9186627+0.0000000i
[3] -0.1253707+0.0975105i -0.1253707-0.0975105i

\$vectors
[,1]           [,2]                    [,3]
[1,] 0.949264118+0i -0.93561508+0i -0.01336028-0.03054433i
[2,] 0.306298338+0i  0.34134741+0i -0.03616819+0.14241169i
[3,] 0.070595039+0i -0.08895451+0i  0.36511901-0.28398118i
[4,] 0.009762363+0i  0.01390883+0i -0.87369452+0.00000000i
[,4]
[1,] -0.01336028+0.03054433i
[2,] -0.03616819-0.14241169i
[3,]  0.36511901+0.28398118i
[4,] -0.87369452+0.00000000i

> eigen(L)\$vectors[,1]/sum(eigen(L)\$vectors[,1])
[1] 0.710569659+0i 0.229278977+0i 0.052843768+0i
[4] 0.007307597+0i
> A <- matrix(c(3,1,4,2),nrow=2)
> A
[,1] [,2]
[1,]    3    4
[2,]    1    2
> kv <- matrix(c(12,8),nrow=2)
> kv
[,1]
[1,]   12
[2,]    8
> solve(A,kv)
[,1]
[1,]   -4
[2,]    6
> integrate(dnorm,0,Inf)
0.5 with absolute error < 4.7e-05
> integrate(dnorm,-Inf,Inf)
1 with absolute error < 9.4e-05
> integrate(function(x) rep(2, length(x)), 0, 1)
2 with absolute error < 2.2e-14
> integrand <- function(x) {1/((x+1)*sqrt(x))}
> integrate(integrand, lower = 0, upper = Inf)
3.141593 with absolute error < 2.7e-05
> xv <- seq(0,10,0.1)
> plot(xv,integrand(xv),type="l")
> phmodel <- function(t,state,parameters){
+     with(as.list(c(state,parameters)),{
+         dv <- r*v*(K-v)/K-b*v*n
+         dn <- c*v*n-d*n
+         result <- c(dv,dn)
+         list(result)
+     })}
> times <- seq(0,500,length=501)
> parameters <- c(r=0.4,K=1000,b=0.02,c=0.01,d=0.3)
> initial <- c(v=50,n=10)
> output <- ode(y=initial,time=times,func=phmodel,parms=parameters)
time        v        n
[1,]    0 50.00000 10.00000
[2,]    1 58.29220 12.75106
[3,]    2 62.99695 17.40172
[4,]    3 60.70065 24.09264
[5,]    4 50.79407 31.32860
[6,]    5 37.68312 36.12636
> plot(output[,1],output[,2],
+      ylim=c(0,60),type="n",ylab="abundance",xlab="time")
> lines(output[,1],output[,2],col="green")
> lines(output[,1],output[,3],col="red")
> plot(output[,3],output[,2],
+      ylim=c(0,70),xlim=c(0,70),type="n",ylab="plant",xlab="herbivore")
> lines(output[,2],output[,3],col="red")```