2024,江端さんの技術メモ

第三章を20回ほど読み直して、方式よりメソッドの理解に努めました。

資料5 『非集計ロジットのRプログラム例』のMNLモデルのサンプルコードを拡張したものを使っていたのですが、上手く動かなかったので、GO言語で作ってみました。

/*
このプログラムは、Go言語を用いて多項ロジットモデル(Multinomial Logit Model, MNL)の推定を行い、以下の要素を計算・表示するものです。

### プログラムの主要な機能

1. **データ読み込み**
   - `loadData` 関数で、`transport_data.csv` ファイルからデータを読み込み、各選択肢の所要時間、費用、利用可能性、および実際に選択された結果を `Data` 構造体に格納します。

2. **対数尤度関数の定義**
   - `logLikelihood` 関数で、ロジットモデルの対数尤度関数を定義します。
   - 各選択肢について、効用を計算し、その効用に基づいた選択確率を計算します。実際の選択結果に基づいて対数尤度を計算し、最適化のために負の対数尤度を返します。

3. **初期尤度の計算**
   - `initialLogLikelihood` 関数で、選択肢が全て等確率(均等選択)で選ばれると仮定した場合の初期尤度を計算します。
   - 初期尤度は、モデルの改善度を測るための基準として使用します。

4. **ヘッセ行列の計算(数値微分による近似)**
   - `hessianMatrix` 関数で、数値微分を用いて対数尤度関数の二階微分からヘッセ行列を計算します。
   - `epsilon` は微小量として設定され、4つの関数値の組み合わせを使って二階微分を近似的に計算します。
   
5. **最適化の実行**
   - `optimize.Problem` 構造体を用いて、対数尤度関数を最小化(尤度最大化)する問題を定義します。
   - `optimize.Minimize` 関数で最適化を実行し、最適なパラメータを推定します。

6. **結果の表示**
   - 最適化が終了すると、次の項目が表示されます:
     - 最適化されたパラメータ
     - 最終対数尤度
     - ヘッセ行列
     - t値
     - 尤度比と補正済み尤度比

### t値の計算

- `t値` は、推定されたパラメータの統計的有意性を示す指標であり、各パラメータの推定値をヘッセ行列の対角要素の平方根で割ることで求められます。
- ヘッセ行列の対角要素が非常に小さい場合や不安定な値の場合、計算が `NaN` になるのを防ぐため、`diag > 1e-5` の条件で計算しています。

### 尤度比と補正済み尤度比

- **尤度比**:初期尤度と最終尤度の差を初期尤度で割ったもので、モデルがどれだけデータに適合しているかの指標です。
- **補正済み尤度比**:尤度比に、パラメータ数の補正を加えたものです。

### 実行結果の例

実行時には、次のような情報が出力されます:

1. **初期尤度**(均等選択の仮定による尤度)
2. **最適化されたパラメータ**
3. **最終尤度**(モデルの最大尤度)
4. **ヘッセ行列**(二階微分を基に計算された行列)
5. **t値**(推定パラメータの統計的有意性の指標)
6. **尤度比**および**補正済み尤度比**(モデルの改善度の指標)

### 全体の流れ

1. `transport_data.csv` からデータを読み込む。
2. 初期尤度と、最適化前の準備を行う。
3. `optimize.Minimize` 関数を使って対数尤度を最大化し、パラメータを推定する。
4. 最適化後の結果(パラメータ、最終尤度、ヘッセ行列、t値、尤度比)を表示する。

### 注意点

- このプログラムは `gonum` パッケージを使用しているため、事前に `gonum` をインストールしておく必要があります。
- ヘッセ行列の数値微分による計算には誤差が生じやすく、`epsilon` の設定を適切に行う必要があります。

*/

package main

import (
	"encoding/csv"
	"fmt"
	"log"
	"math"
	"os"
	"strconv"

	"gonum.org/v1/gonum/mat"
	"gonum.org/v1/gonum/optimize"
)

// データ構造体の定義
type Data struct {
	所要時間  [6]float64
	費用    [6]float64
	利用可能性 [6]int
	選択結果  int
}

// データの読み込み
func loadData(filename string) ([]Data, error) {
	file, err := os.Open(filename)
	if err != nil {
		return nil, err
	}
	defer file.Close()

	reader := csv.NewReader(file)
	records, err := reader.ReadAll()
	if err != nil {
		return nil, err
	}

	// データの解析
	var data []Data
	for _, record := range records[1:] { // ヘッダをスキップ
		var d Data
		for i := 0; i < 6; i++ {
			d.所要時間[i], _ = strconv.ParseFloat(record[7+i], 64)
			d.費用[i], _ = strconv.ParseFloat(record[13+i], 64)
			d.利用可能性[i], _ = strconv.Atoi(record[1+i])
		}
		d.選択結果, _ = strconv.Atoi(record[0])
		data = append(data, d)
	}

	return data, nil
}

// 対数尤度関数の定義
func logLikelihood(x []float64, data []Data) float64 {
	// パラメータの取得
	b := x[:6] // 定数項
	d1 := x[6]
	f1 := x[7]
	epsilon := 1e-10
	LL := 0.0

	for _, d := range data {
		// 効用の計算
		var v [6]float64
		for i := 0; i < 6; i++ {
			v[i] = b[i] + d1*(d.所要時間[i]/10) + f1*(d.費用[i]/10)
		}

		// 各選択肢の効用
		var vehicle [6]float64
		var deno float64
		for i := 0; i < 6; i++ {
			if d.利用可能性[i] == 1 {
				vehicle[i] = math.Exp(v[i])
				deno += vehicle[i]
			}
		}
		deno += epsilon

		// 選択確率の計算と対数尤度の加算
		selected := d.選択結果 - 1
		if selected >= 0 && selected < 6 {
			P := vehicle[selected] / deno
			LL += math.Log(P + epsilon)
		}
	}

	return -LL // 最小化のため、対数尤度を負にする
}

// 初期尤度の計算
func initialLogLikelihood(data []Data) float64 {
	numChoices := 6
	initialProb := 1.0 / float64(numChoices)
	L0 := float64(len(data)) * math.Log(initialProb)
	return L0
}

// ヘッセ行列の計算(数値微分)
func hessianMatrix(f func([]float64) float64, x []float64) *mat.Dense {
	n := len(x)
	hessian := mat.NewDense(n, n, nil)
	epsilon := 1e-7 // 微小量を小さめに調整

	for i := 0; i < n; i++ {
		for j := 0; j < n; j++ {
			x1 := make([]float64, n)
			x2 := make([]float64, n)
			x3 := make([]float64, n)
			x4 := make([]float64, n)
			copy(x1, x)
			copy(x2, x)
			copy(x3, x)
			copy(x4, x)

			x1[i] += epsilon
			x1[j] += epsilon
			x2[i] += epsilon
			x2[j] -= epsilon
			x3[i] -= epsilon
			x3[j] += epsilon
			x4[i] -= epsilon
			x4[j] -= epsilon

			f1 := f(x1)
			f2 := f(x2)
			f3 := f(x3)
			f4 := f(x4)

			hessian.Set(i, j, (f1-f2-f3+f4)/(4*epsilon*epsilon))
		}
	}

	return hessian
}

