r-statistics-fanの日記

統計好きの現場の臨床医の覚書のようなもの

統計推理。どう計算したのだろう。比率の信頼区間

どう計算したのだろう。比率の信頼区間日経メディカルで、コデインの小児への処方が制限されたと言う記事を見た。PDMA(医薬品医療機器総合機構)が運営するMID-NETというデータベースを利用した解析結果が参考にされたらしい。http://www.jpa-web.org/dcms_…

100人を部屋に集めてお金をランダムな相手に渡し続けるとだんだんと貧富の差が生まれる~条件により印象通りの結果に:第2弾

GIGAZINEでの興味深い記事の検証第2弾。 http://gigazine.net/news/20170711-random-people-give-money-to-random-other-people/ 100人を部屋に集めてお金をランダムな相手に渡し続ける」とだんだんと貧富の差が生まれる というもの。記事の中では「45ドルを…

GIZAZINEの記事をなぞってみる:「100人を部屋に集めてお金をランダムな相手に渡し続ける」とだんだんと貧富の差が生まれる

gigazine.net#7/12追記 普通にソートした順でboxplotをすれば分かりやすかったので最下部に図を追加GIGAZINEで興味深い記事が。 「100人を部屋に集めてお金をランダムな相手に渡し続ける」とだんだんと貧富の差が生まれる というもの。記事の中では「45ドル…

#OsakaStan5参加 献本争奪じゃんけんで2回で決着した事象は珍しいものだったのか?

はじめて統計の勉強会に参加した。#7/2一部バージョン違いのコードだったので修正いわゆるアヒル本の勉強会である。 www.amazon.co.jpなんと著者の松浦先生もいらっしゃると聞いて駆けつけた次第である。何故アヒル本と呼ばれるのか? それは、この表紙のブ…

レイド戦後捕獲時のバンギラスCP/固体値と最大強化時CP/レア度順位の表

統計ブログ()だが、わりとポケモン記事へのアクセスが多い。 というわけで追加記事。レイド戦が始まり、バンギラスが身近になった。 前回の記事で http://r-statistics-fan.hatenablog.com/entry/2017/05/22/174305 卵を49個孵化させないと、90%以上の確率…

データの出力時にrownamesを出したくない/ポケモンGO,ラッキーCP777になる固体値全部

ふつうに行列をコンソールに出力すると> lucky777 AT DF HP Lv [1,] 0 2 14 24 [2,] 0 3 11 24 [3,] 0 3 12 24 [4,] 0 4 8 24 [5,] 0 4 9 24 のようにrownamesがNULLであっても、[1,] のように何行目か 分かるようにコンソールに出力される。 これを消したい…

超ムカつくtibbleクラス:バージョンアップでエラー頻発

