sigma - Calcula el área bajo una curva

area bajo la curva pdf

Question

Me gustaría calcular el área bajo una curva para hacer la integración sin definir una función como en integrate() .

Mi información se ve así:

Date Strike Volatility
2003-01-01 20 0.2
2003-01-01 30 0.3
2003-01-01 40 0.4
etc.

Tracé la plot(strike, volatility) para observar la sonrisa de volatilidad. ¿Hay alguna forma de integrar esta "curva" graficada?

Answer 1

De acuerdo, llegué un poco tarde a la fiesta pero, al repasar las respuestas, falta una solución R sencilla para el problema. Aquí va, simple y limpio:

sum(diff(x) * (head(y,-1)+tail(y,-1)))/2

La solución para OP entonces se lee como:

sum(diff(strike) * (head(volatility,-1)+tail(volatility,-1)))/2

Esto efectivamente calcula el área usando el método trapezoidal tomando el promedio de los valores y "izquierdo" y "derecho".

NB: como @Joris ya señaló que podría usar abs(y) si eso tuviera más sentido.

Answer 2

El AUC se aproxima bastante fácilmente al observar muchas figuras de trapecio, cada una de ellas ligada entre x_i , x_{i+1} , y{i+1} y y_i . Usando el paquete rollmean del zoológico, puedes hacer:

library(zoo)
x <- 1:10
y <- 3*x+25
id <- order(x)
AUC <- sum(diff(x[id])*rollmean(y[id],2))

Asegúrese de ordenar los valores x, o su resultado no tendrá sentido. Si tiene valores negativos en algún lugar a lo largo del eje y, tendrá que averiguar cómo exactamente quiere definir el área bajo la curva, y ajustar en consecuencia (por ejemplo, usando abs() )

En cuanto a su seguimiento: si no tiene una función formal, ¿cómo la trazaría? Entonces, si solo tiene valores, lo único que puede aproximar es una integral definida. Incluso si tiene la función en R, solo puede calcular integrales definidas usando integrate() . Trazar la función formal solo es posible si también puedes definirla.

Answer 3

En el mundo de la farmacocinética (PK), el cálculo de diferentes tipos de AUC es una tarea común y fundamental. Hay muchos cálculos diferentes de AUC para farmacokietics, como

• AUC0-t = AUC desde cero hasta el tiempo t
• AUC0-last = AUC desde cero hasta el último punto de tiempo (puede ser el mismo que el anterior)
• AUC0-inf = AUC desde cero hasta el infinito en el tiempo
• AUCint = AUC en un intervalo de tiempo
• AUCall = AUC durante todo el período de tiempo para el que existen datos

Uno de los mejores paquetes que hace estos cálculos es el paquete relativamente nuevo PKNCA de la gente de Pfizer. Echale un vistazo.

Answer 4

Puede usar el paquete ROCR, donde las siguientes líneas le darán el AUC:

pred <- prediction(classifier.labels, actual.labs)
attributes(performance(pred, ''auc''))$y.values[[1]]

Answer 5

Simplemente agregue lo siguiente a su programa y obtendrá el área debajo de la curva:

require(pracma)
AUC = trapz(strike,volatility)

De ?trapz :

Este enfoque coincide exactamente con la aproximación para integrar la función usando la regla trapezoidal con los puntos base x.

Answer 6

Tres opciones más, incluyendo una que usa un método spline y otra que usa la regla de Simpson ...

# get data
n <- 100
mean <- 50
sd <- 50
x <- seq(20, 80, length=n)
y <- dnorm(x, mean, sd) *100
# using sintegral in Bolstad2
require(Bolstad2)
sintegral(x,y)$int
# using auc in MESS
require(MESS)
auc(x,y, type = ''spline'')
# using integrate.xy in sfsmisc
require(sfsmisc)
integrate.xy(x,y)

El método trapezoidal es menos preciso que el método spline, por lo que probablemente se prefiera MESS::auc (usa el método spline) o Bolstad2::sintegral (usa la regla de Simpson). Las versiones de bricolaje de estos (y un enfoque adicional usando la regla de cuadratura) están aquí: http://www.r-bloggers.com/one-dimensional-integrals/

Answer 7

La respuesta de Joris Meys fue genial, pero tuve problemas para eliminar NA de mis muestras. Aquí está la pequeña función que escribí para tratar con ellos:

library(zoo) #for the rollmean function
######
#'' Calculate the Area Under Curve of y~x
#''
#''@param y Your y values (measures ?)
#''@param x Your x values (time ?)
#''@param start : The first x value
#''@param stop : The last x value
#''@param na.stop : returns NA if one value is NA
#''@param ex.na.stop : returns NA if the first or the last value is NA
#''
#''@examples
#''myX = 1:5
#''myY = c(17, 25, NA, 35, 56)
#''auc(myY, myX)
#''auc(myY, myX, na.stop=TRUE)
#''myY = c(17, 25, 28, 35, NA)
#''auc(myY, myX, ex.na.stop=FALSE)
auc = function(y, x, start=first(x), stop=last(x), na.stop=FALSE, ex.na.stop=TRUE){
if(all(is.na(y))) return(NA)
bounds = which(x==start):which(x==stop)
x=x[bounds]
y=y[bounds]
r = which(is.na(y))
if(length(r)>0){
if(na.stop==TRUE) return(NA)
if(ex.na.stop==TRUE & (is.na(first(y)) | is.na(last(y)))) return(NA)
if(is.na(last(y))) warning("Last value is NA, so this AUC is bad and you should feel bad", call. = FALSE)
if(is.na(first(y))) warning("First value is NA, so this AUC is bad and you should feel bad", call. = FALSE)
x = x[-r]
y = y[-r]
}
sum(diff(x[order(x)])*rollmean(y[order(x)],2))
}

Luego lo uso con una aplicación en mi marco de datos: myDF$auc = apply(myDF, MARGIN=1, FUN=auc, x=c(0,5,10,15,20))

Espero que pueda ayudar a noobs como yo :-)

EDITAR: límites agregados