LoginSignup
8
6

More than 5 years have passed since last update.

@mitama

Rで制約付き最適化問題（最適消費計画）

Last updated at Posted at 2017-10-14

はじめに

Rで最適化問題を解くことができるらしいので，学部レベル(初級)のミクロ経済学の序盤で習う，「予算制約の中で効用（満足感）を最大化する購買量の探索」という効用最大化問題（最適消費計画問題）を例にとり，Rで簡単な制約付き最適化問題を解いたりイロイロしてみます。

◇2018年1月20日追記

ここに載せたことをshinyで作動させてみました。
シンプルな効用最大化問題計算機という感じです。
https://nigimitama.shinyapps.io/MicroEconomics_01_UtilityMaximizationProblem/

中身のコードはこのページの下部（おまけ2）にあります。よろしければ御覧ください。

問題

2財モデルのこんな問題があるとします。

第1財の消費量を$x_1$，第2財の消費量を$x_2$と表す。ある家計の効用関数が$U(x_1,x_2)=x_1^{\frac{2}{5}}x_2^{\frac{3}{5}}$で表されるとする。予算が$100$，第1財の価格が$4$，第2財の価格が$6$のとき，効用を最大化する最適消費量はそれぞれいくらか？

$$\max_{x_1,x_2} \hspace{1em} x_1^{\frac{2}{5}}x_2^{\frac{3}{5}}$$
$$s.t. \hspace{1em} 100 = 4x_1+6x_2$$

ラグランジュの未定乗数法で解ける問題ですが，Rで解いてみます。

Rで解く

{Rsolnp}パッケージを使います。

library(Rsolnp)
ObjFunc = function(x) return( - x[1]^(2/5) * x[2]^(3/5) ) #目的関数
#solnp()は最小化をするようになっているので，目的関数にマイナスを掛ける

ConstFunc = function(x) return( x[1]* 4 + x[2] * 6 ) #制約式の右辺
eq.value <- c(100) #制約式の左辺

x0 <- c(1,1) #決定変数を初期化

solution <- solnp(x0, fun = ObjFunc, eqfun = ConstFunc, eqB = eq.value )

# 最適な｢x_1｣と「x_2｣の値
solution$pars

[1] 10 10

最適な消費量$x_1^* , x_2^* $ は $x_1^* = 10, \ x_2^* = 10$であることが示された

（おまけ1）図示

■効用関数の図示

モノ（$x_1,x_2$）をたくさん買うほど効用$U$は上がる

# 3次元プロット:persp()
x_1 <- 1:50
x_2 <- 1:50
u <- function(x_1,x_2) {x_1^(2/5) * x_2^(3/5)} #効用関数を定義
U <- outer(x_1, x_2, u) #outer()はx_1,x_2に対応したf(x_1,x_2)の値を行列で返す
persp(x_1, x_2, U,
      theta = 30, # 横回転の角度
      phi = 30, # 縦回転の角度
      ticktype = "simple", # 線の種類
      lwd = 0.5, # 線の太さ
      col = F,
      border = 8)

■等高線（無差別曲線）と予算制約線と最適消費点

効用関数の曲面を真上から見て等高線にした図。経済学ではこの等高線を無差別曲線という。
予算制約線（青い直線）と最適消費点（藍色の点）も記しておいた。

無差別曲線と予算制約線が接する点が最適消費点となる。

# 等高線:contour()
contour(x_1, x_2, U, method = "edge", labcex = 1,lwd = 2)
abline(a = 100/6, b = -4/6, lwd = 2, col = "blue") #予算制約線
points(x = 10, y = 10, lwd = 3, col = "darkblue", pch = 16) #最適消費点

（おまけ2）Shinyで効用最大化問題計算機を作る

以上の制約付き最適化問題を解く{Rsolnp}の処理や関数の図示などをshinyで行わせてみる。

demo2.PNG

ui.R
# 日本語の文字列が本文中にあるとdeployできないみたい
# server側のRenderUI内の日本語は問題ないのでコレを利用

