Rプログラミングの小ネタ

本日はちょっとPython篇。なぜっていうと、Pythonしか使えないよーって人がいるから。こういうテーブルを加工する系はR（というかdplyr）のほうが断然使いやすいとは思いますが、まあバイリンガルになるのも良いでしょう。 例えば、欲しいのは市区町村別の人口や面積や平均年齢のデータ。で、e-Statに置いてあるのを発見するも（ 国勢調査 都道府県・市区町村別の主な結果 都道府県・市区町村別の主な結果 | ファイル | 統計データを探す | 政府統計の総合窓口 ）、ファイルがエクセルだわ、列名が自治体コードと都道府県名を連結した文字列だわで、なんとも使いづらい。ちゃちゃっと加工しちゃえばいいんでしょうが、全国のいろんなところで、いろんな人が同じような作業をするのが無駄なので、サンプルコードとそれで出力されたCSVファイルを置いておきます。 やっているのは、 "openpyxl"でエクセルファイルを読み込みデータフレームに変換 必要な行や列を取り出す 正規表現で列名を加工 CSVファイルに書き出し あたりです。 import os import openpyxl import pandas as pd os.chdir(r'(エクセルファイルを置いたパス)') # ワークブックを開いた後、ワークシートを指定する # 読み込んでいるのは、e-Statの国勢調査データ wb = openpyxl.load_workbook('major_results_2020.xlsx') ws = wb['第１面事項_2020年'] df = pd.DataFrame(ws.values) # データが入っているのは9行目以降 # 0, 1, 4, 10, 11, 12列にはそれぞれ # 都道府県、市区町村、人口、面積、平均年齢が入っている df = df.loc[9:, [0, 1, 4, 10, 12]] dic = {0:'都道府県', 1:'市区町村', 4:'人口', 10:'面積', 12:'平均年齢'} df = df.rename(columns=dic) # 列名をつける # アンダースコアの前の部分が、 ...

まずは、動作確認用のサンプルコードです。 # サンプルコードのデータフレーム age    <- c( 20,  30,  40) height <- c(170, 168, 175) wait   <- c( 67,  64,  70) df <- data.frame(age, height, wait) df   age height wait 1  20    170   67 2  30    168   64 3  40    175   70 ちょ、待てよ。体重のweightがwaitになってんじゃん。 ってことで、この列名だけ変更したいと。 列名全体のベクトルをまとめて指定して変更する方法↓ names(df) <- c("age", "height", "weight") とか、 インデックス番号を使って、一部を変更する方法↓ names(df)[3] <- "weight" とかがあるんですが、列数がめちゃめちゃ多いデータだと数えるのが大変だし、列の削除や挿入に対してロバストでないし。できれば、列名の指定で変更を行いたい。 ちょっとコードがごちゃごちゃしますが、↓こんな感じでできます。 names(df)[ which( names(df)=="wait" ) ] <- "weight" 内側から順に見ていくと・・・ names(df)=="wait" で、列名が "wait" になっている位置が FALSE FALSE TRUE というベクトルで返ってきます。 それを which() に渡すと、「3」というインデックス番号が返ってきます。 なのでこの「3」を names(df)[] に対するインデックスの指定に使えばOKというわけです。

# サンプルデータの作成 fruits <- c("apple","durian","orange") count  <- c(2, 1, 3) df <- data.frame(fruits, count) df # 中身の確認   fruits count 1  apple     2 2 durian     1 3 orange     3 このデータを使って、ggplotで棒グラフを描いてみると、 library(ggplot2) ggplot(df, aes(x=fruits, y=count))   + geom_bar(stat="identity") 並び順はデータのまま 当然ながら、棒の順番はデータ通りに、２、１、３と並びますね。 これをソートしたい場合は、reorderを使います↓ ggplot(df, aes( x=reorder(fruits, count) , y=count)) +   geom_bar(stat="identity") reorderすると昇順に並ぶ x軸はfruitsなんだけど、reorder関数の第二引数であるcountの値で並べ替えてから使ってね、という指定です。 で、次の課題です。 reorderを使うと、昇順で１、２、３という並びになりましたが、これを降順の３、２，１という並びにしたいときは、どうやればいいか。 データフレームをソートするときなんかに使うorder関数でいうところの「decreasing=TRUE」みたいな指定があればいいのですが、reorder関数にはそのようなオプションはなさそうです。 で、実はごく簡単で、reorderの第二引数（count）の値の正負が逆になれば、順番も逆になるじゃんという理屈です。 ggplot(df, aes(x=reorder(fruits, -count ), y=count)) +   geom_bar(stat="identity") reorder関数の第二引数にマイナスをつければ降順に...

