workingDir="/Users/hito/Dropbox/misc/r/" #ここは手元の環境にあわせて変えて下さい。
inputfilename="test.txt" #ここは手元の環境にあわせて変えて下さい。
setwd(workingDir)
fittingSin <- function(data,index,PDF){
dataNum=colnames(data)[index+1]
filename=paste(dataNum,".pdf",sep="")
print(filename)
x <- data$time
y <- data[,index+1]
x <- x[!is.na(y)]
y <- y[!is.na(y)]
xaxsmin <- ceiling(min(x)/12)
xaxsmax <- floor(max(x)/12)
xaxsvec <- xaxsmin:xaxsmax*12
result <- nls(y ~A*cos(2*pi/T*(x-a))+b, algorithm ="port", start=c(A=(max(y)-min(y))/2,a=2,T=24,b=mean(y)))
result1 <- try(nls(y ~A*cos(2*pi/T*(x-a))+b, algorithm ="port", start=c(A=(max(y)-min(y))/2,a=-6,T=24,b=mean(y))))
result2 <- try(nls(y ~A*cos(2*pi/T*(x-a))+b, algorithm ="port", start=c(A=(max(y)-min(y))/2,a=0,T=24,b=mean(y))))
result3 <- try(nls(y ~A*cos(2*pi/T*(x-a))+b, algorithm ="port", start=c(A=(max(y)-min(y))/2,a=6,T=24,b=mean(y))))
result4 <- try(nls(y ~A*cos(2*pi/T*(x-a))+b, algorithm ="port", start=c(A=(max(y)-min(y))/2,a=12,T=24,b=mean(y))))
result <- result1
if(deviance(result2)<deviance(result1)) {result <- result2}
if(deviance(result3)<deviance(result1)) {result <- result3}
if(deviance(result4)<deviance(result1)) {result <- result4}
par = t(c(coef(result),deviance(result)))
rownames(par) = colnames(data[index+1])
parnames=c(names(coef(result)),"R^2")
predictdata <- predict(result)
if(PDF){pdf(filename)}
plot(x,y,pch=16,xlim=c(min(x),max(x)),ylim=c(0,0.6),yaxt="n",xaxt="n",xlab="Time[h]",ylab="Ratio of P-KaiC",col="grey60")
axis(1,at=xaxsmin:xaxsmax*12,tck=1)
par(new=T)
plot(x,predictdata,col="red",ann=F,type="l",lwd=2,xaxt="n",xlim=c(min(x),max(x)),ylim=c(0,0.6))
sampleNo=sub("[^[:digit:]].","",colnames(data)[index+1])
mtext(paste("No.",sampleNo),cex=0.8)
textx <- (max(x)-min(x))*(-0.05) + min(x)
texty <- -(max(y)-min(y))*0.05 + 0.6
textpx <- (max(x)-min(x))*0.15 + min(x)
textpy <- -(max(y)-min(y))*0.05 + 0.6
options(digits=4)
fittext<- "y=A*cos(2*Pi/T*(t-a))+b\n"
textA <- paste("A=",round(par[1],digits=4),"\n")
texta <- paste("a=",round(par[2],digits=2),"\n")
textT <- paste("T=",round(par[3],digits=2),"\n")
textb <- paste("b=",round(par[4],digits=2),"\n")
text(textx,texty,fittext,pos=4)
text(textpx,textpy,paste(textA,texta,textT,textb),pos=1)
print("A,a,T,b")
print(c(round(par[1],digits=4),round(par[2],digits=2),round(par[3],digits=2),round(par[4],digits=2)))
if(PDF){dev.off()}
}
#get 1PDF
pdf("Fitting.pdf",family="Helvetica")
par(mfrow = c(5,4))
par(oma=c(0,0,0,0))
par(ann=FALSE)
par(mai=c(0.15,0,0.15,0))
par(mgp=c(1,0,0))
fitdata <- read.delim(inputfilename,header=T)
dataNum=length(colnames(fitdata))-1
for(i in 1:dataNum){
PDF=FALSE
fittingSin(fitdata,i,PDF)
}
dev.off()
ここ参考に