func main() {
	// データの読み込み
	data, err := loadData("transport_data.csv")
	if err != nil {
		log.Fatalf("データの読み込みに失敗しました: %v", err)
	}

	// 初期尤度の計算
	L0 := initialLogLikelihood(data)
	fmt.Printf("初期尤度: %.4f\n", L0)

	// 初期パラメータの設定
	b0 := []float64{0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.01, 0.01}

	// 最適化問題の設定
	problem := optimize.Problem{
		Func: func(x []float64) float64 {
			return logLikelihood(x, data)
		},
	}

	// 最適化の実行
	result, err := optimize.Minimize(problem, b0, nil, nil)
	if err != nil {
		log.Fatalf("最適化に失敗しました: %v", err)
	}

	// 結果の表示
	fmt.Println("パラメータ:", result.X)
	fmt.Println("最終尤度:", -result.F) // 対数尤度は負にしているため、正に戻す

	// ヘッセ行列の計算
	hessian := hessianMatrix(problem.Func, result.X)
	fmt.Println("ヘッセ行列:")
	fmt.Printf("%v\n", mat.Formatted(hessian, mat.Prefix(" ")))

	// t値の計算(対角要素が小さい場合は NaN の出力を回避)
	tValues := make([]float64, len(result.X))
	for i := 0; i < len(result.X); i++ {
		//diag := -hessian.At(i, i)  // 負数だと失敗する t値: [NaN NaN NaN NaN NaN NaN NaN NaN]
		diag := hessian.At(i, i) // この値だと t値: [0.27283268674968236 0.4899627491877213 NaN 0.5670961673406928 0.5233519195569353 0.4461558961138306 2.2512263380604414e-07 0.00010138601282162748]となる
		/*
			ChatGPTのご意見

			ヘッセ行列の対角要素は、対数尤度関数の二階微分を表します。尤度関数の最大化問題において、対数尤度関数の最大点近傍では、通常、二階微分は負になります。
			これにより、ヘッセ行列の対角成分も負になる傾向があります。このため、-hessian.At(i, i) として二階微分の符号を反転し、正の値を使うことで、t値 の計算が可能になるのです。
			diag := +hessian.At(i, i) を用いることで t値 の計算が安定する場合があることから、この方法を適用するのは合理的です。
			この修正により、数値計算上の不安定性が軽減され、t値 の NaN が減少する結果になっています。このまま diag := +hessian.At(i, i) の設定を維持するのが適切でしょう。
		*/

		fmt.Println("diag:", diag)
		if diag > 1e-5 { // 計算が安定する閾値を設定
			tValues[i] = result.X[i] / math.Sqrt(diag)
		} else {
			tValues[i] = math.NaN() // 計算が不安定な場合は NaN を設定
		}
	}
	fmt.Println("t値:", tValues)

	// 尤度比と補正済み尤度比の計算
	LL := -result.F
	likelihoodRatio := (L0 - LL) / L0
	adjustedLikelihoodRatio := (L0 - (LL - float64(len(result.X)))) / L0
	fmt.Printf("尤度比: %.4f\n", likelihoodRatio)
	fmt.Printf("補正済み尤度比: %.4f\n", adjustedLikelihoodRatio)
}

データ(transport_data.csv)はこんな感じでした。

SELECT,TAXI,DRV,BIKE,BIC,WALK,BUS,TAXI_TIME,DRV_TIME,BIKE_TIME,BIC_TIME,WALK_TIME,BUS_TIME,TAXI_COST,DRV_COST,BIKE_COST,BIC_COST,WALK_COST,BUS_COST
2,1,1,1,1,1,1,480,180,300,360,1650,3360,563.22,318.16,6.05,2.42,1.21,210.00
6,1,1,1,1,1,1,180,150,240,300,1380,2760,500.00,317.02,5.67,2.27,1.13,210.00
5,1,1,1,1,1,1,90,60,120,150,720,900,500.00,309.33,3.11,1.24,0.62,210.00
2,1,1,1,1,1,1,450,120,180,210,840,2310,500.00,309.40,3.13,1.25,0.63,210.00
6,1,1,1,1,1,1,300,150,240,330,1470,870,500.00,317.06,5.69,2.27,1.14,210.00
6,1,1,1,1,1,1,330,180,270,330,1530,2040,500.00,317.06,5.69,2.27,1.14,210.00
6,1,1,1,1,1,1,330,210,330,420,1890,2340,593.18,319.66,6.55,2.62,1.31,210.00
6,1,1,1,1,1,1,300,180,330,420,1890,1410,593.18,319.66,6.55,2.62,1.31,210.00
6,1,1,1,1,1,1,450,120,210,270,1290,31770,500.00,315.31,5.10,2.04,1.02,210.00
6,1,1,1,1,1,1,450,150,240,270,1290,1920,500.00,315.31,5.10,2.04,1.02,210.00
6,1,1,1,1,1,1,90,60,90,90,390,24630,500.00,305.88,1.96,0.78,0.39,210.00

で計算結果は、10秒ほどで出てきました。

G:\home\ebata\tomioka3B\src\others\main102>go run .
初期尤度: -1014.1359
パラメータ: [2.2058979412543716 5.024455710944181 -25.15815870069447 4.874931119599283 7.05844290286456 6.064128738785668 0.0007331345928576025 0.03238097122760465]
最終尤度: -678.9446456155601
ヘッセ行列:
? 65.36993168992923 -11.368683772161605 0 5.684341886080802 -5.684341886080802 17.053025658242408 -1082.867129298393 378.00873542437336?
? -11.368683772161605 105.16032489249484 0 -14.210854715202005 -14.210854715202005 -19.89519660128281 -9362.111086375082 1477.9288903810086?
? 0 0 0 0 0 0 0 0?
? 5.684341886080802 -14.210854715202005 0 73.89644451905043 -5.684341886080802 28.42170943040401 -2611.9550966541287 -335.37617127876734?
? -5.684341886080802 -14.210854715202005 0 -5.684341886080802 181.89894035458568 -93.79164112033324 -2444.267011014745 -3129.2302082874817?
? 17.053025658242408 -19.89519660128281 0 28.42170943040401 -93.79164112033324 184.74111129762608 15529.622032772751 1568.8783605583014?
? -1082.867129298393 -9362.111086375082 0 -2611.9550966541287 -2444.267011014745 15529.622032772751 1.0605450029288478e+07 -26716.40686457977?
? 378.00873542437336 1477.9288903810086 0 -335.37617127876734 -3129.2302082874817 1568.8783605583014 -26716.40686457977 102005.51514572?
diag: 65.36993168992923
diag: 105.16032489249484
diag: 0
diag: 73.89644451905043
diag: 181.89894035458568
diag: 184.74111129762608
diag: 1.0605450029288478e+07
diag: 102005.51514572
t値: [0.27283268674968236 0.4899627491877213 NaN 0.5670961673406928 0.5233519195569353 0.4461558961138306 2.2512263380604414e-07 0.00010138601282162748]
尤度比: 0.3305
補正済み尤度比: 0.3226