例えば、↓こんなデータがあったとして、 # サンプルデータの作成 name  <- c("Anne", "Bob", "Carl", "Dann", "Eric", "Fred") fruits <- c("orange", "apple", "orange", "durian", "orange", "apple") df <- data.frame(name, fruits) df # 中身を見てみる   name fruits 1 Anne orange 2  Bob  apple 3 Carl orange 4 Dann durian 5 Eric orange 6 Fred  apple みんなの好きな果物のデータだとして、１人しかいないような少数派の行は除去したいと。この例だと、ダンのドリアンを取り除きたいと。 ddplyを使って度数をカウントし、新たに度数（count）の列として追加（transform）。 library(plyr) ddply(df, "fruits", transform, count=length(fruits))   name fruits count 1  Bob  apple     2 2 Fred  apple     2 3 Dann durian     1 4 Anne orange     3 5 Carl orange     3 6 Eric orange     3 で、このcount列を条件として、データフレームをフィルタすればいいかなと。 dd <- ddply(df, "fruits", transform, count=length(fruits)) dd[ dd$count > 1, ] # 少数派の行を削除   name fruits count 1  Bob  apple ...

 ↓この本に載っていた例で、 Rによるデータサイエンス データ解析の基礎から最新手法まで library(MASS) caith        fair red medium dark black blue    326  38    241  110     3 light   688 116    584  188     4 medium  343  84    909  412    26 dark     98  48    403  681    85 ↑このようなデータ（縦に並んでいるのが目の色、横に並んでいるのが髪の色）に対して、↓こんな感じで対応分析を行う、という例が載っていました。 caith.ca <- corresp(caith, nf=4) biplot(caith.ca) Rでのcaithデータの対応分析結果 と、ここまでが前置き。 今、手元にあるのが↓こんなデータだった、としましょう。 gender <- c("M","M","M","M","F","F","F","F","F","F") blood  <- c("A","B","B","O","A","B","A","O","A","AB") gender.blood <- data.frame(gender, blood) gender.blood    gender blood 1       M     A 2       M     B 3       M     B 4       M     O 5  ...

例えば、Rで↓こんなデータを扱っているとします。 # サンプルデータを作る Q1 <- c(1, 1, 2, 2) Q2 <- c(1, 2, 1, 2) Q3 <- c(1, 2, 3, 1) df <- data.frame(Q1, Q2 ,Q3) # 中身の確認 df   Q1 Q2 Q3 1  1  1  1 2  1  2  2 3  2  1  3 4  2  2  1 読み込んだときの都合か何かで、データはinteger型とかnumeric型になっていると。 でも、実は質問 Q1、Q2、Q3に、1:はい、2：いいえ、3:どちらともいえない、とかで答えたもので、因子型として扱いたい。 多重対応分析のmcaとかを使おうとすると、 mca(df = df) でエラー: all variables must be factors が出ちゃうとか、そんなシチュエーション。 因子型に変換したいときには、as.factor関数ですが、これはデータフレームに対しては使えない。 こんな時は、あの一家。そう、applyファミリーの登場です。 lapplyを使って、1列ごとにas.factor関数を適用、リストとして返ってきたものを、またデータフレームに戻してやるという流れです。 df.fctr <- data.frame( lapply(df, as.factor) ) # 型の確認 df.fctr$Q1  [1] 1 1 2 2 Levels: 1 2 無事、因子型になりました。

