[R] Need help for EM algorithm ASAP !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

[David Scott]

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