物理專業(yè)應該掌握的計算機技能No4——運用R做數(shù)據(jù)分析.pdf

物理學家應該掌握的計算機技能No 4 運用R 做數(shù)據(jù)分析 author fame 1 1 什么是什么是 R R 在R的官方教程里是這么給R下注解的 一個數(shù)據(jù)分析和圖形顯示的程序設計環(huán)境 A system for data analysis and visualization which is built based on S language 當然這個解釋還是比較 難解 籠統(tǒng)的講 R是一個可以和SPSS SAS這些商用軟件對抗的統(tǒng)計軟件 我第一次接觸到R就是在 2011年暑假 南京大學和Hopkins大學聯(lián)合舉辦生物統(tǒng)計的暑期課程 在第一節(jié)課上 助教就展示了 關于R的用法和思路 雖然那是一個暑假課程但Hopkins大學的老師一樣很深入的講了R的使用 印象 最深刻的是關于用R來做基因組序列分析的實驗部分 從那以后我對R就很喜歡 甚至很多的實驗都是 借助R來完成的 R 其實是Rick Becker John Chambers 和Allan Wilks 開發(fā)的 S 語言的一種實現(xiàn) 這三個人都在 貝爾實驗室工作 他們合作開發(fā)了S語言來進行統(tǒng)計的運算 并且開發(fā)了商用的SPLUS 關于R的歷史 大家可以google 簡而言之 R是一個開放 開源 共享的軟件系統(tǒng) 借助它你可以實現(xiàn)幾乎所有的統(tǒng)計和計算工作 甚至有人用R來做模擬和分析 R有著開源軟件共有的特點 拓展包多 借助各種拓展包 你可以實現(xiàn) 從生物上的基因序列分析到商業(yè)上網站訪問的聚類分析 同時任何人都可以自由的使用 分享R而不需 要擔心版權的問題 這大大促進了R的發(fā)展 正如Hopkins的老師說的 Not all people can afford the cost of using SPSS or SAS even though they are easy to start with 當然在中國 這點并不是特別重要 但是可怕的是如果你和外國同學交流或者你有機會出國的時候卻發(fā)現(xiàn)沒法和大家 交流因為你使用的是盜版軟件而別人很多在使用開源的軟件 你卻都沒有聽說過 更可怕的是當我們環(huán) 顧圖書館的藏書時候發(fā)現(xiàn)好多都是在講怎么用SPSS和SAS的 老師上課也是這么教同學的 這樣的循 環(huán)是可怕的 因為大家永遠也接觸不到開源軟件的魅力與開源的精神 同時 或許在醫(yī)學 商業(yè)領域SPSS們可以和R相匹敵 但是在生物統(tǒng)計或者大規(guī)模數(shù)據(jù)分析上R 可以說是一騎絕塵 鮮有敵手 因為這些需要高度專門性的領域 R這樣可以編程 自由修改的軟件系 統(tǒng)絕對是適用的 2 2 安裝安裝 R R 2 1windows2 1windows 環(huán)境環(huán)境 在 Google 直接輸入 R 第一個搜索結果就是 R project 的官網 這就是我喜歡google 的地方 客觀公 正而不是所謂的用錢可以買到的鳳巢 或者去 http cran r project org mirrors html 選擇一個鏡像地址 訪問 下載Precompiled Binary Distributions 中的軟件 如果是 Windows 用戶 可以點擊Windows 95 and later 進入base 選擇 rwxxxx exe 下載 很快就可以下載好 然后安裝一般的 windows 下的安裝程序 安裝 除了語言一般選擇默認的配置就可以使用 R 了 然后可以下載 RStudio 這個軟件 它可以把 R 的界面改裝成和 Matlab 類似的界面 工作空間 內存 空間 結果和軟件包展示區(qū) 相當?shù)姆奖?2 2Linux2 2Linux 環(huán)境環(huán)境 在 ubuntu 的環(huán)境下主要是使用RKWard 這樣的類似編輯環(huán)境 當然在 Linux 下也是可以用Rstudio 的只是要去官網上下載deb 的包 我直接展示了同樣的 R 程序在 RKWard 下的默認的效果圖 Illustration 1 RStudio Illustration 2 RKWard 可以看到也可以直接顯示數(shù)據(jù)同時進行操作 當然也可以使用命令行的模式運行 R 3 R3 R 的語法的語法 3 13 1 數(shù)據(jù)結構數(shù)據(jù)結構 首先注意 R 是大小寫敏感的 其次在 R 中基本的數(shù)據(jù)結構 數(shù)值型 numeric 1 1 2 3 1415926 復數(shù)型 complex 1 2i 字符型 character A hello world 邏輯型 logical TRUE FALSE R 中基本的數(shù)據(jù)對象為下列五類 向量 vector 矩陣 matrix 列表 list 數(shù)據(jù)框 data frame 函數(shù) function 3 23 2 基本類型的操作基本類型的操作 賦值操作為 對于任何一個變量名X1 使用下列命令 X1 c 其中c 代表一個數(shù)據(jù)列 比如 X1 c 1 2 3 4 就是將數(shù)組 1 2 3 4 f 賦值給 X1 這個時候我們從RKWard 的內存去就可以看到此時 X1 變量是在全局變量里面 它的概要信息顯示在紅線 框中 最小值為1 00 平均值為 2 50 最大值 4 00 此處簡單驗證數(shù)據(jù)有效性 在數(shù)據(jù)很多時 這個功能是 很實用的 當然 如果我們從 R 的數(shù)據(jù)類型上入手的話 此時我們是生成一個向量 下面演示一個矩陣的生成 方法 這個時候很死板的方法是把矩陣一行一行的輸入 中間用 隔開即可 當然也可以使用比較高 Illustration 3 賦值操作 級的方法 比如 y1 matrix c 1 2 3 4 5 6 7 8 9 nrow 3 ncol 3 byrow FALSE y2 matrix 1 9 nrow 3 ncol 3 byrow FALSE 可以看到兩種方法產生的效果完全一樣 下面介紹列表 列表最大的特點就是可以實現(xiàn)不同類型的數(shù)據(jù)可以用同一個變量表示 從而構成不同數(shù) 據(jù)類型的一個變量集合 比如 x list gender c F M grade c 98 100 90 undergrad TRUE 上面的命令是建立一個list 的列表 可以看到各個對象的屬性并不相同 而且長度也不強行要求一致 你還可以用下列三個命令輸入對于gender 進行操作 如 x gender 和 x 1 最后的一個是要求R 告訴 x 這個變量名有多少變量 如 names x Illustration 4 矩陣操作 Illustration 5 列表操作 下面介紹數(shù)據(jù)框 數(shù)據(jù)框是和列表相同的一類數(shù)據(jù)結構 都可以處理不同類型的數(shù)據(jù) 但是數(shù)據(jù)框嚴格 要求各個組成部分的長度相同 比如 y data frame gender c F M grade c 98 100 undergrad c FALSE TRUE 可以看到上面的操作 我特意選擇了和list 相同的內容 R 直接報錯 而在R 的內存中可以看到正確的 y的 數(shù)據(jù)分布 希望大家可以知道各種不同的數(shù)據(jù)結構在內存中具體的邏輯表示 因為不同的數(shù)據(jù)結構直接決定你 在設計不同的數(shù)據(jù)處理時候的方法 3 33 3 數(shù)據(jù)文件的讀寫操作數(shù)據(jù)文件的讀寫操作 既然是對于統(tǒng)計數(shù)據(jù)進行操作 所以肯定離不開具體的對于數(shù)據(jù)文件的操作 一般而言最最方便的 數(shù)據(jù)讀寫形式就是文本文件 當然很多的時候我們使用的 excel 的格式 這個時候建議大家把 excel 文件 存成cvs 的格式 這樣也可以方便進行操作 當然如果你下載了讀取excel 文件的專門包的話 也是可以 操作的 主要的讀寫操作命令是下面兩個 read table 適合 讀寫 ASCII 文件 read csv 適合讀寫 Excel CSV 文件 要注意的是read table 并不適合讀寫太多太大的矩陣結構 這個時候可以用scan 命令代替 比如我有一個已經存儲好的txt 數(shù)據(jù)文件 名字叫做lorenz txt 下面就有這個文件做一個數(shù)據(jù)讀寫 的展示 則命令是 x read table lorenz txt header FALSE 注意 這里我省略了設置你的工作目錄的工作 這是因為 我把這個文件放在我的 home 文件夾下 如果不是 Illustration 6 數(shù)據(jù)框的操作 的話 運行下列命令 getwd 得到你的工作目錄 一般在Linux 中這個是默認設置在 home 文件夾下的 當然如果你在windows 下可以用如下命令 setwd C Users fame 里面填上你喜歡的路徑即可 這個文件是關于洛淪茲方程的 所以我加了一個 plot 命令展示它的結果 當然如果是 excel 的格式使用下面的命令 y read csv lorenz csv header TRUE 當然 另外一個命令的方向就是從R 向文件中寫數(shù)據(jù)可以使用如下命令 write table x lorenz txt col names FALSE sep t 3 43 4 一個使用一個使用 R R 的例子的例子 下面結合一個實際實驗中例子給大家展示一下用 R 可以做的物理實驗數(shù)據(jù)處理 曾近在一次實驗中 我們需要擬合一條直線 這是在單擺的振動中 我們測量它的不同振幅從而希望可以得到單擺的衰減系 數(shù) 這個實驗設計中包含如何從視頻文件中讀出數(shù)據(jù) 然后進行數(shù)據(jù)處理 這里面包含像素坐標和實際 坐標的轉換 我們假設全部得到這些數(shù)據(jù)了 我下一步需要關注數(shù)據(jù)的擬合操作 首先我們輸入數(shù)據(jù)文件 命令如下 c1 data frame X c 197 224 251 278 306 335 364 396 424 453 Y c 424 381 358 338 327 315 309 300 296 291 這里使用數(shù)據(jù)框從而保證下一步不需要擔心數(shù)據(jù)輸入的時候有缺失或者不全 下面讓我們簡單估計一下數(shù)據(jù)的分布與圖像 從而得到下一步的行動方法 加上一個 plot 命令即可 Illustration 7 lorenz plot c1 圖像如下 下面的任務就是對于這個圖像的擬合 從原理上分析我們已經知道了這個圖像的衰減規(guī)律符合如下的函數(shù)形式 Y a b exp c X 所以所 有的工作就是求解這個擬合的系數(shù) 我們使用如下的命令 nls sol nls Y a b exp c X data c1 start list b 215 8 a 40 18 c 0 01478 nls sum t b 1 347e 03 1 829e 02 7 366 0 000154 a 2 870e 02 2 734e 00 104 965 1 88e 12 c 1 170e 02 7 230e 04 16 182 8 37e 07 Signif codes 0 0 001 0 01 0 05 0 1 1 Residual standard error 2 62 on 7 degrees of freedom Number of iterations to convergence 8 Achieved convergence tolerance 1 686e 06 如果你還是沒有辦法從上面的殘差分析中看出擬合的好壞的話 我們使用如下命令擬合所有的曲線 得 到具體的擬合的圖像結果 xfit seq 180 480 len 2000 Illustration 8 初步plot yfit predict nls sol data frame X xfit plot c1 X c1 Y lines xfit yfit 結果如下 可以看到擬合的結果相當好 最后的結果是 Y 287 0 1347exp 0 0117X 3 53 5 進一步深入提高進一步深入提高 上面的內容充其量只能算是一個導引 大家如果愿意深入的學習 R 的操作有很多的材料可以學習 第一個就是An Introduction to R by Venables and Smith 這個也是 R 的官方簡介 也是全世界學習 R 的入手手冊 幾乎可以說初級中級使用者所有的東西這本書都具備 當然 它的缺點就是例子太過簡單 同時因為簡單所以很多例子缺少實戰(zhàn)效果 推薦薛毅的 統(tǒng)計建模與 R 軟件 這個是相當不錯的書 前 面復習了一下統(tǒng)計的基本知識后面幾乎都是結合例子進行的講解 當然如果你是為了快速查找某些命令 可以看看 R 語言使用筆記 這是從薛毅的 統(tǒng)計建模與 R 軟件 中精簡出來的命令集合 當然網上也 有很多的材料 最著名的就是人大統(tǒng)計之都 COS 論壇 http cos name cn 確實是很有效果