They live in different spaces and mean different things.
A credible interval $[a,b]$ is a subset of the parameter space such that
$$
P(a\leq\Theta\leq b\mid X_1=x_1,\dots,X_n=x_n) = \alpha \, ,
$$
and it means that, after seeing the data, you believe that with probability $\alpha$ the parameter value is inside this interval.
A prediction interval $[u,v]$ is a subset of the sampling space such that
$$
P(u\leq X_{n+1}\leq v\mid X_1=x_1,\dots,X_n=x_n) = \gamma \, ,
$$
and it means that, after seeing the data, you believe that with probability $\gamma$ the value of a future observation $X_{n+1}$ will be inside this interval.
As noted by Henry, you are assuming normal distribution and it's perfectly ok if your data follows normal distribution, but will be incorrect if you cannot assume normal distribution for it. Below I describe two different approaches that you could use for unknown distribution given only datapoints x and accompanying density estimates px.
The first thing to consider is what exactly do you want to summarize using your intervals. For example, you could be interested in intervals obtained using quantiles, but you could also be interested in highest density region (see here, or here) of your distribution. While this should not make much (if any) difference in simple cases like symmetric, unimodal distributions, this will make a difference for more "complicated" distributions. Generally, quantiles will give you interval containing probability mass concentrated around the median (the middle $100\alpha\%$ of your distribution), while highest density region is a region around the modes of the distribution. This will be more clear if you compare the two plots on the picture below -- quantiles "cut" the distribution vertically, while highest density region "cuts" it horizontally.
Next thing to consider is how to deal with the fact that you have incomplete information about the distribution (assuming that we are talking about continuous distribution, you have only a bunch of points rather then a function). What you could do about it is to take the values "as is", or use some kind of interpolation, or smoothing, to obtain the "in between" values.
One approach would be to use linear interpolation (see ?approxfun in R), or alternatively something more smooth like splines (see ?splinefun in R). If you choose such approach you have to remember that interpolation algorithms have no domain knowledge about your data and can return invalid results like values below zero etc.
# grid of points
xx <- seq(min(x), max(x), by = 0.001)
# interpolate function from the sample
fx <- splinefun(x, px) # interpolating function
pxx <- pmax(0, fx(xx)) # normalize so prob >0
Second approach that you could consider is to use kernel density/mixture distribution to approximate your distribution using the data you have. The tricky part in here is to decide about optimal bandwidth.
# density of kernel density/mixture distribution
dmix <- function(x, m, s, w) {
k <- length(m)
rowSums(vapply(1:k, function(j) w[j]*dnorm(x, m[j], s[j]), numeric(length(x))))
}
# approximate function using kernel density/mixture distribution
pxx <- dmix(xx, x, rep(0.4, length.out = length(x)), px) # bandwidth 0.4 chosen arbitrary
Next, you are going to find the intervals of interest. You can either proceed numerically, or by simulation.
1a) Sampling to obtain quantile intervals
# sample from the "empirical" distribution
samp <- sample(xx, 1e5, replace = TRUE, prob = pxx)
# or sample from kernel density
idx <- sample.int(length(x), 1e5, replace = TRUE, prob = px)
samp <- rnorm(1e5, x[idx], 0.4) # this is arbitrary sd
# and take sample quantiles
quantile(samp, c(0.05, 0.975))
1b) Sampling to obtain highest density region
samp <- sample(pxx, 1e5, replace = TRUE, prob = pxx) # sample probabilities
crit <- quantile(samp, 0.05) # boundary for the lower 5% of probability mass
# values from the 95% highest density region
xx[pxx >= crit]
2a) Find quantiles numerically
cpxx <- cumsum(pxx) / sum(pxx)
xx[which(cpxx >= 0.025)[1]] # lower boundary
xx[which(cpxx >= 0.975)[1]-1] # upper boundary
2b) Find highest density region numerically
const <- sum(pxx)
spxx <- sort(pxx, decreasing = TRUE) / const
crit <- spxx[which(cumsum(spxx) >= 0.95)[1]] * const
As you can see on the plots below, in case of unimodal, symmetric distribution both methods return the same interval.
Of course, you could also try to find $100\alpha\%$ interval around some central value such that $\Pr(X \in \mu \pm \zeta) \ge \alpha$ and use some kind of optimization to find appropriate $\zeta$, but the two approaches described above seem to be used more commonly and are more intuitive.
Best Answer
Technically, you've not supplied enough information to answer as we'd need to know your prior.
However, I will assume you're looking at a $p\sim$Beta$(\alpha, \alpha)$ prior. In that case, if you have $N$ trials with a total of $Y$ successes the posterior distribution is $p | \mathcal D \sim$Beta$(\alpha + N - Y, \alpha + Y)$.
You can use the distribution function for Beta to calculate the necessary tails using those parameters.