R語言的一個小范例：數(shù)值模擬和繪圖_數(shù)據(jù)分析師培訓

 下午上課，要給經(jīng)濟工程專業(yè)的本科生簡要介紹一下用于統(tǒng)計分析和繪圖的R語言?？紤]到他們所在的年級并沒有系統(tǒng)的學完統(tǒng)計學和計量經(jīng)濟學，打算不從回歸分析入手，先教他們對R的基本操作有個了解。寫了如下這段代碼，準備課上帶同學們做一下。順便發(fā)到網(wǎng)上，供大家參考。

背景：

       假定是一個由決定的量，滿足關(guān)系       問在的取值空間內(nèi)，值的走勢如何，以及最大值出現(xiàn)在什么地方，并繪圖說明。

talk is cheap. show me the code （廢話少說，放碼過來）

       一點一點寫，并適當做說明

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

       這段代碼的意思是，定義一個名為qiu的function，輸入值alpha之后，會自動求得對應的值。比如

> qiu(0.67)
[1] 0.622206

       現(xiàn)在我們來搞alpha。

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

       定義一個數(shù)組x，取值從0.01開始，每次增加0.001單位，直到0.999。然后看看x這個數(shù)組的情況吧

x
length(x)

       結(jié)果就不貼出來了。現(xiàn)在，把x的值交給qiu來處理，生成的一系列值，我們定義為y。

y <- qiu(x)

       畫個草圖來看看分布情況如何

plot(y~x)

       顯然，隨著x從0到1, y并不是單調(diào)變化的：先增大，后減小。那么，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，供下文使用。

       重新畫張圖，前面那張?zhí)罅?/span>

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的繪圖程序，關(guān)閉此前的圖。第二行的ylim和xlim定義x和y軸的取值范圍。type是告訴plot程序，散點圖以line的形式呈現(xiàn)。ylab和xlab是x和y軸的名稱：由于我們需要使用希臘字母，因此需要使用expression(alpha)和expression(beta)來打出與。

       比原來那張圖看著帥一些了是吧。繼續(xù)調(diào)整。

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="")

       第一行告訴繪圖程序，我下面要在原有那張圖的基礎(chǔ)上，繼續(xù)添加內(nèi)容，不要覆蓋原圖。

       第二行中，把y最大值時的點標出來，pch表示這個點的形狀（1是空心圓，2是空心三角，等等，你們自己試試看從1到20吧）。col="red"，將這個點畫成紅色。ylab和xlab設為空，ylim和xlim和上文的值相同，使得兩張圖x軸、y軸重合。

       差不多行了。如果你想繼續(xù)和我一樣騷包一點的話，跟我向下繼續(xù)設置。

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")
)
)

       收工。這張圖看起來還是有些難看。。。。不過基本信息都已經(jīng)有了。