R語(yǔ)言編程數(shù)學(xué)分析重讀微積分微分學(xué)原理運(yùn)用
1 連續(xù)性
比如下面這個(gè)隨機(jī)函數(shù)
x = seq(0,0.1,0.01) y = runif(11,0,1) plot(x,y) lines(x,y)
x = seq(0,0.01,0.00001) y = runif(1001,0,1) plot(x,y)
無(wú)論我們把區(qū)間縮小到什么程度，這種亂糟糟的仿佛要鋪滿整個(gè)坐標(biāo)圖的點(diǎn)的樣式并不會(huì)發(fā)生變化。
也就是說(shuō)， f(x)從左邊趨近于0的時(shí)候，f(x)在0處是連續(xù)的，而在右側(cè)趨近于0的時(shí)候，卻并不連續(xù)。此即左連續(xù)和右連續(xù)的概念。
2 求導(dǎo)
回到切線的問題，如果曲線 y=f(x)在 x0點(diǎn)并不連續(xù)，那么這點(diǎn)顯然沒有唯一的一條切線。
有的時(shí)候，盡管滿足了連續(xù)性的要求，也不一定存在切線，比如
y = ∣ x ∣
x = seq(-10,10) y = abs(x) plot(x,y,type='l')
在 x=0的位置，我們找到的切線要么和左邊重合，要么和右邊重合，也就是說(shuō)這個(gè)函數(shù)在x=0處存在兩條切線。
同時(shí)也就意味著這一點(diǎn)有兩個(gè)斜率，兩個(gè)導(dǎo)數(shù)。所以，如果把導(dǎo)數(shù)定義為某種映射，則一個(gè)點(diǎn)只能映射為一個(gè)值，所以只能定義這點(diǎn)的導(dǎo)數(shù)不存在。
3 數(shù)值導(dǎo)數(shù)
根據(jù)導(dǎo)數(shù)的定義，當(dāng)函數(shù)的定義域不連續(xù)時(shí)，其不連續(xù)處顯然是不存在導(dǎo)數(shù)的，但圖形可以“欺騙”我們的眼睛。
> x = seq(-1,1,0.1)
> y = sin(x)
> y1 = cos(x)
> xEnd = x+0.1
> yEnd = y+y1*0.1
> plot(x,y)
> for(i in seq_along(x)){
+ lines(c(x[i],xEnd[i]),c(y[i],yEnd[i]),col="red")
+ }
上圖中，圓圈是對(duì)y = sinx 進(jìn)行抽樣的結(jié)果，可以理解為不連續(xù)函數(shù)；紅線則是在每一個(gè)分立的點(diǎn)上，sinx在該點(diǎn)的切線。這兩者看上去如此一致，說(shuō)明連續(xù)函數(shù)的導(dǎo)數(shù)在抽樣之后仍然具備一定的數(shù)學(xué)意義。相應(yīng)地，不連續(xù)的函數(shù)，也應(yīng)該有類似于導(dǎo)數(shù)一樣的存在，從而與連續(xù)函數(shù)的導(dǎo)數(shù)相對(duì)應(yīng)，此即數(shù)值導(dǎo)數(shù)。如果回顧導(dǎo)數(shù)的定義
x = seq(-5,5,0.1) y = cos(x) x1 = seq(-5,5,0.1) end = length(x1) y1 = (sin(x1[2:end])-sin(x1[1:end-1]))/0.1 x5 = seq(-5,5,0.5) end = length(x5) y5 = (sin(x5[2:end])-sin(x5[1:end-1]))/0.5 x10 = seq(-5,5,1) end = length(x10) y10 = (sin(x10[2:end])-sin(x10[1:end-1]))/0.5 plot(x,y,type="l",col="red") lines(x1[2:length(x1)],y1) lines(x5[2:length(x5)],y5) lines(x10[2:length(x10)],y10)
如圖所示
由于我們采用的是后向差分，所以三組差商的值整體右移。此外，隨著 h的增大，其誤差也越來(lái)越明顯。
對(duì)一個(gè)函數(shù)進(jìn)行反復(fù)求導(dǎo)，即可得到高階導(dǎo)數(shù)，可以用數(shù)學(xué)歸納法的方式記為
4 差商與牛頓插值
根據(jù)這個(gè)表達(dá)式，可以通過(guò)一個(gè)簡(jiǎn)單的遞歸程序計(jì)算數(shù)組的差商
# 差商算法,x,y為同等長(zhǎng)度的數(shù)組
diffDiv<-function(x,y){
end = length(x)
ind = 2:end #索引
return(
if(end==1) y[1]
else (diffDiv(x[ind],y[ind])
-diffDiv(x[ind-1],y[ind-1]))/(x[end]-x[1])
)
}
如果要計(jì)算階數(shù)為k的差商，只需重復(fù)調(diào)用diffDiv，
kDiffDiv <- function(x,y,k){
len <- length(x)
if(len<k) return(0)
d<-rep(0,len-k)
for(i in 1:(len-k))
d[i] <- diffDiv(x[i:(i+k)],y[i:(i+k)])
return(d)
}
據(jù)此，繪制出 y=x^5這個(gè)函數(shù)的差商，
> plot(x,y) > k = 1 > lines(x[1:(end-k)],kDiffDiv(x,y,k),col="red") > k = 2 > lines(x[1:(end-k)],kDiffDiv(x,y,k),col="green") > k = 3 > lines(x[1:(end-k)],kDiffDiv(x,y,k),col="blue") > k = 4 > lines(x[1:(end-k)],kDiffDiv(x,y,k),col="purple") > k = 5 > lines(x[1:(end-k)],kDiffDiv(x,y,k),col="yellow") > k = 6 > lines(x[1:(end-k)],kDiffDiv(x,y,k),col="gray")
如圖所示
可見差商與微分在階數(shù)上有著一致的趨勢(shì)。那么，知道了差商之后，可以做點(diǎn)什么呢？
根據(jù)差商定義，可得
polyMulti<-function(x){
n = length(x)
if(n<2) return(c(-x,1))
omega = rep(0,n+1)
omega[1] = - x[1]
omega[2] = 1
for(i in 2:n){
omega[2:i] = -x[i]*omega[2:i]*+omega[1:(i-1)]
omega[1] = -x[i]*omega[1]
omega[i+1] = 1
}
return(omega)
}
Newton<-function(x,y){
n = length(x)
N = rep(0,n+1)
N[1]=y[1]
for(i in 1:n){
omega = polyMulti(x[1:i])
d = kDiffDiv(x[1:i],y[1:i],i-1)
N[1:i] = N[1:i] + d*omega[1:i]
N[i+1] = d*omega[i+1]
}
return(N)
}
驗(yàn)證一下
x = sort(rnorm(10))
y = x^5+2*x^4
N = Newton(x,y)
Y = y*0
for(i in 1:length(Y))
for(j in 1:length(N))
Y[i] = Y[i]+N[j]*x[i]^(j-1)
plot(x,y)
lines(x,Y)
可見效果還是不錯(cuò)的。
5 方向?qū)?shù)
現(xiàn)繪制出100個(gè)隨機(jī)點(diǎn)處x方向的偏導(dǎo)數(shù)
x = matrix(data=seq(-5,4.95,0.05),nrow=200,ncol=200)
y = t(x)
z = 1-(x^2+y^2)
library(rgl)
persp3d(x,y,z,col="red")
N = 1000
x0 = rnorm(N)
y0 = rnorm(N)
z0 = 1-(x0^2+y0^2)
x1 = x0+3
z1 = -2*x0*3+z0
for(i in 1:N)
lines3d(c(x0[i],x1[i]),c(y0[i],y0[i]),c(z0[i],z1[i]),col="green")
從觀感上來(lái)看，綠線的確是沿著 x x x方向。但是這個(gè)切線顯然不是唯一的， y y y軸方向顯然存在另一條切線。推而廣之，一旦坐標(biāo)系旋轉(zhuǎn)，那么曲面上任意一點(diǎn)處的 x和 y方向均會(huì)發(fā)生變化。
那么曲面是否存在一個(gè)只和曲面特征有關(guān)，而與坐標(biāo)系無(wú)關(guān)的參數(shù)？
在解決這個(gè)問題之前，最好先找到曲面某點(diǎn)沿著任意方向的導(dǎo)數(shù)?；仡?x方向偏導(dǎo)數(shù)的定義
如果導(dǎo)數(shù)的方向發(fā)生旋轉(zhuǎn)，則可以寫為
如果寫成矢量形式，則定義梯度
沿這些點(diǎn)梯度方向做一些直線，看看效果如何
x = matrix(data=seq(-5,4.95,0.05),nrow=200,ncol=200)
y = t(x)
z = -(x^2+y^2)
library(rgl)
persp3d(x,y,z,col="red")
theta = seq(-pi,pi,0.01)
x0 = cos(theta)
y0 = sin(theta)
z0 = 1-(x0^2+y0^2)
x1 = x0*0
y1 = y0*0
z1 = z0+2
for(i in 1:N)
lines3d(c(x0[i],x1[i]),c(y0[i],y1[i]),c(z0[i],z1[i]),col="green")
如圖所示，像一頂漂亮的帽子，在某個(gè)投影方向看去，和我們熟知的切線如出一轍。
6 全微分
做圖如下
theta = seq(-pi,pi,0.1)
x = cos(theta)
y = sin(theta)
plot(x,y,type='l',col='red')
x1 = x*0
y1 = (x1-x)/(-2*x)*(-2*y)+y
for(i in 1:length(theta))
lines(c(x[i],x1[i]),c(y[i],y1[i]),col='green')
可見，梯度方向?qū)?yīng)的是圖形的法線方向。對(duì)于二維平面內(nèi)的曲線而言，其法線方向與切線方向垂直。
回顧偏導(dǎo)數(shù)的定義
7 法線
相應(yīng)地最大方向?qū)?shù)的方向即為梯度的歸一化
現(xiàn)隨機(jī)選擇一些點(diǎn)，來(lái)繪制一下這四個(gè)方向的向量
library(rgl)
N = 1500
x<-rnorm(N)
y<-rnorm(N)
z<-1-x^2-y^2
for(i in 1:N){
lines3d(c(x[i],3*x[i]),c(y[i],3*y[i]),c(z[i],z[i]+1),col='red')
if(y[i]>0.1)
lines3d(c(x[i],x[i]),c(y[i],y[i]-1/y[i]/2),c(z[i],z[i]+1),col='green')
if(x[i]>0.1)
lines3d(c(x[i],x[i]-1/x[i]/2),c(y[i],y[i]),c(z[i],z[i]+1),col='green')
lines3d(c(x[i],x[i]*(1-2*z[i])/(2-2*z[i])),c(y[i],y[i]*(1-2*z[i])/(2-2*z[i])),c(z[i],z[i]+1),col='green')
}
可以看到，綠線幾乎重新編織了一遍原函數(shù)，而紅線則刺破了曲面。
8 偏導(dǎo)數(shù)和邊緣檢測(cè)
基于偏導(dǎo)數(shù)的邊緣檢測(cè)
灰度圖像是天然的z=f(x,y)函數(shù)，盡管以一種差分化的形式存在。其中x,y分別代表圖像坐標(biāo)系中的坐標(biāo) z可以表示灰度圖像的灰度值。
那么接下來(lái)我們可以觀察一下偏導(dǎo)數(shù)作用在圖像上是一個(gè)什么效果。圖片當(dāng)然是最經(jīng)典的lena
library(imager)
img = load.image("lena.jpg")
dim(img)
[1] 512 512 1 3
gray = grayscale(img)
par(mfrow=c(1,2))
plot(img)
plot(gray)
可見gray顯然為灰度圖像，從其維度就能看得出來(lái)，然后將其變?yōu)槎S的數(shù)組，接下來(lái)就可以進(jìn)行求導(dǎo)操作了。
dim(gray) [1] 512 512 1 1 mat = array(gray,dim=c(512,512)) mat_x = diff(mat,1) mat_y = t(diff(t(mat),1)) par(mfrow=c(1,2)) image(mat_x) image(mat_y)
由于圖像坐標(biāo)系默認(rèn)是從上向下為 y y y軸，從左向右為 x x x軸，所以在我們熟知的坐標(biāo)系中，圖像是上下顛倒的。而且R語(yǔ)言還非常智能（障）地添加了一層偽彩色，這讓我們更加清晰地看出，對(duì)圖像進(jìn)行差分操作，提取出了邊緣信息。
這很容易理解，所謂“邊緣”，往往意味著變化較大的點(diǎn)，如果我們抽取lena圖的任意一行，
a=mat[256,] par(mfrow=c(1,2)) plot(a,type='l') plot(diff(a,1),type='l')
在變化劇烈處，相應(yīng)地導(dǎo)數(shù)較大。
若將偏導(dǎo)數(shù)在圖像空間中展開，由于任意兩個(gè)像素點(diǎn)之間的差恒為1，則可得到其差分形式
Roberts算子
求導(dǎo)是對(duì)整個(gè)函數(shù)的定義域展開的一次性操作，但在考察其差分形式之后卻發(fā)現(xiàn)，數(shù)值偏導(dǎo)數(shù)可以寫成一種對(duì)局部區(qū)域的反復(fù)操作。
例如，就前向差分而言，可以對(duì)圖像的任意一個(gè)子矩陣
theta = c(pi/3,pi/4,pi/5,pi/6) par(mfrow=c(2,2)) for(i in 1:4) image(mat_x[,1:511]*cos(theta[i])+mat_y[1:511,]*sin(theta[i]))
可以看到，最后一張圖片的法向角度為 30° ，而其右下角正好有一個(gè) 30°附近的清晰的邊緣。
其他算子
如果通過(guò)中心差分來(lái)定義算子，則統(tǒng)一維度后，其x和 y向的梯度算子分別寫為
以上就是R語(yǔ)言編程數(shù)學(xué)分析重讀微積分微分學(xué)原理運(yùn)用的詳細(xì)內(nèi)容，更多關(guān)于R語(yǔ)言數(shù)學(xué)分析微分學(xué)原理的資料請(qǐng)關(guān)注腳本之家其它相關(guān)文章！
相關(guān)文章
R語(yǔ)言繪制corrplot相關(guān)熱圖分析美化示例及詳細(xì)圖解
這篇文章主要為大家介紹了R語(yǔ)言corrplot相關(guān)熱圖分析美化示例及詳細(xì)圖解，有需要的朋友可以借鑒參考下，希望能夠有所幫助，祝大家多多進(jìn)步，早日升職加薪2022-06-06
這篇文章主要介紹了基于R語(yǔ)言for循環(huán)的替換方案，具有很好的參考價(jià)值，希望對(duì)大家有所幫助。一起跟隨小編過(guò)來(lái)看看吧2021-04-04
詳解R語(yǔ)言實(shí)現(xiàn)前向逐步回歸(前向選擇模型)
本文主要介紹了詳解R語(yǔ)言實(shí)現(xiàn)前向逐步回歸，從實(shí)現(xiàn)原理開始，文中通過(guò)示例代碼介紹的非常詳細(xì)，具有一定的參考價(jià)值，感興趣的小伙伴們可以參考一下2021-08-08
R語(yǔ)言常量知識(shí)點(diǎn)總結(jié)
在本篇文章里小編給大家整理了一篇關(guān)于R語(yǔ)言常量知識(shí)點(diǎn)總結(jié)內(nèi)容，有興趣的朋友們可以學(xué)習(xí)分享下。2021-03-03
R包制作后出現(xiàn)not available for錯(cuò)誤問題解決解決
這篇文章主要為大家介紹了R包制作后出現(xiàn)not available for...錯(cuò)誤的問題解決方式，有需要的朋友，可以借鑒參考下，希望能夠有所幫助2021-11-11
這篇文章主要介紹了R語(yǔ)言gsub替換字符工具的具體使用，文中通過(guò)示例代碼介紹的非常詳細(xì)，對(duì)大家的學(xué)習(xí)或者工作具有一定的參考學(xué)習(xí)價(jià)值，需要的朋友們下面隨著小編來(lái)一起學(xué)習(xí)學(xué)習(xí)吧2021-03-03
R語(yǔ)言—自定義函數(shù)求置信區(qū)間的操作
這篇文章主要介紹了R語(yǔ)言—自定義函數(shù)求置信區(qū)間的操作，具有很好的參考價(jià)值，希望對(duì)大家有所幫助。一起跟隨小編過(guò)來(lái)看看吧2021-04-04
R語(yǔ)言學(xué)習(xí)VennDiagram包繪制韋恩圖示例
這篇文章主要為大家介紹了R語(yǔ)言學(xué)習(xí)VennDiagram包繪制韋恩圖示例，有需要的朋友可以借鑒參考下，希望能夠有所幫助，祝大家多多進(jìn)步，早日升職加薪2022-06-06
日期時(shí)間是常用的一種類型，本文主要介紹了R語(yǔ)言日期時(shí)間的使用，文中通過(guò)示例代碼介紹的非常詳細(xì)，具有一定的參考價(jià)值，感興趣的小伙伴們可以參考一下2022-03-03
這篇文章主要介紹了Rstudio 修改工作路徑(三種方法總結(jié))，具有很好的參考價(jià)值，希望對(duì)大家有所幫助。一起跟隨小編過(guò)來(lái)看看吧2021-04-04