陰関数で表された曲線の描画
https://blog.goo.ne.jp/r-de-r/e/ff2399cb40043afca226208e8fefaec9
の例を Octave で描いてみた。
clf;
hold on;
for c = [-50, -30, -20, -15, -10, 0, 1, 2, 5, 10, 20, 30, 50]
ezplot(@(x, y) (x.^2 + y.^2 - 2 .* x).^3- 4*(x .* y).^2 - c, [-1.5, 3.5, -2.5, 2.5], n=1000)
end
clf;
hold on;
for c = -10:10
ezplot(@(x, y) x.^2 - x .* y + y.^2 - c, [-4, 4], n=1000)
end
a の 3 乗根（立方根）は，3乗すると a になる数だ。Julia には cbrt という関数がある。また，定義により a ^ (1/3) でもある。
どちらも同じか？
正の数の場合は同じだ。
julia> cbrt(8)
2.0
julia> 8^(1/3)
2.0
負の数の場合は違う。負の数を 1/3 乗することはできない。
julia> cbrt(-8)
-2.0
julia> (-8)^(1/3)
ERROR:
DomainError with -8.0:
Exponentiation yielding a complex result requires a complex argument.
Replace x^y with (x+0im)^y, Complex(x)^y, or similar.
指示に従って以下のようにする。
julia> (-8 + 0im)^(1/3)
1.0 + 1.7320508075688772im
得られる数を 3 乗すると（誤差の範囲内で）ちゃんと -8 になる。
julia> a = (-8.0 + 0.0im)^(1/3)
1.0 + 1.7320508075688772im
julia> a^3
-7.999999999999998 + 8.881784197001252e-16im
-------
他の言語ではどうなるか?
Octave
cbrt は Julia と同じ結果を返す。
>> cbrt(-8)
ans = -2
負の数を 1/3 乗しようとしても，Julia のようには文句を言わず，答えを返してくれる。
>> a = (-8)^(1/3)
a = 1.000000000000000 + 1.732050807568877i
>> a^3
ans = -7.999999999999995e+00 + 9.797174393178820e-16i
Python
numpy.cbrt も Julia と同じ結果を返す。
>>> import numpy as np
>>> np.cbrt(-8)
-2.0
負の数の 1/3 乗のときは，Octave と同じふるまいをする。
>>> a = (-8)**(1/3)
>>> a
(1.0000000000000002+1.7320508075688772j)
>>> a**3
(-8+3.1086244689504383e-15j)
R
R には cbrt 関数はない(SparkR というライブラリにあるようだが，R 4.2.1 には対応していないようだ）。
負の数の 1/3 乗のときは，Julia のように親切なアドバイスをくれるわけでもなく，ぶっきらぼうに NaN を返す。
> options(digits=16)
> (-8)^(1/3)
[1] NaN
虚数にすれば，答えを出してくれる。
> a = (-8 + 0i)^(1/3)
> a
[1] 1+1.732050807568877i
> a^3
[1] -7.999999999999996e+00+2e-15i
Jupyterlab 側で用意することは唯一，octave kernel のインストール
$ pip install octave_kernel
macOS 用の GNU Octave は
GNU Octave にも手を出してみる
https://blog.goo.ne.jp/r-de-r/e/a59385acd2252786199f2c45663ffaeb
に書いたように，
https://github.com/octave-app/octave-app/releases
で Octave 6.2.0 がインストールされる。
しかし，アプリケーションアイコンは Octave 自身ではなく Octave の IDE 環境（?）のアプレット applet を指している。
Octave 自身は，
/Applications/Octave-6.2.0.app/Contents/Resources/usr/Cellar/octave-octave-app@6.2.0/6.2.0/bin/octave-6.2.0
にある。これに対して /usr/local/bin の中にシンボリックリンク octave を張る。
$ ln -s /Applications/Octave-6.2.0.app/Contents/Resources/usr/Cellar/octave-octave-app@6.2.0/6.2.0/bin/octave-6.2.0 /usr/local/bin/octave
/usr/local/bin にはパスが通っているはず（通っていなければ通す）なので，コマンドライン（プロンプト）から
$ octave
と入力すれば，octave の REPL ウインドウが開く。
GNU Octave, version 6.2.0
Copyright (C) 2021 The Octave Project Developers.
This is free software; see the source code for copying conditions.
There is ABSOLUTELY NO WARRANTY; not even for MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. For details, type 'warranty'.
Octave was configured for "x86_64-apple-darwin18.7.0".
Additional information about Octave is available at https://www.octave.org.
Please contribute if you find this software useful.
For more information, visit https://www.octave.org/get-involved.html
Read https://www.octave.org/bugs.html to learn how to submit bug reports.
For information about changes from previous versions, type 'news'.
octave:1> 1+sqrt(8)
ans = 3.8284
octave:2> quit
これが終わったら Jupyterlab のカーネル Octave を選ぶだけだ。
>> x = [1 2 3 4 5; 6 7 8 9 10]
x =
1 2 3 4 5
6 7 8 9 10
のような場合，総和を求めるのに sum(sum(x)) としたのだが sum(x(:)) とも書けるそうだ。x(:) は行優先で一次元配列（ベクトル）に展開してくれる。
>> x(:)
ans =
1
6
2
7
3
8
4
9
5
10
>> sum(x(:))
ans = 55
>> sum(sum(x))
ans = 55
MATLAB は GNU Octave と部分的な互換性があるようだ。
そこで，MATLAB プログラムを動かすときには GNU Octave をインストールして使えばよいようだ。
https://octave.org/download.html
Windows では 7.1.0 がダウンロードできるようだ。
https://octave.org/download.html#ms-windows
macOS ではちょっと面倒なようだ。
https://github.com/octave-app/octave-app/releases
6.2.0 で我慢する。
====
GNU Octave の実行は他のアプリケーションと同じように，アイコンをクリックするだけ。
IDE 環境ウインドウが開く。コマンドウインドウになにか入力すれば，それが実行され，結果が表示される。
必要なパッケージの準備
例えば，以下のプログラムでは 「chi2cdf は statistics パッケージにあるよ！」 と言われる。
コマンドラインで statistics パッケージの最新バージョンを調べて（たとえば statistics-1.4.3.tar.gz のとき） wget でダウンロードする。
% wget http://downloads.sourceforge.net/project/octave/Octave%20Forge%20Packages/Individual%20Package%20Releases/statistics-1.4.3.tar.gz
パッケージを octave コマンドウインドウでインストールする（1回だけ）
>> pkg install -global statistics-1.4.3.tar.gz
パッケージを使用するとき（セッションごと）
>> pkg load statistics
使用するプログラム。
>> function [chisq, df, p] = ChiSquareTest(x)
[nrow , ncol] = size(x);
df = (nrow - 1) * (ncol - 1);
colsums = sum(x, 1);
rowsums = sum(x, 2);
n = sum(rowsums);
expectations = (rowsums * colsums) / n;
Zscores = (x .- expectations) .^ 2 ./ expectations;
chisq = sum(sum(Zscores));
p = 1 - chi2cdf(chisq, df);
end
>> x = [55 22 16 7; 40 32 24 4];
[chisq, df, p] = ChiSquareTest(x)
chisq = 6.6385
df = 3
p = 0.084359
>> [chisq, df, p] = ChiSquareTest([10 5; 6 11])
chisq = 3.1373
df = 1
p = 0.076523
* プログラム中に日本語があると変なことが起きることがあるようだ（プログラムをコピペできない）。
* MATLAB のすべてが octave で利用できるというものではない。例えば，icdf('chi2', x, df) はないので，1 - chi2cdf(x, df) とするなど，エラー回避が必要になることもある。