library(shinydashboard)
library(shiny)
library(Rsolnp)

header <- dashboardHeader(
  title = "Utility maximization problem"
)

body <- dashboardBody(
    column(width = 3,
           #予算制約
           box(width = NULL,
               uiOutput('budgettext'),
               numericInput("budget", 
                            label = "$$B$$", 
                            value = 100),
               numericInput("p_1", 
                            label = "$$p_1$$", 
                            value = 2),
               numericInput("p_2", 
                            label = "$$p_2$$", 
                            value = 2),
               uiOutput("budgetform")

           ),

           box(width = NULL, 
               #効用関数モジュール
               uiOutput('utilityfunctiontext'),
               withMathJax(),   #TeX表記有効化
               numericInput("a", 
                            label = "$$a$$", 
                            value = 1),
               numericInput("b", 
                            label = "$$b$$", 
                            value = 1),
               numericInput("c", 
                            label = "$$c$$", 
                            value = 1),
               numericInput("d", 
                            label = "$$d$$", 
                            value = 2),
               numericInput("e", 
                            label = "$$e$$", 
                            value = 1),
               numericInput("f", 
                            label = "$$f$$", 
                            value = 2),
               uiOutput("utilityfunction")
               ),
           box(width = NULL,
               #p("変数の範囲（描画範囲）"),
               numericInput("range", 
                            label =  "range of graph", 
                            value = 50)
               )
    ),
    fluidRow(
      column(width = 7,
             #設定された効用最大化問題
             box(width = NULL,
                 uiOutput("Utilitymaximizationproblem")
             ),
             #効用最大化問題の解
             box(width = NULL, status = "warning",
                 uiOutput("maxutility")
                 ),
             #無差別曲線
             box(width = NULL, solidHeader = TRUE,
                 uiOutput('indifferencecurve'),
                 plotOutput("Plot2D", width = 300, height = 300)
                 ),
             box(width = NULL, solidHeader = TRUE,
                 uiOutput('Plot3Dtext'),
                 plotOutput("Plot3D",height = 300)
                 )
             )
      )
)

dashboardPage(
  header,
  dashboardSidebar(disable = TRUE),
  body
)
server.R
library(shinydashboard)
library(shiny)
library(Rsolnp)

