Description
三核苷酸是组成DNA序列的基本片段。具体来说,核苷酸一共有4种,分别用’A’,’G’,’C’,’T’来表示。而三核苷酸就是由3个核苷酸排列而成的DNA片段。三核苷酸一共有64种,分别是’AAA’,’AAG’,…,’GGG’。给定一个长度为L的DNA序列,一共可以分辨出(L-2)个三核苷酸。现在我们想用一些统计学的方法来进行一些分析,步骤如下:
1.对于这(L-2)个三核苷酸,我们从左到右给予编号,分别为1到L-2。
2.从这(L-2)个三核苷酸挑选一对出来,一共有(L-2)*(L-3)/2种可能。如果某一对三核苷酸是一样的,我们就记录他们之间的距离。他们之间的距离定义为他们的编号之差。
3.根据我们所记录的“样本数据”,我们现在需要计算样本数据的方差。方差的计算公式是S2=[(x1-X) 2+(x2-X) 2+…+(xn-X)2]/n, X=(x1+x2+…+xn)/n。如果样本的大小n=0,那么我们认为S2=X=0。
例如,我们要统计DNA序列’ATATATA’:
1. 为三核苷酸编号. L1: ATA, L2:TAT, L3:ATA, L4:TAT, L5:ATA.
2. (L1,L3)=2, (L1,L5)=4, (L3,L5)=2, (L2,L4)=2. 所以样本数据是2,4,2,2.
3. 样本数据平均值X=(2+4+2+2)/4=2.5.
方差S2=[(2-2.5)2+(4-2.5) 2+(2-2.5)2+(2-2.5)2]/4=0.75.
给定一个DNA序列,请你计算出它的方差。
Input
输入包含多组测试数据。第一行包含一个正整数T,表示测试数据数目。每组数据包含一个由’A’,’G’,’C’,’T’组成的字符串,代表要统计的DNA序列。DNA序列的长度大于等于3且不会超过100000。
Output
对每组测试数据,输出一行答案,为一个保留6位精度的实数,代表S2的值。如果你的答案和标准答案的“相对误差”小于1e-8,你的答案会被视为正确的答案。
Sample Input
1
ATATATA
Sample Output
0.750000
分析:设ave为平均数,k为有多少对,t[i]为有每对的值对,答案是
ANS=(ave*ave*k+sum(t[i]^2)-2*ave*sum(t[i]))/k
但是我们不用得到每个t[i],显然只要知道所有t[i]的和以及平方和。
设s[i]为每对的位置,f[i]为3-i,有多少个片段与i的相同,
t[i]=s[i]-s[j] {j为前面和该片段相同的位置}
显然可以记录前面的片段和一次计算,c为上一个位置
t[i]=s[i]*f[c]-sum(s[j])
至于平方和,可以考虑
t[i]^2=(s[i]-s[j])^2
然后套路一波。
注意:
平方和会炸longlong,一定要把k先算出来。
可能一个对都没有,(例如:GGG),此时应该输出0.000000
代码:
const
maxn=100001;
var
t,n:longint;
f,sums,last:array [0..maxn] of int64;
square_sums:array [0..maxn] of extended;
ls:array [0..63] of longint;
x,y,z,c:longint;
sumt,sum:int64;
ave,ans,square_sumt:extended;
s:ansistring;
function getnum(ch:char):longint;
begin
if ch='A' then exit(0);
if ch='G' then exit(1);
if ch='C' then exit(2);
if ch='T' then exit(3);
end;
procedure main;
var i,j,k:longint;
begin
readln(s);
n:=length(s);
fillchar(sums,sizeof(sums),0);
fillchar(square_sums,sizeof(square_sums),0);
fillchar(ls,sizeof(ls),0);
sumt:=0; square_sumt:=0; sum:=0;
x:=getnum(s[1]);
y:=getnum(s[2]);
for i:=3 to n do
begin
z:=getnum(s[i]);
c:=ls[x*16+y*4+z];
f[i]:=f[c]+1;
sums[i]:=sums[c]+i;
if f[c]>0 then
sumt:=sumt+i*f[c]-sums[c];
last[i]:=c;
ls[x*16+y*4+z]:=i;
x:=y; y:=z;
end;
for j:=0 to 63 do
sum:=sum+f[ls[j]]*(f[ls[j]]-1) div 2;
if sum=0 then begin writeln('0.000000'); exit; end;
for i:=3 to n do
begin
c:=last[i];
square_sums[i]:=square_sums[c]+i/sum*i;
if c<>0 then
square_sumt:=square_sumt+i/sum*i*f[c]+square_sums[c]-2*i/sum*sums[c];
end;
ave:=sumt/sum;
ans:=ave*ave-2*ave*sumt/sum+square_sumt;
writeln(ans:0:6);
end;
begin
assign(input,'tri.in');
assign(output,'tri.out');
reset(input);
rewrite(output);
readln(t);
while t>0 do
begin
main;
dec(t);
end;
close(input);
close(output);
end.