ML Topluluk Günü 9 Kasım! TensorFlow, JAX güncellemeler için bize katılın ve daha fazla bilgi edinin

Varyasyon Çıkarımını Kullanarak Genelleştirilmiş Doğrusal Karışık Etki Modellerini Uydurma

TensorFlow.org'da görüntüleyin Google Colab'da çalıştırın Kaynağı GitHub'da görüntüleyin Not defterini indir

Düzenlemek

Düzenlemek

Soyut

Bu ortak çalışmada, TensorFlow Probability'de varyasyon çıkarımını kullanarak genelleştirilmiş bir doğrusal karışık efekt modeline nasıl uyacağımızı gösteriyoruz.

Model Ailesi

Karma etki modelleri doğrusal Genelleştirilmiş (glmM) benzer modelleri doğrusal genelleştirilmiş da tahmin edilen doğrusal bir yanıt bir örnek, belirli bir gürültü dahil dışında (GLM). Bu, kısmen, nadiren görülen özelliklerin daha sık görülen özelliklerle bilgi paylaşmasına izin verdiği için yararlıdır.

Üretken bir süreç olarak, Genelleştirilmiş Doğrusal Karışık Etkiler Modeli (GLMM) şu şekilde karakterize edilir:

$$ \begin{align} \text{for } & r = 1\ldots R: \hspace{2.45cm}\text{# for each random-effect group}\\ &\begin{aligned} \text{for } &c = 1\ldots |C_r|: \hspace{1.3cm}\text{# for each category ("level") of group $r$}\\ &\begin{aligned} \beta_{rc} &\sim \text{MultivariateNormal}(\text{loc}=0_{D_r}, \text{scale}=\Sigma_r^{1/2}) \end{aligned} \end{aligned}\\\\ \text{for } & i = 1 \ldots N: \hspace{2.45cm}\text{# for each sample}\\ &\begin{aligned} &\eta_i = \underbrace{\vphantom{\sum_{r=1}^R}x_i^\top\omega}_\text{fixed-effects} + \underbrace{\sum_{r=1}^R z_{r,i}^\top \beta_{r,C_r(i) } }_\text{random-effects} \\ &Y_i|x_i,\omega,\{z_{r,i} , \beta_r\}_{r=1}^R \sim \text{Distribution}(\text{mean}= g^{-1}(\eta_i)) \end{aligned} \end{align} $$

nerede:

$$ \begin{align} R &= \text{number of random-effect groups}\\ |C_r| &= \text{number of categories for group $r$}\\ N &= \text{number of training samples}\\ x_i,\omega &\in \mathbb{R}^{D_0}\\ D_0 &= \text{number of fixed-effects}\\ C_r(i) &= \text{category (under group $r$) of the $i$th sample}\\ z_{r,i} &\in \mathbb{R}^{D_r}\\ D_r &= \text{number of random-effects associated with group $r$}\\ \Sigma_{r} &\in \{S\in\mathbb{R}^{D_r \times D_r} : S \succ 0 \}\\ \eta_i\mapsto g^{-1}(\eta_i) &= \mu_i, \text{inverse link function}\\ \text{Distribution} &=\text{some distribution parameterizable solely by its mean} \end{align} $$

Başka bir deyişle, bu, her grubun her kategorisinin çok değişkenli bir normalden $\beta_{rc}$ örneğiyle ilişkili olduğunu söyler. $\beta_{rc}$ çekilişleri her zaman bağımsız olsa da, bunlar yalnızca $r$ grubu için girintili olarak dağıtılır: dikkat edin, her $r\in{1,\ldots,R}$ için tam olarak bir $\Sigma_r$ vardır. .

Bir örneğin grubunun özellikleriyle ($z_{r,i}$) yakın bir şekilde birleştirildiğinde, sonuç, $i$-th tahmin edilen doğrusal yanıtta örneğe özgü gürültüdür (aksi halde $x_i^\top\omega$'dır).

${\Sigma_r:r\in{1,\ldots,R}}$ değerini tahmin ettiğimizde, esasen rastgele efektli bir grubun taşıdığı ve aksi takdirde $x_i^\top içinde mevcut sinyali boğabilecek olan gürültü miktarını tahmin ediyoruz. \omega$.

$ \ Metin {Dağıtım} $ ve için çeşitli seçenekler vardır ters bağlantı fonksiyonu {- 1} $, $ g ^. Ortak seçenekler şunlardır:

  • $Y_i\sim\text{Normal}(\text{ortalama}=\eta_i, \text{ölçek}=\sigma)$,
  • $Y_i\sim\text{Binom}(\text{mean}=n_i \cdot \text{sigmoid}(\eta_i), \text{total_count}=n_i)$ ve,
  • $Y_i\sim\text{Poisson}(\text{mean}=\exp(\eta_i))$.

Daha fazla olasılık için bkz tfp.glm modülü.

Varyasyon çıkarımı

Ne yazık ki, $\beta,{\Sigma_r}_r^R$ parametrelerinin maksimum olabilirlik tahminlerini bulmak analitik olmayan bir integral gerektirir. Bu sorunu aşmak için bunun yerine

  1. Ekte $q_{\lambda}$ ile gösterilen parametreli bir dağılım ailesi ("vekil yoğunluk") tanımlayın.
  2. $q_{\lambda}$ gerçek hedef yoğunluğumuza yakın olacak şekilde $\lambda$ parametrelerini bulun.

Dağılımların aile, biz "minimize anlamına gelecektir bağımsız uygun boyutlarda Gauss ve "hedef yoğunluğuna yakın" tarafından olacaktır Kullback-Leibler sapma ". Örnek için, bkz : "İstatistikçiler için Bir İnceleme Varyasyon Çıkarım" Bölüm 2.2 iyi yazılmış türev ve motivasyon. Özellikle, KL sapmasını en aza indirmenin, olumsuz kanıt alt sınırını (ELBO) en aza indirmeye eşdeğer olduğunu göstermektedir.

Oyuncak Problemi

Gelman ve ark. (2007) "radon veri kümesi" bazen regresyon için yaklaşımlar göstermek için kullanılan bir veri kümesi olduğunu. (Örneğin, bu yakından ilişkili PyMC3 blog yazısı .) Radon veri kümesi Birleşik Devletleri boyunca alınan Radon kapalı ölçümlerini içerir. Radon doğal olarak, radyoaktif bir gaz olarak oluşan bir toksik yüksek konsantrasyonlarda.

Gösterimiz için, bodrum içeren evlerde Radon seviyelerinin daha yüksek olduğu hipotezini doğrulamakla ilgilendiğimizi varsayalım. Ayrıca Radon konsantrasyonunun toprak tipiyle, yani coğrafyayla ilgili olduğundan şüpheleniyoruz.

Bunu bir ML problemi olarak çerçevelemek için, okumanın yapıldığı zeminin lineer fonksiyonuna dayalı log-radon seviyelerini tahmin etmeye çalışacağız. Ayrıca ilçeyi rastgele bir etki olarak kullanacağız ve bunu yaparken coğrafyadan kaynaklanan farklılıkları hesaba katacağız. Başka bir deyişle, bir kullanacağız genelleştirilmiş doğrusal karma etki modeli .

%matplotlib inline
%config InlineBackend.figure_format = 'retina'

import os
from six.moves import urllib

import matplotlib.pyplot as plt; plt.style.use('ggplot')
import numpy as np
import pandas as pd
import seaborn as sns; sns.set_context('notebook')
import tensorflow_datasets as tfds

import tensorflow.compat.v2 as tf
tf.enable_v2_behavior()

import tensorflow_probability as tfp
tfd = tfp.distributions
tfb = tfp.bijectors

Ayrıca bir GPU'nun kullanılabilirliği için hızlı bir kontrol yapacağız:

if tf.test.gpu_device_name() != '/device:GPU:0':
  print("We'll just use the CPU for this run.")
else:
  print('Huzzah! Found GPU: {}'.format(tf.test.gpu_device_name()))
We'll just use the CPU for this run.

Veri Kümesini Alın:

Veri kümesini TensorFlow veri kümelerinden yüklüyoruz ve bazı hafif ön işlemler yapıyoruz.

def load_and_preprocess_radon_dataset(state='MN'):
  """Load the Radon dataset from TensorFlow Datasets and preprocess it.

  Following the examples in "Bayesian Data Analysis" (Gelman, 2007), we filter
  to Minnesota data and preprocess to obtain the following features:

  - `county`: Name of county in which the measurement was taken.
  - `floor`: Floor of house (0 for basement, 1 for first floor) on which the
    measurement was taken.

  The target variable is `log_radon`, the log of the Radon measurement in the
  house.
  """
  ds = tfds.load('radon', split='train')
  radon_data = tfds.as_dataframe(ds)
  radon_data.rename(lambda s: s[9:] if s.startswith('feat') else s, axis=1, inplace=True)
  df = radon_data[radon_data.state==state.encode()].copy()

  df['radon'] = df.activity.apply(lambda x: x if x > 0. else 0.1)
  # Make county names look nice. 
  df['county'] = df.county.apply(lambda s: s.decode()).str.strip().str.title()
  # Remap categories to start from 0 and end at max(category).
  df['county'] = df.county.astype(pd.api.types.CategoricalDtype())
  df['county_code'] = df.county.cat.codes
  # Radon levels are all positive, but log levels are unconstrained
  df['log_radon'] = df['radon'].apply(np.log)

  # Drop columns we won't use and tidy the index 
  columns_to_keep = ['log_radon', 'floor', 'county', 'county_code']
  df = df[columns_to_keep].reset_index(drop=True)

  return df

df = load_and_preprocess_radon_dataset()
df.head()

GLMM Ailesini Uzmanlaştırmak

Bu bölümde, GLMM ailesini radon seviyelerini tahmin etme görevi için özelleştiriyoruz. Bunu yapmak için önce bir GLMM'nin sabit etkili özel durumunu ele alıyoruz:

$$ \mathbb{E}[\log(\text{radon}_j)] = c + \text{floor_effect}_j $$

Bu model, $j$ gözlemindeki log radonunun (beklenti olarak) $j$th okumasının alındığı zemin ve ayrıca bir miktar sabit kesişim tarafından yönetildiğini varsayar. Sözde kodda yazabiliriz

def estimate_log_radon(floor):
    return intercept + floor_effect[floor]

Her katta ve evrensel için öğrenilen bir ağırlık var intercept dönem. 0 ve 1. kattaki radon ölçümlerine bakıldığında, bunun iyi bir başlangıç ​​olabileceği görülüyor:

fig, (ax1, ax2) = plt.subplots(ncols=2, figsize=(12, 4))
df.groupby('floor')['log_radon'].plot(kind='density', ax=ax1);
ax1.set_xlabel('Measured log(radon)')
ax1.legend(title='Floor')

df['floor'].value_counts().plot(kind='bar', ax=ax2)
ax2.set_xlabel('Floor where radon was measured')
ax2.set_ylabel('Count')
fig.suptitle("Distribution of log radon and floors in the dataset");

png

Modeli biraz daha karmaşık hale getirmek, coğrafya ile ilgili bir şeyler dahil etmek muhtemelen daha da iyidir: radon, toprakta mevcut olabilecek uranyumun bozunma zincirinin bir parçasıdır, bu nedenle coğrafya hesaba katılması gereken bir anahtar gibi görünüyor.

$$ \mathbb{E}[\log(\text{radon}_j)] = c + \text{floor_effect}_j + \text{county_effect}_j $$

Yine, sözde kodda,

def estimate_log_radon(floor, county):
    return intercept + floor_effect[floor] + county_effect[county]

ilçeye özgü ağırlık dışında öncekiyle aynı.

Yeterince büyük bir eğitim seti verildiğinde, bu makul bir modeldir. Ancak Minnesota'dan aldığımız verilere göre, az sayıda ölçüme sahip çok sayıda ilçe olduğunu görüyoruz. Örneğin, 85 ilçeden 39'unun beşten az gözlemi var.

Bu, ilçe başına gözlem sayısı arttıkça yukarıdaki modele yaklaşan bir şekilde tüm gözlemlerimiz arasında istatistiksel gücün paylaşılmasını motive eder.

fig, ax = plt.subplots(figsize=(22, 5));
county_freq = df['county'].value_counts()
county_freq.plot(kind='bar', ax=ax)
ax.set_xlabel('County')
ax.set_ylabel('Number of readings');

png

Biz bu modeli uyuyorsanız, county_effect vektör olasılıkla belki overfitting ve yoksul genelleme yol açan sadece birkaç antrenman örnekleri vardı ilçeleri için sonuçları ezberlemek sona ereceğini.

GLMM'ler, yukarıdaki iki GLM'ye mutlu bir orta sunar. uydurmayı düşünebiliriz

$$ \log(\text{radon}_j) \sim c + \text{floor_effect}_j + \mathcal{N}(\text{county_effect}_j, \text{county_scale}) $$

Bu model ilk aynıdır, ama biz bir normal dağılım olmak bizim olasılığını düzeltti ve (tek) değişken aracılığıyla tüm ilçeleri varyansı paylaşacak county_scale . Sözde kodda,

def estimate_log_radon(floor, county):
    county_mean = county_effect[county]
    random_effect = np.random.normal() * county_scale + county_mean
    return intercept + floor_effect[floor] + random_effect

Biz üzerinde ortak dağılımını anlaması olacak county_scale , county_mean ve random_effect bizim gözlenen verileri kullanarak. Küresel county_scale birçok gözlemlerle bu birkaç gözlemlerle ilçelerinde uymayacak bir isabet sağlar: bize ilçedeki istatistiksel gücünü paylaşmalarını sağlar. Ayrıca, daha fazla veri topladıkça, bu model havuzlanmış bir ölçek değişkeni olmadan modele yakınsar - bu veri seti ile bile, her iki modelle de en çok gözlemlenen ilçeler hakkında benzer sonuçlara varacağız.

Deney

Şimdi TensorFlow'da varyasyon çıkarımını kullanarak yukarıdaki GLMM'ye uymaya çalışacağız. İlk önce verileri özelliklere ve etiketlere böleceğiz.

features = df[['county_code', 'floor']].astype(int)
labels = df[['log_radon']].astype(np.float32).values.flatten()

Modeli Belirtin

def make_joint_distribution_coroutine(floor, county, n_counties, n_floors):

  def model():
    county_scale = yield tfd.HalfNormal(scale=1., name='scale_prior')
    intercept = yield tfd.Normal(loc=0., scale=1., name='intercept')
    floor_weight = yield tfd.Normal(loc=0., scale=1., name='floor_weight')
    county_prior = yield tfd.Normal(loc=tf.zeros(n_counties),
                                    scale=county_scale,
                                    name='county_prior')
    random_effect = tf.gather(county_prior, county, axis=-1)

    fixed_effect = intercept + floor_weight * floor
    linear_response = fixed_effect + random_effect
    yield tfd.Normal(loc=linear_response, scale=1., name='likelihood')
  return tfd.JointDistributionCoroutineAutoBatched(model)

joint = make_joint_distribution_coroutine(
    features.floor.values, features.county_code.values, df.county.nunique(),
    df.floor.nunique())

# Define a closure over the joint distribution 
# to condition on the observed labels.
def target_log_prob_fn(*args):
  return joint.log_prob(*args, likelihood=labels)

Vekil posterior belirtin

Şimdi, $\lambda$ parametrelerinin eğitilebilir olduğu $q_{\lambda}$ vekil ailesini bir araya getirdik. Bu durumda, ailemiz her parametre için bir tane olmak üzere bağımsız çok değişkenli normal dağılımlardır ve $\lambda = {(\mu_j, \sigma_j)}$, burada $j$ dört parametreyi indeksler.

Yöntem biz vekil aile kullandığı uyacak şekilde kullanmak tf.Variables . Biz de kullanmak tfp.util.TransformedVariable birlikte Softplus pozitif (öğretilebilir) ölçek parametreleri sınırlamak için. Ayrıca, başvuru Softplus tüm için scale_prior olumlu bir parametredir.

Optimizasyona yardımcı olması için bu eğitilebilir değişkenleri biraz titreşimle başlatıyoruz.

# Initialize locations and scales randomly with `tf.Variable`s and 
# `tfp.util.TransformedVariable`s.
_init_loc = lambda shape=(): tf.Variable(
    tf.random.uniform(shape, minval=-2., maxval=2.))
_init_scale = lambda shape=(): tfp.util.TransformedVariable(
    initial_value=tf.random.uniform(shape, minval=0.01, maxval=1.),
    bijector=tfb.Softplus())
n_counties = df.county.nunique()

surrogate_posterior = tfd.JointDistributionSequentialAutoBatched([
  tfb.Softplus()(tfd.Normal(_init_loc(), _init_scale())),           # scale_prior
  tfd.Normal(_init_loc(), _init_scale()),                           # intercept
  tfd.Normal(_init_loc(), _init_scale()),                           # floor_weight
  tfd.Normal(_init_loc([n_counties]), _init_scale([n_counties]))])  # county_prior

Not bu hücre ile değiştirilebilir olduğunu tfp.experimental.vi.build_factored_surrogate_posterior olduğu gibi:

surrogate_posterior = tfp.experimental.vi.build_factored_surrogate_posterior(
  event_shape=joint.event_shape_tensor()[:-1],
  constraining_bijectors=[tfb.Softplus(), None, None, None])

Sonuçlar

Amacımızın izlenebilir parametreli bir dağılım ailesi tanımlamak olduğunu hatırlayın ve ardından hedef dağılımımıza yakın izlenebilir bir dağılıma sahip olmamız için parametreleri seçin.

Yukarıda vekil dağılımı yerleşik ve kullanabilirsiniz tfp.vi.fit_surrogate_posterior negatif Elbo minimize vekil model parametrelerini (arasında Kullback-Liebler sapma en aza indirmek için corresonds bulmak için optimize edici ve adımların belirli bir sayıda kabul hangi vekil ve hedef dağılımı).

Dönüş değeri her adımda negatif ELBO olduğunu ve içinde dağılımları surrogate_posterior parametreleri optimize edici bulduğu güncellendi olacaktır.

optimizer = tf.optimizers.Adam(learning_rate=1e-2)

losses = tfp.vi.fit_surrogate_posterior(
    target_log_prob_fn, 
    surrogate_posterior,
    optimizer=optimizer,
    num_steps=3000, 
    seed=42,
    sample_size=2)

(scale_prior_, 
 intercept_, 
 floor_weight_, 
 county_weights_), _ = surrogate_posterior.sample_distributions()
print('        intercept (mean): ', intercept_.mean())
print('     floor_weight (mean): ', floor_weight_.mean())
print(' scale_prior (approx. mean): ', tf.reduce_mean(scale_prior_.sample(10000)))
intercept (mean):  tf.Tensor(1.4352839, shape=(), dtype=float32)
     floor_weight (mean):  tf.Tensor(-0.6701997, shape=(), dtype=float32)
 scale_prior (approx. mean):  tf.Tensor(0.28682157, shape=(), dtype=float32)
fig, ax = plt.subplots(figsize=(10, 3))
ax.plot(losses, 'k-')
ax.set(xlabel="Iteration",
       ylabel="Loss (ELBO)",
       title="Loss during training",
       ylim=0);

png

Tahmini ortalama ilçe etkilerini, bu ortalamanın belirsizliği ile birlikte çizebiliriz. Bunu, en büyüğü solda olacak şekilde gözlem sayısına göre sıraladık. Belirsizliğin çok sayıda gözlemi olan ilçeler için küçük olduğuna, ancak yalnızca bir veya iki gözlemi olan ilçeler için daha büyük olduğuna dikkat edin.

county_counts = (df.groupby(by=['county', 'county_code'], observed=True)
                   .agg('size')
                   .sort_values(ascending=False)
                   .reset_index(name='count'))

means = county_weights_.mean()
stds = county_weights_.stddev()

fig, ax = plt.subplots(figsize=(20, 5))

for idx, row in county_counts.iterrows():
  mid = means[row.county_code]
  std = stds[row.county_code]
  ax.vlines(idx, mid - std, mid + std, linewidth=3)
  ax.plot(idx, means[row.county_code], 'ko', mfc='w', mew=2, ms=7)

ax.set(
    xticks=np.arange(len(county_counts)),
    xlim=(-1, len(county_counts)),
    ylabel="County effect",
    title=r"Estimates of county effects on log radon levels. (mean $\pm$ 1 std. dev.)",
)
ax.set_xticklabels(county_counts.county, rotation=90);

png

Aslında, bunu, gözlemlerin log sayısını tahmini standart sapmaya karşı çizerek daha doğrudan görebiliriz ve ilişkinin yaklaşık olarak doğrusal olduğunu görebiliriz.

fig, ax = plt.subplots(figsize=(10, 7))
ax.plot(np.log1p(county_counts['count']), stds.numpy()[county_counts.county_code], 'o')
ax.set(
    ylabel='Posterior std. deviation',
    xlabel='County log-count',
    title='Having more observations generally\nlowers estimation uncertainty'
);

png

Karşılaştırarak lme4 R

%%shell
exit  # Trick to make this block not execute.

radon = read.csv('srrs2.dat', header = TRUE)
radon = radon[radon$state=='MN',]
radon$radon = ifelse(radon$activity==0., 0.1, radon$activity)
radon$log_radon = log(radon$radon)

# install.packages('lme4')
library(lme4)
fit <- lmer(log_radon ~ 1 + floor + (1 | county), data=radon)
fit

# Linear mixed model fit by REML ['lmerMod']
# Formula: log_radon ~ 1 + floor + (1 | county)
#    Data: radon
# REML criterion at convergence: 2171.305
# Random effects:
#  Groups   Name        Std.Dev.
#  county   (Intercept) 0.3282
#  Residual             0.7556
# Number of obs: 919, groups:  county, 85
# Fixed Effects:
# (Intercept)        floor
#       1.462       -0.693
<IPython.core.display.Javascript at 0x7f90b888e9b0>
<IPython.core.display.Javascript at 0x7f90b888e780>
<IPython.core.display.Javascript at 0x7f90b888e780>
<IPython.core.display.Javascript at 0x7f90bce1dfd0>
<IPython.core.display.Javascript at 0x7f90b888e780>
<IPython.core.display.Javascript at 0x7f90b888e780>
<IPython.core.display.Javascript at 0x7f90b888e780>
<IPython.core.display.Javascript at 0x7f90b888e780>

Aşağıdaki tablo sonuçları özetlemektedir.

print(pd.DataFrame(data=dict(intercept=[1.462, tf.reduce_mean(intercept_.mean()).numpy()],
                             floor=[-0.693, tf.reduce_mean(floor_weight_.mean()).numpy()],
                             scale=[0.3282, tf.reduce_mean(scale_prior_.sample(10000)).numpy()]),
                   index=['lme4', 'vi']))
intercept   floor     scale
lme4   1.462000 -0.6930  0.328200
vi     1.435284 -0.6702  0.287251

Bu tablo, VI sonuç ~ içinde% 10 olduğunu gösterir lme4 Var. Bu biraz şaşırtıcı çünkü:

  • lme4 dayanmaktadır Laplace yöntemi (değil VI)
  • bu işbirliğinde gerçekten yakınlaşmak için hiçbir çaba gösterilmedi,
  • hiperparametreleri ayarlamak için minimum çaba sarf edildi,
  • verileri düzenlemek veya önceden işlemek için hiçbir çaba gösterilmedi (örneğin, merkez özellikleri, vb.).

Çözüm

Bu ortak çalışmada, Genelleştirilmiş Doğrusal Karışık Etki Modellerini tanımladık ve TensorFlow Olasılığı kullanarak bunlara uyması için varyasyonel çıkarımın nasıl kullanılacağını gösterdik. Oyuncak problemi sadece birkaç yüz eğitim örneğine sahip olsa da, burada kullanılan teknikler ölçekte ihtiyaç duyulanlarla aynıdır.