How to do Math Code R 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 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")SummaryMath code is providing to write code in mathematics.
0 Comments