[原创]Statistics 380 - Statistical Computing

是 否 +2 论坛币

k人 参与回答

经管之家送您一份

应届毕业生专属福利!

求职就业群

赵安豆老师微信：zhaoandou666

经管之家联合CDA

送您一个全额奖学金名额~ !

立即领取

感谢您参与论坛问题回答

经管之家送您两个论坛币！

+2 论坛币

07年3月开始，R的创始人Ross Ihaka 在Auckland大学统计系为3年级学生讲授Statistical Computing课程

本贴将学习Ross讲稿

课程地址：http://www.stat.auckland.ac.nz/~stat380/

该课程的软件计算包当然是R，另外还有一个代码编辑工具Tinn-R，个人觉得十分不错。

R的下载地址（要求能上外网）：

http://cran.cnr.berkeley.edu/bin/windows/base/R-2.4.1-win32.exe

Tinn-R的下载地址（要求能上外网）：

http://www.stat.auckland.ac.nz/~stat380/Tinn-R%201.17.2.4%20setup.exe

[此贴子已经被作者于2007-4-5 10:36:59编辑过]

扫码加我 拉你入群

请注明：姓名-公司-职位

以便审核进群资格，未注明则拒绝

分享0 收藏2 回帖

关键词：Statistical Statistics statistica statistic computing Statistical Statistics computing 原创

课程简介：

This course provides an introduction to statistical computing using the R statistical computing environment. The course will provide an introduction to R programming, including: 

An introduction to R 

Control flow and function writing 

Data structures and their use 

Numerical methods 

Simulation 

Object-oriented programming

Graphics programming

The content of this course differs from that of other courses in the Department of Statistics in that it does not emphasise data analysis and interpretation, but instead concentrates on programming.

We will assume that students have used R or some other similar programming environment, but will not assume that students have had extensive prior programming experience.

[此贴子已经被作者于2007-4-5 20:02:44编辑过]

Tinn-R的主界面为：

[此贴子已经被作者于2007-4-5 10:45:08编辑过]

注意为了使用方便可以对Tinn－R进行适当配置

我个人进行了两点个性化配置：

1、高亮配置。

选择option->main->application

在Active line highlighted前打勾，如图：

2、热键设置。

选择R->hotkey of R，进行如下图的设置

发送代码行：我用alt+L

发送选中代码：我用alt+P

[此贴子已经被作者于2007-4-5 10:58:45编辑过]

Ross的讲义可以到如下地址下载

http://www.stat.auckland.ac.nz/~stat380/lectures.html

共15部分，如果大家不能上外网，我会上传，下面正是开始

基本结构为：执行讲义中的程序代码，给出输出结果

同时我将上传代码Tinn-R文件，方便大家学习

[此贴子已经被作者于2007-4-5 11:59:52编辑过]

01 Introduction 105639.rar (700 Bytes) 本附件包括：

• 01 introduciton.r

2007-4-5 11:07:00 上传

