diff bin/beta_fit.R @ 1:adc0f7765d85 draft

planemo upload

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/bin/beta_fit.R	Thu Sep 07 15:06:58 2017 -0400
@@ -0,0 +1,39 @@
+args = commandArgs();
+input_file=args[4];
+#input_file="control_a.0.1.flt.ratio.tmpca";
+#print (input_file);
+
+
+library(MASS);
+
+p=array(0,dim=1000);
+
+for (i in 0:999)
+{
+	        p[i]=0.001
+}
+
+if ( file.info(input_file)["size"]>0 )
+{
+
+data=read.table(input_file);
+col=1;
+x=data[,col];
+x1=x;
+if (length(x)>10)
+{
+	x1[x==0] <- .Machine$double.eps;
+	x1[x==1] <- (1-.Machine$double.eps);
+	xbar=mean(x1)
+	xvar=var(x1)
+	a <- (xbar*(1-xbar)/xvar - 1)*xbar
+	b <- (1-xbar)*a/xbar
+	(f=fitdistr(x1,"beta",list(shape1=a,shape2=b)))
+	for (i in 0:999)
+	{
+		p[i]=dbeta(i/1000,f[["estimate"]][["shape1"]],f[["estimate"]][["shape2"]])
+	}
+}
+
+}
+write(p,file=paste(input_file,"fit",sep="."),ncolumns=1);