2013-06-07 3 views
7

Хотя я думаю, что это основной вопрос, я не могу показаться, чтобы выяснить, как вычислить это в R:точка пересечения 2 нормальных кривых

точка пересечения (мне нужно х-значение) 2 или более нормальные распределения (установленные на гистограмме), которые имеют, например, следующие параметры:

d=data.frame(mod=c(1,2),mean=c(14,16),sd=c(0.9,0.6),prop=c(0.6,0.4)) 

с среднего и стандартным отклонением моих 2 кривых, а также проп пропорции вклада каждых мод в распределение.

ответ

12

Вы можете использовать uniroot:

f <- function(x) dnorm(x, m=14, sd=0.9) * .6 - dnorm(x, m=16, sd=0.6) * .4 

uniroot(f, interval=c(12, 16)) 

$root 
[1] 15.19999 

$f.root 
[1] 2.557858e-06 

$iter 
[1] 5 

$estim.prec 
[1] 6.103516e-05 


ETA некоторых экспозиции:

uniroot является одномерным корневым искателем, т.е. заданной функцией f одной переменной x, он находит значение x что решает уравнение f(x) = 0.

Чтобы использовать его, вы предоставляете функцию f вместе с интервалом, в пределах которого предполагается значение решения. В этом случае f - это просто разница между двумя плотностями; точка, где они пересекаются, будет равна f. В этом примере я получил интервал (12, 16), создав сюжет и увидев, что они пересекаются вокруг x = 15.

+0

+1, но вы можете добавить некоторое объяснение того, что это делает/как оно работает? Спасибо –

+0

Спасибо за текст. Отлично! –

+0

спасибо, работает отлично !! – Wave

0

Извините, но принятый ответ не является хорошим. Смотрите также: intersection of two curves in matlab

Вы можете получить как корни, используя такую ​​функцию:

intersect <- function(m1, s1, m2, s2, prop1, prop2){ 

B <- (m1/s1^2 - m2/s2^2) 
A <- 0.5*(1/s2^2 - 1/s1^2) 
C <- 0.5*(m2^2/s2^2 - m1^2/s1^2) - log((s1/s2)*(prop2/prop1)) 

(-B + c(1,-1)*sqrt(B^2 - 4*A*C))/(2*A) 
} 

в вашем случае:

> intersect(14,0.9,16,0.6,0.6,0.4) 
[1] 20.0 15.2 
+0

Как вы выбираете наиболее важный пункт? – Spidfire

0

Я согласен с @Flounderer, что вы должны рассчитать как корни, но предлагаемая функция неполна. Когда s1 = s2, A становится нулевым, а Inf и NaN производятся.

Незначительная модификация завершает функцию:

intersect <- function(m1, sd1, m2, sd2, p1, p2){ 

    B <- (m1/sd1^2 - m2/sd2^2) 
    A <- 0.5*(1/sd2^2 - 1/sd1^2) 
    C <- 0.5*(m2^2/sd2^2 - m1^2/sd1^2) - log((sd1/sd2)*(p2/p1)) 

    if (A!=0){ 
    (-B + c(1,-1)*sqrt(B^2 - 4*A*C))/(2*A) 
    } else {-C/B} 
} 
m1=0; sd1=2; m2=2.5; sd2=2; p1=.8; p2=.2 
(is=intersect(m1,sd1,m2,sd2,p1,p2)) 

xs = seq(-6, 6, by=.01) 
plot(xs, p1*dnorm(xs, m1, sd1), type= 'l') 
lines(xs, .2*dnorm(xs, m2,sd2)) 
abline(v=is) 

Общее решение также найти с помощью uniroot.all:

library(rootSolve) 
f <- function(x, m1, sd1, m2, sd2, p1, p2){ 
    dnorm(x, m1, sd1) * p1 - dnorm(x, m2, sd2) * p2 } 
uniroot.all(f, lower=-6, upper=6, m1=m1, sd1=sd1, m2=m2, sd2=sd2, p1=p1, p2=p2) 
Смежные вопросы