view bin/beta_fit.R @ 5:2ebca9da5e42 draft default tip

planemo upload

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);