{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Opgaver til stat4" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "from scipy.stats import norm\n", "from IPython.display import display, Markdown" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Problem 1\n" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "mu = 8.20\n", "sigma = 0.02\n", "samples = np.array([8.18, 8.17, 8.16, 8.15, 8.17, 8.21, 8.22, 8.16, 8.19, 8.18])\n", "n = len(samples)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Part A\n", "\n", "What conclusion can be made if $\\alpha = 0.10$ level of significance.\n", "\n", "We let our $H_0: \\mu = 8.20$.\n", "We can conclude the probability of a *type 1* error.\n", "\n", "We can calculate the stats \n", "$$\n", "v = \\frac {\\sqrt{n}} \\sigma | \\bar{X} - \\mu|\n", "$$\n", "then the *p-value* is the probability that the $H_0$ is falsely rejected when observing $v$.\n", "If the *p-value* is less that $\\alpha$ then $H_0$ is rejected.\n", "\n", "The *p-value* can be calculated by finding the area outside $v$.\n", "$$\n", "p = 2 \\cdot P(Z > v) = 2 \\cdot (1 - \\Phi(V))\n", "$$" ] }, { "cell_type": "code", "execution_count": 25, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "8.179000000000002\n", "5.143928459843998\n", "2.6905202554772245e-07 0.1 False\n" ] }, { "data": { "text/markdown": [ "Thus $H_0$ is rejected" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "x_hat = np.mean(samples)\n", "print(x_hat)\n", "\n", "def assign(alpha):\n", " v = np.sqrt(n) / sigma * np.abs(x_hat - mu)\n", " print(v)\n", " p_val = 2 * (1 - norm.cdf(v))\n", " print(p_val, alpha, p_val >= alpha)\n", " if p_val >= alpha:\n", " display(Markdown(\"Thus H_0 is accepted\"))\n", " else:\n", " display(Markdown(\"Thus $H_0$ is rejected\"))\n", " \n", "assign(0.1)" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "3.320391543176363\n", "0.0008989127881156023 0.05 False\n" ] }, { "data": { "text/markdown": [ "Thus $H_0$ is rejected" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Part b\n", "\n", "assign(0.05)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Problem 2" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Part A\n", "\n", "We want a test where we handle the case of *type II* error where $H_0$ should be rejected with a probability of $0.95$ or larger if $|\\mu - 8.20| \\geq 0.03$.\n" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "sigma_reject = 0.03" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Part B\n", "\n", "The $n$ required can be found with the following formula\n", "$$\n", "n \\approx \\frac {(z_{\\alpha / 2} + z_\\beta)^2 \\sigma^2} {(\\mu_s - \\mu)^2} \\,,\n", "$$\n", "where $n$ is the number of required samples to have $H_0$ rejected with probability $\\beta$ when the real mean is $\\mu_s$." ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "text/markdown": [ "Thus $n$ must be $6.0$" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "beta = 0.05\n", "alpha = 0.05\n", "\n", "def the_z_thing(p):\n", " return norm.ppf(1 - p)\n", "\n", "z_alpha2 = the_z_thing(alpha / 2)\n", "z_beta = the_z_thing(beta)\n", "\n", "n = ((z_alpha2 + z_beta)**2 * sigma**2) / ((mu + sigma_reject - mu)**2)\n", "display(Markdown(f\"Thus $n$ must be ${np.ceil(n)}$\"))" ] }, { "cell_type": "code", "execution_count": 33, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "p_value: 0.0 ts: 13.472193585307627\n", "H_0 is rejected\n" ] } ], "source": [ "## Part C\n", "n = 6\n", "\n", "# We start by calculating a new test statistic\n", "ts = np.sqrt(n) / sigma * (8.31 - mu)\n", "p_value = 2 * (1 - norm.cdf(np.abs(ts)))\n", "print(\"p_value:\", p_value, \"ts:\", ts)\n", "if p_value < alpha:\n", " print(\"H_0 is rejected\")\n", "else:\n", " print(\"H_0 is accepted\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Part D\n", "\n", "The probability of accepting $H_0$ when the true mean is $\\mu$ is:\n", "\n", "$$\n", "\\beta(\\mu) = \\Phi\\left(\\frac {\\mu_0 - \\mu} {\\sigma / \\sqrt{n}} + z_{\\alpha / 2}\\right) - \\Phi\\left(\\frac {\\mu_0 - \\mu} {\\sigma / \\sqrt{n}} - z_{\\alpha / 2}\\right)\n", "$$\n", "\n", "To find the probability of rejecting $H_0$ we subtract this from 1.\n", "\n", "$$\n", "P(reject \\: H_0) = 1 - \\beta(\\mu)\n", "$$" ] }, { "cell_type": "code", "execution_count": 35, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "beta: 1.8420324769867475e-37\n", "There is a 1.0 change that H_0 is rejected\n" ] } ], "source": [ "real_mu = 8.32\n", "blob = (mu - real_mu) / (sigma / np.sqrt(n))\n", "za2 = - norm.ppf(alpha / 2)\n", "\n", "beta = norm.cdf(blob + za2) - norm.cdf(blob - za2)\n", "print(\"beta:\", beta)\n", "\n", "prob = 1 - beta\n", "print(f\"There is a {prob} change that H_0 is rejected\")" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.8.9" } }, "nbformat": 4, "nbformat_minor": 4 }