Bayesian estimation for Pearson r correlation in R I'm trying to do Bayesian estimation for Pearson $r$ correlation using Fisher's r to z transformation. But I'm wondering if I have defined the likelihood accurately in my work below (R code)?
(Note: The posterior will be function(x) prior(x)*likelihood(x) based on Bayes' theorem.)
    m <-  0  # prior mean
    s <- .1  # prior sd
prior <- function(x) dnorm(atanh(x), m, s) * 1/(1-x^2) # `x` denotes `Pearson r`

 n <-  30           # Sample size (i.e., # of pair of datapoints)
SE <- 1/sqrt(n - 3) # Standard Error
 r <- .5            # Observed Correlation

likelihood <- function(x) dnorm(atanh(r), atanh(x), SE)

 A: Bayesian model of correlation with bivariate normal data using the Fisher transformation
There are a few issues here.  Firstly, you need to distinguish between the true correlation $\rho$ in the bivariate normal distribution (which is an unknown parameter) and the Pearson correlation $r$ (which is an observed statistic).  Since the latter is a function of the observed data, there is no "estimation" process applicable to it - once observation occurs, it is known with certainty.  This is evident in your coding, since you have defined r <- 0.5 as an observed value, but then purport to be trying to make an "estimate" of this quantity.
Presumably, what you really mean to do is to use the observed Pearson correlation $r$, together with its approximate distribution based on Fisher's transformation, to estimate the underlying correlation parameter $\rho$.  That would require you to set a prior distribution for $\rho$ and then apply the likelihood that emerges from the approximate distribution of the transformed Pearson correlation.

Simple Bayesian model: To facilitate the analysis, define the quantities:
$$\begin{equation} \begin{aligned}
z &\equiv \text{arctanh}(r), \\[6 pt]
z_\rho &\equiv \text{arctanh}(\rho).
\end{aligned} \end{equation}$$
We use a normal prior for $z_\rho$ with mean $\mu$ and precision $\lambda > 0$.  Using the approximate distribution of $z$ and treating $n>3$ as a fixed value gives us the Bayesian model:
$$\begin{equation} \begin{aligned}
p(z_\rho) &\propto \exp \left( - \frac{\lambda}{2} (z_\rho - \mu)^2 \right), \\[8 pt]
p(z | z_\rho) &\propto \exp \left( - \frac{n-3}{2} (z - z_\rho)^2  \right).
\end{aligned} \end{equation}$$
By completing the square we obtain the resulting posterior:
$$\begin{equation} \begin{aligned}
p(z_\rho | z) 
 &\propto \exp \left( - \frac{\lambda}{2} (z_\rho - \mu)^2 - \frac{n-3}{2} (z - z_\rho)^2 \right) \\[8 pt]
&\propto \exp \left( - \frac{\lambda}{2} z_\rho^2 + \lambda \mu z_\rho - \frac{n-3}{2} z_\rho^2 + (n-3) z z_\rho \right) \\[8 pt]
&= \exp \left( - \frac{\lambda + n-3}{2} z_\rho^2 + (\lambda \mu + (n-3) z) z_\rho \right) \\[8 pt]
&\propto \exp \left( - \frac{\lambda + n-3}{2} \left( z_\rho - \frac{\lambda \mu + (n-3) z}{\lambda + n-3} \right)^2 \right) \\[8 pt]
&\propto \text{N} \left( z_\rho \Bigg| \mu_{\text{post}} = \frac{\lambda \mu + (n-3) z}{\lambda + n-3}, \lambda_{\text{post}} = \lambda + n-3 \right) \\[8 pt]
\end{aligned} \end{equation}$$
The posterior distribution is normal, with updated mean and precision parameters.  The resulting posterior for the untransformed correlation parameter $\rho$ can be obtained from the posterior for $z_\rho$.  The density can be obtained analytically as:
$$\begin{equation} \begin{aligned}
p(\rho | r) &= \frac{1}{1-\rho^2} \cdot \sqrt{\frac{\lambda + n-3}{2 \pi}} \cdot \\[8 pt]
&\text{ } \exp \left( - \frac{\lambda + n-3}{2} \left( \text{arctanh} (\rho) - \frac{\lambda \mu + (n-3)}{\lambda + n-3} \cdot \text{arctanh}(r) \right)^2 \right).
\end{aligned} \end{equation}$$

Special case - diffuse prior: In the special case where $\lambda \rightarrow 0$ we approach a highly diffuse improper prior but we obtain a proper posterior with posterior mean $\mu_{\text{post}} \rightarrow z$ and posterior precision $\lambda_{\text{post}} \rightarrow n-3$.  In this case we have:
$$p(\rho | r) = \frac{1}{1-\rho^2} \cdot \sqrt{\frac{n-3}{2 \pi}} \cdot
\exp \left( - \frac{n-3}{2} \left( \text{arctanh} (\rho) -  \text{arctanh}(r) \right)^2 \right).$$

Coding this model: It is simple to simulate from the posterior distribution by generating values of $z_\rho$ and transforming these to values of $\rho$.  Posterior simulation can be done with the following code:
#Set your prior parameters
mu     <- 0;
lambda <- 0;

#Enter the observed statistics
n      <- 30;
r      <- 0.5;

#Simulate k correlation values from the posterior
lambda_post <- (lambda + (n - 3));
mu_post     <- (lambda*mu + (n - 3)*atanh(r))/lambda_post;
SIM         <- rnorm(k, mean = mu_post, sd = sqrt(1/lambda_post));
SIM_POST    <- tanh(SIM);

#The vector SIM_POST contains k values from the posterior for the correlation parameter

You can easily simulate any number of values from the posterior distribution or plot this distribution analytically.  This should give you everything you need in this kind of problem.
