axjack's blog

axjack is said to be an abbreviation for An eXistent JApanese Cool Klutz.

R

sweep関数の使い方を忘れがちなので自分へのメモ

まとめ sweep関数の MARGIN = 1の時は i 行目のそれぞれの値に対してSTATSの i番目の値をFUNと二項演算する MARGIN = 2の時は j 列目のそれぞれの値に対してSTATSの j番目の値をFUNと二項演算する お気持ち 行平均の引き算のお気持ち sweepで1行目の行平均を…

雑だけどirisを機械学習する

まえがき Data Science: Machine Learningを受講しています。今までのR BasicsやVisualizationに比べて課題が多くてなかなか進まないです。8月中には完了させたいですね。 本題ですが、今回は「雑だけどirisを機械学習する」ということで、とりあえずcaret使…

2次元正規分布のデータから2通りの方法で回帰直線を引く

タイトルの通りです。 lmで単回帰直線 平均と分散共分散からE[Y|X=x] *1を計算 の2通りで回帰直線を引きます。 コード library(MASS) library(scales) options(digits = 3) # 平均 mx <- 10 my <- 22 # 分散 Sx <- sqrt(9) Sy <- sqrt(16) Sxy <- sqrt(9.6)…

標準誤差は統計量の標準偏差であることを確かめる

やること 離散型確率分布 X p 17 0.3 -1 0.7 に於いて、大きさn=100の標本を抽出する時、 和S の標本分布について、Sの平均、標準誤差 平均M の標本分布について、Mの平均、標準誤差 をそれぞれ 計算 モンテカルロ シミュレーション によりそれぞれ求める。 …

Rで文字列→数値に変換する際、NAs introduced by coercionが出て困った時のtips