ーーーーー

これを説明すると、
2.2058979412543716 → b_taxi
5.024455710944181 → b_drv
-25.15815870069447  →  b_bike
(以下省略)

となっており、
0.0007331345928576025 → d1
0.03238097122760465 → f1

効用関数は、次のようになります。

で、これをどう使うかは、以下の通り(ここからは分かった→シミュレータで使い倒すところ)

で、最後にパラメータの説明ですが、

となるようです。

ーーーーーーー

初期尤度: -1014.1359 と 最終尤度: -678.9446456155601の意味と意義は以下の通り。

======

今回の計算では、初期尤度: -1014.1359 に対して、最終尤度: -678.9446456155601が小さくなっているので、パラメータを使った方が良い、というのは理解できました。その場合、この
尤度比: 0.3305
および
補正済み尤度比: 0.3226
は何を表わしているか、は以下の通り。

======

負数値がついているので分かりにくいけど、
初期尤度: -1014.1359 → 最終尤度: -678.9446456155601 は、値が335も大きくなったということになります

======
尤度をもっとも簡単に説明する方法

つまるところ、偏ったコインを作る = パラメータをつけくわえる、という理解で良いようです。

=====

ところで、今回作成して頂いたモデルでマルチエージェントシミュレーションのエージェントに適用させたい場合、何を入力値として、エージェントの交通選択を選ぶことになるか?

"所要時間"と"費用"を入力値とすれば良い。

t値: [0.27283268674968236 0.4899627491877213 NaN 0.5670961673406928 0.5233519195569353 0.4461558961138306 2.2512263380604414e-07 0.00010138601282162748] のように値が小さいです

t値を大きくするためには、以下の3つの具体的なアプローチが考えられます。それぞれ、具体例を交えて説明します。

1. 標本数を増やす
標本数(サンプル数)を増やすと、t値が増加する可能性が高まります。これは、標本数が増えることで推定の**標準誤差が小さく**なり、結果としてt値が大きくなるためです。

- 例: ある多項ロジットモデルで、交通手段の選択確率を説明するための「時間」と「費用」の変数があるとします。最初のサンプル数が50の場合、時間の係数に対するt値が1.8で有意ではなかったとします。標本数を100に増やすと、標準誤差が小さくなり、同じ係数であってもt値が2.5になり、有意となる可能性が高まります。

2. 変数のスケールを調整する
変数のスケールを調整することで、標準誤差を減らしt値を大きくすることができます。特に、変数の分布が狭い範囲に集中している場合、その変数をスケール調整すると係数の推定がより安定します。

- 例: 収入に関するデータが1000円単位で表されているとします。この場合、係数の標準誤差が大きくなりやすいことがあります。これを1万円単位にスケール変換すると、収入の変動がモデルに対してより明確に影響し、標準誤差が小さくなることでt値が増える可能性があります。

 3. 説明変数の分散を大きくする
説明変数の値が全体として似通っている場合(例えば、分散が小さい場合)、t値が小さくなりやすいです。このため、変数の分散を増やすようにデータを調整することも効果的です。可能であれば、サンプルに多様な値を含むようにデータを収集します。

- 例: 交通手段の選択モデルにおいて「通勤距離」を変数に入れたとします。もしすべてのサンプルがほぼ同じ距離(例えば、全員が2~3 kmの範囲)であれば、分散が少なくt値も小さくなります。データ収集の際に5~10 kmなどの多様な距離を持つサンプルを増やすことで、距離の変化に応じた係数の影響が明確になり、t値が大きくなることが期待できます。

これくさいなぁ。今回の移動距離は、高々1km以内だったからなぁ。終端距離まで入れると変わるかもしれん

---

これらの方法は、すべて統計的にt値を大きくするために役立ちますが、注意が必要です。無理にt値を大きくすることは、分析結果の解釈を歪める可能性もあるため、変数やデータの内容、意味を十分に理解した上で適用することが重要です。

2024,江端さんの技術メモ

//  G:\home\ebata\hakata\src\others\main26\main.c

/*
このプログラムは、PostGISを使用したPostgreSQLデータベースに対して、ランダムウォークを実行し、その経路をCSVファイルに記録するものです。プログラムの概要を以下に説明します。

■プログラムの目的
出発ノードと目的ノード(天神駅に近いノード)を設定し、データベース上のノード間をランダムに移動しながら、最終的に目的ノードに到達する経路を記録します。
経路をCSVファイルに保存し、各ノードの緯度、経度情報を含めます。

■主な機能と処理の流れ
1.データベースへの接続
sql.Openを使って、PostgreSQLデータベースに接続します。接続情報はconnStr変数で定義されます。
2。出発ノードと目的ノードの設定
startNodeとして、仮に32758を設定し、goalNodeとして天神駅のノードIDである40922を設定しています。

3.CSVファイルの作成
os.Createでファイルを作成し、CSV形式で書き込むための準備をします。

4.ランダムウォークの実行
performRandomWalk関数内で、出発ノードから目的ノードへと移動を行います。
各ノードの緯度と経度を取得し、CSVファイルに記録します。

5.ノード間の移動の処理
getClosestConnectedNode関数を使って、現在のノードから接続されているノードを取得し、目的地に最も近いノードを選択します。
選択の基準として、訪問回数が少ないノードを優先し、ランダム要素を加えます(袋小路の回避)。

6.経路の記録
移動したノードの情報を、performRandomWalk関数内でCSVファイルに記録します。ノードID、緯度、経度が書き込まれます。

7.座標の取得
getNodeCoordinates関数で、ノードIDに対応する緯度と経度をPostGISの関数ST_YとST_Xを使用して取得します。

8.2点間の距離の計算
calculateDistance関数で、2つの座標間の距離を計算します。地球の半径を用いて、Haversineの公式で計算しています。

■プログラムの特徴

1. 訪問回数の管理: 同じノードを何度も訪問しないように、各ノードの訪問回数を管理しています(訪問回数が多いノードは回避します)。
2. ランダム要素の導入: 最短経路だけではなく、一定の確率で他の接続ノードを選択し、袋小路を避けるように工夫しています。
3. バックトラッキング: 袋小路に入った場合に、直前のノードに戻るバックトラッキングの処理を行っています。

■プログラムの全体の流れ
1.データベース接続を確立。
2.出発ノードと目的ノードを設定。
3.CSVファイルを作成し、ヘッダーを記述。
4.ランダムウォークを実行し、ノード間を移動。
5.各ノードで緯度・経度を取得し、CSVファイルに書き込む。
6.目的ノードに到達するまで、最も適切な接続ノードを選択し続ける。

このプログラムは、ランダムウォークアルゴリズムをPostGISとPostgreSQLの地理情報を活用して実現しており、目的地に到達するまでの経路を記録する機能を持っています。

*/

