Easy Programming: July 2014

b<-rep(0,12)

xstar<-matrix(b,ncol=6)

library(MASS)

a<-c(1.1,1.2,1.3,1.4,1.5,1.6,2.1,2.2,2.3,2.4,2.5,2.6)

x<-matrix(a,ncol=6)

xm<-data.frame(x)

object<-lm.ridge(X1~.,xm,lambda=0.003)

coef4.1<-coef(object)

acoef<-abs(coef4.1)

for (i in 2:6)

{xstar[,i]<-x[,i]/acoef[i]

}

xstar

x[,1]

acoef

Iteration:

g1=function(x)

{a=x/(x+1)

return(a)

}

f1=function(g1,y,nlim)

{interno=0

repeat{

interno=interno+1

if (interno<nlim)

{y=g1(y)

}

else

{

cat("interation=",interno,"\n")

break

}

return(y)

}

f1(g1,1,3)

a<-c(1,1,1,2,2,2,3,3,3)

a[-(4:6)]

1 1 1 3 3 3

g1=function(x,y)

{a=y/(x+y)

return(a)

}

g2=function(x,y)

{b=y*(x+y)

return(b)

}

g3=function(x,y,nlim)

{iterno=0

repeat{

iterno=iterno+1

if (iterno>nlim)

{cat("iterno=",iterno,"\n")

break}

else

{x=g1(x,y)

y=g2(x,y)

}

return(list(x=x,y=y))

}

g3(1,1,2)

iterno= 3

[1] 0.75

[1] 3.375

x=1/2

y=1/2+1

x=(3/2)/(3/2+1/2)

y=(3/2)*(0.75+1.5)

3.375

a=c(1,1,1,1,1,1,1,1,1,1,1,1,1,1)

length(a)

m=matrix(a,ncol=7)

write(t(m),file="C:/Users/Soemone/Desktop/08 Fall-/Research/S Xu/Simulate Large Data/a.txt",ncol=7)

An algorithm that can normalize the columns and rows of a matrix:

m<-matrix(c(1,3,5,4,7,4,6,1,5),3,3)

rowsum<-c(1,2,3)

colsum<-c(1,2,3)

while((max(rowsum)-min(rowsum))>0.0001 && (max(colsum)-min(colsum))>0.0001)

{ for (i in 1:3)

{

rowsum<-apply(m,1,sum)

m[i,]<-m[i,]*12/rowsum[i]

}

for (i in 1:3)

{

colsum<-apply(m,2,sum)

m[,i]<-m[,i]*12/colsum[i]

}

To get the source codes from R

In a package: use methods(package)

In a function: input the name of the function, then return, if the function is in some where, then getAnywhere(function)

fun1<-function(n)

{ if(n==1)

list(success="The computation is successful!")

else

if (n<=0)

list(fail="We need a positive number!")

else while(n!=1)

{

if (n%%2==0) n<-(n/2)

else n<-(n*3+1)

print(n)

}

fun1(3)

[1] 10

[1] 5

[1] 16

[1] 8

[1] 4

[1] 2

[1] 1

fun1<-function(n)

{ if(n==1)

list(success="The computation is successful!")

else

if (n<=0)

list(fail="We need a positive number!")

else while(n!=1)

{

if (n%%2==0) n<-(n/2)

else n<-(n*3+1)

}

print(n)

}

fun1(3)

> fun1(3)

[1] 1

How to export a data frame from R into SAS

Using the package “foreign”, and the function “write.foreign”

Suppose we have a data frame named “b” in R. If the data is not data frame, we can use “as.data.frame” to convert it into a data frame. We set up an empty file name a.txt, then we write a sas code to infile the data from a.txt. At last, we use the following codes to finish the export.

library(foreign)

a<-matrix(c(1,2,3,4,5,6,7,8,9),nrow=3)

b<-as.data.frame(a)

write.foreign(b,"C:\\Users\\Soemone\\Desktop\\a.txt","C:\\Users\\Soemone\\Desktop\\a.sas",package="SAS",dataname="one")

dataname is an optional statement, using that we can convert the name of the data in SAS into the one that we want such as “one” here. May try different kind of data file such as Excel file?

In fact, we do not need to set up the data file a.txt and the sas file a.sas. Just write a path for the data file and the sas codes file, but do not need to put anything there. Then run the codes above, it works.

R can not read Excel file xls directly unless the excel file is converted into the csv file.

Draw multiple graphs in graph.

R codes:

par(mfrow=c(2,2))

a=matrix(rep(0,16),nrow=4)

for (i in 1:4)

{plot(rnorm(a[i,]))

}

Results:

Another codes:

par(mfrow=c(2,3))

a=matrix(rep(0,24),nrow=6)

for (i in 1:6)

{plot(rnorm(a[i,]))

}

par(mfrow=c(3,2))

a=matrix(rep(0,24),nrow=6)

for (i in 1:6)

{plot(rnorm(a[i,]))

}

par(mfcol=c(2,3))

a=matrix(rep(0,24),nrow=6)

for (i in 1:6)

{plot(rnorm(a[i,]))

}

par(mfcol=c(3,2))

a=matrix(rep(0,24),nrow=6)

for (i in 1:6)

{plot(rnorm(a[i,]))

}

Outer:

x<-c(1,1,1,1,1)

y<-c(1,1,1,1,1)

pf<-function(x,y){x+y}

z<-outer(x,y,FUN="pf")

results:

2 2 2 2 2

x<-c(1,1,1,1,1)

> y<-c(1,1,1,1,1)

> pf<-function(x,y){x+y}

> z<-outer(x,y,fun="pf")

> z

[,1] [,2] [,3] [,4] [,5]

[1,] 1 1 1 1 1

[2,] 1 1 1 1 1

[3,] 1 1 1 1 1

[4,] 1 1 1 1 1

[5,] 1 1 1 1 1

Transpose a data frame

Original data frame:

id year_mnth revenue

101 2009-01 1100

101 2009-02 1200

102 2009-02 1900

103 2009-03 2000

R codes:

dat<-read.table(“data”,header=T)

tapply(dat$revenue, list(dat$id, dat$year_mnth), sum)

result:

2009-01 2009-02 2009-03

101 1100 1200 NA

102 NA 1900 NA

103 NA NA 2000

Graphical images for multi-variables

Outline:

outline<-function(x,txt=TRUE)

{if (is.data.frame(x)==TRUE)

x<-ax.matrix(x)

m=nrow(x);n=ncol(x)

plot(c(1,n),c(min(x),max(x)))

for(i in 1:m)

{lines(x[i,],col=i)

if (txt==TRUE)

{k<-dimnames(x)[[1]][i]

text(1+(i-1)%%n,x[i,1+(i-1)%%n],k)

}

Star:

stars()

harmonic curve:

harmonc<-function(x)

{ if (is.data.frame(x)==TRUE)

x<-as.matrix(x)

t<-seq(-pi,pi,pi/30)

m<-nrow(x);n<-ncol(x)

f<-array(0,c(m,length(t)))

for (i in 1:m)

{f[i,]<-x[i,1]/sqrt(2)

for (j in 2:n)

{ if (j%%2==0)

f[i,]<-f[i,]+x[i,j]*sin(j/2*t)

else

f[i,]<-f[i,]+x[i,j]*cos(j%/%2*t)

}

plot(c(-pi,pi),c(min(f),max(f)))

for (i in 1:m) lines(t,f[i,],col=i)

}

data:

x1 x2 x3 x4 x5

1 96 95 98 97 100

2 91 86 83 84 85

3 92 93 94 97 100

4 90 87 88 89 95

5 82 95 94 93 97

6 78 79 76 74 73

7 75 78 76 73 72

codes:

x<-read.table(“C:\\Users\\Soemone\\Desktop\\d.txt”)

source("C:\\Users\\Soemone\\Desktop\\a.R")

harmonc(x)

graph:

copy and paste codes from word document:

When there is a “ which is originally written using word, we need to change it to ", otherwise, there will be an error:

Error: unexpected input in "x<-read.table(“"

Can give a scalar as a value to a vector:

codes:

a<-seq(-1,1,1/10)

f<-array(0,c(4,length(a)))

print(f)

for (i in 1:4)

{f[i,]<-1

print(f[i,])

}

results:

a<-seq(-1,1,1/10)

> f<-array(0,c(4,length(a)))

> print(f)

[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14]

[1,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0

[2,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0

[3,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0

[4,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0

[,15] [,16] [,17] [,18] [,19] [,20] [,21]

[1,] 0 0 0 0 0 0 0

[2,] 0 0 0 0 0 0 0

[3,] 0 0 0 0 0 0 0

[4,] 0 0 0 0 0 0 0

> for (i in 1:4)

+ {f[i,]<-1

+ print(f[i,])

+ }

[1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1

or:

a<-seq(-1,1,1/10)

f<-array(0,c(4,length(a)))

print(f)

for (i in 1:4)

{f[i,]<-sin(a)

print(f[i,])

}

results:

> a<-seq(-1,1,1/10)

> f<-array(0,c(4,length(a)))

> print(f)

[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14]

[1,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0

[2,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0

[3,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0

[4,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0

[,15] [,16] [,17] [,18] [,19] [,20] [,21]

[1,] 0 0 0 0 0 0 0

[2,] 0 0 0 0 0 0 0

[3,] 0 0 0 0 0 0 0

[4,] 0 0 0 0 0 0 0

> for (i in 1:4)

+ {f[i,]<-sin(a)

+ print(f[i,])

+ }

[1] -0.84147098 -0.78332691 -0.71735609 -0.64421769 -0.56464247 -0.47942554

[7] -0.38941834 -0.29552021 -0.19866933 -0.09983342 0.00000000 0.09983342

[13] 0.19866933 0.29552021 0.38941834 0.47942554 0.56464247 0.64421769

[19] 0.71735609 0.78332691 0.84147098

[1] -0.84147098 -0.78332691 -0.71735609 -0.64421769 -0.56464247 -0.47942554

[7] -0.38941834 -0.29552021 -0.19866933 -0.09983342 0.00000000 0.09983342

[13] 0.19866933 0.29552021 0.38941834 0.47942554 0.56464247 0.64421769

[19] 0.71735609 0.78332691 0.84147098

[1] -0.84147098 -0.78332691 -0.71735609 -0.64421769 -0.56464247 -0.47942554

[7] -0.38941834 -0.29552021 -0.19866933 -0.09983342 0.00000000 0.09983342

[13] 0.19866933 0.29552021 0.38941834 0.47942554 0.56464247 0.64421769

[19] 0.71735609 0.78332691 0.84147098

[1] -0.84147098 -0.78332691 -0.71735609 -0.64421769 -0.56464247 -0.47942554

[7] -0.38941834 -0.29552021 -0.19866933 -0.09983342 0.00000000 0.09983342

[13] 0.19866933 0.29552021 0.38941834 0.47942554 0.56464247 0.64421769

[19] 0.71735609 0.78332691 0.84147098

Gibbs sampling:

R codes from Bayesian Statistics: An introduction.

nr<-50

m<-5000

k<-10

n<-16

alpha<-2

beta<-4

lambda<-16

maxn<-24

h<-rep(0,n+1)

for (i in 1:m)

{

pi<-runif(1)

for (j in 1:k)

{

y<-rbinom(1,n,pi)

newalpha<-y+alpha

newbeta<-n-y+beta

pi<-rbeta(1,newalpha,newbeta)

}

for (t in 0:n)

{

if (t==y)

h[t+1]<-h[t+1]+1

}

print(h)

barplot(h)

results:

> print(h)

[1] 200 372 508 573 537 551 470 428 376 300 235 184 141 67 37 16 5

> barplot(h)

R codes translated from computational statistics handbook with MATLAB:

k=1000

m=500

a=2

b=4

x=rep(0,k)

y=rep(0,k)

n=16

x[1]=rbinom(1,n,0.5)

y[1]=rbeta(1,x[1]+a,n-x[1]+b)

for (i in 2:k)

{ x[i]=rbinom(1,n,y[i-1])

y[i]=rbeta(1,x[i]+a,n-x[i]+b)

}

fhat=rep(0,17)

for (j in 1:17)

{ fhat[j]=mean(dbinom(j-1,n,y[500:k]))}

barplot(fhat)

Array a_1,...a_n, want to get sum( a_j: a_j>K) and sum_{1<=j<=n} {max(a_j,K)}

a <- c(3,6,2,5,7,9)
K <- 4
## Question 1
sum(a[a > K])
## Question 2
sum(apply(rbind(K, a), 2, max))

##Question 1:

sum(t(a)%*%(a>k))

##Question 2:

sum(t(a)%*%(a>K))+sum(t(rep(K))%*%(a<K))

Multiple lines put in one plot in R:

a<-scan()

0.6566 2.1185 1.2320

0.647 2.0865 1.2325

0.6532 2.1060 1.2287

0.6487 2.1290 1.2313

0.6594 2.1285 1.2341

0.6577 2.1070 1.2343

0.6579 2.1345 1.2340

0.6734 2.1705 1.2362

0.675 2.1845 1.2372

0.6592 2.1550 1.2340

0.6647 2.1710 1.2305

a<-as.matrix(a)

dim(a)

a<-matrix(a,nrow=11)

dim(a)

colnames(a)<-c("A","B","C")

y<-1:11*5

matplot(y,a,type="l")

Notes: When you read in the multiple data, you should type a<-scan() first, then paste the data in the script, then first submit “a<-scan()”, then select one data line, then submit it. Do this line by line. At last, select a blank line, then submit it. In this way we can read in all the data.

Matplot and legend

R codes:

matplot(a,type="o",pch=c(1,2,3))

legend(1, 9, c("one","two","three"),

pch =c(1,2,3), col =1:3)

results:

Another codes:

matplot(a,type="o",pch=c(1,2,3),col=5:7)

legend(1, 9, c("one","two","three"),

pch =c(1,2,3), col =5:7)

Results:

Another codes:

matplot(a,type="o",pch=c(1,2,3),col=1:3)

legend(1, 9, c("one","two","three"),

pch =c(1,2,3), col =5:7)

Without dots in the plot.

data:

1 2 3

4 5 6

7 8 1

2 3 .

codes:

a<-read.table("C:\\Users\\Soemone\\Desktop\\b.txt")

matplot(a,col=1:3,type="l")

legend(1,7,c("one","two","three"),lty=1:3,col=1:3)

results:

Another codes:

a<-read.table("C:\\Users\\Soemone\\Desktop\\b.txt")

matplot(a,col=4:6,type="l",lty=4:6)

legend(1,7,c("one","two","three"),lty=4:6,col=4:6)

results:

When “.” is used to indicate the missing value, there will be a warning message when we draw the plot. But if “NA” is used to indicate the missing value, there will be no warning message.

Another codes:

a<-read.table("C:\\Users\\Soemone\\Desktop\\b.txt")

matplot(a,col=4:6,type="l",lty=c(1,1,1),lwd=rep(1.8,3))

legend(1,7,c("one","two","three"),lty=c(1,1,1),col=4:6,lwd=rep(1.8,3))

results:

How to calculate a factorial in R:

prod(1:n)=gamma(n+1)

An R code in the paper.

a<-scan()

643.3944 595.0646 577.5084

589.3512 574.781 566.2284

588.3362 577.3106 560.6669

588.9926 578.6509 560.6045

587.5358 578.1017 560.3376

588.4713 578.6467 560.7271

589.0297 578.5529 560.7797

588.5336 578.4009 560.7044

589.0542 579.102 560.6257

588.28 578.4213 560.3738

588.2308 579.052 560.2325

c<-matrix(a,nrow=11,byrow=T)

b<-c(0.95,0.85,0.75,0.65,0.55,0.45,0.35,0.25,0.15,0.05,0)

length(b)

matplot(b,c,type="o",pch=c(1,2,3),col=rep(1,3),xlab="Level",ylab="PE")

legend(0.1,640, c("Sample size=350","Sample size=500","Sample size=650"),

pch =c(1,2,3), col =c(1,1,1))

Program the indicator variables.

ff <- factor(sample(letters[1:5], 25, replace=TRUE))

[1] e b d e a e e a e d a c a e c e a b c e a e e e b

Levels: a b c d e

diag(nlevels(ff))[ff,]

[,1] [,2] [,3] [,4] [,5]

[1,] 0 0 0 0 1

[2,] 0 1 0 0 0

[3,] 0 0 0 1 0

[4,] 0 0 0 0 1

[5,] 1 0 0 0 0

[6,] 0 0 0 0 1

[7,] 0 0 0 0 1

[8,] 1 0 0 0 0

[9,] 0 0 0 0 1

[10,] 0 0 0 1 0

[11,] 1 0 0 0 0

[12,] 0 0 1 0 0

[13,] 1 0 0 0 0

[14,] 0 0 0 0 1

[15,] 0 0 1 0 0

[16,] 0 0 0 0 1

[17,] 1 0 0 0 0

[18,] 0 1 0 0 0

[19,] 0 0 1 0 0

[20,] 0 0 0 0 1

[21,] 1 0 0 0 0

[22,] 0 0 0 0 1

[23,] 0 0 0 0 1

[24,] 0 0 0 0 1

[25,] 0 1 0 0 0

> model.matrix(~ ff - 1)

ffa ffb ffc ffd ffe

1 0 0 0 0 1

2 0 1 0 0 0

3 0 0 0 1 0

4 0 0 0 0 1

5 1 0 0 0 0

6 0 0 0 0 1

7 0 0 0 0 1

8 1 0 0 0 0

9 0 0 0 0 1

10 0 0 0 1 0

11 1 0 0 0 0

12 0 0 1 0 0

13 1 0 0 0 0

14 0 0 0 0 1

15 0 0 1 0 0

16 0 0 0 0 1

17 1 0 0 0 0

18 0 1 0 0 0

19 0 0 1 0 0

20 0 0 0 0 1

21 1 0 0 0 0

22 0 0 0 0 1

23 0 0 0 0 1

24 0 0 0 0 1

25 0 1 0 0 0

attr(,"assign")

[1] 1 1 1 1 1

attr(,"contrasts")

attr(,"contrasts")$ff

[1] "contr.treatment"

How to find the rank of a matrix?

qr()$rank

example:

(a<-diag(c(1,1,1)))

[,1] [,2] [,3]

[1,] 1 0 0

[2,] 0 1 0

[3,] 0 0 1

qr(a)$rank

[1] 3

If you wish to know the index positions of TRUE elements of a logical vector x, then use which(x).

example:

x<-c(1,1,1,2,3,3,3,4,4,4)

which(x%%2==0)

4 8 9 10

x<-c(1,1,1,2,3,3,3,4,4,4)

a=which(x%%2==0)

(x[a])

[1] 2 4 4 4

Textbook: Introduction to scientific programming and simulation with R:

Ex 3.1

x.values <- seq(-2, 2, by = 0.1)

n <- length(x.values)

y.values <- rep(0, n)

for (i in 1:n) {

x <- x.values[i]

if (x<0||x==0)

{y=-x^3}

else if ((0<x && x<1)||x==1)

{y=x^2}

else {y=sqrt(x)}

y.values[i] <- y

}

plot(x.values,y.values,type="l")

Ex 3.2

af<-function(x,n)

{h=0

for (i in 0:n)

{h=h+x^i}

show(h)

}

af(2,3)

Ex 3.4

af<-function(x,n)

{h=0

i=0

while (i<=n)

{h=h+x^i

i=i+1}

show(h)

}

af(2,3)

Using vector operations

x=2; n=3

(sum(rep(x,(n+1))^(0:n)))

Ex. 3.5

a<-c(1,1)

b<-c(cos(1.5),-sin(1.5),sin(1.5),cos(1.5))

theta<-matrix(b,ncol=2,byrow=T)

(c<-theta%*%a)

[,1]

[1,] -0.9267578

[2,] 1.0682322

Ex 3.6

a<-c(1.1,1.2,1.3,1.4,1.5,1.6,1.7)

h=1

for (i in 1:length(a))

{h=h*a[i]}

(ave=h^(1/length(a)))

[1] 1.385527

x<-c(1.1,1.2,1.3,1.4,1.5,1.6,1.7,1.8)

h=0

for ( i in 1:length(x))

{h=h+(1/x[i])

}

a=h^(-1)

print(a)

0.1766348

(b=(sum(x^(-1)))^(-1))

0.1766348

Ex 3.7

x<-c(1.1,1.2,1.3,1.4,1.5,1.6,1.7,1.8)

sum=0

for (i in 1:length(x))

{if (i%%3==0)

{sum=sum+x[i]}

}

print(sum)

2.9

Ex 3.10

x<-c(1.1,0.7,0.4,0.3,1.9,1.8,1.7,2.1)

x.min=x[1]

for (i in 2:length(x))

{if (x[i]<x.min||x[i]==x.min)

{x.min=x[i]}

}

print(x.min)

0.3

Ex 3.13

t<-seq(0,10,0.001)

r<-sqrt(t)

theta<-2*pi*t

xt=r*cos(theta)

yt=r*sin(theta)

plot(xt,yt,type="l")

Ex 2.1.(d)

iff<-function(a)

{

if (a>0)

{show(floor(((a-floor(a))*100)%%10))}

else if(a<0)

{show(floor(((ceiling(a)-a)*100)%%10))}

else {show(0)}

}

a<--2.358

iff(a)

Ex 2.2

(a) a<-c(1:8,7:1)

(b) rep(1:5,1:5)

(d) matrix(c(0,2,3,0,5,0,7,0,0),ncol=3,byrow=T)

Ex.2.3

polar<-function(a)

{b<-rep(0,2)

r=sqrt(a[1]^2+a[2]^2)

if (a[1]==0||a[2]==0)

{theta=0}

else if (x>0||x==0)

{theta=asin(a[2]/r)}

else if (x<0)

{theta=-asin(a[2]/r)+pi}

b[1]=r;b[2]=theta;show(b)

}

a<-c(1,2)

polar(a)

2.236068 1.107149

Ex 2.4

a<-1:100

a[a%%2!=0]

[1] 1 3 5 7 9 11 13 15 17 19 21 23 25 27 29 31 33 35 37 39 41 43 45 47 49

[26] 51 53 55 57 59 61 63 65 67 69 71 73 75 77 79 81 83 85 87 89 91 93 95 97 99

a<-1:100

a[a%%2!=0|a%%3!=0]

[1] 1 2 3 4 5 7 8 9 10 11 13 14 15 16 17 19 20 21 22

[20] 23 25 26 27 28 29 31 32 33 34 35 37 38 39 40 41 43 44 45

[39] 46 47 49 50 51 52 53 55 56 57 58 59 61 62 63 64 65 67 68

[58] 69 70 71 73 74 75 76 77 79 80 81 82 83 85 86 87 88 89 91

[77] 92 93 94 95 97 98 99 100

a<-1:100

a[a%%2!=0|a%%3!=0|a%%7!=0]

[1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19

[20] 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38

[39] 39 40 41 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58

[58] 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77

[77] 78 79 80 81 82 83 85 86 87 88 89 90 91 92 93 94 95 96 97

[96] 98 99 100

Ex 2.5

queue <- c("Steve", "Russell", "Alison", "Liam")

queue[queue!="Alison"]

[1] "Steve" "Russell" "Liam"

Generate a upper triangle matrix:

R codes:

GenerateW <- function(Dimension)

{

result <- matrix(0, nrow=Dimension, ncol=Dimension)

for( i in 1:Dimension)

{

result[i,i:Dimension] <- 1

}

return(result)

}

GenerateW(3)

[,1] [,2] [,3]

[1,] 1 1 1

[2,] 0 1 1

[3,] 0 0 1

How to use sink function in R:

Example:

sink("R001.txt")

x <- c(2,-6,-4,8,5,4,1,3,4,-9,0,1)

A <- matrix(x, ncol=3)

A.prima <- t(A)

A.prima

dim(A)

dim(A.prima)

sink()

Metropolis sampler:

1. Start with some initial value

2. Given this initial value, draw a candidate value

from some proposal distribution.

3. Compute the ratio

of the density at the candidate and initial points,

4. With probability min(

, accept the candidate point, else retain the current point.

5. Return to step 2.

Using R to do a Metropolis sampler:

R codes:

sig1=0.5

sig2=0.1

sig3=10

n=500

x1=rep(0,n)

x2=x1

x3=x1

x1[1]=-10

x2[1]=0

x3[1]=0

for(i in 2:n)

{y=rnorm(1,0,1)*sig1+x1[i-1]

u=runif(1,0,1)

alpha=dnorm(y,0,1)/dnorm(x1[i-1],0,1)

if (u<=alpha)

{x1[i]=y}

else

x1[i]=x1[i-1]

}

plot(x1,type="l")

Plot:

Metropolis-Hastings sampler:

1. Start with some initial value

2. Given this initial value, draw a candidate value

from some proposal distribution.

3. Let

denote the probability of

given the current value is

. Compute the ratio

of the density at the candidate and initial points,

4. With probability min(

, accept the candidate point, else retain the current point.

5. Return to step 2.

R codes which generated a Metropolis-Hastings sampler for a scaled inverse Chi- square distribution:

The MCMC codes deal with generating a Metropolis-Hastings sampler for a scaled inverse Chi-square distribution.

Recall that we can write the probability distribution for a scaled inverse Chi-square as:

Where n is the degrees of freedom, a is the scaling parameter. Here we will assume n=5,a=4. First, we need to define a function f(x,n,a) in R:

f<-function(x,n,a) x^(-n/2)*exp(-a/(2*x))

First we need to specify a vector for the chain. Let’s call it met and give it length 5000.

met<-numeric(5000)

Lets now start the chain, say at value 1. We will call the current value in the chain last,

last<-1

We use a normal Chi-square distribution with 5 degrees of freedom. Then the codes will be:

for( i in 1:5000)

{ cand<-rchisq(1,5)

alpha<-(f(cand,5,4)/f(last,5,4))*(dchisq(cand,5)/dchisq(last,5))

if (runif(1)<min(alpha,1))

last<-cand

met[i]<-last}

plot(met,type=”l”)

Plot a histogram, and control the width of the columns using “lwd”.

a<-c(0,0,1,0,0,1.5,0,0,0,1.3,0,0,2.1)

plot(a,lwd=2,type="h")

a<-c(0,0,1,0,0,1.5,0,0,0,1.3,0,0,2.1)

plot(a,lwd=5,type="h")

Using R codes to write the logistic simulation:

x<-as.matrix(read.table("C:\\Users\\Soemone\\Desktop\\Guelph\\WinBugs\\bino\\gen.csv",sep=",",header=FALSE))

beta<-read.table("C:\\Users\\Soemone\\Desktop\\Guelph\\WinBugs\\bino\\map1.csv",sep=",",header=TRUE)

bet<-beta$effect

n<-numeric(nrow(x))

r<-numeric(nrow(x))

p<-numeric(nrow(x))

for (i in 1:nrow(x))

{n[i]=rpois(1,32)

p[i]=x[i,]%*%bet

p[i]=exp(p[i])/(1+exp(p[i]))

r[i]=rbinom(1,n[i],p[i])

}

out<-matrix(c(n,r),ncol=2)

out

write(t(out),file="\\Users\\Soemone\\Desktop\\Guelph\\WinBugs\\bino\\simulationphe.csv",sep=",",ncolumns=2)

“write” can also output character strings.

Paste Greek

a<-c(1,2,3,2,3,4,5)

plot(a,main=expression(sigma==2))

a<-c(1,2,3,2,3,4,5)

plot(a,main=~sigma~"=2")

Result:

An example:

z1<-as.matrix(read.table("C:\\Users\\Soemone\\Desktop\\Guelph\\example data sets\\highauto.csv",sep=",",header=FALSE))

z2<-as.matrix(read.table("C:\\Users\\Soemone\\Desktop\\Guelph\\example data sets\\lowauto.csv",sep=",",header=FALSE))

par(mfrow=c(2,1))

plot(z2,type="h",lwd=3.7,xlab="Lag",ylab="Autocorrelation",ylim=c(0,1))

abline(h=0)

text(25,0.9,labels="(a)")

plot(z1,type="h",lwd=3.7,xlab="Lag",ylab="Autocorrelation",ylim=c(0,1))

text(25,0.9,labels="(b)")

abline(h=0)

Layout

nf<-layout(matrix(c(1,1,2,3), 2,2, byrow = TRUE))

layout.show(nf)

results:

Curve and histogram

x <- rnorm(1000);

hist(x, freq = FALSE, col = "white");

curve(dnorm, col = 2, add = TRUE)

Density estimate:

plot(density(x.norm),main="Density estimate of data")

Layout:

layout(matrix(c(1,1,2,3), 2,2, byrow = TRUE))

x<-as.matrix(read.table("C:\\Users\\Soemone\\Desktop\\08 Fall-\\Research\\S Xu\\BayesianReview\\damage\\post-sample.csv",sep=",",header=TRUE))

y=x[,2]

par(omi=c(0.02,0.02,0.02,0.02))

plot(y,type="l", main="Diagnostics for beta",ylab="beta",xlab="iteration")

par(omi=c(0.02,0.02,2,2))

acf(y,ylab="Autocorrelation",main="",lag.max=50)

par(omi=c(0.1,1.3,1.3,0.1))

plot(density(y),type="l",ylab="Posterior Density",xlab="beta",main="")

layout(matrix(c(1,1,2,3), 2,2, byrow = TRUE))

x<as.matrix(read.table("C:\\Users\\Soemone\\Desktop\\gamma.csv",sep=",",header=TRUE))

print(x)

y=x[,6]

z=x[,1]

plot(z,y,type="l", main="Diagnostics for beta",ylab="beta",xlab="iteration",axes = FALSE)

axis(1,at=c(10000,20000,30000,40000,50000),lab=c("10000","20000","30000","40000","50000"))

axis(2,at=c(1.0,1.5,2.0),lab=c("1.0","1.5","2.0"))

box()

# par(omi=c(0.02,0.02,2,2))

acf(y,ylab="Autocorrelation",main="",lag.max=50,ylim=c(-1,1))

# par(omi=c(0.1,1.3,1.3,0.1))

plot(density(y),type="l",ylab="Posterior Density",xlab="beta",main="")

hist(y,freq=FALSE,main="")

lines(density(y))

Histogram and density curve:

x<-as.matrix(read.table("C:\\Users\\Soemone\\Desktop\\08 Fall-\\Research\\S Xu\\BayesianReview\\wheat\\gamma22.csv",sep=",",header=TRUE))

y=x[,2]

z=x[,1]

hist(y,prob=TRUE,xlab="QTL effect",ylab="prob",main="",breaks=14)

lines(density(y,adjust=2),lty=1)

box()

text(-0.4,0.63,labels="(b)")

Column:

x<-as.matrix(read.table("C:\\Users\\Soemone\\Desktop\\interval.csv",sep=",",header=TRUE))

y=x[,6]

z=x[,1]

plot(z,y,type="h",lwd=3.7,xlab="",ylab="QTL effect",ylim=c(-1.0,1.5))

text(0,1.5,labels="(a)")

abline(h=0)

abline(v=c(73.51310938,109.4142975,174.3300143,192.46167),lty=3)

require(stats) # normally loaded

x<-as.matrix(read.table("C:\\Users\\Soemone\\Desktop\\figure3-2.csv",sep=",",header=FALSE))

x1=x[,1]

x2=x[,2]

x3=x[,3]

x4=x[,4]

x5=x[,5]

x6=x[,6]

plot(x6,x5,type="h",lwd=3.7,xlab="Genome Position(centiMorgan)",ylab="QTL effect",ylim=c(-3,6.7),xlim=c(0,2400),axes = FALSE)

axis(1,at=c(0,400,800,1200,1600,2000,2400),lab=c("0","400","800","1200","1600","2000","2400"))

axis(2,at=c(-2,0,2,4,6),lab=c("-2","0","2","4","6"))

box()

text(0,6.5,labels="(b)")

abline(h=0)

lines(x6,x1,type = "l",lty=1)

lines(x6,x2,type="l",lty=1)

lines(x6,x3,type="l",lty=2)

lines(x6,x4,type="l",lty=2)

x<-as.matrix(read.table("C:\\Users\\Soemone\\Desktop\\figure1-1.csv",sep=",",header=T))

x1=x[,1]

x2=x[,2]

plot(x2,x1,type="h",lwd=3.7,xlab="Genome Position(centiMorgan)",ylab="QTL effect",ylim=c(-3,6.7),xlim=c(0,2400),axes = FALSE)

axis(1,at=c(0,400,800,1200,1600,2000,2400),lab=c("0","400","800","1200","1600","2000","2400"))

axis(2,at=c(-2,0,2,4,6),lab=c("-2","0","2","4","6"))

box()

text(50,4.66,labels="100%")

text(145,3.36,labels="94%")

text(205,-2.44,labels="73%")

text(290,-1.78,labels="57%")

text(250,2.44,labels="100%")

text(360,3.36,labels="100%")

text(610,1.3,labels="12%")

text(630,-1.3,labels="15%")

text(780,0.97,labels="12%")

text(820,1.93,labels="99%")

text(905,4.01,labels="100%")

text(990,2.44,labels="83%")

text(1100,-1.5,labels="40%")

text(1220,-1.2,labels="34%")

text(1305,-2.44,labels="55%")

text(1335,1.78,labels="40%")

text(1425,1.2,labels="36%")

text(1405,-1.93,labels="35%")

text(1800,0.91,labels="12%")

text(2300,1.09,labels="26%")

abline(h=0)

text(0,6.7,labels="(b)")

x<-as.matrix(read.table("C:\\Users\\Soemone\\Desktop\\figure7.csv",sep=",",header=T))

x1=x[,1]

x2=x[,2]

plot(x1,x2,type="o",pch=2,col=1,xlab="Type I error",ylab="PE")

x<-as.matrix(read.table("C:\\Users\\Soemone\\Desktop\\figure10-1.csv",sep=",",header=T))

x1=x[,1]

x2=x[,2]

x3=x[,3]

x4=x[,4]

x5=x[,5]

x6=x[,6]

plot(x6,x1,type="h",lwd=3.7,xlab="Genome Position(centiMorgan)",ylab="QTL effect",ylim=c(-5,5),xlim=c(0,350),axes = FALSE)

axis(1,at=c(0,50,100,150,200,250,300,350),lab=c("0","50","100","150","200","250","300","350"))

axis(2,at=c(-5,-4,-3,-2,-1,0,1,2,3,4,5),lab=c("-5","-4","-3","-2","-1","0","1","2","3","4","5"))

box()

text(0,5,labels="(a)")

abline(h=0)

abline(v=c(73.6,109.5,174.4,192.5),lty=2)

lines(x6,x2,type = "l",lty=2)

lines(x6,x3,type="l",lty=2)

lines(x6,x4,type="l",lty=1)

lines(x6,x5,type="l",lty=1)

Figure 1

win.graph(width=8,height=5)

x<-as.matrix(read.table("C:\\Users\\Soemone\\Desktop\\figure1-1.csv",sep=",",header=T))

x1=x[,1]

x2=x[,2]

plot(x2,x1,type="h",lwd=3.7,xlab="",ylab="QTL effect",ylim=c(-3,6.7),xlim=c(0,2400),axes = FALSE)

axis(1,at=c(0,400,800,1200,1600,2000,2400),lab=c("0","400","800","1200","1600","2000","2400"))

axis(2,at=c(-2,0,2,4,6),lab=c("-2","0","2","4","6"))

box()

text(0,6.5,labels="(a)")

abline(h=0)

Figure 2

win.graph(width=8,height=5)

x<-as.matrix(read.table("C:\\Users\\Soemone\\Desktop\\figure2.csv",sep=",",header=T))

hist(x,prob=TRUE,xlab="QTL effect",ylab="Density",main="",breaks=37)

lines(density(x),lty=1)

box()

Figure 3

win.graph(width=8,height=5)

x<-as.matrix(read.table("C:\\Users\\Soemone\\Desktop\\figure3-2.csv",sep=",",header=F))

x1=x[,1]

x2=x[,2]

x3=x[,3]

x4=x[,4]

x5=x[,5]

x6=x[,6]

plot(x6,x5,type="h",lwd=3.7,xlab="Genome Position(centiMorgan)",ylab="QTL effect",ylim=c(-3,6.7),xlim=c(0,2400),axes = FALSE)

axis(1,at=c(0,400,800,1200,1600,2000,2400),lab=c("0","400","800","1200","1600","2000","2400"))

axis(2,at=c(-2,0,2,4,6),lab=c("-2","0","2","4","6"))

box()

text(0,6.5,labels="(b)")

abline(h=0)

lines(x6,x1,type = "l",lty=1)

lines(x6,x2,type="l",lty=1)

lines(x6,x3,type="l",lty=2)

lines(x6,x4,type="l",lty=2)

Figure 4

win.graph(width=8,height=5)

x<-as.matrix(read.table("C:\\Users\\Soemone\\Desktop\\figure4-3.csv",sep=",",header=F))

x1=x[,1]

x2=x[,2]

x3=x[,3]

x4=x[,4]

x5=x[,5]

x6=x[,6]

plot(x6,x5,type="h",lwd=3.7,xlab="Genome Position(centiMorgan)",ylab="QTL effect",ylim=c(-3,6.7),xlim=c(0,2400),axes = FALSE)

axis(1,at=c(0,400,800,1200,1600,2000,2400),lab=c("0","400","800","1200","1600","2000","2400"))

axis(2,at=c(-2,0,2,4,6),lab=c("-2","0","2","4","6"))

box()

text(0,6.5,labels="(c)")

abline(h=0)

lines(x6,x1,type = "l",lty=1)

lines(x6,x2,type="l",lty=1)

lines(x6,x3,type="l",lty=2)

lines(x6,x4,type="l",lty=2)

Figure 5 (a)

win.graph(width=8,height=5)

x<-as.matrix(read.table("C:\\Users\\Soemone\\Desktop\\figure1-1.csv",sep=",",header=T))

x1=x[,1]

x2=x[,2]

plot(x2,x1,type="h",lwd=3.7,xlab="Genome Position(centiMorgan)",ylab="QTL effect",ylim=c(-3,6.7),xlim=c(0,2400),axes = FALSE)

axis(1,at=c(0,400,800,1200,1600,2000,2400),lab=c("0","400","800","1200","1600","2000","2400"))

axis(2,at=c(-2,0,2,4,6),lab=c("-2","0","2","4","6"))

box()

text(50,4.66,labels="100%")

text(145,3.36,labels="84%")

text(205,-2.44,labels="48%")

text(290,-1.78,labels="28%")

text(250,2.44,labels="100%")

text(360,3.36,labels="99%")

text(610,1.3,labels="1%")

text(630,-1.3,labels="5%")

text(780,0.97,labels="3%")

text(820,1.93,labels="98%")

text(905,4.01,labels="99%")

text(990,2.44,labels="63%")

text(1100,-1.5,labels="17%")

text(1220,-1.2,labels="12%")

text(1305,-2.44,labels="28%")

text(1335,1.78,labels="18%")

text(1425,1.2,labels="15%")

text(1405,-1.93,labels="17%")

text(1800,0.91,labels="8%")

text(2300,1.09,labels="12%")

abline(h=0)

text(0,6.7,labels="(a)")

Figure 5 (b)

x<-as.matrix(read.table("C:\\Users\\Soemone\\Desktop\\figure1-1.csv",sep=",",header=T))

x1=x[,1]

x2=x[,2]

plot(x2,x1,type="h",lwd=3.7,xlab="Genome Position(centiMorgan)",ylab="QTL effect",ylim=c(-3,6.7),xlim=c(0,2400),axes = FALSE)

axis(1,at=c(0,400,800,1200,1600,2000,2400),lab=c("0","400","800","1200","1600","2000","2400"))

axis(2,at=c(-2,0,2,4,6),lab=c("-2","0","2","4","6"))

box()

text(50,4.66,labels="100%")

text(145,3.36,labels="94%")

text(205,-2.44,labels="73%")

text(290,-1.78,labels="57%")

text(250,2.44,labels="100%")

text(360,3.36,labels="100%")

text(610,1.3,labels="12%")

text(630,-1.3,labels="15%")

text(780,0.97,labels="12%")

text(820,1.93,labels="99%")

text(905,4.01,labels="100%")

text(990,2.44,labels="83%")

text(1100,-1.5,labels="40%")

text(1220,-1.2,labels="34%")

text(1305,-2.44,labels="55%")

text(1335,1.78,labels="40%")

text(1425,1.2,labels="36%")

text(1405,-1.93,labels="35%")

text(1800,0.91,labels="12%")

text(2300,1.09,labels="26%")

abline(h=0)

text(0,6.7,labels="(b)")

Figure 6

win.graph(width=8,height=5)

x<-as.matrix(read.table("C:\\Users\\Soemone\\Desktop\\figure6-2.csv",sep=",",header=T))

x1=x[,1]

x2=x[,2]

plot(x2,x1,type="h",lwd=3.7,xlab="Genome Position(centiMorgan)",ylab="False Positive Rate",ylim=c(0,1),xlim=c(0,2400),axes = FALSE)

axis(1,at=c(0,400,800,1200,1600,2000,2400),lab=c("0","400","800","1200","1600","2000","2400"))

axis(2,at=c(0,0.2,0.4,0.6,0.8,1.0),lab=c("0.0","0.2","0.4","0.6","0.8","1.0"))

box()

abline(h=0)

abline(h=0.1,lty=2)

text(0,1,labels="(b)")

Figure 7

win.graph(width=8,height=5)

x<-as.matrix(read.table("C:\\Users\\Soemone\\Desktop\\figure7.csv",sep=",",header=T))

x1=x[,1]

x2=x[,2]

plot(x1,x2,type="o",pch=2,col=1,xlab="Type I error",ylab="PE",xlim=c(0,1))

Figure 8 (a)

win.graph(width=8,height=5)

x<-as.matrix(read.table("C:\\Users\\Soemone\\Desktop\\figure8-1.csv",sep=",",header=T))

x1=x[,1]

x2=x[,2]

x3=x[,3]

x4=x[,4]

x5=x[,5]

x6=x[,6]

plot(x6,x5,type="h",lwd=3.7,xlab="Genome Position(centiMorgan)",ylab="QTL effect",ylim=c(-6,4),xlim=c(0,360),axes = FALSE)

axis(1,at=c(0,40,80,120,160,200,240,280,320,360),lab=c("0","40","80","120","160","200","240","280","320","360"))

axis(2,at=c(-6,-4,-2,0,2,4),lab=c("-6","-4","-2","0","2","4"))

box()

text(0,4,labels="(a)")

abline(h=0)

abline(v=c(85.6,146.6,212.3,282.2),lty=2)

lines(x6,x1,type = "l",lty=2)

lines(x6,x2,type="l",lty=2)

lines(x6,x3,type="l",lty=1)

lines(x6,x4,type="l",lty=1)

Figure 8 (b)

win.graph(width=8,height=5)

x<-as.matrix(read.table("C:\\Users\\Soemone\\Desktop\\figure8-2.csv",sep=",",header=T))

x1=x[,1]

x2=x[,2]

plot(x1,x2,type="o",pch=2,col=1,xlab="Type I error",ylab="PE",xlim=c(0,1),ylim=c(30,55))

text(0,55,labels="(b)")

Figure 9 (a)

win.graph(width=8,height=5)

x<-as.matrix(read.table("C:\\Users\\Soemone\\Desktop\\figure9-1.csv",sep=",",header=T))

x1=x[,1]

x2=x[,2]

x3=x[,3]

x4=x[,4]

x5=x[,5]

x6=x[,6]

plot(x6,x1,type="h",lwd=3.7,xlab="Genome Position(centiMorgan)",ylab="QTL effect",ylim=c(-1,2),xlim=c(0,1100),axes = FALSE)

axis(1,at=c(0,100,200,300,400,500,600,700,800,900,1000,1100),lab=c("0","100","200","300","400","500","600","700","800","900","1000","1100"))

axis(2,at=c(-1,0,1,2),lab=c("-1","0","1","2"))

box()

text(0,2,labels="(a)")

abline(h=0)

abline(v=c(140,320.8,486.1,632.4,836.2,946.4),lty=2)

lines(x6,x2,type = "l",lty=2)

lines(x6,x3,type="l",lty=2)

lines(x6,x4,type="l",lty=1)

lines(x6,x5,type="l",lty=1)

Figure 9 (b)

win.graph(width=8,height=5)

x<-as.matrix(read.table("C:\\Users\\Soemone\\Desktop\\figure9-2.csv",sep=",",header=T))

x1=x[,1]

x2=x[,2]

plot(x1,x2,type="o",pch=2,col=1,xlab="Type I error",ylab="PE",xlim=c(0,1),ylim=c(18,21),axes = FALSE)

axis(1,at=c(0,0.2,0.4,0.6,0.8,1),lab=c("0.0","0.2","0.4","0.6","0.8","1.0"))

axis(2,at=c(18,19,20,21),lab=c("18","19","20","21"))

box()

text(0,21,labels="(b)")

Figure 10 (a)

win.graph(width=8,height=5)

x<-as.matrix(read.table("C:\\Users\\Soemone\\Desktop\\figure10-1.csv",sep=",",header=T))

x1=x[,1]

x2=x[,2]

x3=x[,3]

x4=x[,4]

x5=x[,5]

x6=x[,6]

plot(x6,x1,type="h",lwd=3.7,xlab="Genome Position(centiMorgan)",ylab="QTL effect",ylim=c(-5,5),xlim=c(0,350),axes = FALSE)

axis(1,at=c(0,50,100,150,200,250,300,350),lab=c("0","50","100","150","200","250","300","350"))

axis(2,at=c(-5,-4,-3,-2,-1,0,1,2,3,4,5),lab=c("-5","-4","-3","-2","-1","0","1","2","3","4","5"))

box()

text(0,5,labels="(a)")

abline(h=0)

abline(v=c(73.6,109.5,174.4,192.5),lty=2)

lines(x6,x2,type = "l",lty=2)

lines(x6,x3,type="l",lty=2)

lines(x6,x4,type="l",lty=1)

lines(x6,x5,type="l",lty=1)

Figure 10 (b)

win.graph(width=8,height=5)

x<-as.matrix(read.table("C:\\Users\\Soemone\\Desktop\\figure10-2.csv",sep=",",header=T))

x1=x[,1]

x2=x[,2]

plot(x1,x2,type="o",pch=2,col=1,xlab="Type I error",ylab="PE",xlim=c(0,1),ylim=c(415,455),axes = FALSE)

axis(1,at=c(0,0.2,0.4,0.6,0.8,1),lab=c("0.0","0.2","0.4","0.6","0.8","1.0"))

axis(2,at=c(415,425,435,445,455),lab=c("415","425","435","445","455"))

box()

text(0,455,labels="(b)")

library(scatterplot3d)

xx<-as.matrix(read.table("C:\\Users\\Soemone\\Desktop\\08 Fall-\\Research\\S Xu\\Compare QTL Methods\\analysis-inter\\sim-poi-blup-03.csv",sep=",",header=T))

x<-xx[,1]

y<-xx[,2]

z<-xx[,3]

scatterplot3d(y,x,z,color="blue",type="h",lwd=3,xlim=c(1,161),ylim=c(1,161))

library(lattice)

xx<-read.table("C:\\Users\\Soemone\\Desktop\\08 Fall-\\Research\\S Xu\\Compare QTL Methods\\analysis-inter\\sim-poi-blup-03.csv",sep=",",header=T)

wireframe(effect~marker1*marker2,data=xx, color="blue",type="h",lwd=3)

Compare two methods:

win.graph(width=8,height=5)

x<-as.matrix(read.table("C:\\Users\\Soemone\\Desktop\\08 Fall-\\Research\\S Xu\\Compare QTL Methods\\analysis-nointer\\simul2\\pseudo\\em-jeffreys-blup-poi-01.csv",sep=",",header=F))

x1=x[,7]

x2=x[,8]

plot(x2,x1,type="h",lwd=3,col=1,xlab="Genome Position(centiMorgan)",ylab="QTL effect",xlim=c(0,2400),ylim=c(-4,6),axes = FALSE)

axis(1,at=c(0,400,800,1200,1600,2000,2400),lab=c("0","400","800","1200","1600","2000","2400"))

axis(2,at=c(-4,-2,0,2,4,6),lab=c("-4","-2","0","2","4","6"))

box()

abline(h=0)

text(0,6,labels="(d)")

Another graph codes:

x<-as.matrix(read.table("C:\\Users\\Soemone\\Desktop\\PRESSBINO.csv",sep=",",header=TRUE))

x1=x[,1]

x2=x[,2]

x3=x[,3]

plot(x1,x2,type="o",pch=1,col=1,xlab="T Values",ylab="PRESS",ylim=c(14000,16500),xlim=c(-2,7),axes = FALSE)

axis(1,at=c(-2,-1.5,-1,-0.5,0,0.5,1,1.5,2,2.5,3,3.5,4,4.5,5,5.5,6,6.5,7),lab=c("-2","-1.5","-1","-0.5","0","0.5","1","1.5","2","2.5","3","3.5","4","4.5","5","5.5","6","6.5","7"))

axis(2,at=c(14000,14500,15000,15500,16000,16500),lab=c("14000","14500","15000","15500","16000","16500"))

lines(x1,x3,type = "o",pch=2,col=2)

box()

legend(5,16250, c("Gneneral","Diag"),

pch =c(1,2), col =1:2)

win.graph(width=8,height=5)

x<-as.matrix(read.table("C:\\Users\\Soemone\\Desktop\\sim-map-1.csv",sep=",",header=T))

x1=as.matrix(x[,3])

x2=as.matrix(x[,5])

plot(x1,x2,type="h",lwd=3.7,xlab="",ylab="QTL effect",ylim=c(-2,2),xlim=c(0,2400),axes = FALSE)

axis(1,at=c(0,400,800,1200,1600,2000,2400),lab=c("0","400","800","1200","1600","2000","2400"))

axis(2,at=c(-2,-1.5,-1,-0.5,0,0.5,1,1.5,2),lab=c("-2","-1.5","-1","-0.5","0","0.5","1.0","1.5","2.0"))

box()

text(0,2,labels="(b)")

abline(h=0)

Time and Date:

x <- as.Date("31DEC2000", format = "%d%b%Y")

y <- as.Date("31DEC1999", format = "%d%b%Y")

diff <- x - y

new.dt <- y + as.difftime(diff, unit="days")

Clean R space

rm(list = ls(all = TRUE))

win.graph(width=8,height=5)

x<-as.matrix(read.table("C:\\Users\\Soemone\\Desktop\\bina.csv",sep=",",header=T))

x1=x[,1]

x2=x[,2]

x3=x[,3]

plot(x1,x2,type="h",lwd=8.0,lty=1,col="grey",xlab="Genome Position(centiMorgan)",ylab="QTL effect",ylim=c(-1.5,3),xlim=c(0,2400),axes = FALSE)

axis(1,at=c(0,400,800,1200,1600,2000,2400),lab=c("0","400","800","1200","1600","2000","2400"))

axis(2,at=c(-1.5,0,1.5,3),lab=c("-1.5","0","1.5","3.0"))

box()

lines(x1,x3,type = "h",lwd=3.7,lty=1,col="black")

Figure 1

win.graph(width=8,height=5)

x<-as.matrix(read.table("C:\\Users\\Soemone\\Desktop\\Figures\\Figure1.csv",sep=",",header=T))

x1=x[,1]

x2=x[,2]

x3=x[,3]

x4=x[,4]

x5=x[,5]

x6=x[,6]

x7=x[,7]

x8=x[,8]

x9=x[,9]

plot(x4,x5,type="h",lwd=3.7,xlab="Genome Position(centiMorgan)",ylab="QTL effect",ylim=c(-2.0,2.0),xlim=c(0,2400),axes = FALSE)

axis(1,at=c(0,200,400,600,800,1000,1200,1400,1600,1800,2000,2200,2400),lab=c("0","200","400","600","800","1000","1200","1400","1600","1800","2000","2200","2400"))

axis(2,at=c(-2.0,-1.5,-1.0,-0.5,0.0,0.5,1.0,1.5,2.0),lab=c("-2.0","-1.5","-1.0","-0.5","0.0","0.5","1.0","1.5","2.0"))

box()

abline(h=0)

abline(h=-2.0)

x10<-(x5!=0)

x11<-rep(-1.8,length(x4))

x12<-x10*x11

x13<-rep(-2.0,length(x4))

x14<-x10*x13

segments(x4,x14,x4,x12)

legend(1000,2,c("True"),lty=1,col=1,lwd=3.7)

Figure 2

win.graph(width=8,height=5)

x<-as.matrix(read.table("C:\\Users\\Soemone\\Desktop\\Figures\\Figure2.csv",sep=",",header=T))

x1=x[,1]

x2=x[,2]

x3=x[,3]

x4=x[,4]

x5=x[,5]

x6=x[,6]

x7=x[,7]

x8=x[,8]

x9=x[,9]

x10=x[,10]

plot(x4,x8,type="h",lwd=3.7,col=6,xlab="Genome Position(centiMorgan)",ylab="QTL effect",ylim=c(-2.0,2.0),xlim=c(0,2400),axes = FALSE)

axis(1,at=c(0,200,400,600,800,1000,1200,1400,1600,1800,2000,2200,2400),lab=c("0","200","400","600","800","1000","1200","1400","1600","1800","2000","2200","2400"))

axis(2,at=c(-2.0,-1.5,-1.0,-0.5,0.0,0.5,1.0,1.5,2.0),lab=c("-2.0","-1.5","-1.0","-0.5","0.0","0.5","1.0","1.5","2.0"))

box()

abline(h=0)

abline(h=-2.0)

x11<-(x5!=0)

x12<-rep(-1.8,length(x4))

x13<-x11*x12

x14<-rep(-2.0,length(x4))

x15<-x11*x14

segments(x4,x15,x4,x13)

legend(1000,2,c("Overdispersion"),lty=1,col=6,lwd=3.7)

Figure 3

win.graph(width=8,height=5)

x<-as.matrix(read.table("C:\\Users\\Soemone\\Desktop\\Figures\\Figure3.csv",sep=",",header=T))

x1=x[,1]

x2=x[,2]

x3=x[,3]

x4=x[,4]

x5=x[,5]

x6=x[,6]

x7=x[,7]

x8=x[,8]

x9=x[,9]

x10=x[,10]

plot(x5,x6,type="h",lwd=3.7,col=4,xlab="Marker location(cM)",ylab="QTL effect",ylim=c(-2.0,2.0),xlim=c(0,360),axes = FALSE)

axis(1,at=c(0,40,80,120,160,200,240,280,320,360),lab=c("0","40","80","120","160","200","240","280","320","360"))

axis(2,at=c(-2.0,-1.5,-1.0,-0.5,0.0,0.5,1.0,1.5,2.0),lab=c("-2.0","-1.5","-1.0","-0.5","0.0","0.5","1.0","1.5","2.0"))

box()

abline(h=-2)

abline(h=0)

abline(v=c(75.8845,114.1891,181.5191,201.6754),lty=2)

x11<-(x3!=0)

x12<-rep(-1.8,length(x3))

x13<-x11*x12

x14<-rep(-2.0,length(x3))

x15<-x11*x14

lines(x5,x8,type="h",lwd=3.7,col=2)

segments(x5,x15,x5,x13)

legend(270,2,c("Expectation","Overdispersion"),lty=1,col=c(4,2),lwd=3.7)

win.graph(width=8,height=5)

x<-as.matrix(read.table("C:\\Users\\Soemone\\Desktop\\Figures\\Figure3.csv",sep=",",header=T))

x1=x[,1]

x2=x[,2]

x3=x[,3]

x4=x[,4]

x5=x[,5]

x6=x[,6]

x7=x[,7]

x8=x[,8]

x9=x[,9]

x10=x[,10]

plot(x5,x7,type="h",lwd=3.7,col=4,xlab="Marker location(cM)",ylab="LOD score",ylim=c(-20,600),xlim=c(0,360),axes = FALSE)

axis(1,at=c(0,40,80,120,160,200,240,280,320,360),lab=c("0","40","80","120","160","200","240","280","320","360"))

axis(2,at=c(0,100,200,300,400,500,600),lab=c("0","100","200","300","400","500","600"))

box()

abline(h=0)

abline(h=-20)

abline(v=c(75.8845,114.1891,181.5191,201.6754),lty=2)

x11<-(x3!=0)

x12<-rep(-20,length(x3))

x13<-x11*x12

x14<-rep(0,length(x3))

x15<-x11*x14

lines(x5,x9,type="h",lwd=3.7,col=2)

segments(x5,x15,x5,x13)

legend(270,600,c("Expectation","Overdispersion"),lty=1,col=c(4,2),lwd=3.7)

Figure 4

win.graph(width=8,height=5)

x<-as.matrix(read.table("C:\\Users\\Soemone\\Desktop\\Figures\\Figure4.csv",sep=",",header=T))

x1=x[,1]

x2=x[,2]

x3=x[,3]

x4=x[,4]

bins=seq(-1.5,31.5,by=3)

hist(x3,breaks=bins,col="gray",main="",xlab="Number of seeded spikelets",axes = FALSE)

axis(1,at=seq(0,30,by=3),lab=c("0","3","6","9","12","15","18","21","24","27","30"))

axis(2,at=c(0,10,20,30,40,50,60,70,80,90),lab=c("0","10","20","30","40","50","60","70","80","90"))

box()

Figure 5

win.graph(width=8,height=5)

x<-as.matrix(read.table("C:\\Users\\Soemone\\Desktop\\Figures\\Figure5.csv",sep=",",header=T))

x1=x[,1]

x2=x[,2]

x3=x[,3]

x4=x[,4]

x5=x[,5]

x6=x[,6]

x7=x[,7]

x8=x[,8]

x9=x[,9]

x10=x[,10]

x11=x[,11]

a=x5-1

b=x5-2.9

plot(b,x6,type="h",lwd=3.7,col=4,xlab="Marker location(cM)",ylab="QTL effect",ylim=c(-0.2,1.4),xlim=c(0,360),axes = FALSE)

axis(1,at=c(0,40,80,120,160,200,240,280,320,360),lab=c("0","40","80","120","160","200","240","280","320","360"))

axis(2,at=c(-0.2,0.0,0.2,0.4,0.6,0.8,1.0,1.2,1.4),lab=c("-0.2","0.0","0.2","0.4","0.6","0.8","1.0","1.2","1.4"))

box()

abline(h=-0.2)

abline(h=0)

abline(v=c(75.8845,114.1891,181.5191,201.6754),lty=2)

x12<-(x3!=0)

x13<-rep(-0.1,length(x3))

x14<-x12*x13

x15<-rep(-0.2,length(x3))

x16<-x12*x15

lines(a,x8,type="h",lwd=3.7,col=2)

lines(x5,x10,type="h",lwd=3.7,col=1)

segments(x5,x16,x5,x14)

legend(270,1.4,c("Expectation","Overdispersion","Mixture"),lty=1,col=c(4,2,1),lwd=3.7)

win.graph(width=8,height=5)

x<-as.matrix(read.table("C:\\Users\\Soemone\\Desktop\\Figures\\Figure5.csv",sep=",",header=T))

x1=x[,1]

x2=x[,2]

x3=x[,3]

x4=x[,4]

x5=x[,5]

x6=x[,6]

x7=x[,7]

x8=x[,8]

x9=x[,9]

x10=x[,10]

x11=x[,11]

a=x5-1

b=x5-2.9

plot(b,x7,type="h",lwd=3.7,col=4,xlab="Marker location(cM)",ylab="LOD score",ylim=c(-2,14),xlim=c(0,360),axes = FALSE)

axis(1,at=c(0,40,80,120,160,200,240,280,320,360),lab=c("0","40","80","120","160","200","240","280","320","360"))

axis(2,at=c(0,2,4,6,8,10,12,14),lab=c("0","2","4","6","8","10","12","14"))

box()

abline(h=-2)

abline(h=0)

abline(v=c(75.8845,114.1891,181.5191,201.6754),lty=2)

x12<-(x3!=0)

x13<-rep(-1,length(x3))

x14<-x12*x13

x15<-rep(-2,length(x3))

x16<-x12*x15

lines(a,x9,type="h",lwd=3.7,col=2)

lines(x5,x11,type="h",lwd=3.7,col=1)

segments(x5,x16,x5,x14)

legend(270,14,c("Expectation","Overdispersion","Mixture"),lty=1,col=c(4,2,1),lwd=3.7)

Figure 6

win.graph(width=8,height=5)

x<-as.matrix(read.table("C:\\Users\\Soemone\\Desktop\\Figures\\Figure6.csv",sep=",",header=T))

x1=x[,1]

x2=x[,2]

x3=x[,3]

x4=x[,4]

x5=x[,5]

x6=x[,6]

x7=x[,7]

x8=x[,8]

x9=x[,9]

x10=x[,10]

x11=x[,11]

a=x5-1.2

plot(a,x6,type="h",lwd=3.7,col=4,xlab="Marker location(cM)",ylab="QTL effect",ylim=c(-0.5,0.8),xlim=c(0,360),axes = FALSE)

axis(1,at=c(0,40,80,120,160,200,240,280,320,360),lab=c("0","40","80","120","160","200","240","280","320","360"))

axis(2,at=c(-0.4,-0.2,0.0,0.2,0.4,0.6,0.8),lab=c("-0.4","-0.2","0.0","0.2","0.4","0.6","0.8"))

box()

abline(h=-0.5)

abline(h=0)

abline(v=c(75.8845,114.1891,181.5191,201.6754),lty=2)

x12<-(x3!=0)

x13<-rep(-0.45,length(x3))

x14<-x12*x13

x15<-rep(-0.5,length(x3))

x16<-x12*x15

lines(x5,x8,type="h",lwd=3.7,col=2)

segments(x5,x16,x5,x14)

legend(270,0.8,c("Expectation","Overdispersion"),lty=1,col=c(4,2,1),lwd=3.7)

win.graph(width=8,height=5)

x<-as.matrix(read.table("C:\\Users\\Soemone\\Desktop\\Figures\\Figure6.csv",sep=",",header=T))

x1=x[,1]

x2=x[,2]

x3=x[,3]

x4=x[,4]

x5=x[,5]

x6=x[,6]

x7=x[,7]

x8=x[,8]

x9=x[,9]

x10=x[,10]

x11=x[,11]

a=x5-1.3

plot(a,x7,type="h",lwd=3.7,col=4,xlab="Marker location(cM)",ylab="QTL effect",ylim=c(-2,70),xlim=c(0,360),axes = FALSE)

axis(1,at=c(0,40,80,120,160,200,240,280,320,360),lab=c("0","40","80","120","160","200","240","280","320","360"))

axis(2,at=c(0,10,20,30,40,50,60,70),lab=c("0","10","20","30","40","50","60","70"))

box()

abline(h=-2)

abline(h=0)

abline(v=c(75.8845,114.1891,181.5191,201.6754),lty=2)

x12<-(x3!=0)

x13<-rep(0,length(x3))

x14<-x12*x13

x15<-rep(-2,length(x3))

x16<-x12*x15

lines(x5,x9,type="h",lwd=3.7,col=2)

segments(x5,x16,x5,x14)

legend(270,70,c("Expectation","Overdispersion"),lty=1,col=c(4,2),lwd=3.7)

Easy Programming

ezoic

Wednesday, July 30, 2014

How I get mysql on mac

Friday, July 25, 2014

The R option to change the width and height of R plot

Monday, July 21, 2014

Some useful R code

looking for a man

Followers