查看: 1417|回复: 0

[个人总经验] DNA解决方案

字体大小: 正常放大

longde

30 主题	8 听众	106 积分

TA的每日心情

	奋斗 2014-12-7 07:58

签到天数: 22 天

[LV.4]偶尔看看III

宣传员

群组: 2014年网络挑战赛交流

群组: 国赛讨论

群组: 2014美赛讨论

群组: 第三届数模基础实训

电梯直达

1^#

发表于 2015-8-17 22:42 |只看该作者 |正序浏览

|招呼Ta 关注Ta

DNA解决方案

摘要  本问题是一个“有人管理分类问题”．首先分别列举出20个学习样本序列中1字符串、2字符串、3字符串出现的频率，构成含41个变量的基本特征集，接着用主成分分析法从中提取出4个特征．然后用Fisher线性判别法进行分类，得出了所求20个人工制造序列及182个自然序列的分类结果如下：
20个人工序列：22, 23，25，27，29，34，35，36，37为A类，其余为B类．
1） 182个自然序列：1，4，8，10，27，29，32，41，43，48，54，63，70，72，75，76，81，86，90，92，102，110，116，119，126，131，144，150，157，159，160，161，162，163，164，165，166，169，170，182为B类，其余为A类．
最后通过检验证明所用的分类数学模型效率较高．

一、问题重述

人类基因组计划中DNA全序列草图是由4个字符A，T，C，G按一定顺序排成的长约30亿的字符序列，其中没有“断句”也没有标点符号．虽然人类对它知之甚少，但也发现了其中的一些规律性和结构．例如，在全序列中有一些是用于编码蛋白质的序列片段，即由这4个字符组成的64种不同的3字符串，其中大多数用于编码构成蛋白质的20种氨基酸．又例如，在不用于编码蛋白质的序列片段中，A和T的含量特别多些，于是以某些碱基特别丰富作为特征去研究DNA序列的结构也取得了一些结果．此外，利用统计的方法还发现序列的某些片段之间具有相关性，等等．这些发现让人们相信，DNA序列中存在着局部的和全局性的结构，充分发掘序列的结构对理解DNA全序列是十分有意义的．目前在这项研究中最普通的思想是省略序列的某些细节，突出特征，然后将其表示成适当的数学对象．
作为研究DNA序列的结构的尝试，提出以下对序列集合进行分类的问题：
1）请从20个已知类别的人工制造的序列（其中序列标号1～10 为A类，11～20为B类）中提取特征，构造分类方法，并用这些已知类别的序列，衡量你的方法是否足够好．然后用你认为满意的方法，对另外20个未标明类别的人工序列（标号21～40）进行分类，把结果用序号（按从小到大的顺序）标明他们的类别（无法分类的不写入）
2）同样方法对182个自然DNA序列(他们都较长)进行分类，像1）一样地给出分类结果．

二、模型的合理假设

1．各序列中DNA碱基三联组（即3字符串）的起始位置和基因表达不影响分类的结果．
2． 64种3字符串压缩为20组后不影响分类的结果．
3．较长的182个自然序列与已知类别的20个样本序列具有共同的特征．

三、模型建立与求解

研究DNA序列具有什么结构，其A，T，C，G 4个碱基排成的看似随机的序列中隐藏着什么规律，是解读人类基因组计划中DNA全序列草图的基础，也是生物信息学（Bioinformatcs）最重要的课题之一．
题目给出了20个已知为两个类别的人工制造的DNA序列，要求我们从中提取特征，构造分类方法，从而对20个未标明类别的人工DNA序列和182个自然DNA序列进行分类．这是模式识别中的“有人管理分类”问题，即事先规定了分类的标准和种类的数目，通过大批已知样本的信息处理找出规律，再用计算机预报未知．给出的已知类别的样本称为学习样本．对于此类问题，我们通过建立分类数学模型（这包括形成和提取特征以及制定分类决策）、考查分类模型的效率、预报未知这几个步骤来进行．

（一）特征的形成和提取

为了有效地实现分类识别，首先要根据被识别的对象产生一组基本特征，并对基本特征进行变换，得到最能反映分类本质的特征．这就是特征形成和提取的过程．在列举了尽可能完备的特征参数集之后，就要借助于数学的方法，使特征参数的数目（在保证分类良好的前提下）减到最小．这是因为：1.多余的特征参数不但没有多少好处，而且会带来噪音，干扰分类和数学模型的建立．2.为了保证样本数和特征参数个数的比值足够大，而又不必要用太多的样本，最好使特征参数的个数降至最少．模式识别计算一般要求样本数至少为变量数的3倍，否则结果不够可靠．本问题的学习样本数为20个，故特征参数的个数以6～8个为宜．
我们通过研究4个字符A,T,C,G在DNA序列中的排列、组合特性，主要是研究字符和字符串的排列在序列中出现的频率，从中提取DNA序列的结构特征参数．

1．特征的形成
分别列举一个字符，2个字符，3个字符的排列在序列中出现的频率，构成基本特征集．
（1）1个字符的出现频率
表1列出了20个样本中A，T，C，G这4个字符出现的频率．由于在不用于编码蛋白质的序列片段中，A和T的含量特别多些，因此我们将A和T是否特别丰富作为一个特征．在表1中，列出了A和T出现的频率之和．(程序见附录一)

                     表 1

            A          C          T          G          A+T
   1. 29.73       17.12    13.51    39.64    43.24
   2. 27.03       16.22    15.32    41.44    42.34
   3. 27.03       21.62       6.31    45.05    33.33
   4. 42.34       10.81    28.83    18.02    71.17
   5. 23.42       23.42    10.81    42.34    34.23
   6. 35.14       12.61    12.61    39.64    47.75
   7. 35.14       9.91    18.92    36.04    54.05
   8. 27.93       16.22    18.92    36.94    46.85
   9. 20.72       20.72    15.32    43.24    36.04
10. 18.18       27.27    13.64    40.91       31.82

   11. 35.45       4.55    50.00    10.00    85.45
   12. 32.73       2.73    50.00    14.55    82.73
   13. 25.45       10.00    51.82    12.73    77.27
   14. 30.00       8.18    50.00    11.82    80.00
   15. 29.09       .00    64.55       6.36    93.64
   16. 36.36       8.18    46.36       9.09    82.73
   17. 35.45       24.55    26.36    13.64    61.82
   18. 29.09       11.82    50.00       9.09    79.09
   19. 21.82       14.55    56.36       7.27    78.18
   20. 20.00       17.27    56.36       6.36    76.36


（2）2字符串的排列出现的频率
A，T，C，G这4个字符组成了16种不同的2字符串．表2列出了20个样本中各2字符串出现的频率．（用“滚动”算法，如ATTCG有AT,TT,TC,CG共4个2字符串）(程序与附录一类似)
表 2
   AA  AC  AT  AG  TA  TC  TG  TT  CA  CT  CC CG  GA  GT GC  GG