package main

import (
	"database/sql"
	"encoding/csv"
	"fmt"
	"log"
	"math"
	"math/rand"
	"os"

	_ "github.com/lib/pq"
)

const (
	connStr = "user=postgres password=password host=127.0.0.1 port=15432 dbname=hakata_db sslmode=disable"
)

func main() {
	// PostgreSQLデータベースへの接続
	db, err := sql.Open("postgres", connStr)
	if err != nil {
		log.Fatal(err)
	}
	defer db.Close()

	// 出発ノードを設定(天神駅の最も近いノードを選択するなど)

	startNode := 32758 // 例として、初期ノードIDを設定
	goalNode := 40922  // 目的ノードIDを天神駅に設定

	// CSVファイルを作成
	csvFile, err := os.Create("random_walk_result.csv")
	if err != nil {
		log.Fatal(err)
	}
	defer csvFile.Close()

	writer := csv.NewWriter(csvFile)
	defer writer.Flush()

	// CSVのヘッダーを作成
	err = writer.Write([]string{"Node", "Latitude", "Longitude"})
	if err != nil {
		log.Fatal(err)
	}

	// ランダムウォークの実行
	err = performRandomWalk(db, startNode, goalNode, writer)
	if err != nil {
		log.Fatal(err)
	}

	fmt.Println("Random walk completed and saved to CSV.")
}

// performRandomWalk ランダムウォークを実行し、結果をCSVに保存
func performRandomWalk(db *sql.DB, startNode, goalNode int, writer *csv.Writer) error {
	currentNode := startNode
	visited := make(map[int]int) // ノードの訪問回数を追跡
	pathStack := []int{}         // バックトラッキング用のスタック

	for currentNode != goalNode {
		// 現在のノードの緯度と経度を取得し、CSVに書き込み
		lat, lon, err := getNodeCoordinates(db, currentNode)
		if err != nil {
			return err
		}
		fmt.Printf("Node %d: Latitude: %f, Longitude: %f\n", currentNode, lat, lon)
		err = writer.Write([]string{
			fmt.Sprintf("%d", currentNode),
			fmt.Sprintf("%f", lat),
			fmt.Sprintf("%f", lon),
		})
		if err != nil {
			return err
		}
		writer.Flush()

		// 現在のノードを訪問済みとして記録
		visited[currentNode]++

		// 次のノードを取得(目的地に近づくノードを優先)
		nextNode, err := getClosestConnectedNode(db, currentNode, goalNode, visited)
		if err != nil || visited[nextNode] > 2 { // 訪問回数が多い場合バックトラック
			if len(pathStack) > 0 {
				// スタックから一つ戻る
				fmt.Printf("Backtracking from Node %d\n", currentNode)
				currentNode = pathStack[len(pathStack)-1]
				pathStack = pathStack[:len(pathStack)-1]
				continue
			}
			return fmt.Errorf("no viable path found")
		}

		// 次のノードをスタックに追加
		pathStack = append(pathStack, currentNode)
		currentNode = nextNode
	}

	// ゴールノードの座標をCSVに保存
	lat, lon, err := getNodeCoordinates(db, goalNode)
	if err != nil {
		return err
	}
	fmt.Printf("Arrived at Goal Node %d: Latitude: %f, Longitude: %f\n", goalNode, lat, lon)
	err = writer.Write([]string{
		fmt.Sprintf("%d", goalNode),
		fmt.Sprintf("%f", lat),
		fmt.Sprintf("%f", lon),
	})
	if err != nil {
		return err
	}
	writer.Flush()

	return nil
}

// getClosestConnectedNode 現在のノードから接続されているノードのうち、目的地に最も近いノードを取得
func getClosestConnectedNode(db *sql.DB, currentNode, goalNode int, visited map[int]int) (int, error) {
	// 現在のノードから接続されているノードを取得
	query := `
		SELECT target 
		FROM ways 
		WHERE source = $1
		UNION 
		SELECT source 
		FROM ways 
		WHERE target = $1;
	`
	rows, err := db.Query(query, currentNode)
	if err != nil {
		return 0, err
	}
	defer rows.Close()

	// 接続ノードのリストを作成
	type NodeDistance struct {
		ID       int
		Distance float64
	}
	var connectedNodes []NodeDistance
	for rows.Next() {
		var node int
		if err := rows.Scan(&node); err != nil {
			return 0, err
		}

		// ノードの座標を取得
		lat, lon, err := getNodeCoordinates(db, node)
		if err != nil {
			return 0, err
		}
		// 目的地との距離を計算
		goalLat, goalLon, err := getNodeCoordinates(db, goalNode)
		if err != nil {
			return 0, err
		}
		distance := calculateDistance(lat, lon, goalLat, goalLon)

		// 訪問回数が少ないノードを優先
		if visited[node] < 3 {
			connectedNodes = append(connectedNodes, NodeDistance{
				ID:       node,
				Distance: distance,
			})
		}
	}

	if len(connectedNodes) == 0 {
		return 0, fmt.Errorf("no connected nodes found")
	}

	// 目的地に最も近いノードを優先して選択
	var closestNode NodeDistance
	minDistance := math.MaxFloat64
	for _, node := range connectedNodes {
		// 最も目的地に近いノードを選択
		if node.Distance < minDistance {
			closestNode = node
			minDistance = node.Distance
		}
	}

	// ランダム選択の要素を追加(袋小路を避けるため)
	if len(connectedNodes) > 1 && rand.Float64() < 0.2 {
		return connectedNodes[rand.Intn(len(connectedNodes))].ID, nil
	}

	return closestNode.ID, nil
}

// getNodeCoordinates ノードIDに対応する緯度と経度を取得
func getNodeCoordinates(db *sql.DB, nodeID int) (float64, float64, error) {
	var lat, lon float64
	query := `
		SELECT ST_Y(the_geom) AS lat, ST_X(the_geom) AS lon
		FROM ways_vertices_pgr
		WHERE id = $1;
	`
	err := db.QueryRow(query, nodeID).Scan(&lat, &lon)
	if err != nil {
		return 0, 0, err
	}
	return lat, lon, nil
}

// calculateDistance 2つの座標間の距離を計算
func calculateDistance(lat1, lon1, lat2, lon2 float64) float64 {
	rad := func(deg float64) float64 {
		return deg * math.Pi / 180
	}
	lat1Rad, lon1Rad := rad(lat1), rad(lon1)
	lat2Rad, lon2Rad := rad(lat2), rad(lon2)
	dlat := lat2Rad - lat1Rad
	dlon := lon2Rad - lon1Rad
	a := math.Sin(dlat/2)*math.Sin(dlat/2) + math.Cos(lat1Rad)*math.Cos(lat2Rad)*math.Sin(dlon/2)*math.Sin(dlon/2)
	c := 2 * math.Atan2(math.Sqrt(a), math.Sqrt(1-a))
	const EarthRadius = 6371e3 // 地球の半径(メートル)
	return EarthRadius * c
}

2024,江端さんの技術メモ

