<jupyter_start><jupyter_text>Day 12: Intro to Probability---- motivation: variable material properties- probability density function- cumulative distribution function- quantile function Important Caveat: This class does **not** satisfy the prob/stat requirement...Before we start, let me give you a warning:> I'm going to teach you useful content in this class. However, I'm going to gloss over the whole subject of statistics. So long as you're dealing not dealing with safety-critical systems, you can get far with what I'm going to teach. If you're working on safety-critical systems someday, make sure to take a class on statistics, or better yet on reliability engineering.<jupyter_code>import pandas as pd 
import grama as gr 
import numpy as np
from scipy.optimize import minimize
from plotnine import *
DF = gr.Intention()

## Dataset of steel alloy specimens
from grama.data import df_ruff<jupyter_output><empty_output><jupyter_text>Motivating Example: Material properties for designWhen doing mechanical design, we need "material properties":| Property | Symbol ||---|---|| Density | $\rho$ || Elasticity | $E$ || Strength (yield) | $\sigma_y$ || Strength (ultimate) | $\sigma_u$ |We can often find "material properties" via Google search; for instance, you can find the following page by searching for "aluminum 6061 strength":[MatWeb aluminum page](http://www.matweb.com/search/datasheet_print.aspx?matguid=1b8c06d0ca7c456694c7777d9e10be5b)But why does this page have the scary warning:> Data points with the AA note have been provided by the Aluminum Association, Inc. and are NOT FOR DESIGN.*Probability* is **key** to understanding why!  Thought Experiment: Tensile testingTensile testing gives us information on the *strength* of materials: ![tensile specimens](./images/Tensile_specimen-round_and_flat.jpg)Via Wikimedia Commons. Questions:**Elasticity**- (Single) Imagine we measure the elasticity $E$ for a *single* specimen multiple times. Do you expect to see the *exact* same value each time? Why or why not?- (Double) Imagine we measure the elasticity $E_i$ for two *different* specimens $i=1, 2$. Do you expect to see the *exact* same value each time? Why or why not?- Are your answers to "why or why not?" the same for (Single) vs (Double)?**Strength**- Imagine we measure the strength (ultimate) $\sigma_u$ for a *single* specimen. Can we do this multiple times?- Imagine we measure the strength (ultimate) $\sigma_{u,i}$ for two *different* specimens $i=1, 2$. Do you expect to see the *exact* same value each time? Why or why not? **Physical mechanisms for variability**- testing equipment  - rounding error- part  - geometry- material  - manufacturing method  - crystal structure is random!  Real Experiments: Real materials dataEvery individual specimen gives us a *different* material property value, for instance the following is a set of tests on cast steel tensile yield strength (TYS, ksi) and tensile ultimate strength (TUS, ksi), along with part thicknesses (in) and part identifiers.<jupyter_code>df_ruff<jupyter_output><empty_output><jupyter_text>Q: What do you observe based on the following plot?<jupyter_code>(
    df_ruff
    >> gr.tf_mutate(idx=DF.index)
    
    >> ggplot(aes(y="idx", x="TYS"))
    + geom_point()
    + labs(
        y="Part Index",
        x="Yield Strength (ksi)"
    )
)<jupyter_output><empty_output><jupyter_text>Let's imagine "stacking" all of the dots:<jupyter_code>(
    df_ruff
    >> ggplot(aes("TYS"))
    + geom_dotplot(bins=30)
    
    + labs(
        y="Count",
        x="Yield Strength (ksi)"
    )
)<jupyter_output><empty_output><jupyter_text>We could instead count the dots in bins; this is a histogram:<jupyter_code>(
    df_ruff
    >> ggplot(aes("TYS"))
    + geom_histogram(aes(y="stat(density)"), bins=30)
    + geom_dotplot(bins=30)
    
    + labs(
        y="Count",
        x="Yield Strength (ksi)"
    )
)<jupyter_output><empty_output><jupyter_text>Using materials data for designImagine we want to design parts that will fail at a rate of 1 in 2.  Q: What *single value* of Yield Strength (`TYS`) should we use to get failures at a rate of 1 in 2, based on the data below?<jupyter_code>(
    df_ruff
    >> ggplot(aes("TYS"))
    + geom_histogram(aes(y="stat(density)"), bins=30)
    + geom_dotplot(bins=30)
    
    + labs(
        y="Count",
        x="Yield Strength (ksi)"
    )
)<jupyter_output><empty_output><jupyter_text>Q: What if we want failure at a rate of 1 in 100? What value of Yield Strength (`TYS`) would you use?  Random Variables---What is a random variable? I'm going to give you two viewpoints:> 1. (Theoretical) A random variable __has__ a *function* that describes *relative frequencies*.> > 2. (Practical) A random variable __is__ a variable that has *no single, well-defined value*.[Theoretically](https://en.wikipedia.org/wiki/Probability_density_functionExample), a random variable is defined by a function that precisely defines how frequently particular values tend to occur. Practically, we will use random variables to model *uncertain* quantities. The Grama function `gr.marg_named(data, dist)` fits a random variable based on two pieces of information:- `data` is a dataset used to inform the model- `dist` is a string used to select a particular distributionIt is common to model strength using a lognormal (`"lognorm"`) distribution; let's do that:<jupyter_code>## Fit a marginal
marg_tys = gr.marg_named(df_ruff.TYS, "lognorm")<jupyter_output><empty_output><jupyter_text>Probability Density Function (PDF)We can access the PDF via dot `.` notation of a marginal with the function `marg.l(x)`. It's `l` because a probability density is sometimes called a *likelihood*.The probability density function describes the relative likelihood of different values. It is a function that looks a lot like a histogram:<jupyter_code>## Evaluate marginal at many points
df_density = gr.df_make(TYS=np.linspace(150, 180))
df_density["l"] = list(map(marg_tys.l, df_density.TYS))

