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

David Scott d.scott at auckland.ac.nz
Tue Sep 21 13:31:07 CEST 2010

```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
> 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

```