One day, I tried to finish a project in which I used some flat file. First, I save the file to .txt tab delimited file. And when I read in the file, there was always an error, like the rows are not of the same length. Then I converted the file to .csv file. There was no more error. Later on, I found, since when it was tab delimited file, some columns were not actually delimited by tab, i.e. two or more columns are collapsed together, and became one column. So for different rows the number of columns are not equal. But when I used the csv file, no such issue, since the commas are always there, and separated the columns.
So maybe the csv files are safer to used than tab delimited file.
Saturday, November 22, 2014
Python's Levenshtein method is much faster than difflib.SquenceMatcher method
I did one project last week. The project is mainly about text mining, comparing the similarity between two strings.
And I used Python's difflib.SequenceMatcher to do the similarity comparison. It took around 120 mins to finish a 400M times computations.
And after I switched to Levenshtein method on Python, it only took 30 mins to finish the same amount of computation.
I was told Levenshtein used C to do the computation and difflib.SequenceMatcher used Python to do the computation. So Levenshtein is much faster.
And I used Python's difflib.SequenceMatcher to do the similarity comparison. It took around 120 mins to finish a 400M times computations.
And after I switched to Levenshtein method on Python, it only took 30 mins to finish the same amount of computation.
I was told Levenshtein used C to do the computation and difflib.SequenceMatcher used Python to do the computation. So Levenshtein is much faster.
Thursday, November 6, 2014
This tutorial greatly helps me understand classes of python language
This tutorial greatly helps me understand classes in python. It explains everything very clear.
http://www.sthurlow.com/python/lesson08/
http://www.sthurlow.com/python/lesson08/
Saturday, October 4, 2014
How to write a new file to a directory
import os
import os.path
import StringIO
import csv
dir = r"C:\Python27"
if not os.path.exists(dir):
os.mkdir(dir)
my_list=[[1,2,3],[4,5,6]]
with open(os.path.join(dir, "filename"+'.csv'), "w") as f:
csvfile=StringIO.StringIO()
csvwriter=csv.writer(csvfile)
for l in my_list:
csvwriter.writerow(l)
for a in csvfile.getvalue():
f.writelines(a)
the file filename.csv, it does not exist on the directory C:\Python27 before. When writing the new file to the directory, we create a file there. Wednesday, August 27, 2014
A good Naive Bayes Classifier tutorial
A good Naive Bayes Classifier tutorial
https://www.youtube.com/watch?v=j3IGd5CjsVA
I think he explained everything very clear, and he also gave a tutorial for the python code. Good one.
https://www.youtube.com/watch?v=j3IGd5CjsVA
I think he explained everything very clear, and he also gave a tutorial for the python code. Good one.
Wednesday, July 30, 2014
How I get mysql on mac
First install mysql workbench here：
http://dev.mysql.com/downloads/workbench/
then install mysql server here：
http://dev.mysql.com/downloads/mysql/
And then , on mac terminal, use command:
cd /usr/local/mysql/bin
sudo ./mysqld_safe --console
Password required
http://dev.mysql.com/downloads/workbench/
then install mysql server here：
http://dev.mysql.com/downloads/mysql/
And then , on mac terminal, use command:
cd /usr/local/mysql/bin
sudo ./mysqld_safe --console
Password required
Friday, July 25, 2014
Monday, July 21, 2014
Some useful R code
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
$x
[1] 0.75
$y
[1] 3.375
x=1/2
y=1/2+1
x=(3/2)/(3/2+1/2)
y=(3/2)*(0.75+1.5)
y
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
2 2
2 2 2
2 2
2 2 2
2 2
2 2 2
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
[1] 1 1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1 1 1 1 1
[1] 1 1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1 1 1 1 1
[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))
K <- 4
## Question 1
sum(a[a > K])
## Question 2
sum(apply(rbind(K, a), 2, max))
Or
##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
y
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:
b
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)
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))
ff
[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)
15
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)
15
Using vector operations
x=2; n=3
(sum(rep(x,(n+1))^(0:n)))
15
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)
5
Ex 2.2
(a)
a<-c(1:8,7:1)
(b) rep(1:5,1:5)
(c)
matrix(rep(1,9),ncol=3)-diag(c(1,1,1))
(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
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.
from some proposal
distribution.
3. Compute the ratio 
of the density at the
candidate and initial points, 
.
of the density at the
candidate and initial points, .
4. With probability min(
, accept the candidate point, else retain the current point.
, 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.
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, 
.
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.
, 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:
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))
or
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))
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))
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
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)")
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)
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)
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)")
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:
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
a)
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)
b)
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
a)
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)
b)
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
a)
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)
b)
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)
Subscribe to:
Posts (Atom)
-
I tried to commit script to bitbucket using sourcetree. I first cloned from bitbucket using SSH, and I got an error, "authentication ... -
https://github.com/boto/boto3/issues/134 import boto3 import botocore client = boto3.client('s3') result = client.list_obje...
-
Previously, I wanted to install "script" on Atom to run PHP. And there was some problem, like the firewall. So I tried atom-runner...