- comp <- sample(c(0, 1), size = 100, prob = c(0.7, 0.3), replace = T)
- rnorm(100, mean = ifelse(comp == 0, 0, 1), sd = ifelse(comp == 0, 1, 2))
猜测分布是两个正态分布的混合，需要估计出函数中的5个参数：p、μ1、μ2、σ1、σ2。
在R中编写对数似然函数时，5个参数都存放在向量para中，由于nlminb()是计算极小值的，因此函数function中最后返回的是对数似然函数的相反数。
|
1
2
3
4
5
6
7
8
|
> l1=
function
(para)
+ {
+ f1=
dnorm
(waiting,para[2],para[3])
+ f2=
dnorm
(waiting,para[4],para[5])
+ f=para[1]*f1+(1-para[1])*f2
+ l1=
sum
(
log
(f))
+
return
(-11)
+ }
|
做参数估计，使用nlminb()之前最大的要点是确定初始值，初始值越接近真实值，计算的结果才能越精确。我们猜想数据的分布是两个正态的混合，概率P直接用0.5做初值即可。通过直方图中两个峰对应的x轴数值(大概为50和80>，就可以将初值设定为μ1和μ2。而概率P处于((0,1)区间内，参数σ1，σ2是正态分布的标准差，必须大于0，所以通过lower和upper两个参数进行一定的约束。
|
1
2
3
4
5
6
7
8
9
10
11
12
|
> geyser.est=
nlminb
(
c
(0.5,50,10,80,10),l1,lower=
c
(0.0001,-
Inf
,0.0001,-
Inf
,0.0001),upper=
c
(0.9999,
Inf
,
Inf
,
Inf
,
Inf
))
>
options
(digits=3)
> geyser.est$par
[1] 0.308 54.203 4.952 80.360 7.508
> p=geyser.est$par[1]
> mu1=geyser.est$par[2];sigma1=geyser.est$par[3]
> mu2=geyser.est$par[4];sigma2=geyser.est$par[5]
> x=
seq
(40,120)
>
#将估计的参凌丈函数代入原密度函数
> f=p*
dnorm
(x,mu1,sigma1)+(1-p)*
dnorm
(x,mu2,sigma2)
>
hist
(waiting,freq=F)
>
lines
(x,f)
|