GNU Octave

自作 m-file

特に家畜育種分野で頻繁に利用する演算を，GNU Octaveのm-file（関数ファイル）として作成しました。利用には制限がありません。再配布・改変も自由に行っていただいて結構です。中にはパフォーマンスが悪いものもあるかと思います。ご指摘いただければ幸いです。

事情により，拡張子を.txtにしています。ダウンロード後，もとの.mに変更してください。

血縁行列関係

additiverel.m - 父と母のリストから分子血縁行列（相加的血縁行列）を計算する [解説]
genomerel.m - 父と母のリストから配偶子血縁行列を計算する [解説]
dominancerel.m - 配偶子血縁行列から優性血縁行列を計算する [解説]
parentmatrix.m - 父と母のリストから両親の位置を表す行列を生成する。 [解説]
mendelianvar.m - 父と母のリストと近交係数からメンデリアンサンプリング分散の期待値ベクトルを計算する。 [解説]
raij.m - 父と母のリストから分子血縁行列のi行j列めの要素を再帰的に計算する。 [解説]
inbcoef.m - 父と母のリストから近交係数を計算する。[解説]

その他の計算

dsum.m - 直和（direct sum）を計算する [解説]
v2r.m - 分散・共分散行列から相関行列を計算する [解説]
stdleg.m - 正規化して標準化したルジャンドル多項式の係数を含む下三角行列を返す [解説]

血縁行列関係のm-file

以下，Henderson（1976）の血縁関係を例にする。

additiverel.m - 父と母のリストから分子血縁行列（相加的血縁行列）を計算する

>> s=[0 0 1 1 3 1 5];
>> d=[0 0 0 2 4 4 6];
>> A=additiverel(s,d)
A =

  1.00000  0.00000  0.50000  0.50000  0.50000  0.75000  0.62500
  0.00000  1.00000  0.00000  0.50000  0.25000  0.25000  0.25000
  0.50000  0.00000  1.00000  0.25000  0.62500  0.37500  0.50000
  0.50000  0.50000  0.25000  1.00000  0.62500  0.75000  0.68750
  0.50000  0.25000  0.62500  0.62500  1.12500  0.56250  0.84375
  0.75000  0.25000  0.37500  0.75000  0.56250  1.25000  0.90625
  0.62500  0.25000  0.50000  0.68750  0.84375  0.90625  1.28125

各個体の近交係数は，それぞれの対角要素から1を減じた値である。

>> f=diag(A)-1
f =

  0.00000
  0.00000
  0.00000
  0.00000
  0.12500
  0.25000
  0.28125

genomerel.m - 父と母のリストから配偶子血縁行列を計算する

>> s=[0 0 1 1 3 1 5];
>> d=[0 0 0 2 4 4 6];
>> G=genomerel(s,d)
G =

  1.00000  0.00000  0.00000  0.00000  0.50000  0.00000  0.50000  0.00000  0.25000  0.25000  0.50000  0.25000  0.25000  0.37500
  0.00000  1.00000  0.00000  0.00000  0.50000  0.00000  0.50000  0.00000  0.25000  0.25000  0.50000  0.25000  0.25000  0.37500
  0.00000  0.00000  1.00000  0.00000  0.00000  0.00000  0.00000  0.50000  0.00000  0.25000  0.00000  0.25000  0.12500  0.12500
  0.00000  0.00000  0.00000  1.00000  0.00000  0.00000  0.00000  0.50000  0.00000  0.25000  0.00000  0.25000  0.12500  0.12500
  0.50000  0.50000  0.00000  0.00000  1.00000  0.00000  0.50000  0.00000  0.50000  0.25000  0.50000  0.25000  0.37500  0.37500
  0.00000  0.00000  0.00000  0.00000  0.00000  1.00000  0.00000  0.00000  0.50000  0.00000  0.00000  0.00000  0.25000  0.00000
  0.50000  0.50000  0.00000  0.00000  0.50000  0.00000  1.00000  0.00000  0.25000  0.50000  0.50000  0.50000  0.37500  0.50000
  0.00000  0.00000  0.50000  0.50000  0.00000  0.00000  0.00000  1.00000  0.00000  0.50000  0.00000  0.50000  0.25000  0.25000
  0.25000  0.25000  0.00000  0.00000  0.50000  0.50000  0.25000  0.00000  1.00000  0.12500  0.25000  0.12500  0.56250  0.18750
  0.25000  0.25000  0.25000  0.25000  0.25000  0.00000  0.50000  0.50000  0.12500  1.00000  0.25000  0.50000  0.56250  0.37500
  0.50000  0.50000  0.00000  0.00000  0.50000  0.00000  0.50000  0.00000  0.25000  0.25000  1.00000  0.25000  0.25000  0.62500
  0.25000  0.25000  0.25000  0.25000  0.25000  0.00000  0.50000  0.50000  0.12500  0.50000  0.25000  1.00000  0.31250  0.62500
  0.25000  0.25000  0.12500  0.12500  0.37500  0.25000  0.37500  0.25000  0.56250  0.56250  0.25000  0.31250  1.00000  0.28125
  0.37500  0.37500  0.12500  0.12500  0.37500  0.00000  0.50000  0.25000  0.18750  0.37500  0.62500  0.62500  0.28125  1.00000

dominancerel.m - 配偶子血縁行列から優性血縁行列を計算する