「R ヒストグラム 一部 色をつける」で検索してみると、hist関数でヒストグラムを描いた後に、polygon関数で色をつける、なんて方法がヒットしました。 polygon使えばなんでもできそうだけど、なんか、ちょっと違うよなあ、とか思ってしまいまして。 で、実はhist関数のcolオプションでも、できるんですよね。 colオプションに1つの値（スカラー）を指定すると、全体が一色で塗りつぶされてしまいますが、ここにベクトルを指定すると、それぞれの棒の色を指定することができます。 例えば、ヒストグラムに10個のビンがあって、それぞれを任意の色で塗りたい場合は、10個の要素を持つベクトルをcolオプションに指定すればOKです。 set.seed(0)       # 再現性のために rd <- rnorm(100) # 100個の乱数 cols <- c("white", "white", "red"  , "white", "white",           "blue" , "white", "white", "white", "white") hist(rd, col=cols) ヒストグラムの一部に色をつける 色を塗りたくない場合は（パワポなどの「塗りつぶしなし」みたいな感じ）、色名の代わりにNAを指定すればいいです。 cols <- c(NA    , NA, "red", NA, NA,           "blue", NA, NA   , NA, NA) hist(rd, col=cols) 最初の例と全く同じ見た目になると思いますが、add=T指定で重ねたときなんかに差がでますね。 階級がいっぱいあって、いちいち全部書き出すのが面倒なときは、下記のような感じで、塗りたいところだけを指定すればいいですね。 cols <- rep("white", 100) # ”white”を詰めた、長めのベクトルを作って...

↓以前、こちらの記事を書いたのですが、 「データの見えざる手」のU分布を、Rでシミュレート よくよく見てみると、書籍と軸の取り方が違ったりして、つっこみどころ満載だったので、悔い改めて、ちゃんとやることにしました。 書籍では、横軸がマスに入っている個数、縦軸が累積確率になっていました。 あと、初期値もちゃんとランダムで設定するようにしました。 それと、対数プロットする際に、軸の目盛ラベルの付けやすさから、ggplotを使ってみました。 library(ggplot2) library(scales) n <- 72000 # 点の個数 m <- 900   # マスの個数 masu <- numeric(m) # 空のマスを用意 # 点をランダムにマスに配置 indices <- sample(1:900, 72000, replace=T) for(i in indices){   masu[i] <- masu[i] + 1 } for( i in 1:100000000 ) {   s <- sample(1:m, 2) # ランダムに2つのマスを選ぶ   if( masu[s[1]] != 0 ) { # 無い袖は振れないケースへの対処     masu[s[1]] <- masu[s[1]] - 1 # 1つ目のマスから取って、     masu[s[2]] <- masu[s[2]] + 1 # 2つ目のマスへ入れる   } } tbl <- table(masu)     # 個数を集計 df <- data.frame(tbl) # データフレームにする # 列名を分かりやすくする colnames(df) <- c("num_of_dots", "freq") # 点の数が因子型なので、整数型に変えておく df$num_of_dots <- as.integer(df$num_of_dots) # 累積確率を計算 df$cum_prob <- rev(cumsum(rev(df$f...

さて、本筋とは関係のないところでケチばっかりつけていた、↓前回と前々回の記事でしたが、 「データの見えざる手」の正規分布の図が一様分布に見えたのでRで試した 「データの見えざる手」の図が分かりにくかったので、Rで一次元プロット 今回は、著者の矢野和男さんの言うところの「U分布」なるものを、計算機シミュレーションで作り出してみましょう。 次に、このようにランダムに玉を分配した後で、マス目間で玉をやりとりさせてみよう。 ランダムにマス目を二つ選んで、一方から他方に玉を１個移す。そして、これを繰り返してみよう。もともと、ランダムに置いた玉なのだから、そこからランダムにマス目を選んで、玉を動かしても、結果は変わらない、と思うだろう。この問題を多くの人に出題してみたが、全員が「結果は変わらない」と答えた。 たしかに、直感的には、ランダムに配置後にランダムに交換しても、マクロな状況は変わらないような気がしますね。でも、そうじゃないのが興味深いところ。 書籍では初期値はランダムとありましたが、手を抜いて、1マス80個の「平等」状態からスタートさせてみました。（結果は同じになりますよね、たぶん） 指定回数の交換を行ったあと、それぞれのマスが持っている個数でソートして、少ない方が左になるようにプロットしています。 1万回、2万回、・・・と実行しながら、プロット結果を画像として出力していきます。jのところのループ回数を変えて、1億回まで実行してみました。 n <- 72000   # 点の個数 m <- 900     # マスの個数 masu <- rep(n/m, m) # 平等に配分 for( i in 1:9 ) {   for( j in 1:10000 ) {     s <- sample(1:m, 2) # ランダムに2つのマスを選ぶ     if( masu[s[1]] != 0 ) { # 無い袖は振れないケースへの対処       masu[s[1]] <- masu[s[1]] - 1 # 1つ目のマスから取って、       masu[s...

