相関を求めて

久し振りにgoshが新版になった。 Gauche A Scheme Implementation 何時ものように下記にて更新。

% gauche-config --reconfigure | sh
% make
% make install

SRFIが貪欲に出てくるのを見ると、schemeを愛してますって熱狂が伝わってき ますよ。

前回、コスタス行列を調べた時、自己相関なんてのが出てきた。何だっけ？ たまには、ChatGPTに質問してみるか。100の説明より、一つのコードだな。世 間で一番使われてる奴でいいだろう。

Pythonで配列 `bldat` の自己相関を求めるには、NumPyというライブラリを使用すると便利です。以下に、NumPyを使用した自己相関を求めるPythonコードの例を示します。

import numpy as np

# サンプルのデータ (bldat) を作成
# ここでは仮のデータを使用します。実際のデータに置き換えてください。
bldat = np.random.rand(100)

# 自己相関を計算
autocorrelation = np.correlate(bldat, bldat, mode='full')

# 自己相関の結果を表示
print("自己相関の結果:")
print(autocorrelation)

# 自己相関のプロット
import matplotlib.pyplot as plt

# ゼロラグの位置を取得
zero_lag_index = len(bldat) - 1

# グラフをプロット
plt.figure(figsize=(10, 5))
plt.plot(range(-zero_lag_index, zero_lag_index + 1), autocorrelation)
plt.xlabel("ラグ")
plt.ylabel("自己相関")
plt.title("自己相関のプロット")
plt.grid()
plt.show()

このコードでは、まずNumPyを使って仮のデータを生成し、その後、`np.correlate()` 関数を使用して自己相関を計算しています。自己相関の結果を表示し、さらにMatplotlibを使用して自己相関のプロットを作成しています。

実際のデータを `bldat` に代入してコードを実行することで、そのデータの自己相関を求めることができます。

numpyを使って説明を省略してるな。これが世間の風潮ですか。プロットもサー ビスしときますんで、今後もご贔屓にして下さい、とな。

ならば、このサービスを受容できる様にしとくのが礼儀ってものだろう。 OpenBSDにはvenv環境が無いので、それを作ってからmatplotlibをコンパイル しようと思った。忍者だとか色々必要だと言われたので、諦めた。素直に py3-matplotlibを入れたよ。numpyもつられて入ってきた。まあ、愛用してた gnuplotよりは、豪華な環境になった。

上のコードをそのまま走らせると、盛んにワーニングが出てくる。CJKのカナ・ フォントが入っていないぞと。そして久し振りに豆腐とご対面したよ。グラフ の説明ぐらいなら英語で十分なんだけどな。一般人向けなコードを生成してる んだな。

オイラーみたいに臍曲りは、numpyにサラリと移譲しちゃうのは嫌い。 pure-pythonで提示を要求。

NumPyを使用せずに、Python標準ライブラリのみを使用して、配列 `bldat` の自己相関を計算する方法を示します。以下がそのコードです：

# サンプルのデータ (bldat) を作成
# ここでは仮のデータを使用します。実際のデータに置き換えてください。
bldat = [0.2, 0.5, 0.8, 0.6, 0.3, 0.1, 0.7, 0.4, 0.9, 0.2]

# 自己相関を計算する関数を定義
def autocorrelation(data, lag):
    n = len(data)
    mean = sum(data) / n
    c0 = sum((x - mean) ** 2 for x in data) / n

    num = sum((data[i] - mean) * (data[i + lag] - mean) for i in range(n - lag))
    den = sum((x - mean) ** 2 for x in data)

    return num / den

# 自己相関の結果を計算
max_lag = len(bldat) - 1
autocorrelation_values = [autocorrelation(bldat, lag) for lag in range(max_lag + 1)]

# 自己相関の結果を表示
print("自己相関の結果:")
print(autocorrelation_values)

このコードは、自己相関を計算するための関数 `autocorrelation` を定義し、それを使用して自己相関を計算しています。データの平均値と分散を使用して自己相関を計算し、計算結果をリスト autocorrelation_values に格納しています。最後に、自己相関の結果を表示します。

