[1] 3
四分位数(Quartile): 用于描述任何类型的数据,尤其是偏态数据的离散程度,即将全部数据从小到大排列,正好排列在上1/4位置叫上四分位数,下1/4位置上的数就叫做下四分位数.
R程序:计算样本(1,2,3,4,5,6,7,8,9)的四分位数
> S<-c(1,2,3,4,5,6,7,8,9) > quantile(S)
0% 25% 50% 75% 100% 1 3 5 7 9 > fivenum(S) [1] 1 3 5 7 9
通用的计算统计函数:
R程序:计算样本(1,2,3,4,5,6,7,8,9)的统计函数
> S<-c(1,2,3,4,5,6,7,8,9) > summary(S)
Min. 1st Qu. Median Mean 3rd Qu. Max. 1 3 5 5 7 9
6). 协方差(Covariance)
协方差用于衡量两个变量的总体误差。而方差是协方差的一种特殊情况,即当两个变量是相同的情况。设X,Y为两个随机变量,称E{[X-E(X)][Y-E(Y)]}为X和Y的协方差,记录Cov(X,Y)。
R程序:计算X(1,2,3,4)和Y(5,6,7,8)的协方差
> X<-c(1,2,3,4) > Y<-c(5,6,7,8) > cov(X,Y) [1] 1.666667
7). 相关系数(Correlation coefficient)
相关系数是用以反映变量之间相关关系密切程度的统计指标。相关系数是按积差方法计算,同样以两变量与各自平均值的离差为基础,通过两个离差相乘来反映两变量之间相关程度。当Var(X)>0, Var(Y)>0时,称Cov(X,Y)/sqrt(Var(X)*Var(Y))为X与Y的相关系统。
R程序:计算X(1,2,3,4)和Y(5,7,8,9)的相关系数
> X<-c(1,2,3,4) > Y<-c(5,7,8,9) > cor(X,Y) [1] 0.9827076
8). 矩
原点矩(moment about origin): 对于正整数k,如果E|X^k|存在,称V^k=E(X^k)为随机变量X的k阶原点矩。X的数学期望是X的一阶原点矩,即E(x)=v1.
R程序:计算S(1,2,3,4,5)的一阶原点矩(均值)
> S<-c(1,2,3,4,5) > mean(S) [1] 3
中心矩(moment about centre): 对于正整数k,如果EX存在,且E(|X - EX|^k)也存在,则称E[X-EX]^k为随机变量X的k阶中心矩。如X的方差是X的二阶中心矩,即D(X)=E{[X-E(X)]^2}
R程序:计算S(1,2,3,4,5)的二阶中心矩(方差)
> S<-c(1,2,3,4,5) > var(S) [1] 2.5
距是广泛应用的一类数学特征,均值和方差分别就是一阶原点矩和二阶中心矩。
偏度(skewness): 是统计数据分布偏斜方向和程度的度量,是统计数据分布非对称程度的数字特征。设分布函数F(x)有中心矩u2=E(X ?E(X))^2, u3 = E(X ?E(X))^3,则Cs=u3/u2^(3/2)为偏度系数。
当Cs>0时,概率分布偏向均值右则,Cs<0时,概率分布偏向均值左则。 R语言:计算10000个正态分布的样本的偏度
> library(PerformanceAnalytics) > S<-rnorm(10000) > skewness(S) [1] -0.00178084 > hist(S,breaks=100)
峰度(kurtosis): 又称峰态系数。表征概率密度分布曲线在平均值处峰值高低的特征数。峰度刻划不同类型的分布的集中和分散程序。设分布函数????F(x)有中心矩u2=E(X ?E(X))^2, u4=E(X ?E(X))^4,则Ck=u4/(u2^2-3)为峰度系数。
R语言:计算10000个正态分布的样本的峰度,(同偏度的样本数据)
> library(PerformanceAnalytics) > kurtosis(S) [1] -0.02443549 > hist(S,breaks=100)
8). 协方差矩阵(covariance matrix)
协方差矩阵是一个矩阵,其每个元素是各个向量元素之间的协方差。是从标量随机变量到高维度随机向量的自然推广。设X = (X1,X2, ... ,Xn), Y = (Y1, Y2, ..., Ym) 为两个随机变量,则Cov(X,Y)为X,Y的协方差矩阵.
R语言:计算协方差矩阵
> x=as.data.frame(matrix(rnorm(10),ncol=2)) > x
V1 V2 1 -2.11315384 -2.55189840 2 -0.96631271 -1.36148355 3 -0.02835058 -0.82328774 4 -1.86669567 -0.07201353 5 0.27324957 -2.23835218
> var(x)
V1 V2 V1 1.13470650 -0.09292042 V2 -0.09292042 1.03172261
> cov(x)
V1 V2 V1 1.13470650 -0.09292042 V2 -0.09292042 1.03172261
3. 极限定理
? ?
大数定律 中心极限定理
1). 大数定律
大数定律(law of large numbers),又称大数定理,是判断随机变量的算术平均值是否向常数收敛的定律,是概率论和数理统计学的基本定律之一。 设X1,X2,...,Xk, 是随机变量序列且E(Xk)存在(k=1,2,3...), Yn = 1/n * (X1 +X2+ ... + Xk),对于任意给定的ε > 0, 有
则称随机变量序列{Xk}服从大数定律。 三个重要定律
? ? ?
Bernoulli大数定律
Chebyshev(切比雪夫)大数定律 Khintchin(辛钦)大数定律
Bernoulli(贝努力)大数定律
设Na是n次独立重复试验中A发生的次数,p是事件A在每次试验中发生的概率,则对任意的正数ε > 0,有
Bernoulli大数定律揭示了“频率稳定于概率”说法的实质。 Chebyshev(切比雪夫)大数定律
设随机变量X1,X2,...Xk相互独立,且具有相同的期望与方差:E(Xk)=μ, Var(Xk) = σ^2, (k = 1, 2, ...), 则对于任意的正数ε > 0, ??有
Khintchin(辛钦)大数定律
设随机变量X1,X2...Xk相互独立,服从相同的分布,且其期望E(Xk) = μ , (k = 1, 2,...), 则对于任意的正数ε > 0, ??有
若对随机变量序列X1, X2, ...Xk存在常数a, 使得对于任意的正数ε > 0, ??有
成立,则称Xk依概率收敛于a,则Chebyshev大数定律和Khintchin大数定律有 大数定律定理
设随机变量X具有期望E(X)=μ,方差Var(X) = σ2, 则对于任意ε > 0, ??有
R语言:假设投硬币,正面概率是0.5,投4次时,计算得到2次正面的概率?根据大数定律,如果投是10000次,计算5000次正面的概率?
#计算2次正面的的概率
> choose(4,2)/2^4 #choose组合数的计算:从4中选择2个 [1] 0.375
#计算5000次正面的的概率
> pbinom(5000, 10000, 0.5) #pbinom二向分布,5000为分位数,产生10000个随机数,每个概率0.5 [1] 0.5039893
2). 中心极限定理(central limit theorem)
中心极限定理是判断随机变量序列部分和的分布是否渐近于正态分布的一类定理。在自然界及生产科学实践中,一些现象受到许多相互独立的随机因素的影响,如果每个因素的影响都很小,那么部的影响可以看作是服从正太分布。中心极限定理正是从数学上论证了这一现象。
设从均值为μ、方差为σ^2;(有限)的任意一个总体中抽取样本量为n的样本,当n充分大时,样本均值的抽样分布近似服从均值为μ、方差为σ^2/n的正态分布。 两个最著名的中心极限宣
? ?
列维定理(Lindburg-Levy)
拉普拉斯定理(de Movire - Laplace)
列维定理(Lindburg-Levy)
即独立同分布随机变量序列的中心极限定理。它表明,独立同分布、且数学期望和方差有限的随机变量序列的标准化和以标准正态分布为极限。
设随机变量X1,X2,......Xn,......相互独立,服从同一分布,且具有数学期望和方差:E(Xk)=μ,D(Xk)=σ^2>0(k=1,2....),则随机变量之和的标准化变量的分布函数Fn(x)对于任意x满足limFn(x)=Φ(x),n→∞ 其中Φ(x)是标准正态分布的分布函数。 拉普拉斯定理(de Movire - Laplace)
即服从二项分布的随机变量序列的中心极限定理。它指出,参数为n, p的二项分布以np为均值、np(1-p)为方差的正态分布为极限。
R的极客理想系列文章,涵盖了R的思想,使用,工具,创新等的一系列要点,以我个人的学习和体验去诠释R的强大。
R语言作为统计学一门语言,一直在小众领域闪耀着光芒。直到大数据的爆发,R语言变成了一门炙手可热的数据分析的利器。随着越来越多的工程背景的人的加入,R语言的社区在迅速扩大成长。现在已不仅仅是统计领域,教育,银行,电商,互联网….都在使用R语言。
要成为有理想的极客,我们不能停留在语法上,要掌握牢固的数学,概率,统计知识,同时还要有创新精神,把R语言发挥到各个领域。让我们一起动起来吧,开始R的极客理想。
前言
覆盖R基础知识,快速上手,RHadoop环境的搭建基础课。 目录 1. 2. 3. 4. 5.
背景知识 开发环境 R语法 R基本函数 R的扩展包
1. 背景知识
R起源
R 是一个有着统计分析功能及强大作图功能的软件系统,是由奥克兰大学统计学系的Ross Ihaka和Robert Gentleman 共同创立。由于R 受Becker, Chambers & Wilks 创立的S 和Sussman 的Scheme两种语言的影响,所以R 看起来和S 语言非常相似。
R 是一个世界范围统计工作者共同协作的产物,至2013 年2 月共计近5000 个包可在互联网上自由下载,这些都是各行业数据分析同行的工作结晶。 R的特点 1. 2. 3. 4. 5. 6. 7. 8. 9.
有效的数据处理和保存机制。
拥有一整套数组和矩阵的操作运算符。 一系列连贯而又完整的数据分析中间工具。
图形统计可以对数据直接进行分析和显示,可用于多种图形设备。
一种相当完善、简洁和高效的程序设计语言。它包括条件语句、循环语句、用户自定义的递归函数以及输入输出接口。 R语言是彻底面向对象的统计编程语言。
R语言和其它编程语言、数据库之间有很好的接口。
R语言是自由软件,可以放心大胆地使用,但其功能却不比任何其它同类软件差。 R语言具有丰富的网上资源
R的下载和安装
R是一个免费的自由软件,它有UNIX、LINUX、MacOS和WINDOWS版本,都是可以免费下载和使用的,在那儿可以下载到R的安装程序、各种外挂程序和文档。在R的安装程序中只包含了8个基础模块,其他外在模块可以通过CRAN获得。 R的官方网站: http://www.r-project.org/ Linux Ubuntu的R安装
~ sudo vi /etc/apt/sources.list
deb http://mirror.bjtu.edu.cn/cran/bin/linux/ubuntu precise/
~ sudo apt-get update
~ sudo apt-get install r-base-core=2.15.3-1precise0precise1
2. 开发环境
R命令行环境:
R默认的界面环境:
RStudio的IDE开发环境: http://www.rstudio.com/
3. R语法
R是一种语法非常简单的表达式语言(expression language),大小写敏感。
可以在R 环境下使用的命名字符集依赖于R 所运行的系统和国家(系统的locale 设置),允许数字,字母,“.”,“_” 1). 命名
命名必须以”.”或者字母开头,以”.”开头时第二个字符不允许是数字。
2). 基本命令
基本命令要么是表达式(expressions),要么就是赋值(assignments)。
? ?
表达式,命令将被解析,并将结果显示在屏幕上,同时清空该命令所占内存。 赋值,命令将被解,并把值传给变量,但结果不会自动显示在屏幕上。
命令可以被”;”隔开或者另起一行。基本命令可以通过大括弧{},放在一起构成一个复合表达式。 注释:一行中以井号”#”开头
换行:如果一条命令在一行结束的时候在语法上还不完整,换行提示符,默认是+ 3). 基本的对象
R创建和控制的实体被称为对象。它们可以是变量,数组,字符串,函数,或者其他通过这些实体定义的一般性的结构。
? ? ? ? ?
矩阵(matrix)或者更为一般的数组(array)是多维的广义向量。实际上,它们就是向量,而且可以同时被两个或者更多个索引引用,并且以特有的方式显示出来。 因子(factor)为处理分类数据提供的一种有效方法。
列表(list)是一种泛化(general form)的向量。它没有要求所有元素是同一类型,许多时候它本身就是向量和列表类型。列表为统计计算的结果返回提供了一种便利的方法。
数据框(data frame)是和矩阵类似的一种结构。在数据框中,列可以是不同的对象。可以把数据框看作是一个行表示观测个体并且(可能)同时拥有数值变量和分类变量的`数据矩阵’ 。许多实验数据都可以很好的用数据框描述:处理方式是分类变量而响应值是数值变量。
函数(function)是可以保存在项目工作空间的R 对象。该对象为R 提供了一个简单而又便利的功能扩充方法。见编写你自己的函数
在R会话过程中,对象是通过名字创建和保存的。objects(), ls()可以显示当前会话的对象名字。rm()可以删除对象。 对象持久化
R 会话中创建的所有对象可以永久地保存在一个文件中以便于以后的R 会话调用。在每一次R 会话结束的时候,你可以保存当前所有可用的对象。如果你想这样做,这些对象将会写入当前目录下一个叫.RData的文件中,并且所有在这次会话中用过的命令行都会被保存在.Rhistory 的文件中。当R 再次在同一目录下启动,这些对象将从这个文件中重新导入工作空间。同时,相关的历史命令文件也会被导入。 4). 向量和赋值
向量是由一串有序数值构成的序列
x <- c(10.4, 5.6, 3.1, 6.4, 21.7)
函数c()完成的赋值语句。这里的函数c() 可以有任意多个参数,而它返回的值则是一个把这些参数首尾相连形成的向量。 赋值也可以用函数assign()实现。
assign(\
赋值符<-,->可以看作是该命令一个语义上的缩写。
c(10.4, 5.6, 3.1, 6.4, 21.7) -> x
向量运算
在算术表达式中使用向量将会对该向量的每一个元素都进行同样算术运算。
出现在同一个表达式中的向量最好是长度一致。如果他们的长度不一样,该表达式的值将是一个和其中最长向量等长的向量。 表达式中短的向量会被循环使用以达到最长向量的长度。 对于一个常数就是简单的重复。
v <- 2*x + y + 1
逻辑向量
逻辑向量元素可以被赋予的值,有TRUE,FALSE 和NA 逻辑向量可以由条件式(conditions)产生 temp <- x > 13 字符向量
字符向量就是字符串,可以用双引号和单引号作分割符。
paste():可以把单独的字符连成字符串,可以有任意多的参数。参数中的任何数字都将被显式地强制转换成字符串,而且以同样的方式在终端显示。默认的分隔符是单个的空格符。 修改分隔符换成”“
labs <- paste(c(\
> dx<-deriv(y~a*cos(a*x),\对一阶导函数求导 > dx(pi/3) [1] -1
attr(,\ x
[1,] -3.464102 # 导函数计算结果y'= -a^2*sin(a*x)=-2^2*sin(2*pi/3)=-3.464102
上面二阶导数的计算,我们是动手划分为两次求导进行计算的,利用deriv3()函数其实合并成一步计算。
> dx<-deriv3(y~sin(a*x),\生成二阶导数公式 > dx(pi/3) # 计算导数 [1] 0.8660254 attr(,\ x
[1,] -1 # 一阶导数结果 attr(,\, , x
x
[1,] -3.464102 # 二阶导数结果
我们再计算另外一个二阶导数,计算y=a*x^4+b*x^3+x^2+x+c,其中a,b,c为常数a=2,b=1,c=3,
根据导数计算公式,用于手动计算的变形结果为一阶导数为y’=2*x^4+x^3+x^2+x+3=4*2*x^3+3*x^2+2*x+1,当x=2时,y’=81, 对y’再求导公式变形为,y”=3*4*2*x^2+2*3*x+2,当x=2时,y”=110。
> dx<-deriv3(y~a*x^4+b*x^3+x^2+x+c,\通过func参数,指定常数值 > dx(2) [1] 49
attr(,\ x
[1,] 81 # 一阶导数结果 attr(,\, , x x
[1,] 110 # 二阶导数结果
这样就直接完成了二阶导数的计算,在R语言中二阶导数是可以直接求出的,想计算更高阶的导数就需要其他的数学工具包了。
4. 偏导数计算
在一元函数中,我们已经知道导数就是函数的变化率。对于二元函数我们同样要研究它的“变化率”。然而,由于自变量多了一个,情况就要复杂的多。在数学中,一个多变量的函数的偏导数,就是它关于其中一个变量的导数而保持其他变量恒定(相对于全导数,在其中所有变量都允许变化)。 偏导数的算子符号为:?。记作?f/?x 或者 f’x。偏导数反映的是函数沿坐标轴正方向的变化率,在向量分析和微分几何中是很有用的。
在xOy平面内,当动点由P(x0,y0)沿不同方向变化时,函数f(x,y)的变化快慢一般说来是不同的,因此就需要研究f(x,y)在(x0,y0)点处沿不同方向的变化率。在这里我们只学习函数f(x,y)在x0y平面沿着平行于x0y轴和平行于y轴两个特殊方位变动时,f(x,y)的变化率。 x方向的偏导:
设有二元函数z=f(x,y),点(x0,y0)是其定义域D内一点.把y固定在y0而让x在x0有增量△x,相应地函数z=f(x,y)有增量(称为对x的偏增量)△z=f(x0+△x,y0)-f(x0,y0)。如果△z与△x之比当△x→0时的极限存在,那么此极限值称为函数z=f(x,y)在(x0,y0)处对x的偏导数(partial derivative)。记作f’x(x0,y0)。 y方向的偏导:
函数z=f(x,y)在(x0,y0)处对x的偏导数,实际上就是把y固定在y0看成常数后,一元函数z=f(x,y0)在x0处的导数。同样,把x固定在x0,让y有增量△y,如果极限存在那么此极限称为函数z=(x,y)在(x0,y0)处对y的偏导数。记作f’y(x0,y0)
同样地,我们可以通过R语言的 deriv()函数进行偏导数的计算。下面我们计算一个二元函数f(x,y)=2*x^2+y+3*x*y^2的偏导数,由于二元函数曲面上每一点都有无穷多条切线,描述这个函数的导数就会相当困难。如果让其中的一个变量y取值为常数,那么就可以求出关于另一个自变量x的偏导数了,即?f/?x。
下面我们分别对x,y两个自变量求偏导数,设变量y为常数,计算x的偏导数?f/?x=4*x+3*y^2,当x=1,y=1时,x的偏导数?f/?x=4*x+3*y^2=7。设变量x为常数,计算y的偏导数?f/?y=1+6*x*y,当x=1,y=1时,y的偏导数?f/?x=1+6*x*y=7。 用R语言程序实现。
> fxy = expression(2*x^2+y+3*x*y^2) # 二元函数公式 > dxy = deriv(fxy, c(\> dxy
function (x, y) {
.expr4 <- 3 * x .expr5 <- y^2
.value <- 2 * x^2 + y + .expr4 * .expr5
.grad <- array(0, c(length(.value), 2L), list(NULL, c(\ .grad[, \ .grad[, \ attr(.value, \ .value }
> dxy(1,1) # 设置自变量 [1] 6
attr(,\
x y # 计算结果,x的偏导数为7,y的偏导数为7 [1,] 7 7
偏导数的程序计算结果与手动计算结果是一致的。下面我们再求一个复杂函数偏导数,计算一个二元函数f(x,y)=x^y + exp(x * y) + x^2 – 2 * x * y + y^3 + sin(x*y)在点(1,3)和点(0,0)的偏导数。 R语言程序实现。
> fxy = expression(x^y + exp(x * y) + x^2 - 2 * x * y + y^3 + sin(x*y)) > dxy = deriv(fxy, c(\> dxy(1,3) # 设置自变量 [1] 43.22666 attr(,\ x y
[1,] 56.28663 44.09554 # 计算结果,x的偏导数为56.28663,y的偏导数为 44.09554 > dxy(0,0) [1] 2
attr(,\ x y
[1,] NaN -Inf # 计算结果,x的偏导数无意义,y的偏导数负无穷大
对于计算的结果,有异议的同学,可以尝试动手计算。
本文我们掌握了R语言对于高等数学的导数计算方法,真的是非常方便,这下更有动力学习高数了。
概率基础和R语言
前方
R语言是统计语言,概率又是统计的基础,所以可以想到,R语言必然要从底层API上提供完整、方便、易用的概率计算的函数。让R语言帮我们学好概率的基础课。 目录 1. 2. 3.
随机变量
随机变量的数字特征 极限定理
1. 随机变量
? ? ?
什么是随机变量? 离散型随机变量 连续型随机变量
1). 什么是随机变量?
随机变量(random variable)表示随机现象各种结果的实值函数。随机变量是定义在样本空间S上,取值在实数载上的函数,由于它的自变量是随机试验的结果,而随机实验结果的出现具有随机性,因此,随机变量的取值具有一定的随机性。 R程序:生成一个在(0,1,2,3,4,5)的随机变量
> S<-1:5 > sample(S,1) [1] 2
> sample(S,1) [1] 3
> sample(S,1) [1] 5
2). 离散型随机变量
如果随机变量X的全部可能的取值只有有限多个或可列无穷多个,则称X为离散型随机变量。 R程序:生成样本空间为(1,2,3)的随机变量X,X的取值是有限的
> S<-1:3
> X<-sample(S,1);X [1] 2
3). 连续型随机变量
随机变量X,取值可以在某个区间内取任一实数,即变量的取值可以是连续的,这随机变量就称为连续型随机变量 R程序:生成样本在空间(0,1)的连续随机函数,取10个值
> runif(10,0,1)
[1] 0.3819569 0.7609549 0.6692581 0.6314708 0.5552201 0.8225527 0.7633086 0.4667188 0.1883553 [10] 0.3741653
2. 随机变量的数字特征
? ? ?
数学期望 方差 标准差
? ? ? ? ? ?
各种分步的期望和方差
常用统计量(最大,最小,中位数,四分位数) 协方差 相关系数
矩(原点矩,中心矩,偏度,峰度) 协方差矩阵
1). 数学期望(mathematical expectation)
离散型随机变量:的一切可能的取值xi与对应的概率Pi(=xi)之积的和称为该离散型随机变量的数学期望,记为E(x)。数学期望是最基本的数学特征之一。它反映随机变量平均取值的大小。
R程序:计算样本(1,2,3,7,21)的数学期望
> S<-c(1,2,3,7,21) > mean(S) [1] 6.8
连续型随机变量:若随机变量X的分布函数F(x)可表示成一个非负可积函数f(x)的积分,则称X为连续性随机变量,f(x)称为X的概率密度函数,积分值为X的数学期望,记为E(X)。
2). 方差(Variance)
方差是各个数据与平均数之差的平方的平均数。在概率论和数理统计中,方差用来度量随机变量和其数学期望(即均值)之间的偏离程度。 设X为随机变量,如果E{[X-E(X)]^2}存在,则称E{[X-E(X)]^2}为X的方差,记为Var(X)。
R程序:计算样本(1,2,3,7,21)的方差
> S<-c(1,2,3,7,21) > var(S) [1] 68.2
3). 标准差(Standard Deviation)
标准差是方差的算术平方根sqrt(var(X))。标准差能反映一个数据集的离散程度。平均数相同的,标准差未必相同。 R程序:计算样本(1,2,3,7,21)标准差
> S<-c(1,2,3,7,21) > sd(S) [1] 8.258329
4). 各种分步的期望和方差
? ?
离散型分布:两点分布,二项分布,泊松分布等
连续型分布:均匀分布,指数分布,正态分布,伽马分布等
对于某一特定场景,其所符合的分布规律一般先验给出 请参考文章:http://blog.fens.me/r-density/ 5). 常用统计量
众数(Mode): 一组数据中出现次数最多的数值,叫众数,有时众数在一组数中有好几个。 R程序:计算样本(1,2,3,3,3,7,7,7,7,9,10,21)的众数
> S<-c(1,2,3,3,3,7,7,7,7,9,10,21) > names(which.max(table(S))) [1] \
最小值(minimum): 在给定情形下可以达到的最小数量或最小数值 R程序:计算样本(2,3,3,3,7,7,7,7,9,10,21)的最小值
> S<-c(2,3,3,3,7,7,7,7,9,10,21) #最小值 > min(S) [1] 2 #最小值的索引 > which.min(S) [1] 1
最大值(maximum): 在给定情形下可以达到的最大数量或最大数值 R程序:计算样本(2,3,3,3,7,7,7,7,9,10,21)的最大值
> S<-c(2,3,3,3,7,7,7,7,9,10,21) #最大值 > max(S) [1] 21 #最大值的索引 > which.max(S) [1] 11
中位数(Medians): 是指将统计总体当中的各个变量值按大小顺序排列起来,形成一个数列,处于变量数列中间位置的变量值就称为中位数。 R程序:计算样本(1,2,3,4,5)的中位数
> S<-c(1,2,3,4,5) > median(S)
# 创建数据点
> x<-seq(-5,5,by=0.01) > y<-f1(x,a,b)
> df<-data.frame(x,y)
# 用ggplot2来画图
> g<-ggplot(df,aes(x,y))
> g<-g+geom_line(col='red') #红色直线
> g<-g+geom_point(aes(result$root,0),col=\点 > g<-g+geom_hline(yintercept=0)+geom_vline(yintercept=0) #坐标轴 > g<-g+ggtitle(paste(\> g
4.2 一元二次方程
一元二次方程:a*x^2+b*x+c=0,设a=1,b=5,c=6,求x?
> f2 <- function (x, a, b, c) a*x^2+b*x+c > a<-1;b<-5;c<-6
> result <- uniroot(f2,c(0,-2),a=a,b=b,c=c,tol=0.0001) > result$root [1] -2
把参数带入方程,用uniroot()函数,我们就解出了方程的一个根,改变计算的区间,我们就可以得到另一个根。
> result <- uniroot(f2,c(-4,-3),a=a,b=b,c=c,tol=0.0001) > result$root [1] -3
方程的两个根,一个是-2,一个是-3。
由于uniroot()函数,每次只能计算一个根,而且要求输入的区间端值,必须是正负号相反的。如果我们直接输入一个(-10,0)这个区间,那么uniroot()函数会出现错误。
> result <- uniroot(f2,c(-10,0),a=a,b=b,c=c,tol=0.0001)
Error in uniroot(f2, c(-10, 0), a = a, b = b, c = c, tol = 1e-04) : 位于极点边的f()值之正负号不相反
这应该是uniroot()为了统计计算对一元多次方程而设计的,所以为了使用uniroot()函数,我们需要取不同的区别来获得方程的根。 以图形展示方程:y = x^2 + 5*x + 6
# 创建数据点
> x<-seq(-5,1,by=0.01) > y<-f2(x,a,b,c) > df<-data.frame(x,y)
# 用ggplot2来画图
> g<-ggplot(df,aes(x,y))
> g<-g+geom_line(col='red') #红色曲线
> g<-g+geom_hline(yintercept=0)+geom_vline(yintercept=0) #坐标轴 > g<-g+ggtitle(paste(\> g
我们从图,并直接的看到了x的两个根取值范围。 4.3 一元三次方程
一元二次方程:a*x^3+b*x^2+c*x+d=0,设a=1,b=5,c=6,d=-11,求x?
> f3 <- function (x, a, b, c,d) a*x^3+b*x^2+c*x+d > a<-1;b<-5;c<-6;d<--11
> result <- uniroot(f3,c(-5,5),a=a,b=b,c=c,d=d,tol=0.0001) > result$root [1] 0.9461458
如果我们设置对了取值区间,那么一下就得到了方程的根。 以图形展示方程:y = x^2 + 5*x + 6
# 创建数据点
> x<-seq(-5,5,by=0.01) > y<-f3(x,a,b,c,d) > df<-data.frame(x,y)
# 用ggplot2画图
> g<-ggplot(df,aes(x,y))
> g<-g+geom_line(col='red') # 3次曲线
> g<-g+geom_hline(yintercept=0)+geom_vline(yintercept=0) #坐标轴 > g<-g+ggtitle(paste(\> g
4.4 二元一次方程组
R语言还可以解二次的方程组,当然计算方法,其实是利用于矩阵计算。 假设方程组:是以x1,x2两个变量组成的方程组,求x1,x2的值
以矩阵形式,构建方程组
# 左矩阵
> lf<-matrix(c(3,5,1,2),nrow=2,byrow=TRUE)
# 右矩阵
> rf<-matrix(c(4,1),nrow=2)
# 计算结果
> result<-solve(lf,rf) > result [,1] [1,] 3 [2,] -1
得方程组的解,x1, x2分别为3和-1。
接下来,我们画出这两个线性方程的图。设y=X2, x=X1,把原方程组变成两个函数形式。
# 定义2个函数
> fy1<-function(x) (-3*x+4)/5 > fy2<-function(x) (-1*x+1)/2
# 定义数据
> x<-seq(-1,4,by=0.01) > y1<-fy1(x) > y2<-fy2(x)
> dy1<-data.frame(x,y=y1,type=paste(\> dy2<-data.frame(x,y=y2,type=paste(\> df <- rbind(dy1,dy2)
# 用ggplot2画图
> g<-ggplot(df,aes(x,y))
> g<-g+geom_line(aes(colour=type,stat='identity')) #2条直线 > g<-g+geom_hline(yintercept=0)+geom_vline(yintercept=0) #坐标轴 > g
我们看到两条直线交点的坐标,就是方程组的两个根。多元一次方程,同样可以用这种方法来解得。
通过R语言,我们实现了对于初等数学的各种计算,真的是非常方便!下一篇文章将介绍,用R语言来解决高级数学中的计算问题。
前言
高等数学是每个大学生都要学习的一门数学基础课,同时也可能是考完试后最容易忘记的一门知识。
我在学习高数的时候绞尽脑汁,但始终都不知道为何而学。生活和工作基本用不到,就算是在计算机行业和金融行业,能直接用到高数的地方也少之又少,学术和实际应用真是相差太远了。
不过,R语言为我打开了一道高数应用的大门,R语言不仅能方便地实现高等数学的计算,还可以很容易地把一篇论文中的高数公式应用于产品的实践中。 因为R语言我重新学习了高数,让生活中充满数学,生活会变得更有意思。 本节并不是完整的高数计算手册,仅介绍了导数计算和偏导数计算的R语言实现。 目录 1. 2. 3. 4.
导数计算
初等函数的导数公式 二阶导数计算 偏导数计算
1. 导数计算
导数(Derivative)是微分学的基本概念,用于计算函数的极值。导数的定义为,当函数y=f(x)在x0的某个领域内有定义,当自变量x在x0处取得增加Δx(点x0+Δx仍在该邻域内)时,相应的函数取得增量Δy=f(x0+Δx)-f(x0);如果Δy与Δx之比当Δx趋于0时的极限存在,则称函数y=f(x)在点x0处可导,并称这个极限为函数y=f(x)在点x0处的导数,记为f`(x0),即
也记作 y’|x=x0 ,dy/dx|x=x0 或 df(x)/dx|x=x0。
通过R语言可以使用deriv()函数直接进行导数的计算,比如要计算 y=x^3 的导数,根据导数计算公式,用于手动计算的变形结果为 y’=3x^2,当x=1时,y’=3,当x=2时,y’=12。 本节的系统环境
? ?
Win7 64bit
R: 3.1.1 x86_64-w64-mingw32/x64 (64-bit)
用R语言程序实现,代码如下。
> dx <- deriv(y ~ x^3, \生成导数公式 expression({ .value <- x^3
.grad <- array(0, c(length(.value), 1L), list(NULL, c(\ .grad[, \
attr(.value, \ .value })
> mode(dx) # 查看dx变量类型 [1] \
> x<-1:2 # 给自变量x赋值 > eval(dx) # 运行求导计算 [1] 1 8 # 原函数的计算结果
attr(,\使用梯度下降法,导函数的计算结果 x
[1,] 3 # x=1,dx=3*1^2=3 [2,] 12 # x=2,dx=3*2^2=12
用R语言程序计算的结果,与我们手动计算的结果是一致的。但计算过程其实是有很大区别的,我们手动计算时是通过给定的导数计算公式,变成后完成的计算。而用计算机程序计算时,是使用梯度下降法来计算一阶导数,是一种最优化的近似算法。对于手动计算导数时,如果函数比较复杂而且比较难应用可变形的公式,那么手动计算就会有非常大的困难,而计算机程序的方法是一般地导数计算方法,不会受到公式难于变形的影响。
我们使用deriv(expr, name)函数时通常要传2个参数,第一参数expr就是原函数公式,用~号来分隔公式的两边,第二参数name用于指定函数的自变量。deriv()函数会返回一个表达式expression类型变量,再用eval()函数运行这个表达式得到就可得到计算结果,如上面的代码实现。 如果希望以函数的形式调用计算公式,那么你还需要传第三个参数func,并让func参数为TRUE,参考下面的代码实现。
计算正弦函数y=sin(x)的导数,根据导数计算公式,用于手动计算的变形结果为 y’=cos(x),当x=pi时,y’=-1,当x=4*pi时,y’=1,其中pi=π表示圆周率。
> dx <- deriv(y ~ sin(x), \生成导数公式的调用函数 function (x) {
.value <- sin(x)
.grad <- array(0, c(length(.value), 1L), list(NULL, c(\ .grad[, \
attr(.value, \ .value }
> mode(dx) # 检查dx的类型 [1] \
> dx(c(pi,4*pi)) # 以参数作为自变量,进行函数调用
[1] 1.224606e-16 -4.898425e-16 attr(,\
x # 导函数的计算结果 [1,] -1 # x=pi,dx=cos(pi)=-1 [2,] 1 # x=4*pi,dx=cos(4*pi)=1
2. 初等函数的导数公式
对于基本的初等函数求导数,通过导数计算公式是可以直接手动完成计算的,下面为一元初等函数的导数计算公式。
函数 原函数 导函数 常数函数 y=C y'=0
幂函数 y=x^n y'=n*x^(n-1) 指数函数 y=a^x y'=a^x*ln(a) y=exp(1)^x y'=exp(1)^x
对数函数 y=log(x,base=a) y'=1/(x*ln(a)) (a>0,且a!=1,x>0) y=ln(x) y'=1/x 正弦函数 y=sin(x) y'=cos(x) 余弦函数 y=cos(x) y'=-sin(x)
正切函数 y=tan(x) y'=sec(x)^2=1/cos(x)^2 余切函数 y=cot(x) y'=-csc(x)^2=1/sin(x)^2 正割函数 y=sec(x) y'=sec(x)*tan(x) 余割函数 y=csc(x) y'=-csc(x)*cot(x) 反正弦函数 y=arcsin(x) y'=1/sqrt(1-x^2) 反余弦函数 y=arccos(x) y'=-1/sqrt(1-x^2) 反正切函数 y=arctan(x) y'=1/(1+x^2) 反余切函数 y=arccot(x) y'=-1/(1+x^2) 反正割函数 y=arcsec(x) y'=1/abs(x)*(x^2-1) 反余割函数 y=arccsc(x) y'=-1/abs(x)*(x^2-1)
公式的注释:
? ? ? ? ? ? ? ? ? ?
y是原函数,x是y函数的自变量,y’是y函数的导函数。 C,n,a为常数。
ln表示以自然常数e为底的对数 exp(1)表示自然常数e
log(x,base=a)表示,以常数a为底的对数 sqrt表示开平方 abs表示绝对值
正割函数sec,计算方法为 sec=1/cos(x) 余割函数csc,计算方法为 csc=1/sin(x) 余切函数cot,计算方法为 cot=1/tan(x)
注: 以上公式不完全匹配于R语言函数
接下来,我们分别对这些一元初等函数进行一阶导数的计算。设y为原函数,x是y函数的自变量,且只有一个自变量。 常数函数
计算 y=3+10*x 函数的导数,根据导数计算公式,用于手动计算的变形结果为y’=0+10*x ,常数项3的导数为0,当x=1时,y’=10。
> dx<-deriv(y~ 3+10*x,\以函数形式生成导数公式
> dx(1) # 传入自变量,并计算
[1] 13 # 原函数计算结果y=3+10*1=13 attr(,\ x
[1,] 10 # 导函数计算结果y'=10*1=10
幂函数
计算 y=x^4 函数的导数,根据导数计算公式,用于手动计算的变形结果为y’=4*x^3,当x=2时,y’=32。
> dx<-deriv(y~x^4,\> dx(2) [1] 16
attr(,\ x
[1,] 32 # 导函数计算结果y'=4*x^3=4*2^3=32
指数函数
计算 y=4^x 函数的导数,根据导数计算公式,用于手动计算的变形结果为y’=4^x*ln(4),当x=2时,y’=22.18071。
> dx<-deriv(y~4^x ,\> dx(2) [1] 16
attr(,\ x
[1,] 22.18071 # 导函数计算结果y'=4^x*log(4)=4*2^3=22.18071
计算 y=exp(1)^x 函数的导数,根据导数计算公式,用于手动计算的变形结果为y’=exp(1)^x,当x=2时,y’=y=7.389056。
> dx<-deriv(y~exp(1)^x ,\> dx(2) [1] 7.389056 attr(,\ x
[1,] 7.389056 # 导函数计算结果y'=exp(1)^x=exp(1)^2=7.389056
对数函数
计算 y=ln(x) 函数的导数,根据导数计算公式,用于手动计算的变形结果为y’=1/x,当x=2时,y’=0.5。
> dx<-deriv(y~log(x),\> dx(2) [1] 0.6931472 attr(,\ x
[1,] 0.5 # 导函数计算结果y'=1/x=1/2=0.5
计算 y=log2(x) 函数的导数,根据导数计算公式,用于手动计算的变形结果为y’=1/(x*log(2)),当x=3时,y’=0.4808983。
但用R语言编程时,只能计算以自然常数为底的对数的导数,对于原函数不是以自然常数为底的对数,首先要变换成以自然常数为底的对数再进行导数计算,根据对数的换底公式,把以2为底的对数转换为以自然常数为底的对数 y=log2(x)=log(x)/log(2),
> dx<-deriv(y~log(x)/log(2),\> dx(3) [1] 1.584963 attr(,\ x
[1,] 0.4808983 # 导函数计算结果y'=1/(x*log(2)=1/(3*log(2)=0.4808983
正弦函数
计算 y=sin(x) 函数的导数,根据导数计算公式,用于手动计算的变形结果为y’=cos(x),当x=pi时,y’=-1,其中pi=π表示圆周率。
> dx<-deriv(y~sin(x),\> dx(pi)
[1] 1.224606e-16 attr(,\ x
[1,] -1 # 导函数计算结果y'=cos(x)=cos(pi)=-1
余弦函数
计算 y=cos(x) 函数的导数,根据导数计算公式,用于手动计算的变形结果为y’=-sin(x),当x=pi/2时,y’=-1。
> dx<-deriv(y~cos(x),\> dx(pi/2) [1] 6.123032e-17 attr(,\ x
[1,] -1 # 导函数计算结果y'=-sin(x)=-sin(pi/2)=-1
正切函数
计算 y=tan(x) 函数的导数,根据导数计算公式,用于手动计算的变形结果为 y’=sec(x)^2=1/cos(x)^2,当x=pi/6时,y’=1.333333。
> dx<-deriv(y~tan(x),\> dx(pi/6) [1] 0.5773503 attr(,\ x
[1,] 1.333333 # 导函数计算结果y'=1/cos(x)^2=1/cos(pi/6)^2=1.333333
余切函数
计算 y=cot(x) 函数的导数,由于R语言没有cot()函数,所以根据三角公式我们动手变形原函数为y=cot(x)=1/tan(x)后再进行导数计算,根据导数计算公式,用于手动计算的变形结果为y’=-csc(x)^2=-1/sin(x)^2,当x=pi/6时,y’=-4。
> dx<-deriv(y~1/tan(x),\> dx(pi/6) [1] 1.732051 attr(,\ x
[1,] -4 # 导函数计算结果y'=-1/sin(x)^2=-1/sin(pi/6)^2=-4
反正弦函数
计算 y=asin(x) 函数的导数,根据导数计算公式,用于手动计算的变形结果为y’=1/sqrt(1-x^2),当x=pi/6时,y’=1.173757。
> dx<-deriv(y~asin(x),\> dx(pi/6) [1] 0.5510696 attr(,\ x
[1,] 1.173757 # 导函数计算结果y'=1/sqrt(1-x^2)=1/sqrt(1-(pi/6)^2)=1.173757
反余弦函数
计算 y=acos(x) 函数的导数,根据导数计算公式,用于手动计算的变形结果为y’=-1/sqrt(1-x^2),当x=pi/8时,y’=-1.08735。
> dx<-deriv(y~acos(x),\> dx(pi/8) [1] 1.167232 attr(,\ x
[1,] -1.08735 # 导函数计算结果y'=-1/sqrt(1-x^2)=-1/sqrt(1-(pi/8)^2)=-1.08735
反正切函数
计算 y=atan(x) 函数的导数,根据导数计算公式,用于手动计算的变形结果为y’=1/(1+x^2),当x=pi/6时,y’=0.7848335。
> dx<-deriv(y~atan(x),\> dx(pi/6) [1] 0.4823479 attr(,\ x
[1,] 0.7848335 # 导函数计算结果y'= 1/(1+x^2) = 1/(1+(pi/6)^2)=0.7848335
3. 二阶导数计算
当我们对一个函数进行多次接连的求导计算,会形成高阶导数。
一般的,函数y=f(x)的导数y’=f'(x)仍然是x的函数,我们就把y’=f'(x)的导数叫做函数y=f(x)的二阶导数,记作y”,即
一阶导数的导数叫做二阶导数,二阶导数的导数叫做三阶导数,N-1阶导数的导数叫做N阶导数,习惯上把二阶以上的导数称之为高阶导数,
比如,计算 y=sin(a*x) 函数的二阶导数导数y”,其中a为常数,根据导数计算公式,用于手动计算的变形结果为一阶导数为y’=a*cos(a*x),对y’再求导公式变形为,y”=-a^2*sin(a*x) 用R语言进行程序实现
> a<-2 # 设置a的值
> dx<-deriv(y~sin(a*x),\生成一阶导数公式 > dx(pi/3) # 计算一阶导数 [1] 0.8660254 attr(,\ x
[1,] -1 # 导函数计算结果y'= a*cos(a*x)=2*cos(2*pi/3)=-1