• Escrevendo suas próprias funções

• Explorando funções para cálculo de Seguros

• É muito simples escrever sua própria função no R.

• Por exemplo, considere a seguinte função:

\[S = \sum_{i=1}^n \frac{p_i \cdot x_i}{(1+d_i)^i} \qquad \mbox{para } x_i \in \mathbb{R}; \,\, p_i \in [0,1]; \,\, d_i \in \mathbb{R}_{+} \]

# Exemplo - função personalizada

f = function(x,p,d){
  n = length(x)
  s = sum( p*x/(1+d)^(1:n) )
  return(s)
}

• Se digitarmos f no console:

f

## function(x,p,d){
##   n = length(x)
##   s = sum( p*x/(1+d)^(1:n) )
##   return(s)
## }

• Para avaliar a função com alguns dados:

f(x=c(100,200,100), p=c(.4,.5,.3), d=0.05)

## [1] 154.7133

• Se executarmos f( c(100,200,100), c(.4,.5,.3), 0.05), sem o nome dos argumentos, vamos ter o mesmo resultado acima. O R assume que os argumentos estão na ordem em que foram definidos na função.

• Também podemos definir um valor default para algum argumento da função usando:

f = function(x,p,d=0.05)

• Nesse caso, a função pode ser avaliada sem que um outro valor seja especificado para d, e o R irá usar o valor default.

• Essas propriedades são similares para funções do R.

• Por exemplo, a função qnorm() computa quantis para a distribuição normal.

qnorm(.95)  # quantil 95% para normal(0,1)

## [1] 1.644854

qnorm(.95, mean=1, sd=2)  # quantil 95% para normal(1,4)

## [1] 4.289707

• Também é possível usar funções recursivamente. Por exemplo, considere a função:

\[h(n) = \begin{cases} 2\,.h(n-1) + 1 & \mbox{se } n \geq 2 \\ 1 & \mbox{se } n=1 \end{cases} \]

h = function(n){
  if(n<=1) 
    return(1) 
  else
    return(2*h(n-1)+1)
}
h(4)

## [1] 15

• A função pode ter variáveis locais com o mesmo nome de variáveis globais.

beta = 0  # criando uma variavel global

slope = function(X,Y){  # criando a funcao 'slope'
  beta = coefficients(lm(Y~X))[2]  # variavel local
  return(as.numeric(beta))
}

attach(cars)
slope(speed,dist)
## [1] 3.932409

beta  # continua com o valor original
## [1] 0

• Podemos usar a função cat() para imprimir comentários ou valores:

slope = function(X,Y){
  beta = coefficients(lm(Y~X))[2]
  cat("A inclinação da reta é ",beta,"\n")
  return(as.numeric(beta))
}
slope(speed,dist)

## A inclinação da reta é  3.932409

## [1] 3.932409

• Vamos usar uma tabela de vida da França do pacote lifecontingencies:

# Funções para cálculo de seguros
require(lifecontingencies)

## Warning: package 'lifecontingencies' was built under R version 4.1.3

data("demoFrance")

• Vamos obter o número de pessoas vivas na idade \(x\):

# carregando a tabela de vida francesa
alive = demoFrance$TV88_90
head(alive)   # l_x

## [1] 100000  99352  99294  99261  99236  99214

• Para construir o vetor com o número de mortes a cada idade:

death = -diff(alive)  # d_x
head(death)

## [1] 648  58  33  25  22  20

• Função de Mortalidade de Makeham (função de sobrevivência):

\[S(x) = \exp \left( -ax - \frac{b}{\log c} [c^x -1] \right) \qquad \forall x \geq 0 \] para alguns parâmetros \(a \geq 0\), \(b \geq 0\) e \(c > 1\).

# função de sobrevivência de Makeham
sMakeham = function(x,a,b,c){
  ifelse( x<0, 1, exp(-a*x-b/log(c)*(c^x-1)) )
}

• Função de probabilidade para idades inteiras:

# função de probabilidade para idades inteiras
dMakeham = function(x,a,b,c){
  ifelse( x>floor(x), 0, sMakeham(x,a,b,c)-sMakeham(x+1,a,b,c) )
}

\[S(x) - S(x+1) = P(T_0>x) - P(T_0>x+1) = P(x<T_0 \leq x+1) \]

• Podemos usar essa função dMakeham para estimar os parâmetros \(a\), \(b\) e \(c\) usando Máxima Verossimilhança.

  (para isso, precisamos remover as mortes ao nascer, e as mortes acima de 105 anos, por restrições da lei de Makeham)

death = death[-c(1,107:111)]  # removendo algumas idades
ages = 1:(length(death))

# função da log verossimilhança
loglikMakeham = function(abc){  # abc: vetor com paramêtros
  - sum( log( dMakeham(ages,abc[1],abc[2],abc[3]))*death[ages] )
}

• Podemos usar a função optim() para obter estimadores de máxima verossimilhança para os parâmetros da função de mortalidade de Makeham.

# estimador de máxima verossimilhança
mlEstim = optim( c(1e-5,1e-4,1.1), loglikMakeham)
abcml = mlEstim$par
abcml

## [1] 4.022473e-04 4.366993e-06 1.123549e+00

• Baseado nas idades observadas de morte, podemos calcular a idade esperada de morte:

sum( (ages+.5)*death )/sum(death)

## [1] 81.19984

• Comparando com o estimado de acordo com o modelo de Makeham:

integrate( sMakeham, 0, Inf, abcml[1], abcml[2], abcml[3])

## 81.15647 with absolute error < 9.2e-05

• Lidando com erros (Seção 1.3.4)

• Funções eficientes (Seção 1.3.5)

• Integração Numérica (Seção 1.3.6)