実際のデータを `bldat` に代入してコードを実行することで、そのデータの自己相関を求めることができます。

どうもChatGPTって、内包表記が好きみたいだな。その方がコードがhaskellチッ クになって、オイラーも好きであります。

>>> import numpy as np
>>> help(np.correlate)

で、豊富なコメントを参照できるんだな。忘れておったわい。

それから、少し高度な説明として、自己相関関数とパワースペクトラム こう いうのが見付かった。

import statistics

x = [1, 2, 3, 4, 5]
y = [5, 4, 3, 2, 1]

correlation = statistics.correlation(x, y)
print("corr:", correlation)

def correlation(x, y, /):
    n = len(x)
    xbar = fsum(x) / n
    ybar = fsum(y) / n
    sxy = fsum((xi - xbar) * (yi - ybar) for xi, yi in zip(x, y))
    sxx = fsum((xi - xbar) ** 2.0 for xi in x)
    syy = fsum((yi - ybar) ** 2.0 for yi in y)
    return sxy / sqrt(sxx * syy)

セカンド・オピニオンでRにも聞いてみる。語源はChatGPTが何気に使っていた 用語。とっかかりが必要だ。

> ??correlation
  :
nlme::corSymm           General Correlation Structure
nlme::VarCorr           Extract variance and correlation components
stats::acf              Auto- and Cross- Covariance and -Correlation
                        Function Estimation
stats::cancor           Canonical Correlations
stats::cor.test         Test for Association/Correlation Between Paired
                        Samples
  Concepts: Kendall correlation coefficient, Pearson correlation
    coefficient, Spearman correlation coefficient
stats::plot.acf         Plot Autocovariance and Autocorrelation
                        Functions
stats::var              Correlation, Variance and Covariance (Matrices)

色々と案内してくれてるけど、目指すは、stats::acfかな。

> help(acf)
Description:

     The function ‘acf’ computes (and by default plots) estimates of
     the autocovariance or autocorrelation function.  Function ‘pacf’
     is the function used for the partial autocorrelations.  Function
     ‘ccf’ computes the cross-correlation or cross-covariance of two
     univariate series.
     :

やっぱりソース祭を開催。~/src/R-4.3.1/src/library/stats/R/acf.R

