BSTE_dice_sumscore_squared_python.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Simulation of squared sum score $S^2$ of $n$ dice"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1. Import `numpy` and `matplotlib` library"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np  # for basic math and random number generation\n",
    "import matplotlib.pyplot as plt # for plotting"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2. Set $n$ (no. of dice) and $M$ (no. of samples)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "n = 12 # no of dice\n",
    "M = int(1e6) # no of samples (\"Monte-Carlo runs\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3. Generate $M$ random samples of $n$ dice"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[5, 4, 1, ..., 3, 5, 4],\n",
       "       [2, 6, 1, ..., 1, 6, 3],\n",
       "       [5, 3, 3, ..., 2, 1, 5],\n",
       "       ...,\n",
       "       [6, 1, 3, ..., 6, 6, 1],\n",
       "       [1, 3, 2, ..., 5, 1, 5],\n",
       "       [4, 6, 2, ..., 4, 2, 3]])"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "dat = np.random.randint(6, size=(M, n)) + 1 # add 1, because random numbers are in {0, 1, ..., 5}\n",
    "display(dat)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4. Calculate sum score $S$ and its square $S^2$ for each sample"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "S = dat.sum(axis=1)\n",
    "S2 = np.multiply(S, S)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 5. Plot distribution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 1.0, 'Histogram of sum score S')"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAEICAYAAABRSj9aAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAXk0lEQVR4nO3df5BdZ33f8ffHAvPDxijEGxCSQQYUiGHqHxGyqTsUbEMk80NpZzK1E2LshioOdgqZJMShaRuS0kk6LcFuHGtUMMWF4DhOmGpcDYZCXEonNpLB2Mi2wmJMtMjG4odtjAuO4Ns/zlFzWd/dPZJ2tdaz79fMnb3nOc+553l2pc997nPPj1QVkqR2HbXYDZAkLSyDXpIaZ9BLUuMMeklqnEEvSY0z6CWpcQa9ZpVkZ5JXLXY7FlOSf5Jkd5JHkpy62O2RDpRBv4QluTfJOdPKLkzymf3LVfXSqrppjtdZnaSSPGmBmrrY/iNwaVUdW1WfX+zGPFEleWeSr/RviFNJ/myx26SOQa8nvCfAG8jzgZ2L3IbD6kB/50neDPwicE5VHQusBT65EG3TgTPoNavRUX+SdUl2JHk4ydeTvKev9un+54P9aO4VSY5K8jtJvprkgSTXJHnmyOte0K/7ZpJ/PW0/v5vk+iQfSvIwcGG/779O8mCS+5L8cZKjR16vkrw1yZeSfCfJ7yd5Yb/Nw0muG60/rY9j25rkKUkeAZYBX0jy5THbJskf9ds9lOT2JC/r192U5C0jdX/k09IhtvlFSf5Xv89vjI6ek7w0ySeSfKv/O72zL39Kkvcm2dM/3pvkKf26V/Wj8N9Kcj/wgf73clmSL/d/p+uSPGuGfyovB26sqi8DVNX9VbVlhro6zAx6HYjLgcur6jjghcB1ffkr+5/L++mNvwYu7B+vBl4AHAv8MUCSk4A/AX4BWAE8E1g5bV8bgeuB5cCHgR8AvwYcD7wCOBt467Rt1gM/DZwBvAPY0u/jBOBlwPkz9GtsW6vq+/3oFODkqnrhmG1f2/f/J/u2/jPgmzPsZ5yDbfPvAx8HfgxYBfxngCTPAP4n8DHgucCL+PuR9b/q93MKcDKwDvidkdd8DvAsuk8wm4B/Cfws8I/71/o2cOUM7bkZuCDJbyZZm2TZsO7rsKgqH0v0AdwLPAI8OPJ4FPjMtDrn9M8/DbwLOH7a66wGCnjSSNkngbeOLL8Y+DvgScC/AT4ysu7pwGMj+/ld4NNztP3twEdHlgs4c2T5VuC3Rpb/E/DeGV5rxraOvPaLZtj2LOBv6AL0qGnrbgLeMrJ84bTf7aG0+Rq6N4VV08rPBz4/wzZfBs4dWf4Z4N7++av6v8FTR9bfBZw9srxi9Pcy5vV/ge5N5rt0b3aXLfa/cR/dwxG9fraqlu9/8PhR8qhfohu53p1ke5LXz1L3ucBXR5a/Shfyz+7X7d6/oqoe5fGj4N2jC0l+MskNSe7vp3P+Pd3oftTXR57/3zHLxzLebG2dVVV9iu6TypXA15NsSXLcXNvNQ5vfAQT4bLojo/55X34CXaCPM66fzx1Z3ltV3xtZfj7w0X667EG64P8BM/xequrDVXUO3Sebi4HfS/IzM7RFh5FBr8Gq6ktVdT7wE8AfAtcnOYZuZDrdHrqg2O95wD66ILuPbroBgCRPA358+u6mLV8F3A2sqW7q6J10QTcfZmvrnKrqiqr6aeCldG+Ev9mv+i7dp5X9nnPoTf3/+7y/qv5FVT0X+GXgT5K8iO4NctwUE4zv557Rl51WfzewYXQgUFVPraqvzdG2v6uqPwdup5t+0iIz6DVYkjclmaiqH9JN80A3wtsL/JBufnu/jwC/luTEJMfSjcD/rKr20c29vyHJP+y/bHwXc4f2M4CHgUeSvAT4lfnq1xxtnVWSlyc5PcmT6YL9e3S/E4DbgH+a5Ol9CP/SfDU4yc8l2f9m+W26kP4BcAPwnCRv7798fUaS00f6+TtJJpIcTzeF9qFZdrMZeHeS5/f7nEiycYb2XJjkdf3+jkqyge6N75ZD7qwOmUGvA7Ee2NkfiXI5cF5Vfa+fenk38H/6j/lnAFcD/41uXv8rdAH4qwBVtbN/fi3d6P47wAPA92fZ928AP9/X/S/AfB6jPWNbBziub8+36aZCvkl33D3AH9HNe38d+CDdl8rz5eXALf3fYivwtqr6SlV9B3gN8AbgfuBLdF8yA/w7YAfdSPsO4HN92Uwu71/740m+Q/eF6+kz1H2Y7lPW39INAv4D8CtV9ZkZ6uswSpU3HtHi6kfRD9JNy3xlkZsjNccRvRZFkjf0UxrH0I2A76A7wkfSPDPotVg20n0RuAdYQzcN5MdLaQEMCvok65PsSjKZ5LIx65Pkin797UlO68tfnOS2kcfDSd4+z33QEaiq3tIfxfHMqjq7qnYtdpukVs05R9+f4fY3dF/wTAHbgfOr6s6ROufSfXl1Lt2XNZdX1eljXudrwOlVNXosryRpAQ25cNE6YLKq7gFIci3dx+47R+psBK7pP3rfnGR5khVVdd9InbOBLw8J+eOPP75Wr149tA+StOTdeuut36iqiXHrhgT9Sn70LMUpHn+I1bg6K+kOndvvPLrjeMdKsonu+ho873nPY8eOHQOaJkkCSDLjIHrIHP24E1mmz/fMWqc/KeaNwJ/PtJOq2lJVa6tq7cTE2DclSdJBGBL0U3TXz9hvFT962vSQOhuAz1XVoFPKJUnzZ0jQbwfW9KeHH003BbN1Wp2tdJcoTX9W5EPT5ufPZ5ZpG0nSwplzjr6q9iW5FLiR7gYMV1fVziQX9+s3A9vojriZpLvM7UX7t0/ydLojdn55/psvSZrLoNuFVdU2ujAfLds88ryAS2bY9lEef2VCSdJh4pmxktQ4g16SGmfQS1LjDHpJatygL2MldVZf9j/Glt/7B687zC2RhnNEL0mNM+glqXEGvSQ1zjl6aYyZ5uIPtL5z93oicEQvSY0z6CWpcQa9JDXOoJekxhn0ktQ4g16SGufhlVrSDvQwyvl6fQ+71OHkiF6SGmfQS1LjDHpJapxBL0mNM+glqXEGvSQ1blDQJ1mfZFeSySSXjVmfJFf0629PctrIuuVJrk9yd5K7krxiPjsgSZrdnEGfZBlwJbABOAk4P8lJ06ptANb0j03AVSPrLgc+VlUvAU4G7pqHdkuSBhoyol8HTFbVPVX1GHAtsHFanY3ANdW5GVieZEWS44BXAu8HqKrHqurB+Wu+JGkuQ4J+JbB7ZHmqLxtS5wXAXuADST6f5H1Jjhm3kySbkuxIsmPv3r2DOyBJmt2QoM+YshpY50nAacBVVXUq8F3gcXP8AFW1parWVtXaiYmJAc2SJA0x5Fo3U8AJI8urgD0D6xQwVVW39OXXM0PQS0uJ18DR4TQk6LcDa5KcCHwNOA/4+Wl1tgKXJrkWOB14qKruA0iyO8mLq2oXcDZw57y1XhpooS9eJj2RzRn0VbUvyaXAjcAy4Oqq2pnk4n79ZmAbcC4wCTwKXDTyEr8KfDjJ0cA909ZJkhbYoMsUV9U2ujAfLds88ryAS2bY9jZg7cE3UZJ0KDwzVpIaZ9BLUuMMeklqnEEvSY0z6CWpcQa9JDXOoJekxhn0ktQ4g16SGmfQS1LjDHpJapxBL0mNM+glqXEGvSQ1btBliqUjxZF+gxHvPKWF4Ihekhpn0EtS4wx6SWqcQS9JjTPoJalxBr0kNc6gl6TGDQr6JOuT7EoymeSyMeuT5Ip+/e1JThtZd2+SO5LclmTHfDZekjS3OU+YSrIMuBJ4DTAFbE+ytaruHKm2AVjTP04Hrup/7vfqqvrGvLVakjTYkBH9OmCyqu6pqseAa4GN0+psBK6pzs3A8iQr5rmtkqSDMCToVwK7R5an+rKhdQr4eJJbk2yaaSdJNiXZkWTH3r17BzRLkjTEkKDPmLI6gDpnVtVpdNM7lyR55bidVNWWqlpbVWsnJiYGNEuSNMSQoJ8CThhZXgXsGVqnqvb/fAD4KN1UkCTpMBkS9NuBNUlOTHI0cB6wdVqdrcAF/dE3ZwAPVdV9SY5J8gyAJMcArwW+OI/tlyTNYc6jbqpqX5JLgRuBZcDVVbUzycX9+s3ANuBcYBJ4FLio3/zZwEeT7N/Xn1bVx+a9F5KkGQ26Hn1VbaML89GyzSPPC7hkzHb3ACcfYhslSYfAM2MlqXEGvSQ1zlsJSkcAbzGoQ2HQ64h0pN8bVjqcnLqRpMYZ9JLUOINekhpn0EtS4wx6SWqcQS9JjTPoJalxBr0kNc6gl6TGGfSS1DiDXpIaZ9BLUuMMeklqnEEvSY0z6CWpcQa9JDXOoJekxhn0ktS4QUGfZH2SXUkmk1w2Zn2SXNGvvz3JadPWL0vy+SQ3zFfDJUnDzHnP2CTLgCuB1wBTwPYkW6vqzpFqG4A1/eN04Kr+535vA+4CjpundmsJ8L6wc5vtd+SNw7XfkBH9OmCyqu6pqseAa4GN0+psBK6pzs3A8iQrAJKsAl4HvG8e2y1JGmhI0K8Edo8sT/VlQ+u8F3gH8MPZdpJkU5IdSXbs3bt3QLMkSUMMCfqMKashdZK8Hnigqm6daydVtaWq1lbV2omJiQHNkiQNMSTop4ATRpZXAXsG1jkTeGOSe+mmfM5K8qGDbq0k6YANCfrtwJokJyY5GjgP2Dqtzlbggv7omzOAh6rqvqr67apaVVWr++0+VVVvms8OSJJmN+dRN1W1L8mlwI3AMuDqqtqZ5OJ+/WZgG3AuMAk8Cly0cE2WJB2IOYMeoKq20YX5aNnmkecFXDLHa9wE3HTALZQkHRLPjJWkxhn0ktQ4g16SGmfQS1LjDHpJapxBL0mNM+glqXEGvSQ1zqCXpMYZ9JLUOINekhpn0EtS4wZd1ExaSN4bdmHM9Hv1XrJLjyN6SWqcQS9JjTPoJalxBr0kNc6gl6TGGfSS1DiDXpIaZ9BLUuMMeklq3KCgT7I+ya4kk0kuG7M+Sa7o19+e5LS+/KlJPpvkC0l2JnnXfHdAkjS7OYM+yTLgSmADcBJwfpKTplXbAKzpH5uAq/ry7wNnVdXJwCnA+iRnzE/TJUlDDBnRrwMmq+qeqnoMuBbYOK3ORuCa6twMLE+yol9+pK/z5P5R89V4SdLchgT9SmD3yPJUXzaoTpJlSW4DHgA+UVW3jNtJkk1JdiTZsXfv3oHNlyTNZUjQZ0zZ9FH5jHWq6gdVdQqwCliX5GXjdlJVW6pqbVWtnZiYGNAsSdIQQ4J+CjhhZHkVsOdA61TVg8BNwPoDbaQk6eANCfrtwJokJyY5GjgP2Dqtzlbggv7omzOAh6rqviQTSZYDJHkacA5w9/w1X5I0lzlvPFJV+5JcCtwILAOurqqdSS7u128GtgHnApPAo8BF/eYrgA/2R+4cBVxXVTfMfzckSTMZdIepqtpGF+ajZZtHnhdwyZjtbgdOPcQ2SpIOgWfGSlLjvGesDhvvDSstDoNeWmK8afjS49SNJDXOoJekxhn0ktQ4g16SGmfQS1LjDHpJapxBL0mNM+glqXEGvSQ1zqCXpMYZ9JLUOINekhpn0EtS4wx6SWqcQS9JjTPoJalx3nhE8847SUlPLAa9JMA7T7XMqRtJatygoE+yPsmuJJNJLhuzPkmu6NffnuS0vvyEJH+V5K4kO5O8bb47IEma3ZxBn2QZcCWwATgJOD/JSdOqbQDW9I9NwFV9+T7g16vqp4AzgEvGbCtJWkBDRvTrgMmquqeqHgOuBTZOq7MRuKY6NwPLk6yoqvuq6nMAVfUd4C5g5Ty2X5I0hyFBvxLYPbI8xePDes46SVYDpwK3jNtJkk1JdiTZsXfv3gHNkiQNMSToM6asDqROkmOBvwDeXlUPj9tJVW2pqrVVtXZiYmJAsyRJQwwJ+inghJHlVcCeoXWSPJku5D9cVX958E2VJB2MIUG/HViT5MQkRwPnAVun1dkKXNAffXMG8FBV3ZckwPuBu6rqPfPacknSIHOeMFVV+5JcCtwILAOurqqdSS7u128GtgHnApPAo8BF/eZnAr8I3JHktr7snVW1bV57IUma0aAzY/tg3jatbPPI8wIuGbPdZxg/fy9JOkw8M1aSGmfQS1LjDHpJapxBL0mN8zLFkmbl5YuPfAa9Dpo3GJGODE7dSFLjDHpJapxBL0mNM+glqXEGvSQ1zqCXpMYZ9JLUOINekhpn0EtS4wx6SWqcQS9JjTPoJalxXtRMc/LiZRrHq1oeORzRS1LjDHpJapxBL0mNGxT0SdYn2ZVkMsllY9YnyRX9+tuTnDay7uokDyT54nw2XJI0zJxBn2QZcCWwATgJOD/JSdOqbQDW9I9NwFUj6/4rsH4+GitJOnBDRvTrgMmquqeqHgOuBTZOq7MRuKY6NwPLk6wAqKpPA9+az0ZLkoYbEvQrgd0jy1N92YHWmVWSTUl2JNmxd+/eA9lUkjSLIUGfMWV1EHVmVVVbqmptVa2dmJg4kE0lSbMYEvRTwAkjy6uAPQdRR5K0CIYE/XZgTZITkxwNnAdsnVZnK3BBf/TNGcBDVXXfPLdVknQQ5rwEQlXtS3IpcCOwDLi6qnYmubhfvxnYBpwLTAKPAhft3z7JR4BXAccnmQL+bVW9f747okPnpQ40H7w0whPPoGvdVNU2ujAfLds88ryAS2bY9vxDaaAk6dB4ZqwkNc6gl6TGGfSS1DiDXpIaZ9BLUuMMeklqnLcSXII8Xl6LwePrF48jeklqnEEvSY0z6CWpcQa9JDXOoJekxnnUTcM8ukZHAo/GWXiO6CWpcQa9JDXOoJekxjlH3wDn4iXNxqCX9ITkl7Tzx6kbSWqcQS9JjXPq5gjiXLzklM7BcEQvSY0bNKJPsh64HFgGvK+q/mDa+vTrzwUeBS6sqs8N2VY/ylG7dHBm+7+z1Ef7cwZ9kmXAlcBrgClge5KtVXXnSLUNwJr+cTpwFXD6wG2XJANdOnyW+nTPkBH9OmCyqu4BSHItsBEYDeuNwDVVVcDNSZYnWQGsHrBtEwxu6chzoP9vj9Q3hiFBvxLYPbI8RTdqn6vOyoHbApBkE7CpX3wkya4BbVtIxwPfWOQ2HG72eWmwzwcpfzgPLVk4z59pxZCgz5iyGlhnyLZdYdUWYMuA9hwWSXZU1drFbsfhZJ+XBvu89AwJ+inghJHlVcCegXWOHrCtJGkBDTm8cjuwJsmJSY4GzgO2TquzFbggnTOAh6rqvoHbSpIW0Jwj+qral+RS4Ea6QySvrqqdSS7u128GttEdWjlJd3jlRbNtuyA9mX9PmGmkw8g+Lw32eYlJd6CMJKlVnhkrSY0z6CWpcUs+6JOckOSvktyVZGeSt/Xlz0ryiSRf6n/+2GK3db4keWqSzyb5Qt/nd/XlzfZ5vyTLknw+yQ39ctN9TnJvkjuS3JZkR1/Wep+XJ7k+yd39/+tXtN7nuSz5oAf2Ab9eVT8FnAFckuQk4DLgk1W1Bvhkv9yK7wNnVdXJwCnA+v5oqZb7vN/bgLtGlpdCn19dVaeMHEfeep8vBz5WVS8BTqb7e7fe59lVlY+RB/Df6a7NswtY0ZetAHYtdtsWqL9PBz5Hd8Zy032mO4/jk8BZwA19Wet9vhc4flpZs30GjgO+Qn+gyVLo85CHI/oRSVYDpwK3AM+u7lwA+p8/sYhNm3f9FMZtwAPAJ6qq+T4D7wXeAfxwpKz1Phfw8SS39pcZgbb7/AJgL/CBforufUmOoe0+z8mg7yU5FvgL4O1V9fBit2ehVdUPquoUulHuuiQvW+QmLagkrwceqKpbF7sth9mZVXUa3RVmL0nyysVu0AJ7EnAacFVVnQp8l6U2TTOGQQ8keTJdyH+4qv6yL/56fwVO+p8PLFb7FlJVPQjcBKyn7T6fCbwxyb3AtcBZST5E232mqvb0Px8APkp3NdqW+zwFTPWfUAGupwv+lvs8pyUf9P1NU94P3FVV7xlZtRV4c//8zXRz901IMpFkef/8acA5wN003Oeq+u2qWlVVq+kuxfGpqnoTDfc5yTFJnrH/OfBa4Is03Oequh/YneTFfdHZdJdFb7bPQyz5M2OT/CPgfwN38Pdzt++km6e/Dnge8LfAz1XVtxalkfMsyT8APkh3WYqjgOuq6veS/DiN9nlUklcBv1FVr2+5z0leQDeKh25K40+r6t0t9xkgySnA++guqngP3SVZjqLhPs9lyQe9JLVuyU/dSFLrDHpJapxBL0mNM+glqXEGvSQ1zqCXpMYZ9JLUuP8HfU4LPHDt+/MAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.hist(S, density = True, bins=np.arange(S.min(), S.max()+1))\n",
    "plt.title('Histogram of sum score S')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 1.0, 'Histogram of squared sum score $S^2$')"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYcAAAEMCAYAAAAvaXplAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAa7ElEQVR4nO3df7RlZX3f8fcnww8NSToQLwSGQUgdXR2jK6FTwBUTLWpkSJOLRiJqAlKyJqwwTdI2lSFNUl3LtGNWG1MiYUoMCcTGCU1inMQx1JASl1aUQRAZyNQrIoxMYCBCRBIR+PaP84yeOfvcc/e9c+HOj/drrbPO3s9+nn2e57kz93P3Pvvsk6pCkqRh37LUHZAk7X8MB0lSh+EgSeowHCRJHYaDJKnDcJAkdRgOkqQOw0HSISfJy5J8IslfJ3l/ksOXuk/7G8NB0qHoi8CZVfUK4G5geon7s98xHA5ySbYneeVS92OpJHlRkluTfCXJzy51fxZDkt9L8s6l7seBrKrur6p/aKtPAk8vZX/2R4bDASzJPUlePVL21iQf27NeVS+uqhvnu5+DyNuAG6vq26vq8qXujJ49Sd6S5JYkjyT5uyQ3JlkxUucUYC3w50vTy/2X4aBnVJLDlrgLzwe2L3EfZrUfzM9+byFzlOStwDuAdcDRwAuB9wFfHqrzHcA1wE9W1ROL0tmDiOFwkBs+KkhyaZIvtVMsO5K8KsnvAycBf5bksSRva3X/WftL65F2aupHh/Z56tCpmv+V5A+HT3O017w0ye3AV5MclmRDks+3Nncmed1I/f+Q5PYkX03yO0mOS/LhVv8vkxw9YYxj+5rkr4B/Cbynje2FY9p25qSVf1+ST7fyP0yyeWSMleQFQ+vfONUzaawT5ueEJH+cZHeSLwyfAhvtC/CcOX7ms41pUp/39Wcw22uuTPInbVwPJ3nPXD+3hczRGBcB/6OqbqmBh6rqvVX1eNv/YcD7gbdX1Y5J83nIqiofB+gDuAd49UjZW4GPjdYBXgTcB5zQyk8G/um4/QCHAzPALwJHAGcCX2n7OILBm3k/1+q9HngCeOfIa94GrASe28rOBU5g8AfJG4GvAscP1b8JOA5YATwIfBr4PuBI4K+A/zTLHMza17b9RuCnZmk7dk6Gxvhv2/7fAHx9ZIwFvGBo/ff2bJ801nHz0+rdAvxKe+3vZvAm6Wv79KXPmHr0eV9+BrPN4zLgM8C7gaMYhNrLe/7ces/RLH3638BdwI8Dzxuz/SeBh9q/jxuBNy71/+f97bHkHfCxDz+8wX+gx4BHhh6PMz4cXtD+w78aOHzMfobD4QeAvwW+Zajs/cDbgR8EvgRkaNvH6IbDv56j77cB00P13zK07Y+BK4fW/w3wp7PsZ9a+tuUbmT0cxs5JG+P9I2P8v/QMh0ljHTc/wOnAvSNtLgN+t09f+oxprj7v489gtnl8GbAbOGwBP7feczRLn74L+A3gXuAp4M+AY5+t/5sHw8PTSge+c6pq+Z4H8DPjKlXVDPDzDH7BP9hOk5wwyz5PAO6rquErOL7I4C/KE4AvVfsf2Nw3Zh97lSU5P8lt7RTCI8D3AM8bqvLA0PI/jFn/tgX0daIJczJujF+ca3979Bgr7D0/zwdO2FO/tflFBn/Fz6sv8/w5j1rQz2DCa64EvlhVT45p1ufn1neOxvXpb6vq56vqJOA04KXApePqajzD4RBSVX9QVS9n8B+tgHft2TRS9X5gZZLhfx8nMThi2AWsSJKhbSvHvdyehSTPB34bWA98ZwuxO4CMaTdfk/o6p1nmZNwYTxpp+jjwrUPr3wXzGutouH5hOORrcHXV2T370mdMs/Z5MczymvcBJ2X8G8p9fm5952iuvt0CfJbBqS31ZDgcIjK43v/MJEcC/8jgL8Gn2uYHGJzD3eOTDM6Tvy3J4Rl8TuJHgM3AJ1q79e1NwmkGf5lNchSD/+i7W18uZPDX9GKY1NeJJszJJxhc+/6zbYyvpzvG24A3J1mW5CzgFa18IWP9FPD37Q3Y57Z9fk+Sf9GzL33GNKnP+2TCa36KQbhtTHJUkuck+f7WbL4/t0lzNNqfDUm+P8mR7fFW4JUMTtOpJ8Ph0HEksJHBm3B/CxzL4LAc4L8Av9QO13+hBpf1/SiD678fAn4LOL+q/qZtez2Dq0EeAX6CwTXiX5vthavqTuC/MfhF9wDwEuDjizGoSX3t0XzsnAyN8a0MLn18I/AnI21/jsEvs0eAtwB/2voz77FW1VNtX98LfKH1573AP+nZlznHNKnPi2C2edwzrhcwOPe/s/V/3j+3SXM0pvp3MAiCh9vrvhF4VVV9ct+HeujI3qcypflL8klgU1UdtH+ZJfk9YGdV/dJS90V6NnjkoHlL8ook39VOc1zA4M2+v1jqfklaPH46UwvxIuA6BlevfB54Q1XtWtouSVpMnlaSJHX0Oq2U5Kz2kfiZJBvGbE+Sy9v225OcOlfbJOe2j8w/nWTNmH2elMEtD35hoYOTJC3MnOGQZBlwBYOrClYDb0qyeqTaWmBVe6wDruzR9g4GV2F8dJaXfjfw4fkMRpK0OPq853AaMFNVdwMk2czgizHuHKozDVzbPsV5U5LlSY5ncI+VsW2r6q5W1nnBJOcwuG/KV/sM4nnPe16dfPLJfapKkppbbrnloaqaGretTzisYO+Pse9kcJ+Tueqs6Nl2L0mOYvAx99cAs55SSrKOwVEKJ510Etu2bZs4CEnS3pLMeiuWPu85jLvFwei72LPV6dN21DuAd1fVY5MqVdVVVbWmqtZMTY0NPknSAvU5ctjJ3vfOOZHBfVH61DmiR9tRpwNvSPJrwHLg6ST/WFXvmdxMkrRY+oTDzcCqDL5O70vAecCbR+psYXCvnc0Mfrk/WlW7kuzu0XYvVfUDe5aTvB14zGCQpGfXnOFQVU8mWQ9cz+DLO66uqu1JLm7bNwFbgbMZfHnH48CFk9oCZPDtWL8JTAEfSnJbVb12sQcoSZq/g+JDcGvWrCnfkJak+UlyS1V1PmcG3ltJkjSG4SBJ6jAcJEkdhoMkqcNbduuAcPKGD31j+Z6NP7yEPZEODR45SJI6DAdJUofhIEnq8D0H7ReG31OA+b2v4PsR0uLzyEGS1GE4SJI6DAdJUofhIEnqMBwkSR2GgySpw0tZtWRGL199Jvbppa3SwnjkIEnqMBwkSR2GgySpw3CQJHUYDpKkDsNBktTRKxySnJVkR5KZJBvGbE+Sy9v225OcOlfbJOcm2Z7k6SRrhspfk+SWJJ9tz2fu6yAlSfMzZzgkWQZcAawFVgNvSrJ6pNpaYFV7rAOu7NH2DuD1wEdH9vUQ8CNV9RLgAuD35z8sSdK+6PMhuNOAmaq6GyDJZmAauHOozjRwbVUVcFOS5UmOB06erW1V3dXK9nqxqrp1aHU78JwkR1bV1xYwPknSAvQ5rbQCuG9ofWcr61OnT9tJfgy4dVwwJFmXZFuSbbt3757HLiVJc+kTDhlTVj3r9Gk7/kWTFwPvAn563Paquqqq1lTVmqmpqT67lCT11Oe00k5g5dD6icD9Pesc0aNtR5ITgQ8A51fV53v0UerFrxSV+ukTDjcDq5KcAnwJOA9480idLcD69p7C6cCjVbUrye4ebfeSZDnwIeCyqvr4fAaj/Y+/jKUD05ynlarqSWA9cD1wF3BdVW1PcnGSi1u1rcDdwAzw28DPTGoLkOR1SXYCLwM+lOT6tq/1wAuAX05yW3scuzjDlST10euW3VW1lUEADJdtGlou4JK+bVv5BxicOhotfyfwzj79kiQ9M/yEtCSpw3CQJHUYDpKkDsNBktRhOEiSOgwHSVKH4SBJ6jAcJEkdhoMkqcNwkCR1GA6SpA7DQZLUYThIkjoMB0lSR69bdkt9DX+5z4FmtO9+OZEOZR45SJI6DAdJUofhIEnqMBwkSR2GgySpw3CQJHUYDpKkjl7hkOSsJDuSzCTZMGZ7klzett+e5NS52iY5N8n2JE8nWTOyv8ta/R1JXrsvA5Qkzd+c4ZBkGXAFsBZYDbwpyeqRamuBVe2xDriyR9s7gNcDHx15vdXAecCLgbOA32r7kSQ9S/ocOZwGzFTV3VX1BLAZmB6pMw1cWwM3AcuTHD+pbVXdVVU7xrzeNLC5qr5WVV8AZtp+JEnPkj7hsAK4b2h9ZyvrU6dP24W8HknWJdmWZNvu3bvn2KUkaT76hEPGlFXPOn3aLuT1qKqrqmpNVa2ZmpqaY5eSpPnoc+O9ncDKofUTgft71jmiR9uFvJ4k6RnU58jhZmBVklOSHMHgzeItI3W2AOe3q5bOAB6tql09247aApyX5MgkpzB4k/tT8xiTJGkfzXnkUFVPJlkPXA8sA66uqu1JLm7bNwFbgbMZvHn8OHDhpLYASV4H/CYwBXwoyW1V9dq27+uAO4EngUuq6qlFHbUkaaJe3+dQVVsZBMBw2aah5QIu6du2lX8A+MAsbX4V+NU+fZMkLT4/IS1J6jAcJEkdfk2o1INfIapDjeGgeRv+RekvSeng5GklSVKH4SBJ6jAcJEkdhoMkqcNwkCR1GA6SpA7DQZLUYThIkjoMB0lSh+EgSeowHCRJHYaDJKnDcJAkdRgOkqQOw0GS1GE4SJI6DAdJUofhIEnq6BUOSc5KsiPJTJINY7YnyeVt++1JTp2rbZJjknwkyefa89Gt/PAk1yT5bJK7kly2GAOVJPU353dIJ1kGXAG8BtgJ3JxkS1XdOVRtLbCqPU4HrgROn6PtBuCGqtrYQmMDcClwLnBkVb0kybcCdyZ5f1XdszhD1nwNf2e0ukbnx+/V1sGgz5HDacBMVd1dVU8Am4HpkTrTwLU1cBOwPMnxc7SdBq5py9cA57TlAo5KchjwXOAJ4O8XNDpJ0oL0CYcVwH1D6ztbWZ86k9oeV1W7ANrzsa38j4CvAruAe4H/WlV/N9qpJOuSbEuybffu3T2GIUnqq084ZExZ9azTp+2o04CngBOAU4B/n+S7Ozupuqqq1lTVmqmpqTl2KUmajz7hsBNYObR+InB/zzqT2j7QTj3Rnh9s5W8G/qKqvl5VDwIfB9b06KckaZH0CYebgVVJTklyBHAesGWkzhbg/HbV0hnAo+1U0aS2W4AL2vIFwAfb8r3AmW1fRwFnAH+zwPFJkhZgzquVqurJJOuB64FlwNVVtT3JxW37JmArcDYwAzwOXDipbdv1RuC6JBcxCIRzW/kVwO8CdzA4LfW7VXX7YgxWktTPnOEAUFVbGQTAcNmmoeUCLunbtpU/DLxqTPljfDMoJElLwE9IS5I6DAdJUofhIEnqMBwkSR2GgySpw3CQJHUYDpKkDsNBktRhOEiSOgwHSVKH4SBJ6jAcJEkdvW68p4Pf8Pcg+x3Ii8u51YHIIwdJUofhIEnqMBwkSR2GgySpw3CQJHUYDpKkDsNBktRhOEiSOgwHSVJHr3BIclaSHUlmkmwYsz1JLm/bb09y6lxtkxyT5CNJPteejx7a9tIkn0iyPclnkzxnXwcqSepvznBIsgy4AlgLrAbelGT1SLW1wKr2WAdc2aPtBuCGqloF3NDWSXIY8D7g4qp6MfBK4OsLH6Ikab76HDmcBsxU1d1V9QSwGZgeqTMNXFsDNwHLkxw/R9tp4Jq2fA1wTlv+IeD2qvoMQFU9XFVPLWx4kqSF6BMOK4D7htZ3trI+dSa1Pa6qdgG052Nb+QuBSnJ9kk8nedu4TiVZl2Rbkm27d+/uMQxJUl99wiFjyqpnnT5tRx0GvBx4S3t+XZJXdXZSdVVVramqNVNTU3PsUpI0H33CYSewcmj9ROD+nnUmtX2gnXqiPT84tK+/rqqHqupxYCtwKpKkZ02fcLgZWJXklCRHAOcBW0bqbAHOb1ctnQE82k4VTWq7BbigLV8AfLAtXw+8NMm3tjenXwHcucDxSZIWYM4v+6mqJ5OsZ/BLexlwdVVtT3Jx276JwV/3ZwMzwOPAhZPatl1vBK5LchFwL3Bua/PlJL/OIFgK2FpV3/y2FEnSM67XN8FV1VYGATBctmlouYBL+rZt5Q8DnfcS2rb3MbicVZK0BPyEtCSpw++QPkQNf6+xJI0yHKQlNBrS92z84SXqibQ3TytJkjoMB0lSh+EgSeowHCRJHYaDJKnDcJAkdRgOkqQOw0GS1GE4SJI6DAdJUofhIEnqMBwkSR2GgySpw3CQJHUYDpKkDsNBktThl/0cQvz2N0l9GQ7SfspvidNS8rSSJKmjVzgkOSvJjiQzSTaM2Z4kl7fttyc5da62SY5J8pEkn2vPR4/s86QkjyX5hX0ZoCRp/uYMhyTLgCuAtcBq4E1JVo9UWwusao91wJU92m4AbqiqVcANbX3Yu4EPL2BMkqR91OfI4TRgpqrurqongM3A9EidaeDaGrgJWJ7k+DnaTgPXtOVrgHP27CzJOcDdwPYFjUqStE/6hMMK4L6h9Z2trE+dSW2Pq6pdAO35WIAkRwGXAu+Y1Kkk65JsS7Jt9+7dPYYhSeqrTzhkTFn1rNOn7ah3AO+uqscmVaqqq6pqTVWtmZqammOXkqT56HMp605g5dD6icD9PescMaHtA0mOr6pd7RTUg638dOANSX4NWA48neQfq+o9PfoqSVoEfY4cbgZWJTklyRHAecCWkTpbgPPbVUtnAI+2U0WT2m4BLmjLFwAfBKiqH6iqk6vqZOA3gP9sMEjSs2vOI4eqejLJeuB6YBlwdVVtT3Jx274J2AqcDcwAjwMXTmrbdr0RuC7JRcC9wLmLOjJJ0oL1+oR0VW1lEADDZZuGlgu4pG/bVv4w8Ko5XvftffonSVpcfkJaktRhOEiSOgwHSVKH4SBJ6vCW3dIByNt565lmOBxE/IUhabF4WkmS1GE4SJI6DAdJUofhIEnqMBwkSR2GgySpw3CQJHUYDpKkDsNBktRhOEiSOgwHSVKH4SBJ6vDGewe40ZvtSbD3vwtvwKiF8MhBktRhOEiSOgwHSVJHr3BIclaSHUlmkmwYsz1JLm/bb09y6lxtkxyT5CNJPteej27lr0lyS5LPtuczF2OgkqT+5gyHJMuAK4C1wGrgTUlWj1RbC6xqj3XAlT3abgBuqKpVwA1tHeAh4Eeq6iXABcDvL3h0kqQF6XPkcBowU1V3V9UTwGZgeqTONHBtDdwELE9y/Bxtp4Fr2vI1wDkAVXVrVd3fyrcDz0ly5MKGJ0laiD7hsAK4b2h9ZyvrU2dS2+OqahdAez52zGv/GHBrVX1tdEOSdUm2Jdm2e/fuHsOQJPXVJxwypqx61unTdvyLJi8G3gX89LjtVXVVVa2pqjVTU1N9dilJ6qlPOOwEVg6tnwjc37POpLYPtFNPtOcH91RKciLwAeD8qvp8jz5KkhZRn3C4GViV5JQkRwDnAVtG6mwBzm9XLZ0BPNpOFU1qu4XBG8605w8CJFkOfAi4rKo+vvChSZIWas7bZ1TVk0nWA9cDy4Crq2p7kovb9k3AVuBsYAZ4HLhwUtu2643AdUkuAu4Fzm3l64EXAL+c5Jdb2Q9V1TeOLA5l3i5D+2r035C319A4ve6tVFVbGQTAcNmmoeUCLunbtpU/DLxqTPk7gXf26Zck6ZnhJ6QlSR2GgySpw3CQJHUYDpKkDsNBktRhOEiSOvya0P2cn2vQs8nPQGgPjxwkSR2GgySpw3CQJHUYDpKkDsNBktTh1Ur7Aa8Q0YHAf6eHFo8cJEkdhoMkqcNwkCR1+J7DEvGTz5L2Z4aDpEXnm9cHPk8rSZI6DAdJUoenlZ4lvscgfdPw/wdPOe2fPHKQJHX0OnJIchbw34FlwHurauPI9rTtZwOPA2+tqk9PapvkGOAPgZOBe4Afr6ovt22XARcBTwE/W1XX79Mol4BHCtLCeFSxf5gzHJIsA64AXgPsBG5OsqWq7hyqthZY1R6nA1cCp8/RdgNwQ1VtTLKhrV+aZDVwHvBi4ATgL5O8sKqeWpwhP3MMBOnZZZA8c/ocOZwGzFTV3QBJNgPTwHA4TAPXVlUBNyVZnuR4BkcFs7WdBl7Z2l8D3Ahc2so3V9XXgC8kmWl9+MTCh7k4/OUvHTjmczmtIdOVwe/zCRWSNwBnVdVPtfWfBE6vqvVDdf4c2FhVH2vrNzD4RX/ybG2TPFJVy4f28eWqOjrJe4Cbqup9rfx3gA9X1R+N9GsdsK6tvgjYscA5ONg8D3hoqTuxn3OOJnN+JjuY5uf5VTU1bkOfI4eMKRtNlNnq9Gm7kNejqq4CrppjX4ecJNuqas1S92N/5hxN5vxMdqjMT5+rlXYCK4fWTwTu71lnUtsH2qkn2vOD83g9SdIzqE843AysSnJKkiMYvFm8ZaTOFuD8DJwBPFpVu+ZouwW4oC1fAHxwqPy8JEcmOYXBm9yfWuD4JEkLMOdppap6Msl64HoGl6NeXVXbk1zctm8CtjK4jHWGwaWsF05q23a9EbguyUXAvcC5rc32JNcxeNP6SeCSA+FKpf2Ip9rm5hxN5vxMdkjMz5xvSEuSDj1+QlqS1GE4SJI6DIf9XJKrkzyY5I6hsmOSfCTJ59rz0UPbLksyk2RHktcOlf/zJJ9t2y5vtzw5KCRZmeT/JLkryfYkP9fKnScgyXOSfCrJZ9r8vKOVOz9DkixLcmv73JbzU1U+9uMH8IPAqcAdQ2W/BmxoyxuAd7Xl1cBngCOBU4DPA8vatk8BL2PwOZIPA2uXemyLOEfHA6e25W8H/l+bC+dpMKYA39aWDwc+CZzh/HTm6d8BfwD8eVs/pOfHI4f9XFV9FPi7keJpBrccoT2fM1S+uaq+VlVfYHD12GntcyTfUVWfqMG/4GuH2hzwqmpXtRs9VtVXgLuAFThPANTAY2318PYonJ9vSHIi8MPAe4eKD+n5MRwOTMfV4HMktOdjW/kK4L6hejtb2Yq2PFp+0ElyMvB9DP46dp6adsrkNgYfNv1IVTk/e/sN4G3A00Nlh/T8GA4Hl8W8jckBJ8m3AX8M/HxV/f2kqmPKDup5qqqnqup7Gdxx4LQk3zOh+iE1P0n+FfBgVd3St8mYsoNufgyHA9N8bz2ysy2Plh80khzOIBj+Z1X9SSt2nkZU1SMM7oB8Fs7PHt8P/GiSe4DNwJlJ3schPj+Gw4FpXrceaYfEX0lyRrt64vyhNge8NqbfAe6qql8f2uQ8AUmmkixvy88FXg38Dc4PAFV1WVWdWFUnM7jFz19V1U9wqM/PUr8j7mPyA3g/sAv4OoO/TC4CvhO4Afhcez5mqP5/ZHD1xA6GrpQA1gB3tG3voX06/mB4AC9ncPh+O3Bbe5ztPH1jTC8Fbm3zcwfwK63c+enO1Sv55tVKh/T8ePsMSVKHp5UkSR2GgySpw3CQJHUYDpKkDsNBktRhOEiSOgwHSVLH/weKczN00gDDTAAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.hist(S2, density = True, bins=np.arange(S2.min(), S2.max()+1, 50))\n",
    "plt.title('Histogram of squared sum score $S^2$')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 6. Calculate mean values $\\bar{S}$ and $\\bar{S^2}$ and compare to expected values $\\mathbb{E}(S)$ and $\\mathbb{E}(S^2)$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Mean of S   = 42.001402; cp. with expected value = 42.0\n",
      "Mean of S^2 = 1799.079366; cp. with expected value = 1799.0\n"
     ]
    }
   ],
   "source": [
    "mean_S  = np.mean(S)\n",
    "mean_S2 = np.mean(S2)\n",
    "E_S  = n*7/2\n",
    "E_S2 = n*91/6 + n*(n-1)*49/4\n",
    "Var_S = E_S2 - E_S**2\n",
    "print('Mean of S   = {}; cp. with expected value = {}'.format(mean_S , E_S))\n",
    "#print('Standard error of mean (of S) = {}'.format(np.sqrt(Var_S/M))\n",
    "print('Mean of S^2 = {}; cp. with expected value = {}'.format(mean_S2, E_S2))"
   ]
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "01da3e92184f855baee79aee4095f9eaf286d7f3cd3041cddfb50dc3e6722963"
  },
  "kernelspec": {
   "display_name": "R",
   "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.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}