genuinearray=Import[StringJoin[path,"FILES2012/a_genuinearray.dat"]];
queryonly=Import[StringJoin[path,"FILES2012/a_queryonly.dat"]];
doubleko=Import[StringJoin[path,"FILES2012/b_doublekoC.dat"]];
doublekosigma=Import[StringJoin[path,"FILES2012/b_doublekosigmaC.dat"]];

(*  *)

blocksize=0.1;   (* Divide genes in blocks according to their growth rate effect *)
binsize=0.0001;   (* For each block, bins are used to plot the chi-squared distribution *)
maxrange=0.03;   (* Chi-squared distributions are plotted from 0 to maxrange *)
xs=Table[x,{x,0.5,(maxrange/binsize),1}];
tabbies={};
plots={};
trend={};
For[i=1,i≤10,i++,
lista={};
For[genei=1,genei<Length[queryonly],genei++,
For[genej=1,genej<Length[genuinearray],genej++,
If[doubleko[[genei,genej]]<blocksize i&&doubleko[[genei,genej]]≥blocksize (i-1),
AppendTo[lista,((doublekosigma[[genei,genej]])^2)];
];
];
];
Print[i,"  ",Length[lista],"  ",Mean[lista]];
AppendTo[trend,{blocksize (i-0.5),Mean[lista],Quartiles[lista]}];
binned=BinCounts[lista,{0,maxrange,binsize}];
tabbed={xs,binned};
tabbed=Transpose[tabbed];
AppendTo[tabbies,tabbed];
AppendTo[plots,ListPlot[tabbed,Filling->Axis,PlotRange->{{0,(maxrange/binsize)},All}]];
]

(* Plot sqrt[meanvariance] (blue); 75% interval (green); median (brown); 25% interval (purple) *)

ListPlot[{Sqrt[Transpose[trend][[2]]],Sqrt[Transpose[Transpose[trend][[3]]][[1]]],Sqrt[Transpose[Transpose[trend][[3]]][[2]]],Sqrt[Transpose[Transpose[trend][[3]]][[3]]]},Joined->True]

plots

