R is free software

 R is free software and comes with ABSOLUTELY NO WARRANTY.

You are welcome to redistribute it under certain conditions.

Type 'license()' or 'licence()' for distribution details.

R code


R is a collaborative project with many contributors.

Type 'contributors()' for more information and

'citation()' on how to cite R or R packages in publications.


Type 'demo()' for some demos, 'help()' for on-line help, or

'help.start()' for an HTML browser interface to help.

Type 'q()' to quit R.


> library("Deriv", lib.loc="~/R/win-library/3.6")

Warning message:

package ‘Deriv’ was built under R version 3.6.3 

> D(expression(2*xˆ3),"x")

Error: unexpected input in "D(expression(2*xˆ"

> D(expression(log(x)),"x")

1/x

> D(expression(a*exp(-b * x)),"x")

-(a * (exp(-b * x) * b))

> D(expression(a/(1+b*exp(-c * x))),"x")

a * (b * (exp(-c * x) * c))/(1 + b * exp(-c * x))^2

> trig.exp <- expression(sin(cos(x + yˆ2)))

Error: unexpected input in "trig.exp <- expression(sin(cos(x + yˆ"

> trig.exp <- expression(sin(cos(x + y^2)))

> D(trig.exp,"x")

-(cos(cos(x + y^2)) * sin(x + y^2))

> D(expression(2*x^3),"x")

2 * (3 * x^2)

> library("deSolve", lib.loc="~/R/win-library/3.6")

> x <- seq(0,10,0.1)

> x

  [1]  0.0  0.1  0.2  0.3  0.4  0.5  0.6  0.7  0.8  0.9  1.0

 [12]  1.1  1.2  1.3  1.4  1.5  1.6  1.7  1.8  1.9  2.0  2.1

 [23]  2.2  2.3  2.4  2.5  2.6  2.7  2.8  2.9  3.0  3.1  3.2

 [34]  3.3  3.4  3.5  3.6  3.7  3.8  3.9  4.0  4.1  4.2  4.3

 [45]  4.4  4.5  4.6  4.7  4.8  4.9  5.0  5.1  5.2  5.3  5.4

 [56]  5.5  5.6  5.7  5.8  5.9  6.0  6.1  6.2  6.3  6.4  6.5

 [67]  6.6  6.7  6.8  6.9  7.0  7.1  7.2  7.3  7.4  7.5  7.6

 [78]  7.7  7.8  7.9  8.0  8.1  8.2  8.3  8.4  8.5  8.6  8.7

 [89]  8.8  8.9  9.0  9.1  9.2  9.3  9.4  9.5  9.6  9.7  9.8

[100]  9.9 10.0

> y1 <- 2+5*x-0.2*x^2

> y1

  [1]  2.000  2.498  2.992  3.482  3.968  4.450  4.928

  [8]  5.402  5.872  6.338  6.800  7.258  7.712  8.162

 [15]  8.608  9.050  9.488  9.922 10.352 10.778 11.200

 [22] 11.618 12.032 12.442 12.848 13.250 13.648 14.042

 [29] 14.432 14.818 15.200 15.578 15.952 16.322 16.688

 [36] 17.050 17.408 17.762 18.112 18.458 18.800 19.138

 [43] 19.472 19.802 20.128 20.450 20.768 21.082 21.392

 [50] 21.698 22.000 22.298 22.592 22.882 23.168 23.450

 [57] 23.728 24.002 24.272 24.538 24.800 25.058 25.312

 [64] 25.562 25.808 26.050 26.288 26.522 26.752 26.978

 [71] 27.200 27.418 27.632 27.842 28.048 28.250 28.448

 [78] 28.642 28.832 29.018 29.200 29.378 29.552 29.722

 [85] 29.888 30.050 30.208 30.362 30.512 30.658 30.800

 [92] 30.938 31.072 31.202 31.328 31.450 31.568 31.682

 [99] 31.792 31.898 32.000