## Visualize
(
    df_ruff
    >> ggplot(aes("TYS"))
    + geom_histogram(aes(y="stat(density)"), bins=30)
    + geom_line(
        data=df_density,
        mapping=aes(y="l")
    )
    
    + labs(
        y="Density",
        x="Yield Strength (ksi)"
    )
)<jupyter_output><empty_output><jupyter_text>- Values of TYS where the Density is smaller will tend to occur less frequently- Values of TYS where the Density is larger will tend to occur more frequently- The *shape* of this distribution is asymmetric; the right tail (extreme right side) stretches much further than the left Q: How frequently would you expect `TYS == 150` to occur? Q: *Roughly*, what is the most-likely TYS value?<jupyter_code>## Point of max PDF
res = minimize(
    lambda x: -1 * marg_tys.l(x),
    155
)
tys_max_pdf = res.x[0]
tys_max_pdf<jupyter_output><empty_output><jupyter_text>Q: How often should we expect to see `TYS == 156.28305687066222` in practice?  ProbabilityWeirdly, for a continuous distribution (a PDF), we don't expect to see any *specific* value. Instead, we tend to think in terms of **areas** of a PDF.<jupyter_code>## Compute median
tys_median = marg_tys.q(0.5)
print("TYS median ~ {0:4.3f}".format(tys_median))

## Evaluate marginal at many points
df_density = gr.df_make(TYS=np.linspace(150, 180, num=200))
df_density["l"] = list(map(marg_tys.l, df_density.TYS))

