解の範囲に制限を付けて連立方程式を解く【備忘録,R,重回帰,最小二乗法】

ほんとに備忘録なんで期待無しで

でもこういう細かいスクリプトとかって共有されないから検索しても引っかからないのよね。

目的

library(tidyverse)
library(dplyr)
library(RCurl)

# 行為に対して所要時間がついている
# ローカルにDLしたなら
# prob <- read.csv("0.csv",encoding="UTF-8",stringsAsFactors=F)

github_data <- getURL("https://raw.githubusercontent.com/amazongodman/open_folder/master/Simultaneous_equations_2020_06_10.csv")
prob <- read.csv(text = github_data,encoding="UTF-8",stringsAsFactors=F)

prob

  m1 m2 m3 mean
1  A          1
2  B          3
3  A  B       4
4  A  B  C    9
5  C  D      13
6  B  D      11
7  A  D       6
8  A  D      10
9  A  B  D   12

## 空白はNODATAに置き換える
prob$m1[nchar(prob$m1)==0] <- "NODATA"
prob$m2[nchar(prob$m2)==0] <- "NODATA"
prob$m3[nchar(prob$m3)==0] <- "NODATA"

## 要素をマトリックスにする
mat4 <-matrix(c(prob$m1,prob$m2,prob$m3),nrow = length(prob$m1),ncol = 3)

mat4

      [,1] [,2]     [,3]    
 [1,] "A"  "NODATA" "NODATA"
 [2,] "B"  "NODATA" "NODATA"
 [3,] "A"  "B"      "NODATA"
 [4,] "A"  "B"      "C"     
 [5,] "C"  "D"      "NODATA"
 [6,] "B"  "D"      "NODATA"
 [7,] "A"  "D"      "NODATA"
 [8,] "A"  "D"      "NODATA"
 [9,] "A"  "B"      "D" 

## 連立方程式になおす
df <- data.frame(rows=seq_len(nrow(mat4))[row(mat4)], values=c(mat4)) %>%
  mutate(a=1) %>%
  pivot_wider(id_cols="rows", names_from="values", values_from="a", values_fn=list(a=length)) %>%
  mutate_all(~ +!is.na(.)) %>%
  select(-rows) %>%
  select(sort(colnames(.)))

      A     B     C     D NODATA
  <int> <int> <int> <int>  <int>
1     1     0     0     0      1
2     0     1     0     0      1
3     1     1     0     0      1
4     1     1     1     0      0
5     0     0     1     1      1
6     0     1     0     1      1
7     1     0     0     1      1
8     1     0     0     1      1
9     1     1     0     1      0

# 行為の無かったNodata列に数字が入っているのは厄介なので列ごと消す
df2 <-df[,colnames(df)[colnames(df) != "NODATA"] ]
df2$y = prob$mean

連立を解くために最小二乗法を使う。

f<-function(param, x, y){
  sum((y - t(param) %*% t(x))^2)
}

# 重回帰で解けないような場合や、重回帰で係数が0以下になってほしくない場合にはoptimを使う

# 係数の最低値を0.1にする制約を付ける
para <- optim(par=rep(0, 4), lower = 0.1,method = "L-BFGS-B",fn = f, x=df2[,-5], y=as.numeric(df2$y) )

[1] 0.6585317 3.3089433 5.1951257 7.6422791

# 係数の最低値を1にする制約を付ける
para2 <- optim(par=rep(0, 4), lower = 1,method = "L-BFGS-B",fn = f, x=df2[,-5], y=as.numeric(df2$y) )

[1] 1.000000 3.166669 5.166660 7.499996

制約を付けないと重回帰と結果は同じになる

para3 <- optim(par=rep(0, 4),method = "L-BFGS-B",fn = f, x=df2[,-5], y=as.numeric(df2$y) )

[1] 0.6585322 3.3089432 5.1951255 7.6422789

# 切片なしの重回帰で係数を探す
summary(lm(data=df2,y~.+0,conf.int = TRUE))

Coefficients:
  Estimate Std. Error t value Pr(>|t|)    
A   0.6585     0.7096   0.928 0.395970    
B   3.3089     0.7194   4.600 0.005842 ** 
C   5.1951     1.0035   5.177 0.003534 ** 
D   7.6423     0.7194  10.624 0.000128 ***

以上