> y2<- 2+5*x-0.4*x^2

> y2

  [1]  2.000  2.496  2.984  3.464  3.936  4.400  4.856

  [8]  5.304  5.744  6.176  6.600  7.016  7.424  7.824

 [15]  8.216  8.600  8.976  9.344  9.704 10.056 10.400

 [22] 10.736 11.064 11.384 11.696 12.000 12.296 12.584

 [29] 12.864 13.136 13.400 13.656 13.904 14.144 14.376

 [36] 14.600 14.816 15.024 15.224 15.416 15.600 15.776

 [43] 15.944 16.104 16.256 16.400 16.536 16.664 16.784

 [50] 16.896 17.000 17.096 17.184 17.264 17.336 17.400

 [57] 17.456 17.504 17.544 17.576 17.600 17.616 17.624

 [64] 17.624 17.616 17.600 17.576 17.544 17.504 17.456

 [71] 17.400 17.336 17.264 17.184 17.096 17.000 16.896

 [78] 16.784 16.664 16.536 16.400 16.256 16.104 15.944

 [85] 15.776 15.600 15.416 15.224 15.024 14.816 14.600

 [92] 14.376 14.144 13.904 13.656 13.400 13.136 12.864

 [99] 12.584 12.296 12.000

> y3 <- 2+4*x-0.6*xˆ2+0.04*xˆ3

Error: unexpected input in "y3 <- 2+4*x-0.6*xˆ"

> y3 <- 2+4*x-0.6*x^2+0.04*x^3

> plot(y3)

> plot(y2)

> y4 <- 2+4*x+2*x^2-0.6*x^3+0.04*x^4

> plot(y4)

> 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")

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

> y1 <- x/(2+5*x)

> y2 <- 1/(x-2+4/x)

> y3 <- 1/(xˆ2-2+4/x)

Error: unexpected input in "y3 <- 1/(xˆ"

> y3 <- 1/(x^2-2+4/x)

> plot(x,y1,type="l",ylab="y",main="Michaelis-Menten")

> plot(x,y2,type="l",ylab="y",main="shallow hump")

> plot(x,y3,type="t",ylab="y",main="steep hump")

Error in plot.xy(xy, type, ...) : invalid plot type 't'

> plot(x,y3,type="l",ylab="y",main="steep hump")

> #Gamma function

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

> t <- seq(0.2,4,0.01)

> t

  [1] 0.20 0.21 0.22 0.23 0.24 0.25 0.26 0.27 0.28 0.29 0.30

 [12] 0.31 0.32 0.33 0.34 0.35 0.36 0.37 0.38 0.39 0.40 0.41

 [23] 0.42 0.43 0.44 0.45 0.46 0.47 0.48 0.49 0.50 0.51 0.52

 [34] 0.53 0.54 0.55 0.56 0.57 0.58 0.59 0.60 0.61 0.62 0.63

 [45] 0.64 0.65 0.66 0.67 0.68 0.69 0.70 0.71 0.72 0.73 0.74

 [56] 0.75 0.76 0.77 0.78 0.79 0.80 0.81 0.82 0.83 0.84 0.85

 [67] 0.86 0.87 0.88 0.89 0.90 0.91 0.92 0.93 0.94 0.95 0.96

 [78] 0.97 0.98 0.99 1.00 1.01 1.02 1.03 1.04 1.05 1.06 1.07

 [89] 1.08 1.09 1.10 1.11 1.12 1.13 1.14 1.15 1.16 1.17 1.18

[100] 1.19 1.20 1.21 1.22 1.23 1.24 1.25 1.26 1.27 1.28 1.29

[111] 1.30 1.31 1.32 1.33 1.34 1.35 1.36 1.37 1.38 1.39 1.40

[122] 1.41 1.42 1.43 1.44 1.45 1.46 1.47 1.48 1.49 1.50 1.51

[133] 1.52 1.53 1.54 1.55 1.56 1.57 1.58 1.59 1.60 1.61 1.62

