DNA

开发平台：
TEXT

R.txt：源码内容
							R语言 
buffon=function(n,l,d){
#win.graph()
par(mfrow=c(2,1))
Y<-runif(n,0,d/2);H<-runif(n,0,pi/2)
x<-as.integer(Y<(l*sin(H))/2)
phat<-cumsum(x)/seq(n)
pihat<-(2*l)/(d*phat)
ts.plot(phat,xlab="N",ylab="Estimated Proportion",ylim=c(0,1))
abline(h=(2*l)/(pi*d),col="red")
title("Proportion of Crossings")
ts.plot(pihat,xlab="N",ylab="Estimated Pi")
title("Running Estimate of Pi")
abline(h=pi,col="red")
out<-list(x,pihat[n],phat[n])
names(out)<-c("Outcomes","Pihat,","Phat")
out
}
BallandX<-function(a,b,c,d,n){
print("BallandX(a,b,c,d,n);a<b<c<d,n is an integer")
Allball<-c("W","W","W","B","B")
x<-NULL;x[1]<-0;NoS1<-0;NoS2<-0
for (i in 1:n){
randomBall<-sample(Allball,3)
S1<-runif(n,a,b);S2<-runif(n,c,d)
if(sum(randomBall=="W")>sum(randomBall=="B")){
x[i]<-sample(S1,1);NoS1=NoS1+1
}else {x[1]<-sample(S2,1);NoS2=NoS2+1
}#if
}#for
plot(x);list(NoS1/n,NoS2/n)
}#function
BallandX(1,2,3,4,5000)
par(mfrow=c(3,2))
x<-seq(-0.01,5,0.01);
plot(x,ppois(x,1),type="s",ylab="F(x)",main="Poisson(l)CDF")
Xrp<-rpois(100,1);plot(Xrp,main="Poisson")
plot(function(x) pnorm(x,log=TRUE),-50,10,main="log{Normal Cumulative}")
Xrn<-rnorm(1000);plot(Xrn,main="Normal")
plot(function(x) plnorm(x,meanlog=0,sdlog=1),0,100,main="longnormal CDF")
Xrln<-rlnorm(1000);plot(Xrln,main"Lognormal")
#Rejection sampling for f(x)=x^2*exp(-x)
rej.sim<-function(n){
r<-NULL
for (i in 1:n){
t<--1
while (t<0){
x<-rexp(1,1);
y<-runif(1,0,exp(-x))
if(x>1) t<--y
else t<-x^2*exp(-x) -y
#step3 accept when f(x)-y(x)>0,i.e t>0
}
r[i]<-x;
};
r;
}
par(mfrow=c(1,2));plot(rej.sim(1000));hist(rej.sim(1000))
bvsim<-function(n,me,sd,ro){
x<-rnorm(n,me[1],sd[1])
y<-rnorm(n,me[2]+ro*sd[2]/sd[1]*(x-me[1]),sd[2]*sqrt(1-ro*ro))
cbind(x,y)
}
y<-bvsim(1000,c(0,0),c(1,1),0.5)
sum(y[,1]<1&y[,2]<1)/1000
par(mfrow=c(3,1))
plot(y[,1])
plot(y[,2])
plot(y[,1],y[,2])
bootstrap<-function(x,nboot,theta,...){
z<-list()
data<-matrix(sample(x,size=length(x)*nboot,replace=T),nrow=nboot)
bd<-apply(data,1,theta,...)
est<-theta(x,...)
z$est<-est
z$distn<--bd
z$bias<-mean(bd)-est
z$se<-sqrt(var(bd))
z
}
xx<-c(1:10)
bootstrap(xx,200,mean)$est
hist(bootstrap(xx,200,mean)$distn)
law.dat<-matrix(c(576.,635.,558.,578.,666.,580.,555.,661.,651.,605.,653.,575.,545.,572.,594.,3.39,3.30,2.81,3.03,3.44,3.07,3.00,3.43,3.36,3.13,3.12,2.74,2.76,2.88,2.96),
nrow=15,ncol=2,dimnames=list(character(0),c("LSAT","GPA")))
theta1<-
grb<-read.table("Lab02RCompIntSta_GRB.R",header=T,skip=1#第一行不要读)
flux<-grb[,2#第二列];hist(flux)#做图;
n<-length(flux);xx=sort(flux)#由高到低排列;yy=(1:n)/n;plot(yy,xx);plot(xx,log(1-yy+1/n),xlab="flux",ylab="log(1-F(flux))")
mean(flux)#this is the MLE
reslm<-lm(log(1-yy+1/n)~xx)#模型
fittedvalue<--1/as.numeric(reslm$coefficients[2]#取第二个参数);fittedvalue
source("Lab02RCompIntSta_motor.R",echo=TRUE)
plot(motor.dat[,1],motor.dat[,2])
plot(motor2.dat[,1],motor2.dat[,2])
em.reg(motor2.dat,30)
library("PROcess")
fdat<-system.file("Test",package="PROcess")
fs<-list.files(fdat,full.names=TRUE)
fl<-read.files(fs[1])
plot(fl,type="l",xlab="m/z")
title(basename(fs[1]))
bseoff<-bslnoff(fl,method="loess",bw=0.1,xlab="m/z",plot=TRUE);title(basename(fs[1]))
pkgobj<-isPeak(bseoff,span=81,sm.span=11,plot=TRUE,zerothrsh=2,area.w=0.003,ratio=0.2,main="a)")
specZoom(pkgobj,xlim=c(5000,10000),main="b)")

						