2017-06-08

長方形に円を並べるという高校生の問題が難しかった

以下のツイートを見たことからはじまった。高校生の問題のはずなのに中々解けないというか、どのように手をつけて良いかわからず少しの間はまってしまった。わからなかったので、勉強会で知り合いに聞いたりして、最後にはほぼ答えを教えてもらってなんとかとけました。よかった。。

この問題って有名な問題でしょうか？
「偏差値72の高校の試験問題」だそうです pic.twitter.com/aYyy73Sel6
— アラクー (@Arakur65536) 2017年5月31日

問題
試行錯誤
- ３つずつ並べてみたり
- 対角線に一列に並べてみたり
回答
pythonのコード
twitter

問題

一辺の長さがそれぞれ4,2000である長方形に敷き詰めることができる単位円（半径が１の円）の個数は2011個以上であることを証明せよ。ただし、敷き詰める円は互いに重なり合わないものとする。必要があれば1.9819< $\sqrt{4\sqrt{3}-3}$ <1.982を用いても良い。

試行錯誤

という問題で、普通に格子状に綺麗にならべたら2000個は並べることができるってのはすぐわかるんだけど、これよりももっと隙間なく並べることができるかっている試行錯誤がはじまるわけ。
f:id:hideisu:20170605230505p:plain

３つずつ並べてみたり

３つの正三角形を組みにしてならべてはどうかと思うけど、2000個のやつよりもならばないし
f:id:hideisu:20170605231751p:plain

対角線に一列に並べてみたり

もしかしたら対角線にならべて、その上と下にうまいことならべたら綺麗にならぶのかとおもったんだけど、1200個くらいしか並びません。
そんなわけで、最終手段として、少しずつ上下にずらしながら、ジグザグに並べたらどうなのかと思い浮かびます。

回答

f:id:hideisu:20170605232207p:plain
こんな風にならべると、例えば下の行を見ると５の円を少しだけ上に並べることができるので、上の2000個の並べ方よりも少しだけ詰めて並べることができます。これはすごい。
直線に並べるよりも、どれだけ詰めてならべることができるか（=2-αとおく）は、上の図の５の球を少し上におくことでどれだけ詰めることができるかに帰着できるから

３と５の円の中心の位置を計算するために３と５を抜き出してして図にすると以下のように書くことができ
f:id:hideisu:20170608024956p:plain
αを三平方の定理を使って求めると $\sqrt{4\sqrt{3}-3}$ となります。

三平方の定理を使った計算はこんな感じ。
2² = α² + (2- $\sqrt{3}$ )²]
α² = 4 - (4 - 4 $\sqrt{3}$ + 3)
α² = 4 $\sqrt{3}$ - 3
α = $\sqrt{4\sqrt{3}-3}$

これで、問題文で与えられていた「必要があれば1.9819< $\sqrt{4\sqrt{3}-3}$ <1.982を用いても良い。」という条件と結びつくんですね。
一直線に並べると、円ある円の中心から隣の円の中心までの距離は２となりますが、５番の円を少しだけ上におくことで、隣の円とのx軸方向の距離は約1.982となり0.018くらいつめて置けるわけです。

このパターンで何個並べることができるか計算すると、2011個となります。並べた画像は以下。
f:id:hideisu:20170605232207p:plain
f:id:hideisu:20170608130349p:plain

回答としては、以上となります。
ここからは、これまでの画像を作るコードとtwitterでのつぶやきのコピー

pythonのコード

円を描くコード

f:id:hideisu:20170608213210p:plain

0から60を描くコード

f:id:hideisu:20170608213439p:plain

1999から2011までを描くコード

f:id:hideisu:20170608130349p:plain

twitter

f:id:hideisu:20170608220300j:plain

やばい、高校の問題とけない笑、どうやって解こうかなー https://t.co/ClwcWUQIO1
— カレーちゃん (@currypurin) 2017年5月31日

