第 3 章 数学

Code

library(tidyverse)
library(plotly)

3.1 距离最短

Code

(A=matrix(c(1,3,6,-1,2,5,2,5,9,6,5.5,11),byrow=F,ncol=2))

##      [,1] [,2]
## [1,]    1  2.0
## [2,]    3  5.0
## [3,]    6  9.0
## [4,]   -1  6.0
## [5,]    2  5.5
## [6,]    5 11.0

Code

b=A[,1]%*%matrix(1,ncol=6,nrow=1)
c=A[,2]%*%matrix(1,ncol=6,nrow=1)
(d=sqrt(((b-t(b))^2+(c-t(c))^2)))

##       [,1]  [,2]  [,3]  [,4]  [,5]  [,6]
## [1,] 0.000 3.606 8.602 4.472 3.640 9.849
## [2,] 3.606 0.000 5.000 4.123 1.118 6.325
## [3,] 8.602 5.000 0.000 7.616 5.315 2.236
## [4,] 4.472 4.123 7.616 0.000 3.041 7.810
## [5,] 3.640 1.118 5.315 3.041 0.000 6.265
## [6,] 9.849 6.325 2.236 7.810 6.265 0.000

Code

dist(A)

##       1     2     3     4     5
## 2 3.606                        
## 3 8.602 5.000                  
## 4 4.472 4.123 7.616            
## 5 3.640 1.118 5.315 3.041      
## 6 9.849 6.325 2.236 7.810 6.265

Code

n=nrow(A)
ax=A[,1]%*%matrix(1,ncol=6,nrow=1)
ay=A[,2]%*%matrix(1,ncol=6,nrow=1)
dx=ax-t(ax)
dy=ay-t(ay)
b=sqrt(dx^2+dy^2)

3.2 牛顿迭代

Code

library(Deriv)  # Version 3.8.5 
f=function(x){
  log(x)/(1+x)
}
f1<- Deriv(f)
f1

## function (x) 
## {
##     .e1 <- 1 + x
##     (1/x - log(x)/.e1)/.e1
## }

Code

f2=Deriv(f1)
x=3
k=1
e=0.0001
for( i in 1:100){
x=x-f1(x)/f2(x)
if(abs(f1(x))<e){break}
k=k+1
}

3.3 微分

3.3.1 高次微分

Code

library(deSolve)
yvini <-c(y=0,v=1,z=8)
derivs3 <- function (t,y,parms){
  with(as.list(c(y,parms)), {
    dy<-v#1 order differentiation
    dv<-z#2 order differentiation
    dz <-(-12*exp(-t)-cc*y-b*v)/a# original equation
    list(c(dy,dv,dz))
  })
}
times <- seq(from = 0, to = 5, by = 0.01)
out3 <- ode(y = c(y=0,v=1,z=8), times = times, func = derivs3,
            parms = c(a=1,b=-3,cc=-2))# solve

f33<-function(t)(2/3*exp(2*t)-2/3*exp(-t)-t*exp(-t)+2*t^2*exp(-t))# Theoretical equation by hand

plot(out3[,1],out3[,2],type='l')# numerical solution
curve(f33,0,5,add=T,col="red",lty=2)# fit perfectly

3.3.2 微分方程组

Code

#(a)
library(deSolve)
model <- function(t, y, p) {
  with(as.list(c(y, p)), {
    dN <-  a*N-g*N^2-b*N*G
    dG <-  -c*G+d*N*G
    list(c(dN, dG))
  })
}

# (b)
y      <- c(N = 10000, G = 60)
parms  <- c(a = 5,b=0.01,c=100,d=0.01,g=0.0001)
times  <- seq(0, 10,0.001)
out <- ode(y, times, model, parms)
plot(out)

Code

out <- ode(y,times, model, parms)


# (c)
a = 5
b=0.01
c=100
d=0.01
g=0.0001
# is the solution when the right-hand side of the original equation is equal to 0
c/d#x

## [1] 10000

Code

(a*d*c-g*c^2)/(b*d*c)#y

## [1] 400

Code

# (d)
for(i in seq(0.2,0.4,0.001)){
  parms  <- c(a = 5,b=i,c=100,d=0.01,g=0.0001)
  out <- ode(y,times, model, parms)
  if(any(out[,3]<1))break
}
i # min b is 0.296

## [1] 0.296

3.4 option

3.4.1 1

Code

knitr::include_graphics('www/R_exam.png')

Code

# 6
Price=readxl::read_xlsx('data/data.xlsx',skip = 11)
Price=as.data.frame(Price)
colnames(Price)[4]="A"
head(Price,6)

##   Year Period    Label     A
## 1 2002    M01 2002 Jan 173.2
## 2 2002    M02 2002 Feb 173.7
## 3 2002    M03 2002 Mar 174.7
## 4 2002    M04 2002 Apr 175.8
## 5 2002    M05 2002 May 175.8
## 6 2002    M06 2002 Jun 175.9

Code

# 7
Price=cbind(T=c(1:240),Price)
plot(Price$T,Price$A,type="l")
lines(Price$T,Price$T*0.355+176.8)

Code

# 8
dat=Price$A#ts data
AM=c()
for(i in 7:length(dat)){
  AM[i-6]=mean(dat[i-6:i])
}
head(AM) #The prediction begins with the 7 phase

## [1] 173.2 173.4 173.9 174.3 174.6 174.8

Code

# 9
xl1=dat[-239:-240]#Yn-2
xl2=dat[c(-1,-240)]#Yn-1
xl=dat[-1:-2]#Y The prediction begins with the third phase
myfun<-function(myvalues){x1<-myvalues[1];x2<-myvalues[2];x3<-myvalues[3]; 
y=x1*xl1+x2*xl2+x3
e=abs(xl-y);sum(e)}
AR<-optim(c(1,1,1), myfun)
AR$par#Model coefficient

## [1] -0.5628  1.5583  1.2166

3.4.2 极值黑塞矩阵

Code

f5=function(x,y){
  l=(x^2+y-11)^2+(x+y^2-7)^2
  l
}
x=y=seq(-4.5,4.5,0.1)

z<- outer(x,y,f5)
persp(x,y,z, theta = 30, phi = 30,expand = 0.5,
      col = "lightblue", ltheta = 120,shade = 0.75,
      ticktype = "detailed",xlab ="x1",
      ylab = "x2",zlab ="f(x1,x2)")

Code

plotly::plot_ly(x = ~x, y = ~y, z = ~z,color = ~z, colors = 'Reds',opacity = 0.8)%>% add_surface()

Code

myfun<-function(myvalues){
  x<-myvalues[1]
  y<-myvalues[2]
  y=(x^2+y-11)^2+(x+y^2-7)^2
  y
}

AR<-optim(c(0,0), myfun,hessian = T)
AR$par#Model coefficient

## [1] 3 2

Code

eigen(AR$hessian)$values

## [1] 82.29 25.72

Code

myfun(AR$par)

## [1] 2.984e-07