AST_402

Profile Likelihood [Page: 17]

x <- c(1,5,10,20,30,40,50,75,100,150,200)
n <- c(31,29,27,25,23,21,19,17,15,15,15)#of trials
y <- c(0,3,6,7,9,13,17,12,11,14,13)# of success
m1 <- glm(y/n~log(x),weights=n,family=binomial)
m1

## 
## Call:  glm(formula = y/n ~ log(x), family = binomial, weights = n)
## 
## Coefficients:
## (Intercept)       log(x)  
##      -4.453        1.296  
## 
## Degrees of Freedom: 10 Total (i.e. Null);  9 Residual
## Null Deviance:       119.3 
## Residual Deviance: 10.63     AIC: 44.43

#the profile likelihood CIs from the built-in R function confint()
confint(m1)

## Waiting for profiling to be done...

##                 2.5 %    97.5 %
## (Intercept) -5.783372 -3.317595
## log(x)       0.978029  1.667461

# Compute profile log-likelihood function manually for beta_2
# logLik is a built-in R function to compute log-likelihood of model

k <- 200
b2 <- seq(.7,2,length=k)
w <- rep(0,k)

for(i in 1:k){
   mm <- glm(y/n~1,offset=b2[i]*log(x),weights=n,family=binomial)
   #logit model: beta1+beta2 * log(xi)
   #y/n~1 indicates constant term beta1, offset is the other part.
   w[i] <- logLik(mm)
}

plot(b2,w,type="l",ylab="Profile log-likelihood",
     cex.lab=1.3,cex.axis=1.3)
abline(h=logLik(m1)-qchisq(.95,1)/2,lty=2)

#qchisq gives the quantile function

f <- function(b2,n,x,y,maxloglik){
   mm <- glm(y/n~1,offset=b2*x,weights=n,family=binomial)
   logLik(mm) - maxloglik + qchisq(.95,1)/2
}

uniroot(f,c(.8,1.1),n=n,x=log(x),y=y,maxloglik=logLik(m1))

## $root
## [1] 0.9779875
## 
## $f.root
## 'log Lik.' 1.084503e-05 (df=1)
## 
## $iter
## [1] 5
## 
## $init.it
## [1] NA
## 
## $estim.prec
## [1] 8.564891e-05

uniroot(f,c(1.5,1.8),n=n,x=log(x),y=y,maxloglik=logLik(m1))

## $root
## [1] 1.667439
## 
## $f.root
## 'log Lik.' -0.0001146866 (df=1)
## 
## $iter
## [1] 4
## 
## $init.it
## [1] NA
## 
## $estim.prec
## [1] 6.103516e-05

Newton Raphson [Page: 21]

library(readr)
dat <- read_csv("G:/My Drive/Turin/4th Year/AST 402/R_code/2. Newton_Raphson/lifetimes.csv")
head(dat)

str(dat)

## spc_tbl_ [48 × 1] (S3: spec_tbl_df/tbl_df/tbl/data.frame)
##  $ 1051: num [1:48] 1337 1389 1921 1942 2322 ...
##  - attr(*, "spec")=
##   .. cols(
##   ..   `1051` = col_double()
##   .. )
##  - attr(*, "problems")=<externalptr>

## Function for score, U
U <- function(y, theta){
  N <- length(y)
  u <- (-(2*N)/theta) + ( (2*sum(y^2)) / theta^3)
  return(u)
}

## Function for  U'
Uprime<- function(y, theta){
  N <- length(y)
  uprime <- ((2*N)/theta^2) - ( (2*3*sum(y^2)) / theta^4)
  return(uprime)
}

#### Function for Newton-Raphson Iterative Method
NR <- function(y, theta.old){
  ### Newton-Raphson Iterative Equation
  theta.new <- theta.old - U(y, theta.old)/Uprime(y, theta.old)
  
  i <- 2
  tol <- 0.0001
  
  while (abs(theta.new-theta.old)>tol){
    
    cat( "iteration",i, ":", "theta=", theta.new,"\n")
    theta.old <- theta.new 
    theta.new <- theta.old - U(y, theta.old)/Uprime(y, theta.old)
    
    i <- i+1
  } ## end of while loop
  
  return(theta.new)
  
} ## end of function NR

