HOME > テクニカルサポート > SYSTAT > モンテカルロ例
ヒューリンクステクニカルサポート
SYSTAT に関する皆様からのご質問のうち、よくある質問を掲載しました。
モンテカルロ法の例
SYSTAT 製品ページ

TOP > モンテカルロ> 例

遺伝子頻度の推定

Rao (1973) には、得点法を使用して血液型 O、A、B の遺伝子頻度の最大尤度の推定を行う方法が示されています。 McLachlan and Krishnan (1997) では、同じ問題に対して EM アルゴリズムが使用されています。 この例では、ギブズ・サンプリングによる、これらの遺伝子頻度のベイズの推定を示します。 4つのセル頻度を持つ次の多項モデルを考えてみます。その確率はパラメータ pqr で、p+q+r=1 です。n = no+ nA+ nB+nAB とします。

データ モデル
no=176 r2
nA=182 p2+2pr
nB=60 q2+2qr
nAB=17 2pq

この問題のために仮説的に大きくしたデータを、多項モデル {n; (1-p-q)2, p2, 2p(1-p-q), q2, 2q(1-p-q), 2pq} nOnAAnAOnBBnBOnAB とします。 後者の完全モデルについては、nAAnBB は欠損値とみなすことができます。

モデル

X ~ Multinomial6 (435; (1-p-q)2, p2, 2p(1-p-q), q2, 2q(1-p-q), 2pq)

事前情報

(p, q, r) ~ Dirichlet (α, β, γ)

完全条件付き密度は次のようになります。

p および q からランダム標本を生成するには、ベータ分布から生成した値と (1-q) および (1-p) とをそれぞれ乗算する必要があります。 これをユーザーのシステムに実装することはできないため、次のようにします。

完全条件付き密度に出現する p および q はすべて、p は (1-q)p に、q は (1-p)q に置換します。α=2、β=2、γ=2 とし、次のように入力します。

FORMAT 10 5
MCMC
GIBBS / SIZE=10000 NSAMP=1 BURNIN=1000 GAP=1 RSEED=1783
FULLCOND / VAR='NAA' DIST=N PAR1='182' ,
PAR2='(((1-Q)*P)^2)/((((1-Q)*P)^2)+(2*((1-Q)*P)*(1-((1-Q)*P)-((1-P)*Q))))' INIT=40
FULLCOND / VAR='NBB' DIST=N PAR1='60',
PAR2='(((1-P)*Q)^2)/((((1-P)*Q)^2)+(2*((1-P)*Q)*(1-((1-P)*Q)-((1-Q)*P))))'INIT=5
FULLCOND / VAR='P' DIST=B PAR1='NAA+182+17+1' PAR2='(2*176)+182+60-NAA-NBB+1',
INIT=0.1
FULLCOND / VAR='Q' DIST=B PAR1='NBB+60+17+1' PAR2='(2*176)+182+60-NAA-NBB+1',
INIT=0.5
SAVE GIBBSGENETIC.SYD
GENERATE
USE GIBBSGENETIC.SYD
LET P=(1-Q1)*P1
LET Q=(1-P1)*Q1
LET R=1-P-Q
LET RBEP= (1-Q)*((NAA1+182+17+2)/((NAA1+182+17+2)+((2*176)+182+60-NAA1-NBB1+2)))
LET RBEQ=(1-P)*((NBB1+60+17+2)/((NBB1+60+17+2)+((2*176)+182+60-NAA1-NBB1+2)))
LET RBER=1-RBEP-RBEQ
STATS
CBSTAT P Q R RBEP RBEQ RBER/ MAXIMUM MEAN MEDIAN MINIMUM SD VARIANCE N ,
PTILE=2.5 50 97.5
DENSITY P RBEP/HIST XMIN=0.20 XMAX=0.35
DENSITY Q RBEQ/HIST XMIN=0.05 XMAX=0.13
DENSITY R RBER/HIST XMIN=0.60 XMAX=0.75
FORMAT

出力は次のようになります。

  P Q R RBEP RBEQ
N of cases 10000 10000 10000 10000 10000
Minimum 0.195 0.058 0.606 0.241 0.087
Maximum 0.309 0.136 0.718 0.293 0.109
Median 0.253 0.090 0.657 0.264 0.096
Mean 0.253 0.090 0.657 0.264 0.096
Standard Dev 0.016 0.010 0.014 0.007 0.003
Variance 0.000 0.000 0.000 0.000 0.000
Method = CLEVELAND          
2.5 % 0.223 0.072 0.628 0.251 0.090
50 % 0.253 0.090 0.657 0.264 0.096
97.5 % 0.284 0.110 0.685 0.278 0.102
 
  RBER        
N of cases 10000        
Minimum 0.612        
Maximum 0.665        
Median 0.640        
Mean 0.640        
Standard Dev 0.007        
Variance 0.000        
Method = CLEVELAND          
2.5 % 0.626        
50 % 0.640        
97.5 % 0.653        

得点法または EM 法で評価した p、q、r の最尤推定値は 0.26444、0.09317、0.64239 です。 使用可能な事後情報を使用すると、p、q、r の推定値はギブズ・サンプリングによって近似されます。 p、q、r の経験的推定値は 0.25294、0.09041、0.65665 です。 Rao-Blackwell 化した推定値は、それぞれ 0.26416、0.09590、0.63994 です。

前のページにもどる