Test_L.Rmd

---
title: "Laplace's test"
author: "Huimin ZHANG"
date: "2022/2023"
output: pdf_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

# Theory of Laplace's test

First, let's define the test statistic. 

Test statistic : $L = \sum_{i=1}^n \frac{T_i}{T^*}$ where $T*$ is a fixed units of time

The interarrival times $W_i$ have an exponential distribution with parameter $\lambda$ and $\mathbb{E}[\mathcal{E(\lambda)}] = \frac{1}{\lambda}$. Thus, if $\lambda$ is constant,  the interarrival times $W_i$ are constant in average, and if $\lambda$ is increasing, the interarrival times are smaller and the $T_i$ are bigger. Therefore, if $\lambda$ is increasing, $\frac{T_i}{T^*}$ is greater than if $\lambda$ is constant. 
So, we reject $\mathcal{H}_0$ when $L$ have big values, i.e. the rejection zone is $\mathcal{R}_{\alpha} = \left\{ L \ge l \right\}$



# Simulation 

```{r}
lambda_exp <- function(alpha, beta, Tstar) {
  lambda <- function(t){
    return(alpha*exp(beta*t))
  }
  LAMBDA <- function(t){
    return((alpha/beta)*(exp(beta*t) - 1))
  }
  
  quantile_function <- function(U){
    return((1/beta) * log(U*(exp(beta*Tstar) - 1) + 1))
  }
  
  return(c(lambda, LAMBDA, quantile_function))
}
```


```{r}
plot_PP <-function(PP){
  # plot the counting process (with jumps = 1): 
  plot(c(0,PP),seq(0,length(PP)),type="s",xlab="time t",ylab="number of events by time t")
  # add the arrival times: 
  points(PP, rep(0,length(PP)),type="p")
  # link the arrival times with the counts:
  lines(PP, seq(1,length(PP)),type="h",lty=2)
}
```



```{r}
expo <- lambda_exp(alpha = 2, beta = 0, Tstar = 5)
lambda <- expo[[1]]
LAMBDA <- expo[[2]]
quantile_function <- expo[[3]]
```

```{r}
L <- function(Ti, T_star){
  return(sum(Ti/T_star))
}
```


## Homogenous Poisson Process

```{r}
# T fixed
simulPPh <- function(alpha, beta, T_star)
{
  lambda <- lambda_fct(alpha, beta, 1)
  Y <- rpois(1, lambda*T_star)
  U <- runif(Y, min = 0, max = T_star)
  return(sort(U))
}
```


```{r}
# simulate a homogeneous Poisson process with T fixed
T_star <- 10
alpha <- 2
beta <- 0
lambda <- lambda_fct(alpha, beta, 1)
PPh = simulPPh(lambda, T_star)

# plot the counting process (with jumps = 1): 
plot(c(0,PPh),seq(0,length(PPh)),type="s",xlab="time t",ylab="number of events by time t")
# add the arrival times: 
points(PPh, rep(0,length(PPh)),type="p")
# link the arrival times with the counts:
lines(PPh, seq(1,length(PPh)),type="h",lty=2)
```

```{r}
test1 <- function(PPh, T_star)
{
  n <- length(PPh)
  L_obs <- sum(PPh/T_star)
  Z_obs <- (L_obs - n/2)/sqrt(n/12)
  return(1 - pnorm(Z_obs))
}

print(test1(PPh, T_star))
```

```{r}
get_size <- function(K, alpha, get_pval, simulHPP, Tstar, alpha_lamb, beta){
  nb.rejects <- 0
  lambda <- lambda_fct(alpha_lamb, beta, 1)
  for(k in 1:K){
    N <- 0
    while(N <= 1){
      PPi <- simulHPP(lambda, Tstar)
      N <- length(PPi)
    }
    pval = get_pval(PPi, Tstar)
    nb.rejects = nb.rejects + (pval<=alpha)  
  }
  return(nb.rejects/K)
}
```

```{r}
alpha <- 0.05
Tstar <- 5
alpha_lamb <- 2
beta <- 0
test_size <- get_size(K = 1000, alpha, test1, simulPPh, Tstar, alpha_lamb, beta)
cat("\nSize of the test : ", test_size)
```



## Inhomogeneous Poisson process

```{r}
lambda_int <- function(alpha, beta, t){
  return((alpha/beta)*(exp(beta*t) - 1))
}

quantile_function <- function(U, Tstar, beta){
  return((1/beta) * log(U*(exp(beta*Tstar) - 1) + 1))
}

simulPPi <- function(Tstar, alpha, beta){
  # simulate the number of event under a Poisson distribution
  N <- rpois(1, lambda_int(alpha, beta, Tstar))
  # simulate a uniform sample
  U <- runif(N,0,1)
  # apply the quantile function (inverse of the cumulative distribution function)
  S <- quantile_function(U, Tstar, beta)
  # sort the sample
  return(sort(S)[1:N])
}
```


```{r}
# simulate an inhomogeneous Poisson process with T fixed
T_star <- 3
alpha <- 1
beta <- 1
PPi = simulPPi(T_star, alpha, beta)

# plot the counting process (with jumps = 1): 
plot(c(0,PPi),seq(0,length(PPi)),type="s",xlab="time t",ylab="number of events by time t")
# add the arrival times: 
points(PPi, rep(0,length(PPi)),type="p")
# link the arrival times with the counts:
lines(PPi, seq(1,length(PPi)),type="h",lty=2)
```

```{r}
print(test1(PPi, T_star))
```

```{r}
get_size <- function(K, alpha, get_pval, simulPPi, Tstar, alpha_lamb, beta){
  nb.rejects <- 0
  for(k in 1:K){
    N <- 0
    while(N <= 1){
      PPi <- simulPPi(T_star, alpha_lamb, beta)
      N <- length(PPi)
      
    }
    pval = get_pval(PPi, Tstar)
    nb.rejects = nb.rejects + (pval<=alpha)  
  }
  return(nb.rejects/K)
}
```



```{r}
alpha <- 0.05
Tstar <- 3
alpha_lamb <- 2
beta <- 1
test_size <- get_size(K = 1000, alpha, test1, simulPPi, Tstar, alpha_lamb, beta)
cat("\nSize of the test : ", test_size)
```