英文:
How can I get the intersection point of two crossing ecdfs in R?
问题
我使用以下代码制作了两个ECDF图:
ecdf1 <- ecdf(data1)
ecdf2 <- ecdf(data2)
这些图相互交叉。我需要获取它们的交点坐标。我应该如何在R中实现这一点?
英文:
I have two ecdf plots using below code:
ecdf1 <- ecdf(data1)
ecdf2 <- ecdf(data2)
These plots are crossing each other. I need to take the crossing point (intersection point) coordinates. How should I do this in R?
答案1
得分: 8
以下是代码的翻译部分:
set.seed(1)
data1 <- rnorm(50) + 1.2
data2 <- rexp(50)
现在我们使用您的代码创建两个ecdf函数:
ecdf1 <- ecdf(data1)
ecdf2 <- ecdf(data2)
如果我们绘制它们,我们会看到这些曲线在两个点处相交:一次在0和1之间,再次在略高于2处。
plot(ecdf1, col = "red")
plot(ecdf2, col = "blue", add = TRUE)
要找到这些确切的点,我们创建一个关于x的函数,该函数是两个ecdf函数在x处的差异。然后,我们使用uniroot
函数确定此差异函数等于0的位置:
diff_func <- function(x) ecdf1(x) - ecdf2(x)
root1 <- uniroot(diff_func, c(0, 1))$root # 找到较低的交点
root2 <- uniroot(diff_func, c(2, 3))$root # 找到较高的交点
我们可以检查我们的结果是否合理:
root1
#> [1] 0.1568627
root2
#> [1] 2.055556
甚至可以绘制线段以演示这些交点是正确的:
segments(root1, y0 = 0, y1 = ecdf1(root1), lty = 2)
segments(root2, y0 = 0, y1 = ecdf1(root2), lty = 2)
创建于2023-02-26,使用 reprex v2.0.2。
英文:
Lets' create a reproducible example to demonstrate:
set.seed(1)
data1 <- rnorm(50) + 1.2
data2 <- rexp(50)
Now we use your code to create the two ecdf functions:
ecdf1 <- ecdf(data1)
ecdf2 <- ecdf(data2)
If we plot them, we will see these curves intersect at two points: once between 0 and 1, and again just above 2.
plot(ecdf1, col = "red")
plot(ecdf2, col = "blue", add = TRUE)
To find these exact points, we create a function of x that is the difference between the two ecdf functions at x. We then use the function uniroot
to determine where this difference function is equal to 0:
diff_func <- function(x) ecdf1(x) - ecdf2(x)
root1 <- uniroot(diff_func, c(0, 1))$root # Finds the lower intersection
root2 <- uniroot(diff_func, c(2, 3))$root # Finds the upper intersection
We can check our results make sense:
root1
#> [1] 0.1568627
root2
#> [1] 2.055556
And even plot segments to demonstrate the intersections are correct:
segments(root1, y0 = 0, y1 = ecdf1(root1), lty = 2)
segments(root2, y0 = 0, y1 = ecdf1(root2), lty = 2)
<sup>Created on 2023-02-26 with reprex v2.0.2</sup>
通过集体智慧和协作来改善编程学习和解决问题的方式。致力于成为全球开发者共同参与的知识库,让每个人都能够通过互相帮助和分享经验来进步。
评论