↓この記事を書いていて思ったのですが、 「データの見えざる手」の正規分布の図が一様分布に見えたのでRで試した 一日の生活の900分は一次元的であるのに、それを30×30の二次元のマスで表現しているところが、そもそも分かりにくい。 人が理解するときのモデルとしても分かりにくいし、二次元になっているせいで、シミュレーションのスクリプトを書くときもいろいろと面倒な処理が必要になったりして（上記のリンクのinteractionのくだりとか）。 素直に、一次元的な図を載せた方が、読者の理解も進むのではと思って、Rで書いてみました。（点やマスの数は見た目がほどよくなるように減らしてあります） n <- 80 # 点の数 m <- 10 # マスの数 x <- runif(n, min=0, max=m) stripchart(x, pch=1, xlab="１分ごとのマス") abline(v=0:m) legend("topright", legend="手の動きのあった時点", pch=1, bg="white")  こういう図の方が分かりやすいと思うけどなあ。

「 データの見えざる手 」は読んでいて、引っ掛かりまくりでした。 読んだかた、「U分布」ってピンときました？ どこで躓いたかというと、こんな感じでコンピュータシミュレーションの結果が紹介されていたのですが、上の「正規分布（ポアソン分布）」と書かれている図↓の見た目が一様分布っぽいんですよね。 正規分布（ポアソン分布）とU分布　（「データの見えざる手」より） 本文を見てみると、 コンピュータシミュレーションでこれを実行するには、玉の位置をランダムに生成すればよい。横方向の位置(x)を決める1～30の乱数と縦方向の場所(y)を決める1～30の乱数を発生させ、(x,y)の位置に玉を置くのだ。 とあります。う～ん、xとyを一様分布に従って発生させて2次元にプロットしている、ってだけだよなー。 ・・・（考え中）・・・ で、しばらく考えてみて、やっと分かりました。こうやって発生させたデータだと、マスの中の点の個数が正規分布になるのね。 実際に、一つずつやってみましょう。 n <- 1000 # 点の個数 m <- 10   # マスの区切りの数 # 座標は一様分布で従って発生させ、プロット x <- runif(n, min=0, max=m) y <- runif(n, min=0, max=m) plot(x, y, pch="・") # マスを書く abline(h=0:m) abline(v=0:m) 一様分布で発生させた点（マス内の個数は正規分布になる） はい、本に載っているのと似たような絵になりました。 マスに入る個数を調べるには、引数を越えない整数を返すceilingが使えますね。 > interaction( ceiling(x), ceiling(y), sep="," )    [1] 4,1   9,3   10,1  1,5   3,1   1,6   8,1   8,3   10,6  5,7   2,6   [12] 2,5   2,7   7,6   6,10 ...

いつもは、↓別のブログの方に載せていたのですが、 「エレガントな問題解決」演習問題 2.1.25の解答（分母が3つの項の積になっている数列の和）: 主張 今回はRを使うので、こちらのブログで。 「エレガントな問題解決」 P29 2.1.27（ひっかけ問題）(c) 偏りのないコインを投げて、表が3回以上出る確率が50%を超えるには、最低何回投げればよいだろうか。 自分が題意を理解できているのか、いまいち自信がありません。私は、ひっかかっているのか？ とにかく解答してみます。 コインを1回投げて、表が3回以上出ることは、もちろんありません。 なので、考慮に値するのは3回以上投げるケースから。 3回投げて、表が3回出る確率をRで求めてみると、 > (1/2)^3 [1] 0.125 なので、12.5%。これは50%を超えていませんね。 では、4回投げる場合。表が3回出る確率と、表が4回出る確率を足せばいいですね。 4回投げて、表が3回出る確率は、（Cは組み合わせの記号だと思ってください） (1/2)^4 × 4 C 3 ですね。Rでは、組み合わせ数を求めるのにchooseコマンドが使えるので、 > (1/2)^4 * choose(4,3) [1] 0.25 4回投げて、表が3回出る確率は、 > (1/2)^4 * choose(4,4) [1] 0.0625 上記の2つを足しても、まだ50%は超えませんね。 いちいち足すのもまどろっこしいので、関数にしてみましょう。 # # n回投げたときに、3回以上でる確率を返す関数 # ProbMoreThan3 <- function(n){   s <- 0    # 組み合わせ数の合計値格納用   for(i in 3:n){     s <- s + choose(n,i)   }   p <- 1/2^n * s   return(p) } 試しに呼んでみましょう。 > ProbMoreThan3(3) [1] 0.125 > ProbMoreThan3(4) ...