acf <-
    function (x, lag.max = NULL,
              type = c("correlation", "covariance", "partial"),
              plot = TRUE, na.action = na.fail, demean = TRUE, ...)
{
    type <- match.arg(type)
    if(type == "partial") {
        m <- match.call()
    :

気合が入ってるんで、福読本が欲しいなあ。起点は、奥村先生の 所。統計・データ解析

アイスクリームの売れ方

相関（と回帰）

概念的には、自己相関関数入門 かな。

で検索すると出てくるのは、ほとんど式の羅列の大学の講義資料ばかり。諦め ないで探していたら、ディジタルデータ処理資料 に出会った。

ソースはMATLABだ。互換のoctaveを使えばいいのか。

ウィキベディ アもお忘れなく。自己相関 そして、C言語の付録付きの 自己相関関数

そもそもの発端になった記事をトレースしてみる。

import matplotlib.pyplot as plt

gol = [0,1,4,13,28,33,47,54,64,70,72]  ## Golomb ruler
idx = [0,2,9,10,12,22,32,42,52,62,72]  ## sync 10 cycle
bldat = [0 for x in range(75)]
for x in gol:       # or idx
    bldat[x] = 1

def acf(dat, lag):
    n = len(dat)
    mean = sum(dat) / n

    num = sum((dat[i] - mean) * (dat[i + lag] - mean) for i in range(n - lag))
    den = sum((x - mean) ** 2 for x in dat)

    return num / den

max_lag = len(bldat) - 1
acf_values = [acf(bldat, lag) for lag in range(max_lag + 1)]

plt.figure(figsize=(10, 5))
plt.stem(range( max_lag + 1), acf_values)
plt.grid()
plt.show()

ゴロム数列(gol)と、ほぼ10置きのインパルス(idx)を(コード上で切り替えて)発生させ、 それの自己相関を表示。ゴロムの方は、どんな組み合わせと取っても、周期が一致す る事はない。従ってlag=0、波形が重なった所でのみ、相関係数が1になり、そ れ以外は(ほぼ)0になる。

一方のインパルスは、ほぼ10おきにインパルスが出てくるので、波形を10シフトしても、そ こでの相関値は高くなる。やっと意味が理解できましたよ。

import matplotlib.pyplot as plt
import math
import random

a = 0.017453292519943295 # pi/180 for degree to radian
bldat = [math.sin(3 * a * t) + random.uniform(-1.0, 1.0)
         for t in range(1024)]

def acf(dat, lag):
    n = len(dat)
    mean = sum(dat) / n
    num = sum((dat[i] - mean) * (dat[i + lag] - mean) for i in range(n - lag))
    den = sum((x - mean) ** 2 for x in dat)
    return num / den

max_lag = len(bldat) - 1
acf_values = [acf(bldat, lag) for lag in range(max_lag + 1)]

fig = plt.figure(figsize=(10, 7))
up = fig.add_subplot(2, 1, 1)      # rows cols loc
dn = fig.add_subplot(2, 1, 2)
up.plot(range( len(bldat) ), bldat)
up.grid()
dn.plot(range(max_lag + 1 ), acf_values)
dn.grid()
plt.show()
fig.savefig("test.svg")

こちらは、3Hzのサイン波に±1のノイズを重畳させ、その元波形と自己相関を 取ってみたものだ。120度毎に相関値が高くなっている事を確認できる。

matplotlibで複数グラフの描画方法とグラフの保存の例題にもなってる。こう してライブラリィーの利用方法を身に付けていくんだなあ。例題の為の例題よ りかは、よっぽど為になるぞ。

import matplotlib.pyplot as plt
import numpy as np

a = 0.017453292519943295 # pi/180 for degree to radian
bldat = np.array([np.sin(3 * a * t) + np.random.uniform(-1.0, 1.0)
                  for t in range(4096)])

autocorrelation = np.correlate(bldat, bldat, mode='full')
zero_lag_index = len(bldat) - 1

fig = plt.figure(figsize=(10, 7))
up = fig.add_subplot(2, 1, 1)      # rows cols loc
dn = fig.add_subplot(2, 1, 2)
up.plot(range( len(bldat) ), bldat)
up.grid()
dn.plot(range(-zero_lag_index, zero_lag_index + 1), autocorrelation)
dn.grid()
#plt.show()
#fig.savefig("test.svg")

vbox$ time python3 wave.py
    0m35.62s real     0m34.60s user     0m00.93s system
vbox$ time python3 wave-np.py
    0m03.60s real     0m02.49s user     0m00.89s system

import numpy as np

csv = np.loadtxt('am.csv',  delimiter=',', dtype='int32', usecols=(0, 1, 2))
hi = csv[:,1]

emacsに付属の高性能RPN式の電卓。 Vectors as Lists (Tutorial) を見ていたら、外部のデータとやり取り出来る ようだ。これで、ちょいとした演算が可能になるな。

For moving data into and out of Calc:
  G  calc-grab-region.  Grab the region defined by mark and point into Calc.
  R  calc-grab-rectangle.  Grab the rectangle defined by mark, point into Calc.
  :  calc-grab-sum-down.  Grab a rectangle and sum the columns.
  _  calc-grab-sum-across.  Grab a rectangle and sum the rows.
  Y  calc-copy-to-buffer.  Copy a value from the stack into the editing buffer.

123  0.22
234  0.33
321  1.55

こんなデータがbufferに有る時、これらをレクタングル指定してから、C-x * r すれば、calcが開いて、データがスタックに取り込まれる。

1:  [ [ 123, 0.22 ]                                 |
      [ 234, 0.33 ]                                 |
      [ 321, 1.55 ] ]

これは便利だな。後は山程あるコマンドをどう使いこなすかの勝負だ。

   The ‘H u C’ (‘calc-vector-correlation’) [‘vcorr’] command computes
the linear correlation coefficient of two vectors.  This is defined by
the covariance of the vectors divided by the product of their standard
deviations.

23090105,140,80,60
  :
23093004,143,76,56

v .
1:  [ [140, 80],                                    |grab [[119], [131], [118],$
      [145, 78],                                    | trn [[119, 131, 118, 131,$
      [132, 74],                                    |grab [[117], [139], [145],$
      ...,                                          | trn [[117, 139, 145, 130,$
      [143, 76] ] 

v t
1:  [ [140, 145, 132, ..., 143],                    |corr 0.386011055468
      [80, 78, 74, ..., 76] ]                       |grab [[140, 80], [145, 78]$
v u
2:  [140, 145, 132, ..., 143]                       | trn [[140, 145, 132, 147,$
1:  [80, 78, 74, ..., 76]                           |unpk [140, 145, 132, 147, $
H u C
1:  0.100366294405                                  |unpk [140, 145, 132, 147, $

The ‘u M’ (‘calc-vector-mean’)
The ‘H u M’ (‘calc-vector-median’)
The ‘u S’ (‘calc-vector-sdev’)
The ‘u X’ (‘calc-vector-max’)
The ‘u N’ (‘calc-vector-min’)
The ‘u #’ (‘calc-vector-count’)
The ‘u +’ (‘calc-vector-sum’)

vbox$ cd /usr/local/share/emacs/29.1/lisp/calc
vbox$ zgrep calc-vector-correlation *el.gz
calc-ext.el.gz: ("calc-stat" calc-vector-correlation calc-vector-count
calc-stat.el.gz:(defun calc-vector-correlation (arg)

 ("calc-stat" calc-vector-correlation calc-vector-count
calc-vector-covariance calc-vector-geometric-mean
calc-vector-harmonic-mean calc-vector-max calc-vector-mean
calc-vector-mean-error calc-vector-median calc-vector-min
calc-vector-pop-covariance calc-vector-pop-sdev
calc-vector-pop-variance calc-vector-product calc-vector-rms
calc-vector-sdev calc-vector-sum calc-vector-variance)

(defun calc-vector-correlation (arg)
  (interactive "P")
  (calc-hyperbolic-func)
  (calc-vector-covariance arg))

(defun calc-vector-covariance (arg)
  (interactive "P")
  (calc-slow-wrapper
   (let ((n (if (eq arg 1) 1 2)))
     (if (calc-is-hyperbolic)
         (calc-enter-result n "corr" (cons 'calcFunc-vcorr
                                           (calc-top-list-n n)))
       (if (calc-is-inverse)
           (calc-enter-result n "pcov" (cons 'calcFunc-vpcov
                                             (calc-top-list-n n)))
         (calc-enter-result n "cov" (cons 'calcFunc-vcov
                                          (calc-top-list-n n))))))))

(defun calcFunc-vpcov (vec1 &optional vec2)
  (math-covariance (list vec1) (list vec2) t 1))

(defun math-covariance (vec1 vec2 pop mode)
  (or (car vec2) (= mode 0)
      (progn
        (if (and (eq (car-safe (car vec1)) 'var)
                 (eq (car-safe (calc-var-value (nth 2 (car vec1)))) 'vec))
            (setq vec1 (symbol-value (nth 2 (car vec1))))
          (setq vec1 (car vec1)))
           :

長年、血圧降下薬を朝に服用していたんだけど、その薬効が低下してきて、 時々、高記録が出るようになった。で、2023-10-13 から、夜用も追加しても らった。これで、朝まで安心。

[sakae@deb blood]$ newstats 231014
          mode(#)   median     mean   N from   to
hi@am      116(3)    119.0   121.27 18 231014 231031
hi@pm      124(3)    127.0   128.38 18 231014 231031
lo@am       71(4)     67.0    66.61 18
lo@pm       68(4)     68.0    68.05 18
[sakae@deb blood]$ newstats 2307
          mode(#)   median     mean   N from   to
hi@am      138(8)    140.0   139.71 100 230701 231008
hi@pm      131(7)    131.0   129.63 100 230701 231008
lo@am       74(13)     74.0     73.2 100
lo@pm       68(14)     67.0    66.41 100

---

This year's Index

Home

---