Rのバージョンアップで、エラーが出ることが多くなった。その多くはどうもtibbleクラスが関わっている! 色んなパッケージがtibbleではエラーを吐きまくる。読み込んだ段階でdeta.frameがtibbleになっていることがしばしばである。たとえばこんな感じ > str(…

ポケモンGO バンギラスの厳選はどこで妥協するか。Rで計算してみる

ポケモンgoで岩祭り開催中。 ヨーギラスがよく出るので、厳選したいが、どの程度で妥協すればよいだろうか。 バンギラスの各攻撃、防御、HP(スタミナ)毎での、最大CPやレア度を知りたい。 Rで計算してみる。#順番をCP強い方が先頭に来るように改定。CP順位…

R3.4.0にアップデート。Rが高速に。

R3.4.0が出たと聞いた。 JITバイトコンパイラになり、loopが高速になるらしい。 適当にフィボナッチ数を計算する関数を作って試してみる。 fibonacci <- function(n){ if(n == 1) 1 else if(n == 2) 2 else (fibonacci(n-2) + fibonacci(n-1)) } fibonacci2 …

探偵ナイトスクープ 得点した選手の背番号の合計は偶数が多い?

ナイトスクープでサッカーの試合で得点した人の背番号の合計が 偶数有利か奇数有利かということを検証していた。偶数背番号の選手も奇数背番号の選手もおなじ得点率なら同じになるだろう。 しかし、実際は異なる。 偶数の選手が得点することがたまたま多かっ…

受験番号と合格率の関係

受験番号と合格率の関係 rs_fan_jp 2016年4月3日 2016年4月3日20:20念のため追記。言うまでもなく受験番号から合格への直接的因果関係はないからね。統計上有意ですが、単なる見かけ上の関係です。間違っても受験番号1ゲットなどという不毛なことを目指さな…

レアなタクシーに遭遇する確率

レアなタクシーに遭遇する確率京都で学会があり、レアなタクシーに遭遇する機会があった。ヤサカタクシーは通常三つ葉なのだが、四つ葉と双葉 が存在する。それぞれ1200台中4台と2台しか無いらしい。今回同じ日に四つ葉と双葉を見た。なんだか良い気分。多分…

ババ抜きで13枚以上のスタートになる確率

ババ抜きで13枚スタートになる確率異動が決まり、週末は長距離移動かあるいはフルに仕事。 サンデープログラマーなので最近記事がかけていない。 秋からは精力的に更新したいものだ。某アニメでババ抜きのシーンが有った4人でゲームを行い、なんと13枚スタ…

RMST生存曲線下面積2

#RMST生存曲線下面積2前回の記事で(surv2sampleComp)パッケージを使ってみた http://r-statistics-fan.hatenablog.com/entry/2014/08/04/225135今回survRM2パッケージが出たことを知った。使ってみる 前回と同じく、イレッサのシミュレーションデータを使用…

RCTで多変量解析#2

http://r-statistics-fan.hatenablog.com/entry/2015/02/19/000248 前回の続きRCTで多変量解析をするとPower検出力が上がる。 しかし、弱い予測因子を投入するとPowerは下がるはずだというのを 前回実証しようとして失敗した。今回は、投入する因子を増やし…

RCTで多変量解析

某社のwebカンファを見てきた。RCTで多変量解析するのは、by-chanceで出た差を調整して 確認するためなどという感じのことを言っていた。しかし、自分の認識は異なる。 RCTで多変量解析するのは、原則的にあらかじめpowerの増加を 期待して計画された時だけ…

頭脳王:三山崩しの「逆型」をRで解く。更に高速化一般化

頭脳王:三山崩しの「逆型」学生時代の旧友よりメールが来た。 数学に関しては私にとっては神様のような方で、学生時代には よく自作の数学問題を披露してくれたりして、楽しんだものだ。懐かしい。旧友によると、三山崩しに関しては非「逆型」の必勝法は数…

C++導入

visual studioの無料版を導入した。10GBとかでて、気が遠くなった。とりあえず、ちょっと触ってみた。 #include <stdio.h> #include <limits.h> #include <iostream> int main(){ std::cout << INT_MAX << "INT_MAX\n"; std::cout << LONG_MAX << "LONG_MAX\n"; std::cout << CHAR_BIT << </iostream></limits.h></stdio.h>…

Rcppを使う:その2

C言語による乱数生成 ここにメルセンヌツイスターのソースコードがあるので利用してみる。あらかじめMT.hをincludeフォルダにおいておく必要がある。以下のファイルをpi4.cppとして保存 #include <Rcpp.h> #include <stdio.h> #include <MT.h> using namespace Rcpp; // [[Rcpp::exp</mt.h></stdio.h></rcpp.h>…

Rcppを試す

Rcppを試すhttps://rs-fan-jp.shinyapps.io/surv_jp/ H22 survival data in Japan前に自作したWebアプリで、すでに人生の折り 返し時点を過ぎてしまっていることに気がついた。 親の平均余命もなんか具体的すぎてグッと来る。時間は無限ではないので、ここら…

暗殺教室の数学問題をRで描画する

暗殺教室の問題をRで描画する追記 hoxo_m さんが、ブラウザでグリグリうごかせるファイルをアップしてくれました。 WebGLファイルの作成方法は知っていましたが、なるほどドロップボックスに 出来たファイルを置けばみんなハッピーにうごかせるのか~。RGL m…

SAS proc meansの重み付けの挙動

SASの%stddiffを使って、解析をしたい。 しかし、Rで完結させたい。そのため、前回翻訳をした。今回、年末年始を利用して最終解析をしようとしたが、その前に SASで検証してみた。げげー。なんだか生データから検証すると一部一致しない。ところが生データか…

日本人の生命曲線をshinyでインタラクティブに表示する

#追記;日本語化してみた。文字によって相性があるようだ。 最終的に何とか日本語でも通った。なお、タバコ吸っている人は実年齢+10才すると大体正しい生存曲線になりますとりあえず、練習でshinyを使ってみた超単純な内容なのに、何故か単独の関数では…

頭脳王2014 三山崩しをRで解く(必敗型スタートの人と必勝型スタートの人がいた件)

#追記 #録画を見なおしたら衝撃的事実が発覚 #なんと敗者の初期値は必敗の組合せだった! #ひでえ。放送事故やん追記2 Nim - sigma425のブログ すでに言及されていた。しかも、コメント欄が熱い!!!追記3 本人のコメント? https://twitter.com/sou_mizu…

ロジスティック回帰から計算される確率の信頼区間

最近忙しくて更新ができないー。でも、ついさっきデスマコロシアムには挑戦した。 正解なら今回もRで最短のはず。ワクワク →12/13 0:30現在最短だ~ 12/15盛大に抜かれた。どうやっても勝てる気がしない。CodeIQ | 『第8回デスマコロシアム』問題 集計報告 …

グローバル環境のデータたちをまとめてコピペできる形にtext化

f1 <- function(){ var.list <- ls(.GlobalEnv) for(i in var.list){ cat(i, "=") eval(parse(text = paste0("dput(", i, ")"))) } } f1() まとめて変数も関数もデータフレームもコピペできる形に出力

LED電球逝く~ワイブル分布で解析する

LEDライトが1年で逝った。 超損した気分。悔しいので、ワイブル分布で平均寿命10年の際に1年で 逝く可能性がどんなもんか計算してみる。 ###LEDライト sh <- seq(1, 10, by = 1) #shape sc <- 10 / gamma(1 + 1 / sh) #scale 平均生存期間 10 year pv <-…

データフレームの内容をテキスト化

自分のRのデータフレームの内容を、メールやblogなどで伝え、他人のRで再現させる関数 - r-statistics-fanの日記 自分のRのデータフレームの内容を、メールやblogなどで伝え、他人のRで再現させる関数 - r-statistics-fanの日記こちらで、書いた関数だが、吉…

シミュレーションで使ったワイブル分布

生存曲線の95%信頼区間のカバー率の記事シングルアームの生存率解析デザイン(2)信頼区間の真値カバー率 - r-statistics-fanの日記 シングルアームの生存率解析デザイン(2)信頼区間の真値カバー率 - r-statistics-fanの日記で使った分布における、理…

シングルアームの生存率解析デザイン(2)信頼区間の真値カバー率

信頼区間の真値カバー率をシミュレーションする。ここ http://www.sascom.jp/download/pdf/usergroups11_A-02.pdf で、シミュレーションされているが、前提によってデコボコとカバー率が変わる。 平均的に、loglogが良さそうという話であった。今回は、指数…

シングルアームの生存率解析デザイン

生存曲線のシングルアーム試験知り合いが試験を計画しているようだ。 なお、結局関係各所を通すには大学の統計家のcheckが 必要とのことで、素人の私など出る幕はない。しかし、興味があるのでやってみる。(もともと興味だけで動いている)2アームならば、…

自分のRのデータフレームの内容を、メールやblogなどで伝え、他人のRで再現させる関数

#10/31追記 データフレームの内容をテキスト化 - r-statistics-fanの日記なんとdput()関数で目的のことは出来るようです。 苦労する必要はなかったんや~ ##追記終了メールやblogなどで、自分の環境上のデータフレームを記載したいことがある。 しかし、a…

慶応がベスト1()~何の雑誌だコレ

https://twitter.com/dankogai/status/522047534022156289算数をなめてる> @kirik: うわ、なんだこれw 算数やんけ RT @fukusanity: 作った人間が統計をなめてることはわかった pic.twitter.com/mcDF5k3jzA— Dan Kogai (@dankogai) October 14, 2014ツイッタ…

デング熱グラフ10/15

デング熱グラフを更新した。10/15またひとり増えた模様 最近新規報告がない模様。そろそろ収束か?

Rで標準化効果量を計算する(Rewriting the code of SAS proc %stddiff to R)

どうしてもstandardized difference(以下stdiff)を計算 する必要が出た。いろいろ調べたが、どうもSASがスタンダードになっている みたい。しかし、どうしても順序変数のstdiffの計算方法が分からない。 いきなりメールを送って良いものか迷いながら、知り合…

9/26現在デング熱グラフ

グラフ更新。代々木公園は感染がおさまってきたように見える。 寒くなって来たしね。一方代々木公園以外は数は少ないものの、くすぶっている。 代々木公園以外は推定日がわからない人が多いので、個々の 潜伏期の確率に従って、微妙な期待値の人数のグラフに…

gmpのisprimeはどこまで正確か~pseudoprime

前回の記事でgmpのisprimeがあまりにも正確というのを書いた。デフォルトのrep=40だと、何年回してもprobably primeなのに実際は そうでない例(pseudoprime)を見つけられないと思う。 しかしreps = を少ない数にすればpseudoprimeを見つけられそう。 rでや…

9/24現在のデング熱のグラフ~厚労省もグラフを作ってきたぞ

http://www.mhlw.go.jp/bunya/kenkou/kekkaku-kansenshou19/dl/20140924-01.pdfhttp://www.mhlw.go.jp/bunya/kenkou/kekkaku-kansenshou19/dl/20140924-01.pdf ついに厚労省でもグラフを表示するようになったようだ。良い傾向である。自分のはこれに、感染日…

スコットランド独立賛成票が素数~Rで素数をチェックする関数色々

某つぶやきで、スコットランド独立賛成票が素数だと知った#2014.09.23中澤さんの別解など追加しましたhttp://business.nikkeibp.co.jp/article/topics/20140919/271500/目的 #本当に素数か確かめる。 #その前後にどのくらい素数があるか。 #ついでに、素…

多倍長出力をキレイに桁を揃えて出力する(改良版

R

多倍長精度計算結果を見やすく出力 - 裏 RjpWiki さんの所で、多倍長出力をキレイに出力する 記事があった。自分の前の記事だと 小数点以下が出力されて、無駄に長くなってダメダメですね。 反省です。今回もコードを勉強させていただいた。 ありがとうござ…

デング熱のグラフ9/19

今回はグラフだけ。代々木公園以外が目立ってきていますね。#一部入力データミス発見し修正。他にもあるかもしれんので参考程度に。

デスマコロシアム:素数を出力

R

library(conf.design) #この行は抜け道でカウントされないのだ。 cat(primes(997),sep=":") 今回もRで最短だったが、都合の良いideone依存の パッケージを発掘しただけのこと。 自慢にはならないわなー。でも2連続R最短で嬉しい。 mamekinさんの正攻法最短…

デング熱グラフ更新2014/9/18分

今回は感染推定日や発症日から発表までにかかる日数をグラフ化してみた。感染から発表までの日数は17日でほぼ安定。 発症から発表までの日数は12日でほぼ安定。ということで、感染日は9/1まで、発症日は9/6までは 大体固定しており、今後大きくは変わらない…

メタアナリシス reiv()

R

メタアナリシス某つぶやきで、reiv()の存在を知った。 http://www.youtube.com/watch?v=0XFlPB0mEzA http://zanet.biz/R/fct/ma_metafor_unibb.zipよりダウンロードできるらしい。そして、自分自身は使うことはないものの、他人に教えることが多いEZR と結果…

9/12現在のデング熱のグラフ。蚊に刺された人の割合の信頼区間も追加

2014/9/12の18時現在のデータに更新 #18:56年齢階級別の症例数も追加今回今までになかったパターンの出現があり。 代々木公園周辺+新宿中央公園どちらも行っている例が出た。 ただ圧倒的に代々木公園が多いこともあり、今回は代々木公園の 例として扱うこ…

デング熱のグラフ:感染日と発症日9/11

#22:20発症日がずれていたのを手直し。本日9/11の21時現在の最新データに更新 http://www.mhlw.go.jp/bunya/kenkou/kekkaku-kansenshou19/dl/20140911-01.pdf代々木公園以外が増えすぎたので、代々木公園以外にまとめちゃいました。 代々木公園かそれ以外の…

2014/09/10 updte デング熱感染日と発症日:Rで同じグラフ内に表示する

本日19時現在のの最新のデータに更新 http://www.mhlw.go.jp/bunya/kenkou/kekkaku-kansenshou19/dl/20140910-01.pdf代々木公園近辺以外にも広がっている模様。http://www.who.int/mediacentre/factsheets/fs117/en/ http://www.mhlw.go.jp/bunya/kenkou/kek…

デング熱の推定感染日をbarplotしてみた。Rで欠損値のあるbarplotをうまく表示する

#追記 ぐぬぬ、夕方は無かったのに、今見たら新データが出てる。 新データに更新した。 発症日も訪問日も未定なヒトがいて、コードを書き換える必要が。 一つのデータのせいで、修正に20分もかかったよ。本日夕方以降に発表のデータにて、8/23-26あたりが…

某疾患患者の某公園や周辺訪問日をカレンダープロット:今回は潜伏期の確率密度で調整

#追記 昨日と同じカラーのも追加。一部バグ修正。 #感染推定日barplotも追加 #発症日のbarplotも追加参考にした潜伏期 http://www.plosone.org/article/info%3Adoi%2F10.1371%2Fjournal.pone.0050972PLOS ONE: The Incubation Periods of Dengue Viruses …

カレンダープロットで某疾患の某公園の訪問をプロット

カレンダープロット - 東京で尻を洗う声優の誕生日のカレンダープロット - 驚異のアニヲタ社会復帰への道カレンダープロッットが流行しているので、やってみた。某疾患の某公園の訪問を表示する。複数候補日がある人は、候補日数で割って 重みとする。つまり…