最近知ったんですが、Rで順列（Permutation）や、組み合わせ（Combination）を計算するデフォルトの関数ってないんですね。 【2016年11月18日追記 開始】 勘違いしていました。訂正します。 組み合せの数を計算するchooseという関数があります。 例えば、4個から2個選ぶ組み合わせの数は > choose(4,2) [1] 6 みたいな感じで求めることができます。 【2016年11月18日追記 終了】 そういうパッケージはあるみたいですけど、わざわざパッケージを使うまでもないような気もしないでもないです。 高校数学を復習する意味でも、自分で実装してみましょう。 まず、順列。 n個の中からk個を取り出して、順番を意識して並べていくような場合の数ですね。 nから初めて、1ずつ減らしながら、r個分かければ計算できます。     n × (n-1) × (n-2) × ... × {n-(r-2)} × {n-(r-1)} 上記をやりたい場合、受け取った引数を全部かけてくれるprodという関数が使えます。１ずつ減っていく引数は、 ":" （コロン）を使えばよさそうですね。 ということで、     prod(n:(n-r+1)) とやればいいのですが、これだと r が 0 のときにうまくいきません。0個を選び出す順列の数は 1 なので、関数にしてしまって、↓こんな感じで、どうでしょうか？ # 順列を計算する関数定義 Permu <- function(n, r){   ifelse( r==0, 1, prod(n:(n-r+1)) ) } 動作確認してみましょう。 > Permu(4, 0) [1] 1 > Permu(4, 1) [1] 4 > Permu(4, 2) [1] 12 > Permu(4, 3) [1] 24 > Permu(4, 4) [1] 24 うん、よさそうです。 次は、組み合わせ。 nからr個選ぶ順列を、rの階乗で割れば（r個の順番を無視する）、OKですね。 階乗はfactorialという関数があるので、そのまま使えます。 ...

例えば↓こんな感じで、各市町村の人口データがあったとします。 city   <- c("A市", "B町", "C村", "D市") male   <- c(95, 49, 26, 45) female <- c(99, 51, 24, 55) df <- data.frame(city, male, female) df   city male female 1  A市   95     99 2  B町   49     51 3  C村   26     24 4  D市   45     55 で、人口や男女比を比較するために、それぞれの市町村に対してグラフを描いてみようとしたとします。 じゃあ、１行ごとにループで回しながら、barplotあたりを使って、グラフを描いてみようとしまして、↓こんな風にスクリプトを書くと（iがループのカウンタだと思ってね）、 # i行目の男女人口を棒グラフにしたい・・・ barplot(df[i, c("male", "female")])  barplot.default(df[i, c("male", "female")]) でエラー:    'height' はベクトルか行列でなければなりません なんでだ？と思って、切り出される部分をチェックしてみると、 df[i, c("male", "female")]   male female 1   95     99 class(df[i, c("male", "female")]) [1] "data.frame" 構成要素としては、２つの数値の並びになっているもの、型がデータフレームなもんだから、barplotにそのまま渡すとエラーになってしまうのね。ちなみに、円グラフのpie...

