{ "cells": [ { "cell_type": "markdown", "metadata": { "id": "7Of9cveNTErL" }, "source": [ "\n", "HEADS AND TAILS\n", "\n", "Abstract\n", "\n", "Science is built upon the Scientific Method, which offers a procedure to formulate and accept or reject propositions about Nature. Randomness is in Science the first explanation of whatever phenomenon. Here we show how this operates when we have two events or outputs produced by a coin, which we call heads and tails. Heads may correspond to being male, or to belonging to the right wing in politics, or to being older than 30 years, or to being indifferent to the cries for help from a country under unjust attack, or to a favorable mutation. In the literature, this theme is known as the Binomial Distribution. If your challenge is to learn or teach the concepts behind statistical formulae and procedures, try this material. It can also be used to learn Python 3, a widespread programming language that is very popular among scientists, and how to work with AI to improve and develop software. The concepts shown here are universal, and the whole of science is more of the same.\n", "\n", " # INTRODUCTION\n", "\n", "I visited a family yesterday, and in that house there were a cat and a dog. Both were boosted by curiosity, but the cat directed its silent attention to the things brought by me, while the dog was continuously barking at me until we finally became friends. From their behavior, we see that they have a neural machinery to distinguish usual and familiar from unusual and unfamiliar. Human babies can exercise this capability as early as their first year of life. So, we train ourselves to compare things of whatever nature and to make ensuing decisions. Science is a specialization of this capability.\n", "\n", "***The purpose of Science is to compare what we think with what we observe to polish our thoughts about nature.***\n", "\n", "> **Example.** The first explanation in Science of whatever phenomenon is randomness. If it does not function, we try to invent another explanation. So, randomness explains gravity as follows:\n", ">>\n", ">> Things move randomly by their very nature, and they happen to fall towards the Earth because by mere casualty their random movement happens to be aligned downwards. But this may change at whatever moment, and things could begin to levitate, or to slide horizontally.\n", ">>\n", ">> Now we, human beings, do not believe in this explanation because we expect from random movement to be random, i.e., to vary without law making it impossible to predict what will be next. But, this is not what happens while bodies fall or when they stand still. In fact, we predict with infallible success that free bodies fall, move towards the Earth. More to the point, we have gathered more data and believe that in the nearness of the Earth, the theory of Newton is compatible with observation: bodies attract one to another with a force that is proportional to the product of their masses and inversely to their distance squared. For planetary data and beyond, as in the nearness of a black hole, the General Relativity Theory of Einstein is much better, in which acceleration results from the curvature of space-time.\n", "\n", "\n", "Special attention is paid to the case in which the difference between observed and expected can be measured. If it is low, Science advises in first instance to continue believing what we think but if that difference is high, Science advises considering other thoughts, other explanations. Notice that **Science is an advisor not a dogmatic dictator. It is human beings that accept or reject the advises of Science and take decisions under their sole responsibility**.\n", "\n", "An immediate question arises: *How can we compare what we think, something subjective, with the objective reality observed by Science?*\n", "\n", "Thoughts are organized in the form of a theory and this is forced to produce observable predictions. So, we end comparing predicted results with observed events or, simply, predicted with observed. For instance, we have said that at the time of a visit to a household, dogs pay more attention to human beings than to things while cats are the other way around. From this we derive a prediction: thefts prefer to robber apartments where there are cats than where there are dogs. We have a theory and a reasonable prediction. The next step is to compare predicted with observed. If we find that robbery of houses happens much more frequently where there were cats and no dogs, we can continue believing that dogs serve as guardians much better than cats.\n", "\n", "Decision-making is easier when things are clear than when data are blurred. In fact, a crucial problem with observed data is that they are not, in general, deterministic, clearly established. Instead, they come with variations that could be difficult to explain and that can be treated as coming at random, i.e., difficult to fit into a deterministic pattern. But we need to make decisions every day, maybe every hour. Think just of the weather that shows us that decision-making in the face of randomness is a truly difficult enterprise that is stained with undesirable failure.\n", "\n", "That is why our purpose is to work out the simplest case of decision-making under random data and to show the underlined formalism that is usual in Science. We consider that, regarding basic principles and procedures, all Science is more of the same, of what we will see here. We will work with a coin, and its behavior will be simulated in Python, a very popular programming language among scientists. Running programs is a duty here because to contrast our beliefs with what we observe, we need to calculate the predictions of our beliefs, a task that is fulfilled by our simulations. Anyway, all Python programs can be run with one click from this notebook thanks to the free Colab service from Google. One click is all you need to run a program.\n", "\n", "***Let us discuss now an important question: how can we make decisions in the face of randomness, given that we cannot understand its behavior? The answer is: although we do not know what will come next, we can estimate what will be better in the long run, what tendencies will arise.***\n", "\n", "\n", "> **Example.** Let us consider a game with a fair die, one whose geometry is perfectly symmetric. The game consists of betting for the best of two options. The first is to bet for Six, while the second is to bet for Five or Less. We immediately choose the second bet because its probability is 5 times that of the first. Notice that we do this in spite of the fact that we do not know what will fall next. But we hope that in the long run, after betting many times, the second option will pay off.\n", "\n", "In general, we can say: to make decisions in the face of randomness, we must look for tendencies towards deterministic observables and exploit them to our profit (an observable is a variable that can be measured or a result of their mathematical manipulation). The tendency in the example of dice is that after many trials, the relative proportion of Five or Less will be 5 times larger than that of Six.\n", "\n", "A very important example of tendencies follows.\n", "\n", "## Theorem of the self-annihilation of deviations from the mean.\n", "\n", "Making decisions always imposes the need to represent a population of data by a sample of maybe by an abstract number, and this is done according to the purposes in mind.\n", "\n", "> **Example.** We need an image of a group of pupils, and all we have is their performance in a math quiz. If our purpose is to impress, we will take the best score. If our purpose is to point to the deficient work done by the teacher, we will take the lowest score. And if we look for the best representative, we will report the mean, which is an abstract number.\n", "\n", "The mean has two fundamental virtues: it represents the population democratically, giving the same weight to all data, and it minimizes the inner tensions created by being different from the mean. This means that there is a discrepancy between the mean and each data point, but all discrepancies cancel one another. Formally:\n", "\n", "\n", "> **Theorem.** Let us consider a data set consisting of numbers $a_1, a_2,...,a_n$ with mean $\\bar{a}$. Then, the deviations from the mean, $\\bar{a} - a_i$, add up to zero.\n", ">\n", ">>Proof. The mean is found by adding all terms and dividing the result by the number of terms:\n", ">>\n", ">>>$\\bar{a} = \\sum^n_{i=1} a_i = \\frac{a_1 + a_2 + ... + a_n}{n}$\n", ">>\n", ">>The deviations from the mean are $\\bar{a} - a_1 $, $\\bar{a} - a_2 $ ... $\\bar{a} - a_n $. So, their sum $S$ is:\n", ">>\n", ">>>$ S = \\bar{a} - a_1 + \\bar{a} - a_2 + ... +\\bar{a} - a_n $\n", ">>\n", ">>Let us associate the means, $n$ at all, and put the other summands inside a parenthesis:\n", ">>>\n", ">>>$ S = n\\bar{a} - (a_1 + a_2+...a_n) $\n", ">>>\n", ">>>$ = n \\frac{a_1 + a_2 + ... + a_n}{n}- (a_1 + a_2+...a_n)$\n", ">>>\n", ">>>$ = (a_1 + a_2+...a_n) - (a_1 + a_2+...a_n) = 0$.\n", ">>\n", "> We have shown that the deviations of the data from the mean add up to 0.\n", "\n", "This theorem is valid for every kind of numeric data, in particular for random data. Now, when the source of data has some degree of stability or homogeneity, the mean is also stable in the sense that its variations from place to place or from moment to moment can be ignored, to first approximation. In that case, the stability grows when more data are added, and the mean becomes more and more determined: the mean of a sample has a tendency to get closer to the mean of the population. All this explains why the mean appears everywhere in Science and decision-making. And this is also why everyone can commit serious errors, say, if one gathers biased data that represents a subpopulation but takes decisions that affect the whole inhomogeneous population.\n", "\n", "> **Example.** An airline promises that, as a strong commitment to its service, the flight will last 1 hour and 40 minutes. How can such an exact promise be formulated, when everyone knows that the duration of the flight has variations? This happens so: the airline finds the mean duration, and the deviations from the mean are found to be equal to or less than 5 minutes. Thus, no one protests against the propaganda given by the airline. The mean has here a great predictive power. But not always does the mean represent the population so good. For instance, if one travels through mountains by car, durations are given in another format. We prefer to say: the trip takes from 8 to 12 hours, and the mean is not mentioned because it has low predictive power. \n", "\n", "\n", "Dealing with coins, we will find means and relative frequencies as our principal observables with a tendency to stability, to determinism, to population values. \n", "\n", "# COINS AS A MODEL\n", "\n", "A coin has two faces that are generically called heads and tails. People are used to the belief that it is fair to bet on any face. In our words, we have a belief: the probability of heads equals that of tails equals 0.5. Probability is here an indication of how fair it is to bet for a given face, so when we say that the probability of heads equals that of tails, we are simply saying that it is fair to bet for any one of the faces of a coin. We refer to heads and tails as events. Betting to flipping coins is a popular game because one cannot predict with certainty what will land next, whether heads else tails. We say that for whatever flip, the coins land at random and with the same probability. We are also assuming that flips are independent, that one flip does not change future relative frequencies, anyhow. The *relative frequency* of heads is the number of heads divided by the total number of trials or flips.\n", "\n", "**Let us frame our discussion inside an aforementioned intriguing question: if coins land at random, and uncertainty is dominant, how can we make Science, how can we formulate verbal propositions that enlighten decision-making?**\n", "\n", "To begin with, we notice that we have the preconception that it is fair to bet on heads or tails. Let us formulate this in the official jargon: our preconception is called a null hypothesis and denoted $H_o$, and it says that the probability of heads equals that of tails.\n", "\n", "The second step in the way of Science is to relate the abstract, subjective concept of probability, based on fairness, with objective observations. The immediate link is the relative frequency of heads, when a given number of independent repetitions is allowed. We hope, according to our $H_o$, that the relative frequency of heads falls near 0.5 and that this is the more certain the larger the number of repetitions. With this, we are claiming that *repetitions create tendencies in the relative frequencies that are reaffirmed when the number of repetitions increases*. Concretely, our expectancy reads:\n", "\n", "When independent repetitions are allowed, a tendency of the relative frequency of heads appears that tends to a limit that is its probability. Does this expectancy hold? Are there other tendencies to a limit? How can we use them in Science?\n", "\n", "> **Tendencies are useful in Science** to characterize an experimental setting. This allows us to compare differing settings, different experiments, and to conclude that they produce the same else different output, that they have the same or different behavior. But, how will we resolve the accusation that results always include randomness, uncertainty, and that therefore our conclusions could be banal or even misleading?\n", "\n", "Our answer to this accusation is to work in detail a complete example. Our case study is just that of a coin, and our belief to be tested, our $H_o$, is just that the probability of heads is 0.5. \n", "\n", "## RELATIVE FREQUENCY AS A TENDENCY\n", "\n", "*In this section, we investigate what happens to the relative frequency of heads when the number of flips of a coin increases.* Does it observe a tendency to a limit that coincides with the probability?\n", "\n", "This question must be answered experimentally. But, attention, we test beliefs that only exist in our minds. Our beliefs refer to experimental settings: a coin with two faces, one with a drawing and the other with a number, and both faces land upwards with the same probability 1/2 = 0.5. Moreover, events are independent, meaning that anyone flipping could be treated as the initial one, or equivalently, that a given tossing or result does not interfere with the output of future ones.\n", "\n", "*It is instructive but extremely difficult to finely control experiments with real coins. That is why we appreciate the possibility to make experiments within the memory of a computer with almost exactly the demanded conditions, the same preparation. Such experiments are called simulations.* **In the rest of this section, we will enlighten the concept of probability as an idealized relative frequency**. We will use the Python language for our simulations because Google Colab allows us to effortlessly run programs in this language. This is possible because in this Jupyterlab notebook, text and programs are in different cells, and Colab allows us to run program cells with one click.\n", "\n", "## Simulation of a coin\n", "\n", "Our first program simulates tossing a coin various times and registering the output.\n", "\n", "> **Example.** Let us simulate flipping a coin 7 times." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "id": "bBUW4XuYTimE" }, "outputs": [], "source": [ "# Program CoinSimulation\n", "import random\n", "\n", "coin = random.choice([\"Heads\", \"Tails\"])\n", "print(coin)\n", "\n", "coin = random.choice([\"Heads\", \"Tails\"])\n", "print(coin)\n", "\n", "coin = random.choice([\"Heads\", \"Tails\"])\n", "print(coin)\n", "\n", "coin = random.choice([\"Heads\", \"Tails\"])\n", "print(coin)\n", "\n", "coin = random.choice([\"Heads\", \"Tails\"])\n", "print(coin)\n", "\n", "coin = random.choice([\"Heads\", \"Tails\"])\n", "print(coin)\n", "\n", "coin = random.choice([\"Heads\", \"Tails\"])\n", "print(coin)" ] }, { "cell_type": "markdown", "metadata": { "id": "y3LDTQbKTqdg" }, "source": [ "\n", "> **Exercise.** Run the code, observe its output, and analyze it. *To run the code, open this document in Colab, drag the cursor over the program cell, and then click on the popping-up icon at the upper-left. To use Colab, type in Google: Colab. This will open a web page in which you can open and see Jupyter and JupiterLab documents and run with one click the Python programs inside it.*\n", ">>\n", ">> **Answer**. The code contains 7 repeats of the same block of two instructions. The second instruction prints the value of the coin, so the first one must be an instruction to simulate a flip.\n", "\n", "Let us see in more detail the meaning of the code.\n", "\n", "> **Example.** What does this code mean?\n", ">> The first instruction `# Program CoinSimulation` begins with the symbol # that is an indicator of a comment in Python. A comment is a message for the human client and for the AI. In our programs, the first comment contains the name of the program.\n", ">>\n", ">> The instruction `import random` commands to add a library with the code to generate random numbers. *A library is a toolbox of code with programs that have been thoroughly tested for truthfulness, doing what they are designed for, style, and efficiency.* These numbers are generated by mini-programs that work deterministically, which is why they are *pseudo-random*, but their output sequences are so chaotic that they produce statistically acceptable results.\n", ">>\n", ">> The instruction `coin = random.choice([\"Heads\", \"Tails\"])` commands to simulate flipping of a coin with two faces with their usual names and with equal probability. The result of the simulation of each flip is written to the variable coin.\n", ">>\n", ">> The instruction `print(coin)` commands to print the value that at the moment of printing has been assigned to the variable coin.\n", "\n", "It is important to notice that, just as in real life, the generator of random numbers produces results that may be different for each new call.\n", "\n", "> **Exercise.** Find a way to repeat the experiment and verify that the answer can vary from run to run. Occasionally, no change is produced. Is that problematic?\n", ">>\n", ">> **Answer.** To rerun a program, we must drag the cursor to the code cell and then click on the Run icon as many times as desired. No change after a new run is no problem because this also may happen by mere randomness, just as in real life.\n", "\n", "*Programs are functional portions of code.* They can be easily modified.\n", "\n", "> **Exercise.** Find a way to modify the code to simulate the tossing of a coin 14 times. Test your new code by running it.\n", ">>\n", ">> **Answer.** The code to simulate flipping a coin 7 times is copied to the clipboard. We use `CTRL + C`. It is then pasted at the end of the existing program with `CTRL + V`, or pasted twice into a new code cell, observing the same indentation as the previous lines. The result is a program that simulates 14 flips of a coin, a fact that can be verified by running the program. The code follows:\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "id": "clVFzVW925JY" }, "outputs": [], "source": [ "# Program CoinSimulation\n", "import random\n", "\n", "coin = random.choice([\"Heads\", \"Tails\"])\n", "print(coin)\n", "\n", "coin = random.choice([\"Heads\", \"Tails\"])\n", "print(coin)\n", "\n", "coin = random.choice([\"Heads\", \"Tails\"])\n", "print(coin)\n", "\n", "coin = random.choice([\"Heads\", \"Tails\"])\n", "print(coin)\n", "\n", "coin = random.choice([\"Heads\", \"Tails\"])\n", "print(coin)\n", "\n", "coin = random.choice([\"Heads\", \"Tails\"])\n", "print(coin)\n", "\n", "coin = random.choice([\"Heads\", \"Tails\"])\n", "print(coin)\n", "\n", "coin = random.choice([\"Heads\", \"Tails\"])\n", "print(coin)\n", "\n", "coin = random.choice([\"Heads\", \"Tails\"])\n", "print(coin)\n", "\n", "coin = random.choice([\"Heads\", \"Tails\"])\n", "print(coin)\n", "\n", "coin = random.choice([\"Heads\", \"Tails\"])\n", "print(coin)\n", "\n", "coin = random.choice([\"Heads\", \"Tails\"])\n", "print(coin)\n", "\n", "coin = random.choice([\"Heads\", \"Tails\"])\n", "print(coin)\n", "\n", "coin = random.choice([\"Heads\", \"Tails\"])\n", "print(coin)" ] }, { "cell_type": "markdown", "metadata": { "id": "SW19nB6I25ud" }, "source": [ "A wonderful surprise to enjoy: *Colab can explain the code.*\n", "\n", "\n", "> **Exercise.** Right-click the mouse over the cell with the code. Select \"Explain code\". Compare the offered explanation with that presented above. Give feedback to Colab: rate the offered answer as perfect else imperfect. If the Explain Code option isn't visible, look for the three-dot menu in the upper right corner of the active code cell. One of the options in this menu is Explain Code. Use it and enjoy!\n", ">>\n", ">> **Answer.** The explanation given by Colab is much better, more understandable, and nicer than the one offered in this text. Therefore, it is a good idea to use the gentleness of Colab routinely.\n", "\n", "\n", "To repeat the experiment 14 times, the code was repeated 14 times. Clearly, this childish strategy does not function for simulations with hundreds of repetitions. In all programming languages, repetitions are implemented by loops.\n", "\n", "> **Example.** The following code uses a loop to repeat a task 7 times. The task is to simulate flipping a coin." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "id": "lx-tkHcrT6iD" }, "outputs": [], "source": [ "# Program CoinSimulationWithLoop.py\n", "import random\n", "\n", "# How many flips\n", "num_flips = 7; # the symbol underscore _ is here a separator.\n", "\n", "\"\"\"Loop to simulate multiple coin flips \"\"\"\n", "for i in range(num_flips):\n", " coin = random.choice([\"Heads\", \"Tails\"])\n", " # Print the values of i and coin\n", " print(i, coin)" ] }, { "cell_type": "markdown", "metadata": { "id": "VbXROFaQT76X" }, "source": [ "\n", "\n", "Comments:\n", "\n", "1. Observe the alternative way to mark a commentary, which is surrounded by triple quote marks \"\"\".\n", "2. The instruction `for i in range(num_flips):` defines the loop. The variable `i` is the loop counter that counts the number of executed loops from `0` to `num_flips`. The variable is incremented at the beginning of each loop, and the instruction says: execute the code in the block below, beginning with the counting variable with 0 until the variable gets equal to its upper limit, which is `num_flips = 7`; the variable `i` counts 0,1,2,...,6 and stops.\n", "\n", "> **Exercise.** Run the code and verify that the `i` variable runs from 0 to 6 and not to 7 because the variable is incremented before each loop, and if the variable attains its upper limit, looping must stop.\n", "\n", "> **Exercise.** Find a way to have the program `CoinSimulationWithLoop.py` simulate 14 flips.\n", ">>\n", ">> **Answer.** The variable that defines the number of flips must be set to 14: `num_flips = 14`;\n", "\n", "\n", "## The role of indentation in Python\n", "\n", "*Indentation refers to the blank spaces at the beginning of a code line.* Python structures the code using indentation. This means that indented code forms a block that is slaved to the first line above with less indentation. Thus, in the program `CoinSimulationWithLoop.py`, the for-loop repeatedly executes a block that consists of the indented code:\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "id": "qj-sWxFbUDEO" }, "outputs": [], "source": [ "coin = random.choice([\"Heads\", \"Tails\"])\n", "print(i, coin)" ] }, { "cell_type": "markdown", "metadata": { "id": "CTtjEsM4UGgO" }, "source": [ "The variable `num_flip` dictates the number of times that this block must be executed.\n", "\n", "\n", "> **Exercise.** Change at will the indentation of each line of this block inside the complete code to see the effect and explain it.\n", ">>\n", ">> **Answer.** The following change is viable and prints the result of flipping the coin once." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "id": "cSUEtr05ULu5" }, "outputs": [], "source": [ "for _ in range(num_flips):\n", " coin = random.choice([\"Heads\", \"Tails\"])\n", "print(coin)" ] }, { "cell_type": "markdown", "metadata": { "id": "eaEoodzlUdtg" }, "source": [ "\n", "> The reason is that the block that is under the scope of the loop consists of just one line, `coin = random.choice([\"Heads\", \"Tails\"])`, while the print instruction below is outside the loop. Therefore, Python executes 7 loops without printing and then prints the result of the last one. Other indentations produce an indentation error because the resultant code does not match the predefined syntax.\n", "\n", "## Underscore as a name of a variable\n", "\n", "The symbol `_` is named *underscore.* It is another letter in Python. It is used as the name of the counting variable for loops when the variable does not appear in the loop block. \n", "\n", "> **Example** The following code uses underscore as the loop variable because it does not appear inside the loop block." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "id": "OeJ-WnTLUetQ" }, "outputs": [], "source": [ "# Program CoinSimulationUnderscore\n", "import random\n", "\n", "# How many flips\n", "num_flips = 7;\n", "\n", "\"\"\"Loop to simulate multiple coin flips \"\"\"\n", "# The symbol _ is the name of the loop variable,\n", "# it is used when the loop variable does not appear in the slave code\n", "for _ in range(num_flips):\n", " coin = random.choice([\"Heads\", \"Tails\"])\n", " print(coin)" ] }, { "cell_type": "markdown", "metadata": { "id": "q_G28f5fUni8" }, "source": [ "\n", "> **Exercise.** Run this snippet or piece of code. What does this code mean?\n", ">>\n", ">> **Answer.** The code simulates flipping a coin 7 times. The result of each flip is printed because the corresponding instruction has the same indentation level as the previous one.\n", "\n", "> **Research.** Let us observe that some lines end with a semicolon and others do not. What is the use of semicolons in Python, and how recommended is it? Hint: type this question as a prompt into your browser.\n", "\n", "\n", "## Interactive code\n", "\n", " *An interactive program interacts with the User*.\n", "\n", " > **Example** The following example is provided by a member of the community. It runs perfectly over Colab + Chrome. With more difficulty over Colab + Visual Studio Code because the input-box appears at the top of the active window. " ] }, { "cell_type": "code", "execution_count": null, "metadata": { "id": "aiSx5e5SUsp7" }, "outputs": [], "source": [ "# Program CoinSimulationInteractive\n", "\n", "# By SuperElephantX\n", "# https://www.reddit.com/r/learnpython/comments/1ci995u/comment/l27y5ca/\n", "# Retrieved 17/VII/2025\n", "\n", "import random\n", "\n", "# Use \" \" to print verbatim\n", "print(\"Welcome to the Python coin flip\")\n", "\n", "while input(\"Press enter to flip the coin or type -1 to exit: \") != \"-1\":\n", " coin = random.choice([\"Heads\", \"Tails\"])\n", " print(coin)\n", "\n", "print(\"Goodbye!\")" ] }, { "cell_type": "markdown", "metadata": { "id": "ZAv_5TOBU1At" }, "source": [ "> Let us explain the code:\n", "\n", ">> The instruction `while input(\"Press enter to flip the coin or type -1 to exit: \") != \"-1\"`: commands to present an input box where the User must print `enter` else `-1`. If the input is `enter`, Python executes the indented code below the `while` instruction. If the input is `-1`, Python leaves the `while block` and exclaims `Goodbye!`. The equal sign `=` means equal to, but the expression `!=` means not equal to, so `!` is a negation. Hence, Python leaves the while loop when the User types something different from `-1`.\n", "\n", "\n", "> **Exercise.** Run the previous code and ask Colab to explain it. Run the code anew, but this time, type any letter. Explain the reaction of the program. What could happen if one inputs a random expression?\n", ">>\n", ">> **Answer.** If one types any single letter into the input box, the program flips a coin. This happens because the program only checks whether the input is different from `-1` to execute the indented code below the while instruction. In the contrary case, the program terminates.\n", "\n", "\n", "## Troubles with Colab?\n", "\n", "If Colab commands you to come later, this means that Google is inviting you to buy a Plus service. So, try some simple recommendations to enjoy the free version of Colab and get ready to discover much better advice:\n", "\n", "- Colab runs best in Chrome.\n", " \n", "- Work continuously, not intermittently.\n", "\n", "- Do your best to run a program only once. Modify parameters for a new run.\n", " \n", "- Use Colab preferentially to get explanations of the code.\n", "To run the code and to play with it by making changes at will,\n", "use a local solution, say,\n", "Visual Studio Code + JupyterLab.\n", "\n", "- Accept complaints and continue working.\n", "\n", "- Colab is very friendly with developers.\n", "It makes suggestions that almost always\n", "are opportune and correct.\n", "Give feedback on Colab when you have a clear evaluation.\n", "\n", "- The free version of Colab is a real help that must be used wisely.\n", "It is slow with heavy-duty computing.\n", "Use it to develop your code, but\n", "consider running it somewhere else.\n", "\n", "- Give a title to you work, in the corresponding box, and save your work with CTRL + s. \n", "It will be saved to Google Drive. \n", "To recover your work if you close it, \n", "type Colab as a prompt to Google and select Drive. \n", "Colab has kept it for you.\n", "\n", "\n", "## Encapsulation\n", "\n", "Blocks can be used various times in the same program. To avoid repeating and repeating, they are encapsulated and receive a name, and then they are referenced by name.\n", "\n", "> **Example.** The procedure for encapsulation follows the following model:" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "id": "50g7PlHAVw2V" }, "outputs": [], "source": [ "# Program encapsulation\n", "import random\n", "\n", "# Encapsulation of a flip.\n", "# It simulates a flip and report the outcome.\n", "# When the capsule is called somewhere else,\n", "# the returned value is the only information that is reported.\n", "def coin_flip():\n", " \"\"\"Simulates a single coin flip.\"\"\"\n", " return random.choice([\"Heads\", \"Tails\"])\n", "\n", "# Using the encapsulated code\n", "def main():\n", " # Encapsulated code is called\n", " # and result is assigned to the variable coin.\n", " coin = coin_flip()\n", " print(coin)\n", "\n", "if __name__ == \"__main__\":\n", " main()" ] }, { "cell_type": "markdown", "metadata": { "id": "t0jA5F_rVxzG" }, "source": [ "Encapsulation relies on functions. *Functions are processors of information that establish a correspondence between input and output. An encapsulated function is called a `method`*. They may receive or not receive an input that is placed inside the final pair of parentheses () of the definition. In this case, it is `coin_flip()` and there is no input. Each method has a name and contains a block to be executed, and the result may or may not be reported. To report a result, it must be preceded by the reserved word `return`. This word is reserved because it can be used with this one and only purpose. To invoke a method, one calls upon its name. If the function returns a value, it can be assigned to a variable, in this case, `coin`.\n", "\n", "> **Exercise.** Run the code and let Colab explain it. Is the last if-block unnecessary because one might replace it with `main()`?. Notice the double equal symbol and explain it.\n", ">>\n", ">> **Answer.** One can replace the last if-block by `main()` and the program runs and produces the expected result:\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "id": "S59mEUvRV4mK" }, "outputs": [], "source": [ "# Program encapsulation2\n", "import random\n", "\n", "def coin_flip():\n", " \"\"\"Simulates a single coin flip.\"\"\"\n", " return random.choice([\"Heads\", \"Tails\"])\n", "\n", "def main():\n", " coin = coin_flip()\n", " print(coin)\n", "\n", "main()" ] }, { "cell_type": "markdown", "metadata": { "id": "Mnl-rZhBV9wY" }, "source": [ "> The if-block" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "id": "makFAgBqWIlB" }, "outputs": [], "source": [ "if __name__ == \"__main__\":\n", " main()" ] }, { "cell_type": "markdown", "metadata": { "id": "q1nTipOLWL2v" }, "source": [ "> has, according to Colab, a purpose: ``This is a standard Python construct. It means that the code inside this if block will only run when the script is executed directly (not when it's imported as a module into another script).''\n", "This means that one can chain diverse encapsulated functions into a module that can be used as a toolbox to fulfill tasks in other programs. So, this if-block says that this script must be executed right now and that it is not a module to be imported somewhere else. The symbol == is the equality operator that is used for comparisons. By contrast, the symbol = is used for assignment of values to variables.\n", "\n", "> **Exercise.** Go to the webpage\n", ">\n", "> https://www.w3schools.com/python/default.asp\n", ">\n", "> which offers a free interactive Python tutorial. In the left column, there is a list of lessons. Select Python functions, and enjoy running in situ the proposed programs that come accompanied by clear and concise explanations. This site will be referred to as w3.\n", "\n", "\n", "# THE RELATIVE FREQUENCY\n", "\n", "*The relative frequency of heads is the total number of heads divided by the number of trials.* Therefore, the relative frequency can be calculated only at the end of the experiment. This task can be simulated in different ways. Let us choose one that primes the organization: we gather the results of each flip in a box, and when we are done, we count the number of heads and then calculate the relative frequency.\n", "\n", "Python has various types of boxes. A general-purpose box is the *Array* that is denoted by `[]`, which keeps an ordered list of objects that are separated by a comma. Objects can be numbers, letters, words, or sentences, or more complex entities.\n", "\n", "> **Example.** In the following program, results from flipping a coin are kept in the Array `results`:" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "id": "MFvPatO1WT7U" }, "outputs": [], "source": [ "# Program headsProportion.py\n", "import random\n", "\n", "num_flips = 7;\n", "\n", "def coin_flip():\n", " \"\"\"Simulates a single coin flip.\"\"\"\n", " return random.choice([\"Heads\", \"Tails\"])\n", "\n", "def simulate_flips(num_flips):\n", " \"\"\"Simulates multiple coin flips, \"\"\"\n", " \"\"\" prints the number of heads and their proportion.\"\"\"\n", " # Semantic respectful printing:\n", " # begin with f and suround variables with curly braces\n", " print(f\"Number of flips = {num_flips}\")\n", " results = [] # This is a void array\n", " for _ in range(num_flips):\n", " # append a new flip at the end of the array,\n", " results.append(coin_flip())\n", " print(f\"Results = {results}\")\n", " number_of_heads = results.count(\"Heads\")\n", " # Semantic respectful printing:\n", " # begin with f and suround variables with curly braces\n", " print(f\"head_count = {number_of_heads}\")\n", " proportion_of_heads = number_of_heads / num_flips\n", " # Print the proportion rounding to 2 decimal places\n", " print(f\"proportion of heads = + {proportion_of_heads:.2f}\")\n", "\n", "def main():\n", " simulate_flips(num_flips)\n", "\n", "if __name__ == \"__main__\":\n", " main()" ] }, { "cell_type": "markdown", "metadata": { "id": "ID5MyN4iWZfK" }, "source": [ "\n", "> **Exercise.** Run the program. The expression .2f is a format specifier. To test what it does, change it to .3f or .4f.\n", ">>\n", ">> **Answer.** The specifier .7f after a floating or real number says that it must be rounded to a precision of 7 decimal places.\n", "\n", "Python uses various types of numbers, distinguished by the memory cost to keep them. We have *integers*, like 3 and 7, and we have *reals* as 4.5 or -3.231. *Reals* are also called *float numbers*.\n", "\n", "# THE TENDENCY TOWARD A PROBABILITY\n", "\n", "Our purpose in this section is to test the idea that the relative frequency must observe a tendency of approaching a limit, and we need to know if that limit is the prescribed probability.\n", "\n", "> **Example.** We flip a coin a predefined number of times, and we calculate the relative frequency of heads. We increase the number of flips by one and find the relative frequency of heads. This is repeated for a certain number of times to catch a tendency, if any. The code, in test mode, follows:\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "id": "OxJ9xruWLbZU" }, "outputs": [], "source": [ "# Program headsProportionTendency.py\n", "import random\n", "\n", "num_experiments = 5\n", "\n", "def coin_flip():\n", " \"\"\"Simulates a single coin flip.\"\"\"\n", " return random.choice([\"Heads\", \"Tails\"])\n", "\n", "def simulate_flips(num_flips):\n", " \"\"\"Simulates multiple coin flips, \"\"\"\n", " \"\"\" prints the number of heads and their proportion.\"\"\"\n", " # Semantic respectful printing:\n", " # begin with f and suround variables with curly braces\n", " print(f\"Number of flips = {num_flips}\")\n", " results = [] #This is an array\n", " for _ in range(num_flips):\n", " # append a new flip at the end of the array,\n", " results.append(coin_flip())\n", " return results\n", "\n", "def calculate_tendency(num_experiments):\n", " \"\"\"The number of flips is increased one by one \"\"\"\n", " \"\"\"from 1 to num_experiments and \"\"\"\n", " \"\"\"the relative frequency of heads is calculated.\"\"\"\n", " num_flips = 1;\n", " for _ in range(num_experiments):\n", " results = simulate_flips(num_flips)\n", " print(f\"Results = {results}\")\n", " number_of_heads = results.count(\"Heads\")\n", " # Semantic respectful printing:\n", " # begin with f and suround variables with curly braces\n", " print(f\"head_count = {number_of_heads}\")\n", " proportion_of_heads = number_of_heads / num_flips\n", " # Print the proportion rounding to 2 decimal places\n", " print(f\"proportion of heads = + {proportion_of_heads:.2f}\")\n", " num_flips += 1\n", "\n", "def main():\n", " calculate_tendency(num_experiments)\n", "\n", "if __name__ == \"__main__\":\n", " main()" ] }, { "cell_type": "markdown", "metadata": { "id": "6gOtvHuoLTs2" }, "source": [ "> **Exercise.** Run the code and test whether it is correct.\n", ">>\n", ">> **Answer.** The code works as promised: it calculates the relative frequency of heads for an increasing number of flips from 1 to 5, and reports results.\n", "\n", "## Displaying results\n", "\n", "After seeing how the previous program works and deciding that it is correct, we can program a report of a list of successive values of the relative frequency when the number of trials grows. The nice news is that this list can be drawn in Python very easily. To learn how that can be done, we went to the W3 page, looked at the left column for the lesson on Graphics, discovered there a section entitled `Python Matplotlib`, looked at its material, and we found very soon what we needed: how to draw the listed values in an array.\n", "\n", "\n", "> **Example** The lesson offered by w3 was implemented as follows:" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "id": "eRs1V12jUgjf" }, "outputs": [], "source": [ "# Program relativeFrequencyTendency.py\n", "# Calculates the tendency of the relative frequency of heads\n", "# when a coin is flipped an increasing number of times.\n", "# Results are plotted.\n", "\n", "\n", "import random\n", "\n", "num_experiments = 5\n", "\n", "def coin_flip():\n", " \"\"\"Simulates a single coin flip.\"\"\"\n", " return random.choice([\"Heads\", \"Tails\"])\n", "\n", "def simulate_flips(num_flips):\n", " \"\"\"Simulates multiple coin flips, \"\"\"\n", " \"\"\" prints the number of heads and their proportion.\"\"\"\n", " # Semantic respectful printing:\n", " # begin with f and suround variables with curly braces\n", " print(f\"Number of flips = {num_flips}\")\n", " results = [] #This is an array\n", " for _ in range(num_flips):\n", " # append a new flip at the end of the array,\n", " results.append(coin_flip())\n", " return results\n", "\n", "def calculate_tendency(num_experiments):\n", " # Calculates the tendency of\n", " # the realtive frequency of heads\n", " # when the number of flips grows from 1 to num_experiments\n", " num_flips = 1;\n", " print(f\"Number of experiments = {num_experiments}\")\n", " print(f\"Number of flips per experiment = {num_flips}\")\n", " relative_frequencies = []\n", " for _ in range(num_experiments):\n", " results = simulate_flips(num_flips)\n", " print(f\"Results = {results}\")\n", " number_of_heads = results.count(\"Heads\")\n", " # Semantic respectful printing:\n", " # begin with f and suround variables with curly braces\n", " print(f\"head_count = {number_of_heads}\")\n", " proportion_of_heads = number_of_heads / num_flips\n", " relative_frequencies.append(proportion_of_heads)\n", " # Print the proportion rounding to 2 decimal places\n", " # The formatting expression \\n creates a void line to easy reading\n", " print(f\"proportion of heads = {proportion_of_heads:.2f} \\n\")\n", " num_flips += 1\n", " return relative_frequencies\n", "\n", "def draw_tendency(relative_frequencies):\n", " import matplotlib.pyplot as plt\n", " plt.plot(relative_frequencies)\n", "\n", "def main():\n", " relative_frequencies = calculate_tendency(num_experiments)\n", " draw_tendency(relative_frequencies)\n", "\n", "if __name__ == \"__main__\":\n", " main()" ] }, { "cell_type": "markdown", "metadata": { "id": "K48NpCZWUV-v" }, "source": [ "> **Exercise.** Run the previous program and verify that everything is sound.\n", ">>\n", ">> **Answer.** Everything is sound.\n", "\n", "> **Exercise.** Since we have every reason to believe that the program is correct, adapt the previous program to draw the results of a simulation with many experiments, say, 7000. Do you consider that the relative frequency observes a tendency, and to what extent?\n", ">>\n", ">> **Answer.** The instructions to print most messages in the previous program were silenced. We researched at w3 how to draw a red horizontal line representing the expected limit of the relative frequency of heads that must be 0.5. The code is the following:" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "id": "bmImw9A5Y9if" }, "outputs": [], "source": [ "# Program relativeFrequencyTendency2.py\n", "# Draws the tendency of the relative frequency of heads\n", "# when a coin is flipped an increasing number of times.\n", "\n", "import random\n", "\n", "num_experiments = 7000\n", "\n", "def coin_flip():\n", " \"\"\"Simulates a single coin flip.\"\"\"\n", " return random.choice([\"Heads\", \"Tails\"])\n", "\n", "def simulate_flips(num_flips):\n", " \"\"\"Simulates multiple coin flips, \"\"\"\n", " \"\"\" prints the number of heads and their proportion.\"\"\"\n", " # Semantic respectful printing:\n", " # begin with f and suround variables with curly braces\n", " # print(f\"Number of flips = {num_flips}\")\n", " results = [] #This is an array\n", " for _ in range(num_flips):\n", " # append a new flip at the end of the array,\n", " results.append(coin_flip())\n", " return results\n", "\n", "def calculate_tendency(num_experiments):\n", " \"\"\"Calculates the tendency of\"\"\"\n", " \"\"\"the realtive frequency of heads\"\"\"\n", " \"\"\"when the number of flips grows from 1 to num_experiments\"\"\"\n", " num_flips = 1;\n", " relative_frequencies = []\n", " for _ in range(num_experiments):\n", " results = simulate_flips(num_flips)\n", " # print(f\"Results = {results}\")\n", " number_of_heads = results.count(\"Heads\")\n", " # Semantic respectful printing:\n", " # begin with f and suround variables with curly braces\n", " # print(f\"head_count = {number_of_heads}\")\n", " proportion_of_heads = number_of_heads / num_flips\n", " relative_frequencies.append(proportion_of_heads)\n", " # Print the proportion rounding to 2 decimal places\n", " # The formatting expression \\n creates a void line to easy reading\n", " # print(f\"proportion of heads = + {proportion_of_heads:.2f} \\n\")\n", " num_flips += 1\n", " return relative_frequencies\n", "\n", "def draw_tendency(relative_frequencies):\n", " import matplotlib.pyplot as plt\n", " import numpy as np\n", " # The tendency of relative frequencies is drawn\n", " plt.plot(relative_frequencies)\n", " plt.xlabel(\"Number of flips\")\n", " plt.ylabel(\"Relative frequency\")\n", " plt.title(\"Tendency of the relative frequency\")\n", " plt.grid(True)\n", "\n", " # Red horizontal line representing the expected limit = 0.5\n", " xpoints = np.array([0, 7000]) # x-coordinates of all points\n", " ypoints = np.array([0.5, 0.5]) # y-coordinates of all points\n", " plt.plot(xpoints, ypoints, c = 'r')\n", " plt.show()\n", "\n", "def main():\n", " print(f\"Number of experiments = {num_experiments}\")\n", " print(\"Tendency of the relative frequency\")\n", " print(f\"when the number of flips grows from 1 to {num_experiments}\")\n", " relative_frequencies = calculate_tendency(num_experiments)\n", " draw_tendency(relative_frequencies)\n", "\n", "if __name__ == \"__main__\":\n", " main()" ] }, { "cell_type": "code", "source": [ "nding" ], "metadata": { "id": "WI5K1MJK8ZFC" }, "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "metadata": { "id": "lENz-Q_DY76G" }, "source": [ "> The graphic clearly shows that the uncertainty diminishes when the number of tries increases. Results tend to accumulate near the red line that marks the probability of heads as 0.5. Nevertheless, converging to this limit is very slow, and uncertainty seems to never end. In spite of this, we are very glad because the tendency to the expected value shows us that Python has a good generator of pseudo-random numbers.\n", "\n", "> **Conclusion.** We have every reason to believe that the relative frequency of heads has a tendency that is reaffirmed when the number of trials per round increases. It seems to very slowly converge to the probability. Anyway, let us identify the probability with the relative frequency when the number of trials is nearly infinite, say, 10000.\n", "\n", "\n", "> **Challenge.** Program `relativeFrequencyTendency2.py` is extremely demanding because it makes 7000 rounds, in which the number of flips for each round increases from 1 to 7000. Research or invent a way to make this work realizable at low computing effort, for instance, by calculating only a sample of values. Implement your idea.\n", "\n", "> **Research.** Resolve the following objection: the relative frequency does not converge to the probability. Instead, the relative frequency dances around the probability and very slowly gets tired but never dead.\n", "\n", "# THE TENDENCY TOWARDS A DISTRIBUTION\n", "\n", "We have observed that the relative frequency has the tendency to converge to a limit. When the model is correct and the source of pseudo-random numbers is acceptable, the limit matches the probability. If we accept this, we can go one step further.\n", "\n", "Let's imagine that we work with a constant number of coins or trials per round, say, 20. If we make an experiment, we might get 14 heads. But if we repeat the experiment, we will almost surely not get 14. Can we do Science under these circumstances? Yes, we can. The reason is that there is a tendency of the following nature:\n", "\n", "Let us count the number of times that we get 14 heads when we repeat the experiment several times, which could be 10000, each time with 20 coins. From that number, we can calculate the relative frequency of 14. But relative frequencies tend to a limit. Hence, we predict that for each k within 0 and 20, the relative frequency of getting k heads tends to a limit. Thus, we get a great prediction that reads: The overall resulting chart that depicts the tendency of each relative frequency of getting k heads must tend to a limit shape that is called the *Binomial Distribution.*\n", "\n", "Our next goal is to estimate that distribution. We will use a large but finite number of rounds, and so we hope that our estimation is accurate. To begin with, let us get acquainted with one astonishing possibility of modern technology.\n", "\n", "## Automatic code development\n", "\n", "While we have been developing our previous programs, we have been continuously assisted by Codey, one of the AIs behind Colab. It makes suggestions, and almost all are opportune and correct. In fact, today it is possible to ask AI services to develop code with specific requirements.\n", "\n", "> **Example.** We prompted Google with\n", ">\n", ">> chatgpt python code for a simulation with coins\n", ">\n", " >ChatGPT did not appear; instead, Gemini answered with the following code:" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "id": "TOzAqW6uWjmz" }, "outputs": [], "source": [ "# Program coinsByGemini.py\n", "\n", "import random\n", "\n", "def coin_flip():\n", " \"\"\"Simulates a single coin flip.\"\"\"\n", " return random.choice([\"Heads\", \"Tails\"])\n", "\n", "def simulate_flips(num_flips):\n", " \"\"\"Simulates multiple coin flips and returns a list of results.\"\"\"\n", " results = []\n", " for _ in range(num_flips):\n", " results.append(coin_flip())\n", " return results\n", "\n", "def count_outcomes(results):\n", " \"\"\"Counts the number of heads and tails in a list of results.\"\"\"\n", " heads_count = results.count(\"Heads\")\n", " tails_count = len(results) - heads_count\n", " return heads_count, tails_count\n", "\n", "def print_results(results):\n", " \"\"\"Prints the individual results of coin flips.\"\"\"\n", " print(\"Individual coin flip results:\")\n", " for result in results:\n", " print(result)\n", "\n", "def display_summary(heads_count, tails_count):\n", " \"\"\"Displays the summary of the simulation.\"\"\"\n", " print(\"\\nSimulation Summary:\")\n", " print(f\"Number of heads: {heads_count}\")\n", " print(f\"Number of tails: {tails_count}\")\n", " total_flips = heads_count + tails_count\n", " if total_flips > 0:\n", " print(f\"Percentage of heads: {heads_count / total_flips * 100:.2f}%\")\n", " print(f\"Percentage of tails: {tails_count / total_flips * 100:.2f}%\")\n", "\n", "def main():\n", " \"\"\"Main function to run the simulation.\"\"\"\n", " try:\n", " num_flips = int(input(\"Enter the number of coin flips: \"))\n", " if num_flips <= 0:\n", " print(\"Please enter a positive number of flips.\")\n", " return\n", "\n", " results = simulate_flips(num_flips)\n", " print_results(results)\n", " heads_count, tails_count = count_outcomes(results)\n", " display_summary(heads_count, tails_count)\n", "\n", " except ValueError:\n", " print(\"Invalid input. Please enter an integer.\")\n", "\n", "if __name__ == \"__main__\":\n", " main()\n" ] }, { "cell_type": "markdown", "metadata": { "id": "urvfbxUBWuVR" }, "source": [ "> Gemini added at the end: \"AI responses may include mistakes\". It also asked to assess its answer. This implies that we must run the code to see what it does.\n", "\n", "\n", "> **Exercise.** Run the previous code and ask Colab to explain it.\n", ">>\n", ">> **Answer.** The code runs spectacularly well, reporting the results of an experiment with the number of coins that is specified by the User in an input box. The output is structured, and each result is correct and is precisely identified.\n", "\n", "> **Exercise.** Run the previous code again but this time, type the word \"seven\" instead of the number 7. What happens and why?\n", ">>\n", ">> **Answer.** The program outputs: \"Invalid input. Please enter an integer\". This is the work of the block of code that begins with \"except ValueError:\". This is a style of programming that preaches that programs must accept any type of Users, specially beginners.\n", "\n", "> **Example.** We went to the web page of Gemini, and we asked it to develop the code for a simulation with coins, including some statistics. To that aim, we gave it the following prompt:\n", ">>\n", ">> Hey, Gemini! Please, help me with the Python code for a simulation with coins, repetitions, and statistics.\n", ">>\n", "> Gemini answered, asking for more details, and gave some ideas. So, we used them to be more specific:\n", ">>\n", ">> Number of coins = 10, number of repetitions = 100, statistics for the total number of heads across all repetitions, the proportion of heads, and the distribution of heads.\n", ">>\n", "> After two seconds, Gemini reported the following code:\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "id": "70nETlyTWzNj" }, "outputs": [], "source": [ "# Program headsCountByGemini.py\n", "# Simulates the binomial Distribution with\n", "# p = 0.5, n = num_coins_per_round.\n", "\n", "\n", "import random\n", "from collections import Counter\n", "\n", "# Parameters for the simulation\n", "num_coins_per_round = 10\n", "num_repetitions = 100\n", "\n", "# Simulate coin flips\n", "results = []\n", "for _ in range(num_repetitions):\n", " heads_in_round = 0\n", " for _ in range(num_coins_per_round):\n", " # 0 for tails, 1 for heads\n", " if random.choice([0, 1]) == 1:\n", " heads_in_round += 1\n", " results.append(heads_in_round)\n", "\n", "# Calculate statistics\n", "total_heads = sum(results)\n", "proportion_heads = total_heads / (num_coins_per_round * num_repetitions)\n", "heads_distribution = Counter(results)\n", "\n", "# Print the results\n", "print(f\"Simulation of {num_repetitions} rounds with {num_coins_per_round} coins per round:\")\n", "print(f\"\\nTotal number of heads across all repetitions: {total_heads}\")\n", "print(f\"Proportion of heads across all repetitions: {proportion_heads:.4f}\")\n", "print(\"Heads distribution\") #Added by the User\n", "print(heads_distribution) #Added by the User\n", "print(\"\\nDistribution of heads per round:\")\n", "# Sort the distribution for better readability\n", "sorted_distribution = dict(sorted(heads_distribution.items()))\n", "for num_heads, count in sorted_distribution.items():\n", " print(f\" {num_heads} heads: {count} times\")\n" ] }, { "cell_type": "markdown", "metadata": { "id": "9rtZdn3oW5gr" }, "source": [ "Must we worry about running the code? Not, absolutely! Gemini itself ran the code and reported the following results:\n", "\n", "\n", ">Simulation of 100 trials with 10 coins per trial:\n", ">\n", ">Total number of heads across all trials: 502\n", ">\n", ">Proportion of heads across all trials: 0.5020\n", ">\n", ">Distribution of heads per trial:\n", ">>\n", ">>1 heads: 2 times\n", ">>\n", ">>2 heads: 2 times\n", ">>\n", ">>3 heads: 11 times\n", ">>\n", ">>4 heads: 23 times\n", ">>\n", ">>5 heads: 21 times\n", ">>\n", ">>6 heads: 25 times\n", ">>\n", ">>7 heads: 12 times\n", ">>\n", ">>8 heads: 4 times\n", ">\n", "> Interestingly enough, Gemini did not add the warning that it can commit errors. Instead, it asked the User to assess its work.\n", "\n", "> **Exercise.** Summarize what the code does and assess it. \n", ">>\n", ">> **Answer.** This code simulates flipping 10 coins per round and making 100 rounds. Then it calculates the desired statistics: the total number of heads, the proportion of heads, and the distribution of how often a given number of heads happens. The output is structured, each item is identified, and the number of heads in the table adds up to 100. The code lacks the main()- formalism.\n", "\n", "\n", "> A portion of the code was not clear to us, so we copied it aside:\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "id": "ovTV-b6TQxGa" }, "outputs": [], "source": [ "# Sort the distribution for better readability\n", "sorted_distribution = dict(sorted(heads_distribution.items()))\n", "for num_heads, count in sorted_distribution.items():\n", " print(f\" {num_heads} heads: {count} times\")" ] }, { "cell_type": "markdown", "metadata": { "id": "4Ex2g5pvQ05U" }, "source": [ "> **Exercise.** Let Colab explain this *snippet.* What does this code do?\n", ">>\n", ">> **Answer.** The program `coinsStatisticsByGemini.py` repeats an experiment 100 times. Each experiment consists in flipping a coin 10 times and reporting the number of heads. Results range from 0 to 10. The program next counts the number of times that each specific number of heads happens. This produces an unordered table with pairs of numbers, saying, for instance, 5 heads happened 12 times. To order the resultant table, the `heads_distribution` is sorted by the number of heads in increasing order, and the result is located in a `dictionary` that is a table with pairs, whose first value is the number of heads and the second is the number of times that this number of heads happened.\n", "\n", "\n", "## Adding Graphs\n", "\n", "Statistics enter through the eyes. So, let us add graphics to the previous code.\n", "\n", "> **Example.** It is mandatory to draw the table of distribution of heads reported by the previous program. So, we asked Google: How to graph a table in Python? Google answered with an explanation and a piece of code, we observed it and decided to make a new try. We asked Google: how to graph in Python a math function? Google answered with a new variant of code that we preferred. The code was:\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "id": "SzcsW1qoI657" }, "outputs": [], "source": [ "# Program graphOfMathFunctionByGemini.py\n", "# how to graph in Python a math function?\n", "\n", "import matplotlib.pyplot as plt\n", "import numpy as np\n", "\n", "# 1. Define the mathematical function\n", "def my_function(x):\n", " return np.sin(x) + np.cos(x) # Example: a trigonometric function\n", "\n", "# 2. Generate x-values\n", "x_values = np.linspace(0, 2 * np.pi, 200) # 200 points from 0 to 2*pi\n", "\n", "# 3. Calculate corresponding y-values\n", "y_values = my_function(x_values)\n", "\n", "# 4. Plot the function\n", "plt.plot(x_values, y_values, label=\"sin(x) + cos(x)\")\n", "\n", "# 5. Add labels, title, and legend\n", "plt.xlabel(\"X-axis\")\n", "plt.ylabel(\"Y-axis\")\n", "plt.title(\"Graph of a Trigonometric Function\")\n", "plt.grid(True)\n", "plt.legend()\n", "\n", "# Display the plot\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": { "id": "LQtoSeDMMIKi" }, "source": [ "> **Exercise.** Run the code, and let Colab explain it.\n", ">>\n", ">> **Answer.** We ran the code and got admired. The basic idea to draw a function or a table is to put the input values in a vector and the output values in another vector, and use imported libraries to plot the graphic. Notice that input data is not given by an array of pairs of coordinates, but all X-coordinates go to one vector and all Y-coordinates to another. This means that the graph facility must build the pairs itself.\n", "\n", "\n", "\n", "> **Example.** Using what Gemini developed, we added code for a graph to the program `coinsStatisticsByGemini.py` as follows:" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "id": "uhl_7b2AOC48" }, "outputs": [], "source": [ "# Program headsCountByGemini2.py\n", "\n", "# Libraries and local rename\n", "import matplotlib.pyplot as plt\n", "import numpy as np #Support for scientific computing\n", "import random\n", "from collections import Counter\n", "\n", "# Parameters for the simulation\n", "num_coins_per_round = 10\n", "num_repetitions = 10000\n", "\n", "# Simulate coin flips\n", "results = []\n", "for _ in range(num_repetitions):\n", " heads_in_round = 0\n", " for _ in range(num_coins_per_round):\n", " # 0 for tails, 1 for heads\n", " if random.choice([0, 1]) == 1:\n", " heads_in_round += 1\n", " results.append(heads_in_round)\n", "\n", "# Calculate statistics\n", "total_heads = sum(results)\n", "proportion_heads = total_heads / (num_coins_per_round * num_repetitions)\n", "heads_distribution = Counter(results)\n", "\n", "x_values = []\n", "y_values = []\n", "\n", "# Print the results\n", "print(f\"Simulation of {num_repetitions} rounds \")\n", "print(f\"with {num_coins_per_round} coins per round:\")\n", "print(f\"\\nTotal number of heads across all repetitions: {total_heads}\")\n", "print(f\"Proportion of heads across all repetitions: {proportion_heads:.4f}\")\n", "print(\"\\nDistribution of heads per round:\")\n", "# Sort the distribution for better readability\n", "sorted_distribution = dict(sorted(heads_distribution.items()))\n", "for num_heads, count in sorted_distribution.items():\n", " x_values.append(num_heads)\n", " y_values.append(count)\n", " print(f\" {num_heads} heads: {count} times\")\n", "\n", "#Plot the function\n", "plt.plot(x_values, y_values, label=\"Heads_distribution\")\n", "\n", "# Add labels, title, and legend\n", "plt.xlabel(\"Number of heads\")\n", "plt.ylabel(\"Count\")\n", "plt.title(\"Graph of heads counts\")\n", "plt.grid(True)\n", "plt.legend(\"Heads count\")\n", "\n", "# Display the plot\n", "plt.show()\n", "\n", "\n", "\n" ] }, { "cell_type": "markdown", "metadata": { "id": "H1mToTT0Xlnc" }, "source": [ "> **Exercise.** Run the code with suitable modifications to see a tendency when the number of experiments increases.\n", ">>\n", ">> **Answer.** With 100 experiments, the data gathered itself around a mountain-like distribution. To test whether this is a tendency, we ran the code with 1000 experiments. The tendency was reaffirmed: an almost symmetric mountain appeared with a mode at five. This agrees with the calculated proportion of heads, that is, 0.493. Extreme values 1, 9, and 10 appeared, but not 0. Almost perfectly symmetric mountains appeared when the number of repetitions was set to 10000. In this case, all extreme values appeared. The observed proportion of heads was 0.4995. This is pretty close to 0.5, which is the expected proportion of heads when we programmed a random choice between 1 and 0, heads else tails. In conclusion, we observed two tendencies. First: the observed proportion of heads tends to 0.4995, and second, the silhouette of the mountain that describes the distribution of the numbers of heads tends to be more symmetric and with the form of a bell that fades away at the extreme values at any side of the bell, but not reaching `0`.\n", "\n", "> **Conclusion.** When we count the number of heads when we toss `n` coins, given that the probability of heads is `p`, the relative frequency of getting `k` heads tends to a limit. The overall graph with all `k` from 0 to n tends to a limit shape that is a bell.\n", "\n", "\n", "Our next goal is to estimate that distribution. The problem is that we seem to be frozen at `p = 0.5`. So, let us show how to get freedom.\n", "\n", "## Working with probabilities other than 0.5\n", "\n", "If we only know how to simulate a fair coin, we can twist the probability of heads with a simple trick.\n", "\n", "\n", "> **Example.** In the previous programs, we have identified heads with 1 and tails with 0, so a fair coin was simulated by making a random choice between heads and tails, else between 1 and 0. To twist the probability of heads, instead of making a random choice between 0 and 1, let us make a random choice between 0, 1, and 1. This is our way to assign to heads, represented by 1, a probability of 2/3 = 0.666 and to tails 1/3 = 0.333... The code follows:" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "id": "cO4rB-4pFlqY" }, "outputs": [], "source": [ "# Program coinsProbOfHeadsTwoThirds.py\n", "\n", "# Libraries and local rename\n", "import matplotlib.pyplot as plt\n", "import numpy as np #Support for scientific computing\n", "import random\n", "from collections import Counter\n", "\n", "# Parameters for the simulation\n", "num_coins_per_round = 20\n", "num_repetitions = 10000\n", "\n", "# Simulate coin flips\n", "results = []\n", "for _ in range(num_repetitions):\n", " heads_in_round = 0\n", " for _ in range(num_coins_per_round):\n", " # 0 for tails, 1 for heads\n", " if random.choice([0, 1, 1]) == 1:\n", " heads_in_round += 1\n", " results.append(heads_in_round)\n", "\n", "# Calculate statistics\n", "total_heads = sum(results)\n", "proportion_heads = total_heads / (num_coins_per_round * num_repetitions)\n", "heads_distribution = Counter(results)\n", "\n", "# Arrays to keep values to be graphed\n", "x_values = []\n", "y_values = []\n", "\n", "# Print the results\n", "print(\"Heads is 1 in a random choice among 0, 1, 1 \")\n", "print(\"So, the probability of heads is 2/3 = 0.666... \")\n", "print(\"and that of tails is 1/3 = 0.333...\")\n", "print(f\"Simulation of {num_repetitions} rounds \")\n", "print(f\"with {num_coins_per_round} coins per round:\")\n", "print(f\"\\nTotal number of heads across all repetitions: {total_heads}\")\n", "print(f\"Proportion of heads across all repetitions: {proportion_heads:.4f}\")\n", "print(\"\\nDistribution of heads per round:\")\n", "# Sort the distribution for better readability\n", "sorted_distribution = dict(sorted(heads_distribution.items()))\n", "for num_heads, count in sorted_distribution.items():\n", " x_values.append(num_heads)\n", " y_values.append(count)\n", " print(f\" {num_heads} heads: {count} times\")\n", "\n", "# Plot the function\n", "plt.plot(x_values, y_values, label=\"Heads_distribution\")\n", "\n", "# Add labels, title, and legend\n", "plt.xlabel(\"Number of heads\")\n", "plt.ylabel(\"Count\")\n", "plt.title(\"Graph of a heads distribution\")\n", "plt.grid(True)\n", "plt.legend()\n", "\n", "# Display the plot\n", "plt.show()\n" ] }, { "cell_type": "markdown", "metadata": { "id": "r-yeVOv7FCoY" }, "source": [ "> Since the probability of heads is 0.666..., we expect a bit more than 13 heads among 20 flips. The mode of the graphic fairly fits this expectation.\n", "\n", "> **Exercise.** Run the program varying the number of experiments to uncover a tendency.\n", ">>\n", ">> **Answer.** The tendency was an asymmetric mountain with one peak or mode at 7 for 10 coins and 13.5 for 20.\n", "\n", "> **Research.** Our method to run simulations with the probability of heads different from 0.5 is almighty but very boring. So, the question is: how are arbitrary probabilities simulated in Python?\n", "\n", "# THE BINOMIAL DISTRIBUTION\n", "\n", "With previous explorations, we got the impression that the silhouette of the distribution of heads for fair coins tends to fit into a bell that is symmetric and that tends to zero sideways. So, let us see what we can do to convince ourselves that this impression is right, and that it is true that the distribution of the number of heads indeed follows a tendency to a limit when the number of repetitions grows. That limit is called the ***Binomial Distribution.***\n", "\n", "> **Example.** Let us propose the following code to tackle the problem of the convergence of the distribution of k heads:\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "id": "wFIWnko6dmao" }, "outputs": [], "source": [ "# Program tendencyToBinomial.py\n", "\n", "# Libraries and local rename\n", "import matplotlib.pyplot as plt\n", "import numpy as np #Support for scientific computing\n", "import random\n", "from collections import Counter\n", "\n", "# Parameters for the simulation\n", "num_coins_per_round = 2\n", "flipping_numbers_to_be_sampled = [1, 10, 100, 1000, 10000]\n", "\n", "\n", "# Simulate coin flips\n", "# Returns an array with the number of heads per round\n", "def simulate_coin_flips(num_coins_per_round, num_repetitions):\n", " results = []\n", " for _ in range(num_repetitions):\n", " heads_in_round = 0\n", " for _ in range(num_coins_per_round):\n", " # 0 for tails, 1 for heads\n", " if random.choice([0, 1]) == 1:\n", " heads_in_round += 1\n", " results.append(heads_in_round)\n", " return results\n", "\n", "\n", "\n", "# Calculate statistics\n", "# Returns the sorted distribution of heads\n", "# as an array with num_heads and another with the corresponding count\n", "def calculate_statistics(results, num_coins_per_round, num_repetitions):\n", " total_heads = sum(results)\n", " proportion_heads = total_heads / (num_coins_per_round * num_repetitions)\n", " heads_distribution = Counter(results)\n", "\n", " # Arrays to keep values\n", " num_heads_values = []\n", " count_values = []\n", " # Print the results\n", " print(f\"Simulation of {num_repetitions} rounds \")\n", " print(f\"with {num_coins_per_round} coins per round:\")\n", " print(f\"\\nTotal number of heads across all repetitions: {total_heads}\")\n", " print(f\"Proportion of heads across all repetitions: {proportion_heads:.4f}\")\n", " print(\"\\nDistribution of heads:\")\n", " # Sort the distribution for better readability\n", " sorted_distribution = dict(sorted(heads_distribution.items()))\n", " for num_heads, count in sorted_distribution.items():\n", " relative_frequency = {count/num_repetitions}\n", " print(f\"{num_heads} heads: {count} times, {relative_frequency} relative frequency\")\n", " num_heads_values.append(num_heads)\n", " count_values.append(count)\n", " return num_heads_values, count_values\n", "\n", "\n", "# Find the observed distribution of the number of heads\n", "# for various flip values for a fix number of coins.\n", "def find_hyper_distribution(num_coins_per_round):\n", " for num_repetitions in flipping_numbers_to_be_sampled:\n", " print(\"\\n=============================\\n\")\n", " print(f\"Flipping number = {num_repetitions}\")\n", " results = simulate_coin_flips(num_coins_per_round, num_repetitions)\n", " num_heads_values, count_values = calculate_statistics(results, num_coins_per_round, num_repetitions)\n", "\n", "def main():\n", " find_hyper_distribution(num_coins_per_round)\n", "\n", "if __name__ == \"__main__\":\n", " main()\n" ] }, { "cell_type": "markdown", "metadata": { "id": "MtM7vodWbfVs" }, "source": [ "> **Exercise.** Run the program `tendencyToBinomial.py` and decide by inspection whether the distribution of the count of heads when normalized as relative frequencies approaches a limit. Compare your results with the mathematical expectancy for the case with 2 coins. Warning: the results of a program are displayed in a window of limited size. So, one must scroll the window to see all published results if the output is too long.\n", ">>\n", ">> **Answer.** The output shows that when the number of repetitions grows, the relative frequencies of 0, 1, and 2 heads when we toss two coins tend to a limit that coincides with the mathematical prediction given by the probability of heads. Let us recall that we are considering the case when heads and tails have equal probability.\n", ">>\n", ">> In fact, we have for two coins, when the probability of heads equals that of tails equals $1/2$, the following 4 possibilities hh, ht, th, tt where h and t denote heads and tails, respectively. Since both h and t have the same probability and events are independent, all 4 possibilities have the same probability, which must be $1/4$. If we count the number of heads, we see that 2 heads fall in the case hh with probability $1/4 = 0.25$, that 1 head happen in ht and th with total probability $1/4 + 1/4 = 1/2 = 0.5$, and we get 0 heads in the case tt with probability $1/4$. Ordering by the number of heads, we expect a table with relative frequencies 0.25, 0.5, 0.25. We see that the tables of relative frequencies of the simulation tend to this expectancy when the number of trials grows. For 10000 trials, the final table fits quite good this mathematical expectancy. So, we can say that the observed distribution of heads tends to a limit that coincides with the mathematical expectancy, the prediction established by the probabilities.\n", "\n", "Our results confirm that the observed distribution of heads tends to a limit that coincides with the mathematical expectancy. This was calculated assuming that heads and tails have the same probability and that events are independent, i.e., that the output of an experiment does not modify the relative frequency of future or parallel trials. In the official jargon, we say that the mathematical expectancy of the distribution of relative frequencies is the *Binomial Distribution* and that this is the limit to which observed distributions tend, assuming that we have an infinite number of trials and that the source of random numbers is perfect or at least acceptable.\n", "\n", "\n", "To better appreciate the tendency of the bells, let us graph them. So, we need to understand how we make 3D bar graphs in Python. Especially important for us is to know how data must be provided to the plotting procedures. It happens that to make 3D graphics, we need to work with 2D arrays.\n", "\n", "\n", "> **Example.** We asked Google about 2D arrays, and Gemini answered with the following code:\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "id": "HyGVy3vsp8By" }, "outputs": [], "source": [ "# Program 2D_array.py\n", "# By Gemini\n", "\n", "import numpy as np\n", "\n", "# Creating a 2D array using NumPy\n", "numpy_2d_array = np.array([[1, 2, 3],\n", " [4, 5, 6],\n", " [7, 8, 9]])\n", "\n", "# Accessing elements\n", "print(numpy_2d_array[0, 1]) # Accesses the element at row 0, column 1 (value: 2)\n", "\n", "# Performing operations (example: element-wise addition)\n", "array_b = np.array([[10, 11, 12],\n", " [13, 14, 15],\n", " [16, 17, 18]])\n", "result_array = numpy_2d_array + array_b\n", "print(result_array)" ] }, { "cell_type": "markdown", "metadata": { "id": "qpkQU0uqqTeh" }, "source": [ "\n", "> **Exercise.** Run the code and pay attention to the form as elements are accessed.\n", ">>\n", ">> **Answer.** With the `numpy` library, elements of a 2x2 matrix are accessed as in ordinary mathematics, but they are declared as vectors of vectors.\n", "\n", " With this, let us continue to 3D graphics.\n", "\n", "\n", "> **Example.** We prompted Google with\n", ">\n", "> examples of 3d bar graphics in Python\n", ">\n", "> Google corrected our prompt to \n", "> examples of 3d bar graphs in Python\n", ">\n", "> Indeed, graph refers to mathematical data, while graphics is more general. Gemini showed the following code and added suitable explanations.\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "id": "2LiterY-h8Cp" }, "outputs": [], "source": [ "# Program 3dBarChartByGemini.py\n", "\n", "import matplotlib.pyplot as plt\n", "import numpy as np\n", "from mpl_toolkits.mplot3d import Axes3D\n", "\n", "# Sample data\n", "x_data = np.arange(5)\n", "print(f\"x_data = {x_data} \\n\")\n", "y_data = np.arange(5)\n", "print(f\"y_data = {y_data} \\n\")\n", "# Random heights for the bars\n", "# Real values with max = 10\n", "# It is a 5X5 matrix\n", "z_heights = np.random.rand(5, 5) * 10\n", "print(f\"z_heights = {z_heights} \\n\")\n", "\n", "# Create a meshgrid for the x and y coordinates of the bars\n", "_xx, _yy = np.meshgrid(x_data, y_data)\n", "print(f\"_xx = {_xx} \\n\")\n", "print(f\"_yy = {_yy} \\n\")\n", "x_pos, y_pos = _xx.ravel(), _yy.ravel()\n", "print(f\"x_pos = {x_pos} \\n\")\n", "print(f\"y_pos = {y_pos} \\n\")\n", "z_pos = np.zeros_like(x_pos) # All bars start at z=0\n", "print(f\"z_pos = {z_pos} \\n\")\n", "\n", "# Define the dimensions of the bars\n", "dx = dy = 0.8 # Width and depth of the bars\n", "dz = z_heights.ravel() # Heights of the bars\n", "print(f\"dz = {dz} \\n\")\n", "\n", "# Create the 3D plot\n", "fig = plt.figure()\n", "ax = fig.add_subplot(111, projection='3d')\n", "\n", "# Plot the 3D bars\n", "ax.bar3d(x_pos, y_pos, z_pos, dx, dy, dz, shade=True)\n", "\n", "# Set labels and title\n", "ax.set_xlabel('X-axis')\n", "ax.set_ylabel('Y-axis')\n", "ax.set_zlabel('Z-axis (Height)')\n", "ax.set_title('Basic 3D Bar Plot')\n", "\n", "# Display the plot\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": { "id": "t-E4hC4hjw9z" }, "source": [ "> **Exercise.** Run the code and pay attention to the form of encoding data for the graphic procedures.\n", ">>\n", ">> **Answer.** The code is complex. We expect to suffer a lot to adapt it to our purpose. In the meantime, we see that the code invents data and then uses it as raw data for the graphs. The form as data must be packed to be understandable by the drawing function is an extension of that of 2D-graphs: the points where bars begin are not encoded as a triple. Instead, all x-coordinates of all points go into the vector, or one dimensional array, `x_pos`, all y-coordinates go into `y_pos`, all z-coordinates into `z_pos`. Besides, each bar has a height that must be declared as `dz`. The width `dx` and depth `dy` of the bars are constants.\n", "\n", "> **Example.** By following the example offered by Gemini, we propose the following program that graphs the tendency of the distribution of the number of heads towards the Binomial Distribution:\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "id": "ELmwAfFQjxxQ" }, "outputs": [], "source": [ "# Program tendencyToBinomialGraph.py\n", "\n", "# Libraries and local rename\n", "import matplotlib.pyplot as plt\n", "import numpy as np # Support for scientific computing\n", "import random\n", "from collections import Counter\n", "\n", "\n", "# Parameters for the simulation\n", "num_coins_per_round = 5\n", "flipping_numbers_to_be_sampled = [10, 100, 1000, 10000]\n", "\n", "\n", "# Data for the graph\n", "# Sample data\n", "x_data = np.arange(num_coins_per_round +1)\n", "print(f\"x_data = {x_data} \\n\")\n", "numb_flipping_values = len(flipping_numbers_to_be_sampled)\n", "print(f\"numb_flipping_values = {numb_flipping_values} \\n\")\n", "y_data = np.arange(numb_flipping_values)\n", "print(f\"y_data = {y_data} \\n\")\n", "\n", "\n", "\n", "# Create a meshgrid for the x and y coordinates of the bars\n", "_xx, _yy = np.meshgrid(x_data, y_data)\n", "x_pos, y_pos = _xx.ravel(), _yy.ravel()\n", "print(f\"x_pos = {x_pos} \\n\")\n", "print(f\"y_pos = {y_pos} \\n\")\n", "z_pos = np.zeros_like(x_pos) # All bars start at z=0\n", "z_heights = 1.0 *np.zeros_like(_xx) #z_heights initialization\n", "print(f\"_xx = {_xx} \\n\")\n", "print(f\"_yy = {_yy} \\n\")\n", "print(f\"z_heights = {z_heights} \\n\")\n", "\n", "\n", "# Define the dimensions of the bars\n", "dx = dy = 0.6 # Width and depth of the bars\n", "\n", "# Simulate coin flips\n", "# Returns an array with the number of heads per round\n", "def simulate_coin_flips(num_coins_per_round, num_repetitions):\n", " results = []\n", " for _ in range(num_repetitions):\n", " heads_in_round = 0\n", " for _ in range(num_coins_per_round):\n", " # 0 for tails, 1 for heads\n", " if random.choice([0, 1]) == 1:\n", " heads_in_round += 1\n", " results.append(heads_in_round)\n", " return results\n", "\n", "\n", "# Calculate statistics\n", "# Returns the sorted distribution of heads\n", "# as an array with num_heads and another with the corresponding count\n", "def calculate_statistics(results, num_coins_per_round, index_num_repetitions):\n", " print(\"\\n calculate_statistics\")\n", " total_heads = sum(results)\n", " num_repetitions = flipping_numbers_to_be_sampled[index_num_repetitions]\n", " proportion_heads = total_heads / (num_coins_per_round * num_repetitions)\n", " heads_distribution = Counter(results)\n", "\n", " # Arrays to keep values\n", " num_heads_values = []\n", " count_values = []\n", " relativeFreq_values = []\n", "\n", " # Print the results\n", " print(f\"Simulation of {num_repetitions} rounds \")\n", " print(f\"with {num_coins_per_round} coins per round:\")\n", " print(f\"\\nTotal number of heads across all repetitions: {total_heads}\")\n", " print(f\"Proportion of heads across all repetitions: {proportion_heads:.4f}\")\n", " print(\"\\nDistribution of heads:\")\n", " # Sort the distribution for better readability\n", " sorted_distribution = dict(sorted(heads_distribution.items()))\n", " for num_heads, count in sorted_distribution.items():\n", " relative_frequency = count/num_repetitions\n", " print(f\"{num_heads} heads: {count} times, {relative_frequency} relative frequency\")\n", " # z_heights[num_heads, index_num_repetitions] = relative_frequency\n", " z_heights[index_num_repetitions, num_heads] = relative_frequency\n", " print(f\"num_heads = {num_heads} index_num_repetitions = {index_num_repetitions}\")\n", " print(f\"z_heights[{index_num_repetitions},{num_heads}] = {z_heights[index_num_repetitions, num_heads]}\")\n", " print(f\"z_heights {z_heights}\")\n", " return z_heights\n", "\n", "\n", "# Find and plot the observed distribution of the number of heads\n", "# for various flip values for a fix number of coins.\n", "def plot_hyper_distribution(num_coins_per_round):\n", " print(\"plot_hyper_distribution(num_coins_per_round):\")\n", " print(f\"num_coins_per_round = {num_coins_per_round}\")\n", " print(f\"flipping_numbers_to_be_sampled = {flipping_numbers_to_be_sampled}\")\n", " print(f\"len(flipping_numbers_to_be_sampled) = {len(flipping_numbers_to_be_sampled)}\")\n", "\n", "for index_num_repetitions in range(len(flipping_numbers_to_be_sampled)):\n", " print(\"\\n=============================\\n\")\n", " print(f\"index_num_repetitions = {index_num_repetitions}\")\n", " num_repetitions = flipping_numbers_to_be_sampled[index_num_repetitions]\n", " print(f\"Flipping number = {num_repetitions}\")\n", " results = simulate_coin_flips(num_coins_per_round, num_repetitions)\n", " z_heights = calculate_statistics(results, num_coins_per_round, index_num_repetitions)\n", " # print(f\"z_heights = {z_heights }\" )\n", " index_num_repetitions += 1\n", "\n", "# Create a meshgrid for the x and y coordinates of the bars\n", "_xx, _yy = np.meshgrid(x_data, y_data)\n", "print(f\"_xx = {_xx} \\n\")\n", "print(f\"_yy = {_yy} \\n\")\n", "x_pos, y_pos = _xx.ravel(), _yy.ravel()\n", "print(f\"x_pos = {x_pos} \\n\")\n", "print(f\"y_pos = {y_pos} \\n\")\n", "z_pos = np.zeros_like(x_pos) # All bars start at z=0\n", "print(f\"z_pos = {z_pos} \\n\")\n", "dz = z_heights.ravel() # Heights of the bars\n", "print(f\"dz = {dz} \\n\\n\")\n", "\n", "print(\"Is there a tendency of the sampled distributions\")\n", "print(\"when the number of repetitions increases?\")\n", "for i in range(4):\n", " print(f\"Number of repetitions {flipping_numbers_to_be_sampled[i]}: sampled distribution{z_heights[i]}\")\n", "print(\"Mathematically expected distribution: [0.25, 0.5, 0.25]\")\n", "\n", "# Create the 3D plot\n", "fig = plt.figure()\n", "ax = fig.add_subplot(111, projection='3d')\n", "\n", "# Plot the 3D bars\n", "ax.bar3d( y_pos, x_pos, z_pos, dx, dy, dz, shade=True)\n", "\n", "# Set labels and title\n", "ax.set_xlabel('num_repetitions = 10^x')\n", "ax.set_ylabel('num_heads')\n", "ax.set_zlabel('relativeFreq')\n", "ax.set_title('Tendency towards the Binomial Distribution')\n", "\n", "# Display the plot\n", "plt.show()\n", "\n", "\n", "\n", "def main():\n", " plot_hyper_distribution(num_coins_per_round)\n", "\n", "if __name__ == \"__main__\":\n", " main()\n", "\n" ] }, { "cell_type": "markdown", "metadata": { "id": "SoIEq2Gfhgw_" }, "source": [ "> **Exercise.** Run the previous program. Is there a stable tendency of the sampled distributions to a limit? Does the limit coincide with that mathematically expected?\n", ">>\n", ">> **Answer.** A stable tendency to a limit is visible, and the limit coincides with the mathematically expected distribution.\n", "\n", "> **Exercise.** The previous program works perfectly for the case in which we toss 2 coins. Generalize it to cover any number of coins.\n", ">>\n", ">> **Answer** We tested the code for 6 coins. It worked perfectly, but it lacked the mathematically expected distribution. So we open a code cell in this document. Colab invited us to ask for a code, so we clicked on the word `generate` and prompted:\n", ">>\n", ">>> the general binomial distribution.\n", "\n", ">> In response, Gemini presented the following code:\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "id": "S0tzQs371T3G" }, "outputs": [], "source": [ "# Program tendencyToBinomialGraphGeneral.py\n", "# Code developed by Gemini inside Colab.\n", "# It reuses and improves our code,\n", "# specially in regard with the graphs.\n", "\n", "# Libraries and local rename\n", "import matplotlib.pyplot as plt\n", "import numpy as np # Support for scientific computing\n", "import random\n", "from collections import Counter\n", "from mpl_toolkits.mplot3d import Axes3D\n", "\n", "# Parameters for the simulation\n", "flipping_numbers_to_be_sampled = [10, 100, 1000, 10000]\n", "\n", "\n", "# Simulate coin flips\n", "# Returns an array with the number of heads per round\n", "def simulate_coin_flips(num_coins_per_round, num_repetitions):\n", " results = []\n", " for _ in range(num_repetitions):\n", " heads_in_round = 0\n", " for _ in range(num_coins_per_round):\n", " # 0 for tails, 1 for heads\n", " if random.choice([0, 1]) == 1:\n", " heads_in_round += 1\n", " results.append(heads_in_round)\n", " return results\n", "\n", "\n", "# Calculate statistics\n", "# Returns the sorted distribution of heads\n", "# as an array with num_heads and another with the corresponding relative frequency\n", "def calculate_statistics(results, num_coins_per_round, index_num_repetitions, num_repetitions, z_heights):\n", " print(\"\\n calculate_statistics\")\n", " total_heads = sum(results)\n", " proportion_heads = total_heads / (num_coins_per_round * num_repetitions)\n", " heads_distribution = Counter(results)\n", "\n", " # Print the results\n", " print(f\"Simulation of {num_repetitions} rounds \")\n", " print(f\"with {num_coins_per_round} coins per round:\")\n", " print(f\"\\nTotal number of heads across all repetitions: {total_heads}\")\n", " print(f\"Proportion of heads across all repetitions: {proportion_heads:.4f}\")\n", " print(\"\\nDistribution of heads:\")\n", " # Sort the distribution for better readability\n", " sorted_distribution = dict(sorted(heads_distribution.items()))\n", " for num_heads, count in sorted_distribution.items():\n", " relative_frequency = count/num_repetitions\n", " print(f\"{num_heads} heads: {count} times, {relative_frequency} relative frequency\")\n", " if num_heads <= num_coins_per_round and index_num_repetitions < len(flipping_numbers_to_be_sampled):\n", " z_heights[index_num_repetitions, num_heads] = relative_frequency\n", " print(f\"z_heights[{index_num_repetitions},{num_heads}] = {z_heights[index_num_repetitions, num_heads]}\")\n", "\n", " print(f\"z_heights {z_heights}\")\n", " return z_heights\n", "\n", "\n", "# Find and plot the observed distribution of the number of heads\n", "# for various flip values for a fix number of coins.\n", "def plot_hyper_distribution(num_coins_per_round):\n", " print(\"plot_hyper_distribution(num_coins_per_round):\")\n", " print(f\"num_coins_per_round = {num_coins_per_round}\")\n", " print(f\"flipping_numbers_to_be_sampled = {flipping_numbers_to_be_sampled}\")\n", " numb_flipping_values = len(flipping_numbers_to_be_sampled)\n", " print(f\"len(flipping_numbers_to_be_sampled) = {numb_flipping_values}\")\n", "\n", " # Data for the graph\n", " x_data = np.arange(num_coins_per_round + 1)\n", " print(f\"x_data = {x_data} \\n\")\n", " y_data = np.arange(numb_flipping_values)\n", " print(f\"y_data = {y_data} \\n\")\n", "\n", " # Create a meshgrid for the x and y coordinates of the bars\n", " _xx, _yy = np.meshgrid(x_data, y_data)\n", " print(f\"_xx = {_xx} \\n\")\n", " print(f\"_yy = {_yy} \\n\")\n", " x_pos, y_pos = _xx.ravel(), _yy.ravel()\n", " print(f\"x_pos = {x_pos} \\n\")\n", " print(f\"y_pos = {y_pos} \\n\")\n", " z_pos = np.zeros_like(x_pos) # All bars start at z=0\n", "\n", " z_heights = 1.0 * np.zeros((numb_flipping_values, num_coins_per_round + 1)) # z_heights initialization\n", " print(f\"z_heights initialized = {z_heights} \\n\")\n", "\n", "\n", " for index_num_repetitions in range(numb_flipping_values):\n", " print(\"\\n=============================\\n\")\n", " print(f\"index_num_repetitions = {index_num_repetitions}\")\n", " num_repetitions = flipping_numbers_to_be_sampled[index_num_repetitions]\n", " print(f\"Flipping number = {num_repetitions}\")\n", " results = simulate_coin_flips(num_coins_per_round, num_repetitions)\n", " z_heights = calculate_statistics(results, num_coins_per_round, index_num_repetitions, num_repetitions, z_heights)\n", "\n", "\n", " dz = z_heights.ravel() # Heights of the bars\n", " print(f\"dz = {dz} \\n\\n\")\n", "\n", " print(\"Is there a tendency of the sampled distributions\")\n", " print(\"when the number of repetitions increases?\")\n", " for i in range(numb_flipping_values):\n", " print(f\"Number of repetitions {flipping_numbers_to_be_sampled[i]}: sampled distribution{z_heights[i]}\")\n", "\n", " # Create the 3D plot\n", " fig = plt.figure(figsize=(10, 8))\n", " ax = fig.add_subplot(111, projection='3d')\n", "\n", " # Define the dimensions of the bars\n", " dx = dy = 0.6 # Width and depth of the bars\n", "\n", " # Plot the 3D bars\n", " ax.bar3d( y_pos, x_pos, z_pos, dx, dy, dz, shade=True)\n", "\n", " # Set labels and title\n", " ax.set_xlabel('Number of Repetitions (10^x)')\n", " ax.set_ylabel('Number of heads')\n", " ax.set_zlabel('Relative Frequency')\n", " ax.set_title('Tendency towards the Binomial Distribution')\n", "\n", " # Set the tick labels for the y-axis (number of repetitions)\n", " ax.set_yticklabels([f'10^{i+1}' for i in range(numb_flipping_values)])\n", "\n", " # Display the plot\n", " plt.show()\n", "\n", "\n", "# Example usage:\n", "# plot_hyper_distribution(num_coins_per_round=5)\n", "# plot_hyper_distribution(num_coins_per_round=10)\n", "plot_hyper_distribution(num_coins_per_round=20) # Example with 20 coins per round" ] }, { "cell_type": "markdown", "metadata": { "id": "9UwPCwuG1LqI" }, "source": [ "> **Exercise.** Run the code by Gemini and evaluate it.\n", ">>\n", ">> **Answer.** The code runs spectacularly well for 20 coins, although the Number of heads in the graphic is incorrectly numbered. This code improves that of the Author.\n", "\n", "\n", "We have observational data, but we lack the mathematically expected distribution. This void hinders us from relating theory to observation. We must fill in this gap to do something with a scientific flavor.\n", "\n", "# THE EXPECTED BINOMIAL DISTRIBUTION\n", "\n", "Science compares what one observes with what one thinks. Here, observation corresponds to the output of our simulations, while what we think is derived from mathematical models.\n", "\n", "Let us show now the main idea of the mathematical model that predicts the shape of the Binomial Distribution for all possible values of $n$, the number of coins. We follow an induction method, beginning with simple observation and inferring what could be more general.\n", "\n", "To proceed further, we need two important rules of probability theory.\n", "\n", "- The first rule is\n", "\n", "> *For independent events, $A$ and $B$, the probability of their intersection $p(A \\cap B)$, is the product of probabilities*:\n", ">>\n", ">> $p(A \\cap B) = p(A)p(B)$\n", ">>\n", ">\n", "> **Example.** If the probability of heads is $1/2$ and that of tails is $1/2$, which is the probability of getting two heads if we throw 2 coins? Because we are assuming that no event influences any other, the probability of heads in the first coin and heads in the second is $(1/2) \\times (1/2) = 1/4$. This means one option in four when all are equiprobable. In fact, and denoting head as h, and tail as t, for two coins we get hh, ht, th, tt and all four options have the same probability because of the independence and the equiprobability of h and t. Hence, the probability of each one is $1/4$.\n", ">\n", "> This rule of the multiplication of probabilities for independent events is too abstract to be accepted by all persons. Let's derive it from a stand that might be easier to accept.\n", ">\n", "> Let us define that the event $A$ is independent from the event $B$ when $B$ plays no role over $A$. More specifically, when the probability of $A$ given $B$ is equal to the probability of $A$. But, how shall we define the probability of $A$ given $B$? We are not considering cause-effect, but only counting procedures. In this vein, $B$ has an influence over $A$ if they have something in common, and the influence is greater the larger is that intersection with respect to $B$. That is why we define the probability of $A$ given $B$, denoted $P(A/B)$, as\n", ">>\n", ">> $P(A/B) = \\frac{p(A \\cap B )}{p(B)}$\n", ">>\n", "> where $A \\cap B $ is the intersection, i.e., the set with all common events.\n", "Thus, the event $A$ is independent of the event $B$ when\n", ">>\n", ">> $P(A/B) = \\frac{p(A \\cap B )}{p(B)} = P(A)$\n", ">>\n", "> Therefore, the event $A$ is independent from the event $B$ when\n", ">>\n", ">> $ p(A \\cap B ) = P(A)p(B) $\n", ">>\n", "> By the same token, the event $B$ is independent from the event $A$ when\n", ">>\n", ">> $ p(B \\cap A ) = P(B)p(A) $\n", ">>\n", "> We say that **events $A$ and $B$ are independent** when neither depends on the other. That is, when\n", ">>\n", ">> $ p(A \\cap B ) = p(A)p(B) $\n", ">>\n", "\n", "- The second rule is:\n", "\n", ">> *When two events $A$ and $B$ are disjoint, the probability of their union, $p(A \\cup B)$, is the sum of probabilities,* i.e., when\n", ">>>\n", ">>> $p(A \\cup B) = p(A) + p(B)$\n", ">>\n", "> where the union lists all events from $A$ and adds those from $B$ that appear for the first time. Thus, if $A$ and $B$ have nothing in common, the probability of their union counts all numbers from both sets without any subtraction because there is no double counting.\n", ">\n", "> **Example.** If we throw two coins, what is the probability of getting 1 head? To obtain 1 head, we must get ht or th. More explicitly, the first coin could be h and the second t, with probability $1/4$, or the first could be t and the other h, with probability $1/4$. These two options have nothing in common. They are disjoint. Thus, the probability of 1 head is the sum of these two probabilities: $1/4 + 1/4 = 1/2$ that means that the half of the events shall contain 1 head. This is correct: all events are tt, ht, th, tt and the half of them have 1 head.\n", "\n", "\n", "\n", "We will apply these rules everywhere in the sequel. The first rule is that of the multiplication of probabilities for independent events, and the second is for the addition of probabilities for disjoint events.\n", "\n", "\n", "How shall we proceed when two events are not disjoint?\n", "\n", "> **Exercise.** If we throw two coins, what is the probability of getting exactly 1 head or 1 tail?\n", ">>\n", ">> **Answer.** If we throw two coins, getting exactly 1 head can happen if we get ht or th, and getting exactly 1 tail can happen in exactly the same situation. Therefore, the probability of getting exactly 1 head or 1 tail is $1/2$.\n", "\n", "With these two rules, we can make suitable predictions.\n", "\n", "We work with a generalized coin whose probability of heads is $p$ and that of tails is $q = 1-p$. Heads is denoted h and tails t. The likelihood of t plus that of h is $q + p = p + q = 1$. The task is to predict the shape of the Distribution of heads, that is, the expected probability of getting $k$ heads when we throw $n$ coins. Thus, we have a function that assigns a probability to each $k$ in between $0$ and $n$. This function receives a name: it is a **random variable** and is denoted $X$. In short, $X$ counts the number of heads and assigns a probability.\n", "\n", "Let us work out some concrete predictions about the shape of the Binomial Distribution.\n", "\n", "\n", "> **Example.** If we throw one coin, we can get t or h. In the case we get t, $X$ is $0$ with probability $q$. If we get h, $X$ is $1$ with probability $p$. Therefore, the random variable $X$, that counts heads, is the function defined by:\n", ">>\n", ">> $t \\rightarrow X = 0 \\rightarrow q$\n", ">>\n", ">> $h \\rightarrow X = 1 \\rightarrow p$\n", ">>\n", ">> All probabilities must add up to $1$. Is that true? Yes:\n", ">>\n", ">> $p + q = p + (1-p) = p +1 -p = 1$\n", "\n", "\n", "> **Example.** If we throw two coins, the output could be tt, th, ht, hh. The corresponding values of $X$ are $0$, $1$, $1$, $2$ with corresponding probabilities $qq = q^2$, $qp$, $pq$, $pp = p^2$. Therefore, the random variable $X$, which counts heads, is the function defined by\n", ">>\n", ">> $tt \\rightarrow X = 0 \\rightarrow q^2$ \n", ">> $th, ht \\rightarrow X = 1 \\rightarrow 2qp = 2 pq$ \n", ">> $hh \\rightarrow X = 2 \\rightarrow p^2$ \n", ">> \n", ">> All probabilities must add up to $1$. Is that true? Yes: \n", ">> $S = q^2 + 2pq + p^2 = (p + q)^2 = 1^2 = 1 $\n", "\n", "\n", "> **Inference.** The shape of the Binomial Distribution for $n = 1$ and $n = 2$ coins is given by $(q+p)^n$. Our inference is that this is valid for every number of coins $n$.\n", "\n", "> **Example.** Is our inference true for $n=3$? Yes, it is. To see that, we need to list all options, but without repetitions. To that aim, let us recur to the binary numbers, or numbers in base $2$. In that notation, $0$ means $0$ heads or t, and $1$ means $1$ head or h:\n", ">>\n", ">> $0 \\rightarrow 000 \\rightarrow ttt \\rightarrow X = 0 \\rightarrow q^3$ \n", ">> $1 \\rightarrow 001 \\rightarrow tth \\rightarrow X = 1 \\rightarrow q^2 p$ \n", ">> $2 \\rightarrow 010 \\rightarrow tht \\rightarrow X = 1 \\rightarrow q^2 p$ \n", ">> $3 \\rightarrow 011 \\rightarrow thh \\rightarrow X = 2 \\rightarrow q p^2$ \n", ">> $4 \\rightarrow 100 \\rightarrow htt \\rightarrow X = 1 \\rightarrow q^2 p$ \n", ">> $5 \\rightarrow 101 \\rightarrow hth \\rightarrow X = 2 \\rightarrow q p^2$ \n", ">> $6 \\rightarrow 110 \\rightarrow hht \\rightarrow X = 2 \\rightarrow q p^2$ \n", ">> $7 \\rightarrow 111 \\rightarrow hhh \\rightarrow X = 3 \\rightarrow p^3$ \n", ">\n", "> Test: all probabilities must add to $1$: Is that true? Let us see:\n", ">>\n", ">> $S = q^3 + q^2 p + q^2 p + q p^2 + q^2 p + q p^2 + q p^2 + p 3$\n", ">>\n", "> Reordering to place together similar terms, we get\n", ">>\n", ">> $ S = q^3 + q^2 p + q^2 p + q^2 p + q p^2 + q p^2 + q p^2 + p 3$\n", ">\n", "> Associating:\n", ">>\n", ">> $ S = q^3 + 3q^2 p + 3q p^2 + p 3 = (q + p)^3 = (p+q)^3 = 1 ^3 = 1$\n", ">\n", "> We passed the test.\n", ">\n", "> *Prediction 1.* The shape of the Binomial Distribution for $n = 3$ is defined by the random variable $X$ given by:\n", ">>\n", ">> $0 \\rightarrow 1q^3 $ \n", ">> $1 \\rightarrow 3q^2p$ \n", ">> $2 \\rightarrow 3q p^2$ \n", ">> $3 \\rightarrow 1p^3 $ \n", ">\n", "> *Prediction 2.* If $p = q = 1/2$, the corresponding bell is symmetric and the probabilities are\n", ">>\n", ">> $0 \\rightarrow q^3 = 1 (1/2)^3 = 1/8$ \n", ">> $1 \\rightarrow 3q^2p = 3 (1/2)^2 (1/2) = 3/8$ \n", ">> $2 \\rightarrow 3q p^2 3 (1/2) (1/2)^2 = 3/8 $ \n", ">> $3 \\rightarrow 1p^3 = (1/2)^3 = 1/8$ \n", ">>\n", ">> So, the heights of the bars are in relation 1:3:3:1.\n", "\n", "The generality of these counting procedures is called ***Combinatorics.***\n", "\n", "\n", "> **Exercise.** Use the trick of binary numbers to predict the shape of the Binomial distribution for $n = 4$. The binomial expansion for $= 4$ is\n", ">>\n", ">>$(q + p)^4 = q^4 + 4 q^3 p + 6 q^2 p^2 + 4 qp^3 + p^4$\n", "\n", ">> **Answer.** All cases with $n= 4$ are:\n", ">>>\n", "\n", ">>> $0 \\rightarrow 0000 \\rightarrow tttt \\rightarrow X = 0 \\rightarrow q^4$ \n", ">>> $1 \\rightarrow 0001 \\rightarrow ttth \\rightarrow X = 1 \\rightarrow q^3p$ \n", ">>> $2\\rightarrow 0010 \\rightarrow ttht \\rightarrow X = 1 \\rightarrow q^3p $ \n", ">>> $3 \\rightarrow 0011 \\rightarrow tthh \\rightarrow X = 2 \\rightarrow q^2 p^2 $ \n", ">>> $4 \\rightarrow 0100 \\rightarrow thtt \\rightarrow X = 1 \\rightarrow q^3 p$ \n", ">>> $5 \\rightarrow 0101 \\rightarrow thth \\rightarrow X = 2 \\rightarrow q^2 p^2$ \n", ">>> $6 \\rightarrow 0110 \\rightarrow thht \\rightarrow X = 2 \\rightarrow q^2 p^2$ \n", ">>> $7 \\rightarrow 0111 \\rightarrow thhh \\rightarrow X = 3 \\rightarrow q p^3$ \n", ">>> $8 \\rightarrow 1000 \\rightarrow httt \\rightarrow X = 1 \\rightarrow q^3 p$ \n", ">>> $9 \\rightarrow 1001 \\rightarrow htth \\rightarrow X = 2 \\rightarrow q^2 p^2$ \n", ">>> $10 \\rightarrow 1010 \\rightarrow htht \\rightarrow X = 2 \\rightarrow q^2 p^2$ \n", ">>> $11 \\rightarrow 1011 \\rightarrow hthh \\rightarrow X = 3 \\rightarrow q p^3$ \n", ">>> $12 \\rightarrow 1100 \\rightarrow hhtt \\rightarrow X = 2 \\rightarrow q^2 p^2$ \n", ">>> $13 \\rightarrow 1101 \\rightarrow hhth \\rightarrow X = 3\\rightarrow q p^3 $ \n", ">>> $14 \\rightarrow 1110 \\rightarrow hhht \\rightarrow X = 3 \\rightarrow q p^3 $ \n", ">>> $15 \\rightarrow 1111 \\rightarrow hhhh \\rightarrow X = 4 \\rightarrow p^4 $ \n", ">>>\n", ">> Test: all probabilities must add to $1$: Is that true? Let us see:\n", ">>>\n", ">>> $S = q^4 + q^3 p + q^3 p + q^2 p^2 + q^3 p + q^2 p^2 + q^2 p^2 + q p^3 + q^3 p + q^2 p^2 + q^2 p^2 + q p^3 + q^2 p^2 + q p^3 + q p^3 + p^4 $\n", ">>\n", ">> Reordering to place together similar terms, we get\n", ">>>\n", ">>> $S = q^4 + ( q^3 p + q^3 p + q^3 p + q^3 p ) + (q^2 p^2 + q^2 p^2 + q^2 p^2 + q^2 p^2 + q^2 p^2 + q^2 p^2 ) + (q p^3 + q p^3 + q p^3 + qp^3 ) + p^4 $\n", ">>\n", ">> Associating:\n", ">>>\n", ">>> $ S = q^4 + 4q^3 p + 6q^2 p^2 + 4 q p^3 + p 4 = (q + p)^4 = (p+q)^4 = 1 ^4 = 1$\n", ">>> such as demanded.\n", "\n", ">> *Prediction 1.* The shape of the Binomial Distribution is defined by the random variable $X$ given by:\n", ">>\n", ">>> $0 \\rightarrow 1q^4 $ \n", ">>> $1 \\rightarrow 4q^3p$ \n", ">>> $2 \\rightarrow 6q^2 p^2 $ \n", ">>> $3 \\rightarrow 4qp^3 $ \n", ">>> $3 \\rightarrow 1p^4 $ \n", ">>\n", ">> Prediction 2: If $p = q = 1/2$, the corresponding bell is symmetric and the probabilities are\n", ">>>\n", ">>> $0 \\rightarrow 1 (1/2)^4 = 1/16$ \n", ">>> $1 \\rightarrow 4q^3p = 4 (1/2)^3 (1/2) = 4 (1/2)^4 = 4/16$ \n", ">>> $2 \\rightarrow 6q^2 p^2 = 6 (1/2)^2 (1/2)^2 = 6 (1/2)^4 = 6/16 $ \n", ">>> $3 \\rightarrow 4 q p^3 =4 (1/2) (1/2)^3 = 4 (1/2)^4 = 4/16$ \n", ">>> $4 \\rightarrow 1p^4 = (1/2)^4 = 1/16$ \n", ">>>\n", ">> So, the heights of the bars are in relation 1:4:6:4:1\n", "\n", "\n", "What does come further? We need to generalize.\n", "\n", "We need a general formula for the expansion of $(a+b)^n$ and a proof thereof.\n", "\n", "The general formula is\n", ">\n", "> $(q + p)^n=\\sum_{k=0}^{n}{n!\\over{(n-k)!k!}}q^kp^{n-k}$\n", "\n", "where the factorial $n!$ is defined recursively:\n", ">\n", "> $0! = 1 $ and $n! = (n-1)! n$,\n", "\n", "> **Example.** $5! = 4! \\times 5 = 3! \\times 4 \\times 5 = 2! \\times 3 \\times 4 \\times 5 = 1! \\times 2 \\times 3 \\times 4 \\times 5 = 0! \\times 1 \\times 2 \\times 3 \\times 4 \\times 5 = 1 \\times 1 \\times 2 \\times 3 \\times 4 \\times 5 = 120$.\n", "\n", "> **Exercise.** Calculate $7!$.\n", "\n", "> **Answer.** $7! = 6! \\times 7 = 5! \\times 6 \\times 7 = 120 \\times 6 \\times 7 \\sim 5000 $\n", "\n", "Let us show that the Binomial Formula reproduces what we already know.\n", "\n", "> **Example.** Let us prove that the formula is correct for $n= 1$. We must show that\n", ">\n", ">> $(q + p)^n=\\sum_{k=0}^{n}{n!\\over{(n-k)!k!}}q^kp^{n-k}$\n", " reduces for $n = 1$ to\n", ">>\n", ">> $ (p+q)^1 = p+q$.\n", ">\n", "> This is achieved by replacing the factorials:\n", ">>\n", ">> $(q + p)^1=\\sum_{k=0}^{1}{1!\\over{(1-k)!k!}}q^kp^{1-k}$\n", ">\n", "> This involves a sum with two terms, the first for $k= 0$ and the second for $k = 1$:\n", ">>\n", ">> $(q + p)^1 = {1!\\over{(1-0)!0!}}q^0p^{1-0} + {1!\\over{(1-1)!1!}}q^1p^{1-1} = $ \n", ">> $ = {1\\over{1}} p^{1} + {1\\over{1}}q^1p^{0} $ \n", ">> $ = p + q $ \n", "\n", "\n", "Thus, we have proved that the formula is correct for $n = 1$.\n", "\n", "\n", ">> **Example.**. Let us verify that the Binomial Expansion is correct for $n = 2$.\n", "\n", "> The formula for $n = 2 $ reads\n", ">>\n", ">> $(q + p)^2=\\sum_{k=0}^{2}{2!\\over{(2-k)!k!}}q^kp^{2-k}$\n", ">>\n", "> and has three terms, for $k = 0, 1, 2$:\n", ">>\n", ">> $(q + p)^2= {2!\\over{(2-0)!0!}}q^0p^{2-0}$ \n", ">> $ + {2!\\over{(2-1)!1!}}q^1p^{2-1}$ \n", ">> $+ {2!\\over{(2-2)!2!}}q^2p^{2-2}$ \n", ">> $= {2\\over{2}}p^2$ \n", ">> $ + {2\\over{1}}q p$ \n", ">> $+ {2\\over{2}}q^2$ \n", ">> $= p^2 + 2pq + q^2 $, such as it was found before.\n", "\n", "\n", "> **Exercise.** Verify that the Binomial Expansion is correct for $n = 3$.\n", ">> **Answer.** We have 4 terms for $k = 0,1,2,3$ \n", ">> $(q + p)^3=\\sum_{k=0}^{3}{3!\\over{(3-k)!k!}}q^kp^{3-k} $ \n", ">> $ = {3!\\over{(3-0)!0!}}q^0p^{3-0} + {3!\\over{(3-1)!1!}}q^1p^{3-1} + {3!\\over{(3-2)!2!}}q^2p^{3-2} + {3!\\over{(3-3)!3!}}q^3p^{3-3} $ \n", ">> $ = {3!\\over{(3)!0!}}q^0p^{3} + {3!\\over{(2)! }}q p^{2} + {3!\\over{1!2!}}q^2p + {3!\\over{0!3!}}q^3p^0 $ \n", ">> $ = {6\\over{6 }} p^{3} + {6\\over{2 }}q p^{2} +{6\\over{2}}q^2p +\n", "{6\\over{6}}q^3 $ \n", ">> $ = p^{3} + 3q p^{2} +3q^2p + q^3 $ \n", "\n", "\n", "The general proof of the Binomial Expansion or Newton Binomial asserts that the formula is valid for each $n$ from 1 to $\\infty$. This must be done by applying the **Mathematical Principle of Induction** that states: to rise up an infinite ladder, it is enough to fulfil two conditions:\n", "- You can reach the first step,\n", "- You can always rise to step $n+1$ whenever you are in step $n$. \n", " \n", "So, to prove the Binomial Theorem, we need to prove it for the first term, something that we already redundantly did for $n= 1, 2, 3 $. In second place, we must show that we can pass from any step to the next one.\n", "\n", "\n", "> **Research.** Assume that we are at step $n$, i.e., that the expansion formula is valid for the power $n$. From that, you must prove that you can rise to the following step, i.e., that the binomial expansion is valid for $n+1$. In other words, you must assume that\n", ">\n", ">> $(q + p)^n=\\sum_{k=0}^{n}{n!\\over{(n-k)!k!}}q^kp^{n-k}$\n", ">\n", "> Then, you must multiply both sides of this equation by $(q+p)$\n", ">>\n", ">> $(q+p) (q + p)^n=(q+p) \\sum_{k=0}^{n}{n!\\over{(n-k)!k!}}q^kp^{n-k}$\n", ">\n", "> and must show that you find:\n", ">>\n", ">> $(q + p)^{n+1}=\\sum_{k=0}^{n+1}{(n+1)!\\over{(n+1-k)!k!}}q^kp^{n+1-k}$\n", ">\n", "> This can be found in Wikipedia:\n", ">\n", "> https://en.wikipedia.org/wiki/Binomial_theorem\n", ">\n", "> together with a prerequisite, the Pascal Triangle,\n", ">\n", "> https://en.wikipedia.org/wiki/Pascal%27s_rule\n", ">\n", "> We also need the general relation of all this algebraic procedures with probability, more exactly with combinatorics:\n", ">\n", "> https://sboyles.github.io/teaching/ce311s/4-combinatorics.pdf\n", "\n", "\n", "\n", "With all this, we can state now what we believe in:\n", "\n", "Let it be an experiment with $n$ coins, whose probability of heads is $p$ and that of tails is $q$, in which we count the number of heads. We repeat the experiment over and over to get an estimation of the Binomial Distribution, which predicts the relative frequency of $k$ heads. Then, the relative frequency of $k$ heads observed in the experiment tends to the predicted probability by the Binomial Expansion: ${n!\\over{(n-k)!k!}}p^kq^{n-k}$.\n", "\n", "If we have a belief that can be tested, which is the next step in Science? To test it.\n", "\n", "# TESTING WHAT WE BELIEVE IN\n", "\n", "We must make our simulations, choose $n$, and for each $k \\le n$ compute the relative frequency of $k$ heads, and compare all those frequencies with the predicted values dictated by our belief.\n", "\n", "Our first step is to see what Python has concerning the Binomial Formula.\n", "\n", "\n", "\n", "> **Example.** We prompted Colab with\n", ">>\n", ">> Code for the binomial expansion (p+q)^n\n", ">\n", "> It answered us with the following code:\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "id": "kd_teeq5OehW" }, "outputs": [], "source": [ "# program binomialExpansionByGemini.py\n", "\n", "import math\n", "\n", "def binomial_expansion(p, q, n):\n", " \"\"\"\n", " Calculates the terms of the binomial expansion (p+q)^n.\n", "\n", " Args:\n", " p: The probability of the first outcome.\n", " q: The probability of the second outcome (q = 1 - p).\n", " n: The number of trials.\n", "\n", " Returns:\n", " A dictionary where keys are the number of successes (k)\n", " and values are the corresponding probabilities.\n", " \"\"\"\n", " if not math.isclose(p + q, 1.0):\n", " print(\"Warning: p + q should be equal to 1.\")\n", "\n", " expansion_terms = {}\n", " for k in range(n + 1):\n", " # Calculate the binomial coefficient (n choose k)\n", " binomial_coeff = math.comb(n, k)\n", "\n", " # Calculate the probability of k successes\n", " probability_k = binomial_coeff * (p ** k) * (q ** (n - k))\n", "\n", " expansion_terms[k] = probability_k\n", " return expansion_terms\n", "\n", "# Example usage for a fair coin (p=0.5, q=0.5) with 10 flips (n=10)\n", "p = 0.5\n", "q = 0.5\n", "n = 10\n", "binomial_probs = binomial_expansion(p, q, n)\n", "print(f\"Binomial expansion for (0.5 + 0.5)^{n}:\")\n", "for k, prob in binomial_probs.items():\n", " print(f\" {k} heads: {prob:.4f}\")\n", "\n", "# Example usage for a biased coin (p=2/3, q=1/3) with 10 flips (n=10)\n", "p = 2/3\n", "q = 1/3\n", "n = 10\n", "binomial_probs_biased = binomial_expansion(p, q, n)\n", "print(f\"\\nBinomial expansion for (2/3 + 1/3)^{n}:\")\n", "for k, prob in binomial_probs_biased.items():\n", " print(f\" {k} heads: {prob:.4f}\")" ] }, { "cell_type": "markdown", "metadata": { "id": "dmFJKMnYOc3S" }, "source": [ "> **Exercise.** Run the code and evaluate it.\n", ">>\n", ">> **Answer.** The program shows two cases. Grosso modo, they seem to be OK. But, to better assess it, we used the code to calculate (0.5 + 0.5)². The answer was: 0.25, 0.5, 0.25, which was OK. So, we provided feedback to Colab with a good point.\n", "\n", "> **Exercise.** Append to the program `tendencyToBinomialGraphGeneral.py`, developed by Gemini, the code for the Binomial Expansion to contrast the last simulated distribution with the mathematically expected limit. \n", ">>\n", ">> **Answer.** The following code shows the tendency of observed distributions of the number of heads for n coins as the number of repetitions increases and contrasts its limit with the mathematically expected distribution. Expectations are fulfilled with flying colors.\n", "\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "id": "CaH7GfFFUrqq" }, "outputs": [], "source": [ "# Program tendencyToBinomialGraphGeneral2.py\n", "# This program shows the tendency of\n", "# the sampled distribution of the number of heads\n", "# and compares its limit with\n", "# the mathematically expected distribution\n", "\n", "# Libraries and local rename\n", "import math\n", "import matplotlib.pyplot as plt\n", "import numpy as np # Support for scientific computing\n", "import random\n", "from collections import Counter\n", "from mpl_toolkits.mplot3d import Axes3D\n", "\n", "# Parameters for the simulation\n", "flipping_numbers_to_be_sampled = [10, 100, 1000, 10000]\n", "num_coins_per_round = 20\n", "\n", "# Simulate coin flips\n", "# Returns an array with the number of heads per round\n", "def simulate_coin_flips(num_coins_per_round, num_repetitions):\n", " results = []\n", " for _ in range(num_repetitions):\n", " heads_in_round = 0\n", " for _ in range(num_coins_per_round):\n", " # 0 for tails, 1 for heads\n", " if random.choice([0, 1]) == 1:\n", " heads_in_round += 1\n", " results.append(heads_in_round)\n", " return results\n", "\n", "\n", "# Calculate statistics\n", "# Returns the sorted distribution of heads\n", "# as an array with num_heads and another with the corresponding relative frequency\n", "def calculate_statistics(results, num_coins_per_round, index_num_repetitions, num_repetitions, z_heights):\n", " print(\"\\n calculate_statistics\")\n", " total_heads = sum(results)\n", " proportion_heads = total_heads / (num_coins_per_round * num_repetitions)\n", " heads_distribution = Counter(results)\n", "\n", " # Print the results\n", " print(f\"Simulation of {num_repetitions} rounds \")\n", " print(f\"with {num_coins_per_round} coins per round:\")\n", " print(f\"\\nTotal number of heads across all repetitions: {total_heads}\")\n", " print(f\"Proportion of heads across all repetitions: {proportion_heads:.4f}\")\n", " print(\"\\nDistribution of heads:\")\n", " # Sort the distribution for better readability\n", " sorted_distribution = dict(sorted(heads_distribution.items()))\n", " for num_heads, count in sorted_distribution.items():\n", " relative_frequency = count/num_repetitions\n", " print(f\"{num_heads} heads: {count} times, {relative_frequency} relative frequency\")\n", " if num_heads <= num_coins_per_round and index_num_repetitions < len(flipping_numbers_to_be_sampled):\n", " z_heights[index_num_repetitions, num_heads] = relative_frequency\n", " print(f\"z_heights[{index_num_repetitions},{num_heads}] = {z_heights[index_num_repetitions, num_heads]}\")\n", "\n", " print(f\"z_heights {z_heights}\")\n", " return z_heights\n", "\n", "\n", "# Find and plot the observed distribution of the number of heads\n", "# for various flip values for a fix number of coins.\n", "def plot_hyper_distribution(num_coins_per_round):\n", " print(\"plot_hyper_distribution(num_coins_per_round):\")\n", " print(f\"num_coins_per_round = {num_coins_per_round}\")\n", " print(f\"flipping_numbers_to_be_sampled = {flipping_numbers_to_be_sampled}\")\n", " numb_flipping_values = len(flipping_numbers_to_be_sampled)\n", " print(f\"len(flipping_numbers_to_be_sampled) = {numb_flipping_values}\")\n", "\n", " # Data for the graph\n", " x_data = np.arange(num_coins_per_round + 1)\n", " print(f\"x_data = {x_data} \\n\")\n", " y_data = np.arange(numb_flipping_values)\n", " print(f\"y_data = {y_data} \\n\")\n", "\n", " # Create a meshgrid for the x and y coordinates of the bars\n", " _xx, _yy = np.meshgrid(x_data, y_data)\n", " print(f\"_xx = {_xx} \\n\")\n", " print(f\"_yy = {_yy} \\n\")\n", " x_pos, y_pos = _xx.ravel(), _yy.ravel()\n", " print(f\"x_pos = {x_pos} \\n\")\n", " print(f\"y_pos = {y_pos} \\n\")\n", " z_pos = np.zeros_like(x_pos) # All bars start at z=0\n", "\n", " z_heights = 1.0 * np.zeros((numb_flipping_values, num_coins_per_round + 1)) # z_heights initialization\n", " print(f\"z_heights initialized = {z_heights} \\n\")\n", "\n", "\n", " for index_num_repetitions in range(numb_flipping_values):\n", " print(\"\\n=============================\\n\")\n", " print(f\"index_num_repetitions = {index_num_repetitions}\")\n", " num_repetitions = flipping_numbers_to_be_sampled[index_num_repetitions]\n", " print(f\"Flipping number = {num_repetitions}\")\n", " results = simulate_coin_flips(num_coins_per_round, num_repetitions)\n", " z_heights = calculate_statistics(results, num_coins_per_round, index_num_repetitions, num_repetitions, z_heights)\n", "\n", "\n", " dz = z_heights.ravel() # Heights of the bars\n", " print(f\"dz = {dz} \\n\\n\")\n", "\n", " print(\"Is there a tendency of the sampled distributions\")\n", " print(\"when the number of repetitions increases?\")\n", " for i in range(numb_flipping_values):\n", " print(f\"Number of repetitions {flipping_numbers_to_be_sampled[i]}: sampled distribution{z_heights[i]}\")\n", "\n", " # Create the 3D plot\n", " fig = plt.figure(figsize=(10, 8))\n", " ax = fig.add_subplot(111, projection='3d')\n", "\n", " # Define the dimensions of the bars\n", " dx = dy = 0.6 # Width and depth of the bars\n", "\n", " # Plot the 3D bars\n", " ax.bar3d( y_pos, x_pos, z_pos, dx, dy, dz, shade=True)\n", "\n", " # Set labels and title\n", " ax.set_xlabel('Number of Repetitions (10^x)')\n", " ax.set_ylabel('Number of heads')\n", " ax.set_zlabel('Relative Frequency')\n", " ax.set_title('Tendency towards the Binomial Distribution')\n", "\n", " # Set the tick labels for the y-axis (number of repetitions)\n", " ax.set_yticklabels([f'10^{i+1}' for i in range(numb_flipping_values)])\n", "\n", " # Display the plot\n", " plt.show()\n", " return z_heights #To compare observed with expected\n", "\n", "\n", "\n", "\n", " # The binomial expansion\n", " def binomial_expansion(p, q, n):\n", " \"\"\"\n", " Calculates the terms of the binomial expansion (p+q)^n.\n", "\n", " Args:\n", " p: The probability of the first outcome.\n", " q: The probability of the second outcome (q = 1 - p).\n", " n: The number of trials.\n", "\n", " Returns:\n", " A dictionary where keys are the number of successes (k)\n", " and values are the corresponding probabilities.\n", " \"\"\"\n", " if not math.isclose(p + q, 1.0):\n", " print(\"Warning: p + q should be equal to 1.\")\n", "\n", " expansion_terms = {}\n", " for k in range(n + 1):\n", " # Calculate the binomial coefficient (n choose k)\n", " binomial_coeff = math.comb(n, k)\n", "\n", " # Calculate the probability of k successes\n", " probability_k = binomial_coeff * (p ** k) * (q ** (n - k))\n", "\n", " expansion_terms[k] = probability_k\n", " return expansion_terms\n", "\n", "# Binomial Expansion for a fair coin (p=0.5, q=0.5)\n", "# This gives the Mathematically Expected Binomial Distribution\n", "p = 0.5\n", "q = 0.5\n", "n = num_coins_per_round\n", "binomial_probs = binomial_expansion(p, q, n)\n", "print(f\"Binomial expansion for (0.5 + 0.5)^{n}:\")\n", "for k, prob in binomial_probs.items():\n", " print(f\" {k} heads: {prob:.4f}\")\n", "\n", "# Simulate the Binomial Distribution for various flipp numbers\n", "# to unveil the tendency towards a limit, if only it exists.\n", "# The matrix with all sampled distribution is returned.\n", "z_heights = plot_hyper_distribution(num_coins_per_round)\n", "\n", "# Contrast of expected with expected\n", "print(\"\\nComparsion between the simulated and expected Binomial Distributions\")\n", "# The last observed distribution is the simulated estimation of the limit\n", "lastRow = len(z_heights)-1\n", "for k, prob in binomial_probs.items():\n", " print(f\"{k} heads, observed: {prob:.4f} expected: {z_heights[lastRow, k]}\")\n", "\n" ] }, { "cell_type": "markdown", "metadata": { "id": "wCPYA893exEo" }, "source": [ "> **Exercise.** Make your best to have the previous program produce a graph with the correct scales for each axis. The number of repetitions must be `[10, 100, 1000, 10000]` that is given by the 1D-array `flipping_numbers_to_be_sampled`. The number of heads must run from `0` to `num_coins_per_round = 20`. The axis of relative frequency correctly shows values between 0 and 1. \n", ">>\n", ">> **Answer.** The documentation for 3D graphics seems to be incomplete. In the meantime, some minor adjustments are proposed in the following code:\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "id": "fx2fh9Tls43L" }, "outputs": [], "source": [ "# Program tendencyToBinomialGraphGeneral3.py\n", "# This program shows the tendency of\n", "# the sampled distribution of the number of heads\n", "# and compares its limit with\n", "# the mathematically expected distribution\n", "\n", "# Libraries and local rename\n", "import math\n", "import matplotlib.pyplot as plt\n", "import numpy as np # Support for scientific computing\n", "import random\n", "from collections import Counter\n", "\n", "\n", "# Parameters for the simulation\n", "flipping_numbers_to_be_sampled = [10, 100, 1000, 10000]\n", "num_coins_per_round = 20\n", "\n", "# Simulate coin flips\n", "# Returns an array with the number of heads per round\n", "def simulate_coin_flips(num_coins_per_round, num_repetitions):\n", " results = []\n", " for _ in range(num_repetitions):\n", " heads_in_round = 0\n", " for _ in range(num_coins_per_round):\n", " # 0 for tails, 1 for heads\n", " if random.choice([0, 1]) == 1:\n", " heads_in_round += 1\n", " results.append(heads_in_round)\n", " return results\n", "\n", "\n", "# Calculate statistics\n", "# Returns the sorted distribution of heads\n", "# as an array with num_heads and another with the corresponding relative frequency\n", "def calculate_statistics(results, num_coins_per_round, index_num_repetitions, num_repetitions, z_heights):\n", " print(\"\\n calculate_statistics\")\n", " total_heads = sum(results)\n", " proportion_heads = total_heads / (num_coins_per_round * num_repetitions)\n", " heads_distribution = Counter(results)\n", "\n", " # Print the results\n", " print(f\"Simulation of {num_repetitions} rounds \")\n", " print(f\"with {num_coins_per_round} coins per round:\")\n", " print(f\"\\nTotal number of heads across all repetitions: {total_heads}\")\n", " print(f\"Proportion of heads across all repetitions: {proportion_heads:.4f}\")\n", " print(\"\\nDistribution of heads:\")\n", " # Sort the distribution for better readability\n", " sorted_distribution = dict(sorted(heads_distribution.items()))\n", " for num_heads, count in sorted_distribution.items():\n", " relative_frequency = count/num_repetitions\n", " print(f\"{num_heads} heads: {count} times, {relative_frequency} relative frequency\")\n", " if num_heads <= num_coins_per_round and index_num_repetitions < len(flipping_numbers_to_be_sampled):\n", " z_heights[index_num_repetitions, num_heads] = relative_frequency\n", " print(f\"z_heights[{index_num_repetitions},{num_heads}] = {z_heights[index_num_repetitions, num_heads]}\")\n", "\n", " print(f\"z_heights {z_heights}\")\n", " return z_heights\n", "\n", "\n", "# Find and plot the observed distribution of the number of heads\n", "# for various flip values for a fix number of coins.\n", "def plot_hyper_distribution(num_coins_per_round):\n", " print(\"plot_hyper_distribution(num_coins_per_round):\")\n", " print(f\"num_coins_per_round = {num_coins_per_round}\")\n", " print(f\"flipping_numbers_to_be_sampled = {flipping_numbers_to_be_sampled}\")\n", " numb_flipping_values = len(flipping_numbers_to_be_sampled)\n", " print(f\"len(flipping_numbers_to_be_sampled) = {numb_flipping_values}\")\n", "\n", " # Data for the graph\n", " x_data = np.arange(num_coins_per_round + 1)\n", " print(f\"x_data = {x_data} \\n\")\n", " y_data = np.arange(numb_flipping_values)\n", " print(f\"y_data = {y_data} \\n\")\n", "\n", " # Create a meshgrid for the x and y coordinates of the bars\n", " _xx, _yy = np.meshgrid(x_data, y_data)\n", " print(f\"_xx = {_xx} \\n\")\n", " print(f\"_yy = {_yy} \\n\")\n", " # Interchange for the graph\n", " x_pos, y_pos = _yy.ravel(),_xx.ravel()\n", " print(f\"x_pos = {x_pos} \\n\")\n", " print(f\"y_pos = {y_pos} \\n\")\n", " z_pos = np.zeros_like(x_pos) # All bars start at z=0\n", "\n", " z_heights = 1.0 * np.zeros((numb_flipping_values, num_coins_per_round + 1)) # z_heights initialization\n", " print(f\"z_heights initialized = {z_heights} \\n\")\n", "\n", "\n", " for index_num_repetitions in range(numb_flipping_values):\n", " print(\"\\n=============================\\n\")\n", " print(f\"index_num_repetitions = {index_num_repetitions}\")\n", " num_repetitions = flipping_numbers_to_be_sampled[index_num_repetitions]\n", " print(f\"Flipping number = {num_repetitions}\")\n", " results = simulate_coin_flips(num_coins_per_round, num_repetitions)\n", " z_heights = calculate_statistics(results, num_coins_per_round, index_num_repetitions, num_repetitions, z_heights)\n", "\n", "\n", " dz = z_heights.ravel() # Heights of the bars\n", " print(f\"dz = {dz} \\n\\n\")\n", "\n", " print(\"Is there a tendency of the sampled distributions\")\n", " print(\"when the number of repetitions increases?\")\n", " for i in range(numb_flipping_values):\n", " print(f\"Number of repetitions {flipping_numbers_to_be_sampled[i]}: sampled distribution{z_heights[i]}\")\n", "\n", " # Create the 3D plot\n", " fig = plt.figure(figsize=(10, 8))\n", " ax = fig.add_subplot(111, projection='3d')\n", "\n", " # Define the dimensions of the bars\n", " dx = dy = 0.6 # Width and depth of the bars\n", "\n", " # Plot the 3D bars\n", " ax.bar3d( x_pos, y_pos, z_pos, dx, dy, dz, shade=True)\n", "\n", " # Set labels and title\n", " ax.set_xlabel('Number of Repetitions = (10^(x+1))')\n", "\n", " ax.set_ylabel('Number of heads - 1')\n", " ax.set_zlabel('Relative Frequency')\n", " ax.set_title('Tendency towards the Binomial Distribution')\n", "\n", " # Set the tick labels for the y-axis (number of repetitions)\n", " # ax.set_xticklabels([f'10^{i+1}' for i in range(numb_flipping_values)])\n", "\n", " # Display the plot\n", " plt.show()\n", " return z_heights\n", "\n", "\n", "\n", "\n", "# The binomial expansion\n", "def binomial_expansion(p, q, n):\n", " \"\"\"\n", " Calculates the terms of the binomial expansion (p+q)^n.\n", "\n", " Args:\n", " p: The probability of the first outcome.\n", " q: The probability of the second outcome (q = 1 - p).\n", " n: The number of trials.\n", "\n", " Returns:\n", " A dictionary where keys are the number of successes (k)\n", " and values are the corresponding probabilities.\n", " \"\"\"\n", " if not math.isclose(p + q, 1.0):\n", " print(\"Warning: p + q should be equal to 1.\")\n", "\n", " expansion_terms = {}\n", " for k in range(n + 1):\n", " # Calculate the binomial coefficient (n choose k)\n", " binomial_coeff = math.comb(n, k)\n", "\n", " # Calculate the probability of k successes\n", " probability_k = binomial_coeff * (p ** k) * (q ** (n - k))\n", "\n", " expansion_terms[k] = probability_k\n", " return expansion_terms\n", "\n", "# Binomial Expansion for a fair coin (p=0.5, q=0.5)\n", "# This gives the Mathematically Expected Binomial Distribution\n", "p = 0.5\n", "q = 0.5\n", "n = num_coins_per_round\n", "binomial_probs = binomial_expansion(p, q, n)\n", "print(f\"Binomial expansion for (0.5 + 0.5)^{n}:\")\n", "for k, prob in binomial_probs.items():\n", " print(f\" {k} heads: {prob:.4f}\")\n", "\n", "# Simulate the Binomial Distribution for various flipp numbers\n", "# to unveil the tendency towards a limit, if only it exists.\n", "# The matrix with all sampled distribution is returned.\n", "z_heights = plot_hyper_distribution(num_coins_per_round)\n", "\n", "# Contrast of expected with expected\n", "print(\"\\nComparsion between the simulated and expected Binomial Distributions\")\n", "# The last observed distribution is the simulated estimation of the limit\n", "lastRow = len(z_heights)-1\n", "for k, prob in binomial_probs.items():\n", " print(f\"{k} heads, Observed: {prob:.4f} Expected: {z_heights[lastRow, k]}\")\n", "\n", "\n", "\n" ] }, { "cell_type": "markdown", "metadata": { "id": "8VD3Ji5D9zW-" }, "source": [ "> **Challenge.** The code of the previous program `tendencyToBinomialGraphGeneral3.py` runs almost perfectly, but looks very unstructured. It causes fatigue. Structure it wisely so that making a change becomes naturally easy. In the jargon, we say: develop an evolvable software.\n", "\n", "\n", "> **PARTIAL CONCLUSION.** We have found that the relative frequency of heads and the distribution of heads tend to a limit that is naturally characterized by a mathematical model: the Binomial Distribution. In that way, we have related theory with observation and have found that their correspondence is quite good.\n", "\n", "\n", "What is the use in Science of the finding of the Binomial Distribution as a limit of observed binomial distributions?\n", "\n", "# THE QUEST FOR THE TRUTH\n", "\n", "One grows with the illusion that Science will show us the truth and that we will always know the true answer to the question: Is this theory correct? In that way, we picture the truth in one of two values: Yes else No. Let us see that when randomness is an important factor, as in every applied discipline, we must change our focus in a way that directly arises from our procedures.\n", "\n", "Let's suppose that we must decide whether the theory $p = 0.5$ is true given that we have tossed a coin 20 times, and we found 14 heads and 6 tails. In general, all scientific studies can be reformulated as a comparison between what one thinks and what one sees, between expected according to a given belief and observed according to experiment, when all important factors are controlled. Here we have:\n", "\n", "> Belief: $H_o$ that says $p= 0.5$ \n", "> Observed: 14 heads out of 20.\n", "\n", "To make a decision, the first task is to calculate the $H_0$. This means gathering too much data to build an estimation of the distribution of heads when we play with 20 coins. Let us recall what we have learned thus far:\n", "\n", "***It is true that randomness interdicts our capability to predict what will be next. But there are tendencies to a limit in ensemble observables such as the proportion of heads, the mode and mean of the numbers of heads, and the silhouette of the distribution of heads. To reaffirm the tendencies, all we must do is to augment the number of experiments. Tendencies to a limit provide us with a kind of stability and robustness, and that is why we say that 10000 means infinite.***\n", "\n", "We have already calculated the $H_o$ that says $p= 0.5$. In fact, we have found that the observed distribution of the number of heads tends to a limit, which, for 10000 repeats, insignificantly differs from the Binomial Distribution.\n", "\n", "\n", "Let us reproduce here the Observed Distribution that was found previously:\n", "\n", "Binomial expansion for (0.5 + 0.5)^20:\n", "\n", " 0 heads: 0.0000\n", "\n", " 1 heads: 0.0000\n", "\n", " 2 heads: 0.0002\n", "\n", " 3 heads: 0.0011\n", "\n", " 4 heads: 0.0046\n", "\n", " 5 heads: 0.0148\n", "\n", " 6 heads: 0.0370\n", "\n", " 7 heads: 0.0739\n", "\n", " 8 heads: 0.1201\n", "\n", " 9 heads: 0.1602\n", "\n", " 10 heads: 0.1762\n", "\n", " 11 heads: 0.1602\n", "\n", " 12 heads: 0.1201\n", "\n", " 13 heads: 0.0739\n", "\n", " 14 heads: 0.0370\n", "\n", " 15 heads: 0.0148\n", "\n", " 16 heads: 0.0046\n", "\n", " 17 heads: 0.0011\n", "\n", " 18 heads: 0.0002\n", "\n", " 19 heads: 0.0000\n", "\n", " 20 heads: 0.0000\n", "\n", "Let us notice that every single event has a low relative frequency that, anyway, fades off towards the extremes 0 and 20. More to the point, our simulation with 10000 repetitions does not produce events with a number of heads in the set 0, 1, 19, 20. This is no miracle, as that is observed in real life:\n", "\n", "> **Example.** Suppose that we make a random sample of 20 people and register their gender, male or female. We do not expect to find all 20 people to be men or to be all women, or to get 19 men and 1 woman, or 19 women and 1 man. Rather, we will find values around 10 men and 10 women. Comment: At birth, there are a bit more males than females, but for adult people, one might find a slight excess of women because the mortality of men is a bit higher than that of females. What do we do when we find 18 women and 2 men out of 20 people? We will always make a new sample, a new observation. If the situation persists, we will ask: Where are the men? To this, they will answer that they are in the army or in the land. \n", "\n", "What will we do once we have calculated the $H_o$ and have a distribution of heads?\n", "\n", "We divide the world into two regions: the first must support the $H_o$, inviting to accept it, while the second must lift suspicion, inviting for a rejection. Is there a rule to divide the world between rejecting and accepting regions? No. There is not. This means that human beings divide the world according to their objectives, or free will. \n", "\n", "To classify an event as suspicious is not an objective decision. It is subjective. We have said that values 0,1, 19, and 20 are suspicious because the relative frequency of these events is null. But this is only one option. People may declare in real life very different proportions of suspicious events.\n", "\n", "> **Example.** Let us think of two candidates, A and B. They appeared on a TV show. They are totally unknown, so neither of them is a favorite. The probability of voting for one of them at that moment is equal to 0.5. But if, after the show, many people are speaking with emotion about candidate B and no one remembers candidate A, we know that the people have a candidate. And everyone would want to know what he says, and after listening to his proposal, they would exclaim: This candidate knows how to tie the ends! The corresponding statistical description would assign 0.95 to the probability of voting for candidate B and 0.05 for candidate A. Thus, the initial hypothesis of indifference is rejected, and a new one is adopted. This implies that candidate A will study the discourse of B, and implement a winner politics. Indeed, candidate A will fire the alarms when polls indicate any advantage for B, so politics appears to work with a proportion of rejecting events equal to 0.99 for the difference in prospects for votes.\n", "\n", "The general behavior is abstracted as expected:\n", "\n", "If we have a supposition that we call the null hypothesis, $H_o$, we use our free will to classify events as normal or accepting, else as abnormal or rejecting. A normal event invites accepting the null hypothesis, while a strange event invites rejecting the null hypothesis and proposing another hypothesis. The exact division is arbitrary, subjective. Every classification of events as normal or else strange is a human decision, a human responsibility. Let us enphasize it:\n", "\n", ">> Science is a human creation, a human responsibility. \n", ">> Science is composed of artworks but not of dogmas.\n", "\n", "**Science has opted for the following default rule: The proportion of ordinary events is defined to be 0.95, and that of strange events is 0.05.** For a central distribution, bell-like, one whose silhouette looks like a one picked mountain that falls away in the distance, strange events are at the far reaches of the universe, and they must add together to 5% or 0.05. If we have many data and a great discovery, this barrier could be lowered to 0.01, or to 0.001, or even less. It is usual to fix the barrier at 0.05 or 5% because it gives a good proportion of events that invite to revise beliefs and that can be found with samples of size 20 or 30. To be more rigorous, the barrier must be lowered to 0.01, which requires samples of size 100, 200. And so on.\n", "\n", "The notion of strange cannot always be applied.\n", "\n", "> **Example.** A roulette is divided into 360 numbered circular sectors of equal probability, 1/360. Is there a sector that deserves to be called strange? No, there is not. Does this mean that we cannot do Science with a roulette? No. Instead, we must create an environment in which it makes sense to divide the world into ordinary and strange. The following does it: The game consists of spinning the roulette twice and adding the two numbers marked by the arrow. The sum will have values that are very rare and others that are very common. For instance, the value 2 has only one possibility, 1 + 1, but the value 360 has too many possibilities, 359 + 1, 358 + 2, ..., 180 + 180, 179 + 181, 178+ 182,... 1 + 359. The value 360 is 359 times more probable than 2. Thus, it makes sense to say that 2 is in the rejection region and the rest are in the acceptance region. But of what $H_o$? Of the $H_o$ that says that all events in the roulette are equiprobable.\n", "\n", "> **Example.** Let us apply the default rule of Science to our case-study of tossing 20 coins with the $H_o$ that says that the probability of heads is 0.5. We define that the proportion of rare events must be 0.05. Because the bell is symmetric, and we do not pay attention to a specific side, we divide 0.05 in both tails, and assign 0.025 to each one. This produces the following table, in which we begin to cumulate from the extremes:\n", "\n", "Binomial expansion for (0.5 + 0.5)^20: \n", " 0 heads: 0.0000 cumulated 0.0 \n", " 1 heads: 0.0000 cumulated 0.0 \n", " 2 heads: 0.0002 cumulated 0.0002 \n", " 3 heads: 0.0011 cumulated 0.0013 \n", " 4 heads: 0.0046 cumulated 0.0059 \n", " 5 heads: 0.0148 cumulated 0.0207 \n", " 6 heads: 0.0370 cumulated 0.0577 \n", " 7 heads: 0.0739 \n", " 8 heads: 0.1201 \n", " 9 heads: 0.1602 \n", " 10 heads: 0.1762 \n", " 11 heads: 0.1602 \n", " 12 heads: 0.1201 \n", " 13 heads: 0.0739 \n", " 14 heads: 0.0370 cumulated 0.0577 \n", " 15 heads: 0.0148 cumulated 0.0207 \n", " 16 heads: 0.0046 cumulated 0.0059 \n", " 17 heads: 0.0011 cumulated 0.0013 \n", " 18 heads: 0.0002 cumulated 0.0002 \n", " 19 heads: 0.0000 cumulated 0.0 \n", " 20 heads: 0.0000 cumulated 0.0 \n", "\n", "The criterion of the measure of the rare events equal to 0.05 renders as strange, extreme, or rare, the events 0,1,2,3,4,5, and 15, 16, 17, 18, 19,20 whose relative frequencies add up together to 0.0414, almost 0.05. It seems like a long list. If this corresponded with the heights of persons, this depicts the behavior of women: they reject small men, at the inferior tail, and prefer tall ones, at the upper one. \n", "\n", "> **Example.** Let us develop a Python program that presents graphically the acceptance and rejection regions when the rejection region has a probability of 0.05.\n", "\n", "Our idea is to use a bar plot, ordinary events will be colored in blue, and strange ones in red. So, we looked at the documentation of MatplotLib for bar plots, and found the following example:\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "id": "nOr4n3QfB3vh" }, "outputs": [], "source": [ "import matplotlib.pyplot as plt\n", "\n", "fig, ax = plt.subplots()\n", "\n", "fruits = ['apple', 'blueberry', 'cherry', 'orange']\n", "counts = [40, 100, 30, 55]\n", "bar_labels = ['red', 'blue', '_red', 'orange']\n", "bar_colors = ['tab:red', 'tab:blue', 'tab:red', 'tab:orange']\n", "\n", "ax.bar(fruits, counts, label=bar_labels, color=bar_colors)\n", "\n", "ax.set_ylabel('fruit supply')\n", "ax.set_title('Fruit supply by kind and color')\n", "ax.legend(title='Fruit color')\n", "\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": { "id": "9eg_BRmCCW4h" }, "source": [ "> **Exercise** Run the code that draws a bar diagram and pay attention to the form as data are given.\n", "\n", ">> **Answer.** Data are of diverse type, and they must be given in separate vectors making sure that they follow the same order.\n", "\n", "\n", "From this example, we propose the following code that graphically represents our expected distribution in which rejection or rare events are distinctly colored form accepting or normal ones:" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "id": "Ch4QO48eCYEt" }, "outputs": [], "source": [ "# Program binomialDistriColoredTails.py\n", "# Binomial Expansion for a fair coin (p=0.5, q=0.5)\n", "# This gives the Mathematically Expected Binomial Distribution\n", "import math\n", "import matplotlib.pyplot as plt\n", "import numpy as np\n", "\n", "# The binomial expansion\n", "def binomial_expansion(p, q, n):\n", " \"\"\"\n", " Calculates the terms of the binomial expansion (p+q)^n.\n", "\n", " Args:\n", " p: The probability of the first outcome.\n", " q: The probability of the second outcome (q = 1 - p).\n", " n: The number of trials.\n", "\n", " Returns:\n", " A dictionary where keys are the number of successes (k)\n", " and values are the corresponding probabilities.\n", " \"\"\"\n", " if not math.isclose(p + q, 1.0):\n", " print(\"Warning: p + q should be equal to 1.\")\n", "\n", " expansion_terms = {}\n", " for k in range(n + 1):\n", " # Calculate the binomial coefficient (n choose k)\n", " binomial_coeff = math.comb(n, k)\n", "\n", " # Calculate the probability of k successes\n", " probability_k = binomial_coeff * (p ** k) * (q ** (n - k))\n", "\n", " expansion_terms[k] = probability_k\n", " return expansion_terms\n", "\n", "# Binomial Expansion for a fair coin (p=0.5, q=0.5)\n", "# This gives the Mathematically Expected Binomial Distribution\n", "p = 0.5\n", "q = 0.5\n", "n = num_coins_per_round = 20\n", "binomial_probs = binomial_expansion(p, q, n)\n", "cumulate_prob = 0.0\n", "print(f\"Binomial expansion for (0.5 + 0.5)^{n}:\")\n", "for k, prob in binomial_probs.items():\n", " cumulate_prob += prob\n", " print(f\" {k} heads: {prob:.4f} cumulate_prob: {cumulate_prob} \")\n", "\n", "####### The graphic with strange and normal events distinctly colored.\n", "\n", "fig, ax = plt.subplots()\n", "\n", "k_class = [] #Horizontal label\n", "counts = [] #vertical values\n", "bar_labels = [] #colors of bars\n", "bar_colors = []\n", "\n", "left_tail = 0.0\n", "right_tail = 0.0\n", "\n", "# Coloring\n", "for k, prob in binomial_probs.items():\n", " k_class.append(k)\n", " counts.append(prob)\n", " left_tail += prob\n", " if left_tail < 0.25:\n", " bar_labels.append('red')\n", " bar_colors.append('tab:red')\n", " if ( (left_tail > 0.25) & (left_tail < 0.75)):\n", " bar_labels.append('blue')\n", " bar_colors.append('tab:blue')\n", " if left_tail > 0.75:\n", " bar_labels.append('red')\n", " bar_colors.append('tab:red')\n", "\n", "\n", "ax.bar(k_class, counts, label=bar_labels, color=bar_colors)\n", "\n", "ax.set_ylabel('probability')\n", "ax.set_title('Binomial Distribution with rejection regions in red')\n", "# ax.legend(title='Fruit color')\n", "\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": { "id": "GJFtuKJzizJv" }, "source": [ "> **Exercise.** Run the previous code and assess it.\n", "\n", ">> **Answer.** The code colors all bars correctly except for one.\n", "\n", "\n", "> **Exercise.** Our previous code produces a graphic with asymmetric red regions. Correct this.\n", "\n", "\n", ">> **Answer.** The following code produces a symmetric graph when the probability of rare events is 0.05. In reality, because of the discreteness of data, the probability of rare events is the maximal possible value less than 0.05. Its value is appended to the title of the graphic.\n", "\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "id": "GN2-RuLCRETa" }, "outputs": [], "source": [ "# Program binomialDistriColoredTails2.py\n", "# Binomial Expansion for a fair coin (p=0.5, q=0.5)\n", "# This gives the Mathematically Expected Binomial Distribution\n", "import math\n", "import matplotlib.pyplot as plt\n", "import numpy as np\n", "\n", "# The binomial expansion\n", "def binomial_expansion(p, q, n):\n", " \"\"\"\n", " Calculates the terms of the binomial expansion (p+q)^n.\n", "\n", " Args:\n", " p: The probability of the first outcome.\n", " q: The probability of the second outcome (q = 1 - p).\n", " n: The number of trials.\n", "\n", " Returns:\n", " A dictionary where keys are the number of successes (k)\n", " and values are the corresponding probabilities.\n", " \"\"\"\n", " if not math.isclose(p + q, 1.0):\n", " print(\"Warning: p + q should be equal to 1.\")\n", "\n", " expansion_terms = {}\n", " for k in range(n + 1):\n", " # Calculate the binomial coefficient (n choose k)\n", " binomial_coeff = math.comb(n, k)\n", "\n", " # Calculate the probability of k successes\n", " probability_k = binomial_coeff * (p ** k) * (q ** (n - k))\n", "\n", " expansion_terms[k] = probability_k\n", " return expansion_terms\n", "\n", "# Binomial Expansion for a fair coin (p=0.5, q=0.5)\n", "# This gives the Mathematically Expected Binomial Distribution\n", "p = 0.5\n", "q = 0.5\n", "n = num_coins_per_round = 20\n", "binomial_probs = binomial_expansion(p, q, n)\n", "cumulate_prob = 0.0\n", "print(f\"Binomial expansion for (0.5 + 0.5)^{n}:\")\n", "for k, prob in binomial_probs.items():\n", " cumulate_prob += prob\n", " print(f\" {k} heads: {prob:.4f} cumulate_prob: {cumulate_prob} \")\n", "\n", "####### The graphic with rare and normal events distinctly colored.\n", "\n", "fig, ax = plt.subplots()\n", "\n", "k_class = []\n", "counts = []\n", "bar_labels = []\n", "bar_colors = []\n", "\n", "left_tail = 0.0\n", "real_alpha_half = -1.0\n", "\n", "epsilon = 1e-7\n", "for k, prob in binomial_probs.items():\n", " k_class.append(k)\n", " counts.append(prob)\n", " left_tail += prob\n", " if left_tail < 0.25:\n", " bar_labels.append('red')\n", " bar_colors.append('tab:red')\n", " real_alpha_half = left_tail\n", " red_area = 2 * real_alpha_half\n", " print(f\"{k} left_tail = {left_tail}\")\n", " if ((left_tail > 0.25) and (1-left_tail + epsilon > real_alpha_half)) :\n", " bar_labels.append('blue')\n", " bar_colors.append('tab:blue')\n", " print(f\" {k} {1-left_tail + epsilon} > {real_alpha_half} blue\")\n", " if 1-left_tail + epsilon <= real_alpha_half:\n", " bar_labels.append('red')\n", " bar_colors.append('tab:red')\n", " print(f\" {k} {1-left_tail + epsilon} <= {real_alpha_half} red\")\n", "\n", "ax.bar(k_class, counts, label=bar_labels, color=bar_colors)\n", "\n", "ax.set_ylabel('Probability')\n", "ax.set_xlabel('Number of heads')\n", "ax.set_title(f\"Binomial Distribution with rejection regions in red, area = {red_area:.4f}\")\n", "# ax.legend(title='Bar color')\n", "\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": { "id": "SfwmXtymBzJs" }, "source": [ "\n", "> **Challenge.** The previous code functions when the proportion of the rejection region is 0.05. This value is officially called the *significance*, which says that the smaller it is, the rarer an event is, and so it is more suspicious of generating a great discovery. Generalize this code for any significance.\n", "\n", "\n", "# Law or advice?\n", "\n", "The subjectivity of the procedure to divide the world of events into two subsets, corresponding to acceptance else rejection, imposes a conclusion: Science cannot be a law, but its verdicts are only advice to be taken into consideration because a result can be explained in many forms, a fact that is endemic to our theme.\n", "\n", "> **Example.** With our accepted criterion to divide normal from strange, many events can be explained by many different null hypotheses. Let us illustrate this fact by showing that a given event can be normal according to the $H_o$ that says p = 0.5 = 1/2 and also according to the hypothesis that declares $p = 2/3$. Our event is given by a number of heads equal to 14 out of 20. According to the previous table, 14 is a normal value when we consider that the probability of heads is $0.5$. Let us prove now that it is also normal when $p = 2/3$. To that aim, let us simulate this probability.\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "id": "jUPYF4S4e4S1" }, "outputs": [], "source": [ "# Program coinsProbOfHeadsTwoThirds.py\n", "\n", "#Libraries and local rename\n", "import matplotlib.pyplot as plt\n", "import numpy as np #Support for scientific computing\n", "import random\n", "from collections import Counter\n", "\n", "# Parameters for the simulation\n", "num_coins_per_round = 20\n", "num_repetitions = 10000\n", "\n", "# Simulate coin flips\n", "results = []\n", "for _ in range(num_repetitions):\n", " heads_in_round = 0\n", " for _ in range(num_coins_per_round):\n", " # 0 for tails, 1 for heads\n", " if random.choice([0, 1, 1]) == 1:\n", " heads_in_round += 1\n", " results.append(heads_in_round)\n", "\n", "# Calculate statistics\n", "total_heads = sum(results)\n", "proportion_heads = total_heads / (num_coins_per_round * num_repetitions)\n", "heads_distribution = Counter(results)\n", "\n", "# Arrays to keep values to be graphed\n", "x_values = []\n", "y_values = []\n", "\n", "# Print the results\n", "print(\"Heads is 1 in a random choice between 0, 1 \")\n", "print(f\"Simulation of {num_repetitions} rounds \")\n", "print(f\"with {num_coins_per_round} coins per round:\")\n", "print(f\"\\nTotal number of heads across all repetitions: {total_heads}\")\n", "print(f\"Proportion of heads across all repetitions: {proportion_heads:.4f}\")\n", "print(\"\\nDistribution of heads per round:\")\n", "# Sort the distribution for better readability\n", "sorted_distribution = dict(sorted(heads_distribution.items()))\n", "for num_heads, count in sorted_distribution.items():\n", " x_values.append(num_heads)\n", " y_values.append(count)\n", " print(f\" {num_heads} heads: {count} times\")\n", "\n", "# Plot the function\n", "plt.plot(x_values, y_values, label=\"Heads_distribution\")\n", "\n", "# Add labels, title, and legend\n", "plt.xlabel(\"Number of heads\")\n", "plt.ylabel(\"Count\")\n", "plt.title(\"Graph of heads distribution\")\n", "plt.grid(True)\n", "plt.legend()\n", "\n", "# Display the plot\n", "plt.show()\n" ] }, { "cell_type": "markdown", "metadata": { "id": "cg05z8BVgRPe" }, "source": [ "> We see that 14 heads is, with almost the highest frequency, at the mode of the distribution. So, it is normal by our criterion. Thus, the observed datum 14 heads can be explained by two hypotheses, the first p = 0.5, the second p = 2/3 = 0.666...\n", "\n", "Therefore, when Science presents a model that cannot be rejected with the data at hand, it does not guarantee that the model is the truth. It only says that the model is compatible with the data. This is a form of saying that more explanations are welcome. Anyway, when too many theories are compatible with the data at hand, we recognize that the field of investigation is rather crude and more research must be done to single out one winning theory.\n", "\n", "\n", "# The Scientific Algorithm\n", "\n", "Our procedure for studying beliefs is almost universal and can be applied to a wide range of situations. Can we abstract it to get a formal method?\n", "\n", "> **Exercise.** By reflecting on our case study of $H_o$ that says $p = 0.05$, formulate a procedure or algorithm for the Scientific Study of null hypotheses in the face of randomness:\n", "\n", ">> **Answer.** The following is a procedure or algorithm to scientifically study beliefs:\n", ">>>\n", ">>> 1. Observe the world to extract possible laws or tendencies that become a belief.\n", ">>>\n", ">>> 2. If possible, formulate the belief to be studied in an operative language that can be simulated or implemented in a laboratory.\n", ">>>\n", ">>> 3. Let randomness be the first possible explanation of your data and belief. This is the null hypothesis, $H_o$.\n", ">>>\n", ">>> 4. Program and simulate your $H_o$ to estimate the distribution of the observable of interest.\n", ">>>\n", ">>> 5. Divide the world of events into two parts: one that supports the belief or $H_o$, the null hypothesis, and the other that rejects it. When the distribution of the measured characteristic has the form of a bell, the region around the mode and covering 95% is the region of acceptance of the $H_o$. The region covering the tails of the bell is the region of rejection. The default criterion of Science assigns 95% to the acceptance region and 5% to the region of rejection.\n", ">>>\n", ">>> 6. Make a verdict: if the observed event is in the region of acceptance, accept the $H_o$. Otherwise, reject it. An acceptance does not mean that the $H_o$ is the explanation of the data, but only that the data and the belief are compatible. By the same token, a rejection does not mean that data is not explained by the $H_o$ but only that it is best explained by other hypotheses that are necessary to concoct.\n", ">>>\n", ">>> 7. Interpret the verdict.\n", "\n", "> **Example.** Our game with 20 coins represents the situation of traffic in a lane and country when the velocity has a bell distribution around 60 km/h. So, 1 head corresponds with a velocity of 6 kilometers per hour, and 14 heads means a velocity of 84km/h, which in this lane, zone, and country is normal, and the client will not pay a fine.\n", "\n", "\n", "# The elegance of the z-distribution\n", "\n", "Statistics is an ancient inheritance consisting of pure mathematical reasoning that assumes that we can a priori define probabilities whose experimental verification demands perfect random generators and many experiments near to infinity. This is exactly what is required when we say that the probability of heads is 0.5. Now, statistics are very useful in real life because valuable conclusions are stable against noise, inexactness, and approximations.\n", "\n", "A characteristic of the ancient inheritance is its elegance. To witness this, let us look at a characterization of the Scientific Method in the presence of noise that causes data to include random components. Randomness appears when various factors that may be deterministic interact in complex ways, making it difficult to predict what will be next. The characterization of the Scientific Method is the following:\n", "\n", "**Science is all about the measurement of the difference between what one thinks, the null hypothesis, and what one observes in an experiment, in a controlled environment. If that difference is too high, one is invited to conclude that what one thinks is wrong, the null hypothesis would be rejected, and that a new explanation must be concocted. But if the difference is small, one declares that we lack reasons to change our thought, and the null hypothesis is accepted.**\n", "\n", "\n", "In relation to our simulation with coins, if we toss a fair coin 20 times, we expect a number of heads that symmetrically dance around 10, and that gives rise to a bell-like distribution. For instance, if we find 14 heads, the difference from the expected, under the null hypothesis of fairness of the coin, is 4. But we must take one more step forward.\n", "\n", "Let us apply the following observation: All bells have the same silhouette. How do we formalize this observation? We must free the shape from the scales, both horizontal, given by the number of coins, and vertical, given by the relative frequencies. Ancient thinkers found the way to do that: The difference between what is observed and what is expected must be divided by the standard deviation of the bell curve, which we denote as $s$, which is the square root of the variance and serves as a measure of its dispersion.\n", "\n", "We end with a measure of discrepancy $d$ between observed and expected as\n", ">\n", "> $d = \\frac{Observed - Expected }{standardDeviation}$\n", ">\n", "In our case, we measure the discrepancy between the observed number of heads, say 14, and the mean, 10, that is the expected average value of the number of heads when the probability of a head is 0.5, and the number of coins per round is 20. The discrepancy between observed and expected in units of standard deviations is\n", ">\n", "> $z = \\frac{14-10}{s}$\n", ">\n", "To decide whether the observed value, 14, is normal else extreme, we must find the distribution of this measure, decide what is normal and what is extreme, and consequently classify our observed result. The distribution of this measure is called the $z-distribution$. The standard deviation $s$ is calculated from a finite number of observations and changes from one experiment to another. But we must fear nothing because the caused error is not mortal, the whole procedure is stable against this variation. In any case, in our simulations we take many events, around 10000, and we hope to be very near to mathematical expectations.\n", "\n", "To advance with our plan, we need the code for the standard deviation of a given dataset. \n", "\n", "> **Example.** We prompted Gemini with:\n", ">\n", ">> Please, help me with the Python code to calculate the standard deviation of an array of numbers called 'y_values'.\n", ">\n", "> In response, Gemini displayed the following code:\n", ">\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true, "id": "eZHIuDvDHCyC" }, "outputs": [], "source": [ "# program standardDeviation.py\n", "# By Gemini\n", "\n", "import statistics\n", "\n", "# Example, array with numbers\n", "y_values = [10, 12, 23, 23, 16, 23, 21, 16]\n", "\n", "# Berechne die Standar Deviation\n", "std_dev = statistics.stdev(y_values)\n", "\n", "# Print results\n", "print(f\"The array of numbers is: {y_values}\")\n", "print(f\"The standard deviation of y_values is: {std_dev:.4f}\")\n" ] }, { "cell_type": "markdown", "metadata": { "id": "XBXrWpLBI4oG" }, "source": [ "> Gemini added the result of the computation\n", "\n", " >Results:\n", ">\n", ">> The array of numbers is: : [10, 12, 23, 23, 16, 23, 21, 16]\n", ">>\n", ">> The standard deviation of y_values is: 5.2372\n", ">>\n", "> We see that Python has a library with the procedure to calculate the standard deviation of an array of data, and of course, the mean. We will apply this to the array of results in our simulations. Warning: the literature distinguishes two kinds of standard deviations, the first for a population when every object in the Universal Set of observation is censed. The second is for a sample, when a random choice of objects is measured. When the number of data is large, greater than 100, the difference between the two is irrelevant. So, we will not pay attention to it.\n", "\n", "\n", "> **Example.** Let us estimate the $z-distribution$ with a simulation." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "id": "OQwYvF-B_9v8" }, "outputs": [], "source": [ "# Program zDistribution.py\n", "# The Binomial Distribution is estimated by a simulation and drawn,\n", "# it is then standardized to produce an estimation of\n", "# the z-distribution.\n", "\n", "\n", "# Libraries\n", "import matplotlib.pyplot as plt\n", "import numpy as np #Support for scientific computing\n", "import random\n", "from collections import Counter\n", "import statistics\n", "\n", "# Parameters for the simulation\n", "num_coins_per_round = 20\n", "num_repetitions = 10000\n", "\n", "# Simulate coin flips\n", "results = []\n", "for _ in range(num_repetitions):\n", " heads_in_round = 0\n", " for _ in range(num_coins_per_round):\n", " # 0 for tails, 1 for heads\n", " if random.choice([0, 1]) == 1:\n", " heads_in_round += 1\n", " results.append(heads_in_round)\n", "\n", "# Calculate statistics\n", "total_heads = 0\n", "for x in results:\n", " total_heads += x\n", "proportion_heads = total_heads / (num_coins_per_round * num_repetitions)\n", "heads_distribution = Counter(results)\n", "\n", "\n", "#Arrays to keep values to be graphed\n", "x_values = []\n", "y_values = []\n", "\n", "\n", "\n", "# Print the results\n", "print(\"Heads corresponds with 1 in a random choice between 0, 1 \")\n", "print(f\"Simulation of {num_repetitions} rounds \")\n", "print(f\"with {num_coins_per_round} coins per round:\")\n", "print(f\"\\nTotal number of heads across all repetitions: {total_heads}\")\n", "print(f\"Proportion of heads across all repetitions: {proportion_heads:.4f}\")\n", "print(\"\\nDistribution of heads per round:\")\n", "# Sort the distribution for better readability\n", "sorted_distribution = dict(sorted(heads_distribution.items()))\n", "for num_heads, count in sorted_distribution.items():\n", " x_values.append(num_heads)\n", " y_values.append(count)\n", " print(f\" {num_heads} heads: {count} times\")\n", "\n", "# Plot the function\n", "plt.plot(x_values, y_values, label=\"Heads_distribution\")\n", "\n", "# Add labels, title, and legend\n", "plt.xlabel(\"Number of heads\")\n", "plt.ylabel(\"Count\")\n", "plt.title(\"Graph of a heads distribution\")\n", "plt.grid(True)\n", "plt.legend()\n", "\n", "# Display the plot\n", "plt.show()\n", "\n", "# values for the z-distribution, z_values and estimated probability\n", "z_values = []\n", "p_values = []\n", "\n", "# the mean number of heads is computed\n", "mean = statistics.mean(results)\n", "mean2 = 0.5*num_coins_per_round\n", "print(f\"Observed mean = {mean}, expected {mean2}\")\n", "# The standard deviation of results is computed\n", "std_dev = statistics.stdev(results)\n", "std_dev_expected = np.sqrt(mean2*0.5) #sqrt(npq)\n", "print(f\"Observed Standard Deviation: {std_dev:.4f}, Expected: {std_dev_expected:.4f}\")\n", "\n", "# The z-distribution\n", "# z-transform to measure the difference between observed and expected\n", "# Observed = num_heads\n", "# Expected = mean2 = np = 20 * 0.5 = 10\n", "total_frequency = 0\n", "for num_heads, count in sorted_distribution.items():\n", " z = (num_heads - mean2) / std_dev # Calculate z-score using the observed mean\n", " z_values.append(z)\n", " p = count/num_repetitions\n", " p_values.append(p)\n", " total_frequency += p\n", " print(f\" {z:.4f} z_value: {p:.4f} probability\")\n", "\n", "\n", "print(f\"Total frequency= {total_frequency}, expected = 1\")\n", "\n", "\n", "# Plot the estimated z-distribution\n", "plt.plot(z_values, p_values, label=\"Estimated z-distribution\")\n", "\n", "# Add labels, title, and legend\n", "plt.xlabel(\"z = standarized discrepancy between observed and expected\")\n", "plt.ylabel(\"Probability\")\n", "plt.title(\"Estimated z-distribution\")\n", "plt.grid(True)\n", "plt.legend()\n", "\n", "\n", "# Display the plot\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": { "id": "1_Jp6Q0j4MxX" }, "source": [ "> **Partial conclusion.** The Binomial distribution is discrete. Nevertheless, we have frequently fit its envelope to a curve that resembles a bell, whose technical name is Gauss-Bell, which is continuous. These bells can be more or less spread, but common sense says that all Gauss-Bells are the same, in the sense that if we measure the events in units of standard deviations, we cannot distinguish one from another. This is made rigorous by the z-transform, which is the name of the $z$-formula. \n", "\n", "> **Exercise.** Develop the code that shows accepting and rejection colored regions for the z-transform generated by the Binomial distribution with $p = 0.5$ and for whatever number of coins.\n", ">\n", ">> **Answer** We asked Colab to devise the corresponding code with the following prompt:\n", ">>\n", ">> Code for the graph of the z-transform showing colored acceptance and rejection regions with confidence 0.05\n", ">>\n", ">> Gemini answered:\n", ">>> *Okay, I will add the code to the selected cell to graph the z-transform of the Binomial distribution, highlighting the acceptance and rejection regions.*\n", ">>\n", ">> And it posted the following code:" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "id": "oXfmKsmZj10A" }, "outputs": [], "source": [ "# Program zDistributionColoredTails.py\n", "# This program calculates and plots the z-transform of the Binomial distribution\n", "# with acceptance and rejection regions highlighted for a given significance level.\n", "# By Gemini\n", "\n", "# Libraries and local rename\n", "import math\n", "import matplotlib.pyplot as plt\n", "import numpy as np # Support for scientific computing\n", "from scipy.stats import binom # Import the binomial distribution from scipy\n", "\n", "# The binomial expansion\n", "def binomial_expansion(p, q, n):\n", " \"\"\"\n", " Calculates the terms of the binomial expansion (p+q)^n.\n", "\n", " Args:\n", " p: The probability of the first outcome.\n", " q: The probability of the second outcome (q = 1 - p).\n", " n: The number of trials.\n", "\n", " Returns:\n", " A dictionary where keys are the number of successes (k)\n", " and values are the corresponding probabilities.\n", " \"\"\"\n", " if not math.isclose(p + q, 1.0):\n", " print(\"Warning: p + q should be equal to 1.\")\n", "\n", " expansion_terms = {}\n", " for k in range(n + 1):\n", " # Calculate the binomial coefficient (n choose k)\n", " binomial_coeff = math.comb(n, k)\n", "\n", " # Calculate the probability of k successes\n", " probability_k = binomial_coeff * (p ** k) * (q ** (n - k))\n", "\n", " expansion_terms[k] = probability_k\n", " return expansion_terms\n", "\n", "# Function to calculate and plot the z-transform with colored regions\n", "def plot_z_distribution_colored(num_coins_per_round, alpha=0.05):\n", " \"\"\"\n", " Calculates and plots the z-transform of the Binomial distribution\n", " with acceptance and rejection regions highlighted.\n", "\n", " Args:\n", " num_coins_per_round: The number of coin flips (n).\n", " alpha: The significance level (proportion of the rejection region).\n", " Defaults to 0.05.\n", " \"\"\"\n", " p = 0.5\n", " n = num_coins_per_round\n", "\n", " # Calculate binomial probabilities using scipy\n", " k_values = np.arange(0, n + 1)\n", " binomial_probs = binom.pmf(k_values, n, p)\n", "\n", " # Calculate the mean and standard deviation of the binomial distribution\n", " mean = n * p\n", " std_dev = np.sqrt(n * p * (1 - p))\n", "\n", " # Calculate z-scores\n", " z_values = (k_values - mean) / std_dev\n", "\n", " # Determine the rejection regions based on the cumulative distribution\n", " cumulative_probs = np.cumsum(binomial_probs)\n", " lower_tail_limit = alpha / 2\n", " upper_tail_limit = 1 - (alpha / 2)\n", "\n", " # Find the k values that define the boundaries of the rejection regions\n", " lower_rejection_k = np.where(cumulative_probs > lower_tail_limit)[0][0] if np.where(cumulative_probs > lower_tail_limit)[0].size > 0 else 0\n", " upper_rejection_k = np.where(cumulative_probs >= upper_tail_limit)[0][0] if np.where(cumulative_probs >= upper_tail_limit)[0].size > 0 else n + 1\n", "\n", " # Prepare data for plotting\n", " bar_colors = []\n", " for k in k_values:\n", " if k < lower_rejection_k or k >= upper_rejection_k:\n", " bar_colors.append('tab:red') # Rejection region\n", " else:\n", " bar_colors.append('tab:blue') # Acceptance region\n", "\n", " # Calculate the actual area of the rejection regions\n", " actual_alpha = np.sum(binomial_probs[k_values < lower_rejection_k]) + np.sum(binomial_probs[k_values >= upper_rejection_k])\n", "\n", "\n", " # Plot the z-distribution with colored regions\n", " fig, ax = plt.subplots(figsize=(10, 6))\n", " ax.bar(z_values, binomial_probs, color=bar_colors, width=z_values[1] - z_values[0] if len(z_values) > 1 else 0.8) # Adjust width for better visualization\n", "\n", " ax.set_ylabel('Probability')\n", " ax.set_xlabel('Z-score (Standardized difference from expected)')\n", " ax.set_title(f'Z-transform of Binomial Distribution (n={n}, p={p}) with Rejection Regions (alpha={alpha}, actual alpha={actual_alpha:.4f})')\n", " ax.grid(True)\n", "\n", "\n", " plt.show()\n", "\n", "# Example usage:\n", "plot_z_distribution_colored(num_coins_per_round=20, alpha=0.05)\n", "plot_z_distribution_colored(num_coins_per_round=10, alpha=0.05)\n", "plot_z_distribution_colored(num_coins_per_round=5, alpha=0.10)" ] }, { "cell_type": "markdown", "metadata": { "id": "dalZYA--jtKQ" }, "source": [ "> **Exercise.** Run the code and evaluate it.\n", ">>\n", ">> **Answer.** The code developed by Gemini shows the z-transform derived from symmetric Binomial Distributions. Resultant graphs are not symmetric.\n", "\n", "Let us add to the z-transform the continuous envelope of the Binomial Distribution. The resultant function is known as the standardized Gauss Bell. So, we prompted Colab with the prompt:\n", ">\n", "> The code that draws the standardized Gauss Bell with rejection regions with significance 0.05.\n", "\n", "Colab answered with the following program:" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "id": "bOVF-WtKmlpb" }, "outputs": [], "source": [ "# Program StandardizedGaussBellColoredTails.py\n", "# This program draws the standardized Gauss Bell (normal distribution)\n", "# with colored acceptance and rejection regions for a given significance level.\n", "# By Gemini\n", "\n", "# Libraries and local rename\n", "import matplotlib.pyplot as plt\n", "import numpy as np # Support for scientific computing\n", "from scipy.stats import norm # Import the normal distribution from scipy\n", "\n", "# Function to plot the standardized Gauss Bell with colored regions\n", "def plot_standard_normal_colored(alpha=0.05):\n", " \"\"\"\n", " Plots the standardized normal distribution with acceptance and rejection regions highlighted.\n", "\n", " Args:\n", " alpha: The significance level (proportion of the rejection region).\n", " Defaults to 0.05.\n", " \"\"\"\n", " # Generate x-values for the plot\n", " x_values = np.linspace(-4, 4, 200)\n", "\n", " # Calculate the probability density function (PDF) for the standard normal distribution\n", " pdf_values = norm.pdf(x_values, 0, 1) # Mean = 0, Standard Deviation = 1\n", "\n", " # Determine the z-scores that define the boundaries of the rejection regions\n", " # Use the percent point function (inverse of CDF) to find the z-score for the given cumulative probability\n", " lower_z = norm.ppf(alpha / 2)\n", " upper_z = norm.ppf(1 - (alpha / 2))\n", "\n", " # Prepare data for plotting - separate into acceptance and rejection regions\n", " x_acceptance = x_values[(x_values >= lower_z) & (x_values <= upper_z)]\n", " pdf_acceptance = pdf_values[(x_values >= lower_z) & (x_values <= upper_z)]\n", "\n", " x_lower_rejection = x_values[x_values < lower_z]\n", " pdf_lower_rejection = pdf_values[x_values < lower_z]\n", "\n", " x_upper_rejection = x_values[x_values > upper_z]\n", " pdf_upper_rejection = pdf_values[x_values > upper_z]\n", "\n", "\n", " # Plot the standardized normal distribution with colored regions\n", " fig, ax = plt.subplots(figsize=(10, 6))\n", "\n", " # Plot the acceptance region\n", " ax.fill_between(x_acceptance, pdf_acceptance, color='tab:blue', alpha=0.7, label='Acceptance Region')\n", "\n", " # Plot the rejection regions\n", " ax.fill_between(x_lower_rejection, pdf_lower_rejection, color='tab:red', alpha=0.7, label='Rejection Region')\n", " ax.fill_between(x_upper_rejection, pdf_upper_rejection, color='tab:red', alpha=0.7) # No label for the second tail\n", "\n", " ax.set_ylabel('Probability Density')\n", " ax.set_xlabel('Z-score')\n", " ax.set_title(f'Standardized Normal Distribution with Rejection Regions (alpha={alpha})')\n", " ax.grid(True)\n", " ax.legend()\n", "\n", " plt.show()\n", "\n", "# Example usage:\n", "plot_standard_normal_colored(alpha=0.05)\n", "plot_standard_normal_colored(alpha=0.01)" ] }, { "cell_type": "markdown", "metadata": { "id": "Y4qVBl-LmfUW" }, "source": [ "> **Exercise.** Run the code and assess it.\n", ">>\n", ">> **Answer.** The code is astonishing. It displays the power of Python that renders things easy, together with a transparent mathematical reasoning. The virtue of this code is that it is extremely general. It presents two examples. The first is with significance 0.05, and the right rejection region begins a bit before 2. This matches the literature: the z-value is $1.96$. So, we consider that the code is correct. The code approximates all symmetric Binomial Distributions that, after the z-transform, are enveloped by this Gaussian Bell. What we need to know is how to pass from a Binomial distribution to this graph. This is done by calculating the mean and the standard deviation. If the Binomial Distribution is defined by $n$ coins and with the probability of heads equal to $p = 0.5$, the mean number of heads is $\\mu = np = n/2$ and the standard deviation is $\\sigma = \\sqrt{np(1-p)} = \\sqrt{n/4}$. With this, the z-transform becomes:\n", ">>\n", ">>> $z(x) = \\frac{x-\\mu}{\\sigma} = \\frac{x-n/2}{\\sqrt{n/4}} = \\frac{x-n/2}{\\sqrt{n}/2} = \\frac{2x-n}{\\sqrt{n}}$\n", ">>\n", ">> For instance, if we have 20 coins, and a number of heads $x = 12$, the z-value becomes: $z(12) = \\frac{2(12) - 20}{\\sqrt{20}} = \\frac{4}{4.5} = 0.9 \\sim 1$\n", ">>\n", ">> This means that 12 is 1 standard deviation upwards from the mean, which is 10. Is that true? Let us see:\n", ">>\n", ">> The standard deviation is\n", ">>\n", ">>> $\\sigma = \\sqrt{npq} = \\sqrt{np(1-p)} = \\sqrt{n/4} = \\sqrt{20/4} = \\sqrt{5} = 2.2$\n", ">>\n", ">> Since the mean is $10$, the mean plus 1 standard deviation is $10 + 2.2 = 12.2 \\sim 12$, such as desired.\n", "\n", "> **Exercise** Superpose a Standard Gauss Bell over the envelope of standardized Binomial Distribution to see how well the fit is.\n", ">>\n", ">> **Answer.** We need to draw two functions on the same graph. So, we prompted Google with:\n", ">>\n", ">>> how to draw two functions in python\n", ">>\n", ">> Gemini answered with an explained procedure, from which we extracted the following code:\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "id": "OuF-TA9SsORl" }, "outputs": [], "source": [ "# Program DrawTwoFunctions.py\n", "# By Gemini\n", "\n", "import matplotlib.pyplot as plt\n", "import numpy as np\n", "\n", "# Functions are defined\n", "def func1(x):\n", " return np.sin(x)\n", "\n", "def func2(x):\n", " return np.cos(x)\n", "\n", "# 400 points along the x-axis are chosen\n", "x_values = np.linspace(0, 2 * np.pi, 400) # From 0 to 2pi, with 400 points\n", "# The y-values f the two functions are calculated\n", "y_values_func1 = func1(x_values)\n", "y_values_func2 = func2(x_values)\n", "# The two functions are drawn\n", "plt.plot(x_values, y_values_func1, label='sin(x)', color='blue', linestyle='-')\n", "plt.plot(x_values, y_values_func2, label='cos(x)', color='red', linestyle='--')\n", "plt.xlabel('X-axis')\n", "plt.ylabel('Y-axis')\n", "plt.title('Plot of Two Functions')\n", "plt.legend() # Displays the labels defined in plt.plot()\n", "\n", "# Display the Plot.\n", "\n", "plt.grid(True) # Optional: adds a grid\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": { "id": "UzSUIlGcsEDu" }, "source": [ "\n", "\n", "> The lesson is: to draw two functions on the same graph, x-values are fabricated and used for the two y-values. Each function is drawn separately. To implement this, we reused programs `zDistribution.py`, which shows the envelope of a standardized Binomial Distribution, and `Program StandardizedGaussBellColoredTails.py` that depicts the Standard Bell. We encountered a problem: the x-values are chosen differently for each function. We tried out the following code:\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "id": "PLc2R1Rin4j-" }, "outputs": [], "source": [ "# Program standardBellFitsBinomialDistribution.py\n", "# The Binomial Distribution is estimated by a simulation and drawn,\n", "# it is then standardized to produce an estimation of\n", "# the z-distribution.\n", "# A Standard Gauss Bell is superposed to check the fitting.\n", "\n", "\n", "# Libraries\n", "import matplotlib.pyplot as plt\n", "import numpy as np #Support for scientific computing\n", "import random\n", "from collections import Counter\n", "import statistics\n", "from scipy.stats import norm # Import the normal distribution from scipy\n", "\n", "# Function to plot the standardized Gauss Bell\n", "def values_standard_normal():\n", " \"\"\"\n", " Plots the standardized normal distribution with acceptance and rejection regions highlighted.\n", "\n", " Args:\n", " alpha: The significance level (proportion of the rejection region).\n", " Defaults to 0.05.\n", " \"\"\"\n", " # Generate x-values for the plot\n", " x_values = np.linspace(-4, 4, 200)\n", "\n", " # Calculate the probability density function (PDF) for the standard normal distribution\n", " pdf_values = norm.pdf(x_values, 0, 1) # Mean = 0, Standard Deviation = 1\n", "\n", " return x_values, pdf_values\n", "\n", "\n", "# Parameters for the simulation of the Binomial\n", "num_coins_per_round = 20\n", "num_repetitions = 10000\n", "\n", "# Simulate coin flips\n", "results = []\n", "for _ in range(num_repetitions):\n", " heads_in_round = 0\n", " for _ in range(num_coins_per_round):\n", " # 0 for tails, 1 for heads\n", " if random.choice([0, 1]) == 1:\n", " heads_in_round += 1\n", " results.append(heads_in_round)\n", "\n", "# Calculate statistics\n", "total_heads = 0\n", "for x in results:\n", " total_heads += x\n", "proportion_heads = total_heads / (num_coins_per_round * num_repetitions)\n", "heads_distribution = Counter(results)\n", "\n", "\n", "#Arrays to keep values to be graphed\n", "x_values = []\n", "y_values = []\n", "\n", "\n", "\n", "# Print the results\n", "print(\"Heads corresponds with 1 in a random choice between 0, 1 \")\n", "print(f\"Simulation of {num_repetitions} rounds \")\n", "print(f\"with {num_coins_per_round} coins per round:\")\n", "print(f\"\\nTotal number of heads across all repetitions: {total_heads}\")\n", "print(f\"Proportion of heads across all repetitions: {proportion_heads:.4f}\")\n", "print(\"\\nDistribution of heads per round:\")\n", "# Sort the distribution for better readability\n", "sorted_distribution = dict(sorted(heads_distribution.items()))\n", "for num_heads, count in sorted_distribution.items():\n", " x_values.append(num_heads)\n", " y_values.append(count)\n", " print(f\" {num_heads} heads: {count} times\")\n", "\n", "# Plot the function\n", "plt.plot(x_values, y_values, label=\"Heads_distribution\")\n", "\n", "# Add labels, title, and legend\n", "plt.xlabel(\"Number of heads\")\n", "plt.ylabel(\"Count\")\n", "plt.title(\"Graph of a heads distribution\")\n", "plt.grid(True)\n", "plt.legend()\n", "\n", "# Display the plot\n", "plt.show()\n", "\n", "\n", "\n", "# Standarization\n", "\n", "# values for the z-distribution, z_values and estimated probability\n", "z_values = []\n", "p_values = []\n", "\n", "# the mean number of heads is computed\n", "mean = statistics.mean(results)\n", "mean2 = 0.5*num_coins_per_round\n", "print(f\"Observed mean = {mean}, expected {mean2}\")\n", "# The standard deviation of results is computed\n", "std_dev = statistics.stdev(results)\n", "std_dev_expected = np.sqrt(mean2*0.5) #sqrt(npq)\n", "print(f\"Observed Standard Deviation: {std_dev:.4f}, Expected: {std_dev_expected:.4f}\")\n", "\n", "# The z-distribution\n", "# z-transform to measure the difference between observed and expected\n", "# Observed = num_heads\n", "# Expected = mean2 = np = 20 * 0.5 = 10\n", "total_frequency = 0\n", "for num_heads, count in sorted_distribution.items():\n", " z = (num_heads - mean2) / std_dev # Calculate z-score using the observed mean\n", " z_values.append(z)\n", " p = count/num_repetitions\n", " p_values.append(p)\n", " total_frequency += p\n", " print(f\" {z:.4f} z_value: {p:.4f} probability\")\n", "\n", "\n", "print(f\"Total frequency= {total_frequency}, expected = 1\")\n", "\n", "# Plot the estimated z-distribution\n", "plt.plot(z_values, p_values, label=\"Estimated z-distribution\")\n", "\n", "x_Bell, y_Bell = values_standard_normal()\n", "# Plot the Gauss Bell\n", "plt.plot(x_Bell, y_Bell, label=\"Standard Gauss Bell\")\n", "\n", "# Add labels, title, and legend\n", "plt.xlabel(\"z_values\")\n", "plt.ylabel(\"Probability\")\n", "plt.title(\"Standard Gauss Bell fits the Binomial\")\n", "plt.grid(True)\n", "plt.legend()\n", "\n", "# Display the plot\n", "plt.show()\n" ] }, { "cell_type": "markdown", "metadata": { "id": "AwZiVNYUmmxb" }, "source": [ "> There was no problem with the different x-values for each function. On the other hand, the resultant functions both have a bell shape, but their heights differ. Why? We see only one cause: the normalization procedures. In fact, the Gauss Bell is by its very nature normalized to represent a probability density function, normalized by area, whose integral is 1. But our estimation of the Binomial from a simulation has been normalized by the height of bars, whose total length must be 1. To test this belief, we have a prediction: there must be a constant that serves as a stretching factor that renders both graphs almost equal. The constant relates the heights of both functions at their maxima. The code follows:\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "id": "bFfxiFuu0xKz" }, "outputs": [], "source": [ "# Program standardBellFitsBinomialDistribution.py\n", "# The Binomial Distribution is estimated by a simulation and drawn,\n", "# it is then standardized to produce an estimation of\n", "# the z-distribution.\n", "# A Standard Gauss Bell is superposed to check the fitting.\n", "\n", "\n", "# Libraries\n", "import matplotlib.pyplot as plt\n", "import numpy as np #Support for scientific computing\n", "import random\n", "from collections import Counter\n", "import statistics\n", "from scipy.stats import norm # Import the normal distribution from scipy\n", "\n", "\n", "##########\n", "\n", "# Function to plot the standardized Gauss Bell with colored regions\n", "def values_standard_normal():\n", " \"\"\"\n", " Plots the standardized normal distribution with acceptance and rejection regions highlighted.\n", "\n", " Args:\n", " alpha: The significance level (proportion of the rejection region).\n", " Defaults to 0.05.\n", " \"\"\"\n", " # Generate x-values for the plot\n", " x_values = np.linspace(-4, 4, 400)\n", "\n", " # Calculate the probability density function (PDF) for the standard normal distribution\n", " pdf_values = norm.pdf(x_values, 0, 1) # Mean = 0, Standard Deviation = 1\n", "\n", " return x_values, pdf_values\n", "\n", "#########\n", "\n", "# Parameters for the simulation of the Binomial\n", "num_coins_per_round = 20\n", "num_repetitions = 10000\n", "\n", "# Simulate coin flips and count heads\n", "results = []\n", "for _ in range(num_repetitions):\n", " heads_in_round = 0\n", " for _ in range(num_coins_per_round):\n", " # 0 for tails, 1 for heads\n", " if random.choice([0, 1]) == 1:\n", " heads_in_round += 1\n", " results.append(heads_in_round)\n", "\n", "# Calculate statistics\n", "total_heads = 0\n", "for x in results:\n", " total_heads += x\n", "total_number_coins = num_coins_per_round * num_repetitions\n", "proportion_heads = total_heads / (total_number_coins)\n", "heads_distribution = Counter(results)\n", "\n", "\n", "# Arrays to keep values to be graphed\n", "x_values = []\n", "y_values = []\n", "\n", "# Print results\n", "print(\"Heads corresponds with 1 in a random choice between 0, 1 \")\n", "print(f\"Simulation of {num_repetitions} rounds \")\n", "print(f\"with {num_coins_per_round} coins per round:\")\n", "print(f\"\\nTotal number of heads across all repetitions: {total_heads}\")\n", "print(f\"Proportion of heads across all repetitions: {proportion_heads:.4f}\")\n", "print(\"\\nDistribution of heads per round:\")\n", "\n", "# Sort the distribution according to human needs\n", "sorted_distribution = dict(sorted(heads_distribution.items()))\n", "for num_heads, count in sorted_distribution.items():\n", " x_values.append(num_heads)\n", " y_values.append(count)\n", " print(f\" {num_heads} heads: {count} times\")\n", "\n", "#Normalization of y_values\n", "\n", "# y_values = [y / total_number_coins for y in y_values]\n", "y_values = [y / num_repetitions for y in y_values]\n", "# Plot the function\n", "plt.plot(x_values, y_values, label=\"Heads_distribution\")\n", "\n", "# Add labels, title, and legend\n", "plt.xlabel(\"Number of heads\")\n", "plt.ylabel(\"Count\")\n", "plt.title(\"Graph of a heads distribution\")\n", "plt.grid(True)\n", "plt.legend()\n", "\n", "# Display the plot\n", "plt.show()\n", "\n", "\n", "\n", "# Standarization\n", "\n", "# values for the z-distribution, z_values and estimated probability\n", "z_values = []\n", "p_values = []\n", "\n", "# the mean number of heads is computed\n", "mean = statistics.mean(results)\n", "mean2 = 0.5*num_coins_per_round\n", "print(f\"Observed mean = {mean}, expected {mean2}\")\n", "# The standard deviation of results is computed\n", "std_dev = statistics.stdev(results)\n", "std_dev_expected = np.sqrt(mean2*0.5) #sqrt(npq)\n", "print(f\"Observed Standard Deviation: {std_dev:.4f}, Expected: {std_dev_expected:.4f}\")\n", "\n", "# The z-distribution\n", "# z-transform to measure the difference between observed and expected\n", "# Observed = num_heads\n", "# Expected = mean2 = np = 20 * 0.5 = 10\n", "total_frequency = 0\n", "for num_heads, count in sorted_distribution.items():\n", " z = (num_heads - mean2) / std_dev # Calculate z-score using the observed mean\n", " z_values.append(z)\n", " p = count/num_repetitions\n", " p_values.append(p)\n", " total_frequency += p\n", " print(f\" {z:.4f} z_value: {p:.4f} relative frequency\")\n", "\n", "\n", "print(f\"Total frequency= {total_frequency}, expected = 1\")\n", "\n", "# values of the Standard Bell\n", "x_Bell, y_Bell = values_standard_normal()\n", "print(\"y_values of the bell \" , y_Bell)\n", "\n", "# Reescaling\n", "bell_0 = y_Bell[200] #max value of the bell with 200 points\n", "print(f\"bell_0 = {bell_0}\")\n", "print(\"p_values : \" , p_values)\n", "\n", "# Maximal value of the Binomial\n", "max = -1\n", "for i in range(len(p_values)):\n", " if p_values[i] > max:\n", " max = p_values[i]\n", "\n", "print(f\"max = {max}\")\n", "\n", "# Rescaling to adjust maxima\n", "k = bell_0/max\n", "print(f\"Streching factor, k = {k}\")\n", "for i in range(len(p_values)):\n", " p_values[i] = p_values[i]*k\n", "\n", "\n", "# Plot the estimated z-distribution\n", "plt.plot(z_values, p_values, label=\"Estimated z-distribution\")\n", "\n", "\n", "# Plot the Gauss Bell\n", "plt.plot(x_Bell, y_Bell, label=\"Standard Gauss Bell\")\n", "\n", "# Add labels, title, and legend\n", "plt.xlabel(\"z_values\")\n", "plt.ylabel(\"Probability\")\n", "plt.title(\"Standard Gauss Bell fits the Standardized Binomial\")\n", "plt.grid(True)\n", "plt.legend()\n", "\n", "# Display the plot\n", "plt.show()\n" ] }, { "cell_type": "markdown", "metadata": { "id": "b9ahgWPNOEhO" }, "source": [ "> After stretching the binomial distribution by an appropriate constant, we find that the graphs of both functions match one another. This result supports our belief that the different modes of normalization cause different maximal heights for each function.\n", "\n", "Let us test our procedures now with another setting.\n", "\n", "> **Exercise.** Superpose a Standard Gauss Bell over the bar diagram of a Binomial Distribution to see how well the fit is.\n", ">>\n", ">> **Answer.** The mean of the Binomial Distribution is $np$ and the deviation is $\\sqrt{np(1-p)}$. When $p = 0.5$ these become $n/2$ and $\\sqrt{n}/2$ respectively. The code of the bar plot of the Binomial was adapted from program `binomialDistriColoredTails.py` and that of the Bell from `binomialDistriColoredTails.py`. The procedure is a variation of the code in the answer to the last exercise. The code follows:\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "id": "IhBNe6BTOMNt" }, "outputs": [], "source": [ "# Program GaussBellFitsBinomialDistribution.py\n", "# The Binomial Distribution comes from the Newton Binomial.\n", "# A Gauss Bell is superposed to check the fitting.\n", "\n", "\n", "# Libraries\n", "import math\n", "import matplotlib.pyplot as plt\n", "import numpy as np #Support for scientific computing\n", "import random\n", "from collections import Counter\n", "import statistics\n", "from scipy.stats import norm # Import the normal distribution from scipy\n", "\n", "###### THE GAUSSIAN BELL ######\n", "\n", "# Function to plot a Gauss Bell with mean m and standard deviation s\n", "# in the range from -1 to n+1, where\n", "# n is the number of coins of a Binomial Distribution.\n", "def values_normal_distribution(m,s, n):\n", " \"\"\"\n", " Plots the standardized normal distribution with acceptance and rejection regions highlighted.\n", "\n", " Args:\n", " alpha: The significance level (proportion of the rejection region).\n", " Defaults to 0.05.\n", " \"\"\"\n", " # Generate x-values for the plot between 0 and n\n", " x_values = np.linspace(0, n, 400) #400 points\n", "\n", " # Calculate the probability density function (PDF) for the standard normal distribution\n", " pdf_values = norm.pdf(x_values, m, s) # Mean = m, Standard Deviation = s\n", "\n", " return x_values, pdf_values\n", "\n", "######## THE BINOMIAL DISTRIBUTION #######\n", "\n", "#The binomial expansion\n", "def binomial_expansion(p, q, n):\n", " \"\"\"\n", " Calculates the terms of the binomial expansion (p+q)^n.\n", "\n", " Args:\n", " p: The probability of the first outcome.\n", " q: The probability of the second outcome (q = 1 - p).\n", " n: The number of trials.\n", "\n", " Returns:\n", " A dictionary where keys are the number of successes (k)\n", " and values are the corresponding probabilities.\n", " \"\"\"\n", " if not math.isclose(p + q, 1.0):\n", " print(\"Warning: p + q should be equal to 1.\")\n", "\n", " expansion_terms = {}\n", " for k in range(n + 1):\n", " # Calculate the binomial coefficient (n choose k)\n", " binomial_coeff = math.comb(n, k)\n", "\n", " # Calculate the probability of k successes\n", " probability_k = binomial_coeff * (p ** k) * (q ** (n - k))\n", "\n", " expansion_terms[k] = probability_k\n", " return expansion_terms\n", "\n", "\n", "# Binomial Expansion for a fair coin (p=0.5, q=0.5)\n", "# This gives the Mathematically Expected Binomial Distribution\n", "p = 0.5\n", "q = 0.5\n", "n = num_coins_per_round = 20\n", "binomial_probs = binomial_expansion(p, q, n)\n", "cumulate_prob = 0.0\n", "print(f\"Binomial expansion for (0.5 + 0.5)^{n}:\")\n", "for k, prob in binomial_probs.items():\n", " cumulate_prob += prob\n", " print(f\" {k} heads: {prob:.4f} cumulate_prob: {cumulate_prob} \")\n", "\n", "####### The bar diagram of the Binomial\n", "\n", "fig, ax = plt.subplots()\n", "\n", "k_class = [] #Horizontal label\n", "p_values = [] #vertical values\n", "bar_labels = [] #colors of bars\n", "bar_colors = []\n", "\n", "for k, prob in binomial_probs.items():\n", " k_class.append(k)\n", " p_values.append(prob)\n", "\n", "ax.bar(k_class, p_values)\n", "ax.set_ylabel('probability')\n", "ax.set_title('Binomial Distribution')\n", "\n", "\n", "############## THE BELL FITS A BINOMIAL ######\n", "\n", "n = 20 # number of coins\n", "p = 0.5 # probability of heads\n", "q = 1 - p\n", "m = n*p # mean\n", "s = math.sqrt(n*p*q) # standard deviation\n", "print(f\"Mean = {m}, Standard Deviation = {s}\")\n", "\n", "# values of the Standard Bell\n", "x_Bell, y_Bell = values_normal_distribution(m,s, n)\n", "\n", "\n", "# Reescaling\n", "\n", "print(\"y_values of the bell \", y_Bell)\n", "\n", "# max value of the bell\n", "bell_0 = y_Bell[200] #max value of the bell with 400 points\n", "print(f\"max value f the bell = bell_0 = {bell_0}\")\n", "\n", "# Maximal value of the Binomial\n", "max = -1\n", "for i in range(len(p_values)):\n", " if p_values[i] > max:\n", " max = p_values[i]\n", "print(f\"max value of the binomial = {max}\")\n", "\n", "# Rescaling to adjust maxima\n", "k = bell_0/max\n", "print(f\"Streching factor, k = {k}\")\n", "for i in range(len(p_values)):\n", " p_values[i] = p_values[i]*k\n", "\n", "\n", "# Envelope of the Binomial\n", "plt.plot(k_class, p_values, label=\"Binomial Distribution\")\n", "\n", "# Plot the Gauss Bell\n", "plt.plot(x_Bell, y_Bell, label=\"Standard Gauss Bell\")\n", "\n", "# Add labels, title, and legend\n", "plt.xlabel(\"z_values\")\n", "plt.ylabel(\"Probability\")\n", "plt.title(\"Standard Gauss Bell fits the Binomial Distribution\")\n", "plt.grid(True)\n", "plt.legend()\n", "\n", "# Display the plot\n", "plt.show()\n" ] }, { "cell_type": "markdown", "metadata": { "id": "SUf3ZHe-ONCI" }, "source": [ "## Interpretation\n", "\n", "To compare the graph of the continuous Gaussian Bell with the bar diagram of the Binomial Distribution, we have to take into account the normalization mode. For the Bell, the normalization arises from areas, while for the Binomial, it arises from the lengths of the bars. Such a problem is largely forgettable for the non-normalized Binomial Distribution because each bar can be interpreted as having an area that is given by its height multiplied by the width, that is, 1. If we accept this, we conclude that the Gauss-Bell and the Binomial Distribution coincide.\n", "\n", "Our result reads: the Gauss Bell is the envelope of the Binomial Distribution, where both distributions have the same mean, $np$, and standard deviations, $\\sqrt{np(1-p)}$. \n", "\n", "This coincidence allows us to propose an explanation of the following fact: the Gauss Bell, or a good approximation to it, is observed in the most varied situations of nature, technology, and society. Why? It is because the Normal Distribution or Gauss Bell arises from the joint operation mode of many independent factors, each one of little importance, but together play a decisive role. And the contribution of each factor is chaotic. Sometimes it can be ignored, but others not. This behavior can be modelled as the random output of a flipping coin, heads or tails.\n", "\n", "# DISCRIMINATION BETWEEN TWO THEORIES\n", "\n", "If we gather enough random data, we can decide between rival theories. Let us illustrate this with a concrete case-study.\n", "\n", "We have shown that 14 heads out of 20 can be explained by two theories: $p= 0.5$ and $p = 2/3$. This means that one single data is insufficient to decide between the two theories. We infer that few data are also insufficient but if the size of data grows, a tendency will possibly appear that we must exploit to discriminate the two theories. Let us devise a simulation to show that gathering enough data enables one to reject one of this null hypotheses. How will we do that?\n", "\n", "Let us analyze the following setting: We will have indeed two simulations. In the first, we throw coins whose heads and tails have equal probabilities to estimate the distribution of heads. This has been done already. In the second, we will throw coins whose heads have a probability of $2/3$ to generate our \"observed data\" from which we will take sets of $n$ data and take the mean. Our question is: Can we use the found distribution in step one to classify the found mean in the second step? The answer is NO! The reason is that the distribution and the observed data refer to different objects. The distribution refers to isolated data, while the observed data refers to means of n data. What shall we correct?\n", "\n", "Our new setting is the following: We fix $n$ as the number of observed data, from which we will calculate the mean that must be classified as ordinary or else strange. So, we will have two simulations. In the first, we throw coins whose heads and tails have equal probabilities, $p = 0.5$. We will throw $n$ coins, take the mean, and estimate the distribution of these means. In the second, we will exercise the rival theory. We will throw coins whose heads have a probability $p =2/3$ to generate our \"observed data,\" from which we will take a random sample of $n$ data and calculate the mean. Our question is: Can we use the found distribution in step one to classify the found mean in the second step as ordinary else strange? The answer is YES! That is because the distribution and the observed mean both have the same type or nature, both refer to means of samples of the same size. In the literature, the found distribution in step one is called the *sample mean with $n$ data*. We must inquire over the value of $n$. All we know is that if $n$ is low, we will not be able of rejecting the rival theory.\n", "\n", "Time to pass to simulations. \n", "\n", "> **Example.** Let us find the distribution of the sample mean for 20 coins and p = 0.5. The code follows. We reuse program `coinsProbOfHeadsTwoThirds.py`." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "id": "SXzTdEFGyM37" }, "outputs": [], "source": [ "# Program sampleMeanDistribution.py\n", "# We estimate the distribution of\n", "# sample means of samples of size n\n", "# generated according to the first theory, p = 0.5.\n", "\n", "\n", "# Libraries and local rename\n", "import matplotlib.pyplot as plt\n", "import numpy as np #Support for scientific computing\n", "import random\n", "from collections import Counter\n", "import statistics\n", "\n", "\n", "# 1. Parameters for the simulation\n", "num_coins_per_round = 20\n", "n = 200 # Size of each sample\n", "num_repetitions = 3000\n", "\n", "\n", "# Simulate coin flips\n", "sample_means = []\n", "for i in range(num_repetitions):\n", " heads_in_round = []\n", "\n", " for j in range(n):\n", " heads = 0\n", " for _ in range(num_coins_per_round):\n", " # 0 for tails, 1 for heads\n", " if random.choice([0, 1]) == 1:\n", " heads += 1\n", " heads_in_round.append(heads)\n", " sample_means.append(statistics.mean(heads_in_round))\n", "\n", "# Calculate distribution\n", "sample_means_distribution = Counter(sample_means)\n", "\n", "# Arrays to keep values to be graphed\n", "x_values = []\n", "y_values = []\n", "\n", "# Print the results\n", "print(\"Heads is 1 in a random choice among 0, 1\")\n", "print(\"So, the probability of heads is 0.5\")\n", "print(\"and that of tails is 0.5\")\n", "print(f\"with {num_coins_per_round} coins per round:\")\n", "print(f\"sample_mean of samples of size {n} rounds\")\n", "print(f\"Simulation of {num_repetitions} repetitions \")\n", "\n", "print(\"\\nDistribution of sample_means:\")\n", "# Sort the distribution for better readability\n", "sorted_distribution = dict(sorted(sample_means_distribution.items()))\n", "for sample_mean, count in sorted_distribution.items():\n", " x_values.append(sample_mean)\n", " y_values.append(count)\n", " print(f\" {sample_mean} heads: {count} times\")\n", "\n", "# Plot the function\n", "plt.plot(x_values, y_values, label=\"Sample_means_distribution\")\n", "\n", "# Add labels, title, and legend\n", "plt.xlabel(\"Sample_mean\")\n", "plt.ylabel(\"Count\")\n", "plt.title(\"Graph of sample means distribution\")\n", "plt.grid(True)\n", "plt.legend()\n", "\n", "# Display the plot\n", "plt.show()\n" ] }, { "cell_type": "markdown", "metadata": { "id": "eRZGsQkMyM37" }, "source": [ "> **Exercise**. Run the previous code and establish the mean and range of the distribution.\n", ">\n", ">> **Answer**: With the probability of heads equal to 0.5, samples of size 200, and 10000 repetitions, the mean of the distribution of sample-means is 10, as expected, and the range of the distribution goes from 9.46 to 10.5. \n", "\n", "> **Exercise**. Estimate the mean and range of the distribution of the sample mean with probability of heads equal to 2/3.\n", ">\n", ">> **Answer**. The previous program was rerun with a minor modification in which the probability of heads was set to 2/3. The modification was included in the line of code `if random.choice([0, 1, 1]) == 1:`. The mean of the distribution of the mean number of heads of samples of size 200 is circa 13.3, as expected. The range of the distribution goes from 12.71 until 13.91.\n", "\n", "\n", "Our example and exercise allow us to claim that when the size of the sample is 200, the observed ranges of the distributions corresponding to the two different theories, p = 0.5 and p = 2/3 = 0.66, do not overlap. As a consequence, the predictions of each theory become strange for the other. This means that a sample size of 200 suffices to discriminate between the two theories.\n", "\n", "> **Exercise**. Run the previous program for increasing values of sample size and observed what happens with the range of the observed distribution.\n", ">>\n", ">> **Answer** For a sample size 50, the observed range goes from 8.8 to 11.12. For a sample size 100, the observed range goes from 9.22 to 10.81\n", "For sample size 200, the observed range goes from 9.46 to 10.5.\n", "For a sample size 500, the observed range goes from 9.64 to 10.37.\n", ">>\n", ">> We observe that when the sample size grows, the observed distribution tends to shrink and gathers itself around the expected mean. \n", "\n", "\n", "> **Challenge**. Resolve, if you can, the following mystery created by two conflicting observations.\n", ">> - First: Our simulations clearly say that when the sample size grows, the observed distribution tends to shrink and gathers itself around the expected mean. From here, we formulate a prediction: for infinitely many data, the observed distribution of sample means has range 0 and every sample mean coincides with the expected mean. Or, the sample mean converges to the expected mean when the sample size grows to infinite.\n", ">> - Second: The number of heads in a Binomial Distribution with 2 coins takes on the values 0, 1 and 2. And for 7 coins, that number runs from 0 to 7. In general, for $n$ coins the range of the distribution goes from 0 to $n$. In short: the range grows as $n$. Therefore, the sample mean has more space to dance the higher is $n$. Hence, the bells of sample means cannot get concentrated the higher is $n$. Convergence, if any, is dubious.\n", "\n", "\n", "Just to reinforce our knowledge, we must dress the previous discussion in the official formalism. But before, let us do something to improve the understandability of our program.\n", "\n", "> **Exercise**. Encapsulate the previous program so that it could be easily invoked and understood.\n", ">>\n", ">> **Answer**. The following program encapsulates two basic procedures of the previous program. To test it, we draw anew the two sample mean distributions we are working with." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "id": "PmxqYfpLyM37" }, "outputs": [], "source": [ "# Program twoSampleMeanDistributions.py\n", "# 1.We estimate the distribution of\n", "# sample means of n data\n", "# generated according to the first theory, p = 0.5.\n", "# 2. We estimate the distribution of\n", "# sample means of n data\n", "# generated according to the second theory, p = 2/3.\n", "\n", "\n", "# Libraries and local rename\n", "import matplotlib.pyplot as plt\n", "import numpy as np #Support for scientific computing\n", "import random\n", "from collections import Counter\n", "import statistics\n", "\n", "\n", "\n", "\n", "# 1. Parameters for the simulation\n", "num_coins_per_round = 20\n", "n = 200 # Size of each sample\n", "num_repetitions = 3000\n", "\n", "\n", "def sample_mean_distribution(sample_size, p):\n", " # Simulate coin flips\n", " sample_means = []\n", " for i in range(num_repetitions):\n", " heads_in_round = []\n", " for j in range(sample_size):\n", " heads = 0\n", " for _ in range(num_coins_per_round):\n", " # 0 for tails, 1 for heads\n", " if p == 0.5:\n", " # Probability of heads = 0.5\n", " if random.choice([0, 1]) == 1:\n", " heads += 1\n", " if p == 0.66:\n", " # Probability of heads = 2/3 = 0.66\n", " if random.choice([0, 1, 1]) == 1:\n", " heads += 1\n", " heads_in_round.append(heads)\n", " sample_means.append(statistics.mean(heads_in_round))\n", "\n", " # Calculate distribution\n", " sample_means_distribution = Counter(sample_means)\n", " return sample_means_distribution\n", "\n", "\n", "\n", "def graph_sample_mean_distribution(sample_means_distribution, n, p):\n", " # Arrays to keep values to be graphed\n", " x_values = []\n", " y_values = []\n", "\n", " # Print the results\n", " print(\"Heads is 1 in a random choice among 0, 1\")\n", " print(\"So, the probability of heads is 0.5\")\n", " print(\"and that of tails is 0.5\")\n", " print(f\"with {num_coins_per_round} coins per round:\")\n", " print(f\"sample_mean of samples of size {n} rounds\")\n", " print(f\"Simulation of {num_repetitions} repetitions \")\n", "\n", " print(\"\\nDistribution of sample_means:\")\n", " # Sort the distribution for better readability\n", " sorted_distribution = dict(sorted(sample_means_distribution.items()))\n", " for sample_mean, count in sorted_distribution.items():\n", " x_values.append(sample_mean)\n", " y_values.append(count)\n", " print(f\" {sample_mean} heads: {count} times\")\n", "\n", " # Plot the function\n", " plt.plot(x_values, y_values, label=\"Sample_means_distribution\")\n", "\n", " # Add labels, title, and legend\n", " plt.xlabel(\"Sample_mean\")\n", " plt.ylabel(\"Count\")\n", " plt.title(f\"Graph of sample means distribution, sample_size = {n}, p = {p}\")\n", " plt.grid(True)\n", " plt.legend()\n", "\n", " # Display the plot\n", " plt.show()\n", "\n", "def main():\n", " # Probability of heads,\n", " # only two possible values: 0.5 and 0.66\n", " p = 0.5\n", " # n is the size of the sample.\n", " sample_means_distribution = sample_mean_distribution(n, p)\n", " graph_sample_mean_distribution(sample_means_distribution, n, p)\n", " p = 0.66\n", " sample_means_distribution = sample_mean_distribution(n, p)\n", " graph_sample_mean_distribution(sample_means_distribution, n, p)\n", "\n", "\n", "main()" ] }, { "cell_type": "markdown", "metadata": { "id": "8h-BnzOwyM38" }, "source": [ "\n", "\n", "> **Exercise**. Discriminate between the two theories, p = 0.5 and p = 2/3 = 0.66, by showing that, if the size of the sample is large enough, a datum, obtained according to the predictions of one theory, is strange according to the distribution dictated by the other.\n", ">>\n", ">> **Answer**. We estimate the distribution of the sample means of size 200 according to the theory that says that p = 0.5. With that distribution, we classify events in ordinary and strange. Then, we find a mean of a sample of size 200 generated by the second theory, p = 2/3 = 0.66, and use the previous classification to decide if its ordinary, making the two theories compatible, or strange, making the two theories incompatible. The code follows. \n", "\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "id": "iCAsydDTyM38" }, "outputs": [], "source": [ "# Program RejectingRivalTheory.py\n", "# We illustrate that two theories can be discriminated\n", "# if we gather enough data.\n", "# First theory: p = 0.5.\n", "# Second theory. p = 2/3.\n", "# Procedure:\n", "# 1. We fix n, the size of observed data.\n", "# 2. We estimate the distribution of\n", "# sample means of n data\n", "# generated according to the first theory, p = 0.5.\n", "# 3. We generate\n", "# $n$ random data according to the second theory, p = 2/3,\n", "# and find the \"observed mean\".\n", "# We classify that \"observed mean\" according to the\n", "# distribution found in the second step.\n", "# We take the sample size n = 200 because\n", "# we already know that this suffices to\n", "# dicrimintae the two theories, i.e.\n", "# the \"observed mean\" becomes strange\n", "# according to the first theory.\n", "\n", "\n", "# Libraries and local rename\n", "import matplotlib.pyplot as plt\n", "import numpy as np #Support for scientific computing\n", "import random\n", "from collections import Counter\n", "import statistics\n", "\n", "\n", "# 1. Parameters for the simulation\n", "num_coins_per_round = 20\n", "n = 200 # Size of each sample\n", "num_repetitions = 1000\n", "\n", "def mean_of_sample(sample_size, p):\n", " heads_in_round = []\n", " for j in range(sample_size):\n", " heads = 0\n", " for _ in range(num_coins_per_round):\n", " # 0 for tails, 1 for heads\n", " if p == 0.5:\n", " # Probability of heads = 0.5\n", " if random.choice([0, 1]) == 1:\n", " heads += 1\n", " if p == 0.66:\n", " # Probability of heads = 2/3 = 0.66\n", " if random.choice([0, 1, 1]) == 1:\n", " heads += 1\n", " heads_in_round.append(heads)\n", " sample_mean = statistics.mean(heads_in_round)\n", " return sample_mean\n", "\n", "\n", "def sample_mean_distribution(sample_size, p):\n", " # Simulate coin flips\n", " sample_means = []\n", " for i in range(num_repetitions):\n", " sample_means.append(mean_of_sample(sample_size, p))\n", " # Calculate distribution\n", " sample_means_distribution = Counter(sample_means)\n", " return sample_means_distribution\n", "\n", "\n", "\n", "def graph_sample_mean_distribution(sample_means_distribution, n, p):\n", " # Arrays to keep values to be graphed\n", " x_values = []\n", " y_values = []\n", "\n", " # Print the results\n", " print(\"Heads is 1 in a random choice among 0, 1\")\n", " print(\"So, the probability of heads is 0.5\")\n", " print(\"and that of tails is 0.5\")\n", " print(f\"with {num_coins_per_round} coins per round:\")\n", " print(f\"sample_mean of samples of size {n} rounds\")\n", " print(f\"Simulation of {num_repetitions} repetitions \")\n", "\n", " print(\"\\nDistribution of sample_means:\")\n", " # Sort the distribution for better readability\n", " sorted_distribution = dict(sorted(sample_means_distribution.items()))\n", " for sample_mean, count in sorted_distribution.items():\n", " x_values.append(sample_mean)\n", " y_values.append(count)\n", " print(f\" {sample_mean} heads: {count} times\")\n", "\n", " range_in = x_values[0]\n", " range_out = x_values[len(x_values)-1]\n", " range = range_out - range_in\n", " print(f\"range_in: {range_in}, range_out: {range_out}, range = {range}\")\n", "\n", " # Plot the function\n", " plt.plot(x_values, y_values, label=\"Sample_means_distribution\")\n", "\n", " # Add labels, title, and legend\n", " plt.xlabel(\"Sample_mean\")\n", " plt.ylabel(\"Count\")\n", " plt.title(f\"Graph of sample means distribution, sample_size = {n}, p = {p}\")\n", " plt.grid(True)\n", " plt.legend()\n", "\n", " # Display the plot\n", " plt.show()\n", "\n", "def main():\n", " # Probability of heads,\n", " # only two possible values: 0.5 and 0.66\n", " p = 0.5\n", " sample_means_distribution = sample_mean_distribution(n, p)\n", " graph_sample_mean_distribution(sample_means_distribution, n, p)\n", " p = 0.66\n", " mean = mean_of_sample(n, p) # n = sample_size\n", " print(f\"Random sample of size = {n}, p = {p}, with mean = {mean}\")\n", "\n", "main()" ] }, { "cell_type": "markdown", "metadata": { "id": "hyMxkK6XyM38" }, "source": [ ">> We clearly see that the mean of a random sample of size 200 generated by the theory p = 0.66 is out of the range of the estimated distribution of sample means with p = 0.5. This means that a sample size of 200 is enough to unveil the incompatibility of the two theories $p = 0.5$ and $p = 0.66$.\n", "\n", "> **Exercise**. Add to the previous program a procedure to graph the two tails of strange events for the portrayed distribution, corresponding to the theory $ p = 0.5$, together with a vertical line pointing to the value of the sample mean generated by the second theory, $p = 0.66$.\n", "\n", ">> **Answer**. The following program draws an estimation of the distribution of sample means with sample means of size 200 generated by the theory that says that $p = 0.5$. The strange events are colored in red given a significance of $0.05$. Next, a sample of size 200 is fabricated from the second theory that says that $p = 0.66$. Its mean is found and drawn to see whether it is normal according to the first theory or abnormal. If this mean is normal, we have not gathered enough data to discriminate them. If this mean is abnormal, we have gathered enough data to establish the discrimination of the two theories.\n", "The code follows. To color the rejection regions, we reuse the code of program `binomialDistriColoredTails2.py`." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "id": "kE_--9AZyM38" }, "outputs": [], "source": [ "# Program RejectingRivalTheory2.py\n", "# We illustrate that two theories can be discriminated\n", "# if we gather enough data.\n", "# First theory: p = 0.5.\n", "# Second theory. p = 2/3.\n", "# Procedure:\n", "# 1. We fix n, the size of observed data.\n", "# 2. We estimate the distribution of\n", "# sample means of n data\n", "# generated according to the first theory, p = 0.5.\n", "# 3. We generate\n", "# $n$ random data according to the second theory, p = 2/3,\n", "# and find the \"observed mean\".\n", "# We classify that \"observed mean\" according to the\n", "# distribution found in the second step.\n", "# We take the sample size n = 200 because\n", "# we already know that this suffices to\n", "# dicriminate the two theories, i.e.\n", "# the \"observed mean\" becomes strange\n", "# according to the first theory.\n", "\n", "\n", "# Libraries and local rename\n", "from matplotlib.collections import LineCollection\n", "import matplotlib.pyplot as plt\n", "import numpy as np #Support for scientific computing\n", "import random\n", "from collections import Counter\n", "import statistics\n", "\n", "\n", "# Parameters for the simulation\n", "num_coins_per_round = 20\n", "n = 200 # Size of each sample\n", "num_repetitions = 5000\n", "# Initialization of the sample mean\n", "# according to the second theory p = 0.66\n", "mean_of_sample_2 = -1\n", "\n", "#Fabricate a sample and take the mean\n", "def mean_of_sample(sample_size, p):\n", " heads_in_round = []\n", " for j in range(sample_size):\n", " heads = 0\n", " for _ in range(num_coins_per_round):\n", " # 0 for tails, 1 for heads\n", " if p == 0.5:\n", " # Probability of heads = 0.5\n", " if random.choice([0, 1]) == 1:\n", " heads += 1\n", " if p == 0.66:\n", " # Probability of heads = 2/3 = 0.66\n", " if random.choice([0, 1, 1]) == 1:\n", " heads += 1\n", " heads_in_round.append(heads)\n", " sample_mean = statistics.mean(heads_in_round)\n", " # rounding to one decimal place\n", " sample_mean = int(sample_mean * 10) / 10\n", " return sample_mean\n", "\n", "\n", "# Simulate data to estimate the distribution of the sample mean.\n", "def sample_mean_distribution(sample_size, p):\n", " # Simulate coin flips\n", " sample_means = []\n", " for i in range(num_repetitions):\n", " sample_means.append(mean_of_sample(sample_size, p))\n", " # Calculate distribution\n", " sample_means_distribution = Counter(sample_means)\n", " return sample_means_distribution\n", "\n", "\n", "\n", "def graph_sample_mean_distribution(sample_means_distribution, n, p):\n", " # Arrays to keep values to be graphed\n", " x_values = []\n", " y_values = []\n", "\n", " # Print the results\n", " print(\"Heads is 1 in a random choice among 0, 1\")\n", " print(\"So, the probability of heads is 0.5\")\n", " print(\"and that of tails is 0.5\")\n", " print(f\"with {num_coins_per_round} coins per round:\")\n", " print(f\"sample_mean of samples of size {n} rounds\")\n", " print(f\"Simulation of {num_repetitions} repetitions \")\n", "\n", " print(\"\\nDistribution of sample_means:\")\n", " # Sort the distribution for better readability\n", " sorted_distribution = dict(sorted(sample_means_distribution.items()))\n", " cumulate_frequency = 0\n", " cumulate_frequencies = []\n", " for sample_mean, count in sorted_distribution.items():\n", " relative_frequency = 1.0 * count / num_repetitions\n", " x_values.append(sample_mean)\n", " y_values.append(relative_frequency)\n", " cumulate_frequency += relative_frequency\n", " cumulate_frequencies.append(cumulate_frequency)\n", " print(f\" mean of heads: {sample_mean}, relative frequency: {relative_frequency} , cumulated : {cumulate_frequency}\")\n", "\n", "\n", " ### The graphic with normal and abnormal events distinctly colored.\n", "\n", " # fig, ax = plt.subplots()\n", "\n", "\n", " bar_labels = []\n", " bar_colors = []\n", "\n", " left_tail = 0.0\n", " real_alpha_half = -1.0\n", "\n", " epsilon = 1e-7\n", " for rel_freq in y_values:\n", " print(f\"rel_freq = y_value = {rel_freq}\")\n", " left_tail += rel_freq\n", " if left_tail < 0.025:\n", " bar_labels.append('red')\n", " bar_colors.append('tab:red')\n", " real_alpha_half = left_tail\n", " red_area = 2 * real_alpha_half\n", " print(f\" left_tail = {left_tail} red\")\n", " if ((left_tail >= 0.025) and (1-left_tail + epsilon > real_alpha_half)) :\n", " bar_labels.append('blue')\n", " bar_colors.append('tab:blue')\n", " print(f\" {1-left_tail + epsilon} > {real_alpha_half} blue\")\n", " if (1-left_tail + epsilon) <= real_alpha_half:\n", " bar_labels.append('red')\n", " bar_colors.append('tab:red')\n", " print(f\" {1-left_tail + epsilon} <= {real_alpha_half} red\")\n", "\n", " # Automatic bar ploting: malfunctions\n", " # ax.bar(x_values, y_values, label=bar_labels, color=bar_colors)\n", " # ax.set_ylabel('Probability')\n", " # ax.set_xlabel('Sample_mean')\n", " # ax.set_title(f\"Sample Mean Distribution with rejection regions in red, area = {red_area:.4f}\")\n", " # ax.legend(title='Bar color')\n", " # plt.show()\n", "\n", " print(\"real_alpha_half = \", real_alpha_half)\n", " print( \"mean relat_Freq cumulate bar_labels bar_colors \")\n", " for i in range(len(y_values)):\n", " print(f\" {x_values[i]} {y_values[i]} {cumulate_frequencies[i]} { bar_labels[i]} {bar_colors[i]}\")\n", "\n", " # We draw the bars as vertical lines\n", " # Format: line = [(x1, y1), (x2, y2) ]\n", " # Collection of bar lines, one for each x value\n", " bar_collection = []\n", " for i in range(len(x_values)):\n", " my_bar = [(x_values[i], 0),(x_values[i], y_values[i])]\n", " bar_collection.append(my_bar)\n", " print(f\"{i} {my_bar}\")\n", " # Make a random sample according to p = 0.66 and take the mean\n", " p = 0.66\n", " mean_of_sample_2 = mean_of_sample(n, p)\n", " print(f\"Random binomial sample of size = {n}, p = {p}, with mean = {mean_of_sample_2}\")\n", "\n", " #Red long line marking a sample mean according to p = 0.66\n", " bar_of_sample_mean_2 = [(mean_of_sample_2, 0),(mean_of_sample_2, 0.4)]\n", " bar_collection.append(bar_of_sample_mean_2)\n", " collection = LineCollection(bar_collection, color=bar_colors)\n", " fig, axes = plt.subplots(1, 1)\n", " axes.add_collection(collection)\n", " axes.autoscale()\n", " axes.set_ylabel('Probability')\n", " axes.set_xlabel('Sample_mean')\n", " axes.set_title(f\"Sample Mean Distribution for p = 0.5, sample size = 20, confidence = {red_area:.4f} \\n and Observed sample mean for p = 0.66, red vertical line at the right.\")\n", " plt.show()\n", "\n", "def main():\n", " # The probability of heads has\n", " # only two possible values: 0.5 and 0.66\n", " p = 0.5 # The distribution of sample means is calculated for this theory\n", " # n = sample_size\n", " sample_means_distribution = sample_mean_distribution(n, p)\n", " graph_sample_mean_distribution(sample_means_distribution, n, p)\n", "\n", "main()" ] }, { "cell_type": "markdown", "metadata": { "id": "wTa7nekSyM39" }, "source": [ ">> The graphic shows clearly that\n", ">>> - The sample means according to the first theory, p = 0.5, fall very close to the expected mean, in this case, 10.\n", ">>> - The unique sample mean according to the second theory, p = 2/3 = 0.66, falls upwards and very far away from the observed distribution of sample means according to the first theory. Thus, that sample mean is very strange, and henceforth, the second theory must be rejected by the first one.\n", "\n", "\n", "\n", "# THE ORIGIN OF THE SHRINKING TENDENCY\n", "\n", "Theories must mature and explain what may initially appear to be a mystery. Otherwise, they are forgotten because they are as blind alleys, from which no progress can be made.\n", "\n", "So, let us explain various observed patterns that are related with a very practical rule: ***To gain precision, augment the size of the sample.*** This result ultimately rests upon the fact that when the sample size grows, the bells of sample means shrink towards their mean. And everything here is built upon the tendency of the distributions of certain observables to reach a definite shape. It is this definiteness that allows us to make predictions that must be fulfilled. The material is almost self-contained and an effort has been made to make it understandable.\n", "\n", "\n", "Concretely, we will pay attention to two facts that we shall explain right down here: the tendency of the relative frequency of a Binomial Distribution towards the probability, and that of the sample mean of the Binomial Distribution to shrink towards the expected mean, a concept that deserves more accurateness.\n", "\n", "# The Expected Value \n", "\n", "The mean of a data set is found by adding all the data and averaging over the number of data. For instance, the mean $\\overline{X}$ of the data set 3, 3, 4, 5, 5, 5 is\n", "\n", "> $\\overline{X} = \\frac{3+3+4+5+5 +5}{6} = \\frac{25}{6}$\n", "\n", "Since some terms are repeated, we can rewrite this fraction using absolute frequencies: that of 3 is 2, that of 4 is 1, and that of 5 is 3. The mean of this data set becomes:\n", "\n", "> $\\overline{X} = \\frac{3 \\times 2 +4 \\times 1 +5\\times 3}{2 + 1 + 3} = \\frac{3 \\times 2 +4 \\times 1 +5\\times 3}{6} = \\frac{25}{6}$\n", "\n", "If we distribute the fraction, we get the mean in terms of relative frequencies:\n", "\n", "> $\\overline{X} = 3 \\times \\frac{2}{6} +4 \\times \\frac{1}{6} +5\\times \\frac{3}{6} = \\frac{25}{6}$.\n", "\n", "This allows us to define the *general formula for the mean* of an enumerable set of data $x_i$ in terms of their relative frequencies $f_i$:\n", "\n", "> $\\overline{X} = \\sum x_i f_i $\n", "\n", "\n", "where $x_i$ is the datum $i$ and $f_i$ is its relative frequency.\n", "\n", "Relative frequencies are observable. They must be compared with what we think. So we ask, what is that ideal item that the relative frequency reflects? It is the probability, which on many occasions is simply an idealized relative frequency. Other times, they correspond with fine engineering work, say, the construction of a perfect die.\n", "\n", "> **Challenge.** Discuss the following rule that is correct: to decree that a relative frequency is a probability is allowed when the sample size is large enough, say, 10 millions.\n", "\n", "If we want the concept of probability to be of some use, we need to make sound definitions, as the following, which is the idealization of the mean of observed data:\n", "\n", "\n", "Let $X$ be a random variable with discrete, enumerable domain. The ***Expected Value*** $E(X)$ of $X$ is defined as:\n", ">\n", "> $E(X) = \\sum x_i p_i $\n", ">\n", "> where $x_i$ is the value $i$ and $p_i$ is its probability. With this, we are requiring that observed relative frequencies converge to the probability. Let notice that this is not a definition directly tied to observation, but to ideas, because probability is, in first instance, a mental concept. The only requisite is that values $x_i$ must be observable. Thus, a discrete random variable must be defined as a function that associates probabilities to events that are results of measurements represented by non-negative integers, 0, 1, 2,....\n", "\n", "> **Example.** For a Binomial Distribution or a Binary Random Variable that counts the number of heads when 20 coins are thrown, $x_i = i$ represents the number of heads and runs from $0$ to $20$. At the same time, $p_i$ is given by the Binomial Expansion that represents what we think.\n", "\n", "\n", "\n", "> **Example.** Let us calculate the Expected Value of the Binomial distribution that counts heads with 1 coin and probability of heads $p$. The corresponding random variable $X_1$ is the function defined by\n", ">\n", "> $0 \\rightarrow q = 1-p$\n", ">\n", "> $1 \\rightarrow p$\n", ">\n", "> The Expected Value of $X_1$ is\n", ">\n", "> $E(X_1) = 0\\times q + 1 \\times p = p$\n", "\n", "\n", "> **Example.** Let us calculate the Expected Value of the Binomial distribution with 2 coins and probability of heads $p$. The corresponding random variable $X_2$ is the function defined by\n", ">\n", "> $0 \\rightarrow q^2$\n", ">\n", "> $1 \\rightarrow 2pq$\n", ">\n", "> $2 \\rightarrow p^2$\n", ">\n", "> and 0 otherwise.\n", ">\n", "> The Expected Value of $X_2$ is the sum of the products of values times their probabilities:\n", ">>\n", ">> $E(X_2) = 0\\times q^2 + 1 \\times 2pq + 2 \\times p^2 = 2pq + 2p^2 = 2p(q +p) = 2p$\n", "\n", "\n", "> **Exercise.** The Expected Value of the Binomial Distribution with 1 coin was $p$ and with 2 coins was $2p$. Is it $3p$ for 3 coins? Explain your result by common sense.\n", "\n", ">> **Answer.** To know if $3p$ is the Expected Value for a Binomial Distribution with 3 coins, we must calculate it. Let $X_3$ be the random variable representing the Binomial Distribution with 3 coins. It is the function defined by\n", ">>>\n", ">>> $0 \\rightarrow q^3$\n", ">>>\n", ">>> $1 \\rightarrow 3q^2p$\n", ">>>\n", ">>> $2 \\rightarrow 3qp^2$\n", ">>>\n", ">>> $3 \\rightarrow p^3$\n", ">>>\n", ">>> and its Expected Value is:\n", ">>>\n", ">>> $E(X_3) = 0 \\times q^3 + 1 \\times 3q^2p + 2 \\times 3qp^2 +3 \\times p^3$\n", ">>>\n", ">>> $ = 3q^2p + 6qp^2 + 3p^3 = 3p (q^2 + 2 pq + p^2) = 3p(p+q)^2 = 3p$.\n", "\n", "Since we have verified that the Expected Value for 3 coins is $3p$, let us postulate that the Expected Value of the Binomial Distribution with parameter $p$ and $n$ coins is $np$. This result says that if we throw $n$ coins, we can expect $np$ heads. For instance, if $p= 1/2$ and we flip the coin 10 times, we expect $10 \\times (1/2) = 5$ heads. In a real experiment, we will rarely get 5 heads, but values will stack around 5 with a distribution given by the Binomial Expansion.\n", "\n", "\n", "To proceed further with the study of our postulate, we need two properties of the Expected Value of random variables.\n", "\n", "- The first is: If all values are multiplied by a constant $c$, the Expected Value is also multiplied by the same constant:\n", ">>\n", ">> $E(cX) = c E(X) $\n", "\n", "> **Example.** Let us think of a tool seller whose monthly income depends on the volume of sales. This is random, but is homogeneous over time. So, it makes sense to speak of probabilities of gaining 1000, 2000, ...... It is also understandable that from time to time the tools, which are indispensable, become more expensive. In consequence, the revenue will increase. Then, this rule says that the expected value of his incomes is multiplied by a constant $c$ if the income of each month is multiplied by the same constant.\n", "\n", "- The second rule refers to a sum of random variables. Let us think of a couple, husband and wife, who work as seller whose income depends on the sales. These are random, but are homogeneous over time. So, it makes sense to speak of probabilities of gaining 1000, 2000,..., 6000. Then, this rule says that the Expected Value of the income of the pair is the sum of the Expected Values of each member of the couple. More specifically: The Expected value of the sum of two random variables is the sum of the expected Values of each variable. Let $X$ and $Y$ be the incomes, then \n", ">>\n", ">> $E(X+Y) = E(X) + E(Y)$\n", "\n", "Let us see a direct application of these rules.\n", "\n", "> **Example.** We have postulated that the Expected Value of the $Bi(n,p)$ is $np$. We might use algebra to prove the postulate for some values, 4, 5, and 7. And, the rest? How can we pass from a postulate to a theorem?\n", ">>\n", ">> When a proposition about numbers 0, 1, 2, ... must be proved, we use the Induction Principle: We prove the proposition for a single case, it might be 1, and we then show that we can always prove that the proposition is valid for $n+1$ assuming that it is valid for $n$. So, let us use this principle to prove that the Expected Value of Binomial distribution is $np$ for every $n$.\n", ">>\n", ">> Step 1: Is it valid for 1? Yes, that was already shown.\n", ">>\n", ">> Step 2: Let us assume that the proposition is valid for $n$, i.e., the Expected value of $Bi(n,p)$ is $np$. From this, we must prove that the expected value of $Bi(n+1,p)$ is $(n+1)p$.\n", ">>\n", ">> Now, $Bi(n+1,p) = Bi(n,p) + Bi(1,p)$, which means that a game with $n+1$ coins is just a game with $n$ coins plus another coin. More clearly, we throw $n$ coins, count heads, and then, we throw another coin, look whether it falls heads, and in that case we add 1 to the number of heads previously found.\n", ">>\n", ">> Since the Expected Value is compatible with the sum of random variables and with the multiplication by a number, we have:\n", ">>>\n", ">>> $E(Bi(n+1,p)) = E(Bi(n,p)) + E(Bi(1,p)= np + p = (n+1)p$\n", ">>>\n", ">> such as desired. We have proved an important result that reads:\n", "\n", "> **Theorem.** The Expected Value of the Binomial distribution $Bi(n,p)$ is $np$.\n", "\n", "\n", "It might be instructive to try to prove by induction a false statement, say\n", ">\n", ">> $E(Bi(n,p)) = np + (n-1)p$:\n", ">\n", "> Step 1: Is it valid for $n= 1$?\n", ">>\n", ">> $E(Bi(1,p)) = 1p + (1-1)p = p$, an affirmation that is correct.\n", ">>\n", "> Step 2: let us assume that the proposition is valid for $n$, i.e., the Expected Value of $Bi(n,p)$ is $np + (n-1)p$. From this, we must prove that the Expected Value of $Bi(n+1,p)$ is\n", ">>\n", ">> $E(Bi(n+1,p)) = (n+1)p + (n+1-1)p = (n+1)p + np = np + p + np = 2np + p$. Is that true? Let us see:\n", ">>\n", ">> $E(Bi(n+1,p)) = E(Bi(n,p)) + E(Bi(1,p))= np + (n-1)p + p = np + np -p + p = 2np $ which is different from $2np + p$. Therefore, the proposition is false.\n", "\n", "\n", "> **Example.** Let $X$ and $Y$ be both random variables that take values on the set 1,2,..., 6, (in thousands) with probabilities equal to $1/6$. These variables represent the random inputs of a couple. So, the Expected Value of each variable is\n", ">>\n", ">> $E(X) = E(Y) = 1 \\times \\frac{1}{6} + 2 \\times \\frac{1}{6} + 3 \\times \\frac{1}{6} + 4 \\times \\frac{1}{6} + 5 \\times \\frac{1}{6} + 6 \\times \\frac{1}{6} = \\frac{21}{6} = 3.5$\n", ">>\n", ">> Therefore, the Expected Value of the familiar income is\n", ">>\n", ">> $E(X+Y) = E(X) + E(Y) = \\frac{21}{6} + \\frac{21}{6} = \\frac{42}{6} = 7$\n", "\n", "It is curious that we have found the Expected Value of the income of the couple $Z = X + Y$, although we have no idea what $Z$ is. So, let us try to understand its nature. The problem is that if the husband gains 2, we cannot say what could be the gain of the wife. All we can assume is that they work in different districts, so they do not compete, thus their gains are random but uncorrelated.\n", "\n", "The following table helps us to decipher what $Z = X + Y$ is:\n", "\n", "Z | 1 & 2 & 3 & 4 & 5 & 6 \n", "1 | 2 & 3 & 4 & 5 & 6 & 7 \n", "2 | 3 & 4 & 5 & 6 & 7 & 8 \n", "3 | 4 & 5 & 6 & 7 & 8 & 9 \n", "4 | 5 & 6 & 7 & 8 & 9 & 10 \n", "5 | 6 & 7 & 8 & 9 & 10 & 11 \n", "6 | 7 & 8 & 9 & 10 & 11 & 12 \n", "\n", "We see in the first column the possible values of $X$, in the first row those of $Y$ and in the other 6 rows, the sum of $X$ and $Y$. The Latex code for this table is the following:\n", "\n", "\n", "\\begin{center}\n", "\\begin{tabular}[c]{|r||rrrrrr|} \\hline \\hline\n", "\\multicolumn{7}{|c|}{ Table of Z = X + Y} \\\\ \\hline \\hline\n", "X ($\\downarrow)$ and Y ($\\rightarrow$)& 1 & 2 & 3 & 4 & 5 & 6\\\\ \\hline \\hline\n", "1 & 2 & 3 & 4 & 5 & 6 & 7\\\\\n", "2 & 3 & 4 & 5 & 6 & 7 & 8\\\\\n", "3 & 4 & 5 & 6 & 7 & 8 & 9\\\\\n", "4 & 5 & 6 & 7 & 8 & 9 & 10\\\\\n", "5 & 6 & 7 & 8 & 9 & 10 & 11\\\\\n", "6 & 7 & 8 & 9 & 10 & 11 & 12\\\\ \\hline\n", "\\end{tabular}\n", "\\end{center}\n", "\n", "So, $Z = X + Y$ is the random variable that takes values on 2, 3,..., 12 with the following probabilities:\n", ">\n", "> $2 \\rightarrow 1/36$\n", ">\n", "> $3 \\rightarrow 2/36$\n", ">\n", "> $4 \\rightarrow 3/36$\n", ">\n", "> $5 \\rightarrow 4/36$\n", ">\n", "> $6 \\rightarrow 5/36$\n", ">\n", "> $7 \\rightarrow 6/36$\n", ">\n", "> $8 \\rightarrow 5/36$\n", ">\n", "> $9 \\rightarrow 4/36$\n", ">\n", "> $10 \\rightarrow 3/36$\n", ">\n", "> $11 \\rightarrow 2/36$\n", ">\n", "> $12 \\rightarrow 1/36$\n", ">\n", "> Test: The sum of all probabilities must add to 1. Is that true?\n", ">\n", "> Yes: the sum of numerators is 1 + 2 + 3 + 4 + 5+ 6 + 5 + 4 + 3 + 2 +1 = 36. Since the denominator is also 36, the sum of probabilities is 1.\n", ">\n", "> Prediction. The Expected Value of $Z$ is the sum of the expected Values of $X$ and $Y$, i.e., $3.5 + 3.5 = 7$. Is that true? Let us see:\n", ">\n", "> $E(Z) = \\frac{2}{36} + \\frac{6}{36} +\\frac{12}{36} +\\frac{20}{36} +\\frac{30}{36} +\\frac{42}{36} +\\frac{40}{36} +\\frac{36}{36} +\\frac{30}{36} +\\frac{22}{36} +\\frac{12}{36} = \\frac{252}{36} = 7$\n", ">\n", "> Thus, it is true that the Expected Value is $7$. The rule of the addition of Expected Values of random variables was verified. And, what about the general proof?\n", "\n", "> *Research.* The general proof of the proposition that the Expected Value of random variables is linear, i.e., compatible with both multiplication by a scalar or number, and with the addition of random variables, already exists. Revise the proof offered by Wikipedia and try to understand it:\n", ">\n", "> https://statproofbook.github.io/P/mean-lin\n", "\n", "> *Challenge.* In the previous example, $X$ and $Y$ abstracted a game with a die, in which all faces fall upwards with the same probability $1/6$. This is an example of a Uniform Distribution, whose graph is a horizontal line. But the shape of the distribution of $X + Y$, abstracting a game with two dice in which the sum is what counts, is not uniform. It is mountain-like. Make a prognostic for three dice, $X + X + X$ and test it. Would you guess what could happen for a sum of many terms $X + X + ... + X$? Is that a happy coincidence or is it general?\n", "\n", "We consider that we have grasped a good idea of what the Expected Value of a random variable is. Now, we can go one step further: to understand why and how the sample mean tends to gather itself around the Expected Value. To prove that, we need to know how to measure the dispersion or spreading of a random variable. Let us proceed.\n", "\n", "# The dispersion of a random variable\n", "\n", "We already know that the sum of the deviations with respect to the mean of a data set balance one with another, they add up to zero. But, the deviations are also indicative of the spreading of data. To recover this meaning, we square all deviations and take the average:\n", ">\n", "> $V (x) = \\frac{\\sum (x_i - \\overline{x})^2}{n}$.\n", "\n", "> **Example.** Let the data set be 4, 5, 6. Then the mean is\n", ">>\n", ">> $\\overline{x} = \\frac{4+5 + 6}{3} = \\frac{15}{3} = 5$\n", ">>\n", ">> Thus, we obtain the variance $V(X)$:\n", ">>\n", ">> $V (x) = \\frac{ (4- 5)^2 + (5-5 )^2 + (6-5 )^2 }{n} = \\frac{ (-1)^2 + (0 )^2 + (1 )^2 }{3} = \\frac{1+1}{3} = \\frac{2}{3}$.\n", "\n", "How must we define the variance of a random variable $X$? We take the deviation from the mean, i.e., from the Expected Value, square them, and take the average. The Expected Value is $E(X)$, the deviations from this value are $(X- E(X))$, whose square is $(X- E(X))^2$, and to finish, we average, i.e., take the Expected Value. Thus, ***the variance $V$ of the random variable $X$*** is defined as\n", "\n", "> $V(X) = E[(X- E(X))^2] = E(X- E(X))^2$.\n", "\n", "The variance is a measure of dispersion, spreading or variation of data sets and random variables. The higher are the deviations from the mean, the higher is the variance. Since, it is an expected value or average, it is very robust against changes in samples but can be distorted by outliers, which are very strange events.\n", "\n", "\n", "Two properties of $V$ will be used right now:\n", "\n", "- If all data are multiplied by a constant $c$, the variance is multiplied by the same constant squared because the deviations are squared:\n", ">>\n", ">> $V(cX) = E(cX- E(cX))^2 = E(cX- cE(X))^2 = E(c(X- E(X)))^2 = E(c^2(X- E(X))^2) = c^2E(X- E(X))^2 = c^2 V(X) $\n", "\n", "- The variance of a sum is the sum of variances plus a residue that might be zero:\n", ">>\n", ">> $V(X+ Y) = E(X + Y - E(X + Y))^2 = E(X + Y - E(X) - E( Y))^2 $\n", ">>\n", ">> $ = E[ (X - E(X))+ (Y - E( Y))]^2 $\n", ">>\n", ">> $ = E[(X - E(X))^2+ (Y - E( Y))^2 + 2(X - E(X)) (Y - E( Y))]$\n", ">>\n", ">> $ = E(X - E(X))^2+ E(Y - E( Y))^2 + 2E[(X - E(X)) (Y - E( Y))]$\n", ">>\n", "> We have gotten:\n", ">>\n", ">> $V(X+ Y) = V(x) + V(Y) + 2E[(X - E(X)) (Y - E( Y))]$\n", "\n", "The residue $2E[(X - E(X)) (Y - E( Y))]$ has a meaning: it is twice the average of the products of the deviations off their corresponding means of the random variables $X$ and $Y$. To get a non-zero value, those deviations must have the same sign in the majority of cases, and the two variables must coordinately dance around their mean values. We say that they are correlated. But if the variables are independent, the terms of the residue will cancel one another because they are uncorrelated. This is precisely the case with our coins: the result of a flip does not influence any other one. In this way, we have shown the following result that will be employed in the sequel:\n", "\n", "> **Theorem.** The variance of a sum of independent random variables is the sum of the variances.\n", "\n", "Let us now pass to specific calculations.\n", "\n", "# Variance of the Binomial Distribution\n", "\n", "Let us calculate the variance of the Binomial Distribution in the simplest cases.\n", "\n", "> **Example.** Let us calculate the variance of the Binomial Distribution for $n=1$. The corresponding random variable $X_1$ is the function defined by\n", ">>\n", ">> $0 \\rightarrow q = 1-p$\n", ">>\n", ">> $1 \\rightarrow p$\n", ">>\n", ">> The Expected Value of $X_1$ is\n", ">>\n", ">>> $E(X_1) = 0\\times q + 1 \\times p = p$\n", ">>\n", ">> The variance of $X_1$ is\n", ">>\n", ">>> $V(X_1) = (0-p)^2 q + (1-p)^2 p = p^2q + (1 -2p + p^2)p = p^2q + p -2p^2 + p^3 = p + p^2(q-2+p) = p + p^2(q +p -2) = p + p^2(1 -2) = p + p^2(-1) = p - p^2 = p (1-p) = pq$\n", "\n", "> **Example.** Let us calculate the variance of the Binomial Distribution for $n=2$. The corresponding random variable $X_2$ is the function\n", ">>\n", ">> $0 \\rightarrow q^2$\n", ">>\n", ">> $1 \\rightarrow 2pq$\n", ">>\n", ">> $2 \\rightarrow p^2$\n", ">>\n", ">> We do not want to calculate this variance from the definition. We prefer to use the properties of the variances. So, $X_2 = X_1 + X_1$, which means that the number of heads in a game with two coins is equal to adding the result of two games, both with 1 coin. We have:\n", ">>> \n", ">>> $V(X_2) = V(X_1) + V(X_1)$\n", ">>\n", ">> where we have used the fact that playing with one coin does not interfere with the results of the game with the second coin. Thus, we are dealing with uncorrelated, independent random variables. Therefore,\n", ">>>\n", ">>> $V(X_2) = pq + pq = 2pq$\n", "\n", "Nevertheless, we have dragged a hole: we believe by instinct that $X_2 = X_1 + X_1$, but in Science, our instinct is best supported by a proof.\n", "\n", "> **Exercise.** Prove that $X_2 = X_1 + X_1$. Notice that, according to our conceptualization, $ X_1 + X_1$ is very different than $2X_1$.\n", ">>\n", ">> **Answer.** Since $X_1$ is the function defined by\n", ">>>\n", ">>> $0 \\rightarrow q = 1-p$\n", ">>>\n", ">>> $1 \\rightarrow p$\n", ">>>\n", ">>> then, we can make the table of the sum of $X_1 + X_1$ if we keep in mind that the number of heads must be added while the probabilties must be multiplied because we are dealing with independent events. The Latex code of the sum is:\n", ">>>\n", ">>>> \\begin{center} \n", "\\begin{tabular}[c]{|c||c|c|} \\hline \n", "\\multicolumn{3}{|c|}{ Table of $Z = X_1 + X_1$} \\\\ \\hline \\hline \n", "$X_1$ ($\\downarrow)$ and $X_1$ ($\\rightarrow$)& 0 $\\rightarrow q$ & 1 $\\rightarrow p$ \\\\ \\hline \\hline \n", "0 $\\rightarrow q$ & 0 $\\rightarrow q^2$ & 1 $\\rightarrow qp$ \\\\ \\hline \n", "1 $\\rightarrow p$ & 1 $\\rightarrow pq$ & 2 $\\rightarrow p^2$ \\\\ \\hline \n", "\\end{tabular} \n", "\\end{center} \n", ">>> \n", ">>> This table roughly looks like this:\n", ">>>\n", ">>> $ ---- $& 0 $\\rightarrow q$ & 1 $\\rightarrow p$ \n", "0 $\\rightarrow q$ & 0 $\\rightarrow q^2$ & 1 $\\rightarrow qp$ \n", "1 $\\rightarrow p$ & 1 $\\rightarrow pq$ & 2 $\\rightarrow p^2$ \n", ">>\n", ">> Therefore, the random variable $X_1 + X_1$ is the function defined by\n", ">>\n", ">>> $0 \\rightarrow q^2$\n", ">>>\n", ">>> $1 \\rightarrow 2pq$\n", ">>>\n", ">>> $2 \\rightarrow p^2$\n", ">>>\n", ">> which coincides with $X_2$.\n", "\n", "\n", "> **Exercise.** Calculate the variance of $X_3$, the random variable whose probability assignments arise from the Binomial expansion with $n = 3$.\n", ">>\n", ">> **Answer.** The random variable $X_3$ represents a game with three coins in which we sum the results of 3 identical games. Thus, we claim that $X_3 = X_1 + X_1 + X_1 = X_2 + X_1$. Let us prove it. The Latex code of the corresponding table is:\n", ">>>\n", ">>> \\begin{center}\n", "\\begin{tabular}[c]{|c||c|c|} \\hline\n", "\\multicolumn{3}{|c|}{ Table of $Z = X_1 + X_1 + X_1 = X_2 + X_1$} \\\\ \\hline \\hline\n", "$X_2$ ($\\downarrow)$ and $X_1$ ($\\rightarrow$)& 0 $\\rightarrow q$ & 1 $\\rightarrow p$ \\\\ \\hline \\hline\n", "0 $\\rightarrow q^2$ & 0 $\\rightarrow q^3$ & 1 $\\rightarrow q^2p$ \\\\ \\hline\n", "1 $\\rightarrow 2qp$ & 1 $\\rightarrow 2q^2p$ & 2 $\\rightarrow 2qp^2$ \\\\ \\hline\n", "2 $\\rightarrow p^2$ & 2 $\\rightarrow qp^2$ & 3 $\\rightarrow p^3$ \\\\ \\hline \n", "\\end{tabular}\n", "\\end{center}\n", ">>\n", ">> This table looks like this:\n", ">>>\n", ">>> $X_2$ y $X_1$ & 0 $\\rightarrow q$ & 1 $\\rightarrow p$ \n", "0 $\\rightarrow q^2$ & 0 $\\rightarrow q^3$ & 1 $\\rightarrow q^2p$ \n", "1 $\\rightarrow 2qp$ & 1 $\\rightarrow 2q^2p$ & 2 $\\rightarrow 2qp^2$ \n", "2 $\\rightarrow p^2$ & 2 $\\rightarrow qp^2$ & 3 $\\rightarrow p^3$ \n", ">> \n", ">> Therefore, the random variable $X_1 + X_1 + X_1$ is the function defined by\n", ">>>\n", ">>> $0 \\rightarrow q^3$\n", ">>>\n", ">>> $1 \\rightarrow 3q^2p$\n", ">>>\n", ">>> $2 \\rightarrow 3qp^2$\n", ">>>\n", ">>> $3 \\rightarrow p^3$\n", ">>>\n", ">> which coincides with $X_3$.\n", ">>\n", ">> We see that the variance of $X_3$ is that of $X_1 + X_1 + X_1$, which is\n", ">>>\n", ">>> $V(X_3) = qp + qp + qp = 3pq$\n", "\n", "> **Exercise.** Generalize: find the variance of the Binomial Distribution in general.\n", "\n", ">> **Answer.** The variance of the Binomial Distribution with $n$ coins is $npq$. Proof by induction:\n", ">>>\n", ">>> Step 1: As it was found before, the variance of the Binomial Distribution with 1 coin is $pq$, so the formula is valid for this case.\n", ">>>\n", ">>> Step 2: Let us pass from step $k$ to step $k+1$. Let us suppose that the variance of $X_k$ is $kpq$. Then\n", ">>>\n", ">>>> $X_{k+1} = X_k + X_ 1$ and they are independent. Hence, the variance is\n", ">>>>\n", ">>>> $V(X_{k+1}) = V(X_k) + V(X_ 1) $\n", ">>>>> $= kpq + pq = (k+1)pq$\n", ">>>>\n", ">>>> such as demanded.\n", ">>>\n", ">>> In summary, the formula is valid for $k = 1$, and if we assume that it is valid for $k$, we can prove that it is valid for $k+1$. As a consequence, it is valid for all $k$.\n", "\n", "\n", "> **Exercise.** Let $X$ be the random variable that counts heads in a Binomial Distribution with $n$ coins and probability of heads $p$. Show that the distribution of its relative frequency $X/n$ shrinks or converges towards the Expected Value when $n$ grows.\n", ">>\n", ">> **Answer.** The Expected Value of $X$ is $np$ and its variance is $V(X) = npq$. The Relative Frequency of a Binomial Distribution $X$ with $n$ coins is $X/n$. Thus, the Expected Value of the Relative Frequency $X/n$ is\n", ">>>\n", ">>> $E(X/n) = \\frac{E(X)}{n} = \\frac{np}{n} = p$\n", ">>>\n", ">> and its variance, or dispersion with respect to its Expected Value $p$, is\n", ">>>\n", ">>> $V(X/n) = \\frac{V(X)}{n^2} = \\frac{npq}{n^2} = \\frac{pq}{n} \\rightarrow 0$\n", ">>>\n", ">> This convergence to zero explains why the Relative Frequency of a Binomial Distribution tends to $p$ and shrinks towards it as the number of coins grows.\n", "\n", "\n", " \n", "> **Exercise.** Explain the meaning of the Sample Mean $\\overline{Y_n}$ of a general random variable $Y$. Give an example.\n", ">>\n", ">> **Answer.** A random variable $Y$ can be imagined as a source of random numbers, with the distribution dictated by $Y$. The sample mean consists of selecting a given number $n$, taking $n$ data from the source, adding them together and averaging over $n$. Next, we must repeat this task as many times as possible to get an estimation of the corresponding distribution.\n", ">>>\n", ">>> For instance, if $Y$ is the random variable that measures the height of adult persons, it takes values from 1m to 2.8m. The sample mean with $n = 30$ consists of taking 30 persons at random, measuring their heights and averaging. This corresponds with 1 value of the sample mean. Next, we must repeat this procedure over and over to gain insight on the shape of the corresponding distribution. That is why we define the *Sample Mean* of $Y$ over random samples of size $n$ as\n", ">>>\n", ">>>> $\\overline{Y_n} = \\frac{\\sum_{i = 1}^n Y_i}{n}$\n", ">>>\n", ">>> where all the $Y_i$ have the same distribution as $Y$ and are independent.\n", ">>>\n", ">>> A plausible theory of height says that beginning from 1m, each additional height of 10cm is caused by an active copy of a gen of growth, which can be on state on else off with equal probability. Thus, if you have 7 active copies, your height is 1.7m. Then, the total height is dictated by $Y = 1 + 0.1X$ where $X$ is the random variable with a binomial distribution with $18$ coins and $p = 0.5$. The sample mean $Y_n$ becomes for this case:\n", ">>>>\n", ">>>> $\\overline{Y_n} = \\frac{\\sum_{i = 1}^n (1 + 0.1X)}{n}$\n", ">>>>\n", ">>> and each $X$ is a Binomial $Bi(18, 0.5)$,\n", "\n", "\n", "> **Example.** Let us prove that the distribution of the Sample Mean of a Binary Random Variable shrinks towards the mean.\n", ">>\n", ">> Proof. We use the variance as a measure of dispersion of distributions. A Binary Random Variable $Y = X_m$, with $m$ coins, has mean $mp$ and variance $mpq$. Let $\\overline{Y_n}$ be its sample mean with $n$ data:\n", ">>>\n", ">>> $\\overline{Y_n} = \\frac{\\sum_{i = 1}^n Y_i}{n}$\n", ">>>> where each $Y_i$ is a copy of $Y$.\n", ">>\n", ">> We assume that the source of random numbers remains the same no matter how many data has generated. This implies that the $Y_i$ are independent. For instance, this happens for a Binary random variable when the coin does not tear out after playing with it. Then\n", ">>>\n", ">>> $V(\\overline{Y_n}) = V(\\frac{\\sum_{i = 1}^n Y_i}{n}) $\n", ">>>>\n", ">>>> $= \\frac{1}{n^2} V( \\sum_{i = 1}^n Y_i ) $\n", ">>>>\n", ">>>> $= \\frac{1}{n^2} \\sum_{i = 1}^n V(Y_i )$\n", ">>>>\n", ">>>> $ = \\frac{1}{n^2} \\sum_{i = 1}^n mpq $\n", ">>>> \n", ">>>> $ = \\frac{n}{n^2} mpq $\n", ">>>>\n", ">>>> $ = \\frac{mpq}{n} \\rightarrow 0$\n", ">>>>\n", ">> This shows that the Sample Mean shrinks towards the Expected Value when $n$ diverges towards infinite.\n", "\n", "> **Exercise.** Generalize the previous result to whatever random variable $X$. Does $X$ must observe some conditions, or is this result unconditional?\n", ">>\n", ">> **Answer.** Let $X$ be a random variable. Let us prove that when $n$ grows, the distribution of every sample mean, $\\overline{X_n}$, shrinks towards the expected Value of $X$.\n", ">>\n", ">> The variance of a random variable measures the averaged quadratic deviation of the variable with respect to its mean or, more exactly, the Expected Value. Thus, to prove that when $n$ grows, the distribution of every sample mean, $\\overline{X_n}$, shrinks towards the Expected Value, it suffices to prove that the variance diminishes and tends to zero as $n$ grows.\n", ">>\n", ">> Let $\\sigma^2$ be the variance of $Y$. The variance $V$ of a Sample Mean $\\overline{Y_n}$ is\n", ">>>\n", ">>> $V(\\overline{Y_n}) = V(\\frac{\\sum_{i = 1}^n Y_i}{n}) $\n", ">>>>\n", ">>>> $ = \\frac{1}{n^2} \\sum_{i = 1}^n V(Y_i ) $\n", ">>>>\n", ">>>> $ = \\frac{1}{n^2} \\sum_{i = 1}^n \\sigma^2 $\n", ">>>>\n", ">>>> $= \\frac{1}{n^2} n \\sigma^2 $\n", ">>>>\n", ">>>> $ = \\frac{\\sigma^2}{n} \\rightarrow 0$\n", ">>>>\n", ">>> This is valid whenever the variance of $X$ is finite. No other condition is required. This reinforces the idea that the higher is the number of random data, much better is the assessment of the mean value of any random variable. How strong is this result! It is an astonishing example of what means to do Science when randomness is the air that we breathe.\n", "\n", "\n", "# GRADUATION\n", "\n", "Choose one of the three following projects:\n", ">\n", "> 1. Reformulate the Scientific Algorithm for comparison between what one thinks and what one observes in terms of the arising discrepancy that must be measured in appropriate units.\n", ">\n", "> 2. Program `relativeFrequencyTendency2.py` draws the tendency of the relative frequency of heads when a coin is flipped an increasing number of times. This tendency appears as a path that dances around the red line that points to the expected value, `0.5`. Now, if we repeat the experiment, the same path will not be found. And if we repeat the experiment many times, for each number of flips or coins, a distribution will appear. Develop the code for a realizable 3D-graphic of that distribution, describe its mode of convergence, and compare your results with mathematical expectations. \n", ">\n", "> 3. We have observed in our simulations that for a Binomial Distribution, the sample mean converges to the expected mean or, more exactly, to the Expected Value. We have been able of explaining this fact with the help of the variance. But other approaches exist that group themselves around a mathematical result with wide coverage that is known as the Law of Large Numbers. Consult its conditions and proof, and decipher it.\n", "\n", "\n", "# CONCLUSIONS\n", "\n", "\n", "> **Exercise.** Express your personal conclusions.\n", ">>\n", ">> **Answer.** Our personal conclusions are:\n", "\n", "- It has been a pleasure to work with Python, Colab, Gemini, Codey and its AI-associates.\n", "\n", "- The purpose of the Scientific Method is\n", "to compare what one thinks with what one observes.\n", "\n", "- The Scientific Method to refine what we think\n", "about nature is based on common sense:\n", "If we find something rare or strange,\n", "it is best explained by a point of view different from ours.\n", "'To be explained' is a terminology used to say that\n", "some data and the offered explanations\n", "seem to be compatible to a good degree.\n", "That does not mean that\n", "the explanation is the correct and\n", "the only one of the given data.\n", "Many theories or explanations can exist.\n", "But if the field of research matures,\n", "gathering more and more data under different settings,\n", "the possible compatible theories tend to diminish.\n", "\n", "- A terrible problem in Science is that\n", "it is rare to see clearly defined, deterministic events.\n", "Rather, we are always surrounded by randomness\n", "that causes every kind of interference and failures.\n", "That is why we face the question:\n", "How can we take decisions when there is randomness?\n", "The answer, based on common sense, is:\n", "look for tendencies. If none exists, create one.\n", "Examples of tendencies:\n", "The relative frequency that tends to a probability\n", "when the number of trials grows.\n", "The mean of a sample that tends to the mean of the population\n", "when the sample is random and grows.\n", "As the number of repetitions increases,\n", "the distribution of a distinguished event\n", "appearing k times tends toward a limit that coincides with the mathematical expectation.\n", "\n", "- Using these tendencies, we can take decisions.\n", "But, beware, decisions are personal and subjective.\n", "This is clearly seen when\n", "the significance, the proportion of rare events, is defined.\n", "The default value is 0.05, and lower values are welcome\n", "because to be rare means a right to formulate new hypotheses.\n", "So, we adhere to the insight that says that\n", "accepting a theory in Science does not convert it into truth,\n", "rather, it only points to a compatibility between\n", "data and explanation, observation and theory.\n", "\n", "- We have seen how the whole machinery functions\n", "over the Binomial Distribution.\n", "This notion is very general.\n", "The very idea is to gather much data to estimate\n", "the distribution of the relative frequency of, say,\n", "heads given a null hypothesis, for instance,\n", "the probability of heads equals that of tails.\n", "Next, we use that distribution to define\n", "what is normal, ordinary, and what is abnormal, extraordinary,\n", "according to the given null hypothesis.\n", "Ordinary events are compatible with the null hypothesis\n", "and invite one to accept it.\n", "Extraordinary events\n", "are incompatible with the null hypothesis\n", "and invite one to reject it.\n", "\n", "- The very foundation of science in the presence of randomness is the existence of perfectly defined theoretical distributions as the limit of the observed distributions driven by randomness as the sample size grows toward infinity.\n", "\n", "- Contrary to the humility of the Scientific Method,\n", "we, human beings, love to fall in love with theories, dogmas,\n", "revelations, intuitions.\n", "Because of this, we deserve to be called irrational.\n", "That is why Real Science is an irrational affair.\n" ] }, { "cell_type": "markdown", "source": [ "# To know more\n", "\n", "- Other ideas like this, from the Author, that cover other topics of statistics:\n", " https://www.ejristos.com/dl/javaevo5Pv3.zip\n", "\n", "- Applications to ordinary life and social sciences, in Spanish:\n", " https://www.ejristos.com/dl/EstVol1V18.pdf\n", "\n", "\n", "## WRITE SOMETHING BETTER THAN THIS\n", "\n", "Modify this material to your liking until you achieve something that truly satisfies you. To do this, reuse this material as much as possible by downloading an editable version.\n", "\n", "This publishing system is based on JupyterLab notebooks, which allow documents containing text cells and Python code cells to be run with one click.\n", "\n", "It's all extremely easy with Google Colab:\n", "\n", "## Google Colab\n", "\n", "\"For Colab, you can simply visit the notebook on the Colab website, https://colab.google/, and click \"Connect.\" The caveat is that it requires an active Google account. Then, just click the \"Connect\" button in the top right corner, and you can start running the code here.\"\n", "\n", "Visit https://colab.google/notebooks/ to find wonderful examples of notebooks written with this technology. These notebooks are not editable and open in read-only mode. This document has editable versions at https://ejristos.com. This version must be uploaded to Colab, and if it opens in read-only mode, editing mode is activated by double-clicking on the document. You can write text or code by selecting a corresponding cell, whose menu is located in the horizontal bar and also at the end of each cell.\n", "\n", "Formatting is done after finishing writing the document. Markdown commands are used for this. https://www.markdownguide.org/basic-syntax/. They are extremely simple and perform poorly if you use many levels of concatenation and also with some few symbols in LaTeX. Use simple formatting and everything will be fine.\n", "\n" ], "metadata": { "id": "jlf4XbY40yFk" } } ], "metadata": { "colab": { "provenance": [] }, "kernelspec": { "display_name": "Python 3", "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.12.3" } }, "nbformat": 4, "nbformat_minor": 0 }