#source("c:/Work07/ShortCourse07/Programs/Students-t-moment.txt",print.eval=TRUE)#
#This does the t moment integral
nsim<-100000
den<-c(1:nsim)
cut<-2.1
h<-function(x)(x^5)*dt(x,df=10)
integrate(h, cut, Inf)
#---------------------Students t candidate------------------
t<-rt(nsim,df=10)
mean((t^5)*(t>cut))
tp<-cumsum((t^5)*(t>cut))/den
#---------------------Normal candidate------------------
z<-rnorm(nsim)
mean((z^5)*(dt(z,df=10)/dnorm(z))*(z>cut))
zp<-cumsum((z^5)*(dt(z,df=10)/dnorm(z))*(z>cut))/den
#---------------------Cauchy candidate------------------
c<-rt(nsim,df=1)
mean((c^5)*(dt(c,df=10)/dt(c,df=1))*(c>cut))
cp<-cumsum((c^5)*(dt(c,df=10)/dt(c,df=1))*(c>cut))/den
#---------------------Uniform candidate------------------
u<-1/runif(nsim,min=0,max=1)
mean((u^7)*dt(u,df=10)*(u>cut))
up<-cumsum((u^7)*dt(u,df=10)*(u>cut))/den
#**********************Plot Pre Processing***************************************
m<-1000		#number of points to plot
step<-round(nsim/m)
tplot<-tp[1];zplot<-zp[1];cplot<-zp[1];uplot<-up[1];xplot<-1
for(i in 1:m) {tplot<-c(tplot,tp[i*step]);zplot<-c(zplot,zp[i*step])}
for(i in 1:m) {cplot<-c(cplot,cp[i*step]);uplot<-c(uplot,up[i*step]);xplot<-c(xplot,i*step)}
#*********************************Plotting*******************************************
#-----------------Plot-----------------------------------
ymax<-15
plot(xplot,tplot,type="l",ylim=c(0,ymax),lty=1,ylab="mean")
par(new=T)
plot(xplot,zplot,type="l",ylim=c(0,ymax),lty=3,ylab="")
par(new=T)
plot(xplot,cplot,type="l",ylim=c(0,ymax),lty=2,ylab="")
par(new=T)
plot(xplot,uplot,type="l",ylim=c(0,ymax),lty=5,ylab="")