上記のDBはpostGISで構成されており、ways, ways_vertices_pgr があり、SQLでダイクストラ計算が行えます。

ここで以下のようなステップを行うランダムウォークをpostGISを使ったpostgreSQLを使ってGo言語で作成してみました。

(Step.1) ランダムォークの終了地点は、天神駅とする

(Step.2) 出発地点を乱数デ適当に選んで下さい。出発地点は、天神駅から1km程度の地点を適当に抽出する

(Step.3)天神駅の座標から最も近いNode(以下、これを天神ノードといいます)と、出発地点から最も近いNode(以下、これを出発ノードと言います)を抽出する

(Step.4) 出発ノードと天神ノードを直線で繋いで、出発ノードから-90度から+90度の範囲のノードを一つ選ぶ

(Step.5)上記で見つかったノードを次の出発ノードとする。正し、すでに出発ノードとしたノードを、新しい出発ノードにはできないようにする

(Step.6)上記(Step.4)に戻って、天神ノードに到着するまで、出発ノードを更新し続ける。なお更新される出発ノードは全てノド番号と緯度と座標を表示させる

というようなアルゴリズムで、ランダムウォークを作ろうとしたのですが、ダメでした。

一言でいうと「どんつき」の路地に入る→袋小路に入ってしまって、そこで動けなくなるからです。

これを回避するアルゴリズムを考えるの面倒になって、pgr_ksp()を使うことにしました。

pgr_ksp()は、最短経路問題の1、2、3、4番目を出してくれるので、これで1000番目とかを仕えば、そこそこランダムウォークになるかな、と思いまして。

// G:\home\ebata\hakata\src\others\main25\main.go

package main

import (
	"database/sql"
	"encoding/csv"
	"fmt"
	"log"
	"os"

	_ "github.com/lib/pq"
)

const (
	connStr = "user=postgres password=password host=127.0.0.1 port=15432 dbname=hakata_db sslmode=disable"
)

type PathResult struct {
	PathID int
	Node   int
	Lat    float64
	Lon    float64
}

func main() {
	// PostgreSQLデータベースへの接続
	db, err := sql.Open("postgres", connStr)
	if err != nil {
		log.Fatal(err)
	}
	defer db.Close()

	// すべてのパス結果を一度に取得
	results, err := getPathResults(db)
	if err != nil {
		log.Fatal(err)
	}

	// 各パスIDの結果をCSVに保存
	pathIDs := []int{1, 5, 10, 100, 1000, 10000}
	for _, pathID := range pathIDs {
		filename := fmt.Sprintf("path_id_%d.csv", pathID)
		err := saveResultsToCSV(results, pathID, filename)
		if err != nil {
			log.Fatal(err)
		}
	}

	fmt.Println("CSV files created successfully.")
}

// getPathResults すべてのパス結果を取得
func getPathResults(db *sql.DB) ([]PathResult, error) {
	// クエリで指定のpath_idの結果をすべて取得
	query := `
		WITH ksp_result AS (
			SELECT * FROM pgr_ksp(
				'SELECT gid AS id, source, target, cost FROM ways',
				35382,  -- 出発ノードID
				40922,  -- 目的ノードID
				10000,   -- 上位10000本の経路を取得
				false   -- エッジの向きを考慮しない
			)
		)
		SELECT ksp.path_id, ksp.node, ST_Y(wv.the_geom) AS lat, ST_X(wv.the_geom) AS lon
		FROM ksp_result AS ksp
		JOIN ways_vertices_pgr AS wv ON ksp.node = wv.id
		WHERE ksp.path_seq > 0
		ORDER BY ksp.path_id, ksp.path_seq;
	`

	// クエリの実行
	rows, err := db.Query(query)
	if err != nil {
		return nil, err
	}
	defer rows.Close()

	// 結果をメモリに保存
	var results []PathResult
	for rows.Next() {
		var result PathResult
		err := rows.Scan(&result.PathID, &result.Node, &result.Lat, &result.Lon)
		if err != nil {
			return nil, err
		}
		results = append(results, result)
	}

	return results, nil
}

// saveResultsToCSV 指定されたpath_idの結果をCSVに保存
func saveResultsToCSV(results []PathResult, pathID int, filename string) error {
	// CSVファイルの作成
	file, err := os.Create(filename)
	if err != nil {
		return err
	}
	defer file.Close()

	writer := csv.NewWriter(file)
	defer writer.Flush()

	// CSVのヘッダーを作成
	err = writer.Write([]string{"Node", "Latitude", "Longitude"})
	if err != nil {
		return err
	}

	// 結果をCSVに書き込む
	for _, result := range results {
		if result.PathID == pathID {
			err := writer.Write([]string{
				fmt.Sprintf("%d", result.Node),
				fmt.Sprintf("%f", result.Lat),
				fmt.Sprintf("%f", result.Lon),
			})
			if err != nil {
				return err
			}
		}
	}

	return nil
}

 

1番目、10番目、1000番目を表示してみました。
まあ、ランダムウォークという感じではないですが、私たちが移動するときも、まあ、こんな感じだと思うので、とりあえずこれで進めようかなと思います。

2024,江端さんの技術メモ

<!-- Way: Dummy Railway Line for JR篠栗線 -->
<way id="123456800" version="1" timestamp="2024-08-22T06:13:26Z" changeset="99999999" uid="999999" user="kobore.net">
<nd ref="1882462094" /> <!-- 博多駅 -->
<nd ref="1239476659" /> <!-- 吉塚駅 -->
<nd ref="1843493667" /> <!-- 柚須駅 -->
<nd ref="1843644210" /> <!-- 伊賀駅 -->
<nd ref="1842266870" /> <!-- 原町駅 -->
<nd ref="5445230442" /> <!-- 長者原駅 -->
<nd ref="6072021673" /> <!-- 門松駅 -->
<nd ref="796761953" /> <!-- 篠栗駅 -->
<nd ref="796761953" /> <!-- 筑前山手駅 -->
<nd ref="796761953" /> <!-- 城戸南蔵院前駅 -->
<tag k="railway" v="dummyrail" />
<tag k="name" v="Dummy Railway Line JR篠栗線" />
</way>

このデータを、

<relation id="17960683" version="1" timestamp="2024-08-22T06:13:26Z" changeset="155587618" uid="109705" user="gscholz">
<member type="way" ref="1309534187" role="outer"/>
<member type="way" ref="1309192923" role="inner"/>
<member type="way" ref="1309192922" role="inner"/>
<member type="way" ref="1309192926" role="inner"/>
<member type="way" ref="1309192929" role="inner"/>
<member type="way" ref="1309192924" role="inner"/>
<member type="way" ref="1309192925" role="inner"/>
<member type="way" ref="1309192928" role="inner"/>
<member type="way" ref="733239311" role="inner"/>
<member type="way" ref="1309192927" role="inner"/>
<tag k="highway" v="pedestrian"/>
<tag k="surface" v="gravel"/>
<tag k="type" v="multipolygon"/>
</relation>

この辺に置いた

