snes1982 at hotmail.com
2010-Sep-21 05:32 UTC
[R] Need help for EM algorithm ASAP !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
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]]
David Scott
2010-Sep-21 11:31 UTC
[R] Need help for EM algorithm ASAP !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
The urgency and the vague description of your problem strongly suggest that this is homework. This list is not for homework---see the posting guide at the bottom of every message. Nonetheless since I know this problem reasonably well I will offer some comments. QRMlib is a package created to accompany a book. If you read that book you would see that it fits the generalized hyperbolic to data using the EM algorithm. If you have QRMlib you have an implementation of the EM algorithm. Also why write code to simulate from the generalized hyperbolic (y in your simulation function below) when you have QRMlib and ghyp, both of which have functions for simulating from the generalized hyperbolic? Your code is pretty difficult to follow, with random indenting and zero comments. The structure of the iteration is totally confused as well. Not too many marks if you handed something like this in to me to grade. David Scott On 21/09/2010 5:32 p.m., snes1982 at hotmail.com wrote:> 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 at 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.-- _________________________________________________________________ David Scott Department of Statistics The University of Auckland, PB 92019 Auckland 1142, NEW ZEALAND Phone: +64 9 923 5055, or +64 9 373 7599 ext 85055 Email: d.scott at auckland.ac.nz, Fax: +64 9 373 7018 Director of Consulting, Department of Statistics