|
1
|
1 args = commandArgs();
|
|
|
2 input_file=args[4];
|
|
|
3 #input_file="control_a.0.1.flt.ratio.tmpca";
|
|
|
4 #print (input_file);
|
|
|
5
|
|
|
6
|
|
|
7 library(MASS);
|
|
|
8
|
|
|
9 p=array(0,dim=1000);
|
|
|
10
|
|
|
11 for (i in 0:999)
|
|
|
12 {
|
|
|
13 p[i]=0.001
|
|
|
14 }
|
|
|
15
|
|
|
16 if ( file.info(input_file)["size"]>0 )
|
|
|
17 {
|
|
|
18
|
|
|
19 data=read.table(input_file);
|
|
|
20 col=1;
|
|
|
21 x=data[,col];
|
|
|
22 x1=x;
|
|
|
23 if (length(x)>10)
|
|
|
24 {
|
|
|
25 x1[x==0] <- .Machine$double.eps;
|
|
|
26 x1[x==1] <- (1-.Machine$double.eps);
|
|
|
27 xbar=mean(x1)
|
|
|
28 xvar=var(x1)
|
|
|
29 a <- (xbar*(1-xbar)/xvar - 1)*xbar
|
|
|
30 b <- (1-xbar)*a/xbar
|
|
|
31 (f=fitdistr(x1,"beta",list(shape1=a,shape2=b)))
|
|
|
32 for (i in 0:999)
|
|
|
33 {
|
|
|
34 p[i]=dbeta(i/1000,f[["estimate"]][["shape1"]],f[["estimate"]][["shape2"]])
|
|
|
35 }
|
|
|
36 }
|
|
|
37
|
|
|
38 }
|
|
|
39 write(p,file=paste(input_file,"fit",sep="."),ncolumns=1);
|
