SLR computation

Calcul du SLR

Génération de quelques données sous-jaçentes

Nous générons des scores respectivement selon les distributions sous Hp et Hd. Pour cela, nous avons pris les estimateurs d’un cas antérieur (155.pdf).

n_within <- 722
shape <- 5.01261
scale <- 14368.05242

n_between <- 10000 # 935550 on peut réduire cette valeur pour gagner un peu de temps (ici 10000 par exemple).
meanlog <- 7.52308
sdlog <- 0.22265

Scores_Hp <- rweibull(n_within,shape,scale)
Scores_Hd <- rlnorm(n_between, meanlog, sdlog)

Estimation des paramètres des distributions

Nous estimons les paramètres des deux distributions sous-jaçentes. Ils seront différents des estimateurs de départ car nous avons un nombre de points limités pour les estimer, respectivement 722 pour la distribution Weibull et 10^{4} pour la distribution lognormale.

#Two required libraries
library(tidyverse)
library(fitdistrplus)

#Within source distribution according to a Weibull
fw <- fitdist(Scores_Hp, "weibull")
summary(fw)

Fitting of the distribution ' weibull ' by maximum likelihood 
Parameters : 
          estimate  Std. Error
shape     5.168741   0.1496142
scale 14327.044435 108.7821517
Loglikelihood:  -6778.191   AIC:  13560.38   BIC:  13569.55 
Correlation matrix:
          shape     scale
shape 1.0000000 0.3163851
scale 0.3163851 1.0000000

plot(fw)

#differences with initial estimates for Weibull
#shape
fw$estimate[1]-shape

    shape 
0.1561312 

#scale
fw$estimate[2]-scale

    scale 
-41.00799 

#Within source distribution according to a LogNormal
fb <- fitdist(Scores_Hd, "lnorm")
summary(fb)

Fitting of the distribution ' lnorm ' by maximum likelihood 
Parameters : 
         estimate  Std. Error
meanlog 7.5225612 0.002227449
sdlog   0.2227449 0.001574901
Loglikelihood:  -74397.71   AIC:  148799.4   BIC:  148813.9 
Correlation matrix:
        meanlog sdlog
meanlog       1     0
sdlog         0     1

plot(fb)

#differences with initial estimates for LogNormal
#shape
fb$estimate[1]-meanlog

      meanlog 
-0.0005187912 

#scale
fb$estimate[2]-sdlog

       sdlog 
9.486877e-05 

Plot des deux distributions selon les estimateurs obtenus et indication du score de la transaction d’intérêt

Pour obtenir un plot analogue à la sortie de XENA, nous construisons des vecteurs sur la base des deux distributions pour les placer sur un même graphe. Nous allons ajouter une ligne verticale représentant le score (evidence_score) obtenu dans le cas d’espère entre la trace et l’empreinte.

# Generate sample data
npoints <- 1000000 #size that allows to obtain a smooth density on the plot
evidence_score <- 11645

weibull_data <- rweibull(npoints, shape = fw$estimate[1], scale = fw$estimate[2])
lognorm_data <- rlnorm(npoints, meanlog = fb$estimate[1], sdlog = fb$estimate[2])

# Create a data frame
data <- data.frame(
  value = c(weibull_data, lognorm_data),
  distribution = factor(c(rep("Weibull", length(weibull_data)), 
                          rep("Lognormal", length(lognorm_data))))
)

# Plot and add evidence score
plotXENA <- ggplot(data, aes(x = value, fill = distribution)) +
  geom_density(alpha = 0.5, linewidth = 1) +
  scale_fill_manual(values = c("Weibull" = "red", "Lognormal" = "blue")) +
  labs(title = "Density Plots for Weibull and Lognormal Distributions",
       x = "Value",
       y = "Density") +
  geom_vline(xintercept = evidence_score, color = "black", linewidth = 1) +
  theme(legend.position = "bottom")  # Positioning legend at the bottom

plotXENA

Calcul du SLR raw

A partir du evidence_score et des paramètres des deux distributions, nous pouvons simplement obtenir le SLR-raw, en obtenant la densité de probabilité au point du *evidence_score” d’intérêt. A noter que les paramètres des deux distributions sont obtenus des estimateurs que nous avons calculés plus haut sur la base des scores originellement générés.

#Get the density under Hp
density_evidence_score_Hp <- dweibull(evidence_score, 
                                      shape = fw$estimate[1], 
                                      scale = fw$estimate[2])
#Get the density under Hd
density_evidence_score_Hd <- dlnorm(evidence_score, 
                                      meanlog = fb$estimate[1], 
                                      sdlog = fb$estimate[2])

#Get the SLR
SLR_raw <- density_evidence_score_Hp/density_evidence_score_Hd
#SLR
log10(SLR_raw)

[1] 14.66487

Calibration

Une régression logistique est utilisée pour calibrer le SLR_raw. Le modèle de calibration est basé sur un modèle dont les paramètres sont contenus dans l’objet: calfunc.Rdata utilisé par XENA. Le code plus bas exploite les paramètres de ce modèle de régression.

#If regression.Rdata is available, we load the regression model
#load("calfunc.Rdata")

#Without the regression.Rdata, we simply put the intercept and the slope in a vector
calfunc <- c(0.4998826, 0.8462218)

#Compute calibrated SLR by applied the logistic regression (in log10)
log10_SLR_calibrated <- as.numeric(log10(SLR_raw) * calfunc[2] + calfunc[1])
SLR_calibrated <- 10^log10_SLR_calibrated

log10_SLR_calibrated

[1] 12.90961

SLR_calibrated

[1] 8.121095e+12