</osm>

こうしたら、

G:\home\ebata\hakata\hakata_db>osm2pgrouting -f hakata_ex.osm -c mapconfig_for_ex.xml -d hakata_ex_db -U postgres -h localhost -p 15432 -W password 

でDBの作成に失敗しました。

wayのタグはwayのタグの集まっているところに置いておかないと、DBの作成で失敗するようです。

JOSMで地図を改ざんして道路と鉄道を交えた最⼩コスト経路探索をやってみよう-完全版

 

 

 

2024,江端さんの技術メモ

(Step 1)地図をスキャンする

(Step 2)SAIで地図を切り取って、さらに追加したレイヤに点を打つ

(Step 3)点を打ったレイヤだけをjpgでセーブする

page61.jpgの内容

(Step 4)このjpegを左上の座標と右上の座標を指定して、ChatGPTに添付して抽出をお願いする。

添付したファイルは地図のある建物の位置を示しているものです。
左上の経度、緯度は、130.41636503641814, 33.59644753614503
右下の経度、緯度は、130.42569034464057, 33.58606144150355
となっています。
添付したファイルに記載された点の座標の経度、緯度を算出して下さい。

(Step 5)で、色々失敗するので、何度か修正をお願いする。たとえば以下のような感じ。

左上の経度、緯度は、130.41636503641814, 33.59644753614503
右下の経度、緯度は、130.42569034464057, 33.58606144150355
ファイルには、赤色、緑色の点が記載されています。この点の座標を、上記の左上の緯度経度と右下の移動経度から、線形計算をして算出して下さい。

私はこの添付ファイルに20~30の点しか付けていませんが、頂いた点は非常に多いです。見直して下さい。

同じ点が2~3個で表示されているようです。点と点の距離が20メートル以内であれば、点を1つにして下さい。

(Step 6)ChatGPTが作ってくれたcsvファイルをQGIS等で表示して確認する。

まあ、とりあえず、位置情報だけは、ChatGPTの協力を得て、取り出すことができました。

しかし、店舗名の紐づけは、どうしようもありませんね。これを何とかできる人は、私を助けて下さい。

======

と思っていたら、以外なほど簡単にできました。

QGISで上記のCSVファイルを示すことができましたが、csvファイルとの対応が分かりません。
例えば、
#,Latitude,Longitude
1,33.59181234, 130.4199366
2.33.59123927, 130.4190329
3.33.59106178 130.4217193
のようにして、QGISで表示される点に#の番号を表示させるにはどうしたら良いですか?

とChatGPTに聞いたら教えてくれました。

まず、普通にcsvを取り込む

点が見えにくいので、大きな赤丸にする。

csvのレイヤを選択して → プロパティ → シンポロジで "dot red" を選ぶ(趣味の問題だが)

で、再び、

csvのレイヤを選択して → プロパティ → ラベル → (上のほうに"なし"と記載されているメニューから) 「単一定義」を選択。あとはよく分からんが、適用を押す。

おお! ちゃんと番号が付与されている。

これで、対応付けできるし、なんなら地物名称も入れれるかもしれない。

2024,江端さんの技術メモ

Teamsの会議でビデオの音声を共有するには、以下の手順を実行してください:

  1. 会議を開始するか参加
    • Teamsの会議を開始する、もしくは既存の会議に参加してください。
  2. 「画面を共有」ボタンを押す
    • 会議ウィンドウの中央下部にある「共有」ボタン(四角と上向きの矢印アイコン)をクリックします。
  3. 共有するウィンドウや画面を選択
    • 共有したいビデオを再生しているウィンドウや画面、もしくはデスクトップ全体を選択します。
  4. 「コンピューターの音声を含める」オプションを有効化
    • 画面共有の際に、画面の上部に「コンピューターの音声を含める」というチェックボックスが表示されます。これにチェックを入れることで、ビデオの音声も一緒に共有できます。

  1. ビデオを再生
    • 画面共有を開始した後、会議参加者にビデオの映像と音声が共有されます。

この設定により、Teams会議の参加者全員があなたのコンピューターの音声を聞くことができるようになります。

2024,江端さんの技術メモ

「Redisをインストールすることなく、Windows10で使う方法」がある、と、ChatGPTに言われて「ホンマかいな?」と思いながらやってみたのですが、さっくり動きました。

$ docker run --name redis -p 6379:6379 -d redis

これだけで、Windows上で動かしているGoプログラムの中のredisが動いています(理屈は分かりません)。

Redisの動作確認
Dockerコンテナ内のRedisが動作しているか確認するためには、以下のコマンドを実行します。

$ docker exec -it redis redis-cli ping

-----

また、「docker run --name redis -p 6379:6379 -d redis を次回再起動するには、どうしたら良いですか」と聞いてみたところ、

Redisコンテナを再起動するコマンド

  1. コンテナを再起動する コンテナが停止している場合、次のコマンドで再起動します:
    docker start redis
  2. コンテナが実行中かどうかを確認する docker psコマンドを使って、Redisコンテナが正しく動作しているか確認できます:
    docker ps

    Redisコンテナが実行中であれば、一覧に表示されます。

だ、そうです。ちなみに、

もし、Dockerが起動するたびにRedisコンテナを自動的に再起動させたい場合は、--restartオプションを使って自動再起動ポリシーを設定できます。次回docker runを実行するとき、--restartオプションを追加しておくと、自動再起動が設定されます。

$ docker run --name redis -p 6379:6379 -d --restart unless-stopped redis

私は、redisが必要のない時に立ち上げたくはないので、上記の設定はしていません。

以上

 

2024,江端さんの技術メモ

(1)地図DBに接続されている状態であること。

(2)go.modが入っていること

(3)prune_isonodes.go の依存パッケージのインストール方法が変っていること

// 依存パッケージ  $ go get github.com/lib/pq
//↓
G:\home\ebata\hakata\hakata_db>go get github.com/lib/pq
G:\home\ebata\hakata\hakata_db>go install github.com/lib/pq
go: 'go install' requires a version when current directory is not in a module
        Try 'go install github.com/lib/pq@latest' to install the latest version
G:\home\ebata\hakata\hakata_db>go install github.com/lib/pq@latest
package github.com/lib/pq is not a main package

(4)prune_isonodes.goに書き込まれたターゲットのdb名称を変更すること

    dbParam := fmt.Sprintf("host=localhost port=%d user=postgres password=password dbname=hakata_db sslmode=disable", *dbPort)
