-3

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 :

overestimated y_max limit of the plot

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 :

under estimated of y_max of distribution

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.

halfer
  • 19,824
  • 17
  • 99
  • 186
  • @akrun `str(y3)` gives `num [1:1000, 1:5] 1.04 1.05 1.04 1.05 1.06 ...` –  Sep 19 '21 at 19:14
  • @akrun . As you can see in my **UPDATE 1**, the y_max is too high : `max_y <- as.double(which.max(density(y3)$y)+20)` gives `429` –  Sep 19 '21 at 19:32
  • I want to do a pretty plot of the 5 distributions. For that, I just want to add a little padding to the y-maximum of all 5 distributions. –  Sep 19 '21 at 19:35
  • @akrun . ok, now, I have `max_y = 23.43461` and the plot on my **UPDATE 2** which is not correct since the y_max is under-estimated. –  Sep 19 '21 at 19:38
  • I am looking for the y-max among the 5 density, I don't know it a-priori : that's what I am looking for to do a pretty plot : after found it, just add a little padding to see clearly all the 5 distributions (or density, we understand ourselves) –  Sep 19 '21 at 19:43
  • @akrun I don't understand : among the 5 distributions, there is a max in height (max_y). If I take this limit and add a small padding, I will have a nice plot where all the other distributions are not over this max_y : is it clearer ? –  Sep 19 '21 at 20:38
  • It looks like you're using base R for the plot, I don't see any ggplot2 used here. Is that intended? – Jon Spring Sep 22 '21 at 06:37
  • @JonSpring yes, I would like to do it with ggplot2. This would be a fine initiation to start with it and secondly, ggplot2 makes beautiful plots from what I have seen. –  Sep 22 '21 at 06:41
  • 1
    Provide example data - https://stackoverflow.com/questions/5963269/how-to-make-a-great-r-reproducible-example – zx8754 Sep 22 '21 at 07:16
  • Is the final update ("hoping this will help") a solution? If so, that needs to be in the answer space. – halfer Sep 23 '21 at 20:53
  • Editor's note: I have rolled this back to an earlier version. Here is my justification: Rev 6 (Sep 23 at 5:56) was made after the accepted answer, and thus was not required in order to produce an answer. I made an enquiry (Sep 23 at 20:53) to see if this edit was an answer, and there was no reply. Rev 7 (Sep 26 at 16:12) seems to confirm my suspicion that rev 6 is a solution and not part of the question ("Thanks all, I managed to produce this following figure"). It also added yet another request, which would need to go in a new question. – halfer Sep 28 '21 at 19:11

1 Answers1

0

No code to generate your sample data was provided, so I used the code @akrun had used previously: y3 <- matrix(rnorm(5000), ncol = 5)

library(tidyverse)
as.data.frame(y3) %>%
    mutate(row = row_number()) %>%     # add row to simplify next step
    pivot_longer(-row) %>%             # reshape long
    ggplot(aes(value, color = name)) + # map x to value, color to name
    geom_density()

enter image description here

Jon Spring
  • 55,165
  • 4
  • 35
  • 53
  • Thanks a lot Jon ! I have provided in my **EDIT** the full R-script. Maybe you could adapt it with your solution ? –  Sep 22 '21 at 08:01
  • what is the signification of symbol `%>%` at the end of lines ? is it a return symbol ? –  Sep 22 '21 at 08:06
  • https://stackoverflow.com/questions/27125672/what-does-function-mean-in-r – Jon Spring Sep 22 '21 at 16:14
  • It seems the first 49 lines of your 57 lines of example code are used to manipulate data that you've provided an external link to. If your question is about how to plot the data as it exists at line 50, please just provide that data in a reproducible format that doesn't rely on external links. See @zx8754's link above for good advice. The `dput()` function is often the easiest way to do this, e.g. by including the output of `dput(y3)` or replacing y3 with a smaller subset of data. – Jon Spring Sep 22 '21 at 20:17
  • I followed your advice, I put in my post the content of input file `Array_total_WITH_Shot_Noise.txt`. Could have please a simple example of plotting with ggplot2 ? And moreover, how to manage the label names of y-axis ("density" how to change it ?) and x-axis ("value" how to change it). Thanks in advance, Best Regards –  Sep 23 '21 at 10:45
  • There are hundreds of resources online for learning ggplot2. Here's one of my favorites: https://r4ds.had.co.nz/data-visualisation.html – Jon Spring Sep 23 '21 at 15:31
  • 2
    Just for general awareness: the OP has a longstanding habit of wasting the time of readers: badly formatted posts, rendering text as images, begging/needy phrasing, not learning from edits made to their work, etc. I expect this is not the first case of asking for free work, but it is not at all surprising. I don't know what happened with the other contributor, but they appeared to have gone to the trouble of answering - and then subsequently deleted their work - my guess is they felt their time had been wasted too. – halfer Sep 23 '21 at 21:07
  • @JonSpring . Thanks, could you take a look please at my **EDIT**, I just need to change the `xlabel` and `ylabel` names. –  Sep 26 '21 at 16:36
  • @youpilat13: please ask a new question for your new problem. You can if you wish ping a helpful answer author and request further help from them, as long as you do it sparingly. – halfer Sep 28 '21 at 19:12