>> s=[0 0 1 1 3 1 5];
>> d=[0 0 0 2 4 4 6];
>> G=genomerel(s,d);
>> D=dominancerel(G)
D =

  1.00000  0.00000  0.00000  0.00000  0.12500  0.25000  0.18750
  0.00000  1.00000  0.00000  0.00000  0.00000  0.00000  0.03125
  0.00000  0.00000  1.00000  0.00000  0.12500  0.00000  0.09375
  0.00000  0.00000  0.00000  1.00000  0.12500  0.25000  0.21875
  0.12500  0.00000  0.12500  0.12500  1.01562  0.15625  0.31641
  0.25000  0.00000  0.00000  0.25000  0.15625  1.06250  0.35156
  0.18750  0.03125  0.09375  0.21875  0.31641  0.35156  1.07910

parentmatrix.m - 父と母のリストから両親の位置を表す行列を生成する。

>> s=[0 0 1 1 3 1 5];
>> d=[0 0 0 2 4 4 6];
>> P=parentmatrix(s,d)
P =

   0   0   0   0   0   0   0
   0   0   0   0   0   0   0
   1   0   0   0   0   0   0
   1   1   0   0   0   0   0
   0   0   1   1   0   0   0
   1   0   0   1   0   0   0
   0   0   0   0   1   1   0

mendelianvar.m - 父と母のリストと近交係数からメンデリアンサンプリング分散の期待値ベクトルを計算する。

>> s=[0 0 1 1 3 1 5];
>> d=[0 0 0 2 4 4 6];
>> f=diag(additiverel(s,d))-1;
>> m=mendelianvar(s,d,f)
m =

  1.00000
  1.00000
  0.75000
  0.50000
  0.50000
  0.50000
  0.40625

以下のようにして，分子血縁行列の逆行列を生成することができる。ただし，この小さな例ではinv(A)を実行した方が速い。

>> s=[0 0 1 1 3 1 5];
>> d=[0 0 0 2 4 4 6];
>> n=length(s);
>> P=parentmatrix(s,d);
>> A=additiverel(s,d);
>> f=diag(A)-1;
>> m=mendelianvar(s,d,f);
>> inva=(eye(n)-0.5*P')*diag(1./m)*(eye(n)-0.5*P)
inva =

   2.33333   0.50000  -0.66667  -0.50000   0.00000  -1.00000   0.00000
   0.50000   1.50000   0.00000  -1.00000   0.00000   0.00000   0.00000
  -0.66667   0.00000   1.83333   0.50000  -1.00000   0.00000   0.00000
  -0.50000  -1.00000   0.50000   3.00000  -1.00000  -1.00000   0.00000
   0.00000   0.00000  -1.00000  -1.00000   2.61538   0.61538  -1.23077
  -1.00000   0.00000   0.00000  -1.00000   0.61538   2.61538  -1.23077
   0.00000   0.00000   0.00000   0.00000  -1.23077  -1.23077   2.46154

>> inva-inv(A)
ans =

  1.0e-16  *

   0.00000   0.00000   0.00000   0.00000   0.00000   0.00000   0.00000
   0.00000   0.00000   0.00000   0.00000   0.00000   0.00000   0.00000
  -1.11022   0.00000   0.00000   0.00000   0.00000   0.00000   0.00000
   0.00000   0.00000   0.00000   0.00000   0.00000   0.00000   0.00000
   0.00000   0.00000   0.00000   0.00000   0.00000   0.00000   0.00000
  -1.11022   0.00000   0.00000   0.00000   0.00000   0.00000   0.00000
   0.00000   0.00000   0.00000   0.00000   0.00000   0.00000   0.00000

raij.m - 父と母のリストから分子血縁行列のi行j列めの要素を再帰的に計算する。
```
>> s=[0 0 1 1 3 1 5];
>> d=[0 0 0 2 4 4 6];
>> raij(1,5,s,d)
ans = 0.50000
>> raij(7,7,s,d)
ans = 1.2812
```
additiverel とは異なり，行列を構築せずにその要素だけを得る。さかのぼる世代が増えると遅くなる。
inbcoef.m - 父と母のリストから近交係数を計算する。
```
>> s=[0 0 1 1 3 1 5];
>> d=[0 0 0 2 4 4 6];
>> f=inbcoef(s,d)
f =

  0.00000
  0.00000
  0.00000
  0.00000
  0.12500
  0.25000
  0.28125
```
計算にはMeuwissen and Luo（1992）のアルゴリズムを使用している。

その他計算向けのm-file

dsum.m - 直和（direct sum）を計算する。引数はいくつでもよい。

>> X=[1 2]';
>> Y=ones(3,1)*5;
>> Z=eye(3);
>> dsum(X,Y,Z)
ans =

    1    0    0    0    0
    2    0    0    0    0
    0    5    0    0    0
    0    5    0    0    0
    0    5    0    0    0
    0    0    1    0    0
    0    0    0    1    0
    0    0    0    0    1

v2r.m - 分散・共分散行列から相関行列を計算する。

>> V=[30 12 4        
      12 32 8
       4  8 40];
>> v2r(V)
ans =

  1.00000  0.38730  0.11547
  0.38730  1.00000  0.22361
  0.11547  0.22361  1.00000

手抜きしているので，何もエラーチェックをしていない。

stdleg.m - 正規化して標準化したルジャンドル多項式の係数を含む下三角行列を返す。

>> L=stdleg(6)
L =

  1.0e+01  *

   0.10000   0.00000   0.00000   0.00000   0.00000   0.00000
   0.00000   0.17321   0.00000   0.00000   0.00000   0.00000
  -0.11180   0.00000   0.33541   0.00000   0.00000   0.00000
   0.00000  -0.39686   0.00000   0.66144   0.00000   0.00000
   0.11250   0.00000  -1.12500   0.00000   1.31250   0.00000
   0.00000   0.62187   0.00000  -2.90205   0.00000   2.61184

これは，正規化したルジャンドル多項式にルート2（切片を1にするための係数）を乗じたものである。

戻る