RTの問題、こんな感じにまず対角線に並べて、その上と下に並べるのが答えかと思ったんだけど違うみたいだ。並ばない(笑 pic.twitter.com/JyOwKjXlVA
— カレーちゃん (@currypurin) 2017年6月3日

となると、微妙に上下にずらしてジグザグに並べるのが答えか。。
— カレーちゃん (@currypurin) 2017年6月3日

横長に長方形を置いて、左端から下辺に接するように上１：下２の円を並べる。その次は上辺に接するように上２：下２の円を並べる。その繰り返しをすると、
3*335+1=2011個置けますね。
— nob (@nob_ymkw) 2017年6月4日

あー、なるほどそういう置き方ですか！！ようやく置き方がわかりました。感謝です。
— カレーちゃん (@currypurin) 2017年6月4日

例の問題、とりあえず60個並べてみた。
PCの性能が悪いのか、pythonのコードが悪いのか500個以上円を描くと処理落ちする。
実際に2000こ並べることに意味はないけど、円に何個目の円かラベルをつけて並べることができたら、解答としてはわかりやすいかも。 pic.twitter.com/Rl8tmH8quk
— カレーちゃん (@currypurin) 2017年6月4日

この問題初見で、短時間で解ける人はなかなかいないとおもう。自分のばあいは1人でやってたら一生とけなかったかも笑
最密充填て分野に分類されるのかしら？ https://t.co/U0fdF2k1QW
— カレーちゃん (@currypurin) 2017年6月4日

最密充填のなかでも、球充填て分野か。https://t.co/iUjeFKyhvf https://t.co/lkbKZ4gl76
— カレーちゃん (@currypurin) 2017年6月4日

最密充填って、そういえば高校の化学で少しやりましたね。
昔すぎて忘れていたわ。でも糸口があれば少しづつ思い出せるんだから、昔の記憶ってすごい。
— カレーちゃん (@currypurin) 2017年6月4日

一気に2000個の円を描く必要はなかったんだな。必要なところだけ描画する方法で、2011個並べることができました。
この方法で2012個目ははみ出してしまって、並べることはできなかったけど、並べることができないっていう証明はできるのかな？ pic.twitter.com/Pz5L2PB8Bn
— カレーちゃん (@currypurin) 2017年6月5日

2011以上である証明に加えて、2012以外である証明もあるらしい。前者は出来たが、後者はどうするのだろう？
— nob (@nob_ymkw) 2017年6月5日

え、まだこの問題の続きがあるの。。。2012以下である証明か面白そう。

2017-03-19

Rでサンプルデータを作成するのが難しかったので、作ってみた

一昨日の勉強会で、参加者全員でモデルに当てはまるサンプルデータをRで作ってみようとなった時になかなか良いサンプルを作れなかったので、落ち着いて作ってみる。 data-refinement.connpass.com

　勉強会の内容と作成するモデル

勉強会の内容は、発表者のSSAS3さんが書いてくれたこれ。ランダム化実験が驚くほどよくわかってすごかった。 qiita.com

作りたいサンプルは以下の条件を満たすことだけで、あとは適当でOKという条件（少し勉強会のときから条件を変えています）

	Y1(任意)	Y2(ランダム)
作業手伝いなし	a	b
作業手伝いあり	A	B

a,b,A,Bはテストの点数の平均。

a<A
- 任意にやらせると、準備が整っている生徒が作業をするので、作業ありの方が点数が高くなる
b>B
- ランダムで作業をやらせると、作業なしグループは準備ができるので点数が高くなる

使う関数の候補

関数	概要
rnorm	正規分布に従う乱数を発生
runif	一様分布に従う乱数を発生
sumple	無作為抽出
round	少数点以下を丸める
floor	少数点以下を繰り下げ
ceiling	小数点以下を繰り上げ

練習

n <- 10 # クラスの人数を10人とする
pre <- runif(n,min=-1,max=1) # 一様分布から乱数を作成
ceiling(pre) # preを切り上げて、0と１にする
sample(c(1,0),size=n,replace=TRUE) #ランダムでトリートメントグループとコントロールグループに分ける 
sample(rep(c(1,0),times=n/2),size=n) # 半分ずつトリートメントグループとコントロールグループに分ける
sample(rep(c(1,0),times=c(10,30)),size=n) #10人、30人のトリートメントグループとコントロールグループに分ける

サンプルデータ作成

library(dplyr)

set.seed(2)

NO <- 1:40 # 出席番号
pre <- runif(40,-1,1) #  準備ができているか
Tvol <- ceiling(pre) #  準備ができていれば１で作業を手伝う、できていなければ０で作業を手伝わない
y1 <- round(pre*20+70) +(1-Tvol)*sample(1:5,40,replace=TRUE) # 1回目の点数

Trand <-sample(c(0,1),40,replace = TRUE) # ランダムに作業手伝いを決める
y2 <- round(pre*20+70)+(1-Trand)*sample(1:5,40,replace=TRUE) # 2回目の点数

sample1 <- data.frame(出席番号=NO,Tvol,y1,Trand,y2) # データフレームを作成


# 1回目の小テストをクロス集計
sample1 %>%
group_by(Tvol) %>%
summarise(n=length(y1),mean=mean(y1))


# 2回目のテストをクロス集計
sample1 %>%
 group_by(Trand) %>%
 summarise(n=length(y2),mean=mean(y2))

jupyterのファイルを下記にあげています。他にもやり方が色々あるだろうから、コメントとか次の勉強会で教えてください。

github.com

2017-03-05

緑本3.4.2（あてはめとあてはまりの良さ）のメモ

緑本を読んでいて、3.4.2のポアソン回帰の最大対数尤度のところで、つまづいたのでメモ。

ある個体iの平均種子数 $λ_i$ が

${ λ_i = exp (β_1+β_2 x_i) }$

であると仮定

リンク関数によって

${log λ_i =β_1+β_2 x_i}$

このモデルの対数尤度は

$\displaystyle{logL(β_1 , β_2) = \sum_i log \frac{{λ_i}^y_i exp(-λ_i)}{y_i !}}$

本のデータで最大対数尤度を計算すると-235.4ぐらいとわかります。
ってところが複雑な式だなぁ、実際にどの数字を使って計算すれば計算できるのかなっていうんで、本のデータで計算してみた。

# データの読み込み
>d <- read.csv("http://hosho.ees.hokudai.ac.jp/~kubo/stat/iwanamibook/fig/poisson/data3a.csv")  

# GLMへの当てはめ
>fit<- glm(y~x,data=d,family=poisson)
Coefficients:
(Intercept)            x  
   1.29172      0.07566  

# 最大対数尤度を計算
> logLik(fit)
'log Lik.' -235.3863 (df=2)

# 最大対数尤度を上記の式で計算してみる
>logL <- function(x,y) log({exp(1.29172+0.07566*x)^y}*exp(-(exp(1.29172+0.07566*x)))/factorial(y))
>sum(mapply(logL,d$x,d$y))
[1] -235.3863

ポアソン回帰で切片と傾きの最尤推定値が得られた後に、最大対数尤度を念のため計算してみて、logLik関数の結果と一致したという結果でした。

データ解析のための統計モデリング入門――一般化線形モデル・階層ベイズモデル・MCMC (確率と情報の科学)

作者: 久保拓弥
出版社/メーカー: 岩波書店
発売日: 2012/05/19
メディア: 単行本
購入: 16人クリック: 163回
この商品を含むブログ (29件) を見る

2017-01-24

2017年1月から3月の統計学学習計画

きっかけ

Niseiさんの以下のつぶやきに触発されて、勉強の計画と結果を書き記していくことにしました。とりあえず読みたい本が決まっている３月までの予定を書きます。４月頭にはどれくらい進んだかちゃんと報告したい。

久々にブログ更新 - 機械学習に本気で取り組むためにやった数学周り前半戦結果 | きのこる庭ブログ https://t.co/r4Hc8c7C5g
— Nisei (@irration) 2017年1月22日

きっかけ
最終的な目標
１月から３月の計画
反省

最終的な目標

統計とか機械学習を数学的なところから理解して、今後もずっと学び続けたい
ガチで競馬の予想をやりたい

１月から３月の計画

１ヶ月間に１冊の本を内容を理解した上で読了することが目標。数学的にわからないところがあれば、なるべく統計学や解析学等の本を読んで理解を目指す。でも、完璧を目指すと時間がかかりすぎるので、どうしてもわからない箇所はメモって飛ばして、Rで何をやっているのみ理解することとする方針。先に進めばわかるかもしれないし、後で数学的な部分のみ復習しても良いので。なお、１月分と２月分は途中まで読み進めているので、全く初めからではない。