Pythonによるハンバーガー統計学 #2

引き続きハンバーガー統計学にようこそ!の第2章。 準備 import numpy as np import pandas as pd import scipy as sp 2.1 平均的ポテトを推定する wakupote = pd.Series([47, 51, 49, 50, 49, 46, 51, 48, 52, 49]) wakupote.describe() 2.3 信頼区間/区間…

Pythonによるハンバーガー統計学 #1

ハンバーガー統計学にようこそ!の第1章をPythonでやってみたもの。 準備 # -*- coding:utf-8 -*- import pandas as pd import numpy as np import matplotlib.pyplot as plt 1.1 ポテトの長さの平均は # ワクワクバーガー waku = pd.Series( [3.5, 4.2, 4.9…

Rによるハンバーガー統計学 #3

(引き続きハンバーガー統計学にようこそ!の第3章をRでやってみたものです。) 3.1 チキンの売上は少ないのか このようなタイプのデータを「とりあえずプロットする」方法として、モザイクプロットというものがある。 sales <- matrix(c(435, 265, 165, 135),…

Rによるハンバーガー統計学 #2

(引き続きハンバーガー統計学にようこそ!の第2章をRでやってみたものです。) 2.1 平均的ポテトを推定する 前回も説明したとおり、標準のvar()、sd()でそれぞれ不偏分散、標本標準偏差が計算される。 wakupote <- c(47, 51, 49, 50, 49, 46, 51, 48, 52, 49)…

Rによるハンバーガー統計学 #1

(ハンバーガー統計学にようこそ!をRでやってみたものです。) 1.1 ポテトの長さの平均は # ワクワクバーガーのポテトの長さデータ waku <- c(3.5, 4.2, 4.9, 4.6, 2.8, 5.6, 4.2, 4.9, 4.4, 3.7, 3.8, 4.0, 5.2, 3.9, 5.6, 5.3, 5.0, 4.7, 4.0, 3.1, 5.8, 3.…

なんでn-1なんですか

不偏分散はなぜnではなくn-1で割るのか、という問題。 以前似たようなことやったけどきちんと証明していなかったのと今ならどう書くか試してみたかったので。 証明をする 確率変数の実現値を、平均をとする。 まず、の期待値を求めよう。 任意の数を足して引…

ユーザーの平均継続期間が解約率の逆数であることを等比数列を使わずに理解する

ユーザーの平均継続期間は「解約率=一定」と仮定すれば「1/解約率」で示せる。これについてはユーザの平均継続期間が「1/解約率」で求められることの数学的証明 - it's an endless world.に等比数列を用いたわかりやすい証明があるが、「解約率=一定」とした…

Project Euler Problem 22

問 Problem 22 - PukiWiki 名前のリストが入ったテキストファイルを読み込み、ソートする。 アルファベットを1から26までの整数に対応付けてその和を計算し、さらに名前の順位を乗ずることで「名前のスコア」を算出する。 例えばCOLINなら3 + 15 + 12 + 9 + …

Project Euler Problem 21

問題の概要 Problem 21 - PukiWiki d(n)をnの「真の約数」の和とする。 「真の約数」はnの約数のうちnを除外したものと定義する。 d(a) == b かつ d(a) == bであるようなaとbをAmicable numbers(友愛数)と呼ぶ。 完全数(すなわちa == b)は除外する。 10000以…

スペインの施設園芸

スペインにはハウスがたくさんある 昔作った資料をもとに話す機会があったので、スペインの施設園芸について手持ちの文献を再度調べていた。 スペインの南部のアルメリア県にはハウスが集中している地域がいくつかあって、特にアルメリア市の西部に位置する…

プランクの法則を綺麗にプロットしたい

技術的に難しい部分とか特に無いのでこちらで。 プランクの法則 Wikipediaで確認しよう。ちなみに英語版のほうが詳しいぞ。 プランクの法則 - Wikipedia あまりきちんと説明するとボロがでるのでいい加減に説明すると、暖かいものほど強い光を出して、暖かい…

有機質資材からの窒素放出量を予測したい

有機質資材からの窒素放出 ほ場に有機質資材、例えば牛ふん堆肥であったり稲わらであったりというものを施用した場合、これらに含まれている窒素というのは必ずしもすぐには使いやすい形にならない。有機質資材中の窒素は基本的には微生物の分解を受けること…

Project Euler Problem 16:20

Problem 16: Power digit sum gmpで解いた。 Problem 17: Number letter counts どの数のときにどれだけの文字数がどこに追加されるかを考えて、1000個のベクトルにまとめて加算していった。 最初綴りを間違えていて解けなかった。恥ずかしい。 Problem 18: …

Project Euler Problem 11:15

Problem 11: Largest product in a grid 問題自体は難しくないし計算時間もかからない。文字列の読み込みと添字操作の練習。 Problem 12: Highly divisible triangular number gmpを使って無理矢理やっても2〜3秒だが、overviewにより効率的な方法が書いてあ…

RGRからゴンペルツ曲線を推定する 3

R

ゴンペルツ曲線では相対成長速度と重量の関係はどうなってるのかと思って適当に回帰してみたらが良く当てはまったので、この仮定(相対成長速度は重量の対数に比例して減少する)のもとに微分方程式を解いてみると、ちゃんとゴンペルツ関数が出て来る。ここで…

RGRからロジスティック曲線を推定する

R

植物の重量に上限があったとして、相対成長速度が上限と今の重量の差に比例するとしてみよう。 これは適当な定数をもってきて書き換えると、要するに相対成長速度が重量の一次式で表現できるということになる。 これを相対成長速度の定義に代入する。 そして…

RGRからゴンペルツ曲線を推定する 2

R

前回はRGRを重量の対数の変化速度として求めたが、そもそもRGRの定義はであったので、複数のデータポイントがあれば何らかの方法で平滑化して微分すれば、より精度良くRGRが推定できるだろう。例えば前後2点を含めた3点をラグランジュ補間して得た二次関数を…

RGRからゴンペルツ曲線を推定する

R

植物の重量について、瞬間的な成長速度は植物の重量に比例すると仮定する。このときを相対成長速度(RGR: Relative Growth Rate)と呼ぶ。もしが変化しないのであれば、植物は指数関数的に増大を続けてしまい、地球は程なくして滅ぶだろう。しかし、植物は無限…

Excelで表形式に入力されてしまったデータを整形するにはINDEX関数が便利

こういう「お前は何がしたかったんだ」みたいなファイルが引き継がれることありますよね?ない?そいつはよかったな!!!!! クソみたいなデータを整形する 普通データはこう入力するだろう。こうでなくては何も始められない。 だが不幸なことに現実はこう…

Project Euler Problem 1:10

久しぶりにProject Eulerをやってみようと思ったが、せっかくなので進捗をリセットして最初からにした。 とりあえず10問やってみたが大分脳がヤバくなってる感あったのでちゃんと継続してボケ防止に努めたい。 解答コードはgithubで管理することにして、ここ…

「データ解析のための統計モデリング入門」第11章メモ

p.254の図11.7をstanを使って再現したかった。stanコードはこんな感じ。 data { int<lower=1> N; // 欠測値を含むデータ数 int<lower=1> N_obs; // 欠測値以外のデータの数 int<lower=1, upper=N> Idx_obs[N_obs]; // 欠測値以外のデータのインデックス int<lower=0> Y_obs[N_obs]; // 欠測値を除いたデー</lower=0></lower=1,></lower=1></lower=1>…

「データ解析のための統計モデリング入門」第9章メモつづき

R

前回。 rion778.hatenablog.com これをrstanarmでやってみた。 準備 install.packages("rstanarm") 割と時間かかった。 library(rstan) library(rstanarm) rstan_options(auto_write = TRUE) options(mc.cores = parallel::detectCores()) # データダウンロ…

geom_barにおけるNAの扱い

geom_barというかstat_countの方だと思うけど、例えばこんなデータがあったとして data <- c(1, 2, 2, 3, 3, 3, NA, NA, NA, NA) 普通に棒グラフを描くとこうだ。 ggplot() + geom_bar(aes(data)) ところがfactorにするとこうだ。 ggplot() + geom_bar(aes(f…

「データ解析のための統計モデリング入門」第9章メモ

緑本9章の内容について、RStanでやってみる。 Stanコード 次の内容で作成し、sample.stanという名前で保存した。 data { int<lower=1> N; vector[N] x; int<lower=0> y[N]; real mean_x; } parameters { real beta1; real beta2; } model { for(i in 1:N){ y[i] ~ poisson(exp(</lower=0></lower=1>…

WindowsにRStanをインストール

基本的にはRStan Getting Started (Japanese) · stan-dev/rstan Wiki · GitHubに従っていけた。 環境 Windows 10 Pro R: 3.3.1 RStudio: 0.99.903 Rtoolsのインストール Building R for WindowsよりRtools34.exeをダウンロード。 セットアップ中、"Edit the …

カテゴリカルデータ解析読書メモ(第10章)

R

決定木 データ全体を説明変数を用いて段階的にグループ分けする手法。 決定木を書く テキストではmvpartパッケージを使用しているが、現在は公開が停止されているので代わりにpartykitパッケージを使用する。 rpart()の代わりに使う関数はctree()であるが、…

Rで逆推定(2)

R

はてなから1年前のブログとかいうメールが来てこんな記事を書いていたのを思い出した。 rion778.hatenablog.com 再び少し調べたところ、これはinvestrというパッケージで実行できるらしい。 データとモデルオブジェクトの作成 データは前回と同じものを使う…

カテゴリカルデータ解析読書メモ(第9章)

対応分析 カテゴリカルデータにおいて各カテゴリーの変数に適当な数字を割り当てると、2つのカテゴリ間で散布図を書いたり相関係数を計算することができる。このとき、相関係数が最大となるような数値をカテゴリー変数に割り当てることで、カテゴリー間の関…

カテゴリカルデータ解析読書メモ(第8章)

R

対数線形モデル 対数線形モデルにおいては、分割表における各セルの期待度数を予測する。すなわち、複数のカテゴリカル変数の中から特定の目的変数を設定するのでなく、全ての変数について相互の関係を調べることを目的としている。 lonlin()とloglm() 対数…

カテゴリカルデータ解析読書メモ(第7章)

R

ポアソン分布 適合度検定はvcdパッケージのgoodfit()関数で行う。 > t <- 0:5 > x <- c(8, 6, 8, 3, 4, 1) > x <- matrix(c(x, t), nr = 6) > library(vcd) > result <- goodfit(x, type = "poisson", method = "MinChisq") > result Observed and fitted va…

カテゴリカルデータ解析読書メモ(第6章)

R

ロジット変換 比率pを変換した値と説明変数xの間に線形関係を仮定する一般化線形モデルにおいて、比率pに対する変換として有名なものにロジット変換とプロビット変換がある。いずれもp→0でf(p)→-∞、p→1でf(p)→∞となるような変換であるが、ロジット変換の方が…

カテゴリカルデータ解析読書メモ(第5章)

シンプソンのパラドックス 層別をしないで解析すると関連が見られるが、層別をして解析すると関連が見られなくなるような現象。層別に用いるカテゴリカル変数が、他のカテゴリカル変数全てに影響を与えているような場合に発生する、見かけ上の相関。連続変数…

カテゴリカルデータ解析読書メモ(第4章)

R

オッズとオッズ比 2つの割合とに対して、オッズの比を考えると、このオッズ比が1より大なら、1より小ならである。 オッズ比は、ロジスティック回帰において回帰係数がオッズ比の対数に一致することや、患者対照研究のようにオッズ比のみが推定できる研究方法…

カテゴリカルデータ解析読書メモ(第3章)

R

第3章の内容は大体知ってる事なのでささっと。 二項検定 > binom.test(405, 765) Exact binomial test data: 405 and 765 number of successes = 405, number of trials = 765, p-value = 0.1116 alternative hypothesis: true probability of success is no…

カテゴリカルデータ解析読書メモ(第2章)

R

気付いたら買って本棚に飾ってあったので読んでいる。カテゴリカルデータ解析 (Rで学ぶデータサイエンス 1)作者: 藤井良宜,金明哲出版社/メーカー: 共立出版発売日: 2010/04/22メディア: 単行本 クリック: 13回この商品を含むブログ (3件) を見る 準備 insta…

Rで複数条件抽出&集計

R

このようなデータフレームがあって、条件c1とc2に基づいて何らかの集計をしたいとする。 また、データフレームはdfというオブジェクトに代入されているとする。 tapply > with(df, tapply(v, list(c1, c2), mean)) 1 2 a 1.689910 2.563432 b 1.780409 2.568…

Excelで複数条件抽出&集計

このようなデータがあって、条件1と条件2に基づいて何らかの集計をしたいとする。 AVERAGEIFS 平均値を計算したいのであれば、AVERAGEIFS関数があるので、例えばこのように入力する。 =AVERAGEIFS(平均対象範囲,条件範囲1,条件1,条件範囲2,条件2)最初の例で…

グループ変数に応じて複数の折れ線グラフを書く

R

緑色の本のp.119に載っている図と似たものを書こうとして、 d <- data.frame(q = rep(c("q0.1", "q0.3", "q0.8"), c(9, 9, 9)), p = c(dbinom(0:8, 8, 0.1), dbinom(0:8, 8, 0.3), dbinom(0:8, 8, 0.8)), y = rep(0:8, 3)) というようなデータを準備したもの…

RStudioのキーバインド変更

Windowsでのお話。 Tools -> Modify Keyboard Shortcuts...で出来るようになっているという話だけど、Ctrl+hをBackSpaceにしたいとかCtrl+zをPgUpにしたいとか思ってもどうもうまく出来ないので結局AutoHotkeyを使うことにした。ほんの少しゴニョゴニョした…

ggplot2で凡例のラベルと項目名を操作する

ラベルと項目名を操作する 以下の説明ではggplot2を含めて3つのパッケージを使う。 library(dplyr) library(tidyr) library(ggplot2) データは次のように準備した。 # テストデータの準備 testdata <- data.frame(x = seq(1, 10, 0.1)) %>% mutate(sin = sin…

facetのstripのtextを変更する

R

ggplot2でfacet_gridやfacet_wrapを使ってグラフを分割描画した際に、各グラフについているラベル(strip)のテキストは、何もしなければ分割に用いたグループ名が自動的に入る。 例: ggplot(iris, aes(x = Petal.Length, y = Petal.Width)) + geom_point()…

飽和水蒸気圧の求め方

何年か前に飽和水蒸気圧の求め方に関して、Goff-Gratchの式やTetensの式やMurrayの式について少し触れたことがあった*1のだが、最近Wikipediaを見ていたらTetensの式だと思っていたものが Tetens(1930)のパラメータ値によるAugust他の式 飽和水蒸気量 - Wiki…

ThinkPad X61 SSD化+その他いろいろ

2万ほど投資したらまだまだ使えそうな感じになった。 用意したもの SSD(A-DATA ASP900S3-256GM-C-7MM ADATA 2.5"SSD 256GB SATA6G A-DATA ASP900S3-256GM-C-7MM) 2.5インチHDケース(【日本正規代理店】 ORICO 2.5インチ HDD/SSD 外付け ドライブ ケース S…

ggplot2で任意の線分

R

baseライブラリ中のlines()に対応するggplot2のレイヤー関数はgeom_line()であると説明される場合が多く、大抵の場合はそれで事足りるのだが、任意の位置に線分を描き込みたいという場合はあまり使い勝手が良くない。 例えば、「統計学:Rを用いた入門書(Mi…

分布の集中度を調べる

R

ランダムか、集中か 水田の調査圃場において、正方形に区切った調査区画毎にニカメイガの卵塊が次のように見つかったとする。 egg <- c(2,4,1,1, 4,0,3,3, 4,3,0,1, 3,2,1,1, 2,1,4,3, 2,2,1,0, 1,0,5,1, 2,1,3,2, 2,2,2,2, 4,1,1,2, 3,2,0,1, 1,2,1,0) もし…

Rで逆推定

R

JMPには逆推定というコマンドがある 次のようなデータを用意して、線形回帰を行うとする。 temp <- rep(c(15, 20, 25), c(5, 5, 5)) day <- c(34, 33, 36, 37, 35, 30, 28, 29, 25, 28, 23, 25, 22, 24, 26) こんなかんじで適当に。 > result <- lm(day~temp…

不審なファイルの存在確認等

事の発端 【注意喚起】潜伏しているかもしれないウイルスの感染検査を今すぐ!:IPA 独立行政法人 情報処理推進機構www.ipa.go.jp 「バッチでの配布でも可能です」とか書くならバッチファイル用意しといてほしい…。 バッチファイル 上からやれといわれたけど…

ggplot2で縦に並べたグラフの横幅を揃える

R

(※2015/06/03 20:27追記あり) 下記の記事に基づく発言に関連してどうも某所で勘違いが発生しているようなので。R Recipe: Aligning Axes in ggplot2 | Exegetic Analyticswww.exegetic.biz 使用データ 下記のように生成したものを用いる。 x <- 0:100 x <-…

RStudioの折りたたみ機能

R

今更ながらコードの折りたたみ機能が付いている事に気付いた。 詳細はCode Folding and Sections – RStudio Supportを参照。 折りたたみ設定されるもの ブレース(波括弧{})で括られた領域。 R SweaveやR Markdownドキュメントにおけるコードチャンク。 R M…

S3

R

RStudioではじめるRプログラミング入門を読んだ。RStudioについて詳しい情報があるかと思っていたのだが、これについてはそこそこの記述だった。 ただ、全体的に見て、Rとプログラミングをゼロから学ぼうと思うのであれば良い情報量だと感じた。リゲス本なん…