如何在R中的能量图中绘制曲线?
Posted
技术标签:
【中文标题】如何在R中的能量图中绘制曲线?【英文标题】:How to draw the curves in an energy diagram in R? 【发布时间】:2015-03-12 22:00:04 【问题描述】:我编写了以下 R 脚本:
#energy diagram
x <- c(0.1, 0.3, 0.5, 0.7, 0.9 ) #chosen randomly, reaction axis
y <- c(-5.057920, -5.057859, -5.057887,-5.057674, -5.057919 ) #energy of the educt, intermediate, transtition states and product
plot(x,y, type="p",
xlim=c(0,1),
ylim=c(-5.058,-5.0575),
xlab="reaction axis",
ylab=expression(paste(E[el] ," / ",10^6," ",kJ/mol)),
xaxt="n" #hide x-axis
)
#h- and v-lines, so i can draw curves by hand
abline(v=seq(0,1,0.1),h=seq(-5.0600,-5.0500,0.00005),col="black",lty=1,lwd=1)
abline(h=c(-5.057920, -5.057859, -5.057887,-5.057674), col="blue", lty=1,lwd=0.7)
是否可以通过看起来像能量图的点绘制曲线。能量图的示例在这里:
【问题讨论】:
lines(spline(x, y))
不幸的是不合适。
所以你提供的是最小/最大 (x,y) 对,这意味着你想要通过那些点的导数为 0 的点拟合曲线。样条线可以用于此,但我不知道任何标准的 R 实现也可以让您指定导数。 bezier
包中的贝塞尔曲线可能会让您更幸运。这将允许您从最小/最大点水平添加“控制点”,这将在最小/最大点处强制执行 0 导数。
再看一点,Hmisc::bezier
似乎更容易使用。我现在没有太多时间,但如果仍然需要,我可以在一天左右的时间内找到解决方案。
没关系,我会把自己读成贝塞尔曲线。感谢您的信息:)
【参考方案1】:
可以做很多事情来简化/矢量化这段代码,但是对于一个小图来说,这很有效:
# get that data
x <- c(0.1, 0.3, 0.5, 0.7, 0.9 ) # reaction axis
y <- c(-5.057920, -5.057859, -5.057887,-5.057674, -5.057919 ) # energies
我将制作一条小贝塞尔曲线将每个点连接到下一个点——这样我们可以确保平滑线穿过数据,而不仅仅是靠近它。我会给每个点一个“控制点”来定义斜率。通过对一个点和它的控制点使用相同的 y 值,该点的斜率将为 0。我将调用该点和控制点之间的偏移量delta
。我们将从一对点开始:
library(Hmisc)
delta = 0.15
bezx = c(0.1, 0.1 + delta, 0.3 - delta, 0.3)
bezy = rep(y[1:2], each = 2)
plot(bezx, bezy, type = 'b', col = "gray80")
lines(bezier(bezx, bezy), lwd = 2, col = "firebrick4")
在这里,我用灰色绘制了点和控制点,用红色绘制了平滑线,这样我们就可以看到发生了什么。 看起来很有希望,让我们把它变成一个我们可以应用到每一对点的函数:
bezf = function(x1, x2, y1, y2, delta = 0.15)
bezier(x = c(x1, x1 + delta, x2 - delta, x2), y = c(y1, y1, y2, y2))
delta
参数你可以玩,我觉得0.1看起来还不错。
plot(x, y, xlab = "Reaction coordinate", ylab = "E", axes = F)
box(bty = "L")
axis(side = 2)
for(i in 1:(length(x) - 1))
lines(bezf(x1 = x[i], x2 = x[i + 1], y1 = y[i], y2 = y[i + 1], delta = 0.1))
您当然可以调整情节、添加标签和ablines
,就像原来的一样。 (使用for
循环和lines
命令仅绘制平滑线。)我留下点以表明我们正在通过它们,而不仅仅是靠近。
我更喜欢在 ggplot2 中绘图,如果你也这样做,你需要将数据提取到 data.frame 中:
bezlist = list()
for (i in 1:(length(x) - 1))
bezlist[[i]] = bezf(x1 = x[i], x2 = x[i + 1], y1 = y[i], y2 = y[i + 1], delta = 0.1)
xx = unlist(lapply(bezlist, FUN = '[', 'y'))
yy = unlist(lapply(bezlist, FUN = '[', 'y'))
bezdat = data.frame(react = xx, E = yy)
library(ggplot2)
ggplot(bezdat, aes(x = react, y = E)) +
geom_line() +
labs(x = "Reaction coordinate")
【讨论】:
Gregor,你这该死的天才,你为我节省了数小时令人沮丧的研究:)。几周前我了解了ggplot2
,但我仍然需要学习如何处理它。无论如何,我真的很感谢你的努力!干杯【参考方案2】:
您可以使用样条拟合。沿着能量图定义一些点,然后使用样条函数拟合它们。您提供的分数越多,您的合身性就越好。您可以查看 stats 包中的 smooth.splines 函数,了解样条拟合的一种实现。
【讨论】:
我在脚本中添加了以下命令:smoothingSpline = smooth.spline(x, y)
和 lines(smoothingSpline)
。虽然不合适。这条线是通过五个给定点定义的,它们是局部最小值和最大值。以上是关于如何在R中的能量图中绘制曲线?的主要内容,如果未能解决你的问题,请参考以下文章
R语言使用pROC包在同一图中绘制两条ROC曲线并通过假设检验检验ROC曲线的AUC或者偏AUC的差异(输出p值)
如何使用 r 中的 ROCR 包绘制 ROC 曲线,*只有分类列联表*