福岡県の市単位（政令指定都市は区単位）の人口密度単位のランキングは、↓こちらにあって、 福岡県の人口密度番付 - 都道府県・市区町村ランキング【日本・地域番付】 　　1位 福岡市 中央区（11770人/㎢） 　　2位 福岡市 城南区（ 8031人/㎢） 　　3位 福岡市 南区　（ 7976人/㎢） となっております。トップは予想通りの福岡市中央区。 ただ、中央区と言ってもそれなりに広いし、人口密度にはムラがあるでしょうし、もしかしたら他の区にもギュギュっと人の密集しているエリアがあるかもしれません。 ということで、Rを使って町丁目単位（博多駅前１丁目とか、天神２丁目とか）でのランキングを出してみました。 とりあえず結果をどうぞ↓ 中央区全体の人口密度は 11770人/㎢ でしたが、中央区荒戸２丁目で人口密度を計算すると 49294人/㎢ と、4倍以上の数字を叩き出しました。 また、上位5位の荒戸、平尾、薬院、高砂はいずれも中央区の町名ですが、6位に入っている愛宕浜は西区にある町です。 西区と言えば福岡市の中では辺境ですし（注：著者の個人的なイメージです）、前述の市区単位のランキングでは17位とかなり下位にランキングされてました、 実は、この愛宕浜2丁目は室見川をはさんで早良区のすぐ対岸にあって、比較的都市部に近い場所なんですよね。そして、ほとんどマンションしかありません（他には中学校とマルキョウがある）。大抵のマンションは10階以上あるので、おのずと人口密度が高くなるわけですね。 ↓Google Earthで見てみるとこんな感じ では、最後にRスクリプトを載せておきます。 シェープファイルは↓こちらに書いてある方法で入手して、 e-Stat(統計局)で公開されているShapeファイルを、Rで表示する - Rプログラミングの小ネタ ↓このやり方でマージしました。 Rで複数のシェープファイルを結合する - Rプログラミングの小ネタ library ( maptools ) library ( ggplot2 )   # シェープファイルを読み込む pj <- CRS ( "+proj=longlat +datum=WGS84" ) shp <-r...

地図で可視化する際の行政界のシェープは、統計局のものをよく利用させてもらっています。↓以下参考。 e-Stat(統計局)で公開されているShapeファイルを、Rで表示する - Rプログラミングの小ネタ 統計局で公開されているシェープは比較的小さい単位に分かれています。例えば、福岡市だと、区単位（東区、博多区、中央区、南区、西区、城南区、早良区）の7つのシェープに分かれています。 こういうのをまとめて１つにしたい時ってありますよね。 QGISのようなGUIのあるGISソフトでやる方法もありますが、数が多くなるとしんどい。ということで、Rのスクリプトでシェープの結合をする方法を紹介します。 ざっくり言うと、readShapePolyで読み込んで、spRbindで結合させるだけなんですけどね。 でも、そう簡単にはいかない時があるんですよね。というのも、ポリゴンのIDがかぶっているとspRbindでエラーが出てしまう。例えば、東区のシェープの中のあるポリゴンに「1」というIDが付いていて、博多区のシェープの中のあるポリゴンにも「1」というIDが付いていると、この２つを結合させようとする際にエラーになってしまいます。 なので、東区のポリゴンはのIDは、"1.1", "1.2", "1.3", ...、博多区のポリゴンのIDは"2.1", "2.2", "2.3", ... みたいにリネームしたあとに、結合させましょうというのが、↓このスクリプトのポイントでございます。 library ( maptools )   # シェープのファイル名を取得 shape_files <- list.files ( path= "shape_folder" , pattern= ".* \\ .shp" , full.names=T )   # 緯度経度、世界測地系 pj <- CRS ( "+proj=longlat +datum=WGS84" )   # ポリゴンのIDが重複しないように前処理したあと、 # リストに格納しておく shape_li...

ベクトル場がxy座標の式で与えられているときに、それがどんな場を表しているかは、具体的にいくつかのベクトルを図示してみると、なんとなく見えてきたりします。 　　f(x, y) = (x, y) みたいな単純なベクトル場でさえ、慣れてないとイメージわかなかったりしますよね。 Rで図示してみましょう。（唐突） # # ベクトル場を表す関数 # f <- function ( x , y ) { X <- x Y <- y return ( c ( X , Y ) ) }   # # 矢印を描画する関数 # write_arrow <- function ( x , y ) { a <- 0.1 # 矢の長さ（比率） b <- 0.05 # 矢尻の長さ arrow_head <- f ( x , y ) * a # aの値で長さを調整   arrows ( x , # 矢のx座標(from) y , # 矢のy座標(from) x + arrow_head [ 1 ] , # 矢のx座標(to) y + arrow_head [ 2 ] , # 矢のy座標(to) length = b ) # 矢尻のサイズ }   # # ここからメイン # n <- 10 # 格子の１辺の数   # 描画する領域を準備 plot ( 0 , 0 , xlim= c ( -n , n ) , ylim= c ( -n , n ) , type= "n" , xlab= "x" , ylab= "y" )   # 各格子点について繰り返す for ( x in -n:n ) { for ( y in -n:n ) { write_arrow ( x , y ) } } f(x, y) = (x, y)のベクトル場 「ぶ...