I created a EM algorithm for Generalized hyperbolic distribution. I want to estimate mutheldaplus, sigmatheldaplus, betasigmaplus in my code. After getting use these value , then my iteration have to be begin of this code. But I can not to do iteration part.
Can you help me use my code and get iteration ? Do know any useful code for EM algorithm for Generalized Hyperbolic library(QRMlib) library(ghyp) ############ simulation part simulation<-function(n,lambda,mu,thelda,gamma,sigma,beta){ set.seed(235) chi<-thelda^2 psi<-gamma^2 W <- rGIG(n, lambda, chi, psi); Z <- rnorm(n,0,1); y<-mu + beta * W + sqrt(W) * Z *gamma; for (i in 1:n){ theldastar<-rep(0,n) zi<-rep(0,n) ti<-rep(0,n) muthelda<-mu gammathelda<-thelda*gamma sigmathelda<-(thelda^2)*sigma betathelda<-(thelda^2)*sigma*beta lambdastar<-lambda-0.5 theldastar[i]<-sqrt(1+((y[i]-muthelda)/sigmathelda)^2) gammastar<-sqrt((gammathelda^2)+((betathelda/sigmathelda)^2)) klambda1<-besselM3(lambdastar+1, x=2, logvalue=FALSE) klambda<-besselM3(lambdastar,x=2,logvalue=FALSE) klambda2<-besselM3(lambdastar-1,x=2,logvalue=FALSE) zi[i]<-((theldastar[i]*klambda1*(theldastar[i]*gammastar))/(gammastar*klambda*theldastar[i]*gammastar)) ti[i]<-((gammastar*klambda2*(theldastar[i]*gammastar))/(theldastar[i]*klambda*theldastar[i]*gammastar)) zimean<-sum(zi)/n timean<-sum(ti)/n mutheldaplus<-(zimean*(1/n)* sum((ti[i]*y[i])-mean(y)))/((zimean*timean)-1) betatheldaplus<- sum(y[i]- mutheldaplus)/(n*zimean) sigmatheldaplus<-((1/n)*sum((ti[i]*((y[i]-mutheldaplus)^2))-(2*betatheldaplus*(y[i]-mutheldaplus))-((betatheldaplus^2)*zi[i]))) print(muthelda) print(mutheldaplus) print(betathelda) print(betatheldaplus) print(sigmathelda) print(sigmatheldaplus) return(ti) } } a<-simulation(20000,-0.5,0,1,1,1,0) [[alternative HTML version deleted]] ______________________________________________ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.