Reproduce a prior density plot in R

问题:

I am still getting used to converting formulas to code in R. I am trying to work through the empirical Bayes chapter (chapter 6) in Computer Age Statistical Inference. I would like to produce plot 6.4, but in order to do so I need to find the marginal distribution of a function of the data.

To obtain the the data and plot it I did the following.

nodes <- read.table("https://web.stanford.edu/~hastie/CASI_files/DATA/nodes.txt",
           header = T)

nodes %>% 
ggplot(aes(x=x/n))+
  geom_histogram(bins = 30)+
  theme_bw()+
  labs(x = "nodes",
       n = "p=x/n")

He then tries the following model

With the following calculations

To get the plot I am trying to reproduce.

I am used to using standard R functions for inference such as lm() and glm(), but now that I need to maximize a mle to find a marginal distribution I am lost.