(* NOTE: tabbed is tebbed for the 13th iteration in the code above if I don't stop it earlier. *)

Transpose[tabbed][[2]];
partial=Transpose[tabbed][[2]];
cumupupu=Accumulate[Reverse[partial]];

pupu=plots[[10]]
logpupu=ListLogPlot[tabbed,Filling->Axis,PlotRange->{{0,(maxrange/binsize)},All}]
logpupu2=ListLogPlot[cumupupu,PlotRange->{{0,300},{0,500}},PlotRange->{0,1000}]

nn=4;
sigma=0.5;

all={};
For[i=1,i≤70000,i++,
nplette=(*RandomReal[NormalDistribution[0,Abs[RandomReal[NormalDistribution[0,1.5]]]],{nn}];*)
nplette=RandomReal[StudentTDistribution[0,0.666,2],{nn}];
(*nplette=RandomReal[StableDistribution[1,1.7,0,0,0.8],{nn}];*)
AppendTo[all,Variance[nplette]];
]

binsize=0.1;
binned=BinCounts[all,{0,Min[{30,Max[all]}],binsize}];
p1=ListPlot[binned,Filling->Axis,PlotRange->{{0,50},{0,6000}},PlotStyle->Green];
(*p2=Plot[(1/(sigma^2))(nn-1) binsize Total[binned]PDF[ChiSquareDistribution[nn-1],(1/(sigma^2))(nn-1) binsize (x-0.5)],{x,0,((Min[{300,Max[all]}])/binsize)},PlotStyle->{Thick,Red},PlotRange->{0,3000}];
Show[p1,p2]*)
Show[pupu,p1]
Show[p1,pupu]

(* Compare the number of experimental data points with the number of points I generate in the simulation above *)

Total[tabbies[[10]]][[2]]
Total[binned]

logp1=ListLogPlot[binned,Filling->Axis,PlotRange->{{0,300},{0,500}},PlotStyle->Green,PlotRange->{0,1000}];
Show[logpupu,logp1]

cumu=Accumulate[Reverse[binned]];
logp2=ListLogPlot[cumu,PlotRange->{{0,300},{0,100000}},PlotStyle->Green];
Show[logp2,logpupu2]

(* StableDistribution[1,1.5,0,0,0.65] *)
(* NormalDistribution[0,2.37] *)
(* CuachyDistribution[0,0.5] *)

nn=1000000;
epsa=Table[0,{nn}];
epsb=Table[0,{nn}];
For[i=1,i≤nn,i++,
mut1=RandomReal[NormalDistribution[0.6,0.1]];   (* the average noise is not 0.1 but 0.04 *)
mut2=RandomReal[NormalDistribution[0.6,0.1]];
mut12=RandomReal[NormalDistribution[0.36,0.1]];
If[mut1≤0||mut2≤0||mut12≤0,Continue[]];
epsa[[i]]=(mut12/(mut1 mut2));
epsb[[i]]=Log[2,(mut12/(mut1 mut2))];
]
pa=ListPlot[BinCounts[epsa,{0.0,2.0,0.02}],PlotRange->All]
pb=ListPlot[BinCounts[epsb,{-10,10,0.05}],PlotRange->All]
pla=ListLogPlot[BinCounts[epsa,{0.0,2.0,0.02}],PlotRange->All]
plb=ListLogPlot[BinCounts[epsb,{-10,10,0.05}],PlotRange->All]

p=LogPlot[1000000PDF[NormalDistribution[200,11],x],{x,1,500},PlotRange->{0.00001,100000}];
Show[plb,p]

(* The width of the bins above is 0.05.. below I discover that the standard deviation of the epsilons is ~0.5, so about 10 bins. In the figure above I see that I can reproduce deviations as large as 10 standard deviations. *)

Mean[epsb]
StandardDeviation[epsb]

nn=4;
npoints=Total[tabbies[[10]]][[2]];
binsize=0.0001; 
sigma=0.03;

pp=ListPlot[tabbies[[10]],Filling->Axis,PlotRange->{{0,30(*(maxrange/binsize)*)},All}];
p=Plot[(1/(sigma^2))(nn-1) binsize npoints PDF[ChiSquareDistribution[nn-1],(1/(sigma^2))(nn-1) binsize (x-0.5)],{x,0,30},PlotRange->{0,9000}]
Show[pp,p]

m=6;
n=100000;
all=Table[0,{n}];
For[i=1,i≤n,i++,
lista=RandomReal[NormalDistribution[0,1],{m}];
all[[i]]=Variance[lista];
]

Mean[all]
StandardDeviation[all]

Median[all]

p1=ListPlot[BinCounts[all,{0,3,0.01}],Filling->Axis,PlotRange->All]

nn=6;
p2=Plot[(nn/100)100000PDF[ChiSquareDistribution[nn],((nn (x-0.5))/100)],{x,0,400},PlotRange->All]

Show[p1,p2]

(4.0/(Sqrt[3]))

StandardDeviation[{0,1}]

m=0.5;
s=0.1;
sigma=(Sqrt[Log[1+((s^2)/(m^2))]])
mu=Log[(Sqrt[((m^4)/((m^2)+(s^2)))])]

p1=Plot[PDF[LogNormalDistribution[mu,sigma],x],{x,0,2},PlotRange->All]
p2=Plot[PDF[NormalDistribution[m,s],x],{x,0,2},PlotRange->All,PlotStyle->Red]
Show[p1,p2]

Mean[LogNormalDistribution[mu,sigma]]
StandardDeviation[LogNormalDistribution[mu,sigma]]

Clear[mu,sigma]

Mean[LogNormalDistribution[mu,sigma]]
StandardDeviation[LogNormalDistribution[mu,sigma]]

Mean[LaplaceDistribution[0,a]]
StandardDeviation[LaplaceDistribution[0,a]]

Mean[CauchyDistribution[mu,sigma]]
StandardDeviation[CauchyDistribution[mu,sigma]]

m=4;
n=10000;
all=Table[0,{n}];
For[i=1,i≤n,i++,
lista=RandomReal[CauchyDistribution[0,1],{m}];
all[[i]]=Variance[lista];
]

ListPlot[BinCounts[all,{0,30,0.1}],Filling->Axis,PlotRange->All]

Mean[all]
Median[all]

ListPlot[Sort[all]]

m=3;
n=1000000;
t=Table[0,{n}];
For[i=1,i≤n,i++,
lista=RandomReal[NormalDistribution[0,1],{m}];
t[[i]]=(((Sqrt[m]) Mean[lista])/(StandardDeviation[lista]));
]

p1=ListLogPlot[BinCounts[t,{-10,10,0.1}],Filling->Axis,PlotRange->All]

p2=LogPlot[1000000PDF[StudentTDistribution[100,10,2],x],{x,0,200},PlotStyle->Red]

Show[p1,p2]

m=4;
n=100000;
all=Table[0,{n}];
For[i=1,i≤n,i++,
lista=RandomReal[NormalDistribution[0,1.0],{m}];
all[[i]]=Variance[lista];
]
Mean[all]
StandardDeviation[all]

m=4;
n=100000;
all=Table[0,{n}];
For[i=1,i≤n,i++,
lista=RandomReal[StudentTDistribution[0,(1.0/(Sqrt[3])),3],{m}];
all[[i]]=Variance[lista];
]
Mean[all]
StandardDeviation[all]

m=1.0;   (* population mean *)
s=0.3;   (* population standard deviation *)
si=(Sqrt[Log[1+((s^2)/(m^2))]])   (* log-normal distribution parameter *)
mu=Log[(Sqrt[((m^4)/((m^2)+(s^2)))])]   (* log-normal distribution parameter *)

m=4;
n=100000;
all=Table[0,{n}];
For[i=1,i≤n,i++,
lista=RandomReal[LogNormalDistribution[mu,sigma],{m}];
all[[i]]=Variance[lista];
]
Mean[all]
StandardDeviation[all]

(* right def
  2.0   22,23
  1.8   33
  1.6   39
  1.4   61
  1.3   68
  1.2   83
  0.9   144
  0.8   171
  0.4   427
*)

(* wrong def
  0.5   19,20
  0.4   28,33
  0.3   86
  0.2   189
  0.1   430
*)

(* Group C *)

singlekor=Import[StringJoin[path,"FILES2012/a_queryonly.dat"]];
singlekoc=Import[StringJoin[path,"FILES2012/a_genuinearray.dat"]];
doubleko=Import[StringJoin[path,"FILES2012/b_doublekoC.dat"]];
doublekosigma=Import[StringJoin[path,"FILES2012/b_doublekosigmaC.dat"]];
eps=Import[StringJoin[path,"FILES2012/c_epsC.dat"]];

threshold=0.4;
nn=100;

stronginteraction=Table[0,{Length[singlekor]},{Length[singlekoc]}];
list={};

For[i=1,i≤Length[singlekor],i++,
For[j=1,j≤Length[singlekoc],j++,
If[doubleko[[i,j]]=="NaN",Continue[]];
If[doubleko[[i,j]]≤0,Continue[]];
gi=RandomReal[StudentTDistribution[singlekor[[i,2]],(1.0/(Sqrt[3]))((singlekor[[i,3]])/2),3],{nn}];
gj=RandomReal[StudentTDistribution[singlekoc[[j,2]],(1.0/(Sqrt[3]))((singlekoc[[j,3]])/2),3],{nn}];
gij=RandomReal[StudentTDistribution[doubleko[[i,j]],(1.0/(Sqrt[3]))((doubleko[[i,j]])/2),3],{nn}];
If[gi≤0||gj≤0||gij≤0,Print["err"];Continue[]];

gi=DeleteCases[gi,x_/;x≤0];
gj=DeleteCases[gj,x_/;x≤0];
gij=DeleteCases[gij,x_/;x≤0];
mm=Min[{Length[gi],Length[gj],Length[gij]}];
gi=Take[gi,mm];
gj=Take[gj,mm];
gij=Take[gij,mm];

(*interaction=gij-gi gj;*)
interaction=Log[2,(gij/(gi gj))];

If[eps[[i,j]]>0,   (* regardless of the def of "interaction" the sign of eps is the same *)
bad=Count[interaction,x_/;x<threshold];
If[bad<0.05 nn,
 stronginteraction[[i,j]]=1;
AppendTo[list,{singlekor[[i,2]],singlekoc[[j,2]],doubleko[[i,j]],Sign[eps[[i,j]]]}];
(*Print[singlekor[[i,2]],"  ",singlekoc[[j,2]],"  ",doubleko[[i,j]],"  ",Sign[eps[[i,j]]]];*)
];
];

If[eps[[i,j]]<0,
bad=Count[interaction,x_/;x>-threshold];
If[bad<0.05 mm,
 stronginteraction[[i,j]]=1;
AppendTo[list,{singlekor[[i,2]],singlekoc[[j,2]],doubleko[[i,j]],Sign[eps[[i,j]]]}];
(*Print[singlekor[[i,2]],"  ",singlekoc[[j,2]],"  ",doubleko[[i,j]],"  ",Sign[eps[[i,j]]]];*)
];
];

];
]

Transpose[list][[1]]

Count[list,{x_,y_,z_,w_}/;x<0.5]
Count[list,{x_,y_,z_,w_}/;y<0.5]
Count[list,{x_,y_,z_,w_}/;x<0.5||y<0.5]
Count[list,{x_,y_,z_,w_}/;x<0.5&&y<0.5]

list[[1]]

(* threshold=0.4 *)
listright=list

(* threshold=0.1 *)
listwrong=list

Length[listright]
Length[listwrong]

(* threshold=0.8 *)
listright=list

(* threshold=0.2 *)
listwrong=list

Length[listright]
Length[listwrong]

(* threshold=1.2 *)
listright=list

(* threshold=0.3 *)
listwrong=list

Length[listright]
Length[listwrong]

(* threshold=1.8 *)
listright=list

(* threshold=0.4 *)
listwrong=list

Length[listright]
Length[listwrong]

(* threshold=2.0 *)
listright=list

(* threshold=0.5 *)
listwrong=list

Length[listright]
Length[listwrong]

Transpose[listright][[1]]
Transpose[listwrong][[1]]

Intersection[Transpose[listright][[1]],Transpose[listwrong][[1]]]

zz=SequenceAlignment[Transpose[listright][[1]],Transpose[listwrong][[1]]]

Table[Depth[zz[[i]]],{i,Length[zz]}]

For[i=1,i≤Length[zz],i++,
Print[zz[[i]]];
]

Count[Flatten[stronginteraction],0]
Count[Flatten[stronginteraction],1]

(* nn=100, threshold 0.50 *)
Count[Flatten[stronginteraction],0]
Count[Flatten[stronginteraction],1]

(* nn=100, threshold 0.08 *)
Count[Flatten[stronginteraction],0]
Count[Flatten[stronginteraction],1]

(* nn=100, threshold 0.10 *)
Count[Flatten[stronginteraction],0]
Count[Flatten[stronginteraction],1]

(* nn=100, threshold 0.12 *)
Count[Flatten[stronginteraction],0]
Count[Flatten[stronginteraction],1]

(* nn=100, threshold 0.14 *)
Count[Flatten[stronginteraction],0]
Count[Flatten[stronginteraction],1]

(* nn=100, threshold 0.16 *)
Count[Flatten[stronginteraction],0]
Count[Flatten[stronginteraction],1]

(* nn=100, threshold 0.18 *)
Count[Flatten[stronginteraction],0]
Count[Flatten[stronginteraction],1]

(* nn=100, threshold 0.26 *)
Count[Flatten[stronginteraction],0]
Count[Flatten[stronginteraction],1]

(* nn=100, threshold 0.40 *)
Count[Flatten[stronginteraction],0]
Count[Flatten[stronginteraction],1]

(* nn=100, threshold 0.50 *)
Count[Flatten[stronginteraction],0]
Count[Flatten[stronginteraction],1]

Count[lista,x_/;x<0.04]

BinCounts[lista,{-0.05,0.05,0.01}]

ListPlot[BinCounts[lista,{-0.1,0.1,0.01}],Filling->Axis]