下载掌阅APP，畅读海量书库
立即打开
如何用R语言进行数学计算？
R语言是统计语言，生来就对数学有良好的支持，用R语言做数学的计算题特别方便。如果计算器中能嵌入R语言的计算函数，那么绝对是一种高科技产品。我真的把R语言当成我的计算器了！
R语言对数学计算有着非常好的支持，本节将完整介绍初等数学中的各种计算操作。
本节的系统环境是：
·Windows 7 64bit
·R：3.1.1 x86_64-w64-mingw32/x64（64-bit）
用R语言实现四则运算操作，包括加、减、乘、除、余数、整除、绝对值、判断正负。
> a<-10；b<-5 # 定义2个变量 > a+b；a-b；a*b；a/b # 加减乘除 [1] 15 [1] 5 [1] 50 [1] 2 > a%%b；a%/%b # 余数，整除 [1] 0 [1] 2 > abs（-a） # 绝对值 [1] 10 > sign（-2：3） # 判断正负 [1] -1 -1 0 1 1 1
用R语言实现数学计算操作，包括幂、自然常数e的幂、平方根、对数。
> a<-10；b<-5；c<-4 # 定义3个变量 > c^b；c^-b；c^（b/10） # 幂运算 [1] 1024 [1] 0.0009765625 [1] 2 > exp（1） # 取自然常数e [1] 2.718282 > exp（3） # 自然常数e的幂 [1] 20.08554 > sqrt（c） # 平方根 [1] 2 > log2（c） # 以2为底的对数 [1] 2 > log10（b） # 以10为底的对数 [1] 0.69897 > log（c，base = 2） # 自定义底的对数 [1] 2 > log（a，base=exp（1）） # 自然常数e的对数 [1] 2.302585 > log（a^b，base=a） # 指数对数操作 [1] 5 > log（exp（3）） [1] 3
用R语言实现比较计算操作，包括==、>、<、！=、<=、>=、isTRUE、identical。
> a<-10；b<-5 # 定义2个变量 > a==a；a！=b；a>b；a<b；a<=b；a>=c # 比较计算 [1] TRUE [1] TRUE [1] TRUE [1] FALSE [1] FALSE [1] TRUE > isTRUE（a） # 判断是否为TRUE [1] FALSE > isTRUE（！a） [1] FALSE > identical（1， as.integer（1）） # 精确比较两个对象 [1] FALSE > identical（NaN， -NaN） [1] TRUE > f <- function（x） x > g <- compiler：：cmpfun（f） > identical（f， g） [1] TRUE
用R语言实现逻辑计算操作，包括&、|、&&、||、xor。
> x<-c（0，1，0，1） # 定义2个向量 > y<-c（0，0，1，1） > x && y；x || y # 只比较向量的第一个元素 &&， || [1] FALSE [1] FALSE > x & y；x | y # S4对象的逻辑运算，比较所有元素 &， | [1] FALSE FALSE FALSE TRUE [1] FALSE TRUE TRUE TRUE > xor（x，y） # 异或比较 [1] FALSE TRUE TRUE FALSE > xor（x，！y） [1] TRUE FALSE FALSE TRUE
用R语言实现约数计算操作，包括ceiling、floor、trunc、round、signif。
> ceiling（5.4） # 向上取整 [1] 6 > floor（5.8） # 向下取整 [1] 5 > trunc（3.9） # 取整数 [1] 3 > round（5.8） # 四舍五入 [1] 6 > round（5.8833， 2） # 四舍五入，并保留2位小数 [1] 5.88 > signif（5990000，2） # 四舍五入，保留前2位整数 [1] 6e+06
用R语言实现数组计算操作，包括求最大值、求最小值、范围、求和、均值、加权平均、连乘、差分、秩、中位数、分位数、任意数、全体数。
> d<-seq（1，10，2）；d # 定义1个向量 [1] 1 3 5 7 9 > max（d）；min（d）；range（d） # 求最大值、最小值、范围range [1] 9 [1] 1 [1] 1 9 > sum（d）；mean（d） # 求和、均值 [1] 25 [1] 5 > weighted.mean（d，rep（1，5）） # 加权平均 [1] 5 > weighted.mean（d，c（1，1，2，2，2）） [1] 5.75 > prod（1：5） # 连乘 [1] 120 > diff（d） # 差分 [1] 2 2 2 2 > rank（d） # 秩 [1] 1 2 3 4 5 > median（d） # 中位数 [1] 5 > quantile（d） # 分位数 0% 25% 50% 75% 100% 1 3 5 7 9 > any（d<5）；all（d<5） # 任意条件any，全体条件all [1] TRUE [1] FALSE
用R语言实现排列组合计算操作，包括阶乘、组合、排列。
> factorial（5） # 阶乘5！
[1] 120
> choose（5， 2） # 组合，从5个中选出2个
[1] 10
> combn（5，2） # 列出从5个中选出2个的所有组合项
[，1] [，2] [，3] [，4] [，5] [，6] [，7] [，8] [，9] [，10]
[1，] 1 1 1 1 2 2 2 3 3 4
[2，] 2 3 4 5 3 4 5 4 5 5
> for （n in 0：10） print（choose（n， k = 0：n）） # 计算0：10的组合个数
[1] 1
[1] 1 1
[1] 1 2 1
[1] 1 3 3 1
[1] 1 4 6 4 1
[1] 1 5 10 10 5 1
[1] 1 6 15 20 15 6 1
[1] 1 7 21 35 35 21 7 1
[1] 1 8 28 56 70 56 28 8 1
[1] 1 9 36 84 126 126 84 36 9 1
[1] 1 10 45 120 210 252 210 120 45 10 1
> choose（5， 2）*factorial（2） # 排列，从5个中选出2个
[1] 20
用R语言实现累积计算操作，包括累加、累乘、最小累积、最大累积。
> cumsum（1：5） # 累加 [1] 1 3 6 10 15 > cumprod（1：5） # 累乘 [1] 1 2 6 24 120 > e<-seq（-3，3）；e # 定义一个向量 [1] -3 -2 -1 0 1 2 3 > cummin（e） # 最小累积cummin [1] -3 -3 -3 -3 -3 -3 -3 > cummax（e） # 最大累积cummax [1] -3 -2 -1 0 1 2 3
用R语言实现两个数组的计算操作，包括取交集、并集、差集、数组是否相等、取唯一、查匹配元素的索引、找重复元素索引。
> x <- c（9：20， 1：5， 3：7， 0：8）；x # 定义两个数组向量 [1] 9 10 11 12 13 14 15 16 17 18 19 20 1 2 3 4 5 [18] 3 4 5 6 7 0 1 2 3 4 5 6 7 8 > y<- 1：10；y [1] 1 2 3 4 5 6 7 8 9 10 > intersect（x，y） # 交集 [1] 9 10 1 2 3 4 5 6 7 8 > union（x，y） # 并集 [1] 9 10 11 12 13 14 15 16 17 18 19 20 1 2 3 4 5 [18] 6 7 0 8 > setdiff（x，y） # 差集，从x中排除y [1] 11 12 13 14 15 16 17 18 19 20 0 > setequal（x， y） # 判断是否相等 [1] FALSE > unique（c（x，y）） # 取唯一 [1] 9 10 11 12 13 14 15 16 17 18 19 20 1 2 3 4 5 [18] 6 7 0 8 > which（x %in% y） # 找到x在y中存在的元素的索引 [1] 1 2 13 14 15 16 17 18 19 20 21 22 24 25 26 27 28 [18] 29 30 31 > which（is.element（x，y）） # 同%in% [1] 1 2 13 14 15 16 17 18 19 20 21 22 24 25 26 27 28 [18] 29 30 31 > which（duplicated（x）） # 找到重复元素的索引 [1] 18 19 20 24 25 26 27 28 29 30
1.三角函数
在直角三角形中仅有锐角（大小在0~90度之间的角）三角函数的定义。给定一个锐角θ，可以做出一个直角三角形，使得其中的一个内角是θ。设这个三角形中，θ的对边、邻边和斜边长度分别是a、b和h，如图1-3所示。
图1-3 直角三角形
三角函数的6种关系：正弦、余弦、正切、余切、正割、余割。
·θ的正弦是对边与斜边的比值：sinθ=a/h
·θ的余弦是邻边与斜边的比值：cosθ=b/h
·θ的正切是对边与邻边的比值：tanθ=a/b
·θ的余切是邻边与对边的比值：cotθ=b/a
·θ的正割是斜边与邻边的比值：secθ=h/b
·θ的余割是斜边与对边的比值：cscθ=h/a
三角函数的特殊值。
用R语言实现三角基本函数计算，包括正弦、余弦、正切。
> sin（0）；sin（1）；sin（pi/2） # 正弦 [1] 0 [1] 0.841471 [1] 1 > cos（0）；cos（1）；cos（pi） # 余弦 [1] 1 [1] 0.5403023 [1] -1 > tan（0）；tan（1）；tan（pi） # 正切 [1] 0 [1] 1.557408 [1] -1.224647e-16
接下来，我们用ggplot2包来画出三角函数的图形。
> library（ggplot2） # 加载ggplot2的库 > library（scales）
三角函数画图，以下代码生成三角函数曲线，如图1-4所示。
> x<-seq（-2*pi，2*pi，by=0.01） # x坐标
> s1<-data.frame（x，y=sin（x），type=rep（'sin'，length（x））） # y坐标，正弦
> s2<-data.frame（x，y=cos（x），type=rep（'cos'，length（x））） # y坐标，余弦
> s3<-data.frame（x，y=tan（x），type=rep（'tan'，length（x））） # y坐标，正切
> s4<-data.frame（x，y=1/tan（x），type=rep（'cot'，length（x））） # y坐标，余切
> s5<-data.frame（x，y=1/cos（x），type=rep（'sec'，length（x））） # y坐标，正割
> s6<-data.frame（x，y=1/sin（x），type=rep（'csc'，length（x））） # y坐标，余割
> df<-rbind（s1，s2，s3，s4，s5，s6）
> g<-ggplot（df，aes（x，y）） # 用ggplot2画图
> g<-g+geom_line（aes（colour=type，stat='identity'））
> g<-g+scale_y_continuous（limits=c（-2， 2））
> g<-g+scale_x_continuous（breaks=seq（-2*pi，2*pi，by=pi），labels=c（"-2*pi"，
"-pi"，"0"，"pi"，"2*pi"））
> g
图1-4 三角函数曲线
2.反三角函数
基本的反三角函数定义如下。
用R语言实现反三角函数的计算，包括反正弦、反余弦、反正切。
> asin（0）；asin（1） # 反正弦asin [1] 0 [1] 1.570796 # pi/2=1.570796 > acos（0）；acos（1） # 反余弦acos [1] 1.570796 # pi/2=1.570796 [1] 0 > atan（0）；atan（1） # 反正切atan [1] 0 [1] 0.7853982 # pi/4=0.7853982
反三角函数画图，以下代码生成反三角函数曲线，如图1-5所示。
> x<-seq（-1，1，by=0.005） # x坐标
> s1<-data.frame（x，y=asin（x），type=rep（'arcsin'，length（x））） # y坐标，反正弦
> s2<-data.frame（x，y=acos（x），type=rep（'arccos'，length（x））） # y坐标，反余弦
> s3<-data.frame（x，y=atan（x），type=rep（'arctan'，length（x））） # y坐标，反正切
> s4<-data.frame（x，y=1/atan（x），type=rep（'arccot'，length（x））） # y坐标，反余弦
> s5<-data.frame（x，y=1/asin（x），type=rep（'arcsec'，length（x））） # y坐标，反正割
> s6<-data.frame（x，y=1/acos（x），type=rep（'arccsc'，length（x））） # y坐标，反余弦
> df<-rbind（s1，s2，s3，s4，s5，s6）
> g<-ggplot（df，aes（x，y）） # 用ggplot2画图
> g<-g+geom_line（aes（colour=type，stat='identity'））
> g<-g+scale_y_continuous（limits=c（-2*pi，2*pi），breaks=seq（-2*pi，2*pi，by=pi），
labels=c（"-2*pi"，"-pi"，"0"，"pi"，"2*pi"））
> g
图1-5 反三角函数曲线
3.三角函数公式
接下来，用单元测试的方式来描述三角函数的数学公式，公式的左边等于公式的右边。通过testthat包进行单元测试，关于testthat包的安装和使用，请参考5.2节。
使用expect_that（right，left）函数，把公式的左右两边表达式，分别以参数形式传入函数中。运行expect_that（）函数，如果没有返回结果则表示两个参数相等，如果有输出则根据输出查看原因。
> library（testthat） # 加载testthat包 > a<-5；b<-10 # 定义变量
平方和公式：
> sin（a）^2+cos（a）^2 [1] 1 > expect_that（1， equals（sin（a）^2+cos（a）^2）） # 用单元测试的方法，判断公式是否两边相等 > expect_that（2， equals（sin（a）^2+cos（a）^2）） # 如果公式两边不相等，会有错误提示 Error： 2 not equal to sin（a）^2 + cos（a）^2 Mean relative difference： 1
和角公式：
和角公式的单元测试如下：
> expect_that（sin（a）*cos（b）+sin（b）*cos（a），equals（sin（a+b））） > expect_that（sin（a）*cos（b）-sin（b）*cos（a），equals（sin（a-b））） > expect_that（cos（a）*cos（b）-sin（b）*sin（a），equals（cos（a+b））） > expect_that（cos（a）*cos（b）+sin（b）*sin（a），equals（cos（a-b））） > expect_that（（tan（a）+tan（b））/（1-tan（a）*tan（b）），equals（tan（a+b））） > expect_that（（tan（a）-tan（b））/（1+tan（a）*tan（b）），equals（tan（a-b）））
2倍角公式：
2倍角公式的单元测试如下：
> expect_that（cos（a）^2-sin（a）^2，equals（cos（2*a））） > expect_that（2*cos（a）^2-1，equals（cos（2*a））） > expect_that（1-2*sin（a）^2，equals（cos（2*a）））
3倍角公式：
3倍角公式的单元测试如下：
> expect_that（4*cos（a）^3-3*cos（a），equals（cos（3*a））） > expect_that（-4*sin（a）^3+3*sin（a），equals（sin（3*a）））
半角公式：
半角公式的单元测试如下：
> expect_that（sqrt（（1-cos（a））/2），equals（abs（sin（a/2）））） > expect_that（sqrt（（1+cos（a））/2），equals（abs（cos（a/2）））） > expect_that（sqrt（（1-cos（a））/（1+cos（a））），equals（abs（tan（a/2）））） > expect_that（abs（sin（a）/（1+cos（a））），equals（abs（tan（a/2）））） > expect_that（abs（（1-cos（a））/sin（a）），equals（abs（tan（a/2））））
和差化积：
和差化积公式的单元测试如下：
> expect_that（（sin（a+b）+sin（a-b））/2，equals（sin（a）*cos（b））） > expect_that（（sin（a+b）-sin（a-b））/2，equals（cos（a）*sin（b））） > expect_that（（cos（a+b）+cos（a-b））/2，equals（cos（a）*cos（b））） > expect_that（（cos（a-b）-cos（a+b））/2，equals（sin（a）*sin（b）））
积化和差：
积化和差公式的单元测试如下：
> expect_that（sin（a）+sin（b），equals（2*sin（（a+b）/2）*cos（（a-b）/2））） > expect_that（sin（a）-sin（b），equals（2*cos（（a+b）/2）*sin（（a-b）/2））） > expect_that（2*cos（（a+b）/2）*cos（（a-b）/2），equals（cos（a）+cos（b））） > expect_that（-2*sin（（a+b）/2）*sin（（a-b）/2），equals（cos（a）-cos（b）））
万能公式：
万能公式的单元测试如下：
> expect_that（sin（2*a），equals（2*tan（a）/（1+tan（a）^2））） > expect_that（（1-tan（a）^2）/（1+tan（a）^2），equals（cos（2*a））） > expect_that（2*tan（a）/（1-tan（a）^2），equals（tan（2*a）））
平方差公式：
平方差公式的单元测试如下：
> expect_that（sin（a）^2-sin（b）^2，equals（sin（a+b）*sin（a-b））） > expect_that（cos（a）^2-sin（b）^2，equals（cos（a+b）*cos（a-b）））
降次升角公式：
降次升角公式的单元测试如下：
> expect_that（（1+cos（2*a））/2，equals（cos（a）^2）） > expect_that（（1-cos（2*a））/2，equals（sin（a）^2））
辅助角公式：
辅助角公式的单元测试如下：
> expect_that（sqrt（a^2+b^2）*sin（a+atan（b/a）），equals（a*sin（a）+b*cos（a）））
复数为实数的延伸，它使任一多项式都有根。复数中的虚数单位i，是-1的一个平方根，即i 2 =-1。任一复数都可表达为x+yi，其中x及y皆为实数，分别称为复数之“实部”和“虚部”。
1.创建一个复数
> ai<-5+2i；ai # 直接创建复数 [1] 5+2i > class（ai） # 查看复数的类型 [1] "complex" > bi<-complex（real=5，imaginary=2）；bi # 通过complex（）函数创建复数 [1] 5+2i > is.complex（bi） [1] TRUE > Re（ai） # 实数部分 [1] 5 > Im（ai） # 虚数部分 [1] 2 > Mod（ai） # 取模 [1] 5.385165 # sqrt（5^2+2^2） = 5.385165 > Arg（ai） # 取辐角 [1] 0.3805064 > Conj（ai） # 取轭 [1] 5-2i
2.复数四则运算
·加法公式：（a+bi）+（c+di）=（a+c）+（b+d）i
·减法公式：（a+bi）-（c+di）=（a-c）+（b-d）i
·乘法公式：（a+bi）（c+di）=ac+adi+bci+bidi=ac+bdi 2 +（ad+bc）i=（ac-bd）+（ad+bc）i
·除法公式：（a+bi）/（c+di）=（（ac+bd）+（bc-ad）i）/（c 2 +d 2 ）
> a<-5；b<-2；c<-3；d<-4 > ai<-complex（real=a，imaginary=b） > bi<-complex（real=c，imaginary=d）
复数四则运算的单元测试如下：
> expect_that（complex（real=（a+c），imaginary=（b+d）），equals（ai+bi）） > expect_that（complex（real=（a-c），imaginary=（b-d）），equals（ai-bi）） > expect_that（complex（real=（a*c-b*d），imaginary=（a*d+b*c）），equals（ai*bi）） > expect_that（complex（real=（a*c+b*d），imaginary=（b*c-a*d））/（c^2+d^2），equals（ai/bi））
3.复数开平方根
> sqrt（-9） # 在实数域，给-9开平方根 [1] NaN Warning message： In sqrt（-9） ： NaNs produced > sqrt（complex（real=-9）） # 在复数域，给-9开平方根 [1] 0+3i
方程计算是数学计算的一种基本形式，R语言也可以很方便地帮助我们解方程，下面将介绍一元多次方程和二元一次方程的解法。解一元多次方程，可以用uniroot（）函数。
1.一元一次方程
一元一次方程：ax+b=0，设a=5，b=10，求x？
> f1 <- function （x， a， b） a*x+b # 定义方程函数
> a<-5；b<-10 # 给a，b常数赋值
> result <- uniroot（f1，c（-10，10），a=a，b=b，tol=0.0001） # 在（-10，10）的区间，精确度
# 为0.0001位，计算方程的根
> result$root # 打印方程的根x
[1] -2
一元一次方程非常容易解得，方程的根是-2！以图形展示函数：y=5x+10，如图1-6所示。
> x<-seq（-5，5，by=0.01） # 创建数据点 > y<-f1（x，a，b） > df<-data.frame（x，y） > g<-ggplot（df，aes（x，y）） # 用ggplot2来画图 > g<-g+geom_line（col='red'） # 红色直线 > g<-g+geom_point（aes（result$root，0），col="red"，size=3） # 点 > g<-g+geom_hline（yintercept=0）+geom_vline（yintercept=0） # 坐标轴 > g<-g+ggtitle（paste（"y ="，a，"* x +"，b）） > g
图1-6 函数y=5x+10
2.一元二次方程
一元二次方程：ax 2 +bx+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 +5x+6，如图1-7所示。 [1]
> x<-seq（-5，1，by=0.01） # 创建数据点 > y<-f2（x，a，b，c） > df<-data.frame（x，y） > g<-ggplot（df，aes（x，y）） # 用ggplot2来画图 > g<-g+geom_line（col='red'） # 红色曲线 > g<-g+geom_hline（yintercept=0）+geom_vline（yintercept=0） # 坐标轴 > g<-g+ggtitle（paste（"y ="，a，"* x ^ 2 +"，b，"* x +"，c）） > g
图1-7 函数y=x 2 +5x+6
我们从图1-7中可以很明显看到x的两个根的取值范围。
3.一元三次方程
一元三次方程：ax 3 +bx 2 +cx+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 3 +5x 2 +6x-11，如图1-8所示。
> x<-seq（-5，5，by=0.01） # 创建数据点 > y<-f3（x，a，b，c，d） > df<-data.frame（x，y） > g<-ggplot（df，aes（x，y）） # 用ggplot2画图 > g<-g+geom_line（col='red'） # 3次曲线 > g<-g+geom_hline（yintercept=0）+geom_vline（yintercept=0） # 坐标轴 > g<-g+ggtitle（paste（"y ="，a，"* x ^ 3 +"，b，"* x ^2 +"，c，"* x + "，d）） > g
图1-8 函数y=x 3 +5x 2 +6x-11
4.二元一次方程组
R语言还可以用来解二元方程组，当然计算方法其实是利用了矩阵计算。
下面是x 1 ，x 2 两个未知变量组成的方程组，求x 1 ，x 2 的值。
以矩阵形式来构建方程组就是
> 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
得方程组的解，x 1 ，x 2 分别为3和-1。
接下来，我们画出这两个线性方程的图，如图1-9所示。设y=x 2 ，x=x 1 ，把原方程组变成两个函数形式。
> fy1<-function（x） （-3*x+4）/5 # 定义2个函数 > 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（"y=（-3*x+4）/5"）） > dy2<-data.frame（x，y=y2，type=paste（"y=（-1*x+1）/2"）） > df <- rbind（dy1，dy2） > g<-ggplot（df，aes（x，y）） # 用ggplot2画图 > g<-g+geom_line（aes（colour=type，stat='identity'）） # 2条直线 > g<-g+geom_hline（yintercept=0）+geom_vline（yintercept=0） # 坐标轴 > g
图1-9 二元一次方程组
我们看到两条直线交点的坐标就是方程组的两个根。多元一次方程同样可以用这种方法来解得。
通过R语言，我们实现了对初等数学的各种计算，真的是非常方便。
[1] 请读者自己运行R语言代码，查看生成的彩色图。