How to modify a logloss scoring rule for forecasted probabilities being a Beta distribution? I've been reading up on scoring rules and posts such as this one Why is LogLoss preferred over other proper scoring rules?
It's clear to me that if I have a set of biased coin and after examining each one I predict that heads is p_i for coin i then after flipping each one I can estimate my performance by taking the average of log(p_outcome) over many flips.
So if I my p_heads forecast is:
0.3, 0.4, 0.7, 0.3 and I observe HHTH then my score will be:
$
\dfrac{1}{4} \left[\log(0.3) + \log(0.4) + \log(0.3) + \log(0.3)\right]
$
However what if my predictions were:
Beta(3,7), Beta(4,6), Beta(7,3), Beta(3,7)
What would my score be then?
 A: Your setup is slightly different from a "standard" probabilistic prediction situation, because you have a probabilistic prediction for $p\sim\text{Beta}(\alpha,\beta)$, but you do not directly observe the quantity $p$ for which you made the prediction. Instead, you observe only the outcome of a Bernoulli experiment with parameter $p$.
So you have a case of a compound distribution, specifically, a Beta-Bernoulli one, which is a very simple case of a Beta-binomial distribution, with $n=1$ and $k\in\{0,1\}$. Fortunately, we can derive the predictive density of the coin toss in a Beta-Bernoulli directly from the predictive density of the parameter $p$. Namely, per the Wikipedia page, if $p$ denotes the probability of throwing heads,
$$ \begin{align*}P(\text{Heads}) &= {1\choose 1} \frac{B(1+\alpha,\beta)}{B(\alpha+\beta+1)} \\
&= \frac{\Gamma(\alpha+1)\Gamma(\beta)}{\Gamma(\alpha+\beta+1)}\cdot
\frac{\Gamma(\alpha+\beta)}{\Gamma(\alpha)\Gamma(\beta)} \\
&= \frac{\Gamma(\alpha+1)}{\Gamma(\alpha)}\cdot
\frac{\Gamma(\alpha+\beta)}{\Gamma(\alpha+\beta+1)} \\
&= \frac{\alpha}{\alpha+\beta},
\end{align*}$$
which of course is just the expectation of the $\text{Beta}(\alpha,\beta)$. Analogously (or simply by subtracting from $1$),
$$ P(\text{Tails})=\frac{\beta}{\alpha+\beta}. $$
To calculate the log score, you now simply take the logarithm of this predictive probability mass for the correct outcome, and average over trials.
So if you have four trials, with predictive densities for $p$ of $\text{Beta}(3,7)$, $\text{Beta}(4,6)$, $\text{Beta}(7,3)$ and $\text{Beta}(3,7)$ and observed outcomes of H, H, T, H, then your score would be
$$\frac{1}{4}\big(\log(0.3)+\log(0.4)+\log(0.3)+\log(0.3)\big) \approx -1.13.$$
We see that this is precisely the same score as for non-probabilistic predictions for $p$ of the expectation of the corresponding Betas.
(Note that by using $\log$ and not $-\log$, you have "positively oriented scores", where larger scores are better. There is also the opposite convention, where smaller scores are better, which is more in keeping with scores as losses.)

EDIT: as you write, since this log-loss only depends on the expectation of the Beta-Bernoulli compound, it does not differentiate between different Beta predictive distributions on $p$ with the same expectation, but different variances. For instance, we could have different density forecasts $p\sim\text{Beta}(3,7)$ and $p\sim\text{Beta}(30,70)$, which have the same expectation, but the second one is much more certain, and the log loss should really include this, by penalizing the more certain one more for incorrect predictions.

xx <- seq(0,1,.01)
opar <- par(mfrow=c(1,2),las=1,mai=c(.5,.5,.5,.5))
    plot(xx,dbeta(xx,3,7),type="l",xlab="",ylab="",main="Beta(3,7)")
    plot(xx,dbeta(xx,30,70),type="l",xlab="",ylab="",main="Beta(30,70)")
par(opar)

However, this is not a problem of the log scoring rule, but one of our sampling scheme. After all, the problem is that the PMF of the resulting Beta-Bernoulli compound is the same for both Betas,
$$P(\text{Heads})=\frac{\alpha}{\alpha+\beta}=\frac{10\alpha}{10\alpha+10\beta},$$
so it's unsurprising that the log loss cannot distinguish between two resulting compounds on the observable number of heads.
The solution is that we need to observe multiple coin tosses with the same $p$. Then instead of the degenerate Beta-Bernoulli for $n=1$, which is simply a Bernoulli again, we see a true Beta-binomial compound, with PMF
$$ P(X=k) = {n \choose k} \frac{B(k+\alpha,n-k+\beta)}{B(\alpha,\beta)}.$$
And this compound distribution indeed distinguishes between two different predictive distributions on $p$ that only differ in the variance. For instance, in the case of $p\sim\text{Beta}(3,7)$ and $p\sim\text{Beta}(30,70)$ above, both forecasts give us an expected probability for heads of $0.3$, but actually observing $k=3$ heads out of $n=10$ trials is more likely under $p\sim\text{Beta}(30,70)$ than under $p\sim\text{Beta}(3,7)$:
> beta(3+3,10-3+7)/beta(3,7)
[1] 0.001547988
> beta(3+30,10-3+70)/beta(30,70)
[1] 0.002119483

And this difference in PMFs then carries directly through to the log loss.
