I am trying to add a padding of 20 for the maximum of 5 curves I am plotting.
library(ggplot2)
... some code to compute 5 (bins) distributions y3(nSample,bin)
color_curve <- c("red", "green", "grey", "black")
max_y <- as.double(which.max(density(y3))+20)
# Save PDF - Chisquare
pdf(file="Chisquare_Distribution.pdf")
for (i in 1:nRed) {
if (i == 1)
plot(density(y3[,i]), col="blue", xlim=c(min(y3),max(y3)), ylim=c(0,max_y), main="z= 0.9595, 1.087, 1.2395, 1.45, 1.688")
else
lines(density(y3[,i]), col=color_curve[i-1])
}
dev.off()
But I get the following error at execution :
Error in which.max(density(y3)) :
'list' object cannot be coerced to type 'double'
Execution halted
I would like to add this padding of 20 to the maximum of all 5 distributions but it fails with this error.
How to fix this ?
Update 1
Thanks for the suggested answer. Unfortunately, now, I get the following figure :
As you can see, the maximum in y_limit is too high for the 5 distributions, I don't know where it could come from.
Update 2
With the new command, I have the following figure :
I get an under-estimated value for searching the max among the 5 distributions.
Edit
I provide the entire code to generate the plots with the input file:
my_data <- read.delim("Array_total_WITH_Shot_Noise.txt", header = FALSE, sep = " ")
array_2D <- array(my_data)
z_ph <- c(0.9595, 1.087, 1.2395, 1.45, 1.688)
b_sp <- c(1.42904922, 1.52601862, 1.63866958, 1.78259615, 1.91956918)
b_ph <- c(sqrt(1+z_ph))
ratio_squared <- (b_sp/b_ph)^2
nRed <- 5
nRow <- NROW(my_data)
nSample_var <- 1000000
nSample_mc <- 1000
Cl<-my_data[,2:length(my_data)]#suppose cl=var(alm)
Cl_sp <- array(0, dim=c(nRow,nRed))
Cl_ph <- array(0, dim=c(nRow,nRed))
length(Cl)
for (i in 1:length(Cl)) {
#(shape/rate) convention :
Cl_sp[,i] <-(Cl[, i] * ratio_squared[i])
Cl_ph[,i] <- (Cl[, i])
}
L <- array_2D[,1]
L <- 2*(array_2D[,1])+1
# Weighted sum of Chi squared distribution
y3_1<-array(0,dim=c(nSample_var,nRed));y3_2<-array(0,dim=c(nSample_var,nRed));y3<-array(0,dim=c(nSample_var,nRed));
for (i in 1:nRed) {
for (j in 1:nRow) {
# Try to summing all the random variable
y3_1[,i] <- y3_1[,i] + Cl_sp[j,i] * rchisq(nSample_var,df=L[j])
y3_2[,i] <- y3_2[,i] + Cl_ph[j,i] * rchisq(nSample_var,df=L[j])
}
y3[,i] <- y3_1[,i]/y3_2[,i]
}
print(paste0('n=',nSample_mc*nSample_var))
for (i in 1:nRed) {
# compute the standard deviation of the ratio by Monte-Carlo
print(paste0('red=',i,',mean_fid = ', ratio_squared[i],',mean_exp = ', mean(y3[,i])))
print(paste0('numerator : var = ', var(y3_1[,i]), ', sigma = ', sd(y3_1[,i])))
print(paste0('denominator : var = ', var(y3_2[,i]), ', sigma = ', sd(y3_2[,i])))
print(paste0('var = ', var(y3[,i]), ', sigma = ', sd(y3[,i])))
}
print('#############################################################')
# par(mfrow=c(2,nRed))
color_curve <- c("red", "green", "grey", "black")
# Save PDF - Chisquare
pdf(file="Chisquare_Distribution.pdf")
for (i in 1:nRed) {
if (i == 1)
plot(density(y3[,i]), col="blue", xlim=c(min(y3),max(y3)), main="z= 0.9595, 1.087, 1.2395, 1.45, 1.688")
else
lines(density(y3[,i]), col=color_curve[i-1])
}
dev.off()
Hoping this will help you to do a version ggplot2 of the classical R plots.
See Question&Answers more detail:
os