In hudie-lab/SC19036: The R Package For StatComp

Question 1

question 1

Answer 1

as for the question, I choose the acceptance-rejection method to solve it. so it is necessary to choose the $Y\sim g(.)$ as the threadshold.

the gamma $X\sim g(x;\beta,\alpha)=\frac{\beta^{\alpha}}{\Gamma(\alpha)}x^{\alpha-1} e^{-\beta x}$ distribution has ths similar distribution. here, I set $\alpha=2$,$\beta=1$, so the $g(x)=x e^{-x}$

with the $ f(x)=\frac{x}{\sigma^2} e^{-\frac{x^2}{2\sigma^2} }$, it is easy to calculate that $f(x)/g(x)<=\sigma^2$ which stands $c$.

##sample generator
raylsamp=function(N,sigma){
    s2=sigma^2
    k=1
  x=0
  while(k<=N){
    ##N : the number of samples
  u=runif(1)
  ##y : the random variable from g(x)
  y=rgamma(1,shape=2,scale=1)
  ##accept-reject step
  if(u<=exp(-y^2/(2*s2)+y-s2/2)) {
    k=k+1
    x[k]=y
  }
  }
  ##print the samples
  x[-1]
  }

##the samples from the acceptance-rejection method above
N=50000
##set the parameter sigma
sigma=2

AR_x=raylsamp(N,sigma)
hist(AR_x,freq = FALSE,main="the acceptance-rejection samples from rayleigh distribution")
y=seq(0,max(AR_x),0.1)
lines(cbind(y,y/sigma^2*exp(-y^2/2/sigma^2)),col="red")
abline(v=sigma,col="blue")
text(2.25,0.25,"theoretical sigma",col="blue")
text(7,0.2,"red lines:theoretical f(x)",col="red")

from the hist diagram, we can find that the mode of the samples is close to the theoretical $\sigma$, and the hists are close to the theoretical pdf $f(x)$.

several sets of $\sigma$

##the samples from the acceptance-rejection method above
N=1000
##set the parameter sigma
sigma=2:5

#par(mfrow=c(2,2))
for(i in 1:length(sigma)){
  s=sigma[i]
  AR_x=raylsamp(N,s)
hist(AR_x,breaks=50,freq = FALSE,main=paste0("the acceptance-rejection samples from rayleigh distribution with sigma=",s))
y=seq(0,max(AR_x),0.1)
lines(cbind(y,y/s^2*exp(-y^2/2/s^2)),col="red")
abline(v=s,col="blue")
text(s,0.1,"theoretical sigma",col="blue")
text(max(AR_x)-5,0.2,"red lines:theoretical f(x)",col="red")
}

we can see from the results, that for each $\sigma$, the samples using the acceptance-rejection method fit the theoretical pdf well because the modes of the samples are all close to the theoretical $\sigma$.

Question 2

question 2

Answer 2

N=1000
####random numbers from the first normal distribution 
x1=rnorm(N)
####random numbers from the second normal distribution 
x2=rnorm(N,mean=3)

p=0.75
####the random variable which stands for the mixing probability
r=sample(c(0,1),N,replace=TRUE,prob=c(p,1-p))
####the mixture random samples
x=(1-r)*x1+r*x2

hist(x,breaks=50,freq = FALSE,ylim=c(0,0.5),main=paste0("mixture of normal distribution with p=","p"))
####generate the density imposed
y=seq(min(x1,x2),max(x1,x2),0.1)
lines(cbind(y,p*dnorm(y)+(1-p)*dnorm(y,mean=3)),col="red")
text(4,0.4,"red lines for the density",col="red")

from the hist diagram, it is evident that the hist of samples is close to the density superimposed. there seems to exist two peaks in the density plot.

different probility $p$

N=1000
x1=rnorm(N)
x2=rnorm(N,mean=3)

#par(mfrow=c(4,3))
####set the different probabilities of p
pi=seq(0.2,0.8,0.05)
for(i in 1:(length(pi)-1)){
  p=pi[i]
  r=sample(c(0,1),N,replace=TRUE,prob=c(p,1-p))
x=(1-r)*x1+r*x2

hist(x,breaks=50,freq = FALSE,ylim=c(0,0.5),main=paste0("mixture of normal distribution with p=",p))
y=seq(min(x1,x2),max(x1,x2),0.1)
lines(cbind(y,p*dnorm(y)+(1-p)*dnorm(y,mean=3)),col="red")
text(max(x)-2,0.4,"red lines for the density",col="red")
}

from the plots above, we can see that the bimodal mixture turns to be more evident with p incresing to 0.5, then the figure vanishes with p increasing from 0.5.

Question 3

question 3

Answer 3

####function for Wishart distribution
rwishart=function(n,sigma){
  ##d is determined by the dimension of sigma
  d=nrow(sigma)
  M=list()
    ##the Bartlett's decomposition
    T=matrix(rep(0,d*d),d,d) 
  for(i in 1:d){
    for(j in 1:i){
        if(i>j) T[i,j]=rnorm(1)
        if(i==j) T[i,j]=sqrt(rchisq(1,df=(n-i+1)))
    }
  }
    ##the choleski factorization
  L=t(chol(sigma))
  Mt=(L%*%T)%*%t(L%*%T)

  ##the random matrix on the wishart distribution
  Mt
}

the random from my own function

library(SC19036)
n=5
d=2
####the sigma needs to be symetric and positive limited.
sigma=matrix(c(5,1,1,3),d,d)
(M=rwishart(n,sigma))

compare the simulating function and the true random matrix

I find the density function of wishart distribution, it is the reference.

from the book, the density function of wishart can get. here, I calculate the densities of every random matrix using my own function $rwishart()$ compared with the densities of the simulations by the R function $rWishart()$. (the two names are different by the W and w.)

##the density function of wishart distribution

dwishart=function(A,n,digma){
  d=nrow(sigma)
  det(A)^((n-d-1)/2)*exp(-0.5*sum(diag(solve(sigma)%*%A)))/2^(d*n/2)/gamma(n/2)/det(sigma)^(d/2)
}

n=5
N=500

##the random generator in stats Rpackage
Ms=rWishart(N,df=nrow(sigma),sigma)
dw=matrix(0,nrow=N,ncol=2)
for(i in 1:N){
  dw[i,1]=dwishart(rwishart(n,sigma),n,sigma)
  dw[i,2]=dwishart(as.matrix(Ms[1:2,1:2,i]),n,sigma)
}

##the random generator in stats Rpackage
Ms=rWishart(n,df=nrow(sigma),sigma)

library(rgl)
ar=matrix(0,ncol=6,nrow=N)
for(i in 1:N){
  m=rwishart(n,sigma)
  mt=rWishart(1,nrow(sigma),sigma)
  ar[i,1:3]=unique(c(m))
  ar[i,4:6]=unique(c((mt[1:2,1:2,1])))
} 
#plot3d(cbind(ar[,1:3],ar[,4:6]),col=c("blue","red"))