Math Coding

How to do Math Coding with Example

mathcoding


 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    1

Differential Equation

detach("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 Error: object 'MASS' not found > head(MASS) Error in head(MASS) : object 'MASS' not found > 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) > head(output) 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")

Post a comment

0 Comments