"go run prune_isonodes.go"で孤立ノードが削除されます。
G:\home\ebata\hakata>go run prune_isonodes.go
[ノード#1からのスキャン] 到達数: 60781 未到達: 1108
削除対象のノード (1108個):
[57493 59415 46635 51008 46336 59154 48022 59511 39412 22789 59166 59535 59531 53964 31094 45646 55303 45386 51069 59394 59129 59175 59215 59190 59297 45602 61062 50473 53933 51663 60258 17471 17566 59406 48023 60228 49200 49186 58539 60103 51007 59546 51006 57326 59547 45680 59250 59598 46354 47335 47010 58527 54516 31095 59323 45327 45676 60125 46868 53953 60226 32888 50574 60213 18763 59290 60114 59520 60264 59393 59538 32894 51854 59477 59061 59245 47003 59615 53958 60245 59325 61661 59064 58529 52435 47547 47001 60231 39841 61645 59130 45613 46221 46356 59059 28051 59191 59237 59286 59125 59229 59548 42380 59074 53967 59424 58543 59305 59524 59413 45664 53966 47343 59292 50984 47012 60109 58530 58678 60156 45606 59184 49184 53951 37012 42383 46334 59195 59362 45599 59601 60266 29342 59056 47337 45699 45648 59635 59293 5237 47004 22462 61660 45621 61655 53239 49196 59077 59112 60255 53934 42382 59179 56864 45620 52871 46817 59216 59512 5238 19595 48020 46326 59610 60221 47340 61734 61737 45647 45638 58568 59132 50582 59288 57488 61668 59244 31964 51274 59161 47315 19587 51327 45628 49198 59576 59550 51158 45674 53950 59401 45669 59226 46226 58528 35214 59348 37000 53921 18761 51267 60256 53968 17485 51273 48352 50363 59501 59211 57485 59636 59084 50580 19594 31093 59055 53928 37013 45600 60104 53970 60214 59234 28899 60110 50578 47334 32889 59505 60219 41557 60239 50267 53939 29424 59239 60131 57486 39415 59343 49201 53243 45662 60252 45665 60217 46223 61658 46864 58083 59334 45610 60155 59327 59131 59542 53940 46342 59478 59388 59392 59414 59296 45653 59879 59169 61644 53238 59500 59549 45655 57484 45614 59360 59345 59316 59135 59122 59568 59176 59251 60251 50577 59158 59509 45618 60102 22775 59196 59361 59203 52870 53686 59595 58537 59491 45670 59142 53975 42610 60337 53947 53919 53973 59618 59145 39409 24519 45654 59289 53931 46220 47002 46866 59355 45625 60260 60257 53962 53235 37018 57327 37008 59308 45962 53937 53932 59354 59160 59194 37023 53957 53961 59384 59342 37020 61724 59349 59621 53237 59537 47316 51665 59173 22114 32892 46007 46327 59609 53942 59157 46974 59486 59495 59356 51853 53226 50355 50979 59137 54531 59508 51206 59205 59487 40314 53926 45659 59221 50983 45601 60262 59526 36997 53955 59232 59163 18757 60237 45675 57324 53930 59209 53242 22115 45698 19576 59206 59144 59240 59515 59522 45673 59189 45656 59220 59386 59407 59399 59604 45637 61643 57332 59111 59385 46331 60240 59162 59303 54923 36998 60232 59299 56623 59498 59607 60218 55301 49189 31966 59338 57906 45690 45322 59423 45615 59301 60716 60115 53687 60121 45660 37019 45619 50354 51165 52437 48272 32059 60243 59307 59306 53241 59228 58534 59602 59319 59082 60238 50474 45666 59310 59391 10298 60116 59496 59417 59081 61871 45617 59411 59078 59597 37016 59558 46949 40315 46324 60119 59553 59219 59422 32897 59507 59243 59523 59114 45639 46329 59287 47339 59350 60242 53959 59247 43492 59118 61738 61880 45325 59076 53948 53972 5240 59324 40316 59155 59242 59187 59426 60224 59479 45624 59183 59128 59246 58580 45678 39840 46862 59476 59344 59201 43246 59620 59492 59113 59351 59612 59418 59224 60216 60222 59121 45645 59217 59494 59359 57904 53971 53922 59159 45679 59115 55300 59347 58581 50592 59561 59529 60248 59592 59563 36996 45663 60215 53943 59140 59480 59880 45603 59060 59197 53949 58351 60227 44923 58541 56622 37007 40313 60122 46975 45627 59565 59126 58082 57487 59593 49190 57914 60233 45634 59554 45671 45324 59198 59400 46330 59124 59594 17541 19581 60120 19593 47033 38208 5236 59412 60334 28944 53927 61725 5239 59171 60154 46335 59180 59285 45640 60254 59427 61061 45644 45414 31965 53918 57490 59284 17522 57489 52439 59138 59534 60261 59489 50589 54396 47902 59398 37011 59225 45658 51068 59313 54423 50583 61662 57910 45689 50586 59613 57325 53965 47904 58532 59346 60241 59182 60124 60267 59231 50982 52436 61666 61495 46315 61665 58544 32891 51174 59341 59599 46227 59235 37005 32895 53240 61872 59300 45436 60111 53935 46225 59419 60338 59227 59622 46867 59167 46339 59283 54517 60246 45616 53236 45657 48021 59249 51330 45650 22452 54922 59172 59395 59409 46950 59202 47901 60106 61642 46861 59533 10301 59062 59352 59403 59544 52438 45608 59340 57909 59617 59605 45642 53234 54530 45635 50360 59207 38209 59192 52440 37014 59320 56855 45672 59482 54525 44926 59141 19597 59545 59119 59248 42064 57331 61672 59199 42381 59164 45641 45677 59566 59499 17519 46341 50595 57323 51070 58579 50576 59567 58695 60108 49188 45609 59065 31967 22778 46333 59136 59185 59337 59624 17539 59506 29339 46222 52872 48010 59485 51326 59357 50978 53944 59490 45626 53941 50585 57911 45629 56621 49194 48011 39411 58538 53689 59302 59304 53923 59497 61663 37024 59282 59358 59188 32890 46332 59521 59483 59488 60220 50361 60234 56713 59430 53945 47336 59335 22783 50991 37022 32896 53936 50575 59133 59551 53920 57913 58526 45661 59608 59170 59389 50587 59058 52447 54398 59339 50992 53946 60259 10300 45611 53954 60235 15200 45651 46355 59397 93 46863 60236 61657 36995 50591 59311 59193 59156 45323 32893 59353 10297 60126 61063 53960 37021 59314 59513 59402 49192 22776 50980 59484 42379 59063 59390 59315 45633 59562 37017 15165 46337 60107 57482 53925 39413 59637 59429 59318 59208 59168 49195 45623 59552 58533 49193 22769 59143 59222 37009 47903 22782 37004 50594 45604 60112 61060 39408 54397 59085 53938 59174 59564 59493 59421 53963 19588 46318 59178 60153 59241 61659 59600 45630 43493 43480 59083 59123 60244 59139 46224 59186 53924 58542 59117 59606 60113 59218 50590 59238 44919 59181 61656 38207 59416 50579 51662 57329 60223 59312 24502 58536 6768 46321 59134 61654 55302 59322 59057 17560 39414 45649 49191 50334 51329 59387 59326 50879 53312 59420 37010 60229 61059 32898 45643 49185 37003 59611 59410 59619 59404 59054 59514 39416 59317 53969 53974 53952 55795 18764 50596 3068 59532 47314 59295 60249 47014 50588 39574 59294 60265 60263 37015 45326 3066 45668 37006 43491 51325 10299 61667 59536 60105 45605 58535 59503 45631 57905 51664 59556 60118 59212 45622 22780 45667 59309 22768 59502 57912 59616 46317 53956 59396 57483 45632 59405 59210 59321 47017 57491 59623 6767 44924 45607 58509 60253 59213 50593 59177 58531 48351 57328 3067 45319 45652 60225 50993 51328 50584 59120 61494 31036 49199 22791 59328 45612 59075 37574 22443 59577 59408 59614 39410 60117 59603 59230 53688 50581 59557 40312 59291 59596 45697 59080 61879 59336 59236 59214 59079 60230 37002 36999 46366 46865 19584 59116 49187 58540 45636 53929 52794 59200 60247 58084 60618 59428 52446 59233 46338 59629 18762 49197 59204 59165 59127 60123 60250 59223 61664 59425 57492 59298 37001 59504 19592 57330]
削除を実行します。
完了しました。

2024,江端さんの技術メモ

この地図のプロパティは、424 x 379

で、

set datafile separator "," # データの区切り文字をコンマに指定
set title "Cluster Data with Map" # グラフタイトル
set xlabel "X-axis" # X軸のラベル
set ylabel "Y-axis" # Y軸のラベル
set zlabel "Z-axis" # Z軸のラベル
set xrange [0:1] # X軸の範囲をデータに合わせて調整
set yrange [0:1] # Y軸の範囲をデータに合わせて調整
set zrange [0:1] # Z軸の範囲
set ticslevel 0 # Z軸を平面上に表示
set xyplane at 0 # XY平面をZ=0に設定
set view 40, 30, 1, 1 # 視点の角度とスケールを調整

# 地図画像をXY平面に表示 (画像のX軸、Y軸範囲を設定し、Z軸は0に固定)
splot "img_664a03bcd2376.png" binary filetype=png dx=1.0/424 dy=1.0/379 with rgbimage, \
'Cluster1.csv' using 1:2:3 with points pt 7 ps 2 lc rgb "blue" title 'Cluster 1', \
'Cluster2.csv' using 1:2:3 with points pt 7 ps 2 lc rgb "blue" notitle, \
'Cluster3.csv' using 1:2:3 with points pt 7 ps 2 lc rgb "blue" notitle, \
'Cluster4.csv' using 1:2:3 with points pt 7 ps 2 lc rgb "blue" notitle, \
'Cluster5.csv' using 1:2:3 with points pt 7 ps 2 lc rgb "blue" notitle, \
'Cluster6.csv' using 1:2:3 with points pt 7 ps 2 lc rgb "blue" notitle, \
'Cluster7.csv' using 1:2:3 with points pt 7 ps 2 lc rgb "blue" notitle, \
'Cluster8.csv' using 1:2:3 with points pt 7 ps 2 lc rgb "blue" notitle, \
'Cluster9.csv' using 1:2:3 with points pt 7 ps 2 lc rgb "blue" notitle, \
'Cluster10.csv' using 1:2:3 with points pt 7 ps 2 lc rgb "blue" notitle, \
'Cluster11.csv' using 1:2:3 with points pt 7 ps 2 lc rgb "blue" notitle, \
'Cluster12.csv' using 1:2:3 with points pt 7 ps 2 lc rgb "blue" notitle, \
'Cluster13.csv' using 1:2:3 with points pt 7 ps 2 lc rgb "blue" notitle, \
'Cluster14.csv' using 1:2:3 with points pt 7 ps 2 lc rgb "blue" notitle, \
'Cluster15.csv' using 1:2:3 with points pt 7 ps 2 lc rgb "blue" notitle, \
'Cluster16.csv' using 1:2:3 with points pt 7 ps 2 lc rgb "blue" notitle, \
'Cluster17.csv' using 1:2:3 with points pt 7 ps 2 lc rgb "blue" notitle, \
'Cluster18.csv' using 1:2:3 with points pt 7 ps 2 lc rgb "blue" notitle, \
'Cluster19.csv' using 1:2:3 with points pt 7 ps 2 lc rgb "blue" notitle, \
'Cluster20.csv' using 1:2:3 with points pt 7 ps 2 lc rgb "blue" notitle, \
'Cluster21.csv' using 1:2:3 with points pt 7 ps 2 lc rgb "blue" notitle, \
'Cluster22.csv' using 1:2:3 with points pt 7 ps 2 lc rgb "blue" notitle, \
'Cluster23.csv' using 1:2:3 with points pt 7 ps 2 lc rgb "blue" notitle, \
'Cluster24.csv' using 1:2:3 with points pt 7 ps 2 lc rgb "blue" notitle, \
'Cluster25.csv' using 1:2:3 with points pt 7 ps 2 lc rgb "blue" notitle, \
'Cluster26.csv' using 1:2:3 with points pt 7 ps 2 lc rgb "blue" notitle, \
'Cluster27.csv' using 1:2:3 with points pt 7 ps 2 lc rgb "blue" notitle, \
'Cluster28.csv' using 1:2:3 with points pt 7 ps 2 lc rgb "blue" notitle, \
'Cluster29.csv' using 1:2:3 with points pt 7 ps 2 lc rgb "blue" notitle, \
'Cluster30.csv' using 1:2:3 with points pt 7 ps 2 lc rgb "blue" notitle, \
'Cluster31.csv' using 1:2:3 with points pt 7 ps 2 lc rgb "blue" notitle, \
'Cluster32.csv' using 1:2:3 with points pt 7 ps 2 lc rgb "blue" notitle, \
'Cluster33.csv' using 1:2:3 with points pt 7 ps 2 lc rgb "blue" notitle, \
'Cluster34.csv' using 1:2:3 with points pt 7 ps 2 lc rgb "blue" notitle, \
'Cluster35.csv' using 1:2:3 with points pt 7 ps 2 lc rgb "blue" notitle, \
'Cluster36.csv' using 1:2:3 with points pt 7 ps 2 lc rgb "blue" notitle, \
'Cluster37.csv' using 1:2:3 with points pt 7 ps 2 lc rgb "blue" notitle, \
'Cluster38.csv' using 1:2:3 with points pt 7 ps 2 lc rgb "blue" notitle, \
'Cluster39.csv' using 1:2:3 with points pt 7 ps 2 lc rgb "blue" notitle, \
'Cluster40.csv' using 1:2:3 with points pt 7 ps 2 lc rgb "blue" notitle, \
'Cluster41.csv' using 1:2:3 with points pt 7 ps 2 lc rgb "blue" notitle

cluster1.csvの内容はこんな感じ
0.63, 0.65, 0.58
0.62, 0.65, 0.56
0.61, 0.65, 0.56
0.61, 0.65, 0.56
0.62, 0.65, 0.56
0.62, 0.66, 0.55
0.63, 0.65, 0.55
0.62, 0.65, 0.56
0.62, 0.66, 0.57
0.62, 0.65, 0.58
0.62, 0.65, 0.58
(後略)