1. 9.01 9.01  3.60  8.11  4.50  .90  4.50  3.60  3.60  3.60 1.80  8.11 11.7 1 2.70  5.41 18.92
2. 9.91 7.21  3.60  5.41  2.70 1.80  5.41  5.41  4.50  1.80  .90  9.01  9.91 4.50  5.41 21.62
3. 5.41 11.71 3.60  5.41  2.70 1.80 .90 .90 5.41  .90  .90  14.41 13.51  .90 7.21 23.42
4.  18.92 5.41 11.71  5.41 10.81 1.80  5.41  10.81  5.41  1.80 .90 2.70  6.31 4.50  2.70 4.50
5. 6.31 8.11  1.80  7.21 1.80  2.70  2.70  3.60  5.41  4.50 2.70 10.81  9.91  .90  9.01 21.62
6.  15.32 2.70  6.31  9.91 3.60 1.80 1.80  5.41  4.50 .00  .00  8.11  10.81  .90 8.11 19.82
7.  15.32 1.80 10.81  7.21 4.50  2.70  6.31  5.41  .90 1.80  .90  6.31  13.51 .90 4.50 16.22
8. 8.11 3.60  6.31 9.91  5.41  3.60  2.70  7.21  2.70  3.60 1.80  8.11  10.81 1.80 7.2116.22
9. 9.01  .90  4.50  6.31  .00 3.60  7.21  4.50  3.60  2.70 2.70  11.71  7.21 3.60 13.5118.02
10.  6.36  3.64 1.82  6.36 1.82 5.45  2.73  3.64  5.45 3.64 4.55 13.64 4.55 3.64 13.64 18.18
11. 15.45  2.73 14.55 2.73 16.36  .91 1.82 30.00  .91 .91  .91  1.82 2.73  4.55 .00  2.73
12. 13.64  .91  10.91 6.36 15.45 1.82 1.82 30.91  .91  .91 .00  .91 2.73  7.27  .00 4.55
13. 6.36 4.55 10.00 4.55  12.73 1.82 2.73  34.55  2.73 2.73  1.82  1.8 2 3.64 4.55 1.82  2.73
14. 8.18 .91  12.73 7.27 13.64  6.36 1.82  28.18  2.73  4.55  .00  .91  5.45  4.55  .91  .91
15.13.64  .00  12.73 1.82 13.64  .00 2.73 48.18  .00  .00 .00  .00  1.82  3.64  .00  .91
16. 16.36  3.64  15.45  .9113.64  4.55 4.55  22.73  1.82  5.45  .00  .91  4.55  2.73  .00  1.82
17.17.27  5.45 10.91 1.82 10.00  6.36 4.55  5.45 4.55  7.27  9.09 2.73  3.64  2.73  3.64 3.64
18.8.18  7.27  11.82 1.82 15.45  1.82  .91  30.91  3.64  3.64 1.82  2.73 1.82 3.64  .91  2.73
19.2.73  2.73 13.64 1.82  14.55  9.09  .913  1.82  1.82  8.18  1.82  2.73  2.73  2.73  .91 .91
20. 6.36  6.36 6.36  .91  9.09  10.00  3.64 32.73  2.73  13.64  .91  .00 1.82  3.64 .00 .91

（3）3字符串的排列出现的频率
A，T，C，G这4个字符组成了64种不同的3字符串．这64种3字符串构成生物蛋白质的20种氨基酸．在参考文献[1]的Figur2中，给出了这20种氨基酸的编码（见图1）．因此，在计算3字符串的出现频率时，我们根据图1将代表同一种氨基酸的3字符串合成一类，只统计20类3字符串的出现频率．（不考虑字符串在序列片段中的起始位置，也采用“滚动”算法．如ACGTCC中就有ACG,CGT,GTC,TCC共4个3字符串）见表3．(程序与附录一类似)

Symmetries of the diamond code sort the 64 codons into 20 classes, indicated here by 20 colors. All the codons in each class specified the same amino acid.
   图1 Brian Hayes 在论文“The Invention of the Genetic Code”中给出的图形
            （注：图中DNA被转录为RNA,“U”代表“T”）

表 3

b1 b2    b3 b4 b5 b6 b7 b8 b9 b10  b11 b12 b13 b14 b15 b16 b17 b18 b19 b20
1 1.77  3.54  2.65  0.88  0.00  0.00 7.96  0.88  4.42 2.65  17.70 10.62 3.54  4.42  4.42  7.08 1.77  3.54  13.27  7.08
2 1.89  1.89  0.94  0.94  0.00  0.94 1.89  0.94  4.72 12.26  7.55  11.32  8.49  3.77 3.77  6.60 9.43  6.60  7.55  2.83
3 0.98  0.00  0.00  5.88  0.98  8.82 2.94  0.00  0.00 2.94 10.78  5.88 13.73  0.00  4.90  3.92 19.61  1.96  8.82  5.88
4 0.00  0.00  0.00  0.87  0.00  0.87 13.04 1.74  6.09 2.61  11.30 13.04 3.48  5.22 3.48  8.70 3.48  1.74 14.78， 7.83
5 2.86  0.00  0.00  3.81  0.95  3.81 3.81  0.00  3.81 3.81  9.52  9.52 12.38  2.86  9.52  4.76 7.62 2.86  7.62  9.52
6 0.00  0.00  0.88  2.63  0.00  1.75 13.16  0.88  4.39 1.75  14.04  9.65 7.02  5.26  4.39 11.40  2.63 1.75  10.53 6.14
7 1.92  0.00  0.00  2.88  0.96  4.81 2.88 0.00  1.92 4.81  12.50  6.73  13.46  1.92  6.73 4.81 10.58  3.85  9.62  7.69
8 2.56  3.42  0.00  0.85  0.85  0.85 12.82  0.85  1.71 0.85  20.51  2.56  3.42 9.40  5.98 11.11  0.85 4.27 11.97  3.42
9 0.00  0.00  0.00  2.97  2.97  9.90 2.97 0.00  0.99 3.96  6.93  1.98 13.86  1.98  2.97  3.96 23.76  2.97  8.91  6.93
10 1.87  0.93  3.74  2.80  0.00  0.00 2.80 0.00  7.48 8.41  9.35  7.48 3.74 14.95 12.15  0.00 2.80 4.67  7.48  7.48

11 0.00  0.89  0.00  0.00  0.00  1.79  8.04 0.00  5.36  4.46 15.18  8.04  8.93 4.46  3.57  8.04 4.46 6.25 13.39  5.36
12 2.73  0.00  0.91  2.73  0.91  3.64  4.55 3.64  3.64  1.82 9.09 5.45  3.64 5.45  6.36  7.27 8.18 5.45 10.91  9.09
13 1.80  0.90  0.90  0.90  0.00  0.90  9.01 0.00  3.60  7.21 14.41  8.11  7.21 6.31  7.21  4.50 1.80 7.21  11.71 4.50
14 2.94  0.00  0.00  5.88  0.00  6.86  1.96 0.00  3.92  6.86 3.92 9.80  13.73 0.98  5.88  2.94 10.78  0.98  1 0.78 9.80
15 2.91  1.94  2.91  1.94  0.00  5.83  1.94 0.00  1.94  9.71 5.83 8.74  10.68  1.94 3.88  3.88 8.74 2.91  11.65 10.68
16 2.86  0.95  0.00 11.43  1.90  1.90  2.86 0.00 4.76  3.81  5.71 8.57 8.57 6.67  9.52 4.76  5.71 2.86  7.62  7.62
17 1.92  0.96  1.92  4.81  1.92  3.85  1.92 0.96  0.96 6.73  4.81 8.65  10.58  2.88  6.73 2.88 9.62 6.73  8.65  7.69
18 1.71  0.85  1.71  0.85  0.85  2.56  16.24  0.85  1.71 0.85  16.24  5.13  6.84 5.98  3.42 11.11  1.71 5.13  11.11 3.42
19 0.94  0.94  1.89  0.94  0.94  0.94  1.89 0.94  10.38  7.55  5.66 9.43  8.49 8.49  7.55 5.66 6.60  11.32  6.60  0.94
20 0.86  0.86  0.00  1.72  0.86  0.86  17.24  0.86  2.59  1.72 15.52  7.76  5.17 3.45  4.31 9.48 5.17 5.17  9.48  5.17

