Hello, You are not closing the pdf device. The only changes I have made to your code are right at the beginning of the plotting instructions and at the end of the code. ## The rest of the code is for plotting the image pdf(file = "power.pdf") op <- par(mfrow = c(4,2), cex = 0.45) [...] par(op) dev.off() ################# The comments only line is your last code line. The result is attached. Hope this helps, Rui Barradas ?s 19:39 de 09/05/21, varin sacha via R-help escreveu:> Dear R-experts, > > I am trying to get the 8 graphs like the ones in this paper : > https://statweb.stanford.edu/~tibs/reshef/comment.pdf > My R code does not show any error message neither warnings but I d'on't get what I would like to get (I mean the 8 graphs), so I am missing something. What's it ? Many thanks for your precious help. > > ################# > set.seed(1) > library(energy) > > # Here we define parameters which we use to simulate the data > # The number of null datasets we use to estimate our rejection reject #regions for an alternative with level 0.05 > nsim=50 > > # Number of alternative datasets we use to estimate our power > nsim2=50 > > # The number of different noise levels used > num.noise <- 30 > > # A constant to determine the amount of noise > noise <- 3 > > # Number of data points per simulation > n=100 > > # Vectors holding the null "correlations" (for pearson, for spearman, for kendall and dcor respectively) for each # of the nsim null datasets at a #given noise level > val.cor=val.cors=val.cork=val.dcor=rep(NA,nsim) > > # Vectors holding the alternative "correlations" (for pearson, for #spearman, for kendall and dcor respectively) #for each of the nsim2 alternative datasets at a given noise level > val.cor2=val.cors2=val.cork2=val.dcor2= rep(NA,nsim2) > > > # Arrays holding the estimated power for each of the 4 "correlation" types, for each data type (linear, #parabolic, etc...) with each noise level > power.cor=power.cors=power.cork=power.dcor= array(NA, c(8,num.noise)) > > ## We loop through the noise level and functional form; each time we #estimate a null distribution based on #the marginals of the data, and then #use that null distribution to estimate power > ## We use a uniformly distributed x, because in the original paper the #authors used the same > > for(l in 1:num.noise) { > > ????? for(typ in 1:8) { > > ## This next loop simulates data under the null with the correct marginals (x is uniform, and y is a function of a #uniform with gaussian noise) > > ??? for(ii in 1:nsim) { > ????? x=runif(n) > > #lin+noise > if(typ==1) { > y=x+ noise *(l/num.noise)* rnorm(n) > } > > #parabolic+noise > if(typ==2) { > y=4*(x-.5)^2+? noise * (l/num.noise) * rnorm(n) > } > > #cubic+noise > if(typ==3) { > y=128*(x-1/3)^3-48*(x-1/3)^3-12*(x-1/3)+10* noise? * (l/num.noise) *rnorm(n) > } > > #sin+noise > if(typ==4) { > y=sin(4*pi*x) + 2*noise * (l/num.noise) *rnorm(n) > } > > #their sine + noise > if(typ==5) { > y=sin(16*pi*x) + noise * (l/num.noise) *rnorm(n) > } > > #x^(1/4) + noise > if(typ==6) { > y=x^(1/4) + noise * (l/num.noise) *rnorm(n) > } > > #circle > if(typ==7) { > y=(2*rbinom(n,1,0.5)-1) * (sqrt(1 - (2*x - 1)^2)) + noise/4*l/num.noise *rnorm(n) > } > > #step function > if(typ==8) { > y = (x > 0.5) + noise*5*l/num.noise *rnorm(n) > } > > # We resimulate x so that we have the null scenario > x <- runif(n) > > # Calculate the 4 correlations > val.cor[ii]=(cor(x,y)) > val.cors[ii]=(cor(x,y,method=c("spearman"))) > val.cork[ii]=(cor(x,y,method=c("kendal"))) > val.dcor[ii]=dcor(x,y) > } > > ## Next we calculate our 4 rejection cutoffs > cut.cor=quantile(val.cor,.95) > cut.cors=quantile(val.cors,.95) > cut.cork=quantile(val.cork,.95) > cut.dcor=quantile(val.dcor,.95) > > ## Next we simulate the data again, this time under the alternative > > ??? for(ii in 1:nsim2) { > ????? x=runif(n) > > #lin+noise > if(typ==1) { > y=x+ noise *(l/num.noise)* rnorm(n) > } > > #parabolic+noise > if(typ==2) { > y=4*(x-.5)^2+? noise * (l/num.noise) * rnorm(n) > } > > #cubic+noise > if(typ==3) { > y=128*(x-1/3)^3-48*(x-1/3)^3-12*(x-1/3)+10* noise? * (l/num.noise) *rnorm(n) > } > > #sin+noise > if(typ==4) { > y=sin(4*pi*x) + 2*noise * (l/num.noise) *rnorm(n) > } > > #their sine + noise > if(typ==5) { > y=sin(16*pi*x) + noise * (l/num.noise) *rnorm(n) > } > > #x^(1/4) + noise > if(typ==6) { > y=x^(1/4) + noise * (l/num.noise) *rnorm(n) > } > > #circle > if(typ==7) { > y=(2*rbinom(n,1,0.5)-1) * (sqrt(1 - (2*x - 1)^2)) + noise/4*l/num.noise *rnorm(n) > } > > #step function > if(typ==8) { > y = (x > 0.5) + noise*5*l/num.noise *rnorm(n) > } > > ## We again calculate our 4 "correlations" > val.cor2[ii]=(cor(x,y)) > val.cors2[ii]=(cor(x,y,method=c("spearman"))) > val.cork2[ii]=(cor(x,y,method=c("kendal"))) > val.dcor2[ii]=dcor(x,y) > } > > ## Now we estimate the power as the number of alternative statistics #exceeding our estimated cutoffs > power.cor[typ,l] <- sum(val.cor2 > cut.cor)/nsim2 > power.cors[typ,l] <- sum(val.cors2 > cut.cor)/nsim2 > power.cork[typ,l] <- sum(val.cork2 > cut.cor)/nsim2 > power.dcor[typ,l] <- sum(val.dcor2 > cut.dcor)/nsim2 > } > } > > save.image() > > ## The rest of the code is for plotting the image > pdf("power.pdf") > par(mfrow = c(4,2), cex = 0.45) > plot((1:30)/10, power.cor[1,], ylim = c(0,1), main = "Linear", xlab = "Noise Level", ylab = "Power", pch = 1, col = "black", type = 'b') > points((1:30)/10, power.cors[1,], pch = 2, col = "green", type = 'b') > points((1:30)/10, power.cork[1,], pch = 3, col = "blue", type = 'b') > points((1:30)/10, power.dcor[1,], pch = 4, col = "red", type = 'b') > legend("topright",c("cor pearson","cor spearman", "cor kendal","dcor"), pch = c(1,2,3), col = c("black","green","blue","red")) > > plot((1:30)/10, power.cor[2,], ylim = c(0,1), main = "Quadratic", xlab = "Noise Level", ylab = "Power", pch = 1, col = "black", type = 'b') > points((1:30)/10, power.cors[2,], pch = 2, col = "green", type = 'b') > points((1:30)/10, power.cork[2,], pch = 3, col = "blue", type = 'b') > points((1:30)/10, power.dcor[2,], pch = 4, col = "red", type = 'b') > legend("topright",c("cor pearson","cor spearman", "cor kendal","dcor"), pch = c(1,2,3), col = c("black","green","blue","red")) > > plot((1:30)/10, power.cor[3,], ylim = c(0,1), main = "Cubic", xlab = "Noise Level", ylab = "Power", pch = 1, col = "black", type = 'b') > points((1:30)/10, power.cors[3,], pch = 2, col = "green", type = 'b') > points((1:30)/10, power.cork[3,], pch = 3, col = "blue", type = 'b') > points((1:30)/10, power.dcor[3,], pch = 4, col = "red", type = 'b') > legend("topright",c("cor pearson","cor spearman", "cor kendal","dcor"), pch = c(1,2,3), col = c("black","green","blue","red")) > > plot((1:30)/10, power.cor[5,], ylim = c(0,1), main = "Sine: period 1/8", xlab = "Noise Level", ylab = "Power", pch = 1, col = "black", type = 'b') > points((1:30)/10, power.cors[5,], pch = 2, col = "green", type = 'b') > points((1:30)/10, power.cork[5,], pch = 3, col = "blue", type = 'b') > points((1:30)/10, power.dcor[5,], pch = 4, col = "red", type = 'b') > legend("topright",c("cor pearson","cor spearman", "cor kendal","dcor"), pch = c(1,2,3), col = c("black","green","blue","red")) > > plot((1:30)/10, power.cor[4,], ylim = c(0,1), main = "Sine: period 1/2", xlab = "Noise Level", ylab = "Power", pch = 1, col = "black", type = 'b') > points((1:30)/10, power.cors[4,], pch = 2, col = "green", type = 'b') > points((1:30)/10, power.cork[4,], pch = 3, col = "blue", type = 'b') > points((1:30)/10, power.dcor[4,], pch = 4, col = "red", type = 'b') > legend("topright",c("cor pearson","cor spearman", "cor kendal","dcor"), pch = c(1,2,3), col = c("black","green","blue","red")) > > plot((1:30)/10, power.cor[6,], ylim = c(0,1), main = "X^(1/4)", xlab = "Noise Level", ylab = "Power", pch = 1, col = "black", type = 'b') > points((1:30)/10, power.cors[6,], pch = 2, col = "green", type = 'b') > points((1:30)/10, power.cork[6,], pch = 3, col = "blue", type = 'b') > points((1:30)/10, power.dcor[6,], pch = 4, col = "red", type = 'b') > legend("topright",c("cor pearson","cor spearman", "cor kendal","dcor"), pch = c(1,2,3), col = c("black","green","blue","red")) > > plot((1:30)/10, power.cor[7,], ylim = c(0,1), main = "Circle", xlab = "Noise Level", ylab = "Power", pch = 1, col = "black", type = 'b') > points((1:30)/10, power.cors[7,], pch = 2, col = "green", type = 'b') > points((1:30)/10, power.cork[7,], pch = 3, col = "blue", type = 'b') > points((1:30)/10, power.dcor[7,], pch = 4, col = "red", type = 'b') > legend("topright",c("cor pearson","cor spearman", "cor kendal","dcor"), pch = c(1,2,3), col = c("black","green","blue","red")) > > plot((1:30)/10, power.cor[8,], ylim = c(0,1), main = "Step function", xlab = "Noise Level", ylab = "Power", pch = 1, col = "black", type = 'b') > points((1:30)/10, power.cors[8,], pch = 2, col = "green", type = 'b') > points((1:30)/10, power.cork[8,], pch = 3, col = "blue", type = 'b') > points((1:30)/10, power.dcor[8,], pch = 4, col = "red", type = 'b') > legend("topright",c("cor pearson","cor spearman", "cor kendal","dcor"), pch = c(1,2,3), col = c("black","green","blue","red")) > > ################# > > ______________________________________________ > R-help at r-project.org mailing list -- To UNSUBSCRIBE and more, see > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide http://www.R-project.org/posting-guide.html > and provide commented, minimal, self-contained, reproducible code. >-------------- next part -------------- A non-text attachment was scrubbed... Name: power.pdf Type: application/pdf Size: 44649 bytes Desc: not available URL: <https://stat.ethz.ch/pipermail/r-help/attachments/20210509/a24b613f/attachment.pdf>
Dear Rui,
I thank you for your response but when I run the code with your few
modifications, I still don't get the 8 graphs but I get the following answer
:
null device
????????? 1
Here below my R code with your modifications. I don't know what I am still
missing ?
##############
set.seed(1)
library(energy)
# Here we define parameters which we use to simulate the data
# The number of null datasets we use to estimate our rejection reject #regions
for an alternative with level 0.05
nsim=50
# Number of alternative datasets we use to estimate our power
nsim2=50
# The number of different noise levels used
num.noise <- 30????????????????????
# A constant to determine the amount of noise
noise <- 3?
# Number of data points per simulation
n=100
# Vectors holding the null "correlations" (for pearson, for spearman,
for kendall and dcor respectively) for each # of the nsim null datasets at a
#given noise level
val.cor=val.cors=val.cork=val.dcor=rep(NA,nsim)
# Vectors holding the alternative "correlations" (for pearson, for
#spearman, for kendall and dcor respectively) #for each of the nsim2 alternative
datasets at a given noise level
val.cor2=val.cors2=val.cork2=val.dcor2= rep(NA,nsim2)
# Arrays holding the estimated power for each of the 4 "correlation"
types, for each data type (linear, #parabolic, etc...) with each noise level
power.cor=power.cors=power.cork=power.dcor= array(NA, c(8,num.noise))
## We loop through the noise level and functional form; each time we #estimate a
null distribution based on #the marginals of the data, and then #use that null
distribution to estimate power
## We use a uniformly distributed x, because in the original paper the #authors
used the same
for(l in 1:num.noise) {??
????? for(typ in 1:8) {
## This next loop simulates data under the null with the correct marginals (x is
uniform, and y is a function of a #uniform with gaussian noise)
??? for(ii in 1:nsim) {?????
????? x=runif(n)
#lin+noise??????????????????????????????????????????????????????
if(typ==1) {??????
y=x+ noise *(l/num.noise)* rnorm(n)?????
}
#parabolic+noise
if(typ==2) {??????
y=4*(x-.5)^2+? noise * (l/num.noise) * rnorm(n)????
}
#cubic+noise
if(typ==3) {??????
y=128*(x-1/3)^3-48*(x-1/3)^3-12*(x-1/3)+10* noise? * (l/num.noise) *rnorm(n)????
}
#sin+noise
if(typ==4) {??????
y=sin(4*pi*x) + 2*noise * (l/num.noise) *rnorm(n)????
}
#their sine + noise
if(typ==5) {??????
y=sin(16*pi*x) + noise * (l/num.noise) *rnorm(n)????
}
#x^(1/4) + noise
if(typ==6) {??????
y=x^(1/4) + noise * (l/num.noise) *rnorm(n)????
}
#circle
if(typ==7) {??????
y=(2*rbinom(n,1,0.5)-1) * (sqrt(1 - (2*x - 1)^2)) + noise/4*l/num.noise
*rnorm(n)????
}
#step function
if(typ==8) {??????
y = (x > 0.5) + noise*5*l/num.noise *rnorm(n)????
}??????
?
# We resimulate x so that we have the null scenario
x <- runif(n)
# Calculate the 4 correlations????? ?????
val.cor[ii]=(cor(x,y))
val.cors[ii]=(cor(x,y,method=c("spearman")))
val.cork[ii]=(cor(x,y,method=c("kendal")))
val.dcor[ii]=dcor(x,y)??? ?????????
}
## Next we calculate our 4 rejection cutoffs???????
cut.cor=quantile(val.cor,.95) ??
cut.cors=quantile(val.cors,.95)
cut.cork=quantile(val.cork,.95)
cut.dcor=quantile(val.dcor,.95)
## Next we simulate the data again, this time under the alternative
??? for(ii in 1:nsim2) {?????
????? x=runif(n)
#lin+noise??????????????????????????????????????????????????????
if(typ==1) {??????
y=x+ noise *(l/num.noise)* rnorm(n)????
}
#parabolic+noise
if(typ==2) {??????
y=4*(x-.5)^2+? noise * (l/num.noise) * rnorm(n)????
}
#cubic+noise
if(typ==3) {??????
y=128*(x-1/3)^3-48*(x-1/3)^3-12*(x-1/3)+10* noise? * (l/num.noise) *rnorm(n)????
}
#sin+noise
if(typ==4) {??????
y=sin(4*pi*x) + 2*noise * (l/num.noise) *rnorm(n)????
}
#their sine + noise
if(typ==5) {??????
y=sin(16*pi*x) + noise * (l/num.noise) *rnorm(n)????
}
#x^(1/4) + noise
if(typ==6) {??????
y=x^(1/4) + noise * (l/num.noise) *rnorm(n)????
}
#circle
if(typ==7) {??????
y=(2*rbinom(n,1,0.5)-1) * (sqrt(1 - (2*x - 1)^2)) + noise/4*l/num.noise
*rnorm(n)????
}
#step function
if(typ==8) {??????
y = (x > 0.5) + noise*5*l/num.noise *rnorm(n)????
}??????
## We again calculate our 4 "correlations"????????????
val.cor2[ii]=(cor(x,y))????
val.cors2[ii]=(cor(x,y,method=c("spearman")))
val.cork2[ii]=(cor(x,y,method=c("kendal")))
val.dcor2[ii]=dcor(x,y)??????????????
}
## Now we estimate the power as the number of alternative statistics #exceeding
our estimated cutoffs???????
power.cor[typ,l] <- sum(val.cor2 > cut.cor)/nsim2??
power.cors[typ,l] <- sum(val.cors2 > cut.cor)/nsim2
power.cork[typ,l] <- sum(val.cork2 > cut.cor)/nsim2
power.dcor[typ,l] <- sum(val.dcor2 > cut.dcor)/nsim2???????
}
}
save.image()?
## The rest of the code is for plotting the image
pdf(file = "power.pdf")
op <- par(mfrow = c(4,2), cex = 0.45)
plot((1:30)/10, power.cor[1,], ylim = c(0,1), main = "Linear", xlab =
"Noise Level", ylab = "Power", pch = 1, col =
"black", type = 'b')
points((1:30)/10, power.cors[1,], pch = 2, col = "green", type =
'b')
points((1:30)/10, power.cork[1,], pch = 3, col = "blue", type =
'b')
points((1:30)/10, power.dcor[1,], pch = 4, col = "red", type =
'b')
legend("topright",c("cor pearson","cor spearman",
"cor kendal","dcor"), pch = c(1,2,3), col =
c("black","green","blue","red"))
?plot((1:30)/10, power.cor[2,], ylim = c(0,1), main = "Quadratic",
xlab = "Noise Level", ylab = "Power", pch = 1, col =
"black", type = 'b')
points((1:30)/10, power.cors[2,], pch = 2, col = "green", type =
'b')
points((1:30)/10, power.cork[2,], pch = 3, col = "blue", type =
'b')
points((1:30)/10, power.dcor[2,], pch = 4, col = "red", type =
'b')
legend("topright",c("cor pearson","cor spearman",
"cor kendal","dcor"), pch = c(1,2,3), col =
c("black","green","blue","red"))
?plot((1:30)/10, power.cor[3,], ylim = c(0,1), main = "Cubic", xlab =
"Noise Level", ylab = "Power", pch = 1, col =
"black", type = 'b')
points((1:30)/10, power.cors[3,], pch = 2, col = "green", type =
'b')
points((1:30)/10, power.cork[3,], pch = 3, col = "blue", type =
'b')
points((1:30)/10, power.dcor[3,], pch = 4, col = "red", type =
'b')
legend("topright",c("cor pearson","cor spearman",
"cor kendal","dcor"), pch = c(1,2,3), col =
c("black","green","blue","red"))
plot((1:30)/10, power.cor[5,], ylim = c(0,1), main = "Sine: period
1/8", xlab = "Noise Level", ylab = "Power", pch = 1,
col = "black", type = 'b')
points((1:30)/10, power.cors[5,], pch = 2, col = "green", type =
'b')
points((1:30)/10, power.cork[5,], pch = 3, col = "blue", type =
'b')
points((1:30)/10, power.dcor[5,], pch = 4, col = "red", type =
'b')
legend("topright",c("cor pearson","cor spearman",
"cor kendal","dcor"), pch = c(1,2,3), col =
c("black","green","blue","red"))
plot((1:30)/10, power.cor[4,], ylim = c(0,1), main = "Sine: period
1/2", xlab = "Noise Level", ylab = "Power", pch = 1,
col = "black", type = 'b')
points((1:30)/10, power.cors[4,], pch = 2, col = "green", type =
'b')
points((1:30)/10, power.cork[4,], pch = 3, col = "blue", type =
'b')
points((1:30)/10, power.dcor[4,], pch = 4, col = "red", type =
'b')
legend("topright",c("cor pearson","cor spearman",
"cor kendal","dcor"), pch = c(1,2,3), col =
c("black","green","blue","red"))
plot((1:30)/10, power.cor[6,], ylim = c(0,1), main = "X^(1/4)", xlab =
"Noise Level", ylab = "Power", pch = 1, col =
"black", type = 'b')
points((1:30)/10, power.cors[6,], pch = 2, col = "green", type =
'b')
points((1:30)/10, power.cork[6,], pch = 3, col = "blue", type =
'b')
points((1:30)/10, power.dcor[6,], pch = 4, col = "red", type =
'b')
legend("topright",c("cor pearson","cor spearman",
"cor kendal","dcor"), pch = c(1,2,3), col =
c("black","green","blue","red"))
plot((1:30)/10, power.cor[7,], ylim = c(0,1), main = "Circle", xlab =
"Noise Level", ylab = "Power", pch = 1, col =
"black", type = 'b')
points((1:30)/10, power.cors[7,], pch = 2, col = "green", type =
'b')
points((1:30)/10, power.cork[7,], pch = 3, col = "blue", type =
'b')
points((1:30)/10, power.dcor[7,], pch = 4, col = "red", type =
'b')
legend("topright",c("cor pearson","cor spearman",
"cor kendal","dcor"), pch = c(1,2,3), col =
c("black","green","blue","red"))
plot((1:30)/10, power.cor[8,], ylim = c(0,1), main = "Step function",
xlab = "Noise Level", ylab = "Power", pch = 1, col =
"black", type = 'b')
points((1:30)/10, power.cors[8,], pch = 2, col = "green", type =
'b')
points((1:30)/10, power.cork[8,], pch = 3, col = "blue", type =
'b')
points((1:30)/10, power.dcor[8,], pch = 4, col = "red", type =
'b')
legend("topright",c("cor pearson","cor spearman",
"cor kendal","dcor"), pch = c(1,2,3), col =
c("black","green","blue","red"))
par(op)
dev.off()
#################
Le dimanche 9 mai 2021 ? 22:44:22 UTC+2, Rui Barradas <ruipbarradas at
sapo.pt> a ?crit :
Hello,
You are not closing the pdf device.
The only changes I have made to your code are right at the beginning of
the plotting instructions and at the end of the code.
## The rest of the code is for plotting the image
pdf(file = "power.pdf")
op <- par(mfrow = c(4,2), cex = 0.45)
[...]
par(op)
dev.off()
#################
The comments only line is your last code line.
The result is attached.
Hope this helps,
Rui Barradas
?s 19:39 de 09/05/21, varin sacha via R-help escreveu:> Dear R-experts,
>
> I am trying to get the 8 graphs like the ones in this paper :
> https://statweb.stanford.edu/~tibs/reshef/comment.pdf
> My R code does not show any error message neither warnings but I
d'on't get what I would like to get (I mean the 8 graphs), so I am
missing something. What's it ? Many thanks for your precious help.
>
> #################
> set.seed(1)
> library(energy)
>
> # Here we define parameters which we use to simulate the data
> # The number of null datasets we use to estimate our rejection reject
#regions for an alternative with level 0.05
> nsim=50
>
> # Number of alternative datasets we use to estimate our power
> nsim2=50
>
> # The number of different noise levels used
> num.noise <- 30
>
> # A constant to determine the amount of noise
> noise <- 3
>
> # Number of data points per simulation
> n=100
>
> # Vectors holding the null "correlations" (for pearson, for
spearman, for kendall and dcor respectively) for each # of the nsim null
datasets at a #given noise level
> val.cor=val.cors=val.cork=val.dcor=rep(NA,nsim)
>
> # Vectors holding the alternative "correlations" (for pearson,
for #spearman, for kendall and dcor respectively) #for each of the nsim2
alternative datasets at a given noise level
> val.cor2=val.cors2=val.cork2=val.dcor2= rep(NA,nsim2)
>?
>
> # Arrays holding the estimated power for each of the 4
"correlation" types, for each data type (linear, #parabolic, etc...)
with each noise level
> power.cor=power.cors=power.cork=power.dcor= array(NA, c(8,num.noise))
>
> ## We loop through the noise level and functional form; each time we
#estimate a null distribution based on #the marginals of the data, and then #use
that null distribution to estimate power
> ## We use a uniformly distributed x, because in the original paper the
#authors used the same
>
> for(l in 1:num.noise) {
>
>? ????? for(typ in 1:8) {
>
> ## This next loop simulates data under the null with the correct marginals
(x is uniform, and y is a function of a #uniform with gaussian noise)
>
>? ??? for(ii in 1:nsim) {
>? ????? x=runif(n)
>
> #lin+noise
> if(typ==1) {
> y=x+ noise *(l/num.noise)* rnorm(n)
> }
>
> #parabolic+noise
> if(typ==2) {
> y=4*(x-.5)^2+? noise * (l/num.noise) * rnorm(n)
> }
>
> #cubic+noise
> if(typ==3) {
> y=128*(x-1/3)^3-48*(x-1/3)^3-12*(x-1/3)+10* noise? * (l/num.noise)
*rnorm(n)
> }
>
> #sin+noise
> if(typ==4) {
> y=sin(4*pi*x) + 2*noise * (l/num.noise) *rnorm(n)
> }
>
> #their sine + noise
> if(typ==5) {
> y=sin(16*pi*x) + noise * (l/num.noise) *rnorm(n)
> }
>
> #x^(1/4) + noise
> if(typ==6) {
> y=x^(1/4) + noise * (l/num.noise) *rnorm(n)
> }
>
> #circle
> if(typ==7) {
> y=(2*rbinom(n,1,0.5)-1) * (sqrt(1 - (2*x - 1)^2)) + noise/4*l/num.noise
*rnorm(n)
> }
>
> #step function
> if(typ==8) {
> y = (x > 0.5) + noise*5*l/num.noise *rnorm(n)
> }
>
> # We resimulate x so that we have the null scenario
> x <- runif(n)
>
> # Calculate the 4 correlations
> val.cor[ii]=(cor(x,y))
> val.cors[ii]=(cor(x,y,method=c("spearman")))
> val.cork[ii]=(cor(x,y,method=c("kendal")))
> val.dcor[ii]=dcor(x,y)
> }
>
> ## Next we calculate our 4 rejection cutoffs
> cut.cor=quantile(val.cor,.95)
> cut.cors=quantile(val.cors,.95)
> cut.cork=quantile(val.cork,.95)
> cut.dcor=quantile(val.dcor,.95)
>
> ## Next we simulate the data again, this time under the alternative
>
>? ??? for(ii in 1:nsim2) {
>? ????? x=runif(n)
>
> #lin+noise
> if(typ==1) {
> y=x+ noise *(l/num.noise)* rnorm(n)
> }
>
> #parabolic+noise
> if(typ==2) {
> y=4*(x-.5)^2+? noise * (l/num.noise) * rnorm(n)
> }
>
> #cubic+noise
> if(typ==3) {
> y=128*(x-1/3)^3-48*(x-1/3)^3-12*(x-1/3)+10* noise? * (l/num.noise)
*rnorm(n)
> }
>
> #sin+noise
> if(typ==4) {
> y=sin(4*pi*x) + 2*noise * (l/num.noise) *rnorm(n)
> }
>
> #their sine + noise
> if(typ==5) {
> y=sin(16*pi*x) + noise * (l/num.noise) *rnorm(n)
> }
>
> #x^(1/4) + noise
> if(typ==6) {
> y=x^(1/4) + noise * (l/num.noise) *rnorm(n)
> }
>
> #circle
> if(typ==7) {
> y=(2*rbinom(n,1,0.5)-1) * (sqrt(1 - (2*x - 1)^2)) + noise/4*l/num.noise
*rnorm(n)
> }
>
> #step function
> if(typ==8) {
> y = (x > 0.5) + noise*5*l/num.noise *rnorm(n)
> }
>
> ## We again calculate our 4 "correlations"
> val.cor2[ii]=(cor(x,y))
> val.cors2[ii]=(cor(x,y,method=c("spearman")))
> val.cork2[ii]=(cor(x,y,method=c("kendal")))
> val.dcor2[ii]=dcor(x,y)
> }
>
> ## Now we estimate the power as the number of alternative statistics
#exceeding our estimated cutoffs
> power.cor[typ,l] <- sum(val.cor2 > cut.cor)/nsim2
> power.cors[typ,l] <- sum(val.cors2 > cut.cor)/nsim2
> power.cork[typ,l] <- sum(val.cork2 > cut.cor)/nsim2
> power.dcor[typ,l] <- sum(val.dcor2 > cut.dcor)/nsim2
> }
> }
>
> save.image()
>
> ## The rest of the code is for plotting the image
> pdf("power.pdf")
> par(mfrow = c(4,2), cex = 0.45)
> plot((1:30)/10, power.cor[1,], ylim = c(0,1), main = "Linear",
xlab = "Noise Level", ylab = "Power", pch = 1, col =
"black", type = 'b')
> points((1:30)/10, power.cors[1,], pch = 2, col = "green", type =
'b')
> points((1:30)/10, power.cork[1,], pch = 3, col = "blue", type =
'b')
> points((1:30)/10, power.dcor[1,], pch = 4, col = "red", type =
'b')
> legend("topright",c("cor pearson","cor
spearman", "cor kendal","dcor"), pch = c(1,2,3), col =
c("black","green","blue","red"))
>
> plot((1:30)/10, power.cor[2,], ylim = c(0,1), main = "Quadratic",
xlab = "Noise Level", ylab = "Power", pch = 1, col =
"black", type = 'b')
> points((1:30)/10, power.cors[2,], pch = 2, col = "green", type =
'b')
> points((1:30)/10, power.cork[2,], pch = 3, col = "blue", type =
'b')
> points((1:30)/10, power.dcor[2,], pch = 4, col = "red", type =
'b')
> legend("topright",c("cor pearson","cor
spearman", "cor kendal","dcor"), pch = c(1,2,3), col =
c("black","green","blue","red"))
>
> plot((1:30)/10, power.cor[3,], ylim = c(0,1), main = "Cubic",
xlab = "Noise Level", ylab = "Power", pch = 1, col =
"black", type = 'b')
> points((1:30)/10, power.cors[3,], pch = 2, col = "green", type =
'b')
> points((1:30)/10, power.cork[3,], pch = 3, col = "blue", type =
'b')
> points((1:30)/10, power.dcor[3,], pch = 4, col = "red", type =
'b')
> legend("topright",c("cor pearson","cor
spearman", "cor kendal","dcor"), pch = c(1,2,3), col =
c("black","green","blue","red"))
>
> plot((1:30)/10, power.cor[5,], ylim = c(0,1), main = "Sine: period
1/8", xlab = "Noise Level", ylab = "Power", pch = 1,
col = "black", type = 'b')
> points((1:30)/10, power.cors[5,], pch = 2, col = "green", type =
'b')
> points((1:30)/10, power.cork[5,], pch = 3, col = "blue", type =
'b')
> points((1:30)/10, power.dcor[5,], pch = 4, col = "red", type =
'b')
> legend("topright",c("cor pearson","cor
spearman", "cor kendal","dcor"), pch = c(1,2,3), col =
c("black","green","blue","red"))
>
> plot((1:30)/10, power.cor[4,], ylim = c(0,1), main = "Sine: period
1/2", xlab = "Noise Level", ylab = "Power", pch = 1,
col = "black", type = 'b')
> points((1:30)/10, power.cors[4,], pch = 2, col = "green", type =
'b')
> points((1:30)/10, power.cork[4,], pch = 3, col = "blue", type =
'b')
> points((1:30)/10, power.dcor[4,], pch = 4, col = "red", type =
'b')
> legend("topright",c("cor pearson","cor
spearman", "cor kendal","dcor"), pch = c(1,2,3), col =
c("black","green","blue","red"))
>
> plot((1:30)/10, power.cor[6,], ylim = c(0,1), main = "X^(1/4)",
xlab = "Noise Level", ylab = "Power", pch = 1, col =
"black", type = 'b')
> points((1:30)/10, power.cors[6,], pch = 2, col = "green", type =
'b')
> points((1:30)/10, power.cork[6,], pch = 3, col = "blue", type =
'b')
> points((1:30)/10, power.dcor[6,], pch = 4, col = "red", type =
'b')
> legend("topright",c("cor pearson","cor
spearman", "cor kendal","dcor"), pch = c(1,2,3), col =
c("black","green","blue","red"))
>
> plot((1:30)/10, power.cor[7,], ylim = c(0,1), main = "Circle",
xlab = "Noise Level", ylab = "Power", pch = 1, col =
"black", type = 'b')
> points((1:30)/10, power.cors[7,], pch = 2, col = "green", type =
'b')
> points((1:30)/10, power.cork[7,], pch = 3, col = "blue", type =
'b')
> points((1:30)/10, power.dcor[7,], pch = 4, col = "red", type =
'b')
> legend("topright",c("cor pearson","cor
spearman", "cor kendal","dcor"), pch = c(1,2,3), col =
c("black","green","blue","red"))
>
> plot((1:30)/10, power.cor[8,], ylim = c(0,1), main = "Step
function", xlab = "Noise Level", ylab = "Power", pch =
1, col = "black", type = 'b')
> points((1:30)/10, power.cors[8,], pch = 2, col = "green", type =
'b')
> points((1:30)/10, power.cork[8,], pch = 3, col = "blue", type =
'b')
> points((1:30)/10, power.dcor[8,], pch = 4, col = "red", type =
'b')
> legend("topright",c("cor pearson","cor
spearman", "cor kendal","dcor"), pch = c(1,2,3), col =
c("black","green","blue","red"))
>
> #################
>
> ______________________________________________
> R-help at r-project.org mailing list -- To UNSUBSCRIBE and more, see
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>
Rui, The created pdf.file is off-screen device. Indeed after dev.off() I should view the pdf file on my computer. But I don't find it. Where do I find the pdf.file ? Regards, Le dimanche 9 mai 2021 ? 22:44:22 UTC+2, Rui Barradas <ruipbarradas at sapo.pt> a ?crit : Hello, You are not closing the pdf device. The only changes I have made to your code are right at the beginning of the plotting instructions and at the end of the code. ## The rest of the code is for plotting the image pdf(file = "power.pdf") op <- par(mfrow = c(4,2), cex = 0.45) [...] par(op) dev.off() ################# The comments only line is your last code line. The result is attached. Hope this helps, Rui Barradas ?s 19:39 de 09/05/21, varin sacha via R-help escreveu:> Dear R-experts, > > I am trying to get the 8 graphs like the ones in this paper : > https://statweb.stanford.edu/~tibs/reshef/comment.pdf > My R code does not show any error message neither warnings but I d'on't get what I would like to get (I mean the 8 graphs), so I am missing something. What's it ? Many thanks for your precious help. > > ################# > set.seed(1) > library(energy) > > # Here we define parameters which we use to simulate the data > # The number of null datasets we use to estimate our rejection reject #regions for an alternative with level 0.05 > nsim=50 > > # Number of alternative datasets we use to estimate our power > nsim2=50 > > # The number of different noise levels used > num.noise <- 30 > > # A constant to determine the amount of noise > noise <- 3 > > # Number of data points per simulation > n=100 > > # Vectors holding the null "correlations" (for pearson, for spearman, for kendall and dcor respectively) for each # of the nsim null datasets at a #given noise level > val.cor=val.cors=val.cork=val.dcor=rep(NA,nsim) > > # Vectors holding the alternative "correlations" (for pearson, for #spearman, for kendall and dcor respectively) #for each of the nsim2 alternative datasets at a given noise level > val.cor2=val.cors2=val.cork2=val.dcor2= rep(NA,nsim2) >? > > # Arrays holding the estimated power for each of the 4 "correlation" types, for each data type (linear, #parabolic, etc...) with each noise level > power.cor=power.cors=power.cork=power.dcor= array(NA, c(8,num.noise)) > > ## We loop through the noise level and functional form; each time we #estimate a null distribution based on #the marginals of the data, and then #use that null distribution to estimate power > ## We use a uniformly distributed x, because in the original paper the #authors used the same > > for(l in 1:num.noise) { > >? ????? for(typ in 1:8) { > > ## This next loop simulates data under the null with the correct marginals (x is uniform, and y is a function of a #uniform with gaussian noise) > >? ??? for(ii in 1:nsim) { >? ????? x=runif(n) > > #lin+noise > if(typ==1) { > y=x+ noise *(l/num.noise)* rnorm(n) > } > > #parabolic+noise > if(typ==2) { > y=4*(x-.5)^2+? noise * (l/num.noise) * rnorm(n) > } > > #cubic+noise > if(typ==3) { > y=128*(x-1/3)^3-48*(x-1/3)^3-12*(x-1/3)+10* noise? * (l/num.noise) *rnorm(n) > } > > #sin+noise > if(typ==4) { > y=sin(4*pi*x) + 2*noise * (l/num.noise) *rnorm(n) > } > > #their sine + noise > if(typ==5) { > y=sin(16*pi*x) + noise * (l/num.noise) *rnorm(n) > } > > #x^(1/4) + noise > if(typ==6) { > y=x^(1/4) + noise * (l/num.noise) *rnorm(n) > } > > #circle > if(typ==7) { > y=(2*rbinom(n,1,0.5)-1) * (sqrt(1 - (2*x - 1)^2)) + noise/4*l/num.noise *rnorm(n) > } > > #step function > if(typ==8) { > y = (x > 0.5) + noise*5*l/num.noise *rnorm(n) > } > > # We resimulate x so that we have the null scenario > x <- runif(n) > > # Calculate the 4 correlations > val.cor[ii]=(cor(x,y)) > val.cors[ii]=(cor(x,y,method=c("spearman"))) > val.cork[ii]=(cor(x,y,method=c("kendal"))) > val.dcor[ii]=dcor(x,y) > } > > ## Next we calculate our 4 rejection cutoffs > cut.cor=quantile(val.cor,.95) > cut.cors=quantile(val.cors,.95) > cut.cork=quantile(val.cork,.95) > cut.dcor=quantile(val.dcor,.95) > > ## Next we simulate the data again, this time under the alternative > >? ??? for(ii in 1:nsim2) { >? ????? x=runif(n) > > #lin+noise > if(typ==1) { > y=x+ noise *(l/num.noise)* rnorm(n) > } > > #parabolic+noise > if(typ==2) { > y=4*(x-.5)^2+? noise * (l/num.noise) * rnorm(n) > } > > #cubic+noise > if(typ==3) { > y=128*(x-1/3)^3-48*(x-1/3)^3-12*(x-1/3)+10* noise? * (l/num.noise) *rnorm(n) > } > > #sin+noise > if(typ==4) { > y=sin(4*pi*x) + 2*noise * (l/num.noise) *rnorm(n) > } > > #their sine + noise > if(typ==5) { > y=sin(16*pi*x) + noise * (l/num.noise) *rnorm(n) > } > > #x^(1/4) + noise > if(typ==6) { > y=x^(1/4) + noise * (l/num.noise) *rnorm(n) > } > > #circle > if(typ==7) { > y=(2*rbinom(n,1,0.5)-1) * (sqrt(1 - (2*x - 1)^2)) + noise/4*l/num.noise *rnorm(n) > } > > #step function > if(typ==8) { > y = (x > 0.5) + noise*5*l/num.noise *rnorm(n) > } > > ## We again calculate our 4 "correlations" > val.cor2[ii]=(cor(x,y)) > val.cors2[ii]=(cor(x,y,method=c("spearman"))) > val.cork2[ii]=(cor(x,y,method=c("kendal"))) > val.dcor2[ii]=dcor(x,y) > } > > ## Now we estimate the power as the number of alternative statistics #exceeding our estimated cutoffs > power.cor[typ,l] <- sum(val.cor2 > cut.cor)/nsim2 > power.cors[typ,l] <- sum(val.cors2 > cut.cor)/nsim2 > power.cork[typ,l] <- sum(val.cork2 > cut.cor)/nsim2 > power.dcor[typ,l] <- sum(val.dcor2 > cut.dcor)/nsim2 > } > } > > save.image() > > ## The rest of the code is for plotting the image > pdf("power.pdf") > par(mfrow = c(4,2), cex = 0.45) > plot((1:30)/10, power.cor[1,], ylim = c(0,1), main = "Linear", xlab = "Noise Level", ylab = "Power", pch = 1, col = "black", type = 'b') > points((1:30)/10, power.cors[1,], pch = 2, col = "green", type = 'b') > points((1:30)/10, power.cork[1,], pch = 3, col = "blue", type = 'b') > points((1:30)/10, power.dcor[1,], pch = 4, col = "red", type = 'b') > legend("topright",c("cor pearson","cor spearman", "cor kendal","dcor"), pch = c(1,2,3), col = c("black","green","blue","red")) > > plot((1:30)/10, power.cor[2,], ylim = c(0,1), main = "Quadratic", xlab = "Noise Level", ylab = "Power", pch = 1, col = "black", type = 'b') > points((1:30)/10, power.cors[2,], pch = 2, col = "green", type = 'b') > points((1:30)/10, power.cork[2,], pch = 3, col = "blue", type = 'b') > points((1:30)/10, power.dcor[2,], pch = 4, col = "red", type = 'b') > legend("topright",c("cor pearson","cor spearman", "cor kendal","dcor"), pch = c(1,2,3), col = c("black","green","blue","red")) > > plot((1:30)/10, power.cor[3,], ylim = c(0,1), main = "Cubic", xlab = "Noise Level", ylab = "Power", pch = 1, col = "black", type = 'b') > points((1:30)/10, power.cors[3,], pch = 2, col = "green", type = 'b') > points((1:30)/10, power.cork[3,], pch = 3, col = "blue", type = 'b') > points((1:30)/10, power.dcor[3,], pch = 4, col = "red", type = 'b') > legend("topright",c("cor pearson","cor spearman", "cor kendal","dcor"), pch = c(1,2,3), col = c("black","green","blue","red")) > > plot((1:30)/10, power.cor[5,], ylim = c(0,1), main = "Sine: period 1/8", xlab = "Noise Level", ylab = "Power", pch = 1, col = "black", type = 'b') > points((1:30)/10, power.cors[5,], pch = 2, col = "green", type = 'b') > points((1:30)/10, power.cork[5,], pch = 3, col = "blue", type = 'b') > points((1:30)/10, power.dcor[5,], pch = 4, col = "red", type = 'b') > legend("topright",c("cor pearson","cor spearman", "cor kendal","dcor"), pch = c(1,2,3), col = c("black","green","blue","red")) > > plot((1:30)/10, power.cor[4,], ylim = c(0,1), main = "Sine: period 1/2", xlab = "Noise Level", ylab = "Power", pch = 1, col = "black", type = 'b') > points((1:30)/10, power.cors[4,], pch = 2, col = "green", type = 'b') > points((1:30)/10, power.cork[4,], pch = 3, col = "blue", type = 'b') > points((1:30)/10, power.dcor[4,], pch = 4, col = "red", type = 'b') > legend("topright",c("cor pearson","cor spearman", "cor kendal","dcor"), pch = c(1,2,3), col = c("black","green","blue","red")) > > plot((1:30)/10, power.cor[6,], ylim = c(0,1), main = "X^(1/4)", xlab = "Noise Level", ylab = "Power", pch = 1, col = "black", type = 'b') > points((1:30)/10, power.cors[6,], pch = 2, col = "green", type = 'b') > points((1:30)/10, power.cork[6,], pch = 3, col = "blue", type = 'b') > points((1:30)/10, power.dcor[6,], pch = 4, col = "red", type = 'b') > legend("topright",c("cor pearson","cor spearman", "cor kendal","dcor"), pch = c(1,2,3), col = c("black","green","blue","red")) > > plot((1:30)/10, power.cor[7,], ylim = c(0,1), main = "Circle", xlab = "Noise Level", ylab = "Power", pch = 1, col = "black", type = 'b') > points((1:30)/10, power.cors[7,], pch = 2, col = "green", type = 'b') > points((1:30)/10, power.cork[7,], pch = 3, col = "blue", type = 'b') > points((1:30)/10, power.dcor[7,], pch = 4, col = "red", type = 'b') > legend("topright",c("cor pearson","cor spearman", "cor kendal","dcor"), pch = c(1,2,3), col = c("black","green","blue","red")) > > plot((1:30)/10, power.cor[8,], ylim = c(0,1), main = "Step function", xlab = "Noise Level", ylab = "Power", pch = 1, col = "black", type = 'b') > points((1:30)/10, power.cors[8,], pch = 2, col = "green", type = 'b') > points((1:30)/10, power.cork[8,], pch = 3, col = "blue", type = 'b') > points((1:30)/10, power.dcor[8,], pch = 4, col = "red", type = 'b') > legend("topright",c("cor pearson","cor spearman", "cor kendal","dcor"), pch = c(1,2,3), col = c("black","green","blue","red")) > > ################# > > ______________________________________________ > R-help at r-project.org mailing list -- To UNSUBSCRIBE and more, see > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide http://www.R-project.org/posting-guide.html > and provide commented, minimal, self-contained, reproducible code. >