From df13a8138d55f668424b66472a609731dcbfeda9 Mon Sep 17 00:00:00 2001 From: Taro Yabuki Date: Mon, 1 Apr 2024 14:48:18 +0900 Subject: [PATCH] Update README.md --- README.md | 31 +- code/python/python.Rmd | 2206 ---------------------------------------- 2 files changed, 17 insertions(+), 2220 deletions(-) delete mode 100644 code/python/python.Rmd diff --git a/README.md b/README.md index dd426cc..d6938a4 100644 --- a/README.md +++ b/README.md @@ -54,33 +54,36 @@ Wolfram|Alpha，Python，R，Mathematicaをフル活用して，大学教養レ 19. 固有値と固有ベクトル 20. 特異値分解と擬似逆行列 -## 実行環境とサンプルコード +## 実行環境とコード 目次を表示させると見やすくなります[^1]． ### 推奨環境 -システム|サンプルコード|実行結果 +システム|コード|実行結果 --|--|-- [Wolfram\|Alpha](https://www.wolframalpha.com/)|[リンク集](code/wolframalpha)|[キャプチャ画像](code/wolframalpha/results/README.md) -Python notebook|[![Binder](https://mybinder.org/badge_logo.svg)](https://mybinder.org/v2/gh/taroyabuki/comath/HEAD?labpath=code%2Fpython%2Fpython-binder.ipynb) [![Open in Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/taroyabuki/comath/blob/master/code/python/python.ipynb)|[![Open in Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/taroyabuki/comath/blob/master/code/python/python-results.ipynb) -R notebook|[![Binder](https://mybinder.org/badge_logo.svg)](https://mybinder.org/v2/gh/taroyabuki/comath/HEAD?labpath=code%2Fr%2Fr-binder.ipynb) [![Open in Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/taroyabuki/comath/blob/master/code/r/r.ipynb)|[![Open in Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/taroyabuki/comath/blob/master/code/r/r-results.ipynb) -[Wolfram Cloud](https://www.wolframcloud.com)|動的生成[^2]| -Mathematica（[Pi版は無料](https://www.wolfram.com/raspberry-pi/)）|[mathematica.nb](code/mathematica/mathematica.nb)[^3]|[PDF](code/mathematica/results/mathematica-results.pdf) +Python (JupyterLab)|[![Binder](https://mybinder.org/badge_logo.svg)](https://mybinder.org/v2/gh/taroyabuki/comath/HEAD?labpath=code%2Fpython%2Fpython-binder.ipynb) +R (JupyterLab)|[![Binder](https://mybinder.org/badge_logo.svg)](https://mybinder.org/v2/gh/taroyabuki/comath/HEAD?labpath=code%2Fr%2Fr-binder.ipynb) +R (RStudio)|[![Binder](https://mybinder.org/badge_logo.svg)](https://mybinder.org/v2/gh/taroyabuki/comath/HEAD?urlpath=rstudio)[^2]|[Knit to HTML](https://taroyabuki.github.io/comath/r.html) +[Wolfram Cloud](https://www.wolframcloud.com)|動的生成[^3]| +Mathematica（[Pi版は無料](https://www.wolfram.com/raspberry-pi/)）|[mathematica.nb](code/mathematica/mathematica.nb)[^4]|[PDF](code/mathematica/results/mathematica-results.pdf) ### そのほかの環境 -システム|サンプルコード ---|-- -JupyterLab[^4] (Wolfram言語)|[wolfram.ipynb](code/mathematica/wolfram.ipynb) -RStudio (R)|[r.Rmd](code/r/r.Rmd) +システム|コード|実行結果 +--|--|-- +Python (Google Colab)|[![Open in Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/taroyabuki/comath/blob/master/code/python/python.ipynb)|[![Open in Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/taroyabuki/comath/blob/master/code/python/python-results.ipynb) +R (Google Colab)|[![Open in Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/taroyabuki/comath/blob/master/code/r/r.ipynb)|[![Open in Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/taroyabuki/comath/blob/master/code/r/r-results.ipynb) +Wolfram言語 (JupyterLab[^5])|[wolfram.ipynb](code/mathematica/wolfram.ipynb) ### そのほかのコード - [Pythonで使用するモジュールの一括読込（1.2.2項）](code/imports.py) -[^1]: JupyterLabとGoogle Colabでは「≡」のようなアイコンをクリック，RStudioでは「Outline」をクリックします．Wolfram CloudとMathematicaでは，Ctrl-`A`（⌘-`A`）で全選択→View→Close All Subgroupsで全てのセルを閉じます． -[^2]: 無料アカウントではノートブックをアップロードできないので，新規作成したノートブックで，`Import["https://raw.githubusercontent.com/taroyabuki/comath/main/code/mathematica/createnb.m", "NB"]`を実行し，コードを読み込みます． -[^3]: ダウンロードボタンをクリックしてダウンロードするか，Mathematicaで新規作成したノートブックで，`NotebookWrite[EvaluationNotebook[], Import["https://raw.githubusercontent.com/taroyabuki/comath/main/code/mathematica/mathematica.nb", "NB"]]`を実行して読み込みます． -[^4]: Wolfram言語用のJupyterLabは，[Wolfram Engine Community Edition](https://www.wolfram.com/engine/)と[Wolfram Language kernel for Jupyter notebooks](https://github.com/WolframResearch/WolframLanguageForJupyter)を組み合わせて実現します（無料）． +[^1]: JupyterLabとGoogle Colabでは「≡」のようなアイコン，RStudioでは「Outline」をクリックします．Wolfram CloudとMathematicaでは，Ctrl-`A`（⌘-`A`）で全選択→View→Close All Subgroupsで全てのセルを閉じます． +[^2]: code/r/r.Rmdを開いてください． +[^3]: 無料アカウントではノートブックをアップロードできないので，新規作成したノートブックで，`Import["https://raw.githubusercontent.com/taroyabuki/comath/main/code/mathematica/createnb.m", "NB"]`を実行し，コードを読み込みます． +[^4]: ダウンロードボタンをクリックしてダウンロードするか，Mathematicaで新規作成したノートブックで，`NotebookWrite[EvaluationNotebook[], Import["https://raw.githubusercontent.com/taroyabuki/comath/main/code/mathematica/mathematica.nb", "NB"]]`を実行して読み込みます． +[^5]: Wolfram言語用のJupyterLabは，[Wolfram Engine Community Edition](https://www.wolfram.com/engine/)と[Wolfram Language kernel for Jupyter notebooks](https://github.com/WolframResearch/WolframLanguageForJupyter)を組み合わせて実現します（無料）． diff --git a/code/python/python.Rmd b/code/python/python.Rmd deleted file mode 100644 index 0d28f0c..0000000 --- a/code/python/python.Rmd +++ /dev/null @@ -1,2206 +0,0 @@ ---- -title: "Python　コンピュータでとく数学" -output: - html_document: - toc: true - toc_depth: 1 - toc_float: true ---- - -[矢吹太朗『コンピュータでとく数学』（オーム社,\ 2024）](https://github.com/taroyabuki/comath) - -```{r setup, include=FALSE} -system("pip install matplotlib IPython pandas seaborn see scikit-image scikit-learn statsmodels sympy") -``` - -```{python} -from IPython.display import display -``` - -# 1 実行環境 - -```{python} -#!python -m pip install see -``` - -```{python} -import matplotlib.pyplot as plt -import numpy as np -import pandas as pd -import seaborn as sns -import statsmodels.api as sm -import statsmodels.formula.api as smf -import sympy as sym -from collections import Counter -from patsy import dmatrices -from scipy import linalg, stats -from scipy.integrate import quad -from scipy.optimize import minimize -from sklearn.linear_model import LinearRegression -from statsmodels.stats.proportion import binom_test -from statsmodels.stats.weightstats import CompareMeans, DescrStatsW, ttest_ind -from sympy import * -from sympy.stats import * -from sympy.plotting import plot3d -``` - -```{python} -np.sqrt(2) -``` - -```{python} -data = pd.DataFrame({'x1': [1, 3, 6, 10], 'y': [7, 1, 6, 14]}) -model = smf.ols('y ~ x1', data).fit() -model.params -``` - -```{python} -from see import see -see(model) -# .outlier_test() .params .predict() .pvalues -# .remove_data() .resid .resid_pearson .rsquared -# .rsquared_adj .save() .scale .ssr -# .summary() .summary2() .t_test() -``` - -# 2 数と変数 - -```{python} -2 * (-3) -``` - -```{python} -(1 + 2) * 3 -``` - -```{python} -2**10 -``` - -```{python} --2 < -1 -``` - -```{python} -2 + 2 == 5 -``` - -```{python} -10 if (7 < 5) else 20 -``` - -```{python} -var('x') # xが変数であることの宣言 -x < 1 -``` - -```{python} -var('x y') -Eq(x, y) -``` - -```{python} -simplify(Eq((x**2 - 1), (x + 1) * (x - 1))) -``` - -```{python} -not(0<1) -``` - -```{python} -(0<1)or(2>3) -``` - -```{python} -(0 < 1) and (2 > 3) -``` - -```{python} -var('x') -Not(10 < x) -``` - -```{python} -x = 5; x == 5 -``` - -```{python} -a = 1 + 2 -b = 9 -a * (b + 1) -``` - -```{python} -a = 1 + 2; b = 9; a * (b + 1) -``` - -```{python} -a = 3 -var('a') # 変数を記号にする． -expand((a + 1)**2) -``` - -```{python} -import keyword -keyword.kwlist -``` - -```{python} -x1 = 2; x2 = 3; x1 + x2 -``` - -```{python} -x = 1; y = x + 1; x = 2; y -``` - -```{python} -var('x') -f = 2 * x + 3 -f.subs(x, 5) -``` - -```{python} -var('a b x y') -g = a + b -g.subs(((a, x), (b, y))) -``` - -```{python} -f = lambda x: 2 * x + 3 -f(5) -``` - -```{python} -def f(x): return 2 * x + 3 -f(5) -``` - -```{python} -f = lambda x: 2 * x + 3; var('a') -g = f(a) -f(5), g.subs(a, 5) -``` - -```{python} -f = lambda x: 1 / x -f(1) -``` - -```{python} -(lambda x: 2 * x + 3)(5) -``` - -```{python} -f = lambda x, y: x + y -f(2, 3) -``` - -```{python} -g = lambda x: x[0] + x[1] -x = (2, 3); g(x) -``` - -```{python} -f(*x) -``` - -```{python} -g((2, 3)) -``` - -```{python} -var('x') -expand((x + 1)**2) -``` - -```{python} -N(sqrt(2), 30) -``` - -```{python} -import struct -tmp = bin(int(struct.pack('>d', 0.1).hex(), 16))[2:].zfill(64) # ビット列 -s, e, f = tmp[0], tmp[1:12], tmp[12:64] # 1, 11, 52桁に分ける． -s, e, f -``` - -```{python} -s, e, f = int(s), sym.S(int(e, 2)), sym.S(int(f, 2)) # 数値に変換する． -(-1)**s * (1 + f / 2**52) * 2**(e - 1023) -``` - -```{python} -import sympy -sympy.sqrt(5)**2 == 5 # True -``` - -```{python} -import math -import numpy -numpy.sqrt(5)**2 == 5, math.sqrt(5)**2 == 5 # いずれもFalse -``` - -```{python} -0.1 + 0.2 == 0.3 -``` - -```{python} -np.isclose(0.1 + 0.2, 0.3) -``` - -```{python} -sym.S(1) / 10 + sym.S(2) / 10 == sym.S(3) / 10 -``` - -```{python} -var('x') -simplify(sin(x)**2 + cos(x)**2) -``` - -```{python} -sqrtdenest(sqrt(5 + 2 * sqrt(6))) -``` - -```{python} -refine(sqrt((x - 1)**2), Q.nonnegative(x - 1)) -``` - -# 3 データ構造 - -```{python} -v = [2, 3, 5]; len(v) -``` - -```{python} -v[2] = 0.5; v -``` - -```{python} -np.arange(5) -``` - -```{python} -np.arange(0, 1.01, 0.1) # 終点（ここでは1）より大きい値を指定する． -``` - -```{python} -np.linspace(0, 100, 5) -``` - -```{python} -v = np.array([2, 3]) -1.1 * v -``` - -```{python} -u = np.array([10, 20]); v = np.array([2, 3]) -u + v -``` - -```{python} -v + 1 -``` - -```{python} -u = np.array([10, 20]); v = np.array([2, 3]) -u.dot(v) -``` - -```{python} -a = [2, 3, 4]; b = a.copy(); b[2] = 0.5; a -``` - -```{python} -v = [2, -1, 3, -2] -[x for x in v if x > 0] # 内包表記 -``` - -```{python} -v = [2, -1, 3, -2] -np.heaviside(v, 0) -``` - -```{python} -v = [2, -1, 3, -2] -n = len(v) # vのサイズ -u = [None] * n # Noneは「値がない」ということ． -for i in range(n): u[i] = 0 if v[i] < 0 else 1 -u -``` - -```{python} -v = [2, -1, 3, -2] -[(0 if x < 0 else 1) for x in v] -``` - -```{python} -v = [2, -1, 3, -2] -f = lambda x: 0 if x < 0 else 1 -[f(x) for x in v] -``` - -```{python} -f = np.vectorize(lambda x: 0 if x < 0 else 1) -f(v) -``` - -```{python} -u = [1, 7, 2, 9]; v = [2, 3, 5, 7] -[(-1 if a < b else 1) for a, b in zip(u, v)] -``` - -```{python} -x = {"apple" : "りんご", "orange" : "みかん"} -x["orange"] -``` - -```{python} -x["grape"] = "ぶどう" -x["grape"] -``` - -```{python} -x.pop("apple") -"apple" in x -``` - -```{python} -df = pd.DataFrame({'name': ['A', 'B', 'C'], - 'english': [60, 90, 70], - 'math': [70, 80, 90], - 'gender': ['f', 'm', 'm']}) -df -``` - -```{python} -df = pd.DataFrame(({'name': 'A', 'english': 60, 'math': 70, 'gender': 'f'}, - {'name': 'B', 'english': 90, 'math': 80, 'gender': 'm'}, - {'name': 'C', 'english': 70, 'math': 90, 'gender': 'm'})) -df -``` - -```{python} -df[['english', 'math']] -``` - -```{python} -df["english"] -``` - -```{python} -df.english -``` - -```{python} -m = df.iloc[:, [1, 2]].values; m -``` - -```{python} -english, math = m.T; english, math -``` - -# 4 可視化と方程式 - -```{python} -var('x') -plot(x**2 + 2 * x - 4, (x, -5, 3)); -``` - -```{python} -x = np.linspace(-5, 3, 101) -y = x**2 + 2 * x - 4 -plt.plot(x, y); -``` - -```{python} -var('x y') -plot3d(x**2 + y**2, (x, -1, 1), (y, -1, 1)); -``` - -```{python} -plotting.plot_contour(x**2 + y**2, (x, -1, 1), (y, -1, 1)); -``` - -```{python} -plot_implicit(Eq(x**2 + y**2, 1)); -``` - -```{python} -plot_implicit(x**2 + y**2 <= 1); -``` - -```{python} -plot_implicit(Or(Eq(2 * x + 3 * y, 8), Eq(5 * x - 7 * y, -9))); -``` - -```{python} -plot_implicit(And(y <= x, y >= x**2), (x, 0, 2), (y, 0, 2)); -``` - -```{python} -var('x') -a, b = sorted(solve(Eq(x, x**2), x)) -integrate(x - x**2, (x, a, b)) -``` - -```{python} -var('x') -solve(Eq(x**2 + 2 * x - 4, 0), x) -``` - -```{python} -a, b = solve(x**2 + 2 * x - 4, x) -a + b -``` - -```{python} -n = 3; simplify(sum(solve(x**n + 2 * x - 4, x))) -``` - -```{python} -var('x y') -sol = solve([Eq(2 * x + 3 * y, 8), Eq(5 * x - 7 * y, -9)], [x, y]); sol -``` - -```{python} -x1, y1 = sol.values() -x1 + y1 -``` - -```{python} -var('x'); f = 2**x + sin(x) -nsolve(f, x, 0) -``` - -```{python} -var('x') -solve(Lt(x**2 + 2 * x - 4, 0), x) -``` - -# 5 論理式 - -# 6 1次元のデータ - -```{python} -a = pd.Series([36, 43, 53, 55, 56, 56, 57, 60, 61, 73]) -b = pd.Series([34, 39, 39, 49, 50, 52, 52, 55, 83, 97]) -a.hist(); -``` - -```{python} -a.value_counts(bins=np.arange(20, 81, 20), sort=False) -``` - -```{python} -x = [7, 3, 1, 3, 4, 7, 7, 7, 10, 3] -f = Counter(x); f -``` - -```{python} -pd.DataFrame({'A':a, 'B':b}).boxplot(); -``` - -```{python} -a = pd.Series([36, 43, 53, 55, 56, 56, 57, 60, 61, 73]) -a.mean() -``` - -```{python} -b = pd.Series([34, 39, 39, 49, 50, 52, 52, 55, 83, 97]) -sum(b) / len(b), b.sum() / b.count() # 二つの方法 -``` - -```{python} -(a - a.mean()).mean() -``` - -```{python} -a.var(ddof=1) -``` - -```{python} -sum((b - b.mean())**2) / (len(b) - 1) -``` - -```{python} -z = stats.zscore(a, ddof=1); z -``` - -```{python} -z.mean(), z.std(ddof=1) -``` - -```{python} -(a - a.mean()) / a.std(ddof=1) -``` - -```{python} -a.std(ddof=1) * z + a.mean() -``` - -```{python} -10 * z + 50 -``` - -# 7 2次元のデータ - -```{python} -x = pd.Series([35, 45, 55, 65, 75]) -y = pd.Series([114, 124, 143, 158, 166]) -plt.scatter(x, y); -``` - -```{python} -x = pd.Series([35, 45, 55, 65, 75]) -y = pd.Series([114, 124, 143, 158, 166]) -x.cov(y, ddof=1), np.cov(x, y, ddof=1)[0, 1] # 二つの方法 -``` - -```{python} -np.cov(x, y, ddof=1) -``` - -```{python} -(x - x.mean()) @ (y - y.mean()) / (len(x) - 1) -``` - -```{python} -x.corr(y), np.corrcoef(x, y)[0, 1] # 二つの方法 -``` - -```{python} -x = pd.Series([35, 45, 55, 65, 75]) -y = pd.Series([114, 124, 143, 158, 166]) -data = pd.DataFrame({'x': x, 'y': y}) -model = smf.ols('y ~ x', data).fit() -model.params -``` - -```{python} -model.predict({'x': 40}) -``` - -```{python} -sns.regplot(x=x, y=y, ci=None) -``` - -```{python} -var('a b'); L = sum((y - (a * x + b))**2); L -``` - -```{python} -vars = var('p q r s t u') -sol = solve(Eq(L, p * (a - q)**2 + r * (b - (s* a + t))**2 + u), vars) -print(sol) -q.subs(sol), (s * q + t).subs(sol) -``` - -```{python} -a = x.cov(y, ddof=1) / x.var(ddof=1) -b = y.mean() - a * x.mean() -a, b -``` - -```{python} -anscombe = sns.load_dataset('anscombe') -data = anscombe[anscombe.dataset == 'I'] -print(data.x.corr(data.y)) -model = smf.ols('y ~ x', data).fit(); print(model.params) -sns.regplot(x=data.x, y=data.y, ci=None); -``` - -# 8 確率変数と確率分布 - -```{python} -X = DiscreteUniform('X', range(1, 7)) # X = Die('X', 6)でもよい． -density(X)(2) -``` - -```{python} -rv = stats.randint(1, 7) -rv.pmf(2) -``` - -```{python} -P(Eq(X, 2)) # P(X == 2)ではない． -``` - -```{python} -data = list(sample_iter(X, numsamples=1000)) -plt.hist(data); # 結果は割愛 -``` - -```{python} -data = rv.rvs(size=1000) -plt.hist(data); # 結果は割愛 -``` - -```{python} -x = range(1, 7); y = [density(X)(x) for x in x] -_, ax = plt.subplots() # 結果のうち，使わない部分を_とする． -ax.hist(data, bins=np.arange(0.5, 7, 1), density=True, alpha=0.3) -ax.scatter(x, y); -``` - -```{python} -x = range(1, 7); y = rv.pmf(x) -_, ax = plt.subplots() # 結果のうち，使わない部分を_とする． -ax.hist(data, bins=np.arange(0.5, 7, 1), density=True, alpha=0.3) -ax.scatter(x, y); -``` - -```{python} -X = Bernoulli('X', p=sym.S(3) / 10) -data = list(sample_iter(X, numsamples=1000)) -np.bincount(data), Counter(data) # 二つの方法 -``` - -```{python} -rv = stats.bernoulli(3 / 10) -data = rv.rvs(1000) -np.bincount(data), Counter(data) # 二つの方法 -``` - -```{python} -X = Binomial('X', 10, sym.S(3) / 10) -density(X)(3) -``` - -```{python} -rv = stats.binom(10, 3 / 10) -rv.pmf(3) -``` - -```{python} -P(Eq(X, 3)) # P(X == 3)ではない． -``` - -```{python} -var('n p x') -X = Binomial('X', n, p) -density(X)(x) -``` - -```{python} -n = 10; p = sym.S(3) / 10; X = Binomial('X', n, p) -data = list(sample_iter(X, numsamples=1000)) -plt.hist(data); # 結果は割愛 -``` - -```{python} -n = 10; p = 3 / 10; rv = stats.binom(n, p) -data = rv.rvs(1000) -plt.hist(data); # 結果は割愛 -``` - -```{python} -x = range(0, n + 1); y = [density(X)(x) for x in x] -_, ax = plt.subplots() -ax.hist(data, bins=np.arange(-0.5, n + 1, 1), density=True, alpha=0.5) -ax.scatter(x, y); -``` - -```{python} -x = range(0, n + 1); y = rv.pmf(x) -_, ax = plt.subplots() -ax.hist(data, bins=np.arange(-0.5, n + 1, 1), density=True, alpha=0.5) -ax.scatter(x, y); -``` - -```{python} -X = Binomial('X', 10, S(3) / 10) -cdf(X)[3] # cdf(X)(x)ではない． -``` - -```{python} -rv = stats.binom(10, 3 / 10) -rv.cdf(3) -``` - -```{python} -P(X <= 3) -``` - -```{python} -sum([density(X)(k) for k in range(4)]) -``` - -```{python} -sum([rv.pmf(k) for k in range(4)]) -``` - -```{python} -X = Uniform('X', 0, 360) -cdf(X)(200), cdf(X)(150), cdf(X)(200) - cdf(X)(150) -``` - -```{python} -rv = stats.uniform(0, 360) -rv.cdf(200), rv.cdf(150), rv.cdf(200) - rv.cdf(150) -``` - -```{python} -# P(150 <= X <= 200)やP(150 <= X and X <= 200)ではない． -P(And(150 <= X, X <= 200)) -``` - -```{python} -var('x'); integrate(density(X)(x), (x, 150, 200)) -``` - -```{python} -quad(rv.pdf, 150, 200) -``` - -```{python} -var('t x'); integrate(density(X)(t), (t, 0, x)) -``` - -```{python} -diff(x / 360, x) -``` - -```{python} -data = list(sample_iter(X, numsamples=1000)) -plt.hist(data); # 結果は割愛 -``` - -```{python} -data = rv.rvs(1000) -plt.hist(data); # 結果は割愛 -``` - -```{python} -x = np.linspace(0, 360, 100); y = [density(X)(x) for x in x] -_, ax = plt.subplots() -ax.hist(data, bins='sturges', density=True, alpha=0.5) -ax.plot(x, y); -``` - -```{python} -data = rv.rvs(1000) -x = np.linspace(0, 360, 100); y = rv.pdf(x) -_, ax = plt.subplots() -ax.hist(data, bins='sturges', density=True, alpha=0.5) -ax.plot(x, y); -``` - -```{python} -X = Normal('X', 6, 2) -N((cdf(X)(6 + 3 * 2) - cdf(X)(6 - 3 * 2))) -``` - -```{python} -rv = stats.norm(6, 2) -rv.cdf(6 + 3 * 2) - rv.cdf(6 - 3 * 2) -``` - -```{python} -#P(6 - 3 * 2 <= X <= 6 + 3 * 2)ではない． -#P(6 - 3 * 2 <= X and X <= 6 + 3 * 2)ではない． -N(P(And(6 - 3 * 2 <= X, X <= 6 + 3 * 2))) -``` - -```{python} -var('x') -N(integrate(density(X)(x), (x, 6 - 3 * 2, 6 + 3 * 2))) -``` - -```{python} -quad(rv.pdf, 6 - 3 * 2, 6 + 3 * 2) -``` - -```{python} -var('mu sigma x'); X = Normal('X', mu, sigma) -a, b = mu - 3 * sigma, mu + 3 * sigma -(N((cdf(X)(b) - cdf(X)(a))), # 方法1 - N(integrate(density(X)(x), (x, a, b)))) # 方法3 -``` - -```{python} -var('mu sigma x') -X = Normal('X', mu, sigma) -density(X)(x) -``` - -```{python} -X1 = Normal('X1', 0, 1); X2 = Normal('X2', 2, 1); var('x') -plot(density(X1)(x), density(X2)(x), (x, -6, 6)); -``` - -```{python} -rv1 = stats.norm(0, 1); rv2 = stats.norm(2, 1); x = np.linspace(-6, 6, 100) -pd.DataFrame({'x': x, 'X1': rv1.pdf(x), 'X2': rv2.pdf(x)}).plot(x='x'); -``` - -```{python} -X3 = Normal('X3', 0, 2); var('x') -plot(density(X1)(x), density(X3)(x), (x, -6, 6)); -``` - -```{python} -rv3 = stats.norm(0, 2); x = np.linspace(-6, 6, 100) -pd.DataFrame({'x': x, 'X1': rv1.pdf(x), 'X3': rv3.pdf(x)}).plot(x='x'); -``` - -```{python} -Xs = [0, 100, 1000, 10000]; tmp = [0.9, 0.08, 0.015, 0.005] -Ps = np.array(tmp) / sum(tmp) # 念のため合計を1にする． -X = FiniteRV('X', dict(zip(Xs, Ps))) # 確率分布の定義 -data = sample_iter(X, numsamples=1000) -Counter(data) -``` - -```{python} -Xs = [0, 100, 1000, 10000]; tmp = [0.9, 0.08, 0.015, 0.005] -Ps = np.array(tmp) / sum(tmp) # 念のため合計を1にする． -rv = stats.rv_discrete(values=(Xs, Ps)) # 確率分布の定義 -data = rv.rvs(size=1000) -Counter(data) -``` - -```{python} -var('x') -f = Lambda(x, Abs(x)); a, b = -1, 1; -s = integrate(f(x), (x, a, b)) # “全確率” -X = ContinuousRV(x, # 確率分布の定義 - f(x) / s, # 念のため“全確率”で割る． - set=Interval(a, b)) -data = list(sample_iter(X, numsamples=1000)) -plt.hist(data); -``` - -```{python} -f = lambda x: abs(x); a, b = -1, 1; -s = quad(f, a, b)[0] # “全確率” -class MyX(stats.rv_continuous): # 確率分布の定義の準備 - def _pdf(self, x): return f(x) / s # 念のため“全確率”で割る． -rv = MyX(a=a, b=b) # 確率分布の定義 -data = rv.rvs(size=1000) -plt.hist(data); -``` - -```{python} -X = DiscreteUniform('X', range(1, 7)) -Y = X**3 % 4 -data = list(sample_iter(Y, numsamples=1000)) -plt.hist(data); -``` - -```{python} -X = Uniform('X', 0, 1) -Y = X**2 -data = list(sample_iter(Y, numsamples=1000)) -plt.hist(data); -``` - -```{python} -simplify(density(Y)) -``` - -```{python} -var('a b mu sigma'); X = Normal('X', mu, sigma); Y = a * X + b -simplify(density(Y)) -``` - -```{python} -Xs = [0, 100, 1000, 10000]; Ps = [0.9, 0.08, 0.015, 0.005] -X = FiniteRV('X', dict(zip(Xs, Ps))) -E(X) -``` - -```{python} -Xs = [0, 100, 1000, 10000]; Ps = [0.9, 0.08, 0.015, 0.005] -rv = stats.rv_discrete(values=(Xs, Ps)) -rv.mean() -``` - -```{python} -sum(x * density(X)[x] for x in Xs) -``` - -```{python} -sum(x * rv.pmf(x) for x in Xs) -``` - -```{python} -np.dot(Xs, Ps) -``` - -```{python} -np.mean(list(sample_iter(X, numsamples=500000))) -``` - -```{python} -rv.rvs(size=500000).mean() -``` - -```{python} -var('x') -X = ContinuousRV(x, Abs(x), set=Interval(-1, 1)) -integrate(x * density(X)(x), (x, -1, 1)) -``` - -```{python} -class MyX(stats.rv_continuous): - def _pdf(self, x): return abs(x) -rv = MyX(a=-1, b=1) -quad(lambda x: x * rv.pdf(x), -1, 1) -``` - -```{python} -Xs = [0, 100, 1000, 10000]; Ps = [0.9, 0.08, 0.015, 0.005] -X = FiniteRV('X', dict(zip(Xs, Ps))) -variance(X) -``` - -```{python} -Xs = [0, 100, 1000, 10000]; Ps = [0.9, 0.08, 0.015, 0.005] -rv = stats.rv_discrete(values=(Xs, Ps)) -rv.var() -``` - -```{python} -E((X - E(X))**2) -``` - -```{python} -sum((x - E(X))**2 * density(X)[x] for x in Xs) -``` - -```{python} -sum((x - rv.mean())**2 * rv.pmf(x) for x in Xs) -``` - -```{python} -np.dot((Xs - np.dot(Xs, Ps))**2, Ps) -``` - -```{python} -var('x') -X = ContinuousRV(x, Abs(x), set=Interval(-1, 1)) -integrate((x - E(X))**2 * density(X)(x), (x, -1, 1)) -``` - -```{python} -class MyX(stats.rv_continuous): - def _pdf(self, x): return abs(x) -rv = MyX(a=-1, b=1) -quad(lambda x: (x - rv.mean())**2 * rv.pdf(x), -1, 1) -``` - -# 9 多次元の確率分布 - -```{python} -X1 = DiscreteUniform('X1', range(1, 7)) -X2 = DiscreteUniform('X2', range(1, 7)) -X, Y = Max(X1, X2), Min(X1, X2) -{(x, y): P(And(Eq(X, x), Eq(Y, y))) for x in range(1, 7) for y in range(1, 7)} -``` - -```{python} -density(X), density(Y) -``` - -```{python} -{(x, y): P(And(Le(X, x), Le(Y, y))) for x in range(1, 7) for y in range(1, 7)} -``` - -```{python} -X1 = DiscreteUniform('X1', range(1, 7)) -X2 = DiscreteUniform('X2', range(1, 7)) -X, Y = Max(X1, X2), Min(X1, X2) - -(E(X), E(Y), # 平均 - variance(X), variance(Y), # 分散 - std(X), std(Y), # 標準偏差 - covariance(X, Y), # 共分散 - correlation(X, Y)) # 相関係数 -``` - -```{python} -uX, uY = E(X), E(Y); sX, sY = std(X), std(Y) -(E(X), E(Y), # 平均 - E((X - uX)**2), E((Y - uY)**2), # 分散 - E((X - uX) * (Y - uY)), # 共分散 - E((X - uX) * (Y - uY) / sX / sY)) # 相関係数 -``` - -```{python} -print(sum((x * P(Eq(X, x)) for x in range(1, 7)))) # 平均 -print(sum((x - uX) * (y - uY) * P(And(Eq(X, x), Eq(Y, y))) # 共分散 - for x in range(1, 7) for y in range(1, 7))) -``` - -```{python} -[P(And(X <= x, Y <= y)) for x in range(1, 7) for y in range(1, 7)] == \ -[P(X <= x) * P(Y <= y) for x in range(1, 7) for y in range(1, 7)] -``` - -```{python} -U = DiscreteUniform('X', range(1, 7)); X = U % 2; Y = U % 3 -[P(And(X <= x, Y <= y)) for x in range(2) for y in range(3)] == \ -[P(X <= x) * P(Y <= y) for x in range(2) for y in range(3)] -``` - -```{python} -X = Binomial('X', 3, sym.S(1) / 2) -Y = Piecewise((1, Or(Eq(X, 0), Eq(X, 3))), (2, True)) -covariance(X, Y) -``` - -```{python} -[P(And(X <= x, Y <= y)) for x in range(4) for y in (1, 2)] == \ -[P(X <= x) * P(Y <= y) for x in range(4) for y in (1, 2)] -``` - -```{python} -X1 = DiscreteUniform('X1', range(1, 7)) -X2 = DiscreteUniform('X2', range(1, 7)) -X, Y = Max(X1, X2), Min(X1, X2) - -(E(X + Y), E(X) + E(Y), # 平均 - variance(X + Y), variance(X) + variance(Y) + 2 * covariance(X, Y)) # 分散 -``` - -```{python} -n = 15; p = sym.S(4) / 10; Y = Binomial('Y', n, p) -mu = E(Y); sigma = std(Y); Z = Normal('Z', mu, sigma) -x1 = range(0, n + 1); y1 = [density(Y)(x) for x in x1] -x2 = np.linspace(0, n, 101); y2 = [density(Z)(t) for t in x2] -_, ax = plt.subplots(); ax.scatter(x1, y1); ax.plot(x2, y2); -``` - -```{python} -n = 15; p = 4 / 10; Y = stats.binom(n, p) -mu = Y.mean(); sigma = Y.std(); Z = stats.norm(mu, sigma) -x1 = range(0, n + 1); y1 = Y.pmf(x1) -x2 = np.linspace(0, n, 101); y2 = Z.pdf(x2) -_, ax = plt.subplots(); ax.scatter(x1, y1); ax.plot(x2, y2); -``` - -```{python} -X = Uniform('X', 0, 1); Z = Normal('Z', 0, 1) -data = [sum(sample_iter(X, numsamples=12)) - 6 for _ in range(10000)] -x = np.linspace(-4, 4, 101); y = [density(Z)(t) for t in x] -_, ax = plt.subplots(); ax.hist(data, density=True, alpha=0.5); ax.plot(x, y); -``` - -```{python} -X = stats.uniform(); Z = stats.norm() -data = [sum(X.rvs(12)) - 6 for _ in range(10000)] -x = np.linspace(-4, 4, 101); y = Z.pdf(x) -_, ax = plt.subplots(); ax.hist(data, density=True, alpha=0.5); ax.plot(x, y); -``` - -```{python} -mu = [3, 6]; Sigma = [[5, 7], [7, 13]]; -rv = stats.multivariate_normal(mu, Sigma) -Y1, Y2 = np.mgrid[-8:14:0.1, -5:17:0.1] -grid = np.dstack((Y1, Y2)) -z = rv.pdf(grid) -fig = plt.figure() -ax = fig.add_subplot(111, projection='3d') -ax.plot_surface(Y1, Y2, z, cmap='viridis') -ax.set_xlabel('Y1'); ax.set_ylabel('Y2'); -``` - -```{python} -plt.contourf(Y1, Y2, z, cmap='viridis'); plt.xlabel('Y1'); plt.ylabel('Y2'); -``` - -# 10 推測統計 - -```{python} -mu = 2; sigma = 3; rv = stats.norm(mu, sigma) -data1 = [rv.rvs(5).mean() for _ in range(10000)] -data2 = [rv.rvs(50).mean() for _ in range(10000)] -((np.mean(data1), np.var(data1)), (np.mean(data2), np.var(data2))) -``` - -```{python} -plt.hist(data1, density=True), plt.hist(data2, alpha=0.7, density=True); -``` - -```{python} -mu = 2; sigma = 3; rv = stats.norm(mu, sigma) -data1 = [rv.rvs(5).var(ddof=1) for _ in range(10000)] -data2 = [rv.rvs(50).var(ddof=1) for _ in range(10000)] -((np.mean(data1), np.var(data1)), (np.mean(data2), np.var(data2))) -``` - -```{python} -plt.hist(data1, density=True), plt.hist(data2, alpha=0.7, density=True); -``` - -```{python} -n = 4; mu = 5; sigma = 7; rv = stats.norm(mu, sigma) -f = lambda x: (n - 1) * x.var(ddof=1) / sigma**2 -data = [f(rv.rvs(n)) for _ in range(10000)] -x = np.linspace(0, 20, 101); y = stats.chi2(n - 1).pdf(x) -_, ax = plt.subplots() -ax.hist(data, bins='sturges', density=True, alpha=0.5) -ax.plot(x, y); -``` - -```{python} -n = 4; mu = 5; sigma = 7; rv = stats.norm(mu, sigma) -t = lambda x: (x.mean() - mu) / np.sqrt(x.var(ddof=1) / n) -data = [t(rv.rvs(n)) for _ in range(10000)] -x = np.linspace(-5, 5, 101); y = stats.t(n - 1).pdf(x) -_, ax = plt.subplots() -ax.hist(data, density=True, bins=np.arange(-5, 5, 0.5), alpha=0.5) -ax.plot(x, y); -``` - -```{python} -m = 5; muX = 2; sigmaX = 3; rvX = stats.norm(muX, sigmaX) -n = 7; muY = 3; sigmaY = 2; rvY = stats.norm(muY, sigmaY) -f = lambda x, y: (x.var(ddof=1) / sigmaX**2) / (y.var(ddof=1) / sigmaY**2) -data = [f(rvX.rvs(m), rvY.rvs(n)) for _ in range(10000)] -x = np.linspace(0, 5, 101); y = stats.f(m - 1, n - 1).pdf(x) -_, ax = plt.subplots() -ax.hist(data, density=True, bins=np.arange(0, 5.2, 0.2), alpha=0.5) -ax.plot(x, y); -``` - -```{python} -n = 15; p0 = 4 / 10; binom_test(2, 15, 4 / 10) -``` - -```{python} -binom_test(2, n, p0, alternative="smaller") -``` - -```{python} -n = 15; p0 = 4 / 10; mu0 = n * p0; sigma0 = np.sqrt((n * p0 * (1 - p0))) -2 * stats.norm.cdf(2, mu0, sigma0) -``` - -```{python} -alpha = 5 / 100; stats.norm.ppf((alpha/2, 1 - alpha/2), mu0, sigma0) -``` - -```{python} -x = [24.2, 25.3, 26.2, 25.7, 24.4, 25.1, 25.6]; d = DescrStatsW(x); mu0 = 25 -d.ttest_mean(mu0) -``` - -```{python} -m = np.mean(x); s2 = np.var(x, ddof=1); n = len(x); -t = (m - mu0) / np.sqrt(s2 / n); c = stats.t.cdf(t, n - 1) -2 * min(c, 1 - c) -``` - -```{python} -alpha = 5 / 100 -stats.t.ppf((alpha / 2, 1 - alpha / 2), n - 1) -``` - -```{python} -d.tconfint_mean() -``` - -```{python} -x = [25, 24, 25, 26]; y = [23, 18, 22, 28, 17, 25, 19, 16] -ttest_ind(x, y, alternative="larger", usevar='unequal') -``` - -```{python} -tmp = CompareMeans(DescrStatsW(x), DescrStatsW(y)) -tmp.tconfint_diff(alternative="larger", usevar='unequal') -``` - -```{python} -x = [25, 24, 25, 26]; y = [23, 18, 22, 28, 17, 25, 19, 16] -f = np.var(x, ddof=1) / np.var(y, ddof=1); m = len(x); n = len(y); -c = stats.f.cdf(f, m - 1, n - 1) -f, 2 * min(c, 1 - c) -``` - -```{python} -alpha = 0.05 -stats.f.ppf((alpha / 2, 1 - alpha / 2), m - 1, n - 1) -``` - -# 11 線形回帰分析 - -```{python} -data = pd.DataFrame({ - 'x1': [1, 1, 2, 3], 'x2': [2, 3, 5, 7], 'y': [3, 6, 3, 6]}) -model = smf.ols('y ~ x1 + x2', data).fit() -print(model.summary2(alpha=0.05)) -# Results: Ordinary least squares -# ================================================================ -# Model: OLS Adj. R-squared: -1.000 -# Dependent Variable: y AIC: 18.9734 -# Date: 2023-01-01 00:00 BIC: 17.1323 -# No. Observations: 4 Log-Likelihood: -6.4867 -# Df Model: 2 F-statistic: 0.2500 -# Df Residuals: 1 Prob (F-statistic): 0.816 -# R-squared: 0.333 Scale: 6.0000 -# ----------------------------------------------------------------- -# Coef. Std.Err. t P>|t| [0.025 0.975] -# ----------------------------------------------------------------- -# const 3.0000 3.0000 1.0000 0.5000 -35.1186 41.1186 -# x1 -4.0000 7.6811 -0.5208 0.6943 -101.5982 93.5982 -# x2 2.0000 3.3166 0.6030 0.6545 -40.1417 44.1417 -# ---------------------------------------------------------------- -# Omnibus: nan Durbin-Watson: 3.167 -# Prob(Omnibus): nan Jarque-Bera (JB): 0.611 -# Skew: -0.816 Prob(JB): 0.737 -# Kurtosis: 2.000 Condition No.: 35 -# ================================================================ -``` - -```{python} -model.predict({'x1': 1.5, 'x2': 4}) -``` - -```{python} -x1 = np.array([1, 3, 6, 10]); y = np.array([7, 1, 6, 14]) -def L(b): - e = y - (b[0] + b[1] * x1) - return e @ e # 内積 -minimize(L, x0=[0, 0]) -``` - -```{python} -data = pd.DataFrame({ - 'x1': [1, 1, 2, 3], 'x2': [2, 3, 5, 7], 'y': [3, 6, 3, 6]}) -y, X = dmatrices('y ~ x1 + x2', data) -linalg.inv(X.T @ X) @ X.T @ y -``` - -```{python} -linalg.pinv(X) @ y -``` - -```{python} -data = pd.DataFrame({ - 'x1': [1, 1, 2, 3], 'x2': [2, 3, 5, 7], 'y': [3, 6, 3, 6]}) -model = smf.ols('y ~ x1 + x2', data).fit() -model.rsquared -``` - -```{python} -model.rsquared_adj -``` - -```{python} -x1 = [1, 3, 6, 10]; y = [7, 1, 6, 14] -data = pd.DataFrame({'x1': x1, 'y': y}) -_, X = dmatrices('y ~ x1', data) -yh = X @ linalg.pinv(X) @ y -eh = y - yh; fh = yh - np.mean(y); g = y - np.mean(y) -R2 = 1 - np.dot(eh, eh) / np.dot(g, g); R2 -``` - -```{python} -(np.allclose(eh.mean(), 0), # 特徴1 - np.allclose(yh.mean(), data.y.mean()), # 特徴2 - np.allclose(np.dot(g, g), np.dot(fh, fh) + np.dot(eh, eh)), # 特徴3 - np.allclose(R2, np.dot(fh, fh) / np.dot(g, g)), # 特徴4 - np.allclose(R2, np.corrcoef(y, yh)[0, 1]**2), # 特徴5 - 0 <= R2 <= 1, # 特徴6 - np.allclose(np.corrcoef(y, yh), np.corrcoef(y, x1))) # 特徴7 -``` - -```{python} -data = pd.DataFrame({ - 'x1': [1, 1, 2, 3], 'x2': [2, 3, 5, 7], 'y': [3, 6, 3, 6]}) -model = smf.ols('y ~ x1 + x2', data).fit() -model.scale -``` - -```{python} -data = pd.DataFrame({ - 'x1': [1, 1, 2, 3], 'x2': [2, 3, 5, 7], 'y': [3, 6, 3, 6]}) -model = smf.ols('y ~ x1 + x2', data).fit() -model.f_pvalue -``` - -```{python} -data = pd.DataFrame({'x1': [35, 45, 55, 65, 75], - 'y': [114, 124, 143, 158, 166]}) -model = smf.ols('y ~ x1', data).fit() -print(model.summary2(alpha=0.05)) # 信頼区間はこのレポートにも表示される． -model.conf_int(alpha=0.05) -``` - -```{python} -data = pd.DataFrame({'x1': [35, 45, 55, 65, 75], - 'y': [114, 124, 143, 158, 166]}) -model = smf.ols('y ~ x1', data).fit() -tmp = pd.DataFrame({'x1': np.linspace(35, 75, 100)}) -ci = model.get_prediction(tmp).summary_frame(alpha=0.05) -df = ci[['mean', 'mean_ci_lower', 'mean_ci_upper']].assign(x1=tmp.x1) -_, ax = plt.subplots() -ax.scatter(data.x1, data.y) -df.plot(x='x1', ax=ax); -``` - -```{python} -df = ci[['mean', 'obs_ci_lower', 'obs_ci_upper']].assign(x1=tmp.x1) -_, ax = plt.subplots() -ax.scatter(data.x1, data.y) -df.plot(x='x1', ax=ax); -``` - -```{python} -sns.regplot(x=data.x1, y=data.y, ci=95); -``` - -# 12 関数の極限と連続性 - -```{python} -f = lambda x: 2 * x - 3; var('x') -limit(f(x), x, 1) -``` - -```{python} -limit(2 * x - 3, x, 1) -``` - -```{python} -g = lambda x: (x**2 - 2) / (x - sqrt(2)) -limit(g(x), x, sqrt(2)) -``` - -```{python} -limit((1 + 1 / x)**x, x, oo) -``` - -```{python} -limit(1 / x**2, x, 0) -``` - -```{python} -limit(abs(x) / x, x, 0, dir='+'), limit(abs(x) / x, x, 0, dir='-') -``` - -# 13 微分 - -```{python} -f = lambda x: x**3; var('h'); a = 1 -limit((f(a + h) - f(a)) / h, h, 0) -``` - -```{python} -f = lambda x: x**3; var('x') -diff(f(x), x) -``` - -```{python} -f = lambda x: x**3; var('x') -f1 = Lambda(x, diff(f(x), x)) -f2 = diff(f(x), x) -f1, f2 -``` - -```{python} -f1(1), f2.subs(x, 1) -``` - -```{python} -diff(x**3, x, 2) -``` - -```{python} -var('a b t x') -f = lambda t: t**2 -g = lambda x: a * x + b -fp = Lambda(t, diff(f(t), t)); gp = Lambda(x, diff(g(x), x)) -(Lambda(x, diff(Lambda(x, f(g(x)))(x), x))(x), # ① - diff(f(g(x)), x), # ② - diff(f(t), t).subs(t, g(x)) * diff(g(x), x), # ③ - fp(g(x)) * gp(x)) # ④ -``` - -```{python} -var('x'); tmp = series(sin(x), x, x0=0, n=6); tmp -``` - -```{python} -plot(*(tmp.removeO(), sin(x)), (x, -pi, pi)); -``` - -```{python} -a = 0; -np.sum([diff(sin(x), x, k).subs(x, a) * (x - a)**k / factorial(k) - for k in range(6)]) -``` - -```{python} -f = lambda x: x**3 - 12 * x; var('x') -sol = solve(diff(f(x), x), x); print(sol) -series(f(x), x0=sol[0], n=3) -``` - -# 14 積分 - -```{python} -var('x') -integrate(-x**2 + 4 * x + 1, (x, 1, 4)) -``` - -```{python} -f = lambda x: -x**2 + 4 * x + 1; var('k n') -a = 1; b = 4; h = (b - a) / n -s = simplify(Sum(f(a + h * k) * h, (k, 1, n)).doit()); print(s) -limit(s, n, oo) -``` - -```{python} -var('a t x'); integrate(-t**2 + 4 * t + 1, (t, a, x)) -``` - -```{python} -integrate(-x**2 + 4 * x + 1, x) -``` - -```{python} -var('x'); y = Function('y') -dsolve(Eq(diff(y(x), x), -x**2 + 4 * x + 1), y(x)) -``` - -```{python} -dsolve(Eq(diff(y(x), x), -x**2 + 4 * x + 1), y(x), ics={y(0): 1}) -``` - -```{python} -tmp = dsolve(Eq(diff(y(x), x), -x * y(x)), y(x)); tmp -``` - -```{python} -solve(Eq(integrate(tmp.rhs, (x, -oo, oo)), 1), 'C1') -``` - -```{python} -f = Function('f'); var('a t x') -Lambda(x, diff(Lambda(x, integrate(f(t), (t, a, x)))(x), x)) -``` - -```{python} -diff(integrate(f(t), (t, a, x)), x) -``` - -```{python} -var('x') -F = integrate(-x**2 + 4 * x + 1, x) -F.subs(x, 4) - F.subs(x, 1) -``` - -```{python} -var('x'); integrate((x**2 + 2 * x + 1 + (3 * x + 1) * sqrt(x + log(x))) / (x * sqrt(x + log(x)) * (x + sqrt(x + log(x)))), x) -``` - -```{python} -var('a x') -simplify(integrate(1/x**a, (x, 0, 1))) -``` - -```{python} -integrate(1/x**a, (x, 1, oo)) -``` - -```{python} -var('x'); integrate(exp(-x**2), (x, -oo, oo)) -``` - -# 15 多変数関数の微分積分 - -```{python} -f = lambda x, y: x**2 * y / (x**4 + y**2); var('x y r theta') -(limit(limit(f(x, y), x, 0), y, 0), # ① - limit(limit(f(x, y), y, 0), x, 0), # ② - limit(f(r * cos(theta), r * sin(theta)), r, 0), # ③ - limit(f(x, x**2), x, 0)) # ④ -``` - -```{python} -f = lambda x, y: 2 - x**2 - y**2; var('x y') -(diff(f(x, y), x), diff(f(x, y), y)) -``` - -```{python} -diff(f(x, y), Matrix([x, y])) -``` - -```{python} -f = lambda x, y: 2 * x**3 + 5 * x * y + 2 * y**2; var('x y') -hessian(f(x, y), Matrix([x, y])) -``` - -```{python} -f = lambda x: sqrt(x[0]**2 + x[1]**2); var('x1 x2 t'); -x = Matrix([x1, x2]); a = Matrix([1, 1]); h = x - a; -F = lambda t: f(a + t * h) -expr = series(F(t), t, x0=0, n=3).removeO().subs(t, 1) -simplify(expr) -``` - -```{python} -gradf = diff(f(x), x).subs(zip(x, a)) -H = hessian(f(x), x).subs(zip(x, a)) -simplify(f(a) + gradf.dot(x - a) + (x - a).dot(H @ (x - a)) / 2) -``` - -```{python} -f = lambda x: 2 * x[0]**3 + x[0] * x[1]**2 + 5 * x[0]**2 + x[1]**2 -var('x1 x2', real=True); x = Matrix([x1, x2]) # x1, x2は実数 -points = solve(diff(f(x), x), x) # 停留点 -H = hessian(f(x), x) # ヘッセ行列 -def check(H, a): - h = H.subs(zip(x, a)) # 停留点でのヘッセ行列 - if h.is_positive_definite: return (a, f(a), -1) # 極小 - if h.is_negative_definite: return (a, f(a), 1) # 極大 - if h.is_indefinite: return (a, f(a), 0) # 極値ではない - else: return (a, f(a), None) # 不明 -[check(H, a) for a in points] -``` - -```{python} -f = lambda x, y: x**2 + y**2; var('x y') -integrate(integrate(f(x, y), (y, 0, x)), (x, 0, 1)) -``` - -```{python} -integrate(integrate(f(x, y), (x, y, 1)), (y, 0, 1)) -``` - -```{python} -f = lambda x, y: x**2 + y**2; var('u v') -x, y = 2 * u, 3 * v -J = Matrix([x, y]).jacobian((u, v)) -detJ = J.det(); print(detJ) -integrate(integrate(f(x, y) * abs(detJ), - (v, 0, 2 * u / 3)), (u, 0, sym.S(1) / 2)) -``` - -```{python} -var('r theta'); x, y = r * cos(theta), r * sin(theta) -J = Matrix([x, y]).jacobian((r, theta)) -simplify(J.det()) -``` - -# 16 ベクトル - -```{python} -a = Matrix([2, 3, 5]); len(a) -``` - -```{python} -10 * Matrix([2, 3]) -``` - -```{python} -u = Matrix([10, 20]); v = Matrix([2, 3]); u + v -``` - -```{python} -t = sym.S(10) -a = Matrix([1 / t + 2 / t, 1 / t + 2 / t - 3 / t]); b = Matrix([3 / t, 0]) -a == b -``` - -```{python} -a = np.array([1/10 + 2/10, 1/10 + 2/10 - 3/10]); b = np.array([3/10, 0]) -np.allclose(a, b) -``` - -```{python} -100 * Matrix([1, 2]) + 10 * Matrix([3, 1]) -``` - -```{python} -100 * np.array([1, 2]) + 10 * np.array([3, 1]) -``` - -```{python} -a = Matrix([3, 4]) -a.norm() -``` - -```{python} -a = np.array([3, 4]) -linalg.norm(a) -``` - -```{python} -var('x y'); a = Matrix([x, y]); sqrt(a.dot(a)) -``` - -```{python} -var('x y', real=True); a = Matrix([x, y]); a.norm() -``` - -```{python} -a = Matrix([3, 4]) -a.normalized() -``` - -```{python} -a = np.array([3, 4]) -a / linalg.norm(a) -``` - -```{python} -a = Matrix([1, 0]); b = Matrix([1, 1]) -acos(a.dot(b) / (a.norm() * b.norm())) -``` - -```{python} -a = np.array([1, 0]); b = np.array([1, 1]) -np.arccos(a @ b / (linalg.norm(a) * linalg.norm(b))) -``` - -# 17 行列 - -```{python} -A = Matrix([[1, 2, 0], [0, 3, 4]]); A -``` - -```{python} -A = np.array([[1, 2, 0], [0, 3, 4]]); A -``` - -```{python} -x = [5, 7]; diag(*x) # diag(x)はNG．diag(5, 7)はOK． -``` - -```{python} -x = [5, 7]; np.diag(x) # np.diag(*x)はNG．np.diag([5, 7])はOK． -``` - -```{python} -Matrix([[1, 2], [2, 3]]).is_symmetric() -``` - -```{python} -A = np.array([[1, 2], [2, 3]]) -np.allclose(A.T, A) -``` - -```{python} -A = Matrix([[11, 12, 13], [21, 22, 23], [31, 32, 33]]); A -``` - -```{python} -An = np.array([[11, 12, 13], [21, 22, 23], [31, 32, 33]]); An -``` - -```{python} -A[0:2, 0:2], A[:2, :2] # 二つの方法 -``` - -```{python} -An[0:2, 0:2], An[:2, :2] # 二つの方法 -``` - -```{python} -A[:, 2] -``` - -```{python} -An[:, 2] -``` - -```{python} -An[:, [2]] -``` - -```{python} -A[1, :] -``` - -```{python} -An[1, :] -``` - -```{python} -An[[1], :] -``` - -```{python} -10 * Matrix([[2, 3], [5, 7]]) -``` - -```{python} -10 * np.array([[2, 3], [5, 7]]) -``` - -```{python} -Matrix([[10, 20], [30, 40]]) + Matrix([[2, 3], [4, 5]]) -``` - -```{python} -np.array([[10, 20], [30, 40]]) + np.array([[2, 3], [4, 5]]) -``` - -```{python} -A = Matrix([[2, 3], [5, 7]]); B = Matrix([[1, 2], [3, 4]]) -A @ B -``` - -```{python} -A = np.array([[2, 3], [5, 7]]); B = np.array([[1, 2], [3, 4]]) -A @ B -``` - -```{python} -A = Matrix([[2, 3], [5, 7]]); B = Matrix([[1, 2, 3], [4, 5, 6]]); S = A @ B -p, q = A.shape; r, s = B.shape -S1 = Matrix([[A[i, :].dot(B[:, j]) for j in range(s)] for i in range(p)]) # ① -S2 = sum((A[:, j] @ B[j, :] for j in range(q)), zeros(p, s)) # ② -S3 = Matrix.hstack(*[A @ B[:, j] for j in range(s)]) # ③ -S4 = Matrix.vstack(*[A[i, :] @ B for i in range(p)]) # ④ -S == S1, S == S2, S == S3, S == S4 -``` - -```{python} -A = np.array([[2, 3], [5, 7]]); B = np.array([[1, 2, 3], [4, 5, 6]]); S = A @ B -p, q = A.shape; r, s = B.shape -S1 = np.array( - [[A[i, :].dot(B[:, j]) for j in range(s)] for i in range(p)]) # ① -S2 = sum((A[:, [j]] @ B[[j], :] for j in range(q)), np.zeros([p, s])) # ② -S3 = np.vstack([A @ B[:, j] for j in range(s)]).T # ③ -S4 = np.vstack([A[i, :] @ B for i in range(p)]) # ④ -np.allclose(S, S1), np.allclose(S, S2), np.allclose(S, S3), np.allclose(S, S4) -``` - -```{python} -var('a1 a2 x1 x2 p q r s') -x = Matrix([x1, x2]); a = Matrix([a1, a2]) -G = Matrix([[p, q], [q, s]]); A = Matrix([[p, q], [r, s]]) -(Eq(diff(a.dot(x), x), a), - Eq(diff(x.dot(G @ x), x), 2 * G @ x), - Eq(simplify(diff((A @ x).dot(A @ x), x) - 2 * A.T @ A @ x), zeros(2, 1))) -``` - -```{python} -Matrix([[3, 2], [1, 2]]).det() -``` - -```{python} -linalg.det(np.array([[3, 2], [1, 2]])) -``` - -```{python} -Matrix([[2, 3], [5, 7]]).inv() -``` - -```{python} -linalg.inv(np.array([[2, 3], [5, 7]])) -``` - -```{python} -A = Matrix([[3, 2], [1, 2]]); b = Matrix([8, 4]) -A.inv() @ b -``` - -```{python} -A = np.array([[3, 2], [1, 2]]); b = np.array([8, 4]) -linalg.inv(A) @ b -``` - -```{python} -Matrix([[4, 2, 8], [2, 1, 4]]).rref()[0] -``` - -```{python} -A = Matrix([[2, 0, 2], [0, 2, -2], [2, 2, 0]]) -A.rank() -``` - -```{python} -A = np.array([[2, 0, 2], [0, 2, -2], [2, 2, 0]]) -np.linalg.matrix_rank(A) -``` - -# 18 ベクトル空間 - -```{python} -a1 = Matrix([3, 1]); a2 = Matrix([2, 2]); var('c1 c2') -solve(c1 * a1 + c2 * a2, (c1, c2)) -``` - -```{python} -A = Matrix([[1, 0, 1], [1, 1, 0], [0, 1, -1]]) -A.columnspace() -``` - -```{python} -A = Matrix([[1, 0, 1], [1, 1, 0], [0, 1, -1]]) -tmp = A.columnspace() -basis = GramSchmidt(tmp, orthonormal=True); basis -``` - -```{python} -Q = Matrix.hstack(*basis) -Q.T @ Q -``` - -```{python} -A = Matrix([[1, 2], [1, 2], [0, 0]]); B = Matrix([[1, 0], [1, 1], [0, 1]]) -Qa, Ra = A.QRdecomposition() -Qb, Rb = B.QRdecomposition() -display(Qa, Ra, A == Qa @ Ra, Qb, Rb, B == Qb @ Rb) -``` - -```{python} -A = np.array([[1, 2], [1, 2], [0, 0]]); B = np.array([[1, 0], [1, 1], [0, 1]]) -Qa, Ra = linalg.qr(A) -Qb, Rb = linalg.qr(B) -Qa, Ra, np.allclose(A, Qa @ Ra), Qb, Rb, np.allclose(B, Qb @ Rb) -``` - -```{python} -def qrd(A): - (m, n) = A.shape - u = [A[:, i].copy() for i in range(n)]; idx = [] - for i in range(n): - for j in range(i): u[i] -= A[:, i].dot(u[j]) * u[j] - s = u[i].norm() - if not np.isclose(np.double(s), 0): u[i] /= s; idx.append(i) - Q = Matrix.hstack(*[u[i] for i in idx]) if len(idx) != 0 else eye(m) - return (Q, Q.T @ A) - -A = Matrix([[1, 2], [1, 2], [0, 0]]); B = Matrix([[1, 0], [1, 1], [0, 1]]) -qrd(A), qrd(B) # 動作確認 -``` - -```{python} -def qrdn(An): - A = An.astype(float) # 成分を浮動小数点数に変換する． - (m, n) = A.shape - u = [A[:, i].copy() for i in range(n)]; idx = [] - for i in range(n): - for j in range(i): u[i] -= A[:, i].dot(u[j]) * u[j] - s = linalg.norm(u[i]) - if not np.isclose(s, 0): u[i] /= s; idx.append(i) - Q = np.array([u[i] for i in idx]).T if len(idx) != 0 else np.eye(m) - return (Q, Q.T @ A) - -A = np.array([[1, 2], [1, 2], [0, 0]]); B = np.array([[1, 0], [1, 1], [0, 1]]) -qrdn(A), qrdn(B) # 動作確認 -``` - -```{python} -B = Matrix([[1, 0], [1, 1], [0, 1]]) -Q, R = qrd(B) # QR分解 -(Q.T @ Q == eye(Q.shape[1]), # ①（厳密値を想定） - R.is_upper, # ②（厳密値を想定） - Q @ R == B) # ③（厳密値を想定） -``` - -```{python} -B = np.array([[1, 0], [1, 1], [0, 1]]) -Q, R = qrdn(B) # QR分解 -(np.allclose(Q.T @ Q, np.eye(Q.shape[1])), # ① - np.allclose(R, np.triu(R)), # ② - np.allclose(Q @ R, B)) # ③ -``` - -```{python} -A = Matrix([[1, 0], [1, 1], [0, 1]]) -A.T.nullspace() -``` - -```{python} -A = np.array([[1, 0], [1, 1], [0, 1]]) -linalg.null_space(A.T) # 正規直交基底 -``` - -```{python} -A = Matrix([[1, 0], [1, 1], [0, 1]]) -basis1 = A.QRdecomposition()[0] # 列空間 -basis2 = GramSchmidt(A.T.nullspace(), orthonormal=True) # 直交補空間 -Q = Matrix.hstack(basis1, *basis2) -Q, Q.T @ Q == eye(3) -``` - -```{python} -A = np.array([[1, 0], [1, 1], [0, 1]]) -basis1 = linalg.qr(A, mode='economic')[0] # 列空間 -basis2 = linalg.null_space(A.T) # 直交補空間 -Q = np.hstack([basis1, basis2]) -Q, np.allclose(Q.T @ Q, np.eye(3)) -``` - -# 19 固有値と固有ベクトル - -```{python} -A = Matrix([[5, 6, 3], [0, 9, 2], [0, 6, 8]]) -eigs = A.eigenvects() -vals = [e[0] for e in eigs for _ in range(e[1])] # 固有値（順序不明） -vecs = [v for e in eigs for v in e[2]] # 固有ベクトル（非正規） -vals, vecs -``` - -```{python} -An = np.array([[5, 6, 3], [0, 9, 2], [0, 6, 8]]) -valsn, vecsn = linalg.eig(An) -valsn, vecsn # 固有値（順序不明），固有ベクトル（正規） -``` - -```{python} -V = Matrix.hstack(*vecs); A @ V == V @ diag(*vals) -``` - -```{python} -np.allclose(An @ vecsn, vecsn @ np.diag(valsn)) -``` - -```{python} -B = Matrix([[6, 7, 0], [7, 2, 0], [0, 0, 8]]) -eigs = B.eigenvects() -vals = [e[0] for e in eigs for _ in range(e[1])] # 絶対値の降順とは限らない． -vecs = [v for e in eigs for v in e[2]] -idx = np.argsort([-abs(N(x)) for x in vals]) # 近似値で比較する． -vals = [vals[i] for i in idx] # 固有値（絶対値の降順） -vecs = [vecs[i] for i in idx] # 固有ベクトルも並べ替える． -vals # 結果の確認 -``` - -```{python} -B = np.array([[6, 7, 0], [7, 2, 0], [0, 0, 8]]) -vals, vecs = np.linalg.eig(B) # 絶対値の降順とは限らない． -idx = np.argsort(-abs(vals)) # 絶対値の降順の位置 -vals, vecs = vals[idx], vecs[:, idx] # 絶対値の降順 -vals # 結果の確認 -``` - -```{python} -A = Matrix([[5, 6, 3], [0, 9, 2], [0, 6, 8]]); n = A.shape[0] -var('x'); solve((x * eye(n) - A).det(), x) -``` - -```{python} -(5 * eye(n) - A).nullspace() -``` - -```{python} -S = Matrix([[2, 2, -2], [2, 5, -4], [-2, -4, 5]]) -Q, L, V = S.singular_value_decomposition() -Q, L, S == Q @ L @ Q.T == V @ L @ V.T -``` - -```{python} -Sn = np.array([[2, 2, -2], [2, 5, -4], [-2, -4, 5]]) -Q, lambdas, tV = linalg.svd(Sn); L = np.diag(lambdas) -Q, L, np.allclose(Sn, Q @ L @ Q.T), np.allclose(Sn, tV.T @ L @ tV) -``` - -```{python} -S = Matrix([[2, 2, -2], [2, 5, -4], [-2, -4, 5]]) -eigs = S.eigenvects() -vals = [e[0] for e in eigs for _ in range(e[1])] # ①-1 -vecs = [v for e in eigs for v in e[2]] # ①-2 -tmp = GramSchmidt(vecs, orthonormal=True) # ② -Q = Matrix.hstack(*tmp) # ③ -L = diag(*vals) # ④ -Q, L, np.allclose(np.double(S), np.double(Q @ L @ Q.T)) # 近似的な比較 -``` - -```{python} -Sn = np.array([[2, 2, -2], [2, 5, -4], [-2, -4, 5]]) -tmp = linalg.eig(Sn); vals = tmp[0]; vecs = tmp[1] # ① -Q = linalg.qr(vecs)[0] # ②, ③ -L = np.diag(vals) # ④ -Q, L, np.allclose(Sn, Q @ L @ Q.T) -``` - -```{python} -Matrix([[4, 2], [2, 1]]).is_positive_semidefinite -``` - -```{python} -A = Matrix([[4, 2], [2, 1]]) -all(e >= 0 for e in A.eigenvals(multiple=True)) -``` - -```{python} -A = np.array([[4, 2], [2, 1]]) -np.all(linalg.eigvals(A) >= 0) -``` - -```{python} -x1 = [1, 3, 6, 10]; y = [7, 1, 6, 14]; X = Matrix([x1, y]).T -n = X.shape[0]; M = ones(n, n) / sym.S(n) -A = X - M @ X -S = A.T @ A; display(S) -eigs = S.eigenvects() -vals = [re(e[0]) for e in eigs for _ in range(e[1])] # 固有値 -vecs = [re(v) for e in eigs for v in e[2]] # 固有ベクトル -N(vecs[np.argmax(vals)].normalized()) # 最大固有値に対応する固有ベクトル -# 最大固有値の位置の探索と，固有ベクトルの正規化が必要 -``` - -```{python} -x1 = [1, 3, 6, 10]; y = [7, 1, 6, 14]; Xn = np.array([x1, y]).T -n = Xn.shape[0]; Mn = np.ones([n, n]) / n -An = Xn - Mn @ Xn -Sn = An.T @ An; print(Sn) -vals, vecs = linalg.eig(Sn) -vecs[:, np.argmax(vals)] # 最大固有値に対応する固有ベクトル -# 最大固有値の位置の探索が必要 -``` - -```{python} -_, L, V = A.singular_value_decomposition() # 特異値分解 -s = re(L.diagonal()) # 特異値 -idx = np.argmax(s) # 最大特異値の位置（SymPyでは0とは限らない） -display(N(V[:, idx])) # 対応するVの列ベクトル（求めるもの） -s2 = np.array(sorted([N(x)**2 for x in s], reverse=True)) # 特異値の2乗（降順） -np.cumsum(s2) / sum(s2) # 累積寄与率（後述） -``` - -```{python} -_, s, tV = linalg.svd(An) # 特異値分解 -print(tV.T[:, 0]) # Vの第1列（求めるもの） -s2 = s**2 # 特異値の2乗 -np.cumsum(s2) / sum(s2) # 累積寄与率（後述） -``` - -```{python} -X = np.array([[1, 3, 6, 10], [7, 1, 6, 14]]).T - -from statsmodels.multivariate import pca # statsmodelsを使う場合 -pca1 = pca.PCA(X, standardize=False, normalize=False) -P = pca1.scores; print(P) # 主成分スコア -V = pca1.loadings; print(V) # 主成分 -print(V[:, 0]) # 第1主成分（求めるもの） -s2 = pca1.eigenvals # 特異値の2乗 -print(np.cumsum(s2) / sum(s2)) # 累積寄与率（後述） - -from sklearn import decomposition # Scikit-learnを使う場合 -pca2 = decomposition.PCA() -P = pca2.fit_transform(X); print(P) # 主成分スコア -tV = pca2.components_; print(tV.T) # 主成分 -print(tV[0]) # 第1主成分（求めるもの） -s2 = pca2.explained_variance_ratio_ # 特異値の2乗 -print(np.cumsum(s2) / sum(s2)) # 累積寄与率（後述） -``` - -# 20 特異値分解と擬似逆行列 - -```{python} -A = Matrix([[1, 0], [1, 1], [0, 1]]) -Ur, Sr, Vr = A.singular_value_decomposition() -Ur, Sr, Vr.T, A == simplify(Ur @ Sr @ Vr.T) -``` - -```{python} -An = np.array([[1, 0], [1, 1], [0, 1]]) -U, s, tV = linalg.svd(An) -S = np.zeros(An.shape); r = len(s); S[:r, :r] = np.diag(s) -U, S, tV, np.allclose(An, U @ S @ tV) -``` - -```{python} -from skimage import io -url = 'https://github.com/taroyabuki/comath/raw/main/images/boy.jpg' -A = io.imread(url, as_gray=True) # 画像の行列への変換 -U, s, tV = linalg.svd(A) # 特異値分解 -k = 52 -Ak = U[:, :k] @ np.diag(s[:k]) @ tV[:k, :] # 近似 -B = (Ak - np.min(Ak)) / (np.max(Ak) - np.min(Ak)) # 数値を0～1にする． -plt.figure(figsize=(10, 5)) -plt.subplot(1, 2, 1); plt.imshow(A, cmap='gray') -plt.subplot(1, 2, 2); plt.imshow(B, cmap='gray') -plt.show(); -``` - -```{python} -def exbasis(A): # tAの零空間（Aの列空間と直交）の基底を列ベクトルとする行列 - return Matrix.hstack(*GramSchmidt(A.T.nullspace(), orthonormal=True)) -def svd2(A, tol=10e-10): - m, n = A.shape; G = A.T @ A; eigs = G.eigenvects() # ① - vals = [re(e[0]) for e in eigs for _ in range(e[1])] # ②-1a - vecs = [re(v) for e in eigs for v in e[2]] # ②-1b - idx = np.argsort([-abs(N(x)) for x in vals]) # ②-2a - vals, vecs = [vals[i] for i in idx], [vecs[i] for i in idx] # ②-2b - s = [sqrt(x) for x in vals if abs(N(x)) > tol]; r = len(s) # ③ - if r != 0: - Sr = diag(*s) # ④ - Vr = Matrix.hstack(*GramSchmidt(vecs[:r], orthonormal=True)) # ⑤ - Ur = A @ Vr @ diag(*[1 / x for x in s]) # ⑥ - S = 0 * A; S[:r, :r] = Sr # ⑦ - V = Vr if n == r else Vr.row_join(exbasis(Vr)) # ⑧ - U = Ur if m == r else Ur.row_join(exbasis(Ur)) # ⑨ - else: - S = 0 * A; V = eye(n); U = eye(m) - Sr = Matrix([[0]]); Vr = V[:, 0]; Ur = U[: ,0] - return Ur, Sr, Vr, U, S, V - -A = Matrix([[1, 0], [1, 1], [0, 1]]); svd2(A) # 動作確認 -``` - -```{python} -def svd2n(A, tol=10e-10): - m, n = A.shape; G = A.T @ A # ① - vals, vecs = linalg.eig(G); vals, vecs = vals.real, vecs.real # ②-1 - idx = np.argsort(-abs(vals)) # ②-2a - vals, vecs = vals[idx], vecs[:, idx] # ②-2b - s = [np.sqrt(x) for x in vals if x > tol]; r = len(s) # ③ - Sr = np.diag(s) # ④ - V = linalg.qr(vecs[:, :r])[0]; Vr = V[:, :r] # ⑤, ⑧ - Ur = A @ Vr / s # ⑥ - S = 0.0 * A; S[:r, :r] = Sr # ⑦ - if (r != 0): U = np.hstack([Ur, linalg.null_space(Ur.T)]) # ⑨ - else: U = np.eye(m) - return Ur, Sr, Vr, U, S, V - -A = np.array([[1, 0], [1, 1], [0, 1]]); svd2n(A) # 動作確認 -``` - -```{python} -isOrtho = lambda A: np.allclose( - np.double(A.T) @ np.double(A), np.eye(A.shape[1])) -isSquare = lambda A: A.shape[0] == A.shape[1] -isDiagDesc = lambda A: ( - list(A.diagonal()) == sorted(abs(A.diagonal()), reverse=True)) - -A = Matrix([[1., 0], [1, 1], [0, 1]]) -Ur, Sr, Vr, U, S, V = svd2(A) # 特異値分解 -[isOrtho(Ur), isOrtho(Vr), isOrtho(U), isOrtho(V), # ① - isSquare(U), isSquare(V), # ② - isDiagDesc(Sr), isDiagDesc(S), # ③ - np.allclose(np.double(A), np.double(Ur @ Sr @ Vr.T)), # ④-1 - np.allclose(np.double(A), np.double(U @ S @ V.T))] # ④-2 -``` - -```{python} -isOrthon = lambda A: np.allclose(A.T @ A, np.eye(A.shape[1])) -isSquaren = lambda A: A.shape[0] == A.shape[1] -isDiagDescn = lambda A: all( - A.diagonal() == sorted(abs(A.diagonal()), reverse=True)) - -A = np.array([[1, 0], [1, 1], [0, 10]]) -Ur, Sr, Vr, U, S, V = svd2n(A) # 特異値分解 -[isOrthon(Ur), isOrthon(Vr), isOrthon(U), isOrthon(V), # ① - isSquaren(U), isSquaren(V), # ② - isDiagDescn(Sr), isDiagDescn(S), # ③ - np.allclose(np.double(A), np.double(Ur @ Sr @ Vr.T)), # ④-1 - np.allclose(np.double(A), np.double(U @ S @ V.T))] # ④-2 -``` - -```{python} -A = Matrix([[1, 0], [1, 1], [0, 1]]); A.pinv() -``` - -```{python} -A = np.array([[1, 0], [1, 1], [0, 1]]); linalg.pinv(A) -``` - -```{python} -A = Matrix([[1, 0], [1, 1], [0, 1]]); b = Matrix([2, 0, 2]) -A.pinv() @ b -``` - -```{python} -A = np.array([[1, 0], [1, 1], [0, 1]]); b = np.array([2, 0, 2]) -linalg.pinv(A) @ b -``` -