как я могу рассчитать 95% достоверный предел площади под кривой в R?

У меня следующая раздача:

x<-c(22.5,28.14285714,33.78571429,39.42857143,45.07142857,50.71428571,56.35714286,62,67.64285714,73.28571429,78.92857143,84.57142857,90.21428571,95.85714286,101.5,107.1428571,112.7857143,118.4285714,124.0714286,129.7142857,135.3571429,141,146.6428571,152.2857143,157.9285714,163.5714286,169.2142857,174.8571429,180.5,186.1428571,191.7857143,197.4285714,203.0714286,208.7142857,214.3571429,220,225.6428571,231.2857143,236.9285714,242.5714286,248.2142857,253.8571429,259.5,265.1428571,270.7857143,276.4285714,282.0714286,287.7142857,293.3571429,299)
y<-c(0.003541755,0.003541755,0.003541755,0.003541755,0.003541755,0.003541755,0.003541755,0.003541755,0.003541755,0.003541755,0.003541755,0.003541755,0.003541755,0.003541755,0.003541755,0.003541755,0.003541755,0.003541755,0.003541755,0.00328839614285714,0.00296425985714286,0.002655899,0.00236187857142857,0.002080895,0.00181184271428571,0.00155376085714286,0.00130578928571429,0.001074706,0.000877193,0.000709397142857142,0.000567189714285714,0.000447254,0.000346858571428571,0.000263689142857143,0.000195768428571429,0.000141427,9.92657142857141e-05,6.77857142857142e-05,4.48571428571428e-05,2.86428571428571e-05,1.75142857142857e-05,1.01357142857143e-05,5.52e-06,2.78857142857142e-06,1.27285714285713e-06,5.00714285714284e-07,1.5742857142857e-07,3.29857142857142e-08,2.78857142857137e-09,1.74e-12)

plot(x,y)

Я хотел бы найти значение x, которое разделяет область 0,95 под распределением слева и 0,05 области справа (односторонний 95% интервал достоверности).

Я предполагаю, что мне нужно подогнать свою эмпирическую кривую к функции, а затем интегрировать функцию, чтобы получить желаемое значение, но я действительно не знаю, с чего начать.

Как это можно сделать в R?


r interpolation integral
person user18441    schedule 07.01.2013    source источник
comment
Ответ GSee - это правильный путь. Но я хотел бы отметить, что численное интегрирование ваших исходных данных не только проще, чем создание функции соответствия и интегрирование, но и будет иметь меньше ошибок вычислений (в целом).   -  person Carl Witthoft    schedule 07.01.2013
comment
@CarlWitthoft, я не совсем уверен в своем ответе (который был quantile(x, 0.95)). Он разделяет xs на 95% и 5%, но не принимает во внимание площадь (ys).   -  person GSee    schedule 07.01.2013
comment
@Gsee - Полагаю, вам все-таки следует сгенерировать квантиль интеграла Симпсона. Я бы все равно не сгенерировал подходящую функцию.   -  person Carl Witthoft    schedule 08.01.2013

Ответы (2)


arrow_upward
4
arrow_downward

Как указывали другие ответы, это интеграция под проблемой кривой, в сочетании с определением, где площадь достигает 95% от общей площади. Я использую более простой подход к интеграции, чем ответ Дэвида. Вместо того, чтобы интерполировать кривую и интегрировать ее, я просто использую правила интегрирования трапеций, чтобы получить площадь, вносимую каждым интервалом. Затем эти отдельные области складываются, чтобы получить общую площадь. Затем определяется индекс, в котором совокупная площадь превышает 95% от общей площади, и по нему можно провести линию.

piece_area <- c(0, (x[-1] - x[-length(x)])*(y[-1] + y[-length(y)]) / 2)
cum_area <- cumsum(piece_area)
total_area <- cum_area[length(cum_area)]
idx095 <- min(which(cum_area > 0.95 * total_area))

abline(v = x[idx095])

введите описание изображения здесь

Более высокое разрешение точной точки, в которой пересекается 95%, можно получить, используя больше точек в исходной выборке распределения.

person Brian Diggs    schedule 07.01.2013

arrow_upward
2
arrow_downward

Это задача интегрирования (сумма под кривой). Вы можете разделить вашу интеграцию на квадратную + кривую. Однако вы можете использовать быстрое и грязное приближение с помощью сплайнов:

y<-c(0.003541755,0.003541755,0.003541755,0.003541755,0.003541755,0.003541755,0.003541755,0.003541755,0.003541755,0.003541755,0.003541755,0.003541755,0.003541755,0.003541755,0.003541755,0.003541755,0.003541755,0.003541755,0.003541755,0.00328839614285714,0.00296425985714286,0.002655899,0.00236187857142857,0.002080895,0.00181184271428571,0.00155376085714286,0.00130578928571429,0.001074706,0.000877193,0.000709397142857142,0.000567189714285714,0.000447254,0.000346858571428571,0.000263689142857143,0.000195768428571429,0.000141427,9.92657142857141e-05,6.77857142857142e-05,4.48571428571428e-05,2.86428571428571e-05,1.75142857142857e-05,1.01357142857143e-05,5.52e-06,2.78857142857142e-06,1.27285714285713e-06,5.00714285714284e-07,1.5742857142857e-07,3.29857142857142e-08,2.78857142857137e-09,1.74e-12)
x<-c(22.5,28.14285714,33.78571429,39.42857143,45.07142857,50.71428571,56.35714286,62,67.64285714,73.28571429,78.92857143,84.57142857,90.21428571,95.85714286,101.5,107.1428571,112.7857143,118.4285714,124.0714286,129.7142857,135.3571429,141,146.6428571,152.2857143,157.9285714,163.5714286,169.2142857,174.8571429,180.5,186.1428571,191.7857143,197.4285714,203.0714286,208.7142857,214.3571429,220,225.6428571,231.2857143,236.9285714,242.5714286,248.2142857,253.8571429,259.5,265.1428571,270.7857143,276.4285714,282.0714286,287.7142857,293.3571429,299)

sp=smooth.spline(x,y)
f = function(t)
{
    predict(sp,t)$y
}   

N=500 # this is an accuracy parameter
xBis=seq(x[1],x[length(x)],length=N)
yBis=sapply(x,f)

J = function (input)
{   # This function takes input in 1:N
    Integral = 0
    dx=(x[length(x)]-x[1])/N

    for ( j in 1: input)
{   z=xBis[j]
    Integral=Integral+ f(x[1]+z)*dx
}
J=Integral
}
######
I=J(N) # This is the value of the sum under the curve
# It should be roughly equal (given the shape of the curve) to:
index=max(which(y==y[1]))
I = (x[index]-x[1])*(y[index])*3/2
######
res=sapply(1:N,J)/I
Index5=max(which(res<=.05))
Index95=min(which(res>=.95))

x5=xBis[Index5] # This is the 5% quantile 
x95=xBis[Index95]

HTH

Сообщите мне, если что-то неясно

P.S Я думаю, что есть гораздо лучшие способы сделать это ...

person DKK    schedule 07.01.2013