changeset 15:3673b244ad94

add logistic.go
author pyon@macmini
date Sun, 21 Jan 2018 17:31:21 +0900
parents 001e2aa380ad
children 38b64afbaf79
files src/logistic.go
diffstat 1 files changed, 1157 insertions(+), 0 deletions(-) [+]
line wrap: on
line diff
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/logistic.go	Sun Jan 21 17:31:21 2018 +0900
@@ -0,0 +1,1157 @@
+/*
+ * ロジスティック回帰分析
+ *                               Last Change: 2018-01-21 Sun 17:15:18. 
+ */
+package main
+
+import (
+	"bufio"
+	"fmt"
+	"image"
+	"image/color"
+	"image/png"
+	"log"
+	"math"
+	"os"
+	"strconv"
+	"strings"
+	"time"
+)
+
+var RMAX int        // 最大繰り返し回数
+var NN, MM int		// 説明変数の次元
+var lr float64		// 学習率
+var vb []float64	// 重み(切片と傾き)ベクトル
+var db float64      // vb の微増分(デルタ)
+var mx [][]float64  // 説明変数
+var vy []float64    // 結果変数
+var start time.Time
+
+func init() {
+	start = time.Now()
+
+    RMAX = 500
+
+	lr = 0.005
+    db = 0.0001
+
+	NN, MM  = 1000, 3 // i, j
+	vb = make( []float64, MM )
+	vy = make( []float64, NN )
+	mx = make( [][]float64, NN )
+	for i := range mx {
+		mx[i] = make( []float64, MM )
+		mx[i][0] = 1.0
+	}
+
+    for j := range vb {
+        vb[j] = 1.0
+    }
+
+	f, err := os.Open( "input.txt" )
+	if err != nil {
+		log.Fatal( err )
+	}
+	scanner := bufio.NewScanner( f )
+	var buf []string
+	i := 0
+	for scanner.Scan() {
+		buf = strings.Split( scanner.Text(), "\t" )
+		for j := range vb {
+			if s, err := strconv.ParseFloat( buf[j], 64 ); err != nil {
+				log.Fatal( err )
+			} else {
+                if j == MM - 1 {
+                    vy[i] = s
+                } else {
+                    mx[i][j+1] = s
+                }
+			}
+		}
+		i++
+	}
+	if err := scanner.Err(); err != nil {
+        log.Fatal( err )
+	}
+	f.Close()
+}
+
+func main() {
+
+	for n := 0; n < RMAX; n++ {
+        dcdb := f_dCdb() // 対数尤度の勾配を計算
+        odds := f_odds() // オッズを計算
+
+        for j := range dcdb {
+			vb[j] -= dcdb[j] * lr // 重み係数を更新
+		}
+
+		fmt.Printf( "%04d> %7.2f  %6.2f  %3.1f\n", n, dcdb, vb, odds ) // 出力
+	}
+	fmt.Println( time.Now().Sub( start ) )
+	//makegraph()
+}
+
+func f_dCdb() []float64 {
+
+	dcdb := make( []float64, MM )
+
+    for p := range vb {
+        for i:= 0; i < NN; i++ {
+            var xb0, xb float64
+            for j:= range vb {
+                d := 0.0
+                if j == p {
+                    d = db
+                }
+                xb0 += mx[i][j] * vb[j]
+                xb  += mx[i][j] * ( vb[j] + d )
+            }
+            s0 := ( vy[i] - 1 ) * xb0 - math.Log1p( math.Exp( - ( xb0 ) ) )
+            s  := ( vy[i] - 1 ) * xb  - math.Log1p( math.Exp( - ( xb  ) ) )
+            dcdb[p] += - ( s - s0 ) / db
+        }
+    }
+	return dcdb
+}
+
+func f_odds() []float64 {
+    odds := make( []float64, MM )
+    for j:= range vb {
+        odds[j] = math.Exp( vb[j] )
+    }
+    return odds
+}
+
+func makegraph() {
+	const width, height = 256, 256
+	img := image.NewNRGBA( image.Rect( 0, 0, width, height ) )
+	for y := 0; y < height; y++ {
+		for x := 0; x < width; x++ {
+			img.Set( x, y, color.NRGBA{
+				R: uint8( ( x + y ) & 255 ),
+				G: uint8( ( x + y ) << 1 & 255 ),
+				B: uint8( ( x + y ) << 2 & 255 ),
+				A: 255,
+			} )
+		}
+	}
+
+	f, err := os.Create( "image.png" )
+	if err != nil {
+		log.Fatal( err )
+	}
+
+	if err := png.Encode( f, img ); err != nil {
+		f.Close()
+		log.Fatal( err )
+	}
+
+	if err := f.Close(); err != nil {
+		log.Fatal( err )
+	}
+}
+
+/* - input.txt - */
+/*
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	0
+0	0	1
+0	0	1
+0	0	1
+0	0	1
+0	0	1
+0	0	1
+0	0	1
+0	0	1
+0	0	1
+0	0	1
+0	0	1
+0	0	1
+0	0	1
+0	0	1
+0	0	1
+0	0	1
+0	0	1
+0	0	1
+0	0	1
+0	0	1
+0	0	1
+0	0	1
+0	0	1
+0	1	0
+0	1	0
+0	1	0
+0	1	0
+0	1	0
+0	1	0
+0	1	0
+0	1	0
+0	1	0
+0	1	0
+0	1	0
+0	1	0
+0	1	0
+0	1	0
+0	1	0
+0	1	0
+0	1	0
+0	1	0
+0	1	1
+0	1	1
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	0
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	0	1
+1	1	0
+1	1	0
+1	1	0
+1	1	0
+1	1	0
+1	1	0
+1	1	0
+1	1	0
+1	1	0
+1	1	0
+1	1	0
+1	1	0
+1	1	0
+1	1	0
+1	1	0
+1	1	1
+1	1	1
+1	1	1
+1	1	1
+1	1	1
+1	1	1
+1	1	1
+1	1	1
+1	1	1
+1	1	1
+*/