抄録(abstract)
最尤法の考え方を視覚的、体感的に理解するために統計ソフトRを使って簡単な計算例で実験してみた。 To gain a
visual and intuitive understanding of maximum likelihood
estimation, I ran a simple example in R
最尤法をよく理解するためには統計学の専門書を熟読する以外に方法はないのかもしれないが、試しに統計ソフトRを使って簡単な計算例で計算実験しつつ体感的に最尤法の考え方を探ってみた。稚拙な感想文のようなものだが、ざっと読み流していただければと思う。（Rコードは文末に掲載）
よく教科書でベルヌーイ分布のような離散型確率分布で尤度を説明することがある。しかし、計算例などを見ても最尤法の仕組みの話につながらないことが多く直感的に理解することも難しかった。
尤度関数は、母集団の未知の母数（パラメータ）を仮定した下で、無作為に独立に、観測された標本データ（標本要素）\(x_{1},\ldots\) ,\(x_{n}\)が得られる同時確率密度を、母数の関数とみなしたものである。
数式で示すと
n個の標本要素を \(x_{1}\),\(x_{2}\),\(\ldots\) ,\(x_{n}\)とし、母数のもとでの確率密度関数をf(x|)
とすると、尤度関数 \(L(\theta|x_{1}\),\(\ldots\) ,\(x_{n}\))は、以下の積の形で表される。
\(L(\theta|x_{1},\ldots,x_{n})=\prod^{n}_{i}f(x_{i}|\theta)=f(x_{1}|\theta)\cdot f(x_{2}|\theta)\cdot\cdot\cdot f(x_{n}|\theta)\)
尤度は離散型確率分布では尤度は確率（同時確率）と同じになるが、尤もらしさを示す指標を意味している。この同時確率密度\(L(\theta|x_{1},\ldots,x_{n})\)は標本観測値\(x_{1},\ldots,x_{n}\)に依存していると同時に母集団のパラメータにも依存している。標本観測値xiは既に標本抽出で得た確定値なので尤度関数はパラメータの関数とも解釈できる。そこで最尤法はを動かして尤度Lが最大となるような統計量を求めようとする。言い換えると、標本観測値を既に与えられたものとしてを変化させ\(L(\theta|x_{1},\ldots,x_{10})\) が最大となるようなを最尤推定量とする考え方である。
100個の電球が入った箱には良品と不良品が含まれているが、その割合は不明である。そこで、無作為に3個を抽出したところ、（良品、良品、不良品）となり1個が不良品だった。良品には0、不良品には1を与えて、この箱の不良品率ｐを推定したい。電球の状態は良品か不良品かの何れかであり、これはベルヌーイ分布に従うと考えられる。この例ではサンプル数が３個といかにも少なく常識的にも信頼性に疑問符がつくが、話を単純化するために３個としている。
ベルヌーイ分布の確率関数(probability mass function)は\(f(x;p)=p^{x}(1-p)^{1-x}\) \(\,\,\,x=0,1\) で表せる。目的とする事象は出現する(成功）か、出現しない（失敗）かの二者択一なので、事象が出現した（成功の）場合には1、出現しない（失敗の）場合には0をとる。パラメータｐは不良品率としているので不良品が見つかれば成功で1、良品であれば失敗で0とする。
3個のサンプル調査で（良品、良品、不良品）の標本が実現したとすれば\(x=(0,0,1)\)
となる。
１番目のサンプル（標本要素）は良品なのでｘ＝0
その確率はf(x;p)＝\(p^{0}(1-p)^{1}-0=1-p\)
2番目のサンプル（標本要素）も良品でｘ＝0
その確率はf(x;p)＝\(p^{0}(1-p)^{1}-0=1-p\)
3番目のサンプルは不良品なのでｘ＝1
その確率はf(x;p)＝\(p^{1}(1-p)^{1-1}=p\)
先ほどの定義で
尤度関数 \(L(\theta|x_{1},\ldots,x_{n}\))は、以下の積の形で表される。
\(L(\theta|x_{1},\ldots,x_{n})=\prod f(x_{i}|\theta)=f(x_{1}|\theta)\cdot f(x_{2}|\theta)\cdot\cdot\cdot f(x_{n}|\theta)\)
従って、3つの事象が同時に出現する確率は尤度関数は
\(L(p)=p\cdot(1-p)^{2}\) で表され、\(f^{\prime\prime}\)
離散型確率モデルではｘ＝（0，0，1）の同時確率で表される。
最尤法はこの同時確率（連続型確率モデルでは同時確率密度）が最大になるようなｐを求めて、これを最尤推定値とする。この尤度関数は\(L\left(p\right)=\left(1-2p+p^{2}\right)p=p-2p^{2}+p^{3}\)で表されるので、この関数をｐで微分して最大値となるｐを\(L^{\prime}\left(p\right)=1-4p+3p^{2}=0\) として求める。これを解くと\((3p-1)(p-1)=0\) よって\(p=(1/3,1\)) を得る。尤度関数L(p)にこれを代入すると\(L(1/3)=0.148..\,\,\)\(L(1)=0\) となるので\(L(1/3)=0.148\)という最大尤度を与える\(\text{ｐ＝1/3}\)をパラメータｐの最尤推定値とする。( ２次導関数で\(L^{\prime\prime}\left(p\right)=\left(1-4p+3p^{2}\right)^{\prime\prime}\)\(=6p-4\)からp=1/3の時の\(L^{\prime\prime}\left(p\right)<0\)となるので\(p=1/3\)で最大と確認)
統計ソフトRを使って、上記のベルヌーイ分布の最尤推定値を求めた結果をグラフに表してみる。
p=1/3 のところで尤度が0.148と最大になっている。最尤法は尤度関数 \(L(p)=p\cdot(1-p)^{2}\) についてｐを0から1まで変化させて最大尤度0.148を与えるｐ＝1/3をその推定値としている。観測された結果（良品、良品、不良品の同時確率）が最も起こりやすくするようなｐを最尤推定値となっている。最尤法は観測データに基づき、その観測結果が最も起こりやすいパラメータを推定している。この例では観測された結果（良品、良品、不良品）が最も起こりやすくなるようなパラメータを推定している。離散型確率分布では尤度はサンプル（標本要素）の同時確率となるので下限はゼロで上限は1となる。このモデルでは尤度関数を微分すると解析的にきれいな形の推定量が導出できるので、得られた最尤推定量の算式に観測値を入れて推定値が計算できる。これについては後で触れることにする。
この最尤法の考え方は正規分布のような連続型確率モデルでも同じであるが、正規分布のような連続型確率分布では尤度はサンプル（標本要素）の確率密度関数の積で定義される。正規分布では確率密度の積は、下限はゼロになるが上限は1という保証はなく無限大となる。例えば正規分布で分散が限りなくゼロに近い分布の確率密度は極端に尖がり（発散し）、これらの積は極端に大きな値となりうる。
つまり上限はなく+\(\infty\)の可能性がある。離散型確率分布のモデルでは同時確率になるので上限は1である点と大きく異なってくる。
正規分布の場合、尤度関数は密度関数の積で表現されているので、各変数ｘを－\(\infty\)から+\(\infty\)までの変域で重積分すれば1となるように思う人もいるかもしれないが、尤度関数ではそのような計算はできない。尤度関数に出てくる変数ｘは標本実現値であり確定値なので動かすことはできない。尤度関数で動かせる変数は母数（パラメータ）であり、例えば、ベルヌーイ分布であればｐ、正規分布であれば平均とか標準偏差といったパラメータが変数となる。
話を単純化するために、無作為、独立に２個のサンプルを正規母集団から抽出したとする。母集団は、平均は未知、標準偏差10で既知、の正規分布に従うと仮定する。標本要素２個からなる一組のサンプル（37，67）が得られた。
尤度関数はパラメータ（平均）の関数となり
\(L\left(\theta|x_{1},x_{2}\right)=\prod^{2}_{i}f\left(x_{i}|\theta\right)\)
\(=f\left(x_{1}=37|\theta\right)\times f\left(x_{2}=67|\theta\right)\)
対数尤度では
\(log(L(\theta|x_{1},x_{2}))=log(\prod f(x_{i}|\theta))\)
\(=log\left(f(x_{1}=37|\theta)\right)+log\left(f(x_{2}=67|\theta)\right)\)
で表される。ここでを動かして対数尤度が最大となるを探索する。それが最尤推定値となる。ベルヌーイ分布の事例と同じ考え方である。
Rでは正規分布の確率密度はdnorm(x , mean = mu, sd = sigma)で計算できる。例えば、正規分布で x=37 平均=50 標準偏差=10の確率密度関数はdnorm(37 , mean = 50､ sd = 10)で0.01713686と計算できる。
同様にして、正規分布で x=67 平均=50 標準偏差=10の確率密度関数はdnorm(67 , mean = 50､ sd = 10)で0.009404908と計算できる。
従って
\(log(L(\theta|x_{1},x_{2}))=log(\prod^{2}_{i}f(x_{i}|\theta))=\)
\(=log\left(f(x_{1}=37|\theta)\right)+log\left(f(x_{2}=67|\theta)\right)\)
から、一つの対数尤度（平均＝50，標準偏差＝10の場合）＝
log(0.01713686)+log(0.009404908)=－8.733047
を得る。この事例では標準偏差は10と既知のため固定されているが、平均は未知のため50以外の数値をとりうる。そこで平均muを0.1から80まで0.1刻みで変化させ、800個の平均について、対数尤度を800個求め、その対数尤度のうち最大値を与えるmuを平均の最尤推定値とする。Rで実験すると800の対数尤度のうち520番目の対数尤度が-8.693047と最大となり、それに対応する平均は52と計算される。この値は標本平均（37+67）/2＝52 と等しくなっている。
グラフでＲコードでの実験結果を示すと以下のようになる
グラフでは読み取りにくいが、赤色の垂直線で対数尤度が-8.693047と最大なり、平均は52となっている。
無作為に２個のサンプルを抽出したとする。母集団は、平均、標準偏差ともに未知の正規分布に従うと仮定する。
標本要素２個からなる一組のサンプル（37，67）が得られた。
この場合も同様に
\(log(L(\theta|x_{1},x_{2}))=log(\prod^{2}_{i}f(x_{i}|\theta))\ =\)
\(=log\left(f(x_{1}=37|\theta)\right)+log\left(f(x_{2}=67|\theta)\right)\)
のを変化させ最大尤度となるを最尤推定値とする。ただし、平均と標準偏差の2つのパラメータが未知なのでは（平均mu ,標準偏差 sd）の２要素を持つベクトルとなる。平均mu ,標準偏差 sdを同時に変化させて最大尤度となるを最尤推定値とする。グラフにすると尤度曲面が描かれ、その最大点に対応する平均mu ,標準偏差 sdが最尤推定値となる。
この事例では標準偏差と平均の２つのパラメータが未知である。そこで平均muを0.1から80まで0.1刻みで変化させ、同時に標準偏差sdを0.1から16まで0.1刻みで変化させる。800個の平均と160個の標準偏差について、対数尤度を128,000（=800*160）個求め、その対数尤度のうち最大値を与える平均と標準偏差のセットも求める。Rで計算実験そしてみるとサンプル数が２個と少ないので綺麗な３次元グラフにならないが赤点の個所が最大尤度となる。
平均mu＝52、標準偏差sd＝15の点で対数尤度は-8.253977と最大となる。最尤法ではこれを最尤推定値とする。
最尤法の考え方を直感的に理解しやすくするために、素朴な総当たり的なグリッドサーチで推定値を計算しているが、当然に最尤法のソフトでは効率的なアルゴリズムが開発されている。ニュートン法やBFGS法などはよく知られている。最尤法は１回で得られたサンプルから最大限の情報を引き出そうとする考え方で実務的にも広く使われている。反復実験を重視する立場では１回限りのデータに依存することを心配するが、経済事象などでは反復実験は不可能なのでブートストラップ法を併用して信頼性を補強する試みも行われている。
対数尤度関数をパラメータで微分して最尤推定量が解析的にきれいな数式(closed form)で得られれば、その推定量にサンプル値を入れてパラメータの推定値が計算できる。ベルヌーイ分布に従うパラメータｐや正規分布に従うパラメータ（平均と標準偏差）は解析的にきれいな数式で最尤推定量が得られる。その計算過程を以下に示す。
Xはベルヌーイ分布に従う確率変数とする。
ベルヌーイ分布の確率関数（probability mass function）は
\(f\left(x;p\right)=\left\{ p^{x}\left(1-p\right)^{1-x}\right\}\) ｘ＝0，1
pは未知の母数
尤度関数\(L(p)\)は
\(L(p)=\ p^{x1}\left(1-p\right)^{1-x1}p^{x2}\left(1-p\right)^{1-x2}\cdot\cdot\cdot p^{xn}\left(1-p\right)^{1-xn}\)
\(=p\left(1-p\right)^{n-\sum^{n}_{i=1}xi}\)
対数尤度関数\(log(L(p)\))は
\(log(L(p))=\ \left(\sum^{n}_{i=1}xi\right)\cdot log(p)+\left(n-\sum^{n}_{i=1}xi\right)\cdot\log(1-p)\)
\(\frac{dlog\left(L\left(p\right)\right)}{dp}=\frac{\sum^{n}_{i=1}xi}{p}+\frac{-1\cdot n}{1-p}-\frac{-1\cdot\left(-\sum^{n}_{i=1}xi\right)}{1-p}\)
\(=\frac{\sum^{n}_{i=1}xi}{p}-\frac{\left(n-\sum^{n}_{i=1}xi\right)}{1-p}=0\)
\(\widehat{p}=\frac{1}{n}\sum^{n}_{i=1}xi\)
正規分布 平均 分散\(\sigma^{2}\)(標準偏差\(\sigma\))
尤度関数は
\(L\left(\theta\right)=\prod^{n}_{i=1}\frac{1}{\sigma\sqrt{2\pi}}e^{\left\{ -\frac{\left(x_{i}-\mu\right)^{2}}{2{\sigma}^{2}}\right\} }\)
\(L\left(\theta\right)=\left(\frac{1}{2\pi\sigma^{2}}\right)^{n/2}e^{\left\{ -\frac{1}{2\sigma^{2}}\sum^{n}_{i=1}\left(x_{i}-\mu\right)^{2}\right\} }\)
両辺の対数をとると
\(log\left[L\left(\theta\right)\right]=log\left[\left(\frac{1}{2\pi\sigma^{2}}\right)^{n/2}\right]+log\left[e^{\left\{ -\frac{1}{2\sigma^{2}}\sum^{n}_{i=1}\left(x_{i}-\mu\right)^{2}\right\} }\right]\)
対数の法則を使って簡略化すると
\(log\left[\left(\frac{1}{2\pi\sigma^{2}}\right)^{n/2}\right]=\frac{n}{2}log(1)-\frac{n}{2}\left(log(2\pi)+log(\sigma^{2})\right)\)
\(=0-\frac{n}{2}log(2\pi)-\frac{n}{2}log(\sigma^{2})\)
\(log\left[e^{\left\{ -\frac{1}{2{\sigma}^{2}}\sum^{n}_{i=1}\left(x_{i}-\mu\right)^{2}\right\} }\right]=-\frac{1}{2\sigma^{2}}\sum^{n}_{i=1}\left(x_{i}-\mu\right)^{2}\)
従って
\(log\left[L\left(\theta\right)\right]=-\frac{n}{2}log(2\pi)-\frac{n}{2}log({\sigma}^{2})\ -\frac{1}{2\sigma^{2}}\sum^{n}_{i=1}\left(x_{i}-\mu\right)^{2}\)
合成関数の微分で
\(\frac{\partial log\left[L\left(\theta\right)\right]}{\partial\mu}=-0-0-\frac{1}{2{\sigma}^{2}}\cdot2\cdot(-1)\cdot\sum^{n}_{i=1}\left(x_{i}-\mu\right)^{\ }=\frac{1}{\sigma^{2}}\sum^{n}_{i=1}\left(x_{i}-\mu\right)^{\ }\)
商の微分を使うが変数は\(\sigma^{2}\) なので表記が紛らわしくなるが
\(\frac{\partial log\left[L\left(\theta\right)\right]}{\partial\sigma^{2}}=-\frac{n}{2\sigma^{2}}+\frac{1}{2\left(\sigma^{2}\right)^{2}}\sum^{n}_{i=1}\left(x_{i}-\mu\right)^{2}\)
平均の最尤推定量は対数尤度関数をで偏微分して
得た式がゼロになるを求める。
\(\frac{\partial log\left[L\left(\theta\right)\right]}{\partial\mu}=\frac{1}{{\sigma}^{2}}\sum^{n}_{i=1}\left(x_{i}-\mu\right)^{\ }=0\)
\(\frac{1}{\sigma^{2}}\sum^{n}_{i=1}\left(x_{i}-\mu\right)^{\ }=0\)
\(\sum^{n}_{i=1}x^{\ }_{i}-n\cdot\mu=0\)
\(\sum^{n}_{i=1}x^{\ }_{i}=n\cdot\mu\)
の最尤推定量をxとすれば
\(\bar{x}=\frac{\sum^{n}_{i=1}x^{\ }_{i}}{n}\)
分散の最尤推定量は対数尤度関数を\(\sigma^{2}\)で偏微分して得た
式にの代わりに\(\bar{x}\)を代入して得た式がゼロになる\(\sigma^{2}\)を求める
\(-\frac{n}{2\sigma^{2}}+\frac{1}{2\left(\sigma^{2}\right)^{2}}\sum^{n}_{i=1}\left(x_{i}-\bar{x}\right)^{2}=0\)
\(\hat{\sigma^{2}}=\frac{1}{n}\sum^{n}_{i=1}\left(x_{i}-\bar{x}\right)^{2}\)
算式から明らかなように分散の最尤推定量は不偏推定量でないが、他に不変性、漸近正規性などの望ましい性質を備えていることが知られている。
平均の最尤推定量を求めたのと同じ要領で、線形の単純回帰モデルの回帰係数を最尤法で推定するには
\(L\left(\theta\right)=\prod^{n}_{i=1}\frac{1}{\sigma\sqrt{2\pi}}exp\left\{ -\frac{\left(y_{i}-\beta-\beta_{1}x_{1}\right)^{2}}{2{\sigma}^{2}}\right\}\)
と尤度関数を定義すれば、この対数尤度関数をそれぞれの係数で偏微分して正規方程式が導出できる。
対数尤度関数を微分して解析的にきれいな数式(closed form)で推定統計量が得られれば、その式にサンプル値を入れてパラメータの推定値が計算できる。しかし複雑なモデルになると尤度関数は解析的に最大化できないために数値計算で最適値を探索することになる。そのためBFGS, Nelder-Mead, L-BFGS-Bなどのアルゴリズムが使われている。
AIが普及し始めてからコーディングの仕方も大きく変わり、RStudioもCopilotと連係して便利であった。 Copilotが有効化されているとプログラムを書き始めると影法師のように薄い色の文字(ghost text)が表示されコードの作成を先読みしてコードの提案をしてくれた。Copilotの提案が気に入ればTabキーを押すと確定した文字に変わり、そのままコードとして使えるのでプログラムが作成しやすかった。
しかし、このＡＩ機能は便利だと思っていたが、2026年6月頃からRStudio2026.04ではＡＩ機能は有料化されたようだ。ＡＩ機能はPosit AIとGitHubCopilotが選択できるが月額料金を支払う必要があるようだ。念のため自分のGitHub アカウントにログインしてみたらGitHub Copilot Free と表示されていたので月ごとのコード補完２０００回まで無料で 使えるようだ。しかし自分の何かの操作ミスで請求書でも通知されるような事態になることは避けたいのでRStudioのAssistantでChatもCode SuggestionsもNoneに設定している。ちょっとしたコードの書き方が分からない場合は無料ＡＩでも丁寧に教えてくれるので自分のプログラミングのレベルでは不自由を感じていない。世の流れとして今まで無料で使えた物が従量課金に変わっていくように感じる。
最初に１か月だけ有料AIをお試しに使ってみて満足がいくものであれば有料バージョンに切り換えようかと思っていたが、まだ無料ＡＩを使い続けている。有料バージョンだと、例えば、ピタゴラスの定理の証明を質問すると幾何、代数、解析学などで多種類の証明を提示できるそうだ。もし、そのような回答が必要になれば有料ＡＩに切換える時期かもしれない。しかし、最近、標準正規分布関数（累積分布）と誤差関数(Erf function)の変換公式の証明を置換積分を使って、無料AIが丁寧に解説してくれたので、自分にとって何時、有料バージョンに切り換えるのか全く想像もできないでいる。
ベルヌーイ分布の最尤推定
# ベルヌーイ分布の最尤推定
# 3個の製品を検査し、不良品が1個
# このときの不良品率の最尤推定値を ##
x<-seq(0,1,by=0.001)
bernlike<-sapply(x,function(p) (1-p)*(1-p)*p)
maxpvalue<-x[which.max(bernlike)]
maxpvalue
maxlike<-bernlike[which.max(bernlike)]
maxlike
plot(x,bernlike,type="l",lwd=2,
xlab="p",ylab="尤度")
abline(v=maxpvalue,col="red",lwd=2)
text(x=0.6,y=max(bernlike),
labels=paste("最尤推定値=",round(maxpvalue,3),
" 最尤尤度=",round(maxlike,3)),cex=0.4)
maxpvalue
尖った正規分布
#尖った正規分布
## 極端に尖った正規分布 平均値=30,標準偏差=0.0001の確率密度は
curve(dnorm(x,mean=30,sd=0.0001),from=29,to=31,lwd=2,col="blue",
xlab="x",ylab="確率密度",
main="平均値=30,標準偏差=0.0001の正規分布",cex.main=0.8)
points(x=30,y=dnorm(x=30,mean=30,sd=0.0001),col="red",pch=19)
text(x=29.4,y=dnorm(x=30,mean=30,sd=0.0001),
labels=paste("x=30の確率密度=",round(dnorm(x=30,mean=30,sd=0.0001),4)),
cex=0.4)
正規分布で x=37 平均=50 標準偏差=10の確率密度関数のグラフ
#正規分布で x=37 平均=50 標準偏差=10の確率密度関数のグラフ
curve(dnorm(x,mean=50,sd=10),from=0,to=100,lwd=3,col="blue",
xlab="x",ylab="確率密度",
main="平均値=50,標準偏差=10の正規分布",cex.main=0.8)
points(x=37,y=dnorm(x=37,mean=50,sd=10),col="red",pch=19)
text(x=20,y=dnorm(x=37,mean=50,sd=10),
labels=paste("x=37の確率密度=",round(dnorm(x=37,mean=50,sd=10),4)),
cex=0.5)
abline(v=37,col="cyan",lwd=1)
正規分布で 平均の推定( 標準偏差は10と既知 )
# 正規分布で 平均の推定
# === サンプルデータ ===
sampledata <- c(37,67)
# 母数 平均の推定 標準偏差は10と既知
# 尤度関数を定義 ===
# 平均と標準偏差のペア(mu, sd)の尤度を計算する関数
calclikelihood <- function(mu, sd, sample) {
densities <- dnorm(x = sample, mean = mu, sd = sd, log = TRUE)
likelihood <- sum(densities)
return(likelihood)
}
# 平均値の候補を0.1から80まで0.1刻みで生成
muvalue <- seq(from =0.1, to = 80, by = 0.1)
length(muvalue)
# for loopを使って尤度データを計算 # sd=10と既知の場合
loglikelihooddata1 <- matrix(0,nrow=length(muvalue),ncol=1)
for (i in 1:length(muvalue)) {
loglikelihooddata1[i, 1] <- calclikelihood(mu = muvalue[i],
sd = 10,
sample = sampledata)
}
length(loglikelihooddata1)
# which()関数で最大値の位置を特定
maxindex <- which(loglikelihooddata1 == max(loglikelihooddata1), arr.ind = TRUE)
maxindex
max(loglikelihooddata1)
# 最大尤度のインデックス（行、列）を取得
muindex <- maxindex[1, 1]
# === 最尤推定値を特定して表示 ===
# インデックスから最尤推定値を特定
mlemu <- muvalue[muindex]
mlemu
plot(muvalue, loglikelihooddata1, type="l",lwd=2,
xlab="平均", ylab="対数尤度",
main="対数尤度と平均（標準偏差＝10と既知",cex.main=0.5)
abline(v=mlemu, col="red", lty=1,lwd=2)
text(mlemu, max(loglikelihooddata1),
labels=paste("MLE Mean =", round(mlemu,2)), pos=1,cex=0.5)
平均と標準偏差の同時推定
# 母数平均と標準偏差の同時推定
# === サンプルデータ ===
sampledata <- c(37,67)
# 尤度関数は上記と同じ、ただし、sdは可変となりsd = sdvalue[j]に変更
# 平均値の候補を0.1から80まで0.1刻みで生成
muvalue <- seq(from =0.1, to = 80, by = 0.1)
length(muvalue)
# 標準偏差の候補を0.1から16まで0.1刻みで生成
sdvalue <- seq(from = 0.1, to =16, by = 0.1)
length(sdvalue)
# グリッドのすべての組み合わせに対して対数尤度を計算
calclikelihood <- function(mu, sd, sample) {
densities <- dnorm(x = sample, mean = mu, sd = sd, log = TRUE)
likelihood <- sum(densities)
return(likelihood)
}
# for loopを使って尤度データを計算
loglikelihooddata2 <- matrix(0,nrow=length(muvalue),ncol=length(sdvalue))
for (i in 1:length(muvalue)) {
for (j in 1:length(sdvalue)) {
loglikelihooddata2[i, j] <- calclikelihood(mu = muvalue[i],
sd = sdvalue[j],
sample = sampledata)
}}
length(loglikelihooddata2)
# === 最大尤度を与えるインデックスを特定 ===
# which()関数で最大値の位置を特定
maxindex <- which(loglikelihooddata2 == max(loglikelihooddata2), arr.ind = TRUE)
maxindex
max(loglikelihooddata2)
# 最大尤度のインデックス（行、列）を取得
muindex <- maxindex[1, 1]
sdindex <- maxindex[1, 2]
# === 最尤推定値を特定して表示 ===
# インデックスから最尤推定値を特定
mlemu <- muvalue[muindex]
mlesd <- sdvalue[sdindex]
# persp3d関数で3Dプロットを作成
library(rgl)
persp3d(x = muvalue, y = sdvalue, z = loglikelihooddata2,
aspect=c(1,1,1),col="lightblue")
# 最尤推定値の点を追加
points3d(x = mlemu, y = mlesd, z = max(loglikelihooddata2), col = "red", size = 10)
text3d(
x = mlemu,
y = mlesd,
z = max(loglikelihooddata2),
texts = paste("MLE mu=", mlemu,"MLE sd=", mlesd),
pos = 1,
adj = c(5, 5), cex=0.8,
color = "black"
)
mlemu
mean(sampledata)
mlesd
参考文献
1. 森田優三 1968 統計数理入門 日本評論社
2.稲垣宣生 1992 数理統計学 裳華房
3.Damodar N Gujarati 2003 Basic econometrics McGraw-Hill Companies Inc.
4.A.M.Mood,F.A.Graybill 1963 Introduction to the theory of statistics McGraw-Hill Inc.