0

使用任何类型的数据，在我的例子中，来自三个混合伽马分布的数据，目标是参数化分布的 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
4

0 回答 0