How can I practically "marginalize away a nuisance parameter"? I generate my data with this model:
$y=ax+b+\nu$
where $\nu$ is a random value from a random variable which follows a normal distribution with mean equal to zero and standard deviation equal to $\sigma$.
The following R script generates my data.
set.seed(42)
a_t=1 #true value for a
b_t=1 #true value for b
sigma_t=1 #true value for sigma
x=1:10
noise=rnorm(length(x),0,sigma_t)
y=a_t*x+b_t+noise
plot(x,y)
abline(a_t,b_t)


Now I would like to "recover" from $x$ and $y$ the values of $a$, $b$ and $\sigma$ using a Bayesian approach. I will follow the grid approximation1,2


*

*Define the grid. This means you decide how many points to use in estimating the posterior, and then you make a list of the parameter
values on the grid.

*Compute the value of the prior at each parameter value on the grid.

*Compute the likelihood at each parameter value.

*Compute the unstandardized posterior at each parameter value, by multiplying the prior by the likelihood. In this step, we are
approximating the posterior density by a discrete probability
distribution on the grid.

*Standardize the posterior, by dividing each value by the sum of all values. Finally, using the R command sample, take a random sample
with replacement from the discrete distribution.


# define the grid
a=seq(0,2,by=.05)
b=seq(0,2,by=.05)
sigma=seq(.1,2,by=.01)

pars <- expand.grid(a=a,b=b,sigma=sigma)

# define the priors
pars$a_prior <- dunif(pars$a,min=0,max=2)
pars$b_prior <- dunif(pars$b,min=0,max=2)
pars$sigma_prior <- dunif(pars$sigma,min=0,max=2)
pars$prior <- pars$a_prior * pars$b_prior * pars$sigma_prior

# compute the likelihood
sz <- nrow(pars)
for(i in 1:sz){
    likelihoods <- dnorm(y,pars$a[i]*x+pars$b[i],pars$sigma[i])
    pars$likelihood[i] <- prod(likelihoods)
    print(i/sz)
}

# compute the unnormalized posterior
pars$unnormalized_posterior <- pars$likelihood*pars$prior
# normalize the posterior
pars$posterior <- pars$unnormalized_posterior/sum(pars$unnormalized_posterior)

# get the sample from the posterior
sample_indices <- sample(1:nrow(pars), size=1e4, replace = TRUE, prob=pars$posterior)
pars_sample <- pars[sample_indices,c('a','b','sigma')]

# draw a 2d histogram for a and b
library(ggplot2)
p <- ggplot(pars_sample,aes(a,b))
p <- p+stat_bin2d(bins=10)
plot(p)

At the end of the procedure I have a three columns matrix (or the pars_sample dataframe in the R script) that is a sample from the posterior distribution $p(a,b,\sigma|x,y)$; now I would like to think at $\sigma$ as a nuisance parameter (is that meaningful?) and I would like to "marginalize it away", i.e. I would like to get
$p(a,b|x,y)=\int p(a,b,\sigma|x,y) d \sigma$
Hogg et al.3 wrote:

This marginalization can be performed by direct numerical integration
(think gridding up the nuisance parameters, evaluating the posterior
probability at every grid point, and summing)

but I have no idea on how to proceed in R.

1MCELREATH, Richard. Statistical rethinking: A Bayesian course with examples in R and Stan. CRC press, 2020. Page 40.
2ALBERT, Jim. Bayesian Computation with R. New York, NY : Springer New York, 2009. Page 27.
3HOGG, David W., BOVY, Jo and LANG, Dustin. Data analysis recipes: Fitting a model to data. arXiv:1008.4686
 A: It is as simple as dropping the $\sigma$ column from your pars_sample dataframe.
The three column dataframe is a sample from the joint distribution of all three parameters.  Dropping $\sigma$ leaves a sample from the joint distribution of the two remaining parameters.  The procedure is to generate a, b, and $\sigma$ while accounting for influence all three values, but then ignore $\sigma$ to compute marginal distribution of a and b.
This is different from the procedure you quote from Hogg et al.  Their method is a modification of your prior step.  Instead of sampling from the full joint posterior you would reduce the three-way joint posterior to a two-way (a and b) marginal posterior and sample from that.
The two methods are both valid.  You would choose the one that is more computationally convenient in terms of desired precision, computation time, complexity, etc.  From where you are it looks like the easiest thing is to generate samples of all three parameters, drop the $\sigma$, and then compute the a and b marginal distribution.
