Cenzúrázott és csonkított Gaussian eloszlások Bayesian paraméterbecslése

A Gaussian eloszlás az egyik leggyakrabban szembejövő eloszlás a gyakorlati életben. Sajnos bizonyos esetekben a mintavételünk torzított ( jellemzően a farok tartományokban ) és ezért a paraméterbecslésünk nehéz lehet. A mai bejegyzésben két ilyen torzítást fogunk körbejárni: a cenzúrázást és a csonkítást.

---

Cenzúrázott Gaussian eloszlás

Cenzúrázott adatokról akkor beszélünk, ha a megfigyeléseink bizonyos értékhatár alatt vagy felett nem a valós értéket, hanem egy állandót rögzítenek. Klasszikus példa erre az esetre a mérlegek. Ha egy mérleg mérési tartománya 1 kg és 230 kg közötti, akkor egy 240 kg tárgyat nem fogunk tudni megmérni pontosan, csak abban lehetünk biztosak, hogy az 230 kg feletti. Egy ilyen Normális eloszlás a minimum és a maximum között szabályos, viszont a szélső értékekben újabb csúcsok láthatók. Hozzunk létre egy ilyen eloszlást a példa kedvért. Először vegyünk mintát egy Gaussian eloszlásból, ebben az esetben az eloszlás átlaga 2 lesz, és a szórása 5. Ezeket az értékeket szeretnénk majd a továbbiakban visszakapni a paraméterbecslés végén:

import numpy as np
import matplotlib.pyplot as plt 

true_mu = 3 # igazi átlag
true_sigma = 5 # igazi szórás
normal = np.random.normal(true_mu, true_sigma, 1000)
normal = np.random.normal(true_mu, true_sigma, 1000)
# ábra
plt.clf()
plt.title("Eredeti Gaussian eloszlás")
plt.axvline(np.mean(normal), label=f"Átlag: {np.round(np.mean(normal), 2)}", color="red")
plt.hist(normal)
plt.legend()
plt.show()

Következő lépésben cenzúrázzuk az adatokat. A bejegyzésben alsó cenzúrát fogok használni, ahol a cenzúrázási pont 0 lesz:

low = 0 # ami alatt minden mérés egyforma
# cenzúrázott adatok
censored = np.array([ i if i > low else low for i in normal ])
# ábra
import scipy.stats as stats
# Histogram 
plt.subplot(1,2,1)
# igazi minta
plt.hist(censored, density=True, label="Minta")
# standard normál pdf sűrűségfüggvénye 
bins = np.linspace(np.min(normal), np.max(normal), 100)
normalpdf = stats.norm.pdf(bins, true_mu, true_sigma)
plt.plot(bins, normalpdf, label=f"N({true_mu},{true_sigma})")
# jelmagyarázat és cím
plt.title(f"A minta histogramja vs. N({true_mu},{true_sigma})")
plt.legend()
 
# Q-Q ábra
plt.subplot(1,2,2)
stats.probplot(censored, dist="norm", sparams=(true_mu,true_sigma), plot=plt)
plt.show()

A fenti ábrán egyértelműen látszik a 0 értékű megfigyelések nagy száma. Szintén nyilvánvaló, hogy az átlag felfelé mozdult el, és őszintén megvallva a hisztogramról ránézésre azt se merném kijelenteni, hogy ezek a megfigyelések egy Gaussian eloszlásból származnak. A Q-Q ábra ennél nagyobb segítségünkre van. A vízszintes vonal nagyon jellegzetesen utal cenzúrázott vagy csonkított adatokra.

Most lássuk a paraméterbecslést. A következőkben a pymc3 Python csomag segítségével fogjuk megoldani a problémát. A megoldás első lépése az, hogy elkülönítjük a Normál eloszlást követő részét a megfigyeléseknek és azt, amit a minimum értékben tapasztalunk:

# elemszám a minimum értékben
n_left_censored = sum(censored <= low)
# elemszám a Gaussian részben
n_observed = len(censored) - n_left_censored
# a nem cenzúrázott megfigyelések 
uncensored = censored[(censored > low)]
assert len(uncensored) == n_observed

