Particle Filters (and Sequential Importance Resampling; SIR): any way to have them "learn" latent variances?

Question

I'm working my way through learning about Sequential Importance Resampling (SIR) particle filters (starting with a relatively simple example), but am a bit stuck in my understanding.

I'm particularly interested in estimating a full posteior distribution (over time), vs. just point estimates. However, the spread (variance) of the distributions I'm getting seems to be only a function the importance weighting likelihood function's assumed standard deviation (which is an input variable to the algorithm).

My question: is there a way that the SIR algorithm can learn the importance weighting likelihood's standard deviations from data? (Where the data's [observations'] variance potentially changes over time, as it does in my code and the plots below.) So that, I'd hope, the posterior will reflect not only changes in the observations' means, over time (which my code handles), but also changes in their variance (which it doesn't)?

My intuition is that the way to do this is to (somehow) model the observational likelihood's variance (the variable sdObs in my code, below) as, itself, a Markov chain, with its own prior and posterior. I tried doing this, but it didn't work.

In the first plot (link) below there are 100 timesteps with widely dispersed data ("observations"). Which are followed by another 100 timesteps with narrowly dispersed data. There's nothing in the algorithm (as coded) that adjusts the model's standard deviations -- either the transition likelihood's sd (sdTrans) or the observation likelihood's (sdObs). So, the Percentile Interval (PI, and shaded in my plots) just stays at pretty much a fixed variance across all timesteps, regardless of the "observed" data. Whereas I'd like for it to narrow after timestep 100 to reflect the narrower distribution of the observed data starting at timestep 100.

(Apologies that I don't yet have enough "reputation points" to post images! But links are below.)

Image (ggplot): Particle Filter Results on data with a fixed mean and changing variance, with a wide observation variance in the model

If I narrow the observation likelihood's standard deviation, sdObs, from 4 to 1, the whole PI gets correspondingly narrower (without regard to the distribution of observations). Basically, I'm "selecting" the posterior's standard deviation via my choice of sdObs, and which is barely influenced by the observed data's actual standard deviation. (It's fine to specify a prior, but I'm essentially just "picking" the posterior; not what I want!)

Image (ggplot): Particle Filter Results on data with a fixed mean and changing variance, with a narrow observation variance in the model

The filter does track changes in the distribution's means. Here's a plot with a change in the underlying mean.

Image (ggplot): Particle Filter Results on data with a changing mean and variance, with a narrow observation variance in the model

My code is based on this post.

library(ggplot2)
library(reshape2) 

############# Parameters ##########

# Generative params for the data (draws from beta distributions, 
# specified by an (alpha, beta) pair, and which is bounded in [0, 1].)
alphaSet = c(5, 50)
betaSet = c(3, 30)
nSet = c(100, 100) # number of timesteps with each (alpha, beta) pair

# Particle filtering parameters
sdTrans = .1 # SD for transition model
sdObs = 4 # SD for observation model
N = 10000 # number of particles 

# Percentile interval
piSpec = 0.8

set.seed(100)

############# Program Body ########

# Generate "observations" (draws from beta distributions)
samp = c(mapply(rbeta, nSet, shape1=alphaSet, shape2=betaSet))
# Convert to a logit scale so can model it with a normal distribution
sampLogit = logit(samp)
# Number of time steps (as specified by nSet)
T = sum(nSet)

### 1. Initialization (t = 0) ###

x <- matrix( nrow=N, ncol=T ) # Matrix of particles at each timestep
weights <- matrix( nrow=N, ncol=T )
x[, 1] <- rnorm(N, 0, sdTrans) # Draw particles for the 1st timestep

### 2. Importance Sampling Step (t = 0) ###

# Calculate weights, i.e. probability of evidence given sample from X
weights[, 1] <- dnorm(sampLogit[1], x[, 1], sdObs)
# Normalise weights 
weights[, 1] <- weights[, 1]/sum(weights[, 1])

### 3. Selection Step (t = 0) ###

# Weighted resampling with replacement. This ensures that X will converge 
# to the true distribution
x[, 1] <- sample(x[, 1], replace = TRUE, size = N, prob = weights[, 1]) 

for (t in seq(2, T)) {
  
  ### 2. Importance Sampling Step ###
  
  # Predict x_{t} from previous time step x_{t-1}
  # based on process (transition) model
  x[, t] <- rnorm(N, x[, t-1], sdTrans)
  # Calculate  and normalise weights
  weights[, t] <- dnorm(sampLogit[t], x[, t], sdObs) 
  weights[, t] <- weights[, t]/sum(weights[, t])
  
  ### 3. Selection Step ###
  
  # Weighted resampling with replacement
  x[, t] <- sample(x[, t], replace = TRUE, size = N, prob = weights[, t]) 
}

# Create a data frame of the particles
dfX = data.frame(t=rep(1:T), x=c(x))

############# Plot ################

# Convert back to [0, 1] space
# Calculate mean
particleMean = inv_logit( apply( x, 2, mean ) )
# Calculate Percentile Interval (PI)
piVec = c( (1-piSpec)/2, 1-((1-piSpec)/2) )
particlePI = inv_logit( apply( x, 2, quantile, piVec) )

# Create data frames for plotting, and plot
particleMeanDf = data.frame(time=1:T, data=samp, 
                            mean=particleMean)
particleMeanDfMelt = melt(particleMeanDf, id.vars="time")
particlePIDf = data.frame( time=1:T, 
                           low=particlePI[1,], high=particlePI[2,])
particlePIDfMelt = melt(particlePIDf, id.vars="time")
particlePlt = ggplot() +
  geom_point(data=particleMeanDfMelt, aes(x=time, y=value, color=variable),
            alpha=0.8, size=1) +
  geom_ribbon(data=particlePIDf, aes(x=time, ymin=low, ymax=high), 
              fill="steelblue", alpha=0.1) +
  geom_line(data=particlePIDf, aes(x=time, y=low), 
            fill="steelblue", alpha=0.2) +
  geom_line(data=particlePIDf, aes(x=time, y=high), 
            fill="steelblue", alpha=0.2) +
  ylim(0,1) + 
  theme_light()
plot(particlePlt)