Newer
Older
{
"cells": [
{
"cell_type": "markdown",
"metadata": {
"toc": true
},
"source": [
"<h1>Table of Contents<span class=\"tocSkip\"></span></h1>\n",
"<div class=\"toc\"><ul class=\"toc-item\"><li><span><a href=\"#Data-sets\" data-toc-modified-id=\"Data-sets-1\"><span class=\"toc-item-num\">1 </span>Data sets</a></span><ul class=\"toc-item\"><li><span><a href=\"#Synthetic-data-with-unobservables\" data-toc-modified-id=\"Synthetic-data-with-unobservables-1.1\"><span class=\"toc-item-num\">1.1 </span>Synthetic data with unobservables</a></span></li><li><span><a href=\"#Data-without-unobservables\" data-toc-modified-id=\"Data-without-unobservables-1.2\"><span class=\"toc-item-num\">1.2 </span>Data without unobservables</a></span></li></ul></li><li><span><a href=\"#Algorithms\" data-toc-modified-id=\"Algorithms-2\"><span class=\"toc-item-num\">2 </span>Algorithms</a></span><ul class=\"toc-item\"><li><span><a href=\"#Contraction-algorithm\" data-toc-modified-id=\"Contraction-algorithm-2.1\"><span class=\"toc-item-num\">2.1 </span>Contraction algorithm</a></span></li><li><span><a href=\"#Causal-approach---metrics\" data-toc-modified-id=\"Causal-approach---metrics-2.2\"><span class=\"toc-item-num\">2.2 </span>Causal approach - metrics</a></span></li></ul></li><li><span><a href=\"#Performance-comparison\" data-toc-modified-id=\"Performance-comparison-3\"><span class=\"toc-item-num\">3 </span>Performance comparison</a></span><ul class=\"toc-item\"><li><span><a href=\"#With-unobservables-in-the-data\" data-toc-modified-id=\"With-unobservables-in-the-data-3.1\"><span class=\"toc-item-num\">3.1 </span>With unobservables in the data</a></span></li><li><span><a href=\"#Without-unobservables\" data-toc-modified-id=\"Without-unobservables-3.2\"><span class=\"toc-item-num\">3.2 </span>Without unobservables</a></span></li></ul></li></ul></div>"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"Our model is defined by the probabilistic expression \n",
"\n",
"P(Y=0 | \\text{do}(R=r)) = \\sum_x \\underbrace{P(Y=0|X=x, T=1)}_\\text{1} \n",
"\\overbrace{P(T=1|R=r, X=x)}^\\text{2} \n",
"\\underbrace{P(X=x)}_\\text{3}\n",
"\\end{equation}\n",
"\n",
"which is equal to \n",
"\n",
"\\begin{equation}\\label{model_cont}\n",
"P(Y=0 | \\text{do}(R=r)) = \\int_x P(Y=0|X=x, T=1)P(T=1|R=r, X=x)P(X=x)\n",
"\\end{equation}\n",
"\n",
"for continuous $x$. In the model Z is a latent, unobserved variable, and can be excluded from the expression with do-calculus by showing that $X$ is admissible for adjustment. Model as a graph:\n",
"![Model as picture](../figures/intervention_model.png \"Intervention model\")\n",
"\n",
"For predicting the probability of negative outcome the following should hold because by Pearl $P(Y=0 | \\text{do}(R=r), X=x) = P(Y=0 | R=r, X=x)$ when $X$ is an admissible set:\n",
"\n",
"\\begin{equation} \\label{model_pred}\n",
"P(Y=0 | \\text{do}(R=r), X=x) = P(Y=0|X=x, T=1)P(T=1|R=r, X=x).\n",
"\\end{equation}\n",
"\n",
"Still it should be noted that this prediction takes into account the probability of the individual to be given a positive decision ($T=1$), see second term in \\ref{model_pred}.\n",
"\n",
"----\n",
"\n",
"### Notes\n",
"\n",
"* Equations \\ref{model_disc} and \\ref{model_cont} describe the whole causal effect in the population (the causal effect of changing $r$ over all strata $X$).\n",
"* Prediction should be possible with \\ref{model_pred}. Both terms can be learned from the data. NB: the probability $P(Y=0 | \\text{do}(R=r), X=x)$ is lowest when the individual $x$ is the most dangerous or the least dangerous. How could we infer/predict the counterfactual \"what is the probability of $Y=0$ if we were to let this individual go?\" has yet to be calculated.\n",
"* Is the effect of R learned/estimated correctly if it is just plugged in to a predictive model (e.g. logistic regression)? **NO**\n",
"* $P(Y=0 | do(R=0)) = 0$ only in this application. <!-- My predictive models say that when $r=0$ the probability $P(Y=0) \\approx 0.027$ which would be a natural estimate in another application/scenario (e.g. in medicine the probability of an adverse event when a stronger medicine is distributed to everyone. Then the probability will be close to zero but not exactly zero.) -->\n",
"\n",
"Imports and settings."
"metadata": {},
"outputs": [],
"source": [
"# Imports\n",
"\n",
"import numpy as np\n",
"import pandas as pd\n",
"from datetime import datetime\n",
"import matplotlib.pyplot as plt\n",
"import scipy.stats as scs\n",
"import scipy.integrate as si\n",
"import seaborn as sns\n",
"import numpy.random as npr\n",
"from sklearn.preprocessing import OneHotEncoder\n",
"from sklearn.linear_model import LogisticRegression\n",
"from sklearn.ensemble import RandomForestClassifier\n",
"from sklearn.model_selection import train_test_split\n",
"\n",
"# Settings\n",
"\n",
"%matplotlib inline\n",
"\n",
"plt.rcParams.update({'font.size': 16})\n",
"plt.rcParams.update({'figure.figsize': (14, 7)})\n",
"\n",
"# Suppress deprecation warnings.\n",
"\n",
"import warnings\n",
"\n",
"def fxn():\n",
" warnings.warn(\"deprecated\", DeprecationWarning)\n",
"\n",
"with warnings.catch_warnings():\n",
" warnings.simplefilter(\"ignore\")\n",
" fxn()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Data sets\n",
"\n",
"### Synthetic data with unobservables\n",
"\n",
"In the chunk below, we generate the synthetic data as described by Lakkaraju et al. The default values and definitions of $Y$ and $T$ values follow their description.\n",
"\n",
"**Parameters**\n",
"\n",
"* M = `nJudges_M`, number of judges\n",
"* N = `nSubjects_N`, number of subjects assigned to each judge\n",
"* betas $\\beta_i$ = `beta_i`, where $i \\in \\{X, Z, W\\}$ are coefficients for the respected variables\n",
"\n",
"**Columns of the data:**\n",
"\n",
"* `judgeID_J` = judge IDs as running numbering from 0 to `nJudges_M - 1`\n",
"* R = `acceptanceRate_R`, acceptance rates\n",
"* X = `X`, invidual's features observable to all (models and judges)\n",
"* Z = `Z`, information observable for judges only\n",
"* W = `W`, unobservable / inaccessible information\n",
"* T = `decision_T`, bail-or-jail decisions where $T=0$ represents jail decision and $T=1$ bail decision.\n",
"* Y = `result_Y`, result variable, if $Y=0$ person will or would recidivate and if $Y=1$ person will or would not commit a crime.\n",
"\n",
"The generated data will have M\\*N rows."
" '''Return value of sigmoid function (inverse of logit) at x.'''\n",
"\n",
"def dataWithUnobservables(nJudges_M=100,\n",
" nSubjects_N=500,\n",
" beta_X=1.0,\n",
" beta_Z=1.0,\n",
"\n",
" # Assign judge IDs as running numbering from 0 to nJudges_M - 1\n",
" df = df.assign(judgeID_J=np.repeat(range(0, nJudges_M), nSubjects_N))\n",
"\n",
" # Sample acceptance rates uniformly from a closed interval\n",
" # from 0.1 to 0.9 and round to tenth decimal place.\n",
" acceptance_rates = np.round(npr.uniform(.1, .9, nJudges_M), 10)\n",
"\n",
" # Replicate the rates so they can be attached to the corresponding judge ID.\n",
" df = df.assign(acceptanceRate_R=np.repeat(acceptance_rates, nSubjects_N))\n",
"\n",
" # Sample the variables from standard Gaussian distributions.\n",
" df = df.assign(X=npr.normal(size=nJudges_M * nSubjects_N))\n",
" df = df.assign(Z=npr.normal(size=nJudges_M * nSubjects_N))\n",
" df = df.assign(W=npr.normal(size=nJudges_M * nSubjects_N))\n",
" probabilities_Y = sigmoid(beta_X * df.X + beta_Z * df.Z + beta_W * df.W)\n",
" df = df.assign(probabilities_Y=probabilities_Y)\n",
"\n",
" # Result is 0 if P(Y = 0| X = x; Z = z; W = w) >= 0.5 , 1 otherwise\n",
" df = df.assign(result_Y=np.where(df.probabilities_Y >= 0.5, 0, 1))\n",
" # For the conditional probabilities of T we add noise ~ N(0, 0.1)\n",
" probabilities_T = sigmoid(beta_X * df.X + beta_Z * df.Z)\n",
" probabilities_T += np.sqrt(0.1) * npr.normal(size=nJudges_M * nSubjects_N)\n",
" df = df.assign(probabilities_T=probabilities_T)\n",
" # Sort by judges then probabilities in decreasing order\n",
" # Most dangerous for each judge are at the top.\n",
" df.sort_values(by=[\"judgeID_J\", \"probabilities_T\"],\n",
" ascending=False,\n",
" inplace=True)\n",
" # Iterate over the data. Subject will be given a negative decision\n",
" # if they are in the top (1-r)*100% of the individuals the judge will judge.\n",
" # I.e. if their within-judge-index is under 1 - acceptance threshold times\n",
" # the number of subjects assigned to each judge they will receive a\n",
" # negative decision.\n",
" df['decision_T'] = np.where((df.index.values % nSubjects_N) <\n",
" ((1 - df['acceptanceRate_R']) * nSubjects_N),\n",
" 0, 1)\n",
" # Halve the data set to test and train\n",
" train, test = train_test_split(df, test_size=0.5)\n",
" train_labeled = train.copy()\n",
" test_labeled = test.copy()\n",
" # Set results as NA if decision is negative.\n",
" train_labeled.loc[train_labeled.decision_T == 0, 'result_Y'] = np.nan\n",
" test_labeled.loc[test_labeled.decision_T == 0, 'result_Y'] = np.nan\n",
" return train_labeled, train, test_labeled, test, df"
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Data without unobservables\n",
"\n",
"In the chunk below, we generate a simplified data. The default values and definitions of $Y$ and $T$ values follow the previous description.\n",
"\n",
"**Parameters**\n",
"\n",
"* M = `nJudges_M`, number of judges\n",
"* N = `nSubjects_N`, number of subjects assigned to each judge\n",
"\n",
"**Columns of the data:**\n",
"\n",
"* `judgeID_J` = judge IDs as running numbering from 0 to `nJudges_M - 1`\n",
"* R = `acceptanceRate_R`, acceptance rates\n",
"* X = `X`, invidual's features observable to all (models and judges), now $X \\sim \\mathcal{N}(0, 1)$\n",
"* T = `decision_T`, bail-or-jail decisions where $T=0$ represents jail decision and $T=1$ bail decision.\n",
"* $p_y$ = `probabilities_Y`, variable where $p_y = P(Y=0)$\n",
"* Y = `result_Y`, result variable, if $Y=0$ person will or would recidivate and if $Y=1$ person will or would not commit a crime. Here $Y \\sim \\text{Bernoulli}(\\frac{1}{1+exp\\{-\\beta_X \\cdot X\\}})$"
"def dataWithoutUnobservables(nJudges_M=100,\n",
" nSubjects_N=500,\n",
" sigma=0.0):\n",
"\n",
" df = pd.DataFrame()\n",
"\n",
" # Assign judge IDs as running numbering from 0 to nJudges_M - 1\n",
" df = df.assign(judgeID_J=np.repeat(range(0, nJudges_M), nSubjects_N))\n",
"\n",
" # Sample acceptance rates uniformly from a closed interval\n",
" # from 0.1 to 0.9 and round to tenth decimal place.\n",
" acceptance_rates = np.round(npr.uniform(.1, .9, nJudges_M), 10)\n",
"\n",
" # Replicate the rates so they can be attached to the corresponding judge ID.\n",
" df = df.assign(acceptanceRate_R=np.repeat(acceptance_rates, nSubjects_N))\n",
"\n",
" # Sample feature X from standard Gaussian distribution, N(0, 1).\n",
" df = df.assign(X=npr.normal(size=nJudges_M * nSubjects_N))\n",
"\n",
" # Calculate P(Y=1|X=x) = 1 / (1 + exp(-beta_X * x)) = sigmoid(beta_X * x)\n",
" # Draw Y ~ Bernoulli(sigmoid(beta_X * x)) = Bin(1, p)\n",
" size=nJudges_M * nSubjects_N)\n",
"\n",
" df = df.assign(result_Y=results)\n",
" # Invert the probabilities. P(Y=0 | X) = 1 - P(Y=1 | X)\n",
" df.probabilities_Y = 1 - df.probabilities_Y\n",
"\n",
" # Assign the prediction probabilities and add some Gaussian noise\n",
" # if sigma is set to != 0.\n",
" df = df.assign(probabilities_T=df.probabilities_Y)\n",
"\n",
" df.probabilities_T += npr.normal(size=nJudges_M * nSubjects_N) * sigma\n",
"\n",
" # Sort by judges then probabilities in decreasing order\n",
" # I.e. the most dangerous for each judge are first.\n",
" df.sort_values(by=[\"judgeID_J\", \"probabilities_T\"],\n",
"\n",
" # Iterate over the data. Subject is in the top (1-r)*100% if\n",
" # his within-judge-index is over acceptance threshold times\n",
" # the number of subjects assigned to each judge. If subject\n",
" # is over the limit they are assigned a zero, else one.\n",
" df.reset_index(drop=True, inplace=True)\n",
"\n",
" df['decision_T'] = np.where((df.index.values % nSubjects_N) <\n",
" ((1 - df['acceptanceRate_R']) * nSubjects_N),\n",
" 0, 1)\n",
"\n",
" # Halve the data set to test and train\n",
" train, test = train_test_split(df, test_size=0.5)\n",
" train_labeled = train.copy()\n",
" test_labeled = test.copy()\n",
" # Set results as NA if decision is negative.\n",
" train_labeled.loc[train_labeled.decision_T == 0, 'result_Y'] = np.nan\n",
" test_labeled.loc[test_labeled.decision_T == 0, 'result_Y'] = np.nan\n",
" return train_labeled, train, test_labeled, test, df"
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Algorithms\n",
"\n",
"### Contraction algorithm\n",
"\n",
"Below is an implementation of Lakkaraju's team's algorithm presented in [their paper](https://helka.finna.fi/PrimoRecord/pci.acm3098066). Relevant parameters to be passed to the function are presented in the description."
]
},
{
"cell_type": "code",
"def contraction(df, judgeIDJ_col, decisionT_col, resultY_col, modelProbS_col,\n",
" accRateR_col, r):\n",
" '''\n",
" This is an implementation of the algorithm presented by Lakkaraju\n",
" et al. in their paper \"The Selective Labels Problem: Evaluating \n",
" Algorithmic Predictions in the Presence of Unobservables\" (2017).\n",
" Parameters:\n",
" df = The (Pandas) data frame containing the data, judge decisions,\n",
" judge IDs, results and probability scores.\n",
" judgeIDJ_col = String, the name of the column containing the judges' IDs\n",
" in df.\n",
" decisionT_col = String, the name of the column containing the judges' decisions\n",
" resultY_col = String, the name of the column containing the realization\n",
" modelProbS_col = String, the name of the column containing the probability\n",
" scores from the black-box model B.\n",
" accRateR_col = String, the name of the column containing the judges' \n",
" acceptance rates\n",
" r = Float between 0 and 1, the given acceptance rate.\n",
" Returns the estimated failure rate at acceptance rate r.\n",
" most_lenient_ID_q = df[judgeIDJ_col].loc[df[accRateR_col].idxmax()]\n",
" # Subset. \"D_q is the set of all observations judged by q.\"\n",
" D_q = df[df[judgeIDJ_col] == most_lenient_ID_q].copy()\n",
" # All observations of R_q have observed outcome labels.\n",
" # \"R_q is the set of observations in D_q with observed outcome labels.\"\n",
" R_q = D_q[D_q[decisionT_col] == 1].copy()\n",
" # Sort observations in R_q in descending order of confidence scores S and\n",
" # assign to R_sort_q.\n",
" # \"Observations deemed as high risk by B are at the top of this list\"\n",
" R_sort_q = R_q.sort_values(by=modelProbS_col, ascending=False)\n",
"\n",
" number_to_remove = int(\n",
" round((1.0 - r) * D_q.shape[0] - (D_q.shape[0] - R_q.shape[0])))\n",
"\n",
" # \"R_B is the list of observations assigned to t = 1 by B\"\n",
" R_B = R_sort_q[number_to_remove:R_sort_q.shape[0]]\n",
"\n",
" return np.sum(R_B[resultY_col] == 0) / D_q.shape[0]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"Generalized performance:\n",
"\n",
"$$\n",
"\\mathbf{gp} = \\sum_{x\\in\\mathcal{X}} f(x)\\delta(F(x) < r)P(X=x)\n",
"$$\n",
"\n",
"and empirical performance:\n",
"\n",
"$$\n",
"\\mathbf{ep} = \\dfrac{1}{n} \\sum_{(x, y) \\in \\mathcal{D}} \\delta(y=0) \\delta(F(x) < r)\n",
"$$\n",
"\n",
"where\n",
"\n",
"$$\n",
"$$\n",
"\n",
"and\n",
"\n",
"$$\n",
"F(x_0) = \\int P(x)~\\delta(P(Y=0|T=1, X=x) > P(Y=0|T=1, X=x_0)) ~ dx = \\int P(x)~\\delta(f(x) > f(x_0)) ~ dx.\n",
"NB: in code the direction of inequality was changed. CDF changed to `bailIndicator` function.\n",
"\n",
"**Rationale for `bailIndicator`:**\n",
"\n",
"* Bail decision is based on prediction $P(Y=0|T=1, X=x)$.\n",
" * Uniform over all judges\n",
"* Judges rationing: \"If this defendant is in the top 10% of \"dangerousness rank\" and my $r = 0.85$, I will jail him.\"\n",
"* Overall: this kind of defendant $(X=x)$ is usually in the $z^{th}$ percentile in dangerousness (sd +- $u$ percentiles). Now, what is the probability that this defendant has $z \\leq 1-r$?\n",
"\n",
"\n",
"\n",
"1. Train model for $P(Y=0|T=1, X=x)$\n",
"* Estimate quantile function for $P(T=1|R=r, X=x)$\n",
"* Calculate $P(Y=0|do(r'), do(x'))=P(Y=0|T=1, X=x') \\cdot P(T=1|R=r', X=x')$ for all instances of the training data\n",
"* Order in ascending order based on the probabilities obtained from previous step\n",
"* Calculate $$\\dfrac{\\sum_{i=0}^{r\\cdot |\\mathcal{D}_{all}|}}{|\\mathcal{D}_{all}|}$$--->"
"def getProbabilityForClass(x, model, class_value):\n",
" Function (wrapper) for obtaining the probability of a class given x and a \n",
" predictive model.\n",
" x = individual features, an array, shape (observations, features)\n",
" model = a trained sklearn model. Predicts probabilities for given x. Should\n",
" accept input of size (observations, features)\n",
" class_value = the resulting class to predict (usually 0 or 1).\n",
" The probabilities of given class label for each x.\n",
" '''\n",
" if x.ndim == 1:\n",
" # if x is vector, transform to column matrix.\n",
" f_values = model.predict_proba(np.array(x).reshape(-1, 1))\n",
" else:\n",
" f_values = model.predict_proba(x)\n",
"\n",
" # Get correct column of predicted class, remove extra dimensions and return.\n",
" return f_values[:, model.classes_ == class_value].flatten()\n",
"\n",
"# def cdf(x_0, model, class_value):\n",
"# '''\n",
"# Cumulative distribution function as described above.\n",
"# '''\n",
"# prediction = lambda x: getProbabilityForClass(\n",
"# np.array([x]).reshape(-1, 1), model, class_value)\n",
"# results = np.zeros(x_0.shape[0])\n",
"\n",
"# for i in range(x_0.shape[0]):\n",
"\n",
"# y_copy = y_values.copy()\n",
"\n",
"# y_copy[x_preds > prediction_x_0[i]] = 0\n",
"\n",
"# results[i] = si.simps(y_copy, x=x_values)\n",
"\n",
"# return results\n",
"\n",
"\n",
"def bailIndicator(r, y_model, x_train, x_test):\n",
" '''\n",
" Indicator function for whether a judge will bail or jail a suspect.\n",
" (1) Calculate recidivism probabilities from training set with a trained \n",
" model and assign them to predictions_train.\n",
" (2) Calculate recidivism probabilities from test set with the trained \n",
" model and assign them to predictions_test.\n",
" (3) Construct a quantile function of the probabilities in\n",
" if pred belongs to a percentile (computed from step (3)) lower than r\n",
" return True\n",
" else\n",
" return False\n",
" Returns:\n",
" --------\n",
" (1) Boolean list indicating a bail decision (bail = True).\n",
" '''\n",
" predictions_train = getProbabilityForClass(x_train, y_model, 0)\n",
" predictions_test = getProbabilityForClass(x_test, y_model, 0)\n",
"\n",
" return [\n",
" scs.percentileofscore(predictions_train, pred, kind='weak') < r\n",
" for pred in predictions_test\n",
"\n",
"def estimatePercentiles(x_train, y_model, N_bootstraps=2000, N_sample=100):\n",
"\n",
" res = np.zeros((N_bootstraps, 101))\n",
"\n",
" percs = np.arange(101)\n",
"\n",
" for i in range(N_bootstraps):\n",
"\n",
" sample = npr.choice(x_train, size=N_sample)\n",
"\n",
" predictions_sample = getProbabilityForClass(sample, y_model, 0)\n",
"\n",
" res[i, :] = np.percentile(predictions_sample, percs)\n",
" return res\n",
"\n",
"\n",
"def calcReleaseProbabilities(r,\n",
" x_train,\n",
" x_test,\n",
" y_model,\n",
" N_bootstraps=2000,\n",
" N_sample=100,\n",
" percentileMatrix=None):\n",
" '''\n",
" Similar to bailIndicator, but calculates probabilities for bail decisions by\n",
" bootstrapping the data set.\n",
" Returns probabilities for positive bail decisions.\n",
" '''\n",
"\n",
" if percentileMatrix is None:\n",
" percentileMatrix = estimatePercentiles(x_train, y_model, N_bootstraps,\n",
" N_sample)\n",
"\n",
" probs = np.zeros(len(x_test))\n",
" for i in range(len(x_test)):\n",
" probs[i] = np.nan\n",
"\n",
" pred = getProbabilityForClass(x_test[i], y_model, 0)\n",
"\n",
" probs[i] = np.mean(pred < percentileMatrix[:, r])\n",
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Performance comparison\n",
"\n",
"Below we try to replicate the results obtained by Lakkaraju and compare their model's performance to the one of ours."
"def fitLogisticRegression(x_train, y_train, x_test, class_value):\n",
" '''\n",
" Fit logistic regression model with given inputs. Checks their shape if \n",
" incompatible.\n",
" \n",
" Parameters:\n",
" \n",
" \n",
" Returns:\n",
" (1) Trained LogisticRegression model\n",
" (2) Probabilities for given test inputs for given class.\n",
" '''\n",
" # Instantiate the model (using the default parameters)\n",
" logreg = LogisticRegression(solver='lbfgs')\n",
"\n",
" # Check shape and fit the model.\n",
" if x_train.ndim == 1:\n",
" logreg = logreg.fit(x_train.values.reshape(-1, 1), y_train)\n",
" else:\n",
" logreg = logreg.fit(x_train, y_train)\n",
" label_probs_logreg = getProbabilityForClass(x_test, logreg, class_value)\n",
" \n",
" return logreg, label_probs_logreg"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### With unobservables in the data\n",
"Lakkaraju says that they used logistic regression. We train the predictive models using only *observed observations*, i.e. observations for which labels are available. We then predict the probability of negative outcome for all observations in the test data and attach it to our data set."
{
"name": "stdout",
"output_type": "stream",
"text": [
"name": "stderr",
"output_type": "stream",
"text": [
"/Users/rikulain/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:107: IntegrationWarning: The maximum number of subdivisions (50) has been achieved.\n",
" If increasing the limit yields no improvement it is advised to analyze \n",
" the integrand in order to determine the difficulties. If the position of a \n",
" local difficulty can be determined (singularity, discontinuity) one will \n",
" probably gain from splitting up the interval and calling the integrator \n",
" on the subranges. Perhaps a special-purpose integrator should be used.\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"/Users/rikulain/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:107: IntegrationWarning: The maximum number of subdivisions (50) has been achieved.\n",
" If increasing the limit yields no improvement it is advised to analyze \n",
" the integrand in order to determine the difficulties. If the position of a \n",
" local difficulty can be determined (singularity, discontinuity) one will \n",
" probably gain from splitting up the interval and calling the integrator \n",
" on the subranges. Perhaps a special-purpose integrator should be used.\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"name": "stderr",
"output_type": "stream",
"text": [
"/Users/rikulain/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:107: IntegrationWarning: The maximum number of subdivisions (50) has been achieved.\n",
" If increasing the limit yields no improvement it is advised to analyze \n",
" the integrand in order to determine the difficulties. If the position of a \n",
" local difficulty can be determined (singularity, discontinuity) one will \n",
" probably gain from splitting up the interval and calling the integrator \n",
" on the subranges. Perhaps a special-purpose integrator should be used.\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"/Users/rikulain/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:107: IntegrationWarning: The maximum number of subdivisions (50) has been achieved.\n",
" If increasing the limit yields no improvement it is advised to analyze \n",
" the integrand in order to determine the difficulties. If the position of a \n",
" local difficulty can be determined (singularity, discontinuity) one will \n",
" probably gain from splitting up the interval and calling the integrator \n",
" on the subranges. Perhaps a special-purpose integrator should be used.\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"4 "
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"/Users/rikulain/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:107: IntegrationWarning: The maximum number of subdivisions (50) has been achieved.\n",
" If increasing the limit yields no improvement it is advised to analyze \n",
" the integrand in order to determine the difficulties. If the position of a \n",
" local difficulty can be determined (singularity, discontinuity) one will \n",
" probably gain from splitting up the interval and calling the integrator \n",
" on the subranges. Perhaps a special-purpose integrator should be used.\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"[2] 0 "
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"/Users/rikulain/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:107: IntegrationWarning: The maximum number of subdivisions (50) has been achieved.\n",
" If increasing the limit yields no improvement it is advised to analyze \n",
" the integrand in order to determine the difficulties. If the position of a \n",
" local difficulty can be determined (singularity, discontinuity) one will \n",
" probably gain from splitting up the interval and calling the integrator \n",
" on the subranges. Perhaps a special-purpose integrator should be used.\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"1 "
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"/Users/rikulain/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:107: IntegrationWarning: The maximum number of subdivisions (50) has been achieved.\n",
" If increasing the limit yields no improvement it is advised to analyze \n",
" the integrand in order to determine the difficulties. If the position of a \n",
" local difficulty can be determined (singularity, discontinuity) one will \n",
" probably gain from splitting up the interval and calling the integrator \n",
" on the subranges. Perhaps a special-purpose integrator should be used.\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"2 "
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"/Users/rikulain/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:107: IntegrationWarning: The maximum number of subdivisions (50) has been achieved.\n",
" If increasing the limit yields no improvement it is advised to analyze \n",
" the integrand in order to determine the difficulties. If the position of a \n",
" local difficulty can be determined (singularity, discontinuity) one will \n",
" probably gain from splitting up the interval and calling the integrator \n",
" on the subranges. Perhaps a special-purpose integrator should be used.\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"3 "
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"/Users/rikulain/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:107: IntegrationWarning: The maximum number of subdivisions (50) has been achieved.\n",
" If increasing the limit yields no improvement it is advised to analyze \n",
" the integrand in order to determine the difficulties. If the position of a \n",
" local difficulty can be determined (singularity, discontinuity) one will \n",
" probably gain from splitting up the interval and calling the integrator \n",
" on the subranges. Perhaps a special-purpose integrator should be used.\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"4 "
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"/Users/rikulain/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:107: IntegrationWarning: The maximum number of subdivisions (50) has been achieved.\n",
" If increasing the limit yields no improvement it is advised to analyze \n",
" the integrand in order to determine the difficulties. If the position of a \n",
" local difficulty can be determined (singularity, discontinuity) one will \n",
" probably gain from splitting up the interval and calling the integrator \n",
" on the subranges. Perhaps a special-purpose integrator should be used.\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"[3] 0 "
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"/Users/rikulain/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:107: IntegrationWarning: The maximum number of subdivisions (50) has been achieved.\n",
" If increasing the limit yields no improvement it is advised to analyze \n",
" the integrand in order to determine the difficulties. If the position of a \n",
" local difficulty can be determined (singularity, discontinuity) one will \n",
" probably gain from splitting up the interval and calling the integrator \n",
" on the subranges. Perhaps a special-purpose integrator should be used.\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"1 "
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"/Users/rikulain/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:107: IntegrationWarning: The maximum number of subdivisions (50) has been achieved.\n",
" If increasing the limit yields no improvement it is advised to analyze \n",
" the integrand in order to determine the difficulties. If the position of a \n",
" local difficulty can be determined (singularity, discontinuity) one will \n",
" probably gain from splitting up the interval and calling the integrator \n",
" on the subranges. Perhaps a special-purpose integrator should be used.\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"2 "
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"/Users/rikulain/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:107: IntegrationWarning: The maximum number of subdivisions (50) has been achieved.\n",
" If increasing the limit yields no improvement it is advised to analyze \n",
" the integrand in order to determine the difficulties. If the position of a \n",
" local difficulty can be determined (singularity, discontinuity) one will \n",
" probably gain from splitting up the interval and calling the integrator \n",
" on the subranges. Perhaps a special-purpose integrator should be used.\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"3 "
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"/Users/rikulain/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:107: IntegrationWarning: The maximum number of subdivisions (50) has been achieved.\n",
" If increasing the limit yields no improvement it is advised to analyze \n",
" the integrand in order to determine the difficulties. If the position of a \n",
" local difficulty can be determined (singularity, discontinuity) one will \n",
" probably gain from splitting up the interval and calling the integrator \n",
" on the subranges. Perhaps a special-purpose integrator should be used.\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"4 "
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"/Users/rikulain/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:107: IntegrationWarning: The maximum number of subdivisions (50) has been achieved.\n",
" If increasing the limit yields no improvement it is advised to analyze \n",
" the integrand in order to determine the difficulties. If the position of a \n",
" local difficulty can be determined (singularity, discontinuity) one will \n",
" probably gain from splitting up the interval and calling the integrator \n",
" on the subranges. Perhaps a special-purpose integrator should be used.\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"[4] 0 "
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"/Users/rikulain/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:107: IntegrationWarning: The maximum number of subdivisions (50) has been achieved.\n",
" If increasing the limit yields no improvement it is advised to analyze \n",
" the integrand in order to determine the difficulties. If the position of a \n",
" local difficulty can be determined (singularity, discontinuity) one will \n",
" probably gain from splitting up the interval and calling the integrator \n",
" on the subranges. Perhaps a special-purpose integrator should be used.\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"1 "
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"/Users/rikulain/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:107: IntegrationWarning: The maximum number of subdivisions (50) has been achieved.\n",
" If increasing the limit yields no improvement it is advised to analyze \n",
" the integrand in order to determine the difficulties. If the position of a \n",
" local difficulty can be determined (singularity, discontinuity) one will \n",
" probably gain from splitting up the interval and calling the integrator \n",
" on the subranges. Perhaps a special-purpose integrator should be used.\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"2 "
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"/Users/rikulain/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:107: IntegrationWarning: The maximum number of subdivisions (50) has been achieved.\n",
" If increasing the limit yields no improvement it is advised to analyze \n",
" the integrand in order to determine the difficulties. If the position of a \n",
" local difficulty can be determined (singularity, discontinuity) one will \n",
" probably gain from splitting up the interval and calling the integrator \n",
" on the subranges. Perhaps a special-purpose integrator should be used.\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"3 "
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"/Users/rikulain/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:107: IntegrationWarning: The maximum number of subdivisions (50) has been achieved.\n",
" If increasing the limit yields no improvement it is advised to analyze \n",
" the integrand in order to determine the difficulties. If the position of a \n",
" local difficulty can be determined (singularity, discontinuity) one will \n",
" probably gain from splitting up the interval and calling the integrator \n",
" on the subranges. Perhaps a special-purpose integrator should be used.\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"4 "
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"/Users/rikulain/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:107: IntegrationWarning: The maximum number of subdivisions (50) has been achieved.\n",
" If increasing the limit yields no improvement it is advised to analyze \n",
" the integrand in order to determine the difficulties. If the position of a \n",
" local difficulty can be determined (singularity, discontinuity) one will \n",
" probably gain from splitting up the interval and calling the integrator \n",
" on the subranges. Perhaps a special-purpose integrator should be used.\n"