### Initial value of theta is the mean of the data, ybar
theta1 <- mean(dat$"1051")
NR(dat$"1051", theta1)

## iteration 2 : theta= 9763.304 
## iteration 3 : theta= 9980.692 
## iteration 4 : theta= 9993.497 
## iteration 5 : theta= 9993.538

## [1] 9993.538

Expectation-Maximization Algorithm

Binormal Mixture Model [Page: 24]

# Function for generating values from binormal mixture distribution
rnormix <- function(n, p, mu, sigma){
  K <- length(p)
  b <- sample.int(K, n, replace=TRUE, prob=p)#sample.int(K, n) is 
           #an integer vector of length n with elements from 1:K.
  
  simulations <- rnorm(n, mu[b], sigma[b])
  return(simulations)
}

probs <- c(0.4, 0.6) #probabilities
mus <- c(52, 80) #mean values
sigmas <- c(5,5) #sigma values
sims <- rnormix(1000, probs, mus, sigmas)
plot(density(sims))

Newton Cotes

Trapezoid [Page: 53]

Use the trapezoidal rule to approximate the integral, \[ \int_0^{\frac{\pi}{4}} x \sin x \, dx \]

#Function for trapezoid
trapezoid <- function(f, a, b){
  h <- (b-a)
  approximate <- (h/2)*(f(a)+f(b))
  return(approximate)
}

func <- function(x){
  return(x*sin(x))
}

trapezoid(func, 0, pi/4)

## [1] 0.2180895

Composite Trapezoidal Rule [Page: 54]

For n=8 subintervals approximate the following integral, \[ \int_0^2 e^{2x} \sin(3x) \, dx \]

area_f <- function(a, b, n, f) {
h <- (b - a) / n
sum_f <- 0

for (j in 1:(n - 1)) {
    x[j] <- a + j * h    # Internal point
    sum_f <- sum_f + f(x[j])
}

area <- h / 2 * (f(a) + f(b) + 2 * sum_f)
return(area)
}

f <- function(x) {
  return(exp(2 * x) * sin(3 * x))
}

result <- area_f(0, 2, 8, f)
print(result)

## [1] -13.57598

Smart one

CTR<-function(a,b,n,f){
  h<-(b-a)/n
  x<-seq(a,b,length.out=n+1)
  (f(x[1])+2*sum(f(x[2:n]))+f(x[n+1]))*h/2
}

f1<-function(x){
  exp(2 * x) * sin(3 * x)
}

CTR(a = 0,b =2,n=8,f=f1)

## [1] -13.57598

Simpson’s Rule

\[ \int_0^{\frac{\pi}{4}} x \sin x \, dx \]

simpsons.rule <- function(f, a, b) {
  h <- (b - a) / 2
  x0 <- a
  x1 <- a + h
  x2 <- b
  s <- (h / 3) * (f(x0) + 4 * f(x1) + f(x2))
  return(s)
}

f <- function(x) {
  return(x * sin(x))
}

simpsons.rule(f, 0, pi/4)

## [1] 0.1513826

Composite Simpson’s Rule

\[ \int_0^2 e^{2x} \sin(3x) \, dx \]

composite.simpson <- function(f, a, b, n) {
  
h <- (b - a) / n
  
xj <- seq.int(from = a, to = b, by = h)
xj <- xj[-1]
xj <- xj[-length(xj)]
#sequence theke 1st and last data deduct
  
approx <- (h / 3) * (f(a) + 2 * sum(f(xj[seq.int(from = 2, to = length(xj), by = 2)])) + 4 * sum(f(xj[seq.int(from = 1, to = length(xj), by = 2)])) + f(b))
  
return(approx)
}

f <- function(x) {
  return(exp(2 * x) * sin(3 * x))
}

composite.simpson(f, 0, 2, 8)

## [1] -14.18334

Basic Monte Carlo Integration

Standard normal [Page: 55]

\[ I[-2, 1] = \int_{-2}^1 \theta^{2} \phi(\theta) \, dx \] \(\phi\) indicates standard normal pdf.

set.seed(1)
sample<-rnorm(100000)

mean(sample[sample > -2 & sample < 1] ^ 2)

## [1] 0.5755548

Cauchy distribution [Page: 56]

\[ I[e, \pi] = \int_e^\pi arctan(\theta^{1/3}) C(\theta|\mu = 3, \sigma = 2) \, d\theta \]