其中 b1 =aaa+ata    b2=aca+aga    b3=cac+ctc    b4=ccc+cgc
      b5 =gag+gtg    b6=gcg+ggg    b7=tat+ttt    b8=tct+tgt
      b9 =aac+caa+atc+cta             b10=aag+gaa+atg+gta
      b11=aat+taa+att+tta             b12=acc+cca+agc+cga
      b13=acg+gac+ctg+gtc             b14=act+tca+agt+tga
      b15=cag+gac+ctt+ttc             b16=cat+tac+ctt+ttc
      b17=ccg+gcc+cgg+ggc             b18=cct+tcc+cgt+tgc
      b19=gat+tag+gtt+ttg             b20=gct+tcg+ggt+tgg

综合起来，形成了有41个变量的基本特征集．

2. 特征的提取
上述基本特征集中有41个变量，即样本处于一个高维空间中．特征的提取就是通过变换的方法用低维空间来表示样本，使得X的大部分特性能由Y来表达，即将p维随机向量X变换成q维随机向量 Y（q<p）．我们用主成分分析法进行特征的提取，其步骤是：
（1）求X的均方差矩阵V的特征根，记为：
λ1≥λ2≥…≥λk＞0 λk+1=…=λP=0
（2）求λ1, λ2…λK对应的标准正交的特征向量r1，r2,…,rk
得到第i个主成分为yi=riX, i=1,2,…,k.
（3）求第i个主成分的贡献率ui=λi/  λj, i=1,2,…,k, 及前m个主成分的累计贡献率vm= ui .
（4）求得q，使得Vq≥V0(V0一般在0.85到1之间)，则取
W=(r1,r2,…,rq)
Y=XW
第3步所求的贡献率，代表主成分表达X的能力，贡献率越大，对应的主成分表达X的能力越强．只要前q个主成分的累计贡献率超过给定的百分比V．就可以用低维特征Y=（y1,y2, …,yq）来反映高维特征（x1,x2,…,xp）的变化特性．
现将反映20个已知类别样本的41个特征的随机向量X进行特征提取．
计算得前4个主成分的累计贡献率为96%，故提取特征为4个变量，取
W=（r1,r2,r3,r4），则Y=XW，Y的4个分量就是从基本特征集提取所得的特征参数向量．(程序及结果见附录二)

（二）分类决策的制定

前面已选取了特征参数，把特征参数张成的多维空间称为特征空间．分类决策就是在特征空间中用统计的方法把被识别对象归为某一类别．基本作法是在学习样本集的基础上确定某个判决规则，使按这种判决规则对被甄别对象进行分类所造成的错误识别率最小或引起的损失最少．
这里，我们的分类决策选取Fisher线性判别法．即选取线性判别函数U(x)，使得：
   U(x)={E1[U(x)]-E2[U(x)]}2/{D1 [U(x)]+D2[U(x)]}=max       (1)
其中Ei与Di分别表示母体i的期望和方差运算，i=1，2．
（1）式的含义是：构造一个线性判别函数U(x)对样本进行分类，使得平均出错概率最小．即应在不同母体下，使U(x)的取值尽量分开．具体地说，要使母体间的差异    (E1(U(x))-E2(U(x)))2相对于母体内的差异D1[U(x)]+D2[U(x)] 为最大．取
      U(x)=( 1- 2)T(∑1+∑2)-1X
就可满足(1)．其中 i为第i类母体的均值矩阵的估计，∑i为第i类母体的方差矩阵的估计．取分类门槛值为：
            U0=U(α* 1+(1-α)* 2)
其中0<α<1，本问题中两类样本的个数相等，可取 α=1/2．若U( 1)>U0,U( 2)<U0 ,  则当U(X)>U0., 就认为X取自母体1；当U(X)<U0, 就认为X取自母体2．
用上面得出的4个主成分构成的特征组和此分类决策，对20个学习样本进行分类，能得出正确的结果．但是，若取W=（r1,r2,r3），求Y=XW，以Y的3个分量作为特征参数向量，再用Fisher线性判别法对20个学习样本进行分类，则第四个样本不能正确分类．
因此，得出分类的数学模型为：
（1）特征选取：取W=（r1,r2,r3,r4），求Y=XW，得出特征参数向量就是Y的4个列向量．其中X是反映20个学习样本的41个特征的随机向量．
（2）分类决策：Fisher线性判别法．

（三）分类模型的有效性考察

前面建立的分类数学模型对20个学习样本进行了正确分类．为了进一步考查分类模型的有效性和可靠性，我们采用的方法是：预先留一部分学习样本不参加训练，然后用分类决策模型对其作预报，将预报成功率作为预报能力的指标．
每次取出一个学习样本，以其余学习样本作训练集，用分类决策模型对取出的一个样本作预报，同时对给出的后20种样本作预报．结果见表4．

表  4
取出样品序号取出样本类别预报后20组样本中A类序号预报
1       A 22，23，25，27，29，34，35，36，37
2       A 22，23，25，27，29，34，35，36，37
3       A 22，23，25，27，29，34，35，36，37
4       A 23，25，27，29，34，35，36，37
5       A 22，23，25，27，29，34，35，36，37
6       A 22，23，25，27，29，34，35，36，37
7       A 22，23，25，27，29，34，35，36，37
8       A 22，23，25，27，29，34，35，36，37
9       A 22，23，25，27，29，34，35，36，37
10       A 22，23，25，27，29，34，35，36，37
11       B 22，23，25，27，29，34，35，36，37
12       B 22，23，25，27，29，34，35，36，37
13       B 22，23，25，27，29，34，35，36，37
14       B 22，23，25，27，29，34，35，36，37
15       B 22，23，25，27，29，34，35，36，37，39
16       B 22，23，25，27，29，34，35，36，37
17       B 22，23，25，27，29，34，35，36，37，30，39
18       B 22，23，25，27，29，34，35，36，37
19       B 22，23，25，27，29，34，35，36，37
20       B 22，23，25，27，29，34，35， 37

从表4可以看出：
1．每次取出一个学习样本，以其余学习样本作训练集，用分类模型对该学习样本的预报的成功率是100%．
2．每次取出一个学习样本，以其余学习样本作训练集，用分类模型对未知类别的第21~40个样本进行预报，其结果有以下特点：
（1）除分别取出4、15、17，20的预报结果不同外，分别取出其余16中一个，预报结果均为：22，23，25，27，29，34，35，36，37，占80%．
（2）分别取出4、15、20的预报结果，与（1）的结果相比，只有一个样本的差异，占15%．
（3）取出17的预报结果，与（1）的结果相比，有两个样本的差异，占5%．
第一种结果和第二种结果非常接近，合计占总数的95%．只有第三组的这一个结果有较大差异，占总数的5%．
由以上检验得出结论：所建立的分类数学模型分类效果很好．

（四）未知样本的预报
现在用前面建立的数学模型对题目所给的未知类型的20个人工序列和182个自然序列进行预报．（程序见附录三）
结果为：
1） 20个人工序列的类别
A类：22，23，25，27，29，34，35，36，37
B类：21、24、26、28、30、31、32、33、38、39、40
2） 182个自然序列的类别
A类：（共142个）2，3，5，6，7，9，11，12，13，14，15，16，17，18，19，20，21，22，23，24，25，26，28，30，31，33，34，35，36，37，38，39，40，42，44，45，46，47，49，50，51，52，53，55，56，57，58，59，60，61，62，64，65，66，67，68，69，71，73，74，77，78，79，80，82，83，84，85，87，88，89，91，93，94，95，96，97，98，99，100，101，103，104，105，106，107，108，109，111，112，113，114，115，117，118，120，121，122，123，124，125，127，128，129，130，132，133，134，135，136，137，138，139，140，141，142，143，145，146，147，148，149，151，152，153，154，155，156，158，167，168，171，172，173，174，175，176，177，178，179，180，181
B类：（共40个）1，4，8，10，27，29，32，41，43，48，54，63，70，72，75，76，81，86，90，92，102，110，116，119，126，131，144，150，157，159，160，161，162，163，164，165，166，169，170，182

四、模型的优缺点分析

