2013-06-07 42 views
7

Mặc dù tôi cho rằng đây là câu hỏi cơ bản, tôi dường như không thể tìm ra cách tính toán điều này trong R:điểm giao nhau 2 đường cong bình thường

điểm giao nhau (Tôi cần giá trị x) 2 hoặc nhiều phân phối chuẩn (lắp trên một biểu đồ) có ví dụ các thông số sau:

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

với giá trị trung bình và độ lệch chuẩn của 2 đường cong của tôi, và chống đỡ các tỷ lệ đóng góp của từng mod để phân phối.

Trả lời

12

Bạn có thể sử dụng 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 một số triển lãm:

uniroot là một công cụ tìm gốc đơn biến, tức là cung cấp một chức năng f của một biến x, nó tìm thấy giá trị của x rằng giải phương trình f(x) = 0.

Để sử dụng nó, bạn cung cấp hàm f, cùng với khoảng thời gian trong đó giá trị giải pháp được giả định nằm. Trong trường hợp này, f chỉ là sự khác biệt giữa hai mật độ; điểm mà họ giao nhau sẽ là nơi f là số không. Tôi nhận được khoảng thời gian (12, 16) trong ví dụ này bằng cách vẽ một cốt truyện và thấy rằng chúng giao nhau quanh x = 15.

+0

+1, nhưng bạn có thể thêm một số giải thích về điều này không/bao nó hoạt động? Cảm ơn –

+0

Cảm ơn văn bản. Điều đó thật tuyệt! –

+0

cảm ơn, hoạt động hoàn hảo !! – Wave

0

Xin lỗi, nhưng câu trả lời được chấp nhận là không tốt. Xem thêm: intersection of two curves in matlab

Bạn có thể lấy cả rễ sử dụng một chức năng như thế này:

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) 
} 

trong trường hợp của bạn:

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

Làm thế nào để bạn chọn điểm quan trọng nhất? – Spidfire

0

Tôi đồng ý với @Flounderer mà bạn nên tính toán cả gốc rễ, nhưng chức năng được cung cấp không đầy đủ. Khi s1 = s2, A trở thành 0 và Inf và NaN được tạo ra.

Một sửa đổi chút ít hoàn thành chức năng:

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) 

Một giải pháp chung cũng được tìm thấy sử dụng 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) 
Các vấn đề liên quan