0

Using any type of data, in my case data from three gamma distributions mixed, the goal is to parameterise the distributions' theta and the distribution weight alpha that are all 0< and sum to 1.

How do I restrict the alphas to sum to one, where no alpha is negative?

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
TylerH
  • 20,799
  • 66
  • 75
  • 101
BDL
  • 1
  • Here's a relevant transformation: http://stackoverflow.com/questions/23835294/converting-optim-to-constroptim-in-r. – MrFlick Feb 23 '15 at 16:19
  • The mixture weights are computed from the "responsibility" p(component | datum), so by construction they are nonnegative and sum to 1. If you are seeing something else, there is a bug in your code; in that case maybe you should post the output you are seeing. – Robert Dodier Feb 23 '15 at 17:35
  • Hi Robert, Thank you for your input. The output I get from printing the alpha parameter is three vectors, one for each distribution, about the size: [,1] 1.866397e+170 [,2] 1.866397e+170 [,3] 4.046550e-02 – BDL Feb 24 '15 at 09:04
  • As you see they are not bound to any restrictions, how do I impose that? – BDL Feb 24 '15 at 09:05

0 Answers0