使用任何类型的数据,在我的例子中,来自三个混合伽马分布的数据,目标是参数化分布的 theta 和分布权重 alpha,它们都是 0< 并且总和为 1。
如何将 alpha 的总和限制为一个,其中没有 alpha 为负数?
func1<-function(x, comp_numb=3, a=c(0,0,0), k=c(0,0,0), th=c(0,0,0), maxit=10, tol=0.3){
a1<-c(1/3,1/3,1/3)
k1<-c(0.1,0.1,0.1)
th1<-c(0.1,0.1,0.1)
m<-2
sumcomp<-numeric(comp_numb)
out<-matrix(,length(x),comp_numb)
the<-matrix(,maxit, comp_numb)
alp<-matrix(,maxit, comp_numb)
if(a==0 && k==0 && th==0){ # controlling initial parameter values
a<-a1
k<-k1
th<-th1
}
while(m<=maxit+1){ # iteration
for (l in 1:length(x)){
for (i in 1:comp_numb){ # estimating the denominator sum in fY|X=x
sumcomp[i]<-(a[i]*x[l]^(k[i]-1)*exp(-x[l]/th[i]))/(gamma(k[i])*(th[i])^k[i])
# estimating the numerator in fY|X=x
out[l,i]<-(a[i]*((x[l]^(k[i]-1)*exp(-x[i]/th[i]))/(gamma(k[i])*(th[i])^k[i]))/sum(sumcomp))
} #close out loop
} #close x loop
for (i in 1:comp_numb){ # estimating parameters
a[i]<- sum(out[,i])/length(x) #replacing alphas
th[i]<- sum(x*out[,i])/(k[i]*sum(out[,i])) # replacing thetas
alp[m-1,i]<-a[i] # stacking parameters
the[m-1,i]<-th[i]
if(m>3){
difalp<-abs(alp[m-1,]-alp[m-2,]) ##change in alpha estimate
difthe<-abs(the[m-1,]-the[m-2,]) ##change in theta estimate
}
} # close parameter estimator
if(m>3){
if(all(difalp<tol) && all(difthe<tol))break} ##breaks while loop if dif<tol
m<-m+1
} # close while loop parameter maximator
print(alp)
print(the)
} # close function