{ "cells": [ { "cell_type": "code", "execution_count": 1, "id": "92d02b64-d5dc-47b5-abbd-7ca1c413c6e7", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "⪅ (generic function with 1 method)" ] }, "execution_count": 1, "metadata": {}, "output_type": "execute_result" } ], "source": [ "using Distributions\n", "using Roots\n", "using Plots\n", "using StatsFuns: logistic, logit\n", "\n", "# \\approx TAB → ≈\n", "# \\lessapprox TAB → ⪅\n", "x ⪅ y = x < y || x ≈ y" ] }, { "cell_type": "code", "execution_count": 2, "id": "8d954a2c-d1b8-42d8-a152-a55d0c3ec235", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "confidence_interval" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Bernoulli分布モデルでの通常の信頼区間\n", "\n", "\"\"\"P値函数(exact版)\"\"\"\n", "function pvalue_exact(n, p, k)\n", " bin = Binomial(n, p)\n", " p0 = pdf(bin, k)\n", " s = sum(pdf(bin, j) for j in support(bin) if pdf(bin, j) ⪅ p0)\n", " min(1, s)\n", "end\n", "\n", "\"\"\"P値函数(正規分布近似)\"\"\"\n", "function pvalue_normal(n, p, k)\n", " bin = Binomial(n, p)\n", " normal = Normal(mean(bin), std(bin))\n", " min(1, 2cdf(normal, k), 2ccdf(normal, k))\n", "end\n", "\n", "\"\"\"信頼区間函数\"\"\"\n", "function confidence_interval(pvalue_func, n, k; α = 0.05)\n", " f(t) = pvalue_func(n, logistic(t), k) - α\n", " CI = logistic.(find_zeros(f, -10.0, 10.0))\n", " if length(CI) < 2\n", " return 2k ≤ n ? [0.0, first(CI)] : [first(CI), 1.0]\n", " else\n", " return [first(CI), last(CI)]\n", " end\n", "end" ] }, { "cell_type": "code", "execution_count": 3, "id": "ca2d2ac3-cb4b-4b28-8f1c-32a56f043489", "metadata": {}, "outputs": [ { "data": { "image/svg+xml": [ "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "n = 20\n", "x = 0:n\n", "θ = 0:0.002:1\n", "z_exact = pvalue_exact.(n, θ, x')\n", "z_normal = pvalue_normal.(n, θ, x')\n", "\n", "P = plot(; size=(500, 500), colorbar=false)\n", "heatmap!(x, θ, z_exact; clim=(-0.1, 1.0))\n", "plot!(; xlabel=\"x\", ylabel=\"θ\")\n", "plot!(; xtick=0:n, ytick=0:0.05:1)\n", "title!(\"pvalue_exact(n=$n, θ, x)\"; titlefontsize=12)\n", "\n", "Q = plot(; size=(500, 500), colorbar=false)\n", "heatmap!(x, θ, z_normal; clim=(-0.1, 1.0))\n", "plot!(; xlabel=\"x\", ylabel=\"θ\")\n", "plot!(; xtick=0:n, ytick=0:0.05:1)\n", "title!(\"pvalue_normal(n=$n, θ, x)\"; titlefontsize=12)\n", "\n", "plot(P, Q; size=(1000, 500))\n", "plot!(; leftmargin=3Plots.mm, bottommargin=3Plots.mm)" ] }, { "cell_type": "code", "execution_count": null, "id": "1b69e787-9a44-4243-a5ff-223296742165", "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "jupytext": { "formats": "ipynb,jl:hydrogen" }, "kernelspec": { "display_name": "Julia 1.6.4", "language": "julia", "name": "julia-1.6" }, "language_info": { "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", "version": "1.6.4" } }, "nbformat": 4, "nbformat_minor": 5 }