Majd ha ez megvan, megpróbáljuk minimalizálni a log-likelihoodot mind a két részre. A két részt pedig a közösen használt és értékek fogják összekötni:

from pymc3.distributions.dist_math import normal_lcdf  # bal oldali cenzúrázához
from pymc3.distributions.dist_math import normal_lccdf # jobb oldali cenzúrázáshoz
import pymc3 as pm
import arviz as az

high = np.max(censored)
# cenzúrázott adatok log-likelihood-ja 
def left_censored_likelihood(mu, sigma, n_left_censored, lower_bound):
    """ Likelihood of left-censored data. """
    return n_left_censored  * normal_lcdf(mu, sigma, lower_bound)

# Model
with pm.Model() as censored_model:
    # kereset paraméterek priorja
    mu = pm.Uniform("mu", low, high)
    sigma = pm.Uniform("sigma", np.sqrt(np.var(uncensored))/2, 2*np.sqrt(np.var(uncensored)))
    
    # Gaussian megfigyelések
    observed = pm.Normal(
        "observed",
        mu=mu,
        sigma=sigma,
        observed=uncensored,
    )
    
    # cenzúrázott megfigyelések
    left_censored = pm.Potential(
        "left_censored", left_censored_likelihood(mu, sigma, n_left_censored, low)
    )

# Mintavétel
with censored_model:
    trace = pm.sample(1000, tune=500, return_inferencedata=True)
    # eredmény ábrázolása
    az.plot_posterior(trace, var_names=["mu", "sigma"], ref_val=[true_mu, true_sigma], round_to=3)

Ami eredménye:

Az eredmény egyértelműen sokkal közelebb van a helyes paraméterekhez, mint amit naiv módon becsültünk volna.

Csonkított Gaussian eloszlás

Csonkított eloszlásról akkor beszélünk ha valamilyen okból levágjuk egy eloszlás szélső értékeit és eldobjuk az ide érkező megfigyeléseket. Egy naiv módon végrehajtott kiugróérték tisztítás tipikusan az az eset, ami ilyen adatsort tud produkálni. Az előbbiekhez hasonlóan kezdjük szintetikus tesztadatok előállításával:

# csonkított adatok
trunceted = np.array([ i for i in normal if i > low ])
# Histogram 
plt.subplot(1,2,1)
# minta
plt.hist(trunceted, density=True, label="Minta")
# standard normál pdf sűrűségfüggvénye 
bins = np.linspace(np.min(normal), np.max(normal), 100)
normalpdf = stats.norm.pdf(bins, true_mu, true_sigma)
plt.plot(bins, normalpdf, label=f"N({true_mu},{true_sigma})")
# jelmagyarázat és cím
plt.title(f"A minta histogramja vs. N({true_mu},{true_sigma})")
plt.legend()
 
# Q-Q ábra
plt.subplot(1,2,2)
stats.probplot(censored, dist="norm", sparams=(true_mu,true_sigma), plot=plt)
plt.show()

A Q-Q ábrán itt is megjelenik a jellegzetes vízszintes vonal, amit a cenzúrázott adatoknál is láttunk. Viszont itt nem látható a minimum megfigyelések nagy tömege, ebből gondolhatjuk, hogy csonkított adatokkal van dolgunk.

Szerencsénkre a PyMC3-nak van csonkított Gaussian eloszlása, így a modellépítés nagyon egyszerű:

with pm.Model() as truncated_model:
    mu = pm.Uniform("mu", low, high)
    sigma = pm.Uniform("sigma", np.sqrt(np.var(uncensored))/2, 2*np.sqrt(np.var(uncensored)))
    N = pm.TruncatedNormal('N', mu=mu,sd=sigma,lower=low, observed = trunceted)
    
with truncated_model:
    trace = pm.sample(1000, tune=500, return_inferencedata=True)
    az.plot_posterior(trace, var_names=["mu", "sigma"], ref_val=[true_mu, true_sigma], round_to=3)

Ami eredménye:

Ez se olyan rossz, főleg ha figyelembe vesszük, hogy a csonkítás során kidobtuk a megfigyelések több mint harmadát.

Irodalom

• Luis Mario Domenzain és mások — Censored Data Models