[144] 1.63 1.64 1.65 1.66 1.67 1.68 1.69 1.70 1.71 1.72 1.73

[155] 1.74 1.75 1.76 1.77 1.78 1.79 1.80 1.81 1.82 1.83 1.84

[166] 1.85 1.86 1.87 1.88 1.89 1.90 1.91 1.92 1.93 1.94 1.95

[177] 1.96 1.97 1.98 1.99 2.00 2.01 2.02 2.03 2.04 2.05 2.06

[188] 2.07 2.08 2.09 2.10 2.11 2.12 2.13 2.14 2.15 2.16 2.17

[199] 2.18 2.19 2.20 2.21 2.22 2.23 2.24 2.25 2.26 2.27 2.28

[210] 2.29 2.30 2.31 2.32 2.33 2.34 2.35 2.36 2.37 2.38 2.39

[221] 2.40 2.41 2.42 2.43 2.44 2.45 2.46 2.47 2.48 2.49 2.50

[232] 2.51 2.52 2.53 2.54 2.55 2.56 2.57 2.58 2.59 2.60 2.61

[243] 2.62 2.63 2.64 2.65 2.66 2.67 2.68 2.69 2.70 2.71 2.72

[254] 2.73 2.74 2.75 2.76 2.77 2.78 2.79 2.80 2.81 2.82 2.83

[265] 2.84 2.85 2.86 2.87 2.88 2.89 2.90 2.91 2.92 2.93 2.94

[276] 2.95 2.96 2.97 2.98 2.99 3.00 3.01 3.02 3.03 3.04 3.05

[287] 3.06 3.07 3.08 3.09 3.10 3.11 3.12 3.13 3.14 3.15 3.16

[298] 3.17 3.18 3.19 3.20 3.21 3.22 3.23 3.24 3.25 3.26 3.27

[309] 3.28 3.29 3.30 3.31 3.32 3.33 3.34 3.35 3.36 3.37 3.38

[320] 3.39 3.40 3.41 3.42 3.43 3.44 3.45 3.46 3.47 3.48 3.49

[331] 3.50 3.51 3.52 3.53 3.54 3.55 3.56 3.57 3.58 3.59 3.60

[342] 3.61 3.62 3.63 3.64 3.65 3.66 3.67 3.68 3.69 3.70 3.71

[353] 3.72 3.73 3.74 3.75 3.76 3.77 3.78 3.79 3.80 3.81 3.82

[364] 3.83 3.84 3.85 3.86 3.87 3.88 3.89 3.90 3.91 3.92 3.93

[375] 3.94 3.95 3.96 3.97 3.98 3.99 4.00

> plot(t,gamma(t),type="l")

> abline(h=1,lty=2)

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

> x <- seq(0,10,0.1)

> plot(x,y,type="l",main="three-parameter logistic")

Error in xy.coords(x, y, xlabel, ylabel, log) : object 'y' not found

> y <- 100/(1+90*exp(-1*x))

> plot(x,y,type="l",main="three-parameter logistic")

> y <- 20+100/(1+exp(0.8*(3-x)))

> plot(x,y,ylim=c(0,140),type="l",main="four-parameter logistic")

> x <- -200:100

> y <- 100*exp(-exp(0.02*x))

> plot(x,y,type="l",main="negative Gompertz")

> x <- 0:100

> y <- 50*exp(-5*exp(-0.08*x))

> plot(x,y,type="l",main="positive Gompertz")

> #Biexponential model

> a <- 10

> b <- -0.8

> c <- 10

> d <- -0.05

> y <- a*exp(b*x)+c*exp(d*x)

> plot(x,y,main="+ - + -",type="l")

> a <- 10

> b <- -0.8

> c <- 10

> d <- 0.05

> y <- a*exp(b*x)+c*exp(d*x)

> plot(x,y,main="+ - + +",type="l")

> a <- 200

> b <- 0.2