# サーバロジックの定義
shinyServer(function(input, output) {

  # 予算制約式box見出し
  output$budgettext <- renderUI({
    withMathJax(helpText("予算制約と価格の設定$$p_1 x_1 + p_2 x_2 \\leq B$$"))
  })
  #指定された予算をTeXでUIに表示
  output$budgetform <- renderUI({
    withMathJax(helpText(paste0('指定された予算制約式:$$',
                                input$p_1,'x_1 + ', input$p_2,'x_2',
                                '\\leq',input$budget,'$$')))
  })


  # 効用関数box見出し
  output$utilityfunctiontext <- renderUI({
    withMathJax(helpText('効用関数のパラメータを設定$$U(x_1,x_2)=a x_1^ \\frac{c}{d} b x_2^ \\frac{e}{f}$$'))
  })
  #指定された効用関数をTeXでUIに表示
  output$utilityfunction <- renderUI({
    withMathJax(helpText(paste0('指定された効用関数:$$U(x_1,x_2)=',
                              input$a,'x_1','^\\frac{',input$c,'}','{',input$d,'}', input$b,'x_2','^\\frac{',input$e,'}','{',input$f,'}','$$')))
  })


  #設定された効用最大化問題
  output$Utilitymaximizationproblem <- renderUI({
    withMathJax(helpText(paste0(
      '設定された効用最大化問題:',

      '$$ \\max_{x_1,x_2} \\hspace{1em}  U(x_1,x_2)=',input$a,'x_1','^\\frac{',input$c,'}','{',input$d,'}', input$b,'x_2','^\\frac{',input$e,'}','{',input$f,'}$$',

      '$$s.t. \\hspace{1em} ', input$p_1,'x_1 + ', input$p_2,'x_2 \\leq',input$budget,'$$'
      )))
  })


  #最適消費計画の解
  output$maxutility <- renderUI({

    ObjFunc = function(x) return( - input$a * x[1]^(input$c/input$d) * input$b * x[2]^(input$e/input$f) ) #目的関数
    #solnp()は最小化をするようになっているので，目的関数にマイナスを掛ける
    ConstFunc = function(x) return( x[1]* input$p_1 + x[2] * input$p_2 ) #制約式の右辺
    eq.value <- c(input$budget) #制約式の左辺
    x0 <- c(1,1) #決定変数を初期化

    solution <- solnp(x0, fun = ObjFunc, eqfun = ConstFunc, eqB = eq.value )

    #return
    withMathJax(helpText(paste0("効用を最大化する消費量:",
                                '$$x_{1}^*=',solution$pars[1],', \\hspace{1em}  x_{2}^*=',solution$pars[2],'$$')))
  })


  #無差別曲線と予算制約線と最適消費点
  output$indifferencecurve <- renderUI({
    withMathJax(helpText('無差別曲線と予算制約線と最適消費点'))
  })
  output$Plot2D <- renderPlot({
    #以下のreactiveな要素をrenderPlotに入れれば使えるっぽい
    #効用関数をfunctionの形で定義
    u <- function(x_1,x_2) {input$a * x_1 ^ (input$c/input$d) * input$b * x_2 ^ (input$e/input$f)} #効用関数を定義

    #xの範囲指定
    x_1 <- 1:input$range
    x_2 <- 1:input$range
    U <- outer(x_1, x_2, u) #outer()はx_1,x_2に対応したf(x_1,x_2)の値を行列で返す

    #plot
    par(mar=c(3,3,1,1))
    contour(x_1, x_2, U, method = "edge", labcex = 1,lwd = 2)
    abline(a = input$budget/input$p_2, b = -input$p_1/input$p_2, lwd = 2, col = "blue") #予算制約線 x_2 = ...

    #最適消費点
    ObjFunc = function(x) return( - input$a * x[1]^(input$c/input$d) * input$b * x[2]^(input$e/input$f) ) #目的関数
    #solnp()は最小化をするようになっているので，目的関数にマイナスを掛ける
    ConstFunc = function(x) return( x[1]* input$p_1 + x[2] * input$p_2 ) #制約式の右辺
    eq.value <- c(input$budget) #制約式の左辺
    x0 <- c(1,1) #決定変数を初期化
    solution <- solnp(x0, fun = ObjFunc, eqfun = ConstFunc, eqB = eq.value )
    points(x = solution$pars[1], y = solution$pars[2], lwd = 3, col = "darkblue", pch = 16) #最適消費点

  })  

  #3D plot of utility function 効用関数の3次元図
  output$Plot3Dtext <- renderUI({
    withMathJax(helpText('効用関数の3次元図'))
  })
  output$Plot3D <- renderPlot({

    #以下のreactiveな要素をrenderPlotに入れれば使えるっぽい
    #効用関数をfunctionの形で定義
    u <- function(x_1,x_2) {input$a * x_1 ^ (input$c/input$d) * input$b * x_2 ^ (input$e/input$f)} #効用関数を定義

    #xの範囲指定
    x_1 <- 1:input$range
    x_2 <- 1:input$range
    U <- outer(x_1, x_2, u) #outer()はx_1,x_2に対応したf(x_1,x_2)の値を行列で返す    
    persp(x_1, x_2, U,
          theta = 30, # 横回転の角度。ここをいじって表示のアングルを変える。
          phi = 20, # 縦回転の角度。ここをいじって表示のアングルを変える。
          ticktype = "simple", # 線の種類
          lwd = 0.5, # 線が太いと見づらかったりするので0.5にしておいた。
          col = F, #塗りつぶしなし
          border = 8 #枠線灰色
          )
  })  
})
8
6
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
8
6