优点：
1．针对`“有人管理分类”问题，成功地建立解决这类难题的数学模型，并可立即运用到实践中去．
2．仅用4个特征参数即圆满解决了较为复杂的分类问题．而且模型假设条件少，因而能准确地反映实际情况，可靠性高．
3．采用模块化分析，逐渐深入，提高了准确性．
4．突出特征，假设合理，避免了在一些细节问题上的纠缠．
缺点：
由于只考虑了DNA样本序列中1字符串、2字符串、3字符串出现的频率作为特征，DNA序列的分类不一定与实际情况完全相符．（可以由科学家用物理的或化学的方法测定，作为补充）．

五、模型的改进方向及推广

模型的改进：因为模型没考虑DNA序列的实际特性，当序列变得很多很长很复杂时，分类的准确性会降低而不可用，因此应增加对DNA序列的生物特性的考虑．
模型的推广：该模型对一般的“有人管理分类”问题的求解有重要意义．对研究DNA序列的规律性和结构提供了一种有效的分类模型．对人类基因组的研究有现实意义，有利于加快科研步伐．

六、参考文献

[１] Brain Hayes(美)．The Invention of the Genetic Code． American Scientist—Computing 　Science，Jan.-Feb., 1998
[２] 萧树铁主编．数学实验．北京：高等教育出版社，1999
[３] 复旦大学．概率论第二册—数理统计．北京：高等教育出版社，1985
[４] William F. Lucas 主编．生命科学模型。长沙：国防科技大学出版社，1996
[５] 徐光煇主编．运筹学基础手册．北京：科学出版社，1999
[６] 姜启源主编．数学模型．北京：高等教育出版社，1993

七、附录



附录一 1个字符出现频率的计算程序]
CHARACTER*121 LINE(40)
integer a,c,t,g,at
READ*,LINE
DO 20 II=1,40
iii=ii+20
A=0
C=0
T=0
G=0
DO 10 I=1,121
IF(LINE(ii)(I:I).EQ.’a’)THEN
A=A+1
else if(line(ii)(I:I).eq.’c’)then
c=c+1
else if(line(ii)(I:I).eq.’t’)then
t=t+1
else if(line(ii)(I:I).eq.’g’)then
g=g+1
END IF
10 continue
at=a+t
actg=a+c+t+g
aa=a/actg*100.
cc=c/actg*100.
tt=t/actg*100.
gg=g/actg*100.
aatt=at/actg*100.
open(5,file='t1.dat',status='old')
write(5,1)aa,cc,tt,gg
1 format(1x,4f7.2)
20 CONTINUE
END
附录二  基本特征量的提取程序及结果
d=[ 27.43  19.47  36.28  16.81  63.72;
      28.85  24.04  22.12  25.00  50.96;
      17.65  25.49  18.63  38.24  36.27;
      20.87  19.13  40.87  19.13  61.74;
      24.76  22.86  21.90  30.48  46.67;
      21.93  21.05  38.60  18.42  60.53;
      23.08  20.19  23.08  33.65  46.15;
      25.64  14.53  44.44  15.38  70.09;
      14.85  21.78  18.81  44.55  33.66;
      28.97  24.30  25.23  21.50  54.21;
      24.11  17.86  35.71  22.32  59.82;
      17.43  22.94  33.03  26.61  50.46;
      27.03  18.92  33.33  20.72  60.36;
      23.53  23.53  16.67  36.27  40.20;
      24.27  21.36  20.39  33.98  44.66;
      22.86  30.48  20.95  25.71  43.81;
      21.36  25.24  20.39  33.01  41.75;
      22.22  17.09  43.59  17.09  65.81;
      27.36  28.30  23.58  20.75  50.94;
      19.83  19.83  43.10  17.24  62.93];
   dd=[ 5.31 4.42 7.96 8.85 9.73 6.19 1.77  18.58 6.19 4.42 4.42 4.42 6.19 4.42 4.42 1.77;
      7.69 9.62 3.85 7.69 9.62 3.85 .96 6.73 2.88 1.92 7.69  11.54 7.69 8.65 2.88 4.81;
      2.94 3.92 5.88 4.90 3.92 2.94 1.96 9.80 .00 1.96  12.75 9.80  10.78 .98 4.90  21.57;
      1.74 4.35 3.48  11.30  13.04 1.74 2.61  22.61 2.61 9.57 4.35 2.61 3.48 4.35 8.70 2.61;
      6.67 3.81 3.81 9.52 5.71 1.90 4.76 9.52 7.62 4.76 7.62 2.86 4.76 3.81 9.52  12.38;
      3.51 3.51 5.26 9.65 7.89 4.39 1.75  24.56 7.89 6.14 1.75 4.39 2.63 2.63  11.40 1.75;
      5.77 4.81 4.81 7.69 6.73 2.88 2.88  10.58 2.88 2.88 7.69 6.73 7.69 4.81 4.81  15.38;
      3.42 5.13 9.40 6.84  11.97 5.13 3.42  23.93 2.56 6.84 2.56 2.56 7.69 3.42 1.71 2.56;
      1.98 1.98 3.96 6.93 3.96 2.97 2.97 8.91 1.98 .99 8.91 8.91 6.93 4.95 7.92  24.75;
      9.35 5.61 2.80  10.28 7.48 5.61 5.61 6.54 8.41 7.48 2.80 5.61 3.74 8.41 9.35 .00;
      2.68 5.36 4.46  11.61  15.18 1.79 .89  16.96 3.57 6.25 3.57 4.46 2.68 7.14 7.14 5.36;
      5.50 2.75 2.75 6.42 6.42 7.34 4.59  13.76 4.59 5.50 6.42 6.42 .92  10.09 6.42 8.26;
      5.41 7.21 7.21 7.21  10.81 1.80 5.41  15.32 3.60 4.50 2.70 7.21 7.21 6.31 6.31 .90;
      7.84 4.90 .98 8.82 4.90 .98 2.94 7.84 2.94 3.92 9.80 6.86 7.84 3.92 6.86  17.65;
      5.83 4.85 3.88 9.71 7.77 3.88 1.94 6.80 3.88 2.91 3.88 9.71 6.80 6.80 8.74  11.65;
      4.76 3.81 1.90  12.38 8.57 5.71 .00 6.67 5.71 3.81  10.48  10.48 3.81 8.57 9.52 2.86;
      3.88 2.91 2.91  10.68 5.83 .97 6.80 5.83 5.83 5.83 9.71 3.88 4.85 5.83  11.65  10.68;
      3.42 9.40 5.98 3.42  10.26 1.71 4.27  27.35 5.13 3.42 4.27 3.42 2.56 6.84 1.71 5.98;
      8.49 5.66 4.72 8.49 4.72 8.49 2.83 6.60  11.32 1.89 9.43 5.66 2.83 9.43 4.72 3.77;
      3.45 7.76 4.31 4.31  10.34 .86 3.45  27.59 1.72 6.03 8.62 3.45 4.31 5.17 1.72 6.03];
  ddd=[ 1.77 3.54 2.65 .88 .00 .00 7.96 .88 4.42 2.65  17.70  10.62 3.54 4.42 4.42 7.08 1.77 3.54  13.27 7.08;
      1.92 1.92 .96 .96 .00 .96 1.92 .96 4.81  12.50 7.69  11.54 8.65 3.85 3.85 6.73 9.62 6.73 7.69 2.88;
         .98 .00 .00 5.88 .98 8.82 2.94 .00 .00 2.94  10.78 5.88  13.73 .00 4.90 3.92  19.61 1.96 8.82 5.88;
         .00 .00 .00 .87 .00 .87  13.04 1.74 6.09 2.61  11.30  13.04 3.48 5.22 3.48 8.70 3.48 1.74  14.78 7.83;
      2.86 .00 .00 3.81 .95 3.81 3.81 .00 3.81 3.81 9.52 9.52  12.38 2.86 9.52 3.81 7.62 2.86 7.62 9.52;
         .00 .00 .88 2.63 .00 1.75  13.16 .88 4.39 1.75  14.04 9.65 7.02 5.26 4.39  11.40 2.63 1.75  10.53 6.14;
      1.92 .00 .00 2.88 .96 4.81 2.88 .00 1.92 4.81  12.50 6.73  13.46 1.92 6.73 4.81  10.58 3.85 9.62 7.69;
      2.56 3.42 .00 .85 .85 .85  12.82 .85 1.71 .85  20.51 2.56 3.42 9.40 5.98  11.11 .85 4.27  11.97 3.42;
         .00 .00 .00 2.97 2.97 9.90 2.97 .00 .99 3.96 6.93 1.98  13.86 1.98 2.97 3.96  23.76 2.97 8.91 6.93;
      1.87 .93 3.74 2.80 .00 .00 2.80 .00 7.48 8.41 9.35 7.48 3.74  14.95  12.15 .00 2.80 4.67 7.48 7.48;
         .00 .89 .00 .00 .00 1.79 8.04 .00 5.36 4.46  15.18 8.04 8.93 4.46 3.57 8.04 4.46 6.25  13.39 5.36;
      2.75 .00 .92 2.75 .92 3.67 4.59 3.67 3.67 1.83 9.17 5.50 3.67 5.50 6.42 7.34 8.26 5.50  11.01 9.17;
      1.80 .90 .90 .90 .00 .90 9.01 .00 3.60 7.21  14.41 8.11 7.21 6.31 7.21 4.50 1.80 7.21  11.71 4.50;
      2.94 .00 .00 5.88 .00 6.86 1.96 .00 3.92 6.86 3.92 9.80  13.73 .98 5.88 2.94  10.78 .98  10.78 9.80;
      2.91 1.94 2.91 1.94 .00 5.83 1.94 .00 1.94 9.71 5.83 8.74  10.68 1.94 3.88 3.88 8.74 2.91  11.65  10.68;
      2.86 .95 .00  11.43 1.90 1.90 2.86 .00 4.76 3.81 5.71 8.57 8.57 6.67 9.52 4.76 5.71 2.86 7.62 7.62;
      1.94 .97 1.94 4.85 1.94 3.88 1.94 .97 .97 6.80 4.85 8.74  10.68 2.91 6.80 2.91 9.71 6.80 8.74 7.77;
      1.71 .85 1.71 .85 .85 2.56  16.24 .85 1.71 .85  16.24 5.13 6.84 5.98 3.42  11.11 1.71 5.13  11.11 3.42;
         .94 .94 1.89 .94 .94 .94 1.89 .94  10.38 7.55 5.66 9.43 8.49 8.49 7.55 5.66 6.60  11.32 6.60 .94;
         .86 .86 .00 1.72 .86 .86  17.24 .86 2.59 1.72  15.52 7.76 5.17 3.45 4.31 9.48 5.17 5.17 9.48 5.17];
x=[ 29.73  17.12  13.51  39.64  43.24;
      27.03  16.22  15.32  41.44  42.34;
      27.03  21.62 6.31  45.05  33.33;
      42.34  10.81  28.83  18.02  71.17;
      23.42  23.42  10.81  42.34  34.23;
      35.14  12.61  12.61  39.64  47.75;
      35.14 9.91  18.92  36.04  54.05;
      27.93  16.22  18.92  36.94  46.85;
      20.72  20.72  15.32  43.24  36.04;
      18.18  27.27  13.64  40.91  31.82;;
      35.45 4.55  50.00  10.00  85.45;
      32.73 2.73  50.00  14.55  82.73;
      25.45  10.00  51.82  12.73  77.27;
      30.00 8.18  50.00  11.82  80.00;
      29.09 .00  64.55 6.36  93.64;
      36.36 8.18  46.36 9.09  82.73;
      35.45  24.55  26.36  13.64  61.82;
      29.09  11.82  50.00 9.09  79.09;
      21.82  14.55  56.36 7.27  78.18;
      20.00  17.27  56.36 6.36  76.36];
xx=[ 9.01 9.01 3.60 8.11 4.50 .90 4.50 3.60 3.60 3.60 1.80 8.11  11.71 2.70 5.41  18.92;
      9.91 7.21 3.60 5.41 2.70 1.80 5.41 5.41 4.50 1.80 .90 9.01 9.91 4.50 5.41  21.62;
      5.41  11.71 3.60 5.41 2.70 1.80 .90 .90 5.41 .90 .90  14.41  13.51 .90 7.21  23.42;
      18.92 5.41  11.71 5.41  10.81 1.80 5.41  10.81 5.41 1.80 .90 2.70 6.31 4.50 2.70 4.50;
      6.31 8.11 1.80 7.21 1.80 2.70 2.70 3.60 5.41 4.50 2.70  10.81 9.91 .90 9.01  21.62;
      15.32 2.70 6.31 9.91 3.60 1.80 1.80 5.41 4.50 .00 .00 8.11  10.81 .90 8.11  19.82;
      15.32 1.80  10.81 7.21 4.50 2.70 6.31 5.41 .90 1.80 .90 6.31  13.51 .90 4.50  16.22;
      8.11 3.60 6.31 9.91 5.41 3.60 2.70 7.21 2.70 3.60 1.80 8.11  10.81 1.80 7.21  16.22;
      9.01 .90 4.50 6.31 .00 3.60 7.21 4.50 3.60 2.70 2.70  11.71 7.21 3.60  13.51  18.02;
      6.36 3.64 1.82 6.36 1.82 5.45 2.73 3.64 5.45 3.64 4.55  13.64 4.55 3.64  13.64  18.18;
      15.45 2.73  14.55 2.73  16.36 .91 1.82  30.00 .91 .91 .91 1.82 2.73 4.55 .00 2.73;
      13.64 .91  10.91 6.36  15.45 1.82 1.82  30.91 .91 .91 .00 .91 2.73 7.27 .00 4.55;
      6.36 4.55  10.00 4.55  12.73 1.82 2.73  34.55 2.73 2.73 1.82 1.82 3.64 4.55 1.82 2.73;
      8.18 .91  12.73 7.27  13.64 6.36 1.82  28.18 2.73 4.55 .00 .91 5.45 4.55 .91 .91;
      13.64 .00  12.73 1.82  13.64 .00 2.73  48.18 .00 .00 .00 .00 1.82 3.64 .00 .91;
      16.36 3.64  15.45 .91  13.64 4.55 4.55  22.73 1.82 5.45 .00 .91 4.55 2.73 .00 1.82;
      17.27 5.45  10.91 1.82  10.00 6.36 4.55 5.45 4.55 7.27 9.09 2.73 3.64 2.73 3.64 3.64;
      8.18 7.27  11.82 1.82  15.45 1.82 .91  30.91 3.64 3.64 1.82 2.73 1.82 3.64 .91 2.73;
      2.73 2.73  13.64 1.82  14.55 9.09 .91  31.82 1.82 8.18 1.82 2.73 2.73 2.73 .91 .91;
      6.36 6.36 6.36 .91 9.09  10.00 3.64  32.73 2.73  13.64 .91 .00 1.82 3.64 .00 .91];
  xxx=[ 5.41 .90 2.70 .90 5.41 3.60 .90 1.80 2.70 8.11 4.50 1.80  25.23 3.60 3.60 5.41  13.51 .00 3.60 4.50;
      2.70 2.70 .00 .00 3.60 6.31 2.70 .90 7.21 7.21 6.31 1.80  18.92 .90 6.31 1.80  14.41 .00 3.60  10.81;
      2.70 2.70 2.70 .00 3.60 6.31 .00 .90 4.50 5.41 1.80 .90  29.73 .00 5.41 4.50  22.52 .00 1.80 2.70;
      15.32 6.31 .00 .00 .00 .90 9.01 1.80 6.31  10.81  12.61 3.60 4.50 1.80 2.70 5.41 1.80 1.80 7.21 6.31;
      3.60 1.80 2.70 .00 5.41 7.21 .90 .00 4.50 1.80 2.70 3.60  20.72 1.80 6.31 4.50  19.82 1.80 1.80 7.21;
      9.01 .90 .90 .00 2.70 5.41 4.50 .00 2.70  13.51 6.31 .00  25.23 .90 1.80 1.80  16.22 .00 2.70 3.60;
      9.01 1.80 .00 .00 1.80 4.50 4.50 .90 3.60  16.22 8.11 .00  17.12 2.70 1.80 1.80  10.81 .90 6.31 6.31;
      2.70 1.80 .90 .90 2.70 3.60 2.70 .90 4.50 9.91 8.11 3.60  18.92 .90 2.70 4.50  12.61 .90 7.21 8.11;
      5.41 .00 .90 1.80 5.41 9.01 1.80 .90 3.60 6.31 1.80 3.60  11.71 2.70 2.70 2.70  20.72 1.80 4.50  10.81;
      3.64 .91 2.73 6.36 3.64  10.91 .91 1.82 3.64 2.73 2.73 .91  17.27 .00 4.55 4.55  17.27 4.55 1.82 7.27;
      9.09 .91 .00 .00 .00 .00  24.55 .00 3.64 6.36  33.64 .91 4.55 1.82 .00 1.82 .00 2.73 5.45 2.73;
      2.73 .91 .00 .00 .00 .00  19.09 .00 1.82 8.18  37.27 .00 4.55 4.55 .00 2.73 .00 .91  10.00 5.45;
         .91 2.73 .00 .00 .00 .00  27.27 1.82 1.82 5.45  26.36 2.73 4.55 2.73 4.55 5.45 1.82 2.73 5.45 1.82;
      6.36 5.45 .00 .00 1.82 .00  20.00 5.45 2.73 2.73  24.55 .00 1.82 3.64 3.64 8.18 .91 .91 9.09 .91;
      11.82 .91 .00 .00 1.82 .00  47.27 1.82 .00 3.64  25.45 .00 .91 .91 .00 .00 .00 .00 2.73 .91;
      10.00 2.73 .91 .00 .00 .00  14.55 4.55 5.45 3.64  31.82 .91 .91 3.64 1.82 6.36 .00 .00 7.27 3.64;
      10.91 .91 3.64 3.64 .00 .91 8.18 2.73  12.73 9.09  11.82 3.64 3.64 6.36 1.82 1.82 6.36 6.36 1.82 1.82;
      4.55 4.55 .00 .00 .91 .91  21.82 .91 4.55 .91  29.09 .00 3.64 1.82 .91  10.91 2.73 4.55 4.55 .91;
      3.64 .91 1.82 .91 .91 .00  25.45 5.45 3.64 .00  21.82 1.82 1.82 3.64 .91  13.64 .91 2.73 5.45 2.73;
      2.73 .91 5.45 .00 .00 .00  23.64  10.00 6.36 1.82  13.64 .00 1.82 8.18 1.82  13.64 .00 1.82 6.36 .00];
   ffx=[x xx xxx];
   ffd=[d dd ddd];
   cx=cov(ffx);
   [vx,ex]=eig(cx);
   ex1=eig(cx);
   e1=mean(ex1)*41;
   ex2=ex1(38：41,

;
   e2=mean(ex2)*7;
   e2/e1
      vx1=[vx(:,38:41)];
   s=ffx*vx1;ss=ffd*vx1;
   x=s(1:10,

;
y=s(11:20,

;
u1=mean(x);u2=mean(y);
u1-u2;
z=8/9*(cov(x)+cov(y));
ux=0.5*(u1-u2)*inv(z);
u12=0.5*u1+0.5*u2;
u0=ux*u12.';
la=0;
for i=1:10
p(i)=ux*ss(i,

.';
tx(i)=ux*x(i,

.';
fy(i)=ux*y(i,

.';
  if p(i)>u0
   pbd(i)=1;
   la=la+1;
      else
         pbd(i)=2 ;
end
      if tx(i)>u0
      lbx(i)=1 ;
      else
         lbx(i)=2;
      end
      if fy(i)>u0
      lby(i)=1 ;
      else
         lby(i)=2 ;
      end
  for n=11:20
   p(n)=ux*ss(n,

';
   if p(n)>u0
      pbd(n)=1 ;
      la=la+1;
      else
         pbd(n)=2;
      end
tx ,fy ,p
pbd,lbx,lby
ans =0.9847
u0 =-2.4812
tx= Columns 1 through 7
8.2471 9.7074 10.8780 3.8672 9.3837 9.7612 9.2014
Columns 8 through 10
6.2700 11.6489 5.4181
fy =Columns 1 through 7
-15.2467  -15.2121  -14.2828 -8.0112  -13.4839  -11.1970  -11.2608
Columns 8 through 10
-15.0827  -14.9635  -15.2662
p =Columns 1 through 7
-6.5147 -3.6869 0.7514 -6.0838 0.3758 -6.7805 0.1074
Columns 8 through 14
-8.1194 5.0825 -6.1039 -7.0908 -2.7297 -6.0715 4.1447
Columns 15 through 20
4.5919 -4.2199 0.9096 -9.2269 -8.1303  -10.7112
pbd =Columns 1 through 12
2    2    1    2    1    2    1    2    1    2    2    2
Columns 13 through 20
2    1    1    2    1    2    2    2
lbx =1    1    1    1    1    1    1    1    1    1
lby = 2    2    2    2    2    2    2    2    2    2

附录三对未知序列进行分类的运算程序
d=[ 27.43  19.47  36.28  16.81  63.72;
      28.85  24.04  22.12  25.00  50.96;
      17.65  25.49  18.63  38.24  36.27;
      20.87  19.13  40.87  19.13  61.74;
      24.76  22.86  21.90  30.48  46.67;
      21.93  21.05  38.60  18.42  60.53;
      23.08  20.19  23.08  33.65  46.15;
      25.64  14.53  44.44  15.38  70.09;
      14.85  21.78  18.81  44.55  33.66;
      28.97  24.30  25.23  21.50  54.21;
      24.11  17.86  35.71  22.32  59.82;
      17.43  22.94  33.03  26.61  50.46;
      27.03  18.92  33.33  20.72  60.36;
      23.53  23.53  16.67  36.27  40.20;
      24.27  21.36  20.39  33.98  44.66;
      22.86  30.48  20.95  25.71  43.81;
      21.36  25.24  20.39  33.01  41.75;
      22.22  17.09  43.59  17.09  65.81;
      27.36  28.30  23.58  20.75  50.94;
      19.83  19.83  43.10  17.24  62.93];
   dd=[ 5.31 4.42 7.96 8.85 9.73 6.19 1.77  18.58 6.19 4.42 4.42 4.42 6.19 4.42 4.42 1.77;
      7.69 9.62 3.85 7.69 9.62 3.85 .96 6.73 2.88 1.92 7.69  11.54 7.69 8.65 2.88 4.81;
      2.94 3.92 5.88 4.90 3.92 2.94 1.96 9.80 .00 1.96  12.75 9.80  10.78 .98 4.90  21.57;
      1.74 4.35 3.48  11.30  13.04 1.74 2.61  22.61 2.61 9.57 4.35 2.61 3.48 4.35 8.70 2.61;
      6.67 3.81 3.81 9.52 5.71 1.90 4.76 9.52 7.62 4.76 7.62 2.86 4.76 3.81 9.52  12.38;
      3.51 3.51 5.26 9.65 7.89 4.39 1.75  24.56 7.89 6.14 1.75 4.39 2.63 2.63  11.40 1.75;
      5.77 4.81 4.81 7.69 6.73 2.88 2.88  10.58 2.88 2.88 7.69 6.73 7.69 4.81 4.81  15.38;
      3.42 5.13 9.40 6.84  11.97 5.13 3.42  23.93 2.56 6.84 2.56 2.56 7.69 3.42 1.71 2.56;
      1.98 1.98 3.96 6.93 3.96 2.97 2.97 8.91 1.98 .99 8.91 8.91 6.93 4.95 7.92  24.75;
      9.35 5.61 2.80  10.28 7.48 5.61 5.61 6.54 8.41 7.48 2.80 5.61 3.74 8.41 9.35 .00;
      2.68 5.36 4.46  11.61  15.18 1.79 .89  16.96 3.57 6.25 3.57 4.46 2.68 7.14 7.14 5.36;
      5.50 2.75 2.75 6.42 6.42 7.34 4.59  13.76 4.59 5.50 6.42 6.42 .92  10.09 6.42 8.26;
      5.41 7.21 7.21 7.21  10.81 1.80 5.41  15.32 3.60 4.50 2.70 7.21 7.21 6.31 6.31 .90;
      7.84 4.90 .98 8.82 4.90 .98 2.94 7.84 2.94 3.92 9.80 6.86 7.84 3.92 6.86  17.65;
      5.83 4.85 3.88 9.71 7.77 3.88 1.94 6.80 3.88 2.91 3.88 9.71 6.80 6.80 8.74  11.65;
      4.76 3.81 1.90  12.38 8.57 5.71 .00 6.67 5.71 3.81  10.48  10.48 3.81 8.57 9.52 2.86;
      3.88 2.91 2.91  10.68 5.83 .97 6.80 5.83 5.83 5.83 9.71 3.88 4.85 5.83  11.65  10.68;
      3.42 9.40 5.98 3.42  10.26 1.71 4.27  27.35 5.13 3.42 4.27 3.42 2.56 6.84 1.71 5.98;
      8.49 5.66 4.72 8.49 4.72 8.49 2.83 6.60  11.32 1.89 9.43 5.66 2.83 9.43 4.72 3.77;
      3.45 7.76 4.31 4.31  10.34 .86 3.45  27.59 1.72 6.03 8.62 3.45 4.31 5.17 1.72 6.03];
  ddd=[ 1.77 3.54 2.65 .88 .00 .00 7.96 .88 4.42 2.65  17.70  10.62 3.54 4.42 4.42 7.08 1.77 3.54  13.27 7.08;
      1.92 1.92 .96 .96 .00 .96 1.92 .96 4.81  12.50 7.69  11.54 8.65 3.85 3.85 6.73 9.62 6.73 7.69 2.88;
         .98 .00 .00 5.88 .98 8.82 2.94 .00 .00 2.94  10.78 5.88  13.73 .00 4.90 3.92  19.61 1.96 8.82 5.88;
         .00 .00 .00 .87 .00 .87  13.04 1.74 6.09 2.61  11.30  13.04 3.48 5.22 3.48 8.70 3.48 1.74  14.78 7.83;
      2.86 .00 .00 3.81 .95 3.81 3.81 .00 3.81 3.81 9.52 9.52  12.38 2.86 9.52 3.81 7.62 2.86 7.62 9.52;
         .00 .00 .88 2.63 .00 1.75  13.16 .88 4.39 1.75  14.04 9.65 7.02 5.26 4.39  11.40 2.63 1.75  10.53 6.14;
      1.92 .00 .00 2.88 .96 4.81 2.88 .00 1.92 4.81  12.50 6.73  13.46 1.92 6.73 4.81  10.58 3.85 9.62 7.69;
      2.56 3.42 .00 .85 .85 .85  12.82 .85 1.71 .85  20.51 2.56 3.42 9.40 5.98  11.11 .85 4.27  11.97 3.42;
         .00 .00 .00 2.97 2.97 9.90 2.97 .00 .99 3.96 6.93 1.98  13.86 1.98 2.97 3.96  23.76 2.97 8.91 6.93;
      1.87 .93 3.74 2.80 .00 .00 2.80 .00 7.48 8.41 9.35 7.48 3.74  14.95  12.15 .00 2.80 4.67 7.48 7.48;
         .00 .89 .00 .00 .00 1.79 8.04 .00 5.36 4.46  15.18 8.04 8.93 4.46 3.57 8.04 4.46 6.25  13.39 5.36;
      2.75 .00 .92 2.75 .92 3.67 4.59 3.67 3.67 1.83 9.17 5.50 3.67 5.50 6.42 7.34 8.26 5.50  11.01 9.17;
      1.80 .90 .90 .90 .00 .90 9.01 .00 3.60 7.21  14.41 8.11 7.21 6.31 7.21 4.50 1.80 7.21  11.71 4.50;
      2.94 .00 .00 5.88 .00 6.86 1.96 .00 3.92 6.86 3.92 9.80  13.73 .98 5.88 2.94  10.78 .98  10.78 9.80;
      2.91 1.94 2.91 1.94 .00 5.83 1.94 .00 1.94 9.71 5.83 8.74  10.68 1.94 3.88 3.88 8.74 2.91  11.65  10.68;
      2.86 .95 .00  11.43 1.90 1.90 2.86 .00 4.76 3.81 5.71 8.57 8.57 6.67 9.52 4.76 5.71 2.86 7.62 7.62;
      1.94 .97 1.94 4.85 1.94 3.88 1.94 .97 .97 6.80 4.85 8.74  10.68 2.91 6.80 2.91 9.71 6.80 8.74 7.77;
      1.71 .85 1.71 .85 .85 2.56  16.24 .85 1.71 .85  16.24 5.13 6.84 5.98 3.42  11.11 1.71 5.13  11.11 3.42;
         .94 .94 1.89 .94 .94 .94 1.89 .94  10.38 7.55 5.66 9.43 8.49 8.49 7.55 5.66 6.60  11.32 6.60 .94;
         .86 .86 .00 1.72 .86 .86  17.24 .86 2.59 1.72  15.52 7.76 5.17 3.45 4.31 9.48 5.17 5.17 9.48 5.17];
x=[ 29.73  17.12  13.51  39.64  43.24;
      27.03  16.22  15.32  41.44  42.34;
      27.03  21.62 6.31  45.05  33.33;
      42.34  10.81  28.83  18.02  71.17;
      23.42  23.42  10.81  42.34  34.23;
      35.14  12.61  12.61  39.64  47.75;
      35.14 9.91  18.92  36.04  54.05;
      27.93  16.22  18.92  36.94  46.85;
      20.72  20.72  15.32  43.24  36.04;
      18.18  27.27  13.64  40.91  31.82;;
      35.45 4.55  50.00  10.00  85.45;
      32.73 2.73  50.00  14.55  82.73;
      25.45  10.00  51.82  12.73  77.27;
      30.00 8.18  50.00  11.82  80.00;
      29.09 .00  64.55 6.36  93.64;
      36.36 8.18  46.36 9.09  82.73;
      35.45  24.55  26.36  13.64  61.82;
      29.09  11.82  50.00 9.09  79.09;
      21.82  14.55  56.36 7.27  78.18;
      20.00  17.27  56.36 6.36  76.36];
xx=[ 9.01 9.01 3.60 8.11 4.50 .90 4.50 3.60 3.60 3.60 1.80 8.11  11.71 2.70 5.41  18.92;
      9.91 7.21 3.60 5.41 2.70 1.80 5.41 5.41 4.50 1.80 .90 9.01 9.91 4.50 5.41  21.62;
      5.41  11.71 3.60 5.41 2.70 1.80 .90 .90 5.41 .90 .90  14.41  13.51 .90 7.21  23.42;
      18.92 5.41  11.71 5.41  10.81 1.80 5.41  10.81 5.41 1.80 .90 2.70 6.31 4.50 2.70 4.50;
      6.31 8.11 1.80 7.21 1.80 2.70 2.70 3.60 5.41 4.50 2.70  10.81 9.91 .90 9.01  21.62;
      15.32 2.70 6.31 9.91 3.60 1.80 1.80 5.41 4.50 .00 .00 8.11  10.81 .90 8.11  19.82;
      15.32 1.80  10.81 7.21 4.50 2.70 6.31 5.41 .90 1.80 .90 6.31  13.51 .90 4.50  16.22;
      8.11 3.60 6.31 9.91 5.41 3.60 2.70 7.21 2.70 3.60 1.80 8.11  10.81 1.80 7.21  16.22;
      9.01 .90 4.50 6.31 .00 3.60 7.21 4.50 3.60 2.70 2.70  11.71 7.21 3.60  13.51  18.02;
      6.36 3.64 1.82 6.36 1.82 5.45 2.73 3.64 5.45 3.64 4.55  13.64 4.55 3.64  13.64  18.18;
      15.45 2.73  14.55 2.73  16.36 .91 1.82  30.00 .91 .91 .91 1.82 2.73 4.55 .00 2.73;
      13.64 .91  10.91 6.36  15.45 1.82 1.82  30.91 .91 .91 .00 .91 2.73 7.27 .00 4.55;
      6.36 4.55  10.00 4.55  12.73 1.82 2.73  34.55 2.73 2.73 1.82 1.82 3.64 4.55 1.82 2.73;
      8.18 .91  12.73 7.27  13.64 6.36 1.82  28.18 2.73 4.55 .00 .91 5.45 4.55 .91 .91;
      13.64 .00  12.73 1.82  13.64 .00 2.73  48.18 .00 .00 .00 .00 1.82 3.64 .00 .91;
      16.36 3.64  15.45 .91  13.64 4.55 4.55  22.73 1.82 5.45 .00 .91 4.55 2.73 .00 1.82;
      17.27 5.45  10.91 1.82  10.00 6.36 4.55 5.45 4.55 7.27 9.09 2.73 3.64 2.73 3.64 3.64;
      8.18 7.27  11.82 1.82  15.45 1.82 .91  30.91 3.64 3.64 1.82 2.73 1.82 3.64 .91 2.73;
      2.73 2.73  13.64 1.82  14.55 9.09 .91  31.82 1.82 8.18 1.82 2.73 2.73 2.73 .91 .91;
      6.36 6.36 6.36 .91 9.09  10.00 3.64  32.73 2.73  13.64 .91 .00 1.82 3.64 .00 .91];
  xxx=[ 5.41 .90 2.70 .90 5.41 3.60 .90 1.80 2.70 8.11 4.50 1.80  25.23 3.60 3.60 5.41  13.51 .00 3.60 4.50;
      2.70 2.70 .00 .00 3.60 6.31 2.70 .90 7.21 7.21 6.31 1.80  18.92 .90 6.31 1.80  14.41 .00 3.60  10.81;
      2.70 2.70 2.70 .00 3.60 6.31 .00 .90 4.50 5.41 1.80 .90  29.73 .00 5.41 4.50  22.52 .00 1.80 2.70;
      15.32 6.31 .00 .00 .00 .90 9.01 1.80 6.31  10.81  12.61 3.60 4.50 1.80 2.70 5.41 1.80 1.80 7.21 6.31;
      3.60 1.80 2.70 .00 5.41 7.21 .90 .00 4.50 1.80 2.70 3.60  20.72 1.80 6.31 4.50  19.82 1.80 1.80 7.21;
      9.01 .90 .90 .00 2.70 5.41 4.50 .00 2.70  13.51 6.31 .00  25.23 .90 1.80 1.80  16.22 .00 2.70 3.60;
      9.01 1.80 .00 .00 1.80 4.50 4.50 .90 3.60  16.22 8.11 .00  17.12 2.70 1.80 1.80  10.81 .90 6.31 6.31;
      2.70 1.80 .90 .90 2.70 3.60 2.70 .90 4.50 9.91 8.11 3.60  18.92 .90 2.70 4.50  12.61 .90 7.21 8.11;
      5.41 .00 .90 1.80 5.41 9.01 1.80 .90 3.60 6.31 1.80 3.60  11.71 2.70 2.70 2.70  20.72 1.80 4.50  10.81;
      3.64 .91 2.73 6.36 3.64  10.91 .91 1.82 3.64 2.73 2.73 .91  17.27 .00 4.55 4.55  17.27 4.55 1.82 7.27;
      9.09 .91 .00 .00 .00 .00  24.55 .00 3.64 6.36  33.64 .91 4.55 1.82 .00 1.82 .00 2.73 5.45 2.73;
      2.73 .91 .00 .00 .00 .00  19.09 .00 1.82 8.18  37.27 .00 4.55 4.55 .00 2.73 .00 .91  10.00 5.45;
         .91 2.73 .00 .00 .00 .00  27.27 1.82 1.82 5.45  26.36 2.73 4.55 2.73 4.55 5.45 1.82 2.73 5.45 1.82;
      6.36 5.45 .00 .00 1.82 .00  20.00 5.45 2.73 2.73  24.55 .00 1.82 3.64 3.64 8.18 .91 .91 9.09 .91;
      11.82 .91 .00 .00 1.82 .00  47.27 1.82 .00 3.64  25.45 .00 .91 .91 .00 .00 .00 .00 2.73 .91;
      10.00 2.73 .91 .00 .00 .00  14.55 4.55 5.45 3.64  31.82 .91 .91 3.64 1.82 6.36 .00 .00 7.27 3.64;
      10.91 .91 3.64 3.64 .00 .91 8.18 2.73  12.73 9.09  11.82 3.64 3.64 6.36 1.82 1.82 6.36 6.36 1.82 1.82;
      4.55 4.55 .00 .00 .91 .91  21.82 .91 4.55 .91  29.09 .00 3.64 1.82 .91  10.91 2.73 4.55 4.55 .91;
      3.64 .91 1.82 .91 .91 .00  25.45 5.45 3.64 .00  21.82 1.82 1.82 3.64 .91  13.64 .91 2.73 5.45 2.73;
      2.73 .91 5.45 .00 .00 .00  23.64  10.00 6.36 1.82  13.64 .00 1.82 8.18 1.82  13.64 .00 1.82 6.36 .00];
   ffx=[x xx xxx];
   ffx=[ffx(1:16,

;ffx(18:20,

]
   ffd=[d dd ddd];
   cx=cov(ffx);
   [vx,ex]=eig(cx);
   ex1=eig(cx)
   e1=mean(ex1)*41;
   ex2=ex1(36:41,

;
   e2=mean(ex2)*6;
   e2/e1
      vx1=[vx(:,38:41)];
   s=ffx*vx1;ss=ffd*vx1;
   x=s(1:10,

;
y=s(11:19,

;
u1=mean(x);u2=mean(y);
u1-u2;
z=8/9*(cov(x)+cov(y));
ux=0.5*(u1-u2)*inv(z);
u12=0.5*u1+0.5*u2;
u0=ux*u12.';
la=0
for i=1:9
fd(i)=ux*ss(i,

.';
tx(i)=ux*x(i,

.';
fy(i)=ux*y(i,

.';
  if fd(i)>u0
   pbd(i)=1;
   la=la+1;
      else
         pbd(i)=2 ;
end
      if tx(i)>u0
      lbx(i)=1 ;
      else
         lbx(i)=2;
      end
      if fy(i)>u0
      lby(i)=1 ;
      else
         lby(i)=2 ;
      end
  for n=10:19
   fd(n)=ux*ss(n,