Bayes自由エネルギーの精密な評価

Precise estimation of Baysian free energies.ipynb

{
  "cells": [
    {
      "metadata": {},
      "cell_type": "markdown",
      "source": "# Bayes自由エネルギーの精密な評価\n\n黒木玄\n\n2017-11-23"
    },
    {
      "metadata": {
        "trusted": true
      },
      "cell_type": "code",
      "source": "import PyPlot; plt = PyPlot\nusing Distributions\nusing QuadGK\n#using Mamba\nusing KernelDensity\nfunction makefunc_pdfkde(X)\n    local ik = InterpKDE(kde(X))\n    local pdfkde(x) = pdf(ik, x)\n    return pdfkde\nend\nfunction makefunc_pdfkde(X,Y)\n    local ik = InterpKDE(kde((X,Y)))\n    local pdfkde(x, y) = pdf(ik, x, y)\n    return pdfkde\nend\n\nmacro sum(f_k, k, itr)\n    quote\n        begin\n            local s = zero(($(esc(k))->$(esc(f_k)))($(esc(itr))[1]))\n            for $(esc(k)) in $(esc(itr))\n                s += $(esc(f_k))\n            end\n            s\n        end\n    end\nend\n\nmacro prod(f_k, k, itr)\n    quote\n        begin\n            local p = one(($(esc(k))->$(esc(f_k)))($(esc(itr))[1]))\n            for $(esc(k)) in $(esc(itr))\n                p *= $(esc(f_k))\n            end\n            p\n        end\n    end\nend\n\nn = 128",
      "execution_count": 1,
      "outputs": [
        {
          "data": {
            "text/plain": "128"
          },
          "execution_count": 1,
          "metadata": {},
          "output_type": "execute_result"
        }
      ]
    },
    {
      "metadata": {},
      "cell_type": "markdown",
      "source": "## 分散を1に固定した正規分布モデル"
    },
    {
      "metadata": {},
      "cell_type": "markdown",
      "source": "分散を1に固定した正規分布モデル\n\n$$\np(x|\\mu) = \\frac{1}{\\sqrt{2\\pi}}\\exp\\left(-\\frac{1}{2}(x-\\mu)^2\\right)\n$$\n\nと事前分布\n\n$$\n\\varphi(\\mu) = \\frac{1}{\\sqrt{2\\pi}}\\exp\\left(-\\frac{1}{2}\\mu^2\\right)\n$$\n\nに関する x=(x_1,\\ldots,x_n)$ の分配函数を\n\n$$\nZ^\\beta(x) = \\int_{-\\infty}^\\infty (p(x_1|\\mu)\\cdots p(x_n|\\mu))^\\beta \\varphi(\\mu)\\,d\\mu\n$$\n\nと定義する. 函数 $F(\\beta|x)$ を\n\n$$\nF^\\beta(x) = -2\\log Z^\\beta(x)\n$$\n\nと定義する. このとき, 以下の公式が成立している.\n\n$$\nZ^\\beta(x)=\n\\frac{1}{\\sqrt{(2\\pi)^{n\\beta} (n\\beta+1)}}\n\\exp\\left(-\\frac{1}{2}\\left(\n  \\beta\\sum_{k=1}^n x_k^2 - \\frac{\\beta^2}{n\\beta+1}\\left(\\sum_{k=1}^n x_k\\right)^2\n\\right)\\right),\n$$\n\n$$\nF^\\beta(x)=\n\\beta\\sum_{k=1}^n x_k^2 - \\frac{\\beta^2}{n\\beta+1}\\left(\\sum_{k=1}^n x_k\\right)^2 + \n\\log(n\\beta+1) + n\\beta\\log(2\\pi),\n$$\n\n$$\n\\frac{\\partial F^\\beta(x)}{\\partial\\beta} =\n\\sum_{k=1}^n x_k^2 - \\left(\\frac{2\\beta}{n\\beta+1} - \\frac{n\\beta^2}{(n\\beta+1)^2}\\right)\n\\left(\\sum_{k=1}^n x_k\\right)^2 + \\frac{n}{n\\beta+1} + n\\log(2\\pi),\n$$\n\n$$\nF^\\beta(x) = \\int_0^1 \\frac{\\partial F^\\beta(x)}{\\partial\\beta}\\,d\\beta.\n$$\n\nこれらの公式を次のセルで数値的に確認する. \n\n$$\n\\mathrm{quadgk}(f, a, b) = \\text{numerical estimate of}\\ \\int_a^b f(x)\\,dx.\n$$"
    },
    {
      "metadata": {
        "trusted": true
      },
      "cell_type": "code",
      "source": "function Z_exact_normal1(beta, x)\n    local n = length(x)\n    return 1/sqrt( (2*pi)^(n*beta) * (n*beta+1) )*exp(-( beta*sum(x.^2) - beta^2/(n*beta+1)*sum(x)^2 )/2)\nend\n\nfunction F_exact_normal1(beta, x)\n    local n = length(x)\n    return beta*sum(x.^2) - beta^2/(n*beta+1)*sum(x)^2 + log(n*beta+1) + n*beta*log(2*pi)\nend\n\nfunction E2nLn_exact_normal1(beta, x)\n    local n = length(x)\n    return sum(x.^2) - ( 2*beta/(n*beta+1) - n*beta^2/(n*beta+1)^2 )*sum(x)^2 + n/(n*beta+1) + n*log(2*pi)\nend\n\nfunction Z_quadgk_normal1(beta, x)\n    local n = length(x)\n    local f(mu) = 1/sqrt((2*pi)^(n*beta+1))*exp(-( beta*sum((mu.-x).^2) + mu^2 )/2)\n    return quadgk(f, -20, 20)[1]\nend\n\nfunction F_quadgk_normal1(beta, x)\n    local f(beta) = E2nLn_exact_normal1(beta, x)\n    return quadgk(f, 0, beta)[1]\nend\n\nbeta = 0.7\nx = randn(n)\n@time @show -2*log(Z_exact_normal1(beta, x))\n@time @show F_exact_normal1(beta, x)\n@time @show -2*log(Z_quadgk_normal1(beta, x))\n@time @show F_quadgk_normal1(beta, x)\n;",
      "execution_count": 2,
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": "-2 * log(Z_exact_normal1(beta, x)) = 256.00967377255785\n  0.130419 seconds (67.36 k allocations: 3.584 MiB)\nF_exact_normal1(beta, x) = 256.00967377255785\n  0.026722 seconds (18.71 k allocations: 885.968 KiB)\n-2 * log(Z_quadgk_normal1(beta, x)) = 256.00967377255785\n  1.273119 seconds (1.43 M allocations: 77.392 MiB, 2.85% gc time)\nF_quadgk_normal1(beta, x) = 256.00967377255796\n  0.180445 seconds (123.67 k allocations: 6.575 MiB)\n"
        }
      ]
    },
    {
      "metadata": {},
      "cell_type": "markdown",
      "source": "## 正規分布モデル、分散1の正規分布モデル、その混合モデル"
    },
    {
      "metadata": {
        "trusted": true
      },
      "cell_type": "code",
      "source": "normal(w) = Normal(w[1], w[2])\nnormal1(w) = Normal(w[1], 1.0)\nmixnormal(w) = MixtureModel(Normal[Normal(w[2], 1.0), Normal(w[3], 1.0)], [1-w[1], w[1]])\n\npdf_normal(w,x)  = pdf(normal(w),x)\npdf_normal1(w,x) = pdf(normal1(w),x)\npdf_mixnormal(w,x) = (1-w[1])*pdf(normal1(w[2]),x) + w[1]*pdf(normal1(w[3]),x)\n\nfunction rand_prior_normal(M)\n    local mu = randn(M)\n    local sigma = abs.(randn(M))\n    return [mu';sigma']\nend\n\nfunction rand_prior_normal1(M)\n    local mu = randn(M)\n    return mu'\nend\n\nfunction rand_prior_mixnormal(M)\n    local a = rand(M)\n    local b = randn(M)\n    local c = randn(M)\n    return [a';b';c']\nend",
      "execution_count": 3,
      "outputs": [
        {
          "data": {
            "text/plain": "rand_prior_mixnormal (generic function with 1 method)"
          },
          "execution_count": 3,
          "metadata": {},
          "output_type": "execute_result"
        }
      ]
    },
    {
      "metadata": {},
      "cell_type": "markdown",
      "source": "### 高温展開の係数\n\n分配函数は形式的に\n\n$$\nZ^\\beta(x) = \\int\\left(\\prod_{i=1}^n p(x_i|w)\\right)^\\beta\\varphi(w)\\,dw = \\sum_{k=0}^\\infty a_k \\beta^k\n$$\n\nと展開できる. ここで\n\n$$\na_k = \\frac{1}{k!}\\int (-H(w))^k\\varphi(w)\\,dw,\n\\quad\n-H(w) = \\sum_{i=1}^n \\log p(x_i|w).\n$$\n\n次のセルでモンテカルロ法で `a_k = U[k+1]` を計算する."
    },
    {
      "metadata": {
        "trusted": true
      },
      "cell_type": "code",
      "source": "function hightempcoeff(pdf_model, rand_prior, sample; N=10, M=10^5)\n    local n = length(sample)\n    local negH(w) = @sum(pdf_model(w, sample[i]), i, 1:n)\n    local w = rand_prior(M)\n    local h = [negH(w[:,m]) for m in 1:M]\n    local U = Array{Float64}(N+1)\n    U[1] = 1.0\n    for k in 1:N\n        U[k+1] = @sum(h[m]^k, m, 1:M)/M/factorial(k)\n    end\n    return U\nend\n\nsample = randn(n);\n@time begin\n    @show hightempcoeff(pdf_normal1,   rand_prior_normal,    sample)\n    @show hightempcoeff(pdf_mixnormal, rand_prior_mixnormal, sample)\n    @show hightempcoeff(pdf_normal1,   rand_prior_normal,    sample)\nend;",
      "execution_count": 13,
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": "hightempcoeff(pdf_normal1, rand_prior_normal, sample) = [1.0, 29.7762, 471.689, 5150.31, 43066.7, 2.92254e5, 1.67007e6, 8.24555e6, 3.58455e7, 1.39219e8, 4.88675e8]\nhightempcoeff(pdf_mixnormal, rand_prior_mixnormal, sample) = [1.0, 29.8122, 463.034, 4925.53, 40057.1, 2.64388e5, 1.47058e6, 7.07456e6, 3.00004e7, 1.13786e8, 3.90459e8]\nhightempcoeff(pdf_normal1, rand_prior_normal, sample) = [1.0, 29.7987, 472.008, 5152.45, 43073.4, 2.92237e5, 1.6697e6, 8.24267e6, 3.58296e7, 1.39147e8, 4.8839e8]\n  2.163058 seconds (10.14 M allocations: 198.338 MiB, 5.95% gc time)\n"
        }
      ]
    },
    {
      "metadata": {},
      "cell_type": "markdown",
      "source": "高温展開の係数は急激に増大する."
    },
    {
      "metadata": {},
      "cell_type": "markdown",
      "source": "### 事前分布に関するモンテカルロ法による自由エネルギーの計算\n\n次のセルで自由エネルギーを事前分布が生成したサンプルを用いたモンテカルロ法で計算する."
    },
    {
      "metadata": {
        "trusted": true
      },
      "cell_type": "code",
      "source": "function Z(pdf_model, rand_prior, sample; M=10^5)\n    local n = length(sample)\n    local P(w) = @prod(pdf_model(w, sample[i]), i, 1:n)\n    local w = rand_prior(M)\n    return @sum(P(w[:,m]), m, 1:M)/M\nend\nF(pdf_model, rand_prior, sample; M=10^5) = -2*log(Z(pdf_model, rand_prior, sample, M=M))\n\nfunction Z(beta, pdf_model, rand_prior, sample; M=10^5)\n    local n = length(sample)\n    local P(w) = @prod(pdf_model(w, sample[i]), i, 1:n)\n    local w = rand_prior(M)\n    return @sum(P(w[:,m])^beta, m, 1:M)/M\nend\nF(beta, pdf_model, rand_prior, sample; M=10^5) = -2*log(Z(beta, pdf_model, rand_prior, sample, M=M))\n\nsample = randn(2^8)\nprintln(\"-\"^70)\nbeta = 1.0\n@show beta\n@show F_exact_normal1(beta, sample)\n@time @show F(pdf_normal1,   rand_prior_normal1,   sample)\n@time @show F(beta, pdf_normal1,   rand_prior_normal1,   sample)\n@time @show F(pdf_mixnormal, rand_prior_mixnormal, sample)\n@time @show F(beta, pdf_mixnormal, rand_prior_mixnormal, sample)\n@time @show F(pdf_normal,    rand_prior_normal,    sample)\n@time @show F(beta,pdf_normal,    rand_prior_normal,    sample)\nprintln(\"-\"^70)\nbeta = 0.3\n@show beta\n@show F_exact_normal1(beta, sample)\n@time @show F(beta, pdf_normal1,   rand_prior_normal1,   sample)\n@time @show F(beta, pdf_mixnormal, rand_prior_mixnormal, sample)\n@time @show F(beta,pdf_normal,    rand_prior_normal,    sample)\nprintln(\"-\"^70)\n;",
      "execution_count": 15,
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": "----------------------------------------------------------------------\nbeta = 1.0\nF_exact_normal1(beta, sample) = 732.3089268491906\nF(pdf_normal1, rand_prior_normal1, sample) = 732.302233037779\n  0.900596 seconds (411.78 k allocations: 16.660 MiB, 0.47% gc time)\nF(beta, pdf_normal1, rand_prior_normal1, sample) = 732.3096928655483\n  0.868880 seconds (412.57 k allocations: 16.700 MiB, 0.44% gc time)\nF(pdf_mixnormal, rand_prior_mixnormal, sample) = 733.2295531928136\n  1.645751 seconds (306.93 k allocations: 18.686 MiB, 0.21% gc time)\nF(beta, pdf_mixnormal, rand_prior_mixnormal, sample) = 733.2198856064724\n  1.671460 seconds (307.54 k allocations: 18.716 MiB)\nF(pdf_normal, rand_prior_normal, sample) = 738.0627926146833\n  0.985316 seconds (407.37 k allocations: 19.473 MiB, 0.32% gc time)\nF(beta, pdf_normal, rand_prior_normal, sample) = 737.9994094575254\n  0.980434 seconds (407.98 k allocations: 19.504 MiB, 0.51% gc time)\n----------------------------------------------------------------------\nbeta = 0.3\nF_exact_normal1(beta, sample) = 222.38872292239267\nF(beta, pdf_normal1, rand_prior_normal1, sample) = 222.3851342344374\n  0.812597 seconds (400.50 k allocations: 16.071 MiB, 0.47% gc time)\nF(beta, pdf_mixnormal, rand_prior_mixnormal, sample) = 222.84345450699283\n  1.640887 seconds (300.37 k allocations: 18.357 MiB, 0.21% gc time)\nF(beta, pdf_normal, rand_prior_normal, sample) = 227.01508213299778\n  0.972565 seconds (400.37 k allocations: 19.120 MiB, 0.35% gc time)\n----------------------------------------------------------------------\n"
        }
      ]
    },
    {
      "metadata": {
        "collapsed": true,
        "trusted": true
      },
      "cell_type": "code",
      "source": "",
      "execution_count": null,
      "outputs": []
    }
  ],
  "metadata": {
    "kernelspec": {
      "name": "julia-0.6.0",
      "display_name": "Julia 0.6.0",
      "language": "julia"
    },
    "language_info": {
      "file_extension": ".jl",
      "name": "julia",
      "mimetype": "application/julia",
      "version": "0.6.0"
    },
    "toc": {
      "threshold": 4,
      "number_sections": true,
      "toc_cell": false,
      "toc_window_display": false,
      "toc_section_display": "block",
      "sideBar": true,
      "navigate_menu": true,
      "moveMenuLeft": true,
      "widenNotebook": false,
      "colors": {
        "hover_highlight": "#DAA520",
        "selected_highlight": "#FFD700",
        "running_highlight": "#FF0000",
        "wrapper_background": "#FFFFFF",
        "sidebar_border": "#EEEEEE",
        "navigate_text": "#333333",
        "navigate_num": "#000000"
      },
      "nav_menu": {
        "height": "12px",
        "width": "252px"
      }
    },
    "gist": {
      "id": "",
      "data": {
        "description": "Bayes自由エネルギーの精密な評価",
        "public": true
      }
    }
  },
  "nbformat": 4,
  "nbformat_minor": 2
}