> ############# 01 Introduction#################
> ### Simple R Expression
> 1+2
[1] 3
> 1/2
[1] 0.5
> 17^2
[1] 289
> 1+2*3
[1] 7
> (1+2)*3
[1] 9
> ### R functions
> sqrt(2)
[1] 1.414214
> log(10)
[1] 2.302585
> log10(10)
[1] 1
> sin(1)
[1] 0.841471
> 4*atan(1)
[1] 3.141593
> ### Assignment
> z=17
> z*(z+23)
[1] 680
> z=18
> z*(z+23)
[1] 738
> ### Numeric Vectors
> x=c(1,2,3,4)
> x
[1] 1 2 3 4
> ### Sequences
> 1:50
[1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22
[23] 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44
[45] 45 46 47 48 49 50
> 5:-5
[1] 5 4 3 2 1 0 -1 -2 -3 -4 -5
> ### Combining vectors
> x=c(1,2,3,4)
> c(x,10)
[1] 1 2 3 4 10
> c(x,x)
[1] 1 2 3 4 1 2 3 4
> ### Vector Arithmetic
> x
[1] 1 2 3 4
> 2*x+1
[1] 3 5 7 9
> sqrt(x)
[1] 1.000000 1.414214 1.732051 2.000000
> ### The Recycling Rule
> c(1,2,3,4)+c(1,2)
[1] 2 4 4 6
> ### The Modulo and Integer Division Operators
> 11%%3
[1] 2
> 1:10%%2
[1] 1 0 1 0 1 0 1 0 1 0
> 13.5%%2
[1] 1.5
> 13.5%/%2
[1] 6
> ### Elements and Subsets
> x[3]
[1] 3
> x[c(1,3)]
[1] 1 3
> x[1:3]
[1] 1 2 3
> ### Negative Subscripts
> x[c(-1,-3)]
[1] 2 4
> ### Changing Vector Subsets
> y=1:10
> y[4:6]=0
> y
[1] 1 2 3 0 0 0 7 8 9 10
> ### Special Numerical Values - Infinity
> 1/0
[1] Inf
> -1/0
[1] -Inf
> 1+Inf
[1] Inf
> 1000/Inf
[1] 0
> ### Special Numerical Values - Not a Number
> 0/0
[1] NaN
> Inf-Inf
[1] NaN
> sqrt(-1)
[1] NaN
Warning message:
NaNs produced in: sqrt(-1)
> ### Special Numerical Values - Not a Available
> 1+sin(NA)
[1] NA
> ### Summary Functions - min, max and rang
> max(1:100)
[1] 100
> max(1:100,Inf)
[1] Inf
> range(1:100)
[1] 1 100
> ### Summary Functions - sum and prod
> sum(1:100)
[1] 5050
> prod(1:10)
[1] 3628800
> ### Summary Functions and NA
> min(NA, 100)
[1] NA
> min(10,20,NA,na.rm=T)
[1] 10
> ### Cumulative Summaries
> cumsum(1:10)
[1] 1 3 6 10 15 21 28 36 45 55
> cumprod(1:10)
[1] 1 2 6 24 120 720 5040 40320
[9] 362880 3628800
> cummax(1:10)
[1] 1 2 3 4 5 6 7 8 9 10
> cummin(1:10)
[1] 1 1 1 1 1 1 1 1 1 1
> ### Parallel Summary Functions
> pmin(c(1,10), c(10,2))
[1] 1 2
> pmax(0,c(-1,0,1))
[1] 0 0 1

[此贴子已经被作者于2007-4-5 11:07:50编辑过]

02 vectors 105640.rar (788 Bytes) 本附件包括：

• 02 vectors.r

2007-4-5 11:09:00 上传

> ##################02 Vectors####################
> ###Logical Vectors
> b=c(T,T,F,F)
> b[1:3]
[1] TRUE TRUE FALSE
> b[1]=F
> b
[1] FALSE TRUE FALSE FALSE
> ###Generating Logical Values
> 3<4
[1] TRUE
> 3>4
[1] FALSE
> c(1,2,3,4)<3
[1] TRUE TRUE FALSE FALSE
> ###Logical Conjuctions
> x=c(1,2,3,4)
> x<2|x>3
[1] TRUE FALSE FALSE TRUE
> all(x>0)
[1] TRUE
> any(x>2)
[1] TRUE
> ###Logical Subsetting
> x=c(1,2,3,4)
> x[x>2]
[1] 3 4
> ###Problem
> x=1:10
> x
[1] 1 2 3 4 5 6 7 8 9 10
> x[x%%2==0]
[1] 2 4 6 8 10
> x[x%%2==1]
[1] 1 3 5 7 9
> x[c(F,T)]
[1] 2 4 6 8 10
> x[c(T,F)]
[1] 1 3 5 7 9
> ###Character Vectors
> s=c("first", "second", "third")
> s[1:2]
[1] "first" "second"
> s[1] = "initial"
> s
[1] "initial" "second" "third"
> ###String Manipulation
> s
[1] "initial" "second" "third"
> nchar(s)
[1] 7 6 5
> substring(s,1,2:3)
[1] "in" "sec" "th"
> ###Pasting Strings Together
> paste("First", "Second", "Third")
[1] "First Second Third"
> paste("First", "Second", "Third", sep=":")
[1] "First:Second:Third"
> paste("First", "Second", "Third", sep="")
[1] "FirstSecondThird"
> ###Pasting Vectors
> paste(s,"element",sep="-")
[1] "initial-element" "second-element" "third-element"
> paste(s,collapse="->")
[1] "initial->second->third"
> ###Complex-Valued Vectors
> z=10+20i
> z
[1] 10+20i
> z=-1+0i
> sqrt(z)
[1] 0+1i
> ###Vector Mode and Length
> mode(1:10)
[1] "numeric"
> length(1:10)
[1] 10
> logical(5)
[1] FALSE FALSE FALSE FALSE FALSE
> numeric(5)
[1] 0 0 0 0 0
> ###Creating Vectors
> vector("numeric",5)
[1] 0 0 0 0 0
> vector("logical",5)
[1] FALSE FALSE FALSE FALSE FALSE
> numeric(0)
numeric(0)
> ###Type Coercion
> c(T,17)
[1] 1 17
> c(T,17,"twelve")
[1] "TRUE" "17" "twelve"
> ###Type Coercion Idioms
> x=1:10
> sum(x>5)
[1] 5
> cumsum(x>5)
[1] 0 0 0 0 0 1 2 3 4 5
> ###Explicit coercion
> "1"+"2"
Error in "1" + "2" : non-numeric argument to binary operator
> as.numeric("1")
[1] 1
> ###Named Vectors
> v=1:4
> names(v)=c("A","B","C","D")
> v
A B C D
1 2 3 4
> names(v)
[1] "A" "B" "C" "D"
> ###Subsetting Using Names
> v=1:4
> names(v)=c("A","B","C","D")
> v["A"]
A
1
> v[c("A","D")]
A D
1 4
> ###List
> lst=list(10,"eleven",T)
> ###Printing Lists
> lst
[[1]]
[1] 10

[[2]]
[1] "eleven"

[[3]]
[1] TRUE

> ###Named Lists
> list(a=10,b="eleven",T)
$a
[1] 10

$b
[1] "eleven"

[[3]]
[1] TRUE

> ###Elements and Subsets of Lists
> lst=list(10,"eleven",T)
> lst[[1]]
[1] 10
> lst[1]
[[1]]
[1] 10

> ###Extracting Named Elements
> lst=list(a=10,b="eleven",T)
> lst[["a"]]
[1] 10
> lst$a
[1] 10
> ###The NULL Object
> length(NULL)
[1] 0
> as.numeric(NULL)
numeric(0)
> as.list(NULL)

list()

[此贴子已经被作者于2007-4-5 11:10:06编辑过]

03 Control Flow 105641.rar (682 Bytes) 本附件包括：

• 03 control flow.r

2007-4-5 11:11:00 上传

> ################## 03 Control Flow ###############################
> ### Expressions and Compound Expressions
> x={10;20}
> x
[1] 20
> ### Assignments within Compound Expressions
> z={x=10;y=x^2;x+y}
> x
[1] 10
> y
[1] 100
> z
[1] 110
> ### If then Else
> ### For Loop Example 1
> x=1:10
> s=0
> for(i in 1:length(x))
+ s=s+x
> s
[1] 55
> ### For Loop Example 2
> x=rnorm(10)
> x
[1] 0.09689273 -1.03691591 -1.56250044 1.46877663 0.86120905
[6] 0.25595239 -0.34517049 -0.81650957 -1.11443628 -0.27185452
> s=0
> for(i in x)
+ s=s+i
> s
[1] -2.464556
> ### The "Next" Statement
> ### While loop
> ### While loop Example
> n=0
> s=0
> while(s<=100){
+ n=n+1
+ s=s+n
+ }
> c(n,s)
[1] 14 105
> ### Repeat Loops
> ### Exiting from Loops with "Break"
> n=0
> s=0
> repeat{
+ if(s>100)
+ break
+ n=n+1
+ s=s+n
+ }
> c(n,s)
[1] 14 105
> ### Defining Functions
> ### A Simple Example of a Function
> square=function(x) x*x
> square(10)
[1] 100
> square(1:10)
[1] 1 4 9 16 25 36 49 64 81 100
> ### Functions Defined in Terms of Other Functions
> sumsq=function(x) sum(square(x))
> sumsq(1:10)
[1] 385
> sumsq=function(x) sum(square(x-mean(x)))
> sumsq(1:10)
[1] 82.5
> ### Functions in General
> ### Evaluation of Functions
> hypot=function(a,b) sqrt(a^2 + b^2)
> hypot(3,4)
[1] 5
> ### Optional Arguments to Functions
> sumsq=function(x,about=0) sum((x-about)^2)
> ### The Benefits of Default Arguments
> ### Specifying Particular Optional Arguments
> sumsq(1:10,about=mean(1:10))
[1] 82.5
> sumsq(about=mean(1:10),1:10)
[1] 82.5
> ### Argument Matching Rules
> ### Argument Matching Examples
> sumsq(1:10,mean(1:10))
[1] 82.5
> sumsq(1:10,about=mean(1:10))
[1] 82.5
> sumsq(1:10,a=mean(1:10))
[1] 82.5
> sumsq(x=1:10,mean(1:10))
[1] 82.5
> sumsq(mean(1:10),x=1:10)
[1] 82.5
>

[此贴子已经被作者于2007-4-5 11:11:41编辑过]

04 Functions 105680.rar (1.07 KB) 本附件包括：

• 04 functions.r

2007-4-5 15:13:00 上传

> ################## 04 Functions ###############################
> ### Computationsal Methods
> ### Computing Factorials by Iteration
> factorial = function(n) {
+ f = 1
+ for(i in 1:n)
+ f = f * i
+ f
+ }
> factorial(10)
[1] 3628800
> ### Computing Factorials by Recursion
> factorial =
+ function(n)
+ if (n==1) 1 else n * factorial(n-1)
> factorial(10)
[1] 3628800
> factorial(1000)
[1] Inf
> ### Computing Factorials by Vector Arithmetic
> factorial = function(n) prod(1:n)
> factorial(10)
[1] 3628800
> ### The Gamma Function
> ### Using F(x) to Compute Factorial
> factorial =
+ function(n) gamma(n+1)
> factorial(1:10)
[1] 1 2 6 24 120 720 5040 40320
[9] 362880 3628800
> ### Computing Binomial Probabilities
> bincoef =
+ function(n, k)
+ factorial(n)/(factorial(k) * factorial(n-k))
> binprob =
+ function(n, p)
+ bincoef(n, 0:n) * p^(0:n) * (1-p)^(n:0)
> round(binprob(10, .25), 4)
[1] 0.0563 0.1877 0.2816 0.2503 0.1460 0.0584 0.0162 0.0031 0.0004
[10] 0.0000 0.0000
> ### Cumulative Binomial probabilities
> bincum =
+ function(n, p)
+ cumsum(binprob(n, p))
> ### Factorials and Binomial Coefficients
> factorial = function(n) gamma(n + 1)
> choose =
+ function(n, k)
+ factorial(n)/(factorial(k)*factorial(n-k))
> ### Binomial Coefficient Example
> choose(19,3) * choose(23, 3)
[1] 1716099
> ### The Binomial Distribution
> ### Individual Binomial Probabilities
> # The chance of k successes in n trials
> # with chance of success p.
> binp =
+ function(k, n, p)
+ choose(n, k) * p^k * (1 - p)^(n - k)
> ### Binomial Probability Example
> binp(45:55, 100, 1/2)
[1] 0.04847430 0.05795840 0.06659050 0.07352701 0.07802866
[6] 0.07958924 0.07802866 0.07352701 0.06659050 0.05795840
[11] 0.04847430
> sum(binp(45:55, 100, 1/2))
[1] 0.728747
> ### Binomial Probability Example Continued
> sum(binp(40:60, 100, 1/2))
[1] 0.9647998
> sum(binp(35:65, 100, 1/2))
[1] 0.99821
> ### A Coin Tossing Experiment
> runif(10) < .5
[1] FALSE TRUE FALSE TRUE TRUE TRUE TRUE TRUE TRUE FALSE
> sum(runif(100) < .5)
[1] 53
> ### Coin Tossing Continued
> nheads = numeric(10000)
> for(i in 1:10000)
+ nheads = sum(runif(100) < .5)
> sum(45 <= nheads & nheads <= 55)/10000
[1] 0.728
> sum(40 <= nheads & nheads <= 60)/10000
[1] 0.9652
> sum(35 <= nheads & nheads <= 65)/10000
[1] 0.9979
> ### Packaging the Experiment
> cointoss =
+ function(ntosses, pheads, nreps = 1) {
+ nheads = numeric(nreps)
+ for(i in 1:nreps)
+ nheads = sum(runif(100) < pheads)
+ nheads
+ }
> ### Running the Packaged Experiment
> x = cointoss(100, .5, 10000)
> sum(45 <= x & x <= 55)/length(x)
[1] 0.7249
> ### Statistical Computations
> ### Implementing the t Test in R
> y1=rnorm(20)
> y2=rt(30,1)
> n1 = length(y1)
> m1 = mean(y1)
> s1 = sd(y1)
> n2 = length(y2)
> m2 = mean(y2)
> s2 = sd(y2)
> ### Computing the t Statistic
> sp2 = ((n1 - 1) * s1^2 + (n2 - 1) * s2^2) / (n1 + n2 - 2)
> tval = (m1 - m2) / sqrt((1/n1 + 1/n2) * sp2)
> pval = 2 * pt(- abs(tval), n1 + n2 - 2)
> ### The Complete R Function
> ttest =
+ function(y1, y2)
+ {
+ n1 = length(y1); m1 = mean(y1); s1 = sd(y1)
+ n2 = length(y2); m2 = mean(y2); s2 = sd(y2)
+ sp2 = ((n1 - 1) * s1^2 + (n2 - 1) * s2^2) / (n1 + n2 - 2)
+ tval = (m1 - m2) / sqrt((1/n1 + 1/n2) * sp2)
+ pval = 2 * pt(- abs(tval), n1 + n2 - 2)
+ list(t = tval, df = n1 + n2 - 2, pval = pval)
+ }

[此贴子已经被作者于2007-4-5 15:13:24编辑过]

05 Computation 105686.rar (1.12 KB) 本附件包括：

• 05 computation.r

2007-4-5 16:06:00 上传

> ################### 05 Computation ###########################
> ### Statistical Computations
> ### Key Quantities in the Computation
> ### Implementing the t Test in R
> diff.means =
+ function(y1, y2) mean(y1) - mean(y2)
> pooled.se =
+ function(y1, y2, n1, n2)
+ sqrt((((n1 - 1) * var(y1) +
+ (n2 - 1) * var(y2)) /
+ (n1 + n2 - 2)) * (1/n1 + 1/n2))
> ### Computing the t-Statistic
> tval = diff.means(y1, y2) /
+ pooled.se(y1, y2, length(y1), length(y2))
> ### Computing p Values
> pval = 2 * pt(- abs(tval), n1 + n2 - 2)
> ### The Complete R Function
> ttest =
+ function(y1, y2)
+ {
+ tval = diff.means(y1, y2) /
+ pooled.se(y1, y2, length(y1), length(y2))
+ pval = 2 * pt(- abs(tval), n1 + n2 - 2)
+ list(t = tval, df = n1 + n2 - 2, pval = pval)
+ }
> ### Using The t-Test Function
> y1 = rnorm(10)
> y2 = rnorm(10) + 2
> ttest(y1, y2)
$t
[1] -8.577713

$df
[1] 48

$pval
[1] 2.991950e-11

> ### Rational Approximations
> ### Continued Fractions
> ### Continued Fractions (Cont. . . )
> ### Problem Analysis
> ### Deriving The Continued Fraction
> cf.expand =
+ function(x, n = 5) {
+ cf = numeric(n)
+ for(i in 1:n) {
+ cf = round(x)
+ delta = x - cf
+ x = 1/delta
+ }
+ cf
+ }
> cf.expand(pi, 7)
[1] 3 7 16 -294 3 -4 5
> ### Producing the Rational Approximation
> rat.approx =
+ function(cf) {
+ n = length(cf)
+ num = cf[n]
+ den = 1
+ if (n > 1)
+ for(j in (n - 1):1) {
+ tmp = num
+ num = cf[j] * tmp + den
+ den = tmp
+ }
+ if (den > 0) c(num, den)
+ else c(-num, -den)
+ }
> ### Examples
> rat.approx(cf.expand(pi, 1))
[1] 3 1
> rat.approx(cf.expand(pi, 2))
[1] 22 7
> rat.approx(cf.expand(pi, 3))
[1] 355 113
> print(pi, digits = 15)
[1] 3.14159265358979
> print(22/7, digits = 15)
[1] 3.14285714285714
> print(355/113, digits = 15)
[1] 3.14159292035398
> ### Sample Sizes from Percentages
> ### Problem Analysis
> ### Algorithm
> ### The Algorithm as an R Function
> find.denom =
+ function(f, places) {
+ eps = .5 * 10^-places
+ n = 1
+ repeat {
+ i = round(n * f)
+ if (all(n * (f - eps) <= i &
+ i <= n * (f + eps)))
+ break
+ n = n + 1
+ }
+ n
+ }
> ### Results
> find.denom(c(.146, .122, .073), 3)
[1] 41
> round(41 * c(.146, .122, .073))
[1] 6 5 3
> ### Determining the Number of Digits
> 1/3
[1] 0.3333333
> round(1/3, 4)
[1] 0.3333
> x = .146
> round(x, 2) == x
[1] FALSE
> round(x, 3) == x
[1] TRUE
> eps = 1e-7
> abs(x - round(x, 2)) < eps
[1] FALSE
> abs(x - round(x, 3)) < eps
[1] TRUE
> ### An R Function
> places =
+ function(f, eps = 1e-7) {
+ places = 0
+ repeat {
+ if (all(abs(f - round(f, places)) < eps))
+ break
+ places = places + 1
+ }
+ places
+ }
> places(.146)
[1] 3
> ### An Improved find.denom Function
> find.denom =
+ function(f) {
+ eps = .5 * 10^-places(f)
+ n = 1
+ repeat {
+ i = round(n * f)
+ if (all((f - eps) * n <= i & (f + eps) * n >= i))
+ break
+ n = n + 1
+ }
+ n
+ }
> ### Improving Performance
> n = floor(1/min(f) + eps)
Error in min(f) : object "f" not found
> floor(1/.073 + .0005)
[1] 13
>
>

[此贴子已经被作者于2007-4-5 16:06:57编辑过]