rg255
rg255

Reputation: 4169

Plot line of a function in R

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),")")))

enter image description here 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)

enter image description here

Upvotes: 1

Views: 243

Answers (2)

nicola
nicola

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

Paul Regular
Paul Regular

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

Related Questions