最尤推定

// 依存関係のインストール
:dep image = "0.23"
:dep evcxr_image = "1.1"

:dep rand = "0.8.5"

// プロット用ライブラリ
:dep plotters = { version = "^0.3.5", default_features = false, features = ["evcxr", "all_series", "all_elements"] }

// インポート
use evcxr_image::ImageDisplay;
use image::{GenericImage, imageops::FilterType};

use plotters::prelude::*;
use rand::prelude::*;

use std::{fs::File, io::{self, prelude::*}, path::Path};

// 自作ライブラリ
:dep myml = { path = "../myml/" }

use myml::{utility::linspace, normal::normal};

データの読み込み

// ファイル
let height_file_path = Path::new("./data/height.txt");

// 読み込み
let height_file = File::open(&height_file_path)?;
let reader = io::BufReader::new(height_file);
let height_data = reader
    .lines()
    .filter_map(|line| line.ok())
    .filter_map(|val| val.parse::<f64>().ok())
    .collect::<Vec<_>>();

// 読み取れていることを確認
&height_data[..5]

[167.089607, 181.648633, 176.2728, 173.270164, 172.181037]

let count = height_data.len();

count

25000

ヒストグラムをプロット

evcxr_figure((640, 480), |root| {
    root.fill(&WHITE)?;

    let mut chart = ChartBuilder::on(&root)
        .caption("height", ("Roman", 20).into_font())
        .x_label_area_size(50)
        .y_label_area_size(50)
        .build_cartesian_2d(150..200, 0.0..0.1)?;

    chart.configure_mesh()
        .draw()?;

    // ヒストグラムを描画
    let hist = Histogram::vertical(&chart)
        .style(BLUE.filled())
        .margin(0)
        .data(height_data.iter().map(|&x| (x as i32, 1.0 / count as f64)));

    chart.draw_series(hist)?;

    Ok(())
})

正規分布による生成モデル

やること

• モデル化：身長データが正規分布に従うと仮定

• パラメータ推定：モデル（ここでは正規分布）のパラメータを推定する

すなわち，

\[ f(x; \mu, \sigma) = \frac{1}{\sqrt{2\pi\sigma^2}} \exp\left( -\frac{(x - \mu)^2}{2\sigma^2} \right) \]

の \(\mu,\sigma\) を推定する．

// 身長の平均
let height_mean = height_data.iter().sum::<f64>() / count as f64;

height_mean

172.70250853668063

// 身長の分散
let height_var = height_data
    .iter()
    .map(|h| (h - height_mean).powf(2.0))
    .sum::<f64>() / count as f64;

height_var

23.330517821055636

// 身長の標準偏差
let height_std = height_var.sqrt();

height_std

4.830167473396304

// 対応する正規分布を生成
let x = linspace(150.0, 190.0, 1000);
let y = normal(&x, height_mean, height_std);

let points = x.iter().copied().zip(y.iter().copied()).collect::<Vec<_>>();

プロット

（値が丸められているが，概ね一致していることがわかる）

evcxr_figure((640, 480), |root| {
    root.fill(&WHITE)?;

    let mut chart = ChartBuilder::on(&root)
        .caption("height", ("Roman", 20).into_font())
        .x_label_area_size(50)
        .y_label_area_size(50)
        .build_cartesian_2d(150..200, 0.0..0.1)?;

    chart.configure_mesh()
        .draw()?;

    // ヒストグラムを描画
    let hist = Histogram::vertical(&chart)
        .style(BLUE.filled())
        .margin(0)
        .data(height_data.iter().map(|&x| (x as i32, 1.0 / count as f64)));

    chart.draw_series(hist)?;

    // モデルを描画（同様にxを10倍する）
    chart.draw_series(
        LineSeries::new(
            points.iter().map(|&(x, y)| (x as i32, y)),
            &RED,
        )
    )?;
    
    Ok(())
})

データの生成

推定したモデルをもとに，元のデータに従う新たなデータを生成する．

// 乱数
:dep rand_distr = "0.4.3"

use rand::prelude::{thread_rng, Distribution};
use rand_distr::Normal;

// 乱数生成器
let mut rng = thread_rng();

// 分布
let dist = Normal::<f64>::new(height_mean, height_std)?;

// 試しにサンプリング
dist.sample(&mut rng)

176.0373531958946

// 10000個のサンプリング
let sample_size = 10000;
let generated_heights = dist.sample_iter(&mut rng).take(sample_size).collect::<Vec<_>>();

&generated_heights[..5]

[174.296516478635, 173.45843230673194, 181.0376202621591, 174.63869533511337, 173.8480648420886]

元データに重ねてプロットしてみる

evcxr_figure((640, 480), |root| {
    root.fill(&WHITE)?;

    let mut chart = ChartBuilder::on(&root)
        .caption("height data / generated data", ("Roman", 20).into_font())
        .x_label_area_size(50)
        .y_label_area_size(50)
        .build_cartesian_2d(150..200, 0.0..0.1)?;

    chart.configure_mesh()
        .draw()?;

    // 元データのヒストグラムを描画
    let hist = Histogram::vertical(&chart)
        .style(BLUE.filled())
        .margin(0)
        .data(height_data.iter().map(|&x| (x as i32, 1.0 / count as f64)));

    chart.draw_series(hist)?;

    // 生成されたデータのヒストグラムを描画
    let hist2 = Histogram::vertical(&chart)
        .style(RGBAColor(255, 128, 0, 0.4).filled())
        .margin(0)
        .data(generated_heights.iter().map(|&x| (x as i32, 1.0 / sample_size as f64)));

    chart.draw_series(hist2)?;

    Ok(())
})