# Batches of bayesian updates for gaussian with unknown variance different from computation with all data

I'm working on a project where I continuously (in batches) update the pdf estimation for an event normally distributed. My variance is unknown, so I'm using the equations given in session 4.1.2 of this book:

I would expect that multiple small updates would lead to approximately the same parameters to one single big update. This is obviously true for the updates of the means, but not so clear for the standard deviation. If we simulate it:

from scipy.stats import norm
import numpy as np

np.random.seed(42)
n_0 = 0
m_0 = 0
s_0 = 0
cum_obs = np.array([])
batch_size = 10
n_batches = 100
for i in range(n_batches):
batch_obs = norm.rvs(loc=10, scale=5, size=batch_size)
cum_obs = np.concatenate([cum_obs, batch_obs])

n = batch_size
nu_0 = -1 + n_0
nu = nu_0 + n
x_bar = np.mean(batch_obs)
s = np.std(batch_obs)

var = (
(1 / nu)
* (
(s ** 2) * (n - 1)
+ (s_0 ** 2) * nu_0
+ n_0 * n * ((x_bar - m_0) ** 2) / (n_0 + n)
)
)

mult = (1 / nu) * (1 / (n + n_0))

m_0 = x_bar
n_0 += n
s0 = s

s = np.sqrt(var)
if i == n_batches - 1:
print("...\n")
if i < 5 or i == n_batches - 1:
print("Using update rules:\t\t\t\t", s)
print("Computing the std for all observed data:\t", cum_obs.std())
print("")

Using update rules:                          3.429529651472475
Computing the std for all observed data:     3.429529651472475

Using update rules:                          4.02186514878462
Computing the std for all observed data:     4.67859979446889

Using update rules:                          2.5465800067817965
Computing the std for all observed data:     4.424395957981525

Using update rules:                          2.6167707409764502
Computing the std for all observed data:     4.704117427614569

Using update rules:                          1.83449085034948
Computing the std for all observed data:     4.621424689362122

Using update rules:                          0.4917832154508146
Computing the std for all observed data:     4.893631038736771


I set the priors in a way that the first update in the two computations are equal, but as the batches go, the two get further and further apart and the presented equations get far from the expect value of 5.

I would like to know if there's something fundamentally wrong in what I'm doing or even if I should expect similar results for the two computations.

Edit: I noticed a typo in my code with s0 vs s_0 and the results make more sense now. But the stds computed with the equation still seems more unprecise

The update of the full sample and the sample in batches should yield the same results for all the parameters.

I also suggest you re-think the prior distributions. You don't want a prior standard deviation of zero, and your first round of degrees of freedom is a -1 when the prior sample size is zero.

R code is below:

set.seed(42)

m_true <- 10
s2_true <- 25

data <- rnorm(100, m_true, sqrt(s2_true))
batch_size <- 10
n_batches <- 10

update <- function(m_0, n_0, s_02, v_0, data)
{
ybar <- mean(data)
s2 <- var(data)
n <- length(data)
m_n <- (n * ybar + n_0*m_0) / (n + n_0)
n_n <- n_0 + n
v_n <- v_0 + n
s_n2 <- 1 / v_n * (s2*(n-1) + s_02*v_0 + n_0*n/n_n*(ybar - m_0)^2)
return(c(m_n, n_n, s_n2, v_n))
}

m_0 <- 3 # prior mean
n_0 <- 3 # prior sample size
s_02 <- 20 # prior variance
v_0 <- 1 # prior degrees of freedom

cat("Full Sample\n")
update(m_0, n_0, s_02, v_0, data)
#[1]  10.16257 100.00000  27.11061  99.00000

for (i in 1:n_batches)
{
temp <- update(m_0, n_0, s_02, v_0, data[((i-1)*batch_size + 1):(i*batch_size)])
m_0 <- temp[1]
n_0 <- temp[2]
s_02 <- temp[3]
v_0 <- temp[4]
print(paste("Update", i, "   ", paste(round(temp, 5), collapse=", ")))
}
#[1] "Update 1     12.73648, 10, 17.44937, 9"
#[1] "Update 2     10.9596, 20, 43.07483, 19"
#[1] "Update 3     10.34293, 30, 39.37739, 29"
#[1] "Update 4     9.80232, 40, 37.3558, 39"
#[1] "Update 5     9.82164, 50, 33.14745, 49"
#[1] "Update 6     9.8667, 60, 33.13476, 59"
#[1] "Update 7     10.27079, 70, 30.90093, 69"
#[1] "Update 8     10.10077, 80, 28.8877, 79"
#[1] "Update 9     10.22908, 90, 27.22863, 89"
#[1] "Update 10     10.16257, 100, 27.11061, 99"


EDIT: Removed an incorrect response about resetting the batch size.

• Thanks! But what do you mean about: 1-that I'm reseting the batch size and 2-that they should yield the same results. This depends on the prior, right? With the code you presented I believe std(data)^2 doesn't yield the same value as the update 10.
– jcp
Commented Apr 22, 2021 at 12:16
• 1. I'll take back the batch size comment and edit. I was wrong on that. On 2. I think that I showed the same results are achieved. The "update" function on all the data and the "updated" function on the last 10 results end at the same parameter estimates in the example. In my example std(data)^2 = 27.11061 which is var(data) in R. Does the python "std" function use the population or sample standard deviation? Commented Apr 23, 2021 at 0:00
• All clear! I was actually using your R code but I failed to consider the priors. Thanks a lot!
– jcp
Commented Apr 26, 2021 at 8:06