R语言的一个小范例：数值模拟和绘图_数据分析师培训

 下午上课，要给经济工程专业的本科生简要介绍一下用于统计分析和绘图的R语言。考虑到他们所在的年级并没有系统的学完统计学和计量经济学，打算不从回归分析入手，先教他们对R的基本操作有个了解。写了如下这段代码，准备课上带同学们做一下。顺便发到网上，供大家参考。

背景：

       假定是一个由决定的量，满足关系       问在的取值空间内，值的走势如何，以及最大值出现在什么地方，并绘图说明。

talk is cheap. show me the code （废话少说，放码过来）

       一点一点写，并适当做说明

qiu <- function(alpha){(alpha / (1-alpha))^{-alpha}}

       这段代码的意思是，定义一个名为qiu的function，输入值alpha之后，会自动求得对应的值。比如

> qiu(0.67)
[1] 0.622206

       现在我们来搞alpha。

x <- seq(0.001,0.99, by=0.001)

       定义一个数组x，取值从0.01开始，每次增加0.001单位，直到0.999。然后看看x这个数组的情况吧

x
length(x)

       结果就不贴出来了。现在，把x的值交给qiu来处理，生成的一系列值，我们定义为y。

y <- qiu(x)

       画个草图来看看分布情况如何

plot(y~x)

       显然，随着x从0到1, y并不是单调变化的：先增大，后减小。那么，y的最高值是多少？对应的x值是多大？

max(y)
x[which.max(y)]

       第一行告诉你y的最大值1.3211

       第二行的which.max(y)告诉你当y值最大时，所对应的编号是218。把218带入x中，找到对应的x的值：x[218]=0.218

schx=c(x[which.max(y)])
schy=c(max(y))

       分别将y值最大时的x和y值，定义为schx和schy，供下文使用。

       重新画张图，前面那张太丑了

dev.off()
plot(y~x, ylim=c(0,1.5), xlim=c(0,1), type="l", ylab=expression(beta), xlab=expression(alpha))

       第一行dev.off()告诉R的绘图程序，关闭此前的图。第二行的ylim和xlim定义x和y轴的取值范围。type是告诉plot程序，散点图以line的形式呈现。ylab和xlab是x和y轴的名称：由于我们需要使用希腊字母，因此需要使用expression(alpha)和expression(beta)来打出与。

       比原来那张图看着帅一些了是吧。继续调整。

par(new=TRUE)
plot(max(y)~x[which.max(y)], pch=2, lty=3, col="red",ylim=c(0,1.5), xlim=c(0,1), ylab="", xlab="")

       第一行告诉绘图程序，我下面要在原有那张图的基础上，继续添加内容，不要覆盖原图。

       第二行中，把y最大值时的点标出来，pch表示这个点的形状（1是空心圆，2是空心三角，等等，你们自己试试看从1到20吧）。col="red"，将这个点画成红色。ylab和xlab设为空，ylim和xlim和上文的值相同，使得两张图x轴、y轴重合。

       差不多行了。如果你想继续和我一样骚包一点的话，跟我向下继续设置。

text(schx, schy-0.2, substitute(paste("(", xx ,", ", yy, ")"), list(xx=schx,yy=schy)))

       这行代码的作用是，在途中加入一段文字。语法text(a,b,c)的意思是，a代表横坐标的位置，b代表纵坐标的位置，c代表所需要加入的语言。a和b加在一起，告诉plot()需要把一段文字c放在(a,b)这个地方。

       c需要做进一步的说明。"(", 以及 ")", 表示这其中是有纯文字部分的，plot()直接把它们打出来即可（注意逗号是要保留的）。加在两个小括号中间的是xx和yy，这是两个值。

              substitute(paste(),list())：

       substitute()中包括两部分，一部分是paste(),一部分是list()。paste告诉plot()要加入这段东西了，其中包括两个待赋值的xx和yy。list()为它们分别赋值。

       看看效果吧。

       看起来还可以。要是能把       这个式子也写上去就更完美了。

text(0.4, 0.5, 
substitute(
paste(beta == (frac(alpha, 1-alpha))^{-alpha}, 
",", 
"0<", alpha, "<1")
)
)

       收工。这张图看起来还是有些难看。。。。不过基本信息都已经有了。