set.seed(99)
c_sample<-rcauchy(100000, 3, 2)

mean(atan(c_sample[c_sample > exp(1) & c_sample < pi] ^ (1/3)))

## [1] 0.96073

Importance Sampling

It is well known that the normal distribution is the limiting distribution of the binomial. That is, as the number of binomial experiments becomes large, the histogram of the binomial outcomes increasingly resembles a normal distribution. Using R, generate 25 values from a Binomial distribution:

\[ \theta_i \sim \text{Binomial}(10, 0.5) \] Use importance sampling with a normal approximation density as the proposal distribution to estimate the expected value \(\mathbb{E}(\theta)\).

even though the true expected value is known to be: \(\mathbb{E}(\theta) = np = 10 \times 0.5 = 5\).

theta <- round(rnorm(25, mean=5, sd=sqrt(2.5)), 0)

weight <- dbinom(theta, size=10, prob=0.5)/dnorm(theta, mean=5, sd=sqrt(2.5))

(imp.average <- sum(theta*weight)/sum(weight))

## [1] 5.02628

Rejection Sampling

Proposal Density: Beta

#generate thitaTilda from h
thitaTilda <- rbeta(2500, 2, 6)

#generate U from U(0,1)
U <- runif(2500)

M <- 1.67

#calculate the quantity require to accept that thitaTilda
quantity <- dbeta(thitaTilda, 2.7, 6.3)/(M*dbeta(thitaTilda, 2, 6))

#accept those thitaTilda for which (U<quantity)
theta <- thitaTilda*(U<quantity)
accepted.values <- theta[theta!=0]

#percent accepted
length(accepted.values)/2500

## [1] 0.6144

Proposal Density: uniform

thitaTilda <- runif(2500)
U <- runif(2500)
M <- 2.669

quantity <- dbeta(thitaTilda, 2.7, 6.3)/(M*dunif(thitaTilda))
#dunif(x)=1 actually

theta <- thitaTilda*(U<quantity)
accepted.values <- theta[theta!=0]

length(accepted.values)/2500

## [1] 0.3792

MCMC Algorithm

Gibbs Sampling

# priors
mu0<-1.9  ; t20<-0.95^2
s20<-.01 ; nu0<-1

# data
y<-c(1.64,1.70,1.72,1.74,1.82,1.82,1.82,1.90,2.08)
n<-length(y) ; mean.y<-mean(y) ; var.y<-var(y)

## starting values
set.seed(1)
S<-1000
PHI<-matrix(nrow=S,ncol=2)
PHI[1,]<-phi<-c( mean.y, 1/var.y)

for(s in 2:S) {
  # generate a new theta value from its full conditional
  mun<-  ( mu0/t20 + n*mean.y*phi[2] ) / ( 1/t20 + n*phi[2] )
  t2n<- 1/( 1/t20 + n*phi[2] )
  phi[1]<-rnorm(1, mun, sqrt(t2n) )
  
  # generate a new sigma^2 value from its full conditional
  nun<- nu0+n
  s2n<- (nu0*s20 + (n-1)*var.y + n*(mean.y-phi[1])^2 ) /nun
  phi[2]<- rgamma(1, nun/2, nun*s2n/2)
  
  PHI[s,]<-phi         }

#### Posterior quantiles
quantile(PHI[,1],c(.025,.5,.975))

##     2.5%      50%    97.5% 
## 1.707282 1.804348 1.901129

quantile(PHI[,2],c(.025,.5, .975))

##      2.5%       50%     97.5% 
##  17.48020  53.62511 129.20020

quantile(1/sqrt(PHI[,2]),c(.025,.5, .975))

##       2.5%        50%      97.5% 
## 0.08797701 0.13655763 0.23918408

Metropolis hastings [Page: 62]

nsim <- 10000
theta <- vector()
theta[1] <- 1 #initialize theta value

for(i in 2:nsim){
  theta_tilda <- rnorm(1, mean=theta[i-1], sd=1)
  alpha <- dgamma(theta_tilda, 5, 5)/dgamma(theta[i-1], 5, 5)
  
  if(runif(1)<alpha){
    theta[i] <- theta_tilda
  } else{
    theta[i] <- theta[i-1]
  }
}#end of for loop

hist(theta)