#Galaxy Data Question

pdf('GalaxyDataFits.pdf',paper='USr',width=11,height=9)
K<-5
source('Ex4-2015-EMsub.R')
(BIC.5<--2*Ly+(3*K-1)*log(n))

K<-4
source('Ex4-2015-EMsub.R')
(BIC.4<--2*Ly+(3*K-1)*log(n))

K<-3
source('Ex4-2015-EMsub.R')
(BIC.3<--2*Ly+(3*K-1)*log(n))
dev.off()

###############################################
Task II
###############################################

#EM for the Student t(theta,sigma,nu)
#Location only

set.seed(23842)
n<-10
nu<-5
theta<-5
y<-rt(n,nu)+theta

th<-0
del<-Inf
tol<-1e-10
while(del > tol){
	
	w<-1/(nu+(y-th)^2)
	new.th<-sum(w*y)/sum(w)
	del<-abs(th-new.th)
	th<-new.th
	#print(w)
	print(th)
}

#Check
th.vec<-seq(-10,10,by=0.01)
dmat<--log(nu+t(outer(y,th.vec,"-")^2))
lvec<-apply(dmat,1,sum)
plot(th.vec,lvec,type='l')
abline(v=th,col='red')


#EM for the Student t(theta,sigma,nu)
#Location and scale

set.seed(23842)
n<-10
nu<-5
theta<-5
y<-rt(n,nu)+theta

th<-0
sig<-1
del<-Inf
tol<-1e-10
while(del > tol){
		
	z<-(y-th)/sig
	w<-1/(nu+z^2)
	new.th<-sum(w*y)/sum(w)
	new.sig<-sqrt((nu+1)*sum((y-new.th)^2/(nu+z^2))/n)
	del<-sqrt((th-new.th)^2+(sig-new.sig)^2)
	th<-new.th
	sig<-new.sig
	#print(w)
	print(c(th,sig))
}

#Check
tfunc<-function(pars,yv,nuv){
	th<-pars[1]
	sig<-exp(pars[2])
	like<-sum(dt((yv-th)/sig,nuv,log=T)-log(sig))
	return(-like)
}

tfunc(c(0,1),y,nu)

omax<-optim(c(0,1),fn=tfunc,yv=y,nuv=nu,method='BFGS')
th.hat<-omax$par[1]
sig.hat<-exp(omax$par[2])

print(c(th.hat,sig.hat)) #via optim
print(c(th,sig)) #via EM

##############################################################

#For the galaxies data

set.seed(23842)
nu.vec<-4:30

y<-MASS::galaxies
n<-length(y)

for(i in 4:30){

	nu<-i

	th<-0
	sig<-1
	del<-Inf
	tol<-1e-10
	istep<-0
	while(del > tol){
		z<-(y-th)/sig
		w<-1/(nu+z^2)
		new.th<-sum(w*y)/sum(w)
		new.sig<-sqrt((nu+1)*sum((y-new.th)^2/(nu+z^2))/n)
		del<-sqrt((th-new.th)^2+(sig-new.sig)^2)
		th<-new.th
		sig<-new.sig
		istep<-istep+1
	}
	print(c(i,th,sig,istep))

}