## Visualize
(
    df_ruff
    >> ggplot(aes("TYS"))
    + geom_ribbon(
        data=df_density
        >> gr.tf_filter(DF.TYS <= tys_median),
        mapping=aes(ymax="l", ymin=0),
        fill="salmon",
        alpha=1/3,
    )
    + geom_line(
        data=df_density,
        mapping=aes(y="l")
    )
    
    + labs(
        y="Density",
        x="Yield Strength (ksi)"
    )
)<jupyter_output>TYS median ~ 157.167<jupyter_text>The median is the value that divides the PDF in "half"; we can use the median to give us a strength value at which parts will break at a rate of 1 in 2, aka 50%, aka $0.5$. Since probability is defined in terms of an area, when we think about probability we have to think in terms of **sets** (not single values).$$\mathbb{P}[X \in A] = \int_{x \in A} \text{PDF}(x) dx$$We'll come back to this definition in our next class.  Cumulative Distribution FunctionThe *indefinite integral* of the PDF is called the *cumulative distribution function* (CDF):$$\text{CDF}(x) = \int_{-\infty}^x \text{PDF}(x') dx'$$> The CDF tells us the area to the left of a given point $x$We can access the CDF of a marginal with dot `.` notation via `marg.p(x)`. Here `p` stands for probability; unlike the PDF, the CDF does return a probability:<jupyter_code>## Evaluate the CDF at the median
# marg_tys.p(150)
marg_tys.p(tys_median)<jupyter_output><empty_output><jupyter_text>Above, we saw that the median was the point that cuts the PDF "in half"---this CDF calculation verifies that assertion. Here we see that `tys_median` has `0.5` of the PDF below, and thus `0.5` of the PDF above. Connection to "NOT FOR DESIGN"Remember this scary statement from the [MatWeb aluminum page](http://www.matweb.com/search/datasheet_print.aspx?matguid=1b8c06d0ca7c456694c7777d9e10be5b)?> Data points with the AA note have been provided by the Aluminum Association, Inc. and are NOT FOR DESIGN.The median (also the mean) are sometimes used for *typical values*; the AA values for strength reported are typical values. This means, if you design parts using those values (and load to designed failure), you will tend to get a failure rate of 50%. This probably isn't what you want in most cases!  *Aside*: Upper-side of distributionWe can also get the area for values above a threshold $x$ by taking $1 - \text{CDF}(x)$. For instance, here's an area on the right-side of a distribution:<jupyter_code>## Compute median
tys_q75 = float(marg_tys.q(0.75))

## Evaluate marginal at many points
df_density = gr.df_make(TYS=np.linspace(150, 180, num=200))
df_density["l"] = list(map(marg_tys.l, df_density.TYS))

## Visualize
(
    df_ruff
    >> ggplot(aes("TYS"))
    + geom_ribbon(
        data=df_density
        >> gr.tf_filter(tys_q75 <= DF.TYS),
        mapping=aes(ymax="l", ymin=0),
        fill="salmon",
        alpha=1/3,
    )
    + geom_line(
        data=df_density,
        mapping=aes(y="l")
    )
    
    + labs(
        y="Density",
        x="Yield Strength (ksi)"
    )
)<jupyter_output><empty_output><jupyter_text>We can compute this area with the CDF using 1 - CDF:<jupyter_code>## Evaluate the CDF at the median
# 1 - CDF(q75)
1 - marg_tys.p(tys_q75)<jupyter_output><empty_output><jupyter_text>Apparently this upper shaded region is 25% of the PDF.  Quantile FunctionThe *inverse* of the CDF is called the *quantile function*:$$\text{Q}(p) = \text{CDF}^{-1}(p)$$Note that the quantile function takes in a *probability*; it returns a specific value of the associated random variable. We can use this to find values that bracket a specified fraction of a PDF.  Using the quantile function for designWhat if we want to find a strength value where parts will break at a rate of 1 in 100? We can use the quantile function to compute such a value.We can access the quantile function of a marginal in Grama using `marg.q(p)`; remember that the quantile function takes a probability and returns a random variable value.<jupyter_code>## Compute lower quantile as design value
tys_lower = marg_tys.q(0.01)
print("TYS lower ~ {0:4.3f}".format(tys_lower))

## Evaluate marginal at many points
df_density = gr.df_make(TYS=np.linspace(150, 180, num=200))
df_density["l"] = list(map(marg_tys.l, df_density.TYS))

## Visualize
(
    df_ruff
    >> ggplot(aes("TYS"))
    + geom_ribbon(
        data=df_density
        >> gr.tf_filter(DF.TYS <= tys_lower),
        mapping=aes(ymax="l", ymin=0),
        fill="salmon",
        alpha=1/3,
    )
    + geom_line(
        data=df_density,
        mapping=aes(y="l")
    )
    + geom_vline(
        xintercept=tys_lower,
        linetype="dashed",
        color="salmon",
    )
    
    + labs(
        y="Density",
        x="Yield Strength (ksi)"
    )
)<jupyter_output>TYS lower ~ 154.698
{'fraction_non_alphanumeric': 0.10613831899404368, 'fraction_numerical': 0.0164626075446724, 'mean_word_length': 4.155223880597015, 'pattern_counts': {'":': 2, '<': 58, '<?xml version=': 0, '>': 86, 'https://': 1, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '14', 'dir': 'jupyter-structured-clean-dedup', 'id': '555059', 'license': 'MIT', 'n_tokens_mistral': 4075, 'n_tokens_neox': 3772, 'n_words': 1405, 'path': 'day12.ipynb', 'repo_name': 'zdelrosario/pdo-notes', 'size': '12088'}