Reputation: 4169
I have the following script:
FGM = function (n,r,z){
x = r*sqrt(n)/(2*z)
Px = 1-pnorm(x)
}
re = 10000
data = data.frame(abs(rnorm(re,0,1)), abs(rnorm(re,0,1)), abs(rnorm(re,0,1)))
colnames(data) = c("n","r","z")
data$Px = FGM(data$n,data$r,data$z)
data$x = data$r*sqrt(data$n)/(2*data$z)
par(mar=c(4.5,4.5,1,1))
plot(data$x,data$Px, xlim = c(0,3), pch = 19, cex = 0.1, xaxs="i", yaxs="i",
xlab = expression(paste("Standardized mutational size (",italic(x), ")")),
ylab = expression(paste("P"[a],"(",italic(x),")")))
Which is a recreation of the graph found here (box 2). You can see in this script that I do this by just plotting 10000 small black points with various values of n,z, and r. This seems like an ugly work around, I think I should just be able to give R my function
FGM = function (n,r,z){
x = r*sqrt(n)/(2*z)
Px = 1-pnorm(x)
}
and have it plot a line on a graph. However, a few hours of scouring the web has been unproductive, and I tried a few ways with abline and lines but nothing worked, is there a way of doing it with these functions or another function?
Tried this...
plot(data$x,data$Px, xlim = c(0,3), ylim = c(0,0.5), xaxs="i", yaxs="i",
xlab = expression(paste("Standardized mutational size (",italic(x), ")")),
ylab = expression(paste("P"[a],"(",italic(x),")")), type = "n")
curve(1-pnorm(r*sqrt(n)/(2*z)), add=T)
>Error in curve(1 - pnorm(r * sqrt(n)/(2 * z)), add = T) :
'expr' must be a function, or a call or an expression containing 'x'
>
@PaulRegular offered this solution but it still plots based on data, not the formula itself. I'm looking for a solution which can produce the curve properly without large values of "re" - using the following but with "re" set to 10 you can see what I mean...
data <- data[order(data$x),]
lines(data$x, data$Px, lwd=1)
Upvotes: 1
Views: 243
Reputation: 24480
You can pass a function of just one variable to plot
. I guess that you are looking for:
plot(function(x) 1-pnorm(x),0,3)
Upvotes: 3
Reputation: 518
Try sorting your data by x, then add the line:
data <- data[order(data$x),]
lines(data$x, data$Px, lwd=2)
Upvotes: 0