> c <- -1

> d <- 0.7

> y <- a*exp(b*x)+c*exp(d*x)

> plot(x,y,main="+ + - +",type="l")

> a <- 200

> b <- 0.05

> c <- 300

> d <- -0.5

> y <- a*exp(b*x)+c*exp(d*x)

> plot(x,y,main="+ + + -",type="l")

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

> x <- 0:6

> plot(x,factorial(x),type="s",main="factorial x",log="y")

> choose(8,3)

[1] 56

> plot(0:8,choose(8,0:8),type="s",main="binomial coefficients")

> curve(pnorm(x),-3,3)

> arrows(-1,0,-1,pnorm(-1),col="red")

> arrows(-1,pnorm(-1),-3,pnorm(-1),col="green")

> pnorm(-1)

[1] 0.1586553

> curve(dnorm(x),-3,3)

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

> x <- seq(-3,3,0.01)

> y <- exp(-abs(x))

> plot(x,y,type="l",main= "x")

> y <- exp(-abs(x)ˆ3)

Error: unexpected input in "y <- exp(-abs(x)ˆ"

> y <- exp(-abs(x)^3)

> plot(x,y,type="l",main= "x")

> y <- exp(-abs(x)^8)

> plot(x,y,type="l",main= "x")

> y <- exp(-abs(x)^21)

> plot(x,y,type="l",main= "x")

> x <- seq(-3,3,0.01)

> z <- seq(-3,-1.25,0.01)

> p <- dnorm(z)

> z <- c(z,-1.25,-3)

> p <- c(p,min(p),min(p))

> plot(x,dnorm(x),type="l",xaxt="n",ylab="probability density",xlab="height")

> axis(1,at=-3:3,labels=c("146","154","162","170","178","186","192"))

> polygon(z,p,col="red")

> z <- seq(1.875,3,0.01)

> p <- dnorm(z)

> z <- c(z,3,1.875)

> p <- c(p,min(p),min(p))

> plot(x,dnorm(x),type="l",xaxt="n",ylab="probability density",xlab="height")

> axis(1,at=-3:3,labels=c("146","154","162","170","178","186","192"))

> polygon(z,p,col="red")

> z <- seq(-0.635,1.25,0.01)

> p <- dnorm(z)

> z <- c(z,1.25,-0.635)

> p <- c(p,0,0)

> plot(x,dnorm(x),type="l",xaxt="n",ylab="probability density",xlab="height")

> axis(1,at=-3:3,labels=c("146","154","162","170","178","186","192"))

> polygon(z,p,col="red")

> #The central limit theorem

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

> hist(runif(10000)*10,main="")

> means <- numeric(10000)

> for (i in 1:10000){

+     means[i] <- mean(runif(5)*10)

+ }

> hist(means,ylim=c(0,1600),main="")

> mean(means)

[1] 4.990829

> sd(means)

[1] 1.293516

> xv <- seq(0,10,0.1)

> yv <- dnorm(xv,mean=4.998581,sd=1.28996)*5000

> lines(xv,yv)

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

> hist(sample(1:6,replace=T,10000),breaks=0.5:6.5,main="",xlab="one die")

> a <- sample(1:6,replace=T,10000)

> b <- sample(1:6,replace=T,10000)

> hist(a+b,breaks=1.5:12.5,main="", xlab="two dice")

> c <- sample(1:6,replace=T,10000)

> hist(a+b+c,breaks=2.5:18.5,main="", xlab="three dice")

> d <- sample(1:6,replace=T,10000)

> e <- sample(1:6,replace=T,10000)

> hist(a+b+c+d+e,breaks=4.5:30.5,main="", xlab="five dice")

> mean(a+b+c+d+e)

[1] 17.5012

> sd(a+b+c+d+e)

[1] 3.819877

> lines(seq(1,30,0.1),dnorm(seq(1,30,0.1),17.5937,3.837668)*10000)

> #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="")

Post a comment

0 Comments