次のような変換を考えます。 "1名" → 1 "2名" → 2 "なし" → 0 "調査中" → NA これをdplyrのパイプの中でmutate( case_when(...) )を駆使して実行していたのですが、エラーとなってしまいました。 データフレーム(見栄えのためtibble) mydf <- tibble( 同居…

edXのHarvardX's Data Scienceを受講しています。

R

edXのHarvardX's Data Scieceとは? → HarvardX Data Science Professional Certificate | edX とりあえずR Basicsが終わってVisualizationのIntroの途中まで来ました。2年前も同じコースを受講したのですが、その時は途中でドロップアウトしてました・・・…

一元配置分散分析をRで実装する

手計算でもaovでも一元配置分散分析は出来るので、Rで実装してみようと思った次第です。F分布の累積分布を除いてほぼほぼベクトル演算?を使っています。データは水準の繰り返し数のトータルから適当に生成しています。統計学的観点はほぼ0です。 実装 # デ…

Rで線形計画法

R

lpSolve*1を使って線形計画法を解いてみた。理論は難しいから敬遠していたがソルバー使うと一瞬で解けて気持ち良い。 関数lpに渡すパラメータは以下の通り。第一パラメータは"max"か"min"。詳しいことはrdocumentationなどを参照。 f_obj : ベクトル。目的関…

母比率の信頼区間に含まれる2次不等式を解く

母比率の信頼区間 2次不等式を解く 式変形 具体例で検証 公式を用いる 2次不等式を解いた結果を用いる 参考 母比率の信頼区間 母比率 の母集団からサイズの標本を抽出する。このとき標本割合 について、 は近似的に平均、分散 の正規分布 に従う。ただし、…

行簡約行列をRで

pracmaのrrefを使って行簡約行列を出してみます。これで適当な行列を手計算で簡約化して答え合わせができますね。 ソース # install.packages('pracma') # library('pracma') ## 行ベクトルを4本 a1 <- c(1,2,0) a2 <- c(2,4,0) a3 <- c(0,1,3) a4 <- c(1,3,…

『確率4万分の1、県の入札くじに6回連続当選「奇跡的」』なのかどうかを調べよう

www.asahi.com 動機 面白いニュースを見かけたので計算してみましたなポストです。記事のこの部分に注目しました。 6回のくじにはそれぞれ制限価格で入札した3~8の鑑定業者が参加した。6回連続で当たる確率を計算すると約4万分の1となる。 奇跡ってど…

超幾何分布

何度やっても忘れるのでブログに書いて覚えよう。 超幾何分布の確率質量関数 計算例 その1 男50人女50人から10人を選ぶ。10人のうち男3人女7人となる確率pは、超幾何分布を用いて その2 統計検定準一級2017年問10より引用。表の80人から30人を無作為抽出する…

連の検定(Runs Test for Detecting Non-randomness)

はじめに 連の検定についてよく分からなかったのでRで実装して確かめてみました。 連とは? 2値{A,B}を取る系列ABAABBAAABBBABがあった時、これをA | B | AA | BB | AAA | BBB | A | Bと連続する同じ値ごとに分割できるように見える。この時同じ文字または…

重回帰分析をRで

やること ソースコード 結果の確認 (1)回帰式 (2)自由度調整済み寄与率 (3)同じ地区で広さ=70平米, 築年数=10年, 価格=5.8千万円の提示は妥当か やること 永田『多変量解析法入門』(以下参考書)よりp.2のデータをもとに重回帰分析を行う。 ソースコード # デ…

初期値1から始めたニュートン法による平方根の近似値計算は、7回ほど反復計算しておけば大丈夫そうです。 

はじめに ニュートン法による平方根の近似 例 結論 コード 結果確認 要約統計量 箱ひげ図 おわりに はじめに 今日で令和元年ゴールデンウィークも無事終了です。 前回に引き続きニュートン法で平方根を求めてみました。今回は、N = [1,1000] の範囲でニュー…

ニュートン法で平方根の値を求める

R

きっかけ 電卓のルートってどうやって求めるんだっけ? → たぶんニュートン法だろうなぁ。ということで、Rで簡単に実装して計算してみました。 実装 ニュートン法 fnewton <- function(xn=1,p=2,a=2){ return ( xn - ( xn^p - a )/( p * ((xn)^(p-1)) ) ) } …

ggplot, geom_point, facet_gridの練習

R

はじめに ggplotやqplotを日々練習しています。ggplot, geom_point, facet_gridを組み合わせると、 グラフ内側のx軸(量的) グラフ内側のy軸(量的) グラフ外側の横側(質的) グラフ外側の上側(質的) 点の色(質的・量的) のように5変数くらい同時に、割と分か…

所沢市の平成29年版統計書-2 人口(その1)-1.人口の推移から、人口の総数=男+女となっていることを確認する

はじめに 動機 やること データソース コード sublimeでの加工 1. 総数-男-女の列を取得する 2. 数値表現をカンマ有りからカンマ無しへと置換する 3. ダブルクオーテーションを削除する Rで読み込み 結果確認 おわりに はじめに 動機 昨今不正統計が話題であ…

次の元号をRで予想してみる〜前処理編〜

はじめに コード rvestの準備 rvestでwebスクレイピング→テキスト処理(1文字ずつに分解) 結果確認 終わりに 参考にしたサイト はじめに 平成もあと3ヶ月足らずで終わってしまいますね。色んな方が次の元号を予想している流れがあるので、自分でも予想して…

RStudioにあるTerminalで文字が打てない時の対処(Mac)

掲題の通りです。 対処 キーボード環境設定を開き、「英語」の入力ソースを追加する。(※[A] ABC と表示される) これでRStudioにあるTerminalで文字が打てるようになるはずです。

顎がガクッと

あけましてゴールデンウィーク2018ですね。私は初日のお昼にベーコンエピをもぐもぐしていたら、右顎の奥の方から「ガクッ」というか「パキッ」というかそんな音が聞こえて、すぐさま痛くなりました……。妻には顎関節症だよ、と言われ、ついになってしまった…

axjack is said to be an abbreviation for An eXistent JApanese Cool Klutz.