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()
ここ参考に