1 Introdução

Neste trabalho, implementamos os resultados, obtidos na dissertação de mestrado, para o Custo Mínimo Esperado (Teórico) entre eventos de dois processos de Poisson independentes e simulamos, a partir de amostras aleatórias de tempos de chegada, obtendo o Custo Mínimo Estimado.

A programação foi realizada na linguagem R (https://www.r-project.org/), utilizando o IDE RStudio (link).

IDE é o ambiente de desenvolvimento integrado (Em inglês, Integrated Development Environment).

2 Valores Teóricos

2.1 Distância Esperada (Em termos da Beta Incompleta)

Considere dois processos de Poisson independentes com taxas de chegada \(\lambda_1\) e \(\lambda_2\), cujos tempos de chegadas são, respectivamente, \(X_1,X_2,\cdots\) e \(Y_1,Y_2,\cdots\). Então, a fórmula

\[E[|X_k-Y_k|]=k\bigg(\frac{1}{\lambda_1}-\frac{1}{\lambda_2}\bigg)+2k\bigg(\frac{I_x(k,k+1)}{\lambda_2}-\frac{I_x(k+1,k)}{\lambda_1}\bigg),\] é válida para \(r\geq0\), \(k\geq1\) inteiros e \(\lambda_1,\lambda_2>0\).

\(I_x(a,b)\) é a função beta incompleta regularizada e \(x=\lambda_1/(\lambda_1+\lambda_2)\).

(Demonstração: Ver a dissertação de Mestrado.)

E<-function(k,l1,l2){
  
  x<-l1/(l1+l2) 
  
  res<-k*(1/l1-1/l2)+2*k*((1/l2)*pbeta(x,k,k+1)-(1/l1)*pbeta(x,k+1,k))
  
  return(res)
}

2.2 Custo Mínimo Esperado (Abordagem da Soma)

O Custo é dado pelo soma da Distância Esperada:

\[ C_{opt}(n,\lambda_1,\lambda_2)=\sum\limits_{k=1}^{n}E[|X_k-Y_k|]\]

Custo_teorico<-function(n,l1,l2){
       
                 custo<-numeric(length(n))
       
                      for (j in 1:length(n)){
             
                              custo[j]<-sum(E(1:n[j],l1,l2))
                      
                              }
                return(custo)
}

3 Valores Simulados

3.1 Custo Mínimo Estimado

Sejam as amostras aleatórias de tempos de chegada \(X_i=(X_{i1},X_{i2},\cdots, X_{im})\) e \(Y_i=(Y_{i1},Y_{i2},\cdots, Y_{im})\); \(i\in \{1,2,\ldots,n\}\) de dois processos de Poisson independentes com taxas \(\lambda_1\) e \(\lambda_2\), respectivamente. Temos que
\[X_{ij}\sim Gama(i,\lambda_1) \quad \mbox{e} \quad Y_{ij}\sim Gama(i,\lambda_2),\quad \forall\; (i,j)\in\{1,\ldots,n\}\times \{1,\ldots,m\}.\]

Nas matrizes abaixo, encontram-se as amostras de tamanho \(m\) dos \(n\) primeiros tempos de chegada. \[ \begin{bmatrix} X_{11} & X_{12} & \cdots & X_{1m}\\ X_{21} & X_{22} & \cdots & X_{2m}\\ \vdots & \vdots & \ddots & \vdots\\ X_{n1} & X_{n2} & \cdots & X_{nm} \end{bmatrix} \qquad \mbox{e}\qquad \begin{bmatrix} Y_{11} & Y_{12} & \cdots & Y_{1m}\\ Y_{21} & Y_{22} & \cdots & Y_{2m}\\ \vdots & \vdots & \ddots & \vdots\\ Y_{n1} & Y_{n2} & \cdots & Y_{nm} \end{bmatrix} \]

O Custo Mínimo Estimado é dado por

\[\hat{C}_{opt}(n,m,\lambda_1,\lambda_2)=\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{m}\frac{1}{m}\big|X_{ij}-Y_{ij}\big|,\] em que \(m\) é o tamanho das amostras e \(n\) é a quantidade de amostras.

Custo.est<-function(n,m,l1,l2) {
                    
                       # l1: taxa do processo 1
                       # l2: taxa do processo 2
                       # n: quantidade de amostras de tempo; (n primeiras chegadas)
                       # m: tamanho de cada amostra do tempo X_i 
                       
                       Custo<-numeric(length(n))

                                  for (j in 1:length(n)){
           
                                              mat.X<-matrix(0,nrow=n[j],ncol=m)
           
                                              mat.Y<-matrix(0,nrow=n[j],ncol=m)
                 
                                                      for (i in 1:n[j]){
                     
                                                               mat.X[i,]<-rgamma(m,i,l1) # m tempos X_i do 1o processo
  
                                                               mat.Y[i,]<-rgamma(m,i,l2) # m tempos Y_i do 2o processo
                                                              } 
                  
                                               Custo[j]<-sum(rowMeans(abs(mat.X-mat.Y)))
                                              }
 
  return(Custo)
}

3.1.1 Estimação dos parâmetros da distribuição Gama

Seja \(X_1,X_2,\ldots,X_n\) uma amostra aleatória de uma distribuição Gama(\(\alpha,\beta\)) com parâmetros de forma e de escala desconhecidos. Nesse caso, podemos usar o método de estimação por Máximo-Verossimilhança. Através desse método, mostra-se que

\[ \hat{\alpha}\approx\frac{0.5}{log(\bar{x})-\overline{log(x)}} \quad \mbox{e} \quad \hat{\beta}\approx\frac{\hat{\alpha}}{\bar{x}}. \]

3.1.2 Cálculo da Esperança do Custo Estimado

\[ \begin{aligned} E(\hat{C}_{opt})&=E\bigg[\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{m}\frac{1}{m}\big|X_{ij}-Y_{ij}\big|\bigg]\\ &\overset{(linear)}{=}\frac{1}{m}\sum\limits_{j=1}^{m}\sum\limits_{i=1}^{n}E\big[|X_{ij}-Y_{ij}|\big]\\ &\overset{(id)}{=}\frac{1}{m}\sum\limits_{j=1}^{m}\underbrace{\sum\limits_{i=1}^{n}E\big[|X_i-Y_i|\big]}_{C_{opt}\;\; \text{(independe de $j$)}}\\ &=\frac{1}{m}\sum\limits_{j=1}^{m}C_{opt}= C_{opt}. \end{aligned} \]

3.1.3 Cálculo da Variância do Custo Estimado

\[ \begin{aligned} Var(\hat{C}_{opt})&=Var\bigg(\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{m}\frac{1}{m}\big|X_{ij}-Y_{ij}\big|\bigg)\\ &=\frac{1}{m^2}\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{m}Var(X_{ij}-Y_{ij})\\ &=\frac{1}{m^2}\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{m}\bigg(Var(X_{ij})+Var(Y_{ij})\bigg)\\ &= \frac{1}{m}\sum\limits_{i=1}^{n}\bigg(\frac{i}{\lambda_{1}^2}+\frac{i}{\lambda_{2}^2}\bigg)=\frac{n(n+1)}{2m}\bigg(\frac{1}{\lambda_1^2}+\frac{1}{\lambda_2^2}\bigg). \end{aligned} \]

3.1.4 Convergência pela Lei dos Grandes Números

A Lei dos Grandes Números garante a seguinte convergência em probabilidade

\[\hat{C}_{opt} \to C_{opt}, \quad m\to \infty \]

3.1.5 Quantidade Pivotal baseada na estatística \(\hat{C}_{opt}\)

Defina-se a variável aleatória \(Z\) como segue

\[ Z:=\frac{\hat{C}_{opt}-E\big[\hat{C}_{opt}\big]}{\sqrt{Var(\hat{C}_{opt})}}, \]

em que \(E\big[\hat{C}_{opt}\big]=C_{opt}\), conforme mostrado em 3.1.2 .

Assintoticamente, \(Z\sim N(0,1)\). Logo, \(Z\) é uma quantidade pivotal (ver Casella and Berger (2002)).

Utilizaremos esta quantidade na construção de um Intervalo de Confiança para \(C_{opt}\).

3.1.6 Construção do Intervalo de Confiança para \(C_{opt}\)

A partir da Quantidade Pivotal, pode-se mostrar que o seguinte Intervalo de Confiança (IC) para \(C_{opt}\). \[IC_{100(1-\alpha)\%}(C_{opt})=\bigg[\;\hat{C}_{opt}-z_{(\alpha/2)}\sqrt{Var(\hat{C}_{opt})}\;,\;\hat{C}_{opt}+z_{(1-\alpha/2)}\sqrt{Var(\hat{C}_{opt})}\;\bigg]\] é válido.

## Variancia do Custo Estimado

Var_Custo.est<-function(n,m,l1,l2){(n*(n+1)/(2*m))*(1/l1^2+1/l2^2)}

## Intervalo de Confiança ao nível de 95%

Intervalo<-function(n,m,l1,l2){
                
                int<-Custo.est(n,m,l1,l2)+c(-1,1)*qnorm(0.975)*sqrt(Var_Custo.est(n,m,l1,l2))
               
           return( int )
}

set.seed(23/12/20)

n<-1200:1250

Int<-sapply(n,Intervalo,m=100,l1=0.95,l2=0.9)
LI<-Int[1,]
LS<-Int[2,]

CUSTO<-sapply(n,Custo_teorico,l1=0.95,l2=0.90)

plot(n,CUSTO,type="l")
lines(n,LI,col="red" )
lines(n,LS,col="blue" )
title(main="Custo Mínimo Esperado e Intervalo de Confiança a 95%")
legend("topleft", inset=.02, legend=c("Limite Superior", "Custo Esperado","Limite Inferior"),
       col=c("blue", "black","red"), lty=c(1,1,1), cex=0.8)

data <- data.frame(n = n,
                   custo = CUSTO,
                   li   =   LI,
                   ls = LS)

#install.packages("ggplot2")
library(ggplot2)
ggplot(data, aes(n, custo)) +
  geom_point() +
  geom_line() +
  geom_errorbar(aes(ymin = li, ymax = ls)) +
  labs(x = "n",
       y = "Custo",
       title = "Custo Mínimo Esperado de Transporte e Intervalo de Confianca ao nível 95%")

#ggsave(p, filename = "C:/Users/adolfoamds/Documents/test2.pdf",height = 10, width = 5)
par(mfrow=c(1,2))

l1<-0.95 ; l2<-0.90
n<-100:130
CUSTO<-sapply(n,Custo_teorico,l1=l1,l2=l2)
CUSTO.est<-sapply(n,Custo.est,m=100,l1=l1,l2=l2)


## Grafico do Custo 

plot(n,CUSTO.est,type="l",ylab="Custo")
lines(n,CUSTO,col="red")
title(main="Custo Esperado X Custo Estimado")
legend("topleft", inset=.02, expression(lambda[1]~"= 0.95, "~lambda[2]~" = 0.90"))
legend("bottomright",inset=.02,legend=c("Custo Esperado","Custo Estimado"),
       col=c("red", "black"), lty=c(1,1), cex=0.8)


plot(n,CUSTO.est,pch=1,ylab="Custo",main = "Custo Esperado X Custo Estimado")
lines(n,CUSTO, col="red")
legend("topleft", inset=.02, expression(lambda[1]~"= 0.95, "~lambda[2]~" = 0.90"))
legend("bottomright",inset=.02,legend=c("Custo Esperado","Custo Estimado"),
       col=c("red", "black"), lty=c(1,1), cex=0.8)

n0<-10:100
m0<-10:100

CUSTO<-sapply(n0,Custo_teorico,l1=l1,l2=l2)
CUSTO.est<-mapply(Custo.est,n0,m0,MoreArgs = list(l1=l1,l2=l2))
#x11()
#plot(m0,CUSTO.est,pch=1,ylab="Custo",main = "Custo Esperado X Custo Estimado")
#lines(n0,CUSTO, col="red")
#legend("topleft", inset=.02, expression(lambda[1]~"= 0.95, "~lambda[2]~" = 0.90"))
#legend("bottomright",inset=.02,legend=c("Custo Esperado","Custo Estimado"),
       #col=c("red", "black"), lty=c(1,1), cex=0.8)

3.1.7 Função de Distribuição Acumulada Empírica

\[ F_n(x)=\begin{cases} 0, & \text{se}\;\; x<x_{(1)},\\ \frac{i}{n}, & \text{se}\;\; x_{(i)}\leq x < x_{(i+1)}, \;\;\forall\;i: 1 \leq i<n,\\ 1, & \text{se}\;\; x\geq x_{(n)}, \end{cases}, \]

em que \(x_{(1)}, x_{(2)},\cdots, x_{(n)}\) são realizações ordenadas das variáveis aleatórias \(X_1,X_2,\cdots,X_n\).

Alternativamente, temos

\[ F_n(x)=\sum\limits_{i=1}^{n}\frac{I_{\{x_i\leq x\}}(x_i)}{n}. \]

## Grafico da Funcao de Distribuicao Empirica

n<-100:130
l1<-rep(0.95,4)
l2<-c(0.90,0.92,0.94,0.95)


par(mfrow=c(2,2))
for (i in 1:4) {
 
   plot.ecdf(Custo.est(n,m=100,l1[i],l2[i]),main="Distribuição Empírica do Custo Estimado")

   legend("topleft", legend=bquote(lambda[1]==.(l1[[i]]) ~ "," ~lambda[2]==.(l2[[i]])) )
}

4 Tabelas para o Custo de Transporte

#install.packages("knitr")
library(knitr)
n<-100:110
C_est<-Custo.est(n,m=100,0.95,0.90)

C_Teo<-Custo_teorico(n,0.95,0.90)

Int<-sapply(n,Intervalo,m=100,l1=0.95,l2=0.90)

LI<-Int[1,]
LS<-Int[2,]

dados1<-data.frame(n=n,
                  Custo_Estimado=C_est,
                  Custo_Esperado=C_Teo,
                  Limite_Inferior=LI,
                  Limite_Superior=LS
                 )

kable(dados1,digits = 2,caption=" Tabela para o Custo de Transporte (Taxa 1=0.95 e Taxa 2=0.90)")
Tabela para o Custo de Transporte (Taxa 1=0.95 e Taxa 2=0.90)
n Custo_Estimado Custo_Esperado Limite_Inferior Limite_Superior
100 852.28 852.62 822.85 865.48
101 866.01 865.77 848.85 891.91
102 882.71 878.99 855.86 899.34
103 880.78 892.28 862.17 906.08
104 903.10 905.65 880.03 924.36
105 919.14 919.09 902.49 947.25
106 930.17 932.60 901.15 946.33
107 942.68 946.19 922.50 968.11
108 956.82 959.85 935.24 981.27
109 972.87 973.58 949.36 995.82
110 976.36 987.38 962.81 1009.69 
l2<-c(0.90,0.92,0.94,0.95)
C_teo<-sapply(l2,Custo_teorico,n=100,l1=0.95)


C_est<-sapply(l2,Custo.est,m=100,n=100,l1=0.95)

Int<-sapply(l2,Intervalo,n=100,m=100,l1=0.95)

LI<-Int[1,]
LS<-Int[2,]

tx1<-rep(0.95,4) ; tx2<-c(0.90,0.92,0.94,0.95)

dados2<-data.frame(Taxa_1=tx1,
                   Taxa_2=tx2,
                   Custo_Estimado=C_est,
                   Custo_Esperado=C_teo,
                   Limite_Inferior=LI,
                   Limite_Superior=LS
)

kable(dados2,digits = 2,caption=" Tabela para o Custo de Transporte (n=100 e m=100)")
Tabela para o Custo de Transporte (n=100 e m=100)
Taxa_1 Taxa_2 Custo_Estimado Custo_Esperado Limite_Inferior Limite_Superior
0.95 0.90 847.34 852.62 835.68 878.32
0.95 0.92 813.82 820.29 804.52 846.67
0.95 0.94 807.84 800.39 774.90 816.59
0.95 0.95 791.20 794.81 772.42 813.89 
set.seed(23/12/20)
l1<-seq(0.90,2,length.out = 70)
l2<-seq(0.80,1.5,length.out = 70)
Int<-mapply(Intervalo,l1,l2,MoreArgs = list(n=100,m=100))

CUSTO.est<-mapply(Custo.est,m=100,l1,l2,MoreArgs = list(n=100))

CUSTO<-mapply(Custo_teorico,l1,l2,MoreArgs = list(n=100))


dados3<-matrix(c(l1,l2,CUSTO.est,CUSTO,Int[1,],Int[2,]), ncol=6)
colnames(dados3)<-c("Taxa 1", "Taxa 2", "Custo Estimado", "Custo Esperado", "Limite Inferior", "Limite Superior")


library(DT) ; library(dplyr)


dados3 %>%
  datatable(extensions = 'Buttons',
            options = list(dom = 'Blfrtip',
                           buttons = c('copy', 'csv', 'excel', 'pdf', 'print'),
                           lengthMenu = list(c(10,25,50,-1),
                                             c(10,25,50,"All"))),
            caption = "Tabela para o Custo de Transporte (n=100 e m=100)")%>%
                                                           formatRound(columns=all(),digits=3)

Bibliografia

Casella, George, and Roger L Berger. 2002. Statistical Inference. Vol. 2. Duxbury Pacific Grove, CA.

Kranakis, Evangelos. 2014. “On the Event Distance of Poisson Processes with Applications to Sensors.” Discrete Applied Mathematics 179: 152–62.

R Core Team. 2020. R: A Language and Environment for Statistical Computing. Vienna, Austria: R Foundation for Statistical Computing. https://www.R-project.org/.

Ross, Sheldon M. 1996. Stochastic Processes. Vol. 2.

  1. Universidade de Brasília, ↩︎

  2. Universidade de Brasília, ↩︎