19/01/21 23:54:21.71 zCz5KSqt.net
タバコって燃えるゴミかな？
有害物質を発生するからプラゴミと同じかな？
灰皿は不燃ゴミ
電子タバコは家電ゴミでいいのかな？
776:卵の名無しさん
19/01/22 00:16:37.39 7DkisZz5.net
>>734
吸ってる椰子がゴミ
777:卵の名無しさん
19/01/22 17:20:22.09 XyJ4p7WH.net
f <-function(a){
qa=qnorm(a,50,10,lower=F)
ma=integrate(function(x) x*dnorm(x,50,10)/a,qa,Inf)$value
va=integrate(function(x) (x-ma)^2*dnorm(x,50,10)/a,qa,Inf)$value
c(mean=ma,sd=sqrt(va))
}
m1=f(0.27)[1]
s1=f(0.27)[2]
m2=f(0.50)[1]
s2=f(0.50)[2]
g <- function(x){
y= (x-50)/s1*10+m1
50+(y-m2)/s2*10
}
x=40:70
g(x)
778:卵の名無しさん
19/01/22 17:40:34.59 XyJ4p7WH.net
f <-function(a){
qa=qnorm(a,50,10,lower=F)
ma=integrate(function(x) x*dnorm(x,50,10)/a,qa,Inf)$value
va=integrate(function(x) (x-ma)^2*dnorm(x,50,10)/a,qa,Inf)$value
c(mean=ma,sd=sqrt(va))
}
m1=f(0.27)[1]
s1=f(0.27)[2]
m2=f(0.50)[1]
s2=f(0.50)[2]
g <- function(x){
y= (x-50)/s1*10+m1
50+(y-m2)/s2*10
}
g(c(50,55,57.5))
母集団（例えば18歳人口)の成績が平均値50点標準偏差10点の正規分布に従うとする。
上位27%が進学すると
平均値62.2点標準偏差5.01点の集団Aとなる。
(成績が62.2点なら偏差値50に67.21点なら偏差値60になる)
上位50 %が進学すると
平均値58.0標準偏差6.03の集団Bとなる。
集団Aでの偏差値50, 55, 57.5は各々
集団Bでは偏差値57.0, 73.0, 81.9になる。
779:卵の名無しさん
19/01/22 18:22:22.34 XyJ4p7WH.net
>>737
バグ修正
f <-function(a){
qa=qnorm(a,50,10,lower=F)
ma=integrate(function(x) x*dnorm(x,50,10)/a,qa,Inf)$value
va=integrate(function(x) (x-ma)^2*dnorm(x,50,10)/a,qa,Inf)$value
c(mean=ma,sd=sqrt(va))
}
m1=f(0.27)[1]
s1=f(0.27)[2]
m2=f(0.50)[1]
s2=f(0.50)[2]
g <- function(x){
y= (x-50)/10 *s1+m1
50+(y-m2)/s2*10
}
g(c(50,55,57.5))
母集団（例えば18歳人口)の成績が平均値50点標準偏差10点の正規分布に従うとする。
上位27%が進学すると
平均値62.2点標準偏差5.01点の集団Aとなる。
(成績が62.2点なら偏差値50に67.21点なら偏差値60になる)
上位50 %が進学すると
平均値58.0標準偏差6.03の集団Bとなる。
集団Aでの偏差値50, 55, 57.5は各々
集団Bでは偏差値57.0, 61.2 63.3になる。
780:卵の名無しさん
19/01/22 18:32:41.55 XyJ4p7WH.net
f <-function(a){
qa=qnorm(a,50,10,lower=F)
ma=integrate(function(x) x*dnorm(x,50,10)/a,qa,Inf)$value
va=integrate(function(x) (x-ma)^2*dnorm(x,50,10)/a,qa,Inf)$value
c(mean=ma,sd=sqrt(va))
}
p1=0.27
p2=0.50
m1=f(p1)[1]
s1=f(p1)[2]
m2=f(p2)[1]
s2=f(p2)[2]
g <- function(x){
y= (x-50)/10 *s1+m1
50+(y-m2)/s2*10
}
g(c(50,55,57.5))
xx=seq(40,80,0.1)
plot(xx,g(xx))
781:卵の名無しさん
19/01/22 18:47:08.16 XyJ4p7WH.net
f <-function(a){
qa=qnorm(a,50,10,lower=F)
ma=integrate(function(x) x*dnorm(x,50,10)/a,qa,Inf)$value
va=integrate(function(x) (x-ma)^2*dnorm(x,50,10)/a,qa,Inf)$value
c(mean=ma,sd=sqrt(va))
}
sim <- function(x=50,p1=0.3,p2=0.5){
m1=f(p1)[1]
s1=f(p1)[2]
m2=f(p2)[1]
s2=f(p2)[2]
y= (x-50)/10 *s1+m1
50+(y-m2)/s2*10
}
sim()
782:卵の名無しさん
19/01/22 23:57:18.92 poZn4cZ2.net
>>740
仕事は何をしてるのかね？
いまどき正規職員でないものは完全に負け組だぞ
783:卵の名無しさん
19/01/23 00:04:12.53 p2qfG+zh.net
>>741
何のプログラムかわかる？
784:卵の名無しさん
19/01/24 01:16:27.12 sxN7Of5+.net
テロ工作機関
福山友愛病院
なんと朝木明代市議のような他殺だと日本人に騒がれるから
日本国警察が、
朝鮮殺戮殺人教団に乗っとられたことにより警察から与えられていたのは、
犯罪ライセンス！！！！！！
病院で
785:医師から処方された薬は、何ら疑いを持つことなく飲む人も少なくないだろう。 医学や薬学の知識がある人や、持病で何年も投薬治療を行っている人などは、 また違った見方をしているかもしれないが、一般的に「医師からの処方薬」に対する信頼は篤い。 しかし、そんな「信頼」を揺るがすとんでもない事件が話題となっている。 ■「在庫処理」で不要な治療薬を投与広島県の福山友愛病院は、昨年11～12月の間に、 統合失調症などの患者6名に、本来必要のないパーキンソン病の治療薬を投与していたことが判明。 これは、病院を運営する医療法人会長の指示によるもので、病院側は「使用期限が迫った薬の在庫処理がきっかけのひとつ。」と説明しているとのこと。 会長は、同病院で精神科としても勤務しており、パーキンソン病の治療薬である 「レキップ」の錠剤を統合失調症患者など6人に投与するよう看護師らに指示。 投薬は複数回行われ、末丸会長は通常の 8倍の量を 指示していたという。
786:卵の名無しさん
19/01/24 07:39:17.09 YHL/Shf0.net
ところで、タバコって燃えるゴミですか？
有害物質を発生するのでプラゴミと同じ扱いでしょうか？
電子タバコは家電ゴミですよね
今度、小学校に教育実習に行くので
お家のタバコ用具を全てゴミに出すように指導してきます
787:卵の名無しさん
19/01/24 11:08:18.87 EfYrEwvB.net
URLﾘﾝｸ(www.niid.go.jp)
H3N2で 2/21vs0/21で有意差あるか？
n<5だとΧ二乗検定はあてにならない。
788:卵の名無しさん
19/01/24 15:16:53.44 7fT+1bS6.net
It is common knowledge among doctors and patients that Do-Teihen(exclusively bottom-leveled medical school) graduates mean morons who bought their way to Gachi'Ura(currently called by themselves)
According to the experience of entrance exam to medical school in the era of Showa, when the sense of discrimination against
privately-founded medical schools were more intense than it is now,
all such schools but for Keio had been so compared to some specialized institution for educable mentally retarded kids that nobody but imbecile successors of physicians in private practice had applied for admission.
There had been NOT a single classmate who chose willingly against his/her common sense to go to the Do-Teihen(exclusively bottom-leveled medical school, currently also known as Gachi'Ura),
which would have cost outrageous money and its graduates are destined to be called Uraguchi morons who bought thier way into the Do-Teihen, by thier colleagues and even by thier own clients.
Although people won't call them names to their face,
certain 80-90% people of about my age have been yet scorning and sneering at Uraguchi graduates, speaking in the back of our mind,
" Uraguchi morons shall not behave like somebody."
We never speak out face to face in real life.
789:卵の名無しさん
19/01/24 18:59:58.86 oaTB7QH8.net
テロ工作機関が日本にあることが判明
福山友愛病院
日本人への薬物大量投与テロも発覚
また他は、
単純に、朝鮮殺戮殺人革命の為に、
日本人被害者側を、
カルト医師が医師免許を犯罪に使用し、
偽造カルテをでっち上げているにすぎない
でっち上げる理由は、カルト会社の犯罪にごねている、
だから、
警察から犯罪ライセンス与えられているので、
虚偽説明して騙して閉鎖病棟に閉じ込め、
携帯とりあげ、
監禁罪やり精神保健福祉法違反をやり偽造カルテ書きますね
それだけw
だって、警察がグルで共同組織犯罪をしているから、
犯罪を犯罪でないことにしてくれるんだもの
日本人被害者側を拉致していることがバレたしな
それをやらかしていたので、
福山友愛病院には電凸し、
インターネット上に書いていることを通告もしている
ケースワーカーのイソムラという女性曰わく、
カミサカ先生インターネットに書かれて辛いみたいで～・・・
日本国内にいるテロリストが辛くなることは、
日本人の喜びだっつーのww
790:卵の名無しさん
19/01/24 20:31:49.63 H5DL72sh.net
Rnorm<-function(n,m,s){
x=rnorm(n,m,s)
m + s*(x-mean(x))/sd(x)
}
N=1e3
x=replicate(5,Rnorm(N,50,10))
f <- function(){ # simulated test results of 3 or 5 subjects
y=numeric()
for(i in 1:5) y[i]=sample(x[,i],1)
y3=sum(y[1:3])
y5=sum(y[1:5])
c(y3,y5)
}
p2h <- function(x){ # points to hensa-chi
m=mean(x)
s=sd(x)
x=50+10*(x-m)/s
}
re=replicate(1e5,f())
y3=re[1,] # points of 3 subject test
y5=re[2,]
z3=p2h(y3) # hensa-ti of 3 subject test
z5=p2h(y5)
791:卵の名無しさん
19/01/24 20:32:02.08 H5DL72sh.net
hc <- function(h,d=0.5){ # hensa-chi conversion
idx=which(h-d<z5 & z5<h+d) # hensa-chi range in 5 subject test
mean(z3[idx]) # mean of its hensa-chi in 3 subject
}
hc=Vectorize(hc)
hh=40:80
plot(hh,hc(hh),bty='l',asp=1,xlab='5 subject',ylab='3 subject',
pch=19)
abline(a=0,b=1,col='gray')
792:卵の名無しさん
19/01/24 20:51:53.08 H5DL72sh.net
実行結果
URLﾘﾝｸ(i.imgur.com)
793:卵の名無しさん
19/01/27 05:17:22.65 Ba4uprm1.net
札幌ひばりが丘病院
麻薬取締法違反で書類送検
URLﾘﾝｸ(video.fc2.com)
794:卵の名無しさん
19/01/28 11:16:13.00 C5thEPa8.net
ゆっくり霊夢はFランク大学の就職課に就職したようです
URLﾘﾝｸ(www.youtube.com)
795:卵の名無しさん
19/01/28 17:29:00.31 AVoD+Cyq.net
働いたら負け
796:卵の名無しさん
19/01/28 21:53:43.21 +PCpJyrP.net
It is common knowledge among doctors and patients that Do-Teihen(exclusively bottom-leveled medical school) graduates mean morons who bought their way to Gachi'Ura(currently called by themselves)
According to the experience of entrance exam to medical school in the era of Showa, when the sense of discrimination against
privately-founded medical schools were more intense than it is now,
all such schools but for Keio had been so compared to some specialized institution for educable mentally retarded kids that nobody but imbecile successors of physicians in private practice had applied for admission.
There had been NOT a single classmate who chose willingly against his/her common sense to go to the Do-Teihen(exclusively bottom-leveled medical school, currently also known as Gachi'Ura),
which would have cost outrageous money and its graduates are destined to be called Uraguchi morons who bought thier way into the Do-Teihen, by thier colleagues and even by thier own clients.
Although people won't call them names to their face,
certain 80-90% people of about my age have been yet scorning and sneering at Uraguchi graduates, speaking in the back of our mind,
" Uraguchi morons shall not behave like somebody."
We never speak out face to face in real life.
797:卵の名無しさん
19/02/09 22:33:46.30 F1R4fVG9.net
f=function(p,q=0.99){
a=13*q
c=13*(1-q)
b=13000*p
d=13000*(1-p)
(a/(a+b)-c/(c+d)) / (a/(a+b))
}
p=q=seq(0.01,0.99,0.1)
z=outer(p,q,f)
contour(p,q,z)
g=Vectorize(f)
g((1:19)/20)
798:卵の名無しさん
19/02/10 12:37:21.35 9eXx3CuI.net
options(digits=3)
f=function(p,q=0.99){ # p:pyroli in population, q:pyroli in cancer
a=13*q
c=13*(1-q)
b=13000*p
d=13000*(1-p)
x=(a/(a+b)-c/(c+d)) / (a/(a+b)) # attributable to pyroli(+)
y=(c/(c+d)-a/(a+b)) / (c/(c+d)) # qtteibutabel to pyroli(-)
round(ifelse(x>0,x,y),3)
}
g=Vectorize(f)
g((1:19)/20)
p=seq(0.01,0.99,0.01) # pylori in population
q=seq(0.01,0.99,0.01) # pylori in cancer
z=outer(p,q,f)
image(p,q,z,bty='l',nlevels=20,xlab='pyroli in population',
ylab='pyroli in cancer')
contour(p,q,z,bty='l',nlevels=20,xlab='pyroli in population',
ylab='pyroli in cancer',add=T)
contour(p,q,z,bty='l',nlevels=20,xlab='pyroli in population',
ylab='pyroli in cancer')
source('tools.R')
Persp(p,q,z)
Persp3d(p,q,z)
799:卵の名無しさん
19/02/10 12:37:47.54 9eXx3CuI.net
'
がん 非がん
感染 a b
未感染 c d
attributable risk : x
x={a/(a+b)-c/(c+d)} / {a/(a+b)}
'
800:卵の名無しさん
19/02/10 19:40:18.09 9eXx3CuI.net
q1=0.9934
q2=0.9958
f <- function(p,q){
a=13*q
c=13*(1-q)
b=13000*p
d=13000*(1-p)
(a/(a+b)-c/(c+d)) / (a/(a+b))
}
x0=0.985
f1=Vectorize(function(p)f(p,q=q1))
f2=Vectorize(function(p)f(p,q=q2))
curve(f1(x))
abline(h=0.985)
uniroot(function(x)f1(x)-x0,c(0,0.9))$root
curve(f2(x))
abline(h=0.985)
uniroot(function(x)f2(x)-x0,c(0,0.9))$root
801:卵の名無しさん
19/02/11 12:56:16.04 rzg7XY57.net
一般論として医療統計って全然科学じゃないって、研究者の中ではバカにされてますよ。
802:卵の名無しさん
19/02/11 14:44:26.58 Pny0jm5M.net
統計で使われる数学って難し過ぎてついていけないよね。
AICとかわからないままにソフトで出させているだけ。
803:卵の名無しさん
19/02/13 12:53:38.65 hLnbSCPn.net
It is common knowledge among doctors and patients that Do-Teihen(exclusively bottom-leveled medical school) graduates mean morons who bought their way to Gachi'Ura(currently called by themselves)
According to the experience of entrance exam to medical school in the era of Showa, when the sense of discrimination against
privately-founded medical schools were more intense than it is now,
all such schools but for Keio had been so compared to some specialized institution for educable mentally retarded kids that nobody but imbecile successors of physicians in private practice had applied for admission.
There had been NOT a single classmate who chose willingly against his/her common sense to go to the Do-Teihen(exclusively bottom-leveled medical school, currently also known as Gachi'Ura),
which would have cost outrageous money and its graduates are destined to be called Uraguchi morons who bought thier way into the Do-Teihen, by thier colleagues and even by thier own clients.
Although people won't call them names to their face,
certain 80-90% people of about my age have been yet scorning and sneering at Uraguchi graduates, speaking in the back of our mind,
" Uraguchi morons shall not behave like somebody."
We never speak out face to face in real life.
804:卵の名無しさん
19/02/16 13:25:47.75 QaizKjLi.net
>>761
裏口馬鹿にもわかるように原文もつけなきゃ。
原文はこれな。
私は昭和の時代に大学受験したけど、昔は今よりも差別感が凄く、慶応以外の私立医は特殊民のための特殊学校というイメージで開業医のバカ息子以外は誰も受験しようとすらしなかった。
常識的に考えて、数千万という法外な金を払って、しかも同業者からも患者からもバカだの裏口だのと散々罵られるのをわかって好き好んで私立医に行く同級生は一人もいませんでした。
本人には面と向かっては言わないけれど、俺くらいの年代の人間は、おそらくは８－９割は私立卒を今でも「何偉そうなこと抜かしてるんだ、この裏口バカが」と心の底で軽蔑し、嘲笑しているよ。当の本人には面と向かっては絶対にそんなことは言わないけどね。
805:卵の名無しさん
19/02/17 22:18:47.40 rglxpU9N.net
It is common knowledge among doctors and patients that Do-Teihen(exclusively bottom-leveled medical school) graduates mean morons who bought their way to Gachi'Ura(currently called by themselves)
According to the experience of entrance exam to medical school in the era of Showa, when the sense of discrimination against
privately-founded medical schools were more intense than it is now,
all such schools but for Keio had been so compared to some specialized institution for educable mentally retarded kids that nobody but imbecile successors of physicians in private practice had applied for admission.
There had been NOT a single classmate who chose willingly against his/her common sense to go to the Do-Teihen(exclusively bottom-leveled medical school, currently also known as Gachi'Ura),
which would have cost outrageous money and its graduates are destined to be called Uraguchi morons who bought thier way into the Do-Teihen, by thier colleagues and even by thier own clients.
Although people won't call them names to their face,
certain 80-90% people of about my age have been yet scorning and sneering at Uraguchi graduates, speaking in the back of our mind,
" Uraguchi morons shall not behave like somebody."
We never speak out face to face in real life.
806:卵の名無しさん
19/02/19 19:34:53.54 5dzWj+iV.net
俺は医者ではない。
807:卵の名無しさん
19/02/19 22:25:59.08 j+jAYShv.net
医学統計は苦手だ
できなきゃ論文読めんけど
808:卵の名無しさん
19/02/20 09:47:48.92 P0P24ylG.net
>>765
>456の問題点を指摘できないなら不正統計に騙される素養があるね。
809:卵の名無しさん
19/02/27 17:12:27.40 L7SdQ9d/.net
まあ落ち着いて
萌え動画でもどうぞ
URLﾘﾝｸ(m.youtube.com)
810:卵の名無しさん
19/03/01 20:33:40.05 CpB3LdCN.net
「もし君の妻子や友人にいつまでも生きていてもらいたいと思うならば、君は愚かである。というのは、君は、君の力の及びうる範囲内にないものまでも君の力の及びうる範囲内にあることを欲し、よそのものまでも君のものであることを欲しているからである。」
811:卵の名無しさん
19/03/04 13:04:58.69 PhmNbZK3.net
>>767 これでも見て落ちつけ https://www.youtube.com/watch?v=mB3ZVIc4EF4
813:卵の名無しさん
19/03/14 16:45:10.23 JBQs0yPH.net
韓国で大麻が医療解禁した
国際条約では臨床試験や学術研究(輸出輸入など)が認められているんで一度訪ねてみては?
URLﾘﾝｸ(livedoor.blogimg.jp)
URLﾘﾝｸ(livedoor.blogimg.jp)
韓国で使われてる大塚製薬の大麻製剤は日本の医師会員向けに和訳も開始されている。
URLﾘﾝｸ(news.nicovideo.jp)
814:卵の名無しさん
19/03/18 15:05:39.58 pTprsItO.net
まあ落ち着いて
萌え動画でもどうぞ
URLﾘﾝｸ(www.youtube.com)
815:卵の名無しさん
19/03/27 22:03:49.30 Kh2JcAUc.net
統計楽しいね
816:卵の名無しさん
19/03/28 00:45:03.92 8PkLLiIS.net
URLﾘﾝｸ(www3.nhk.or.jp)
インフルエンザの新しい治療薬「ゾフルーザ」を投与されたＡ香港型のインフルエンザ患者30人を調べたところ、
70％余りに当たる22人から、この薬が効きにくい耐性ウイルスが検出されたことが国立感染症研究所の調査で分かりました。
調査件数は多くないものの、専門家は現在のような使用を続けると、耐性ウイルスが広がるおそれがあるとして使用基準を見直すべきだと指摘しています。
URLﾘﾝｸ(www3.nhk.or.jp)
問題：耐性発生率の95%信頼区間は？ 👀
Rock54: Caution(BBR-MD5:1341adc37120578f18dba9451e6c8c3b)
817:卵の名無しさん
19/03/28 00:49:17.10 8PkLLiIS.net
95%CI 0.541 0.877なので過半数に耐性出現と言っていいな。
818:卵の名無しさん
19/03/28 01:30:14.34 uQgP69pk.net
耐性出現率が過半数である確率はいくらか？
819:卵の名無しさん
19/06/25 16:00:08.86 ZxTc5GC6.net
>>775
あんた医科歯科大卒だろ？
学部1～２回生の時に線形代数と微分積分はやった？
820:卵の名無しさん
19/06/25 16:01:39.44 ZxTc5GC6.net
線形代数と微分積分を知らんと統計学の勉強は無理だね。
まあ、上を見れば切りがないけどね。
821:卵の名無しさん
19/06/28 10:53:20.90 30gP48kS.net
>>777
Rが使えれば統計扱えるよ。
今どき分布表で補間値計算なんてしないし。
非心ｔ分布もRが出してくれる。
URLﾘﾝｸ(www3.nhk.or.jp)
インフルエンザの新しい治療薬「ゾフルーザ」を投与されたＡ香港型のインフルエンザ患者30人を調べたところ、70％余りに当たる22人から、この薬が効きにくい耐性ウイルスが検出されたことが国立感染症研究所の調査で分かりました。
調査件数は多くないものの、専門家は現在のような使用を続けると、耐性ウイルスが広がるおそれがあるとして使用基準を見直すべきだと指摘しています。
耐性化率が50%以上である確率は
pbeta(0.5,1+22,1+8,lower=F)
[1] 0.9946631
822:卵の名無しさん
19/06/28 12:35:11.06 30gP48kS.net
こんな質問があったからRを使って答えておいた。
他の人も同じ数値出していたな。
655 名無しさん＠１周年 sage 2019/06/28(金) 12:05:50.17 ID:vuCIc6fU0
偏差値72.5ってIQにするとどれくらいなの？
URLﾘﾝｸ(asahi.2ch.net)
823:卵の名無しさん
19/06/28 20:16:59.55 30gP48kS.net
統計学のための数学
あんまり公式を使うことないなぁ。
定積分も数値積分でやっちゃうから
URLﾘﾝｸ(www.data-arts.jp) <
824:卵の名無しさん
19/10/02 15:20:15.02 nv8GYZZt.net
kk=(88+96)/2/100 #TP: kessei kando
kt=(89+100)/2/100 #TN: kessei tokuido
pLHk=kk/(1-kt) #TP/FP
nLHk=(1-kk)/kt #(1-TP)/(1-FP)=FN/TN
bk=(90+98)/2/100 #ben kando
bt=(87+100)/2/100 #ben tokuido
pLHb=bk/(1-bt)
nLHb=(1-bk)/bt
prep=0.5 # pre-possibility
preo=prep/(1-prep) # pre-odds
posto=preo*nLHk*pLHb # kessei(-)&ben(+)
(postp=posto/(1+posto))
#
posto=preo*nLHk*nLHb # kessei(-)&ben(-)
(postp=posto/(1+posto))
825:卵の名無しさん
19/10/02 17:26:13.36 nv8GYZZt.net
#
# 尿素呼気試験(B) 90～100 80～99
# 血清抗体(S) 88～96 89～100
# 尿中抗体(U) 89～97 77～95
# 便中抗原(F) 90～98 87～100
U=c(90/2+100/2,80/2+99/2)/100
S=c(88/2+96/2,89/2+100/2)/100
U=c(89/2+97/2,77/2+95/2)/100
F=c(90/2+98/2,87/2+100/2)/100
DOR <- function(T){
TP=T[1] ; FN=1-TP
TN=T[2] ; FP=1-TN
pLR=TP/FP
nLR=FN/TN
pLR/nLR
}
DOR(U)
DOR(S)
DOR(U)
DOR(F)
826:卵の名無しさん
19/10/02 17:34:50.26 nv8GYZZt.net
B=c(90/2+100/2,80/2+99/2)/100
S=c(88/2+96/2,89/2+100/2)/100
U=c(89/2+97/2,77/2+95/2)/100
F=c(90/2+98/2,87/2+100/2)/100
DOR <- function(T){
TP=T[1] ; FN=1-TP
TN=T[2] ; FP=1-TN
pLR=TP/FP
nLR=FN/TN
pLR/nLR
}
> DOR(B)
[1] 161.9524
> DOR(S)
[1] 197.5909
> DOR(U)
[1] 81.61224
> DOR(F)
[1] 225.359
>
827:卵の名無しさん
19/11/25 07:37:28 ZvHlCyJI.net
# ゴルゴ13は100発100中
# ゴルゴ14は10発10中
# ゴルゴ15は1発1中
# とする。
# 各々10000発撃ったとき各ゴルゴの命中数の期待値はいくらか？
G13 <- function(N,n,r){
pp=0:N
f <- function(x) choose(x,r)*choose(N-x,n-r)/choose(N,n)
sum(pp*f(pp)/sum(f(pp)))
}
G13(10000,100,100)
G13(10000,10,10)
G13(10000,1,1)
#.n発r中の狙撃手が.N発狙撃するときの命中数を返す
Go13 <- function(.N, .n, r, k=10^3){ # k:シミュレーション回数
f <-function(S,N=.N,n=.n){
y=c(rep(1,S),rep(0,N-S))
sum(sample(y,n))
}
xx=r:.N
SS=NULL
for(i in 1:k){
x=sapply(xx,f)
SS=c(SS,which(x==r)-1+r)
}
print(summary(SS))
invisible(SS)
}
828:卵の名無しさん
19/11/27 20:32:52 s7+UeQPF.net
a=c(1,10,100)
curve(dbeta(x,1+a[1],1),bty='l')
curve(dbeta(x,1+a[2],1),bty='l')
curve(dbeta(x,1+a[3],1),bty='l')
1000*qbeta(0.95,a+1,1,lower=F)
829:卵の名無しさん
19/11/27 20:34:40 s7+UeQPF.net
> a=c(1,10,100)
> curve(dbeta(x,1+a[1],1),bty='l')
> curve(dbeta(x,1+a[2],1),bty='l')
> curve(dbeta(x,1+a[3],1),bty='l')
> 10000*qbeta(0.95,a+1,1,lower=F)
[1] 2236.068 7615.958 9707.748
>
830:卵の名無しさん
19/11/30 09:56:41 bijshzSl.net
# 安倍総理招待枠(60-****)が招待状から
# 254x 2409 357xだという。
# URLﾘﾝｸ(youtu.be)
# 通し番号がふられているとして
# 357xを3570と低めに見積もる。
#>いや、それどころか、今年はそれ以上の数字も確認されている。
#選挙で安倍首相が遊説をおこなうとかなりの頻度で目撃されている熱烈な支持者であるM氏という人物がいるのだが、
#氏が菅官房長官や杉田水脈議員とのツーショット写真などとともにSNSにアップしていた
# 今年の「桜を見る会」の受付票には「60-4908」とナンバリングされているのだ。
#ｽﾚﾘﾝｸ(seijinewsplus板)
# 他�
831:ﾉ今年は254? 2409 357? という番号が国会で発言されている。 #参加者の18200の半数9100が官邸枠の可能性があると仮定して官邸枠招待者数の期待値と95%CIを算出せよ。 Nmin=4908 Nmax=9100 n=Nmin:Nmax m=4 pmf=choose(Nmax-1,m-1)/choose(n,m) #Pr(max=60|n) plot(n,pmf,ylab='probability') pdf=pmf/sum (pmf) plot(n,pdf,ylab='density',bty='l') sum(n*pdf) #E(n) plot(n,cumsum(pdf)) abline(h=0.95,lty=3) idx=min(which(cumsum(pdf)>0.95)) n[idx]
832:卵の名無しさん
19/11/30 10:28:03.80 bijshzSl.net
> Nmin=4908
> Nmax=9100
> n=Nmin:Nmax
> m=4
> pmf=choose(Nmax-1,m-1)/choose(n,m) #Pr(max=60|n)
> plot(n,pmf,ylab='probability')
> pdf=pmf/sum (pmf)
> plot(n,pdf,ylab='density',bty='l')
> sum(n*pdf) #E(n)
[1] 6191.377
> plot(n,cumsum(pdf))
> abline(h=0.95,lty=3)
> idx=min(which(cumsum(pdf)>0.95))
> n[idx]
[1] 8406
>
833:卵の名無しさん
19/12/09 08:25:56.61 zd3bijjh.net
# サイコロ
# 正６面体のサイコロがある．４面は青色、２面は赤色である．
# このサイコロを合計20回振るとき、最も起こりそうな順番はどれか？
# １．赤 青 赤 赤 赤
# ２．青 赤 青 赤 赤 赤
# ３．青 赤 赤 赤 赤 赤
sim <- function(){
a=sample(0:1,20, replace=TRUE, prob=c(4,2))
b=as.character(a)
c=paste(b,collapse="")
s1=paste(c(1,0,1,1,1),collapse="")
s2=paste(c(0,1,0,1,1,1),collapse="")
s3=paste(c(0,1,1,1,1,1),collapse="")
res=c(grepl(s1,c),grepl(s2,c), grepl(s3,c))
return(res)
}
k=1e5
re=replicate(k,sim())
mean(re[1,])
mean(re[2,])
mean(re[3,])
834:卵の名無しさん
19/12/10 08:45:53 k/8kaoYw.net
seqNp <- function(N=100,K=5,p=0.5){
if(p==0) return(0)
if(N==K) return(p^K)
q=1-p
a=numeric(N) # a(n)=P0(n)/p^n , P0(n)=a(n)*p^n
for(i in 1:K) a[i]=q/p^i # P0(i)=q
for(i in K:(N-1)){ # recursive formula
a[i+1]=0
for(j in 0:(K-1)){
a[i+1]=(a[i+1]+a[i-j])
}
a[i+1]=q/p*a[i+1]
}
P0=numeric(N)
for(i in 1:N) P0[i]=a[i]*p^i # P0(n)=a(n)*p^n
MP=matrix(rep(NA,N*K),ncol=K)
colnames(MP)=paste0('P',0:(K-1))
MP[,'P0']=P0
MP[1,'P1']=p
for(i in (K-2):K) MP[1,i]=0
for(k in 2:K){
for(i in 1:(N-1)) MP[i+1,k]=p*MP[i,k-1]
} # Pk(n+1)=p*P(k-1)(n)
ret=1-apply(MP,1,sum)
ret[N]
}
835:卵の名無しさん
19/12/10 08:46:19 k/8kaoYw.net
seqNpJ <- function(N,K,p) seqNp(N,K,p)-seqNp(N,K+1,p)
seqNpJ(100,5,0.5)
seq2pCI <- function(N,K,alpha=0.05,Print=T){
vp=Vectorize(function(p)seqNp(N,K,p)-seqNp(N,K+1,p))
if(Print){curve(vp(x),lwd=2,bty='l',xlab='Pr[head]',ylab=paste('Pr[max',K,'-head repetition]'))
abline(h=alpha,lty=3)}
peak=optimize(vp,c(0,1),maximum=TRUE)$maximum
mean=integrate(function(x)x*vp(x),0,1)$value/integrate(function(x)vp(x),0,1)$value
lwr=uniroot(function(x,u0=alpha) vp(x)-u0,c(0.01,peak))$root
upr=uniroot(function(x,u0=alpha) vp(x)-u0,c(peak,0.99))$root
c(lower=lwr,mean=mean,mode=peak,upper=upr)
}
seq2pCI(100,5,0.05,T)
vs=Vectorize(function(K)seq2pCI(N=100,K,alpha=0.05,Print=F))
x=2:23
y=vs(x)
head(y)
y=y*100
plot(x,y['mean',],bty='l',pch=19,ylim=c(0,100),
xlab="最大連続数",ylab="推定裏口入学者数")
points(2:23,y['mode',],bty='l')
segments(x,y['lower',],x,y['upper',])
legend('right',bty='n',pch=c(19,1),legend=c("期待値","最頻値"))
836:卵の名無しさん
19/12/10 08:46:54 k/8kaoYw.net
# pdfからcdfの逆関数を作ってhdiを表示させて逆関数を返す
# 0,1での演算を回避 ∫[1/nxxx,1-1/nxx]dxで計算
pdf2HDI <- function(pdf,xMIN=0,xMAX=1,cred=0.95,nxx=1001){
xx=seq(xMIN,xMAX,length=nxx)
xx=xx[-nxx]
xx=xx[-1]
xmin=xx[1]
xmax=xx[nxx-2]
AUC=integrate(pdf,xmin,xmax)$value
PDF=function(x)pdf(x)/AUC
cdf <- function(x) integrate(PDF,xmin,x)$value
ICDF <- function(x) uniroot(function(y) cdf(y)-x,c(xmin,xmax))$root
hdi=HDInterval::hdi(ICDF,credMass=cred)
c(hdi[1],hdi[2])
}
# N試行で最大K回連続成功→成功確率ｐの期待値、最頻値と95%HDI
# max K out of N-trial to probability & CI
mKoN2pCI <- function(N=100 , K=4 , conf.level=0.95){
pmf=Vectorize(function(p)seqNp(N,K,p)-seqNp(N,K+1,p))
mode=optimize(pmf,c(0,1),maximum=TRUE)$maximum
auc=integrate(pmf,0,1)$value
pdf=function(x) pmf(x)/auc
mean=integrate(function(x)x*pdf(x),0,1)$value
curve(pdf(x),lwd=3,bty='l',xlab='Pr[head]',ylab='density')
lu=pdf2HDI(pdf,cred=conf.level)
curve(pdf(x),lu[1],lu[2],type='h',col='lightblue',add=T)
c(lu[1],mean=mean,mode=mode,lu[2])
}
mKoN2pCI(100,4) # one minute plz
mKoN2pCI(100,5)
837:卵の名無しさん
19/12/10 08:47:01 k/8kaoYw.net
# N試行で最大K回連続成功→成功確率ｐの期待値、最頻値と95% Quantile
# max K out of N-trial to probability & CI_quantile
mKoN2pCIq <- function(N=100 , K=4 , alpha=0.05){
pmf=Vectorize(function(p)seqNp(N,K,p)-seqNp(N,K+1,p))
mode=optimize(pmf,c(0,1),maximum=TRUE)$maximum
auc=integrate(pmf,0,1)$value
pdf=function(x) pmf(x)/auc
curve(pdf(x),bty='l')
mean=integrate(function(x)x*pdf(x),0,1)$value
cdf=function(x) MASS::area(pdf,0,x)
vcdf=Vectorize(cdf)
lwr=uniroot(function(x)vcdf(x)-alpha/2,c(0,mode))$root
upr=uniroot(function(x)vcdf(x)-(1-alpha/2),c(mode,1))$root
c(lower=lwr,mean=mean,mode=mode,upper=upr)
}
mKoN2pCIq(100,4)
mKoN2pCIq(100,5)
838:卵の名無しさん
19/12/10 08:48:13 2lYn70em.net
## simulation
mKoN2pCIga<-function(N=100,K=4,Print=T){
pmf=Vectorize(function(p)seqNp(N,K,p)-seqNp(N,K+1,p))
auc=integrate(pmf,0,1)$value
pdf=function(x) pmf(x)/auc
(mode=optimize(pmf,c(0,1),maximum=TRUE)$maximum)
(mean=integrate(function(x)x*pdf(x),0,1)$value)
(var=integrate(function(x)(x-mean)^2*pdf(x),0,1)$value)
(sd=sqrt(var))
#正規分布で近似してみる
if(Print){
curve(pdf,bty='l',col='red')
curve(dnorm(x,mean,sd),add=T,col='blue')}
c(lower=qnorm(0.025,mean,sd),mode=mode,mean=mean,upper=qnorm(0.975,mean,sd))
}
mKoN2pCIga(100,5)
vmK=Vectorize(function(K) mKoN2pCIga(100,K))
(y=cbind(0,vmK(2:30)))
plot(1:30,y['mean',],pch=19)
points(1:30,y['mode',])
839:卵の名無しさん
19/12/17 22:26:20 MXFfIsV3.net
# 10円硬貨を投げて5回連続で表が出ると
# 稀な現象（0.5の5乗は 0.03125)なので表の出る確率は0.5ではない、というのが頻度主義統計の結論。
# この硬貨はどの程度歪（いびつ）なのかを推論するのがベイズ推論の面白さである。
# この硬貨で表がでる確率の95%信頼区間を計算してみよう。
#
# 表のでる確率分布をβ分布β(a,b)でa=bとすれば、0.5を頂点とする対象な事前確率分布となる。
# 事後分布はβ(a+5,a+0)となるのだが、問題はパラメータaをいくらに設定するかである。
# aを変数として95％信頼区間下限値をグラフにすると最小値があることがわかる。
# URLﾘﾝｸ(i.imgur.com)
# 最尤法で算出させるとa=11.25のとき0.406となる。
# あとは95%区間を計算して、この硬貨が表が出る確率の95%信頼区間は0.406 0.763となる。
# URLﾘﾝｸ(i.imgur.com)
HDInterval::hdi(qbeta,shape1=1+5,shape2=1)
840:卵の名無しさん
19/12/18 12:46:29 HkAD959b.net
n=5
l <- function(a) HDInterval::hdi(qbeta,shape1=a+n,shape2=a)['lower']
vl=Vectorize(l)
curve(vl(x),0.5,1000,bty='l',lwd=2)
(lo=optimise(vl,(c(0.5,1000))))
HDInterval::hdi(qbeta,shape1=lo$min + n,shape2=lo$min)
841:卵の名無しさん
19/12/18 13:45:30 HkAD959b.net
コインをなげたら５回続けて表が出た。
表が出る確率は、平均値が0.5のベータ分布(一様分布もベータ分布の一つ）に従うように鋳造されているとする。
このコインの表の出る確率を95%信頼区間で推定するとき
区域下限値の最小値はいくらか？
つまり、95%以上の信頼性で表の出る確率はいくつ以上�
842:ﾆ言えるか？
843:卵の名無しさん
19/12/18 16:26:08.99 j32hgbYi.net
rm(list=ls())
n=5
l <- function(a) HDInterval::hdi(qbeta,shape1=a+n,shape2=a)['lower']
vl=Vectorize(l)
curve(vl(x),0.5,100,bty='l',lwd=2)
(lo=optimise(vl,(c(0.5,1000))))
HDInterval::hdi(qbeta,shape1=lo$min + n,shape2=lo$min)
#
binom::binom.confint(5,5)
f <- function(x) binom::binom.bayes(5,5,prior.shape1=x,prior.shape2=x)$lower
curve(f(x),0.5,50,bty='l',lwd=2)
optimize(f,c(0.5,50))
844:卵の名無しさん
19/12/18 19:35:43.82 HkAD959b.net
>>798
95%CIをグラフ化
URLﾘﾝｸ(i.imgur.com)
845:卵の名無しさん
19/12/18 19:59:58.70 HkAD959b.net
一般化
#
"コインをなげたらn回続けて表が出た。
表が出る確率は、平均値が0.5のベータ分布(一様分布もベータ分布の一つ）に従うように鋳造されているとする。
このコインの表の出る確率を95%信頼区間で推定するとき区域下限値の最小値はいくらか？
つまり、95%以上の信頼性で表の出る確率はいくつを下回らないと言えるか？"
seq2lowest <- function(n){
l <- function(a) HDInterval::hdi(qbeta,shape1=a+n,shape2=a)['lower']
vl=Vectorize(l)
opt=optimise(vl,(c(0.5,1000)))
hdi=HDInterval::hdi(qbeta,shape1=lo$min + n,shape2=lo$min)
c('shape'=opt$minimum,opt$objective,hdi['upper'])
}
seq2lowest(5)
graphics.off()
nn=2:50
plot(nn,sapply(nn,seq2lowest)['lower',],pch=19,bty='l')
846:卵の名無しさん
19/12/18 22:08:27 RWipilhS.net
よくわからないけど面白そう
847:卵の名無しさん
19/12/19 07:45:22.67 mtJFYvZ2.net
rm(list=ls())
gambling <- function(A=2,B=1,w=3,p=0.5,k=1e5){ # k : how many simulate
sim <- function(){
while(A < w & B < w){
g = rbinom(1,1,p)
if(g==1){
A=A+1
}else{
B=B+1
}
}
A > B
}
mean(replicate(k,sim())) # Pr[A wins]
}
gambling(1,0,4) # 日本シリーズ
848:卵の名無しさん
19/12/19 07:45:58.37 mtJFYvZ2.net
# 日本シリーズ
f <- function(num, N=2, digit = 7){ # winner sequence of 7 games
r=num%%N
q=num%/%N
while(q > 0 | digit > 1){
r=append(q%%N,r)
q=q%/%N
digit=digit-1
}
return(r)
}
vf=Vectorize(f)
dat=t(vf(0:(2^7-1))) # matrirx of winner sequence of 7 games
(dat1=dat[dat[,1]==1,-1]) # winner sequence of 6 games
g <- function(s,A=1,B=0,w=4){ # A beats B T/F?
i=0
while(A < w & B < w){
i=i+1
r = s[i]
if(r==1){
A=A+1
}else{
B=B+1
}
}
A > B
}
sum(apply(dat1,1,g))/nrow(dat1)
849:卵の名無しさん
19/12/19 08:17:50.32 mtJFYvZ2.net
# 通算勝率に従って次の試合に勝つ確率が決まるという設定
# 先勝したAの日本シリーズ前の対戦成績がa勝ｂ負とする
rm(list=ls())
N_series <- function(A=1,B=0,w=4,a=10,b=10,k=1e5){
sim <- function(){
p=(A+a)/(A+B+a+b)
while(A < w & B < w){
g = rbinom(1,1,p)
if(g==1){
A=A+1
p=(A+a+1)/(A+B+a+b+1)
}else{
B=B+1
p=(A+a)/(A+B+a+b+1)
}
}
A > B
}
mean(replicate(k,sim())) # Pr[A wins]
}
N_series()
850:卵の名無しさん
19/12/19 09:38:09.07 +FZyQ+1V.net
互角と言っても1勝1敗と1000勝1000敗だ結果が変わるな
> rm(list=ls())
> N_series <- function(A=1,B=0,w=4,a=10,b=10,k=1e5){
+ sim <- function(){
+ p=(A+a)/(A+B+a+b)
+ while(A < w & B < w){
+ g = rbinom(1,1,p)
+ if(g==1){
+ A=A+1
+ p=(A+a+1)/(A+B+a+b+1)
+ }else{
+ B=B+1
+ p=(A+a)/(A+B+a+b+1)
+ }
+ }
+ A > B
+ }
+ mean(replicate(k,sim())) # Pr[A wins]
+ }
>
> N_series(a=1,b=1)
[1] 0.76425
> N_series(a=10,b=10)
[1] 0.67035
> N_series(a=100,b=100)
[1] 0.65629
> N_series(a=1000,b=1000)
[1] 0.65802
851:卵の名無しさん
19/12/19 10:00:45.03 +FZyQ+1V.net
> N_series(a=1,b=1)
[1] 0.76266
> N_series(a=10,b=10)
[1] 0.6716
> N_series(a=100,b=100)
[1] 0.65707
> N_series(a=1000,b=1000)
[1] 0.65695
>
852:卵の名無しさん
19/12/19 16:11:38.16 mtJFYvZ2.net
NS <- function(w,A,B,p){ # 先にw点得点した方が勝者、A,B：現在の点数 ,p:甲の勝率
ans=0
for(k in 0:(w-B-1)){ # k: Aが勝者になるまでのBの点数
ans=ans+choose(w-A-1+k,w-A-1)*(1-p)^k*p^(w-A)
}
return(ans)
}
NS(3,2,1,0.5)
853:卵の名無しさん
19/12/21 11:16:13.81 JyjCwOQx.net
日本シリーズで賭けをする。
日本シリーズは先に4勝したチームが優勝。
勝率はそれまでの引き分けを除いた通算勝率に従うとする。
開始前の通算成績はA:2勝、B:4勝であった。
今、シリーズでAが先勝（第一試合に勝利）した。
この時点でA,Bのどちらに賭ける方が有利か？
854:卵の名無しさん
19/12/21 21:55:53.22 NW9TiMNu.net
#
855:小数点付きの数numをN進法で表示する rm(list=ls()) dec2basen <- function(num, N, kmin = 5){ # kmin：最小小数点後桁 int=floor(num) r=int%%N q=int%/%N while(q > 0){ r=append(q%%N,r) q=q%/%N } # rに整数のN進法表示を格納 x=num-floor(num) k=max(nchar(x)-2,kmin) # 同長もしくはkminの長さの小数表示 a=numeric(k) for(i in 1:k){ y=x*N a[i]=floor(y) x=y-a[i] # r . a[1] a[2] a[3] ... a[k] } if(N<=10){ # Nが10以下は数値として表示 cat(r,paste(".",paste(a,collapse = ''),sep=''),'\n',sep='') } else{ # Nが1１以上は整数部分と小数部分を数列で表示 print(list(int=r,decimal=a)) } invisible() }
856:卵の名無しさん
19/12/21 23:15:54.59 NW9TiMNu.net
戻り値がある方がいいな。
# 小数点付きの数numをN進法で表示する
rm(list=ls())
dec2basen <- function(num, N, kmin = 5){ # kmin：最小小数点後桁
int=floor(num)
r=int%%N
q=int%/%N
while(q > 0){
r=append(q%%N,r)
q=q%/%N
} # rに整数のN進法表示を格納
x=num-floor(num)
k=max(nchar(x)-2,kmin) # 同長もしくはkminの長さの小数表示
a=numeric(k)
for(i in 1:k){
y=x*N
a[i]=floor(y)
x=y-a[i] # r . a[1] a[2] a[3] ... a[k]
}
b=list(int=r,decimal=a)
if(N<=10){ # Nが10以下は数値として表示
cat(r,paste(".",paste(a,collapse = ''),sep=''),'\n',sep='')
}
else{ # Nが1１以上は整数部分と小数部分を数列で表示
print(b)
}
invisible(b)
}
857:卵の名無しさん
19/12/21 23:26:15.77 NW9TiMNu.net
> dec2basen(5/9,3)
0.120000000000000
>
> dec2basen(3.14159265359,7)
3.06636514320
> dec2basen(3.14159265359,8)
3.11037552421
> dec2basen(3.14159265359,9)
3.12418812407
> dec2basen(3.14159265359,14)
$int
[1] 3
$decimal
[1] 1 13 10 7 5 12 13 10 8 1 4
>
> re=dec2basen(sqrt(3),16) # √３の１６進法表示
$int
[1] 1
$decimal
[1] 11 11 6 7 10 14 8 5 8 4 12 10 10 0 0
> b=re$decimal
> cat(re$int,'.',paste0(c('0':'9','A','B','C','D','E','F')[b+1],sep='',collapse=''),sep='')
1.BB67AE8584CAA00
858:卵の名無しさん
19/12/22 13:49:45.74 D+p3chog.net
# 小数点付きの数numをN進法で表示する
floor((1.2-1)*5) != floor(1)
rm(list=ls())
decN <- function(num, N, kmin = 5){ # kmin：最小小数点後桁
int=floor(num)
r=int%%N
q=int%/%N
while(q > 0){
r=append(q%%N,r)
q=q%/%N
} # rに整数のN進法表示を格納
x=num-floor(num)
k=max(nchar(x)-2,kmin) # 同長もしくはkminの長さの小数表示
a=numeric(k)
for(i in 1:k){
y=x*N
a[i]=floor(y+2e-13) # n=1:9 ;trunc((1+1/n-1)*n) ;trunc((1+1/n-1)*n+5e-14)対応
x=y-a[i] # r . a[1] a[2] a[3] ... a[k]
}
b=list(integer=r,decimal=a,num=sum(c(int,a)*(1/N)^(0:k)))
fig=c(0:9,letters)[1:N]
if(N<=36){ # Nが36以下は数値として表示
cat(b$integer,'.',paste(fig[b$decimal+1],sep='',collapse=''),sep='')
cat('\n')
}
else{ # Nが1１以上は整数部分と小数部分を数列で表示
print(b[1:2])
}
invisible(b)
}
859:卵の名無しさん
19/12/22 13:50:03.28 D+p3chog.net
> decN(5/9,3)
0.120000000000000
> decN(0.728,5,10)
0.331004444
> decN(3.14159265359,7)
3.06636514320
> decN(3.14159265359,14)
3.1da75cda814
860:卵の名無しさん
19/12/22 14:24:46.39 D+p3chog.net
# 小数点付きの数numをN進法で表示する
rm(list=ls())
decN <- function(num, N, kmin = 5){ # kmin：最小小数点後桁
int=floor(num)
r=int%%N
q=int%/%N
while(q > 0){
r=append(q%%N,r)
q=q%/%N
} # rに整数のN進法表示を格納
x=num-floor(num)
k=max(nchar(x)-2,kmin) # 同長もしくはkminの長さの小数表示
a=numeric(k)
for(i in 1:k){
y=x*N
a[i]=trunc(y+2e-13) # n=1:9 ;trunc((1+1/n-1)*n) ;trunc((1+1/n-1)*n+5e-14)対応
x=y-a[i] # r . a[1] a[2] a[3] ... a[k]
}
b=list(integer=r,decimal=a,num=sum(c(int,a)*(1/N)^(0:k)))
fig=c(0:9,letters)[1:N]
if(N<=36){ # Nが36以下は数値として表示
cat(b$integer,'.',paste(fig[b$decimal+1],sep='',collapse=''),sep='')
cat('\n')
}
else{ # Nが37以上は整数部分と小数部分を数列で表示
print(b[1:2])
}
invisible(b) # b$num:検証用
}
861:卵の名無しさん
19/12/22 16:52:54 D+p3chog.net
# 小数点付きの数numをN進法で表示する
rm(list=ls())
decN <- function(num, N, kmin = 5){ # kmin：最小小数点後桁
int=floor(num)
r=int%%N
q=int%/%N
while(q > 0){
r=append(q%%N,r)
q=q%/%N
} # rに整数のN進法表示数列を格納
k=max(nchar(num)-nchar(floor(num))-1,kmin) # 同長もしくはkminの長さの小数表示
a=numeric(k)
x=round(num-floor(num),k) # e.g. 7.28-floor(7.28)-0.28 != 0に対応
for(i in 1:k){
y=round(x*N,k) # e.g. 0.728*5-3.64 !=0 に対応
a[i]=floor(y)
x=y-a[i] # r . a[1] a[2] a[3] ... a[k]
}
b=list(integer=r,decimal=a,num=sum(c(int,a)*(1/N)^(0:k)))
fig=c(0:9,letters)[1:N]
if(N<=36){ # Nが36以下は数値として表示
cat(paste(fig[b$integer+1],sep='',collapse=''),
'.',paste(fig[b$decimal+1],sep='',collapse=''),sep='')
cat('\n')
}
else{ # Nが37以上は整数部分と小数部分を数列で表示
print(b[1:2])
}
invisible(b) # b$num:検証用
}
862:卵の名無しさん
19/12/22 23:14:45 D+p3chog.net
rm(list=ls())
NS <- function(w,A,B,p){ # 先にw点得点した方が勝者、A,B：現在の点数 ,p:甲の勝率
ans=0
for(k in 0:(w-B-1)){ # k: Aが勝者になるまでのBの点数
ans=ans+choose(w-A-1+k,w-A-1)*(1-p)^k*p^(w-A)
}
return(ans)
}
NS(3,2,1,0.5)
NS(4,1,0,0.5)
curve(NS(4,1,0,x),bty='l',lwd=2)
abline(h=0.5,col='gray')
# 日本シリーズで先勝したときに優勝する確率がx以上であるために必要な勝率vf(x)は？
f <- function(u0) uniroot(function(p) NS(4,1,0,p)-u0,c(0,1))$root
f(0.5)
vf=Vectorize(f)
curve(vf(x),bty='l')
vf(seq(0,1,0.1))
863:卵の名無しさん
19/12/23 13:09:52 WZg+GXSV.net
デバッグした
# 通算勝率に従って次の試合に勝つ確率が決まるという設定
# 先勝したAの日本シリーズ前の対戦成績がa勝ｂ負とする
rm(list=ls())
N_series <- function(A=1,B=0,w=4,a=2,b=4,k=1e5){
sim <- function(){
while(A < w & B < w){
p=(A+a)/(A+B+a+b)
g = rbinom(1,1,p)
if(g==1){
A=A+1
}else{
B=B+1
}
}
A > B
}
mean(replicate(k,sim())) # Pr[A wins]
}
N_series()
864:卵の名無しさん
19/12/23 15:14:31.94 +ZpAQitI.net
日本シリーズは先に4勝したチームが優勝。
勝率はそれまでの通算勝率に従うとする。引き分けはないものとする。
勝負がつくごとに次回の勝率が変化する。
（第二試合にAが勝つ確率は通算勝率の3/7
865:Aが勝ったら第三試合に勝つ確率は4/8 Aが負けたら第三試合に勝つ確率は3/8 になるという設定） シリーズ開始前の通算成績はA:2勝、B:4勝であった。 今シリーズでAが先勝（第一試合に勝利）した。 この時点でどちらが優勝するか賭けをする。 A,Bのどちらに賭ける方が有利か？"
866:卵の名無しさん
19/12/26 13:19:39 Andgk16C.net
poisson.test(c(11, 6+8+7), c(800, 1083+1050+878))
prop.test(c(11, 6+8+7), c(800, 1083+1050+878))
867:卵の名無しさん
19/12/26 13:47:03 Andgk16C.net
poisson.test(c(11, 6+8+7), c(800, 1083+1050+878))
prop.test(c(11, 6+8+7), c(800, 1083+1050+878))
library(fmsb)
ratedifference(11, 6+8+7, 800, 1083+1050+878)
rateratio(11, 6+8+7, 800, 1083+1050+878)
source('tools.R')
jags4prop(11, 6+8+7, 800, 1083+1050+878)
868:卵の名無しさん
19/12/27 17:58:44 pKcbCDhI.net
ZooG <- function(L=1,M=1,N=1,A=30,B=60){
alpha=B/180*pi
beta=A/180*pi
a=M/N # length of 0P
b=L/N # length of PQ
P=a*(cos(alpha)+1i*sin(alpha))
Q=P+b*(cos(alpha+pi-beta)+1i*sin(alpha+pi-beta))
(pi-Arg(Q-1))/pi*180 # degree of Q-I-0
max(Arg(Q-P)/pi*180 - Arg(Q-1)/pi*180, Arg(Q-P)/pi*180 - Arg(Q-1)/pi*180+360) # degree of P-Q-I
abs(Q-1) # length of QI
xlim=ylim=max(a,b)*c(-1,2)
plot(0,xlim=xlim,ylim=ylim,bty='l',type='l',axes=FALSE, ann=FALSE,asp=1)
# draw segment of complex a to complex b
seg <- function(a,b,color=2,...) segments(Re(a),Im(a),Re(b),Im(b),col=color,...)
# draw text y at complex x
pt <- function(x,y=NULL,...) text(Re(x),Im(x), ifelse(is.null(y),'+',y), ...)
seg(0,1) ; seg(0,P) ; seg(P,Q)
pt(0,paste('O',B,'°')) ; pt(1,'I') ; pt(P,paste('P ',A,'°')) ; pt(Q,'Q')
pt((P+Q)/2, paste('L = ',L)) ; pt((P+0)/2, paste('M = ',M)) ; pt(1/2, paste('N = ',N))
seg(1,Q,col=1,lty=2)
angleI=(pi-Arg(Q-1))/pi*180 # degree of Q-1-0
if(angleI<0)angleI=angleI+360
if(angleI>360)angleI=angleI-360
if(angleI>180 & angleI<360) angleI=360-angleI
angleQ=Arg(Q-P)/pi*180 - Arg(Q-1)/pi*180 # degree of P-Q-1
if(angleQ<0)angleQ=angleQ+360
if(angleQ>360)angleQ=angleQ-360
if(angleQ>180 & angleQ<360) angleQ=360-angleQ
length=abs(Q-1)*N # length of Q1
c(angleQ=angleQ,angleI=angleI,QIlength=length)
}
869:卵の名無しさん
19/12/27 17:59:16 pKcbCDhI.net
# 長さL,M,NのZ尺を角度A°（LとMのなす角）、B°（MとNのなす角）で折り曲げたとき
# 先端と終端を結ぶ線とZ尺の作る角度および先端と終端の距離
ZooG(1,1,1,36,24)
ZooG(sample(10,1),sample(10,1),sample(10,1),sample(179,1),sample(179,1))
870:卵の名無しさん
19/12/31 20:45:37.19 caKR9jL7.net
# The ambulance arrives X hours ealier at hospital,
# when the patient leave clinic s hours earlier than planned
# and encounter the ambulance t hours later,
# ambulance runs with the velocity of v0 without patient, and v1 with patient
# clinic car runs with the velocity of w
Earlier <- # s hour earlier departure, X hour earlier arrival
function(s=NULL,X=NULL,t=NULL,v0=60,v1=45,w=30){
if(is.null(s)) re=c(s=X/(v0/(v0+w)*w*(1/v0+1/v1)),t=X/w/(1/v0+1/v1))
if(is.null(X)) re=c(X=s*v0/(v0+w)*w*(1/v0+1/v1),t=s*v0/(v0+w))
if(!is.null(t)) re=c(s=t + t*w/v0, X = t*w*(1/v0+1/v1))
return(re)
}
Earlier(X=10/60)*60
Earlier(s=30/60)*60
Earlier(t=30/60)*60
871:卵の名無しさん
20/01/24 09:01:51.18 p9O7hNNU.net
TE=cbind(rep(0:1,c(3,3)),rep(0:2,2))
colnames(TE)=c('JD','service')
TE
cond1 <- function(P,Q) !(P&!Q)
cond2 <- function(P,Q) cond1(P,Q) & cond1(!P,!Q)
f <- function(x){
all(c(
x[1]==1 & cond2(cond2(x[1]==1,x[2]!=1),x[2]>0) | x[1]==0 & !cond2(cond2(x[1]==1,x[2]!=1),x[2]>0)
))
}
TE[(apply(TE,1,f)),]
872:卵の名無しさん
20/01/24
873:10:26:25.09 ID:p9O7hNNU.net
874:卵の名無しさん
20/01/24 10:40:22.42 9ccIONJj.net
TE=cbind(rep(0:1,c(3,3)),rep(0:2,2))
colnames(TE)=c('JD','service')
TE
cond1 <- function(P,Q) !(P&!Q)
cond2 <- function(P,Q) P==Q
f <- function(x){
all(c(
x[1]==1 & cond2(cond2(x[1]==1,x[2]!=1),x[2]>0) | x[1]==0 & !cond2(cond2(x[1]==1,x[2]!=1),x[2]>0)
))
}
TE[(apply(TE,1,f)),]
875:卵の名無しさん
20/01/24 10:42:15.05 9ccIONJj.net
> TE=cbind(rep(0:1,c(3,3)),rep(0:2,2))
> colnames(TE)=c('JD','service')
> TE
JD service
[1,] 0 0
[2,] 0 1
[3,] 0 2
[4,] 1 0
[5,] 1 1
[6,] 1 2
> cond1 <- function(P,Q) !(P&!Q)
> cond2 <- function(P,Q) P==Q
> f <- function(x){
+ all(c(
+ x[1]==1 & cond2(cond2(x[1]==1,x[2]!=1),x[2]>0) | x[1]==0 & !cond2(cond
<(x[1]==1,x[2]!=1),x[2]>0) | x[1]==0 & !cond2(cond2 (x[1]==1,x[2]!=1),x[2]>0)
+ ))
+ }
> TE[(apply(TE,1,f)),]
JD service
[1,] 0 2
[2,] 1 2
>
URLﾘﾝｸ(o.5ch.net)
876:卵の名無しさん
20/01/24 15:17:02.19 p9O7hNNU.net
"某女子大には決して嘘をつかない女子大生と必ず嘘をつく女子大生がいることがわかっている。
この女子大の学生（嘘つきかどうかは不明）から
「あなたのいうことが正しければ手コキかフェラをしてあげる、間違っていれば何もしてあげない」と言われた。
女子大生にフェラをしてもらうには何と言えばいいか？
"
rm(list=ls())
TE=cbind(rep(0:1,c(4,4)),rep(0:3,2))
colnames(TE)=c('JD','service')
TE # JD 0:嘘つき 1:正直, service 0:何もしない 1:手コキ 2:フェラ 3:手コキ＋フェラ
cond1 <- function(P,Q) !(P&!Q) # ｐならばQ
cond2 <- function(P,Q) cond1(P,Q) & cond1(!P,!Q) # （ｐならばQ） かつ（ｐでないならQでない）
f <- function(x){
x[1]==1 & cond2(cond2(x[1]==1,x[2]!=1),x[2]>0) | x[1]==0 & !cond2(cond2(x[1]==1,x[2]!=1),x[2]>0)
}
TE[(apply(TE,1,f)),]
cond3 <- function(P,Q) P==Q
f <- function(x){
x[1]==1 & cond3(cond3(x[1]==1,x[2]!=1),x[2]>0) | x[1]==0 & !cond3(cond3(x[1]==1,x[2]!=1),x[2]>0)
}
TE[(apply(TE,1,f)),]
877:卵の名無しさん
20/01/25 07:15:52.06 v8nYWWGU.net
TEnr <- function(n,r,one=1){ # n(=5),r(=3)を指定して 0 0 1 1 1から1 1 1 0 0までの順列行列を作る
f=function(x){
re=numeric(n) # 容れ子
re[x]=one # 指定のindexにoneを代入
re
}
t(combn(n,r,f)) # oneを入れる組み合わせに上記関数fを実行して転置
}
878:卵の名無しさん
20/01/25 07:43:20.46 v8nYWWGU.net
0 1 以外の場合もあるからこうした方がいいな
TEnr <- function(n,r,zero=0,one=1){ # n(=5),r(=3)を指定して 0 0 1 1 1から1 1 1 0 0までの順列行列を作る
f=function(x){
re=rep(zero,n) # 容れ子
re[x]=one # 指定のindexにoneを代入
re
}
t(combn(n,r,f)) # oneを入れる組み合わせに上記関数fを実行して転置
}
879:卵の名無しさん
20/01/25 13:21:20 v8nYWWGU.net
"
７段の階段を上るさい、１段上るか、２段上るかのいずれかを組合わせて上るこ
ととする。上り方は何通りあるか。
"
# Σxi=n xi<=r 一度にｒ(=3)段まで上がれる時、総数ｎ(=7)段上がる場合の数
stairway <- function(n,r){
count=numeric(n)
m=ceiling(n/r)
for(i in m:n){
pi=gtools::permutations(r,i,rep=T)
pn=pi[apply(pi,1,sum)==n,]
count[i]=ifelse(is.vector(pn),1,nrow(pn))
}
sum(count)
}
stairway(7,2)
880:卵の名無しさん
20/01/27 15:33:10 54kq03+D.net
ここまで複雑だと推論での結論はかなり困難だと思う。
AからHの8人はそれぞれ正直者か嘘つきであり、誰が正直者か嘘つきかはお互いに知っている。
A,B,C,D,Eは嘘つきなら必ず嘘をつくが、F,G,Hは嘘つきでも正しいことを言う場合がある。
次の証言から誰を確実に正直者と断定できるか？
A「嘘つきの方が正直者より多い」
B「Hは嘘つきである」
C「Bは嘘つきである」
D「CもFも嘘つきである」
E「8人の中に、少なくとも1人嘘つきがいる」
F「8人の中に、少なくとも2人嘘つきがいる」
G「Eは嘘つきである」
H「AもFも正直者である」
プログラム解で解くのはたやすい。
TE=expand.grid(0:1,0:1,0:1,0:1,0:1,0:1,0:1,0:1)
colnames(TE)=LETTERS[1:8]
f <- function(x){
all(c(
x[1]==1 & sum(x==0)>sum(x==1) | x[1]!=1 & !(sum(x==0)>sum(x==1)),
x[2]==1 & x[8]==0 | x[2]!=1 & x[8]!=0,
x[3]==1 & x[2]==0 | x[3]!=1 & x[2]!=0 ,
x[4]==1 & (x[3]==0 & x[6]==0) | x[4]!=1 & !(x[3]==0 & x[6]==0),
x[5]==1 & sum(x==0)>=1 | x[5]!=1 & !(sum(x==0)>=1),
x[6]==1 & sum(x==0)>=2 | x[6]==1 & !(sum(x==0)>=
881:2), x[7]==1 & x[5]==0 | x[7]!=1 & x[5]!=0, x[8]==1 & (x[1]==1 & x[6]==1) | x[8]!=1 & !(x[1]==1 & x[6]==1) )) } TE[apply(TE,1,f),]
882:卵の名無しさん
20/01/27 15:49:58 54kq03+D.net
ここまで複雑だと推論での結論はかなり困難だと思う。
AからHの8人はそれぞれ正直者か嘘つきであり、誰が正直者か嘘つきかはお互いに知っている。
A,B,C,D,Eは嘘つきなら必ず嘘をつくが、F,G,Hは嘘つきでも正しいことを言う場合がある。
次の証言から誰を確実に正直者と断定できるか？
A「嘘つきの方が正直者より多い」
B「Hは嘘つきである」
C「Bは嘘つきである」
D「CもFも嘘つきである」
E「8人の中に、少なくとも1人嘘つきがいる」
F「8人の中に、少なくとも2人嘘つきがいる」
G「Eは嘘つきである」
H「AもFも正直者である」
プログラムで解くのはたやすいが正しいプログラムを書くには注意力が必要。 ,や)のミス一つで誤答がでる
TE=expand.grid(0:1,0:1,0:1,0:1,0:1,0:1,0:1,0:1)
colnames(TE)=LETTERS[1:8]
f <- function(x){
all(c(
x[1]==1 & sum(x==0)>sum(x==1) | x[1]!=1 & !(sum(x==0)>sum(x==1)),
x[2]==1 & x[8]==0 | x[2]!=1 & x[8]!=0,
x[3]==1 & x[2]==0 | x[3]!=1 & x[2]!=0 ,
x[4]==1 & (x[3]==0 & x[6]==0) | x[4]!=1 & !(x[3]==0 & x[6]==0),
x[5]==1 & sum(x==0)>=1 | x[5]!=1 & !(sum(x==0)>=1),
x[6]==1 & sum(x==0)>=2 | x[6]==0,
x[7]==1 & x[5]==0 | x[7]==0,
x[8]==1 & (x[1]==1 & x[6]==1) | x[8]==0
))
}
TE[apply(TE,1,f),]
883:卵の名無しさん
20/01/28 08:12:46.55 xazq0gTe.net
では次の問題です
次の動画をみさくら語で解説しましょう
URLﾘﾝｸ(www.youtube.com)
国試浪人の事務員でなければできるはずだねW
884:卵の名無しさん
20/01/28 09:10:54 G+B9uPjF.net
>>834
みさくら語ってド底辺頭脳の言語か？
RとかPythonなら使えるけど。
885:コンブ薬屋
20/01/28 10:56:41.57 xazq0gTe.net
ほほう、どうやら
底辺の国試浪人には解けないようだね
886:卵の名無しさん
20/01/28 10:59:24.13 G+B9uPjF.net
rm(list=ls())
graphics.off()
seg <- function(a,b,...){
segments(Re(a),Im(a),Re(b),Im(b),col='gray',...)}
plot(0:10,type='n',xlim=c(0,10),ylim=c(0,10),asp=1,bty='l')
seg(10,10+10i,lwd=5) ; seg(0,10,lwd=5) ; seg(0,10i,lwd=5) ; seg(10i,10+10i,lwd=5)
cir <- function(a,b=0,r){
x=seq(a-r,a+r,by=0.01)
x=x[0<x & x<10]
y=b+sqrt(r^2-(x-a)^2)
lines(x,y)
}
for(i in 0:10) cir(i,0, 10-i*0.5)
cir2 <- function(a=10,b,r){
y=seq(b-r,b+r,by=0.01)
y=y[0<y & y<10]
x=a-sqrt(r^2-(y-b)^2)
lines(x,y)
}
for(i in 0:10) cir2(10,i,5-i*0.5)
887:卵の名無しさん
20/01/28 11:00:43.78 G+B9uPjF.net
cir3 <- function(a=0,b,r){
y=seq(b-r,b+r,by=0.01)
y=y[0<y & y<10]
x=a+sqrt(r^2-(y-b)^2)
lines(x,y)
}
for(i in 0:10) cir3(0,i,10-i*0.5)
cir4 <- function(a,b=10,r){
x=seq(a-r,a+r,by=0.01)
x=x[0<x & x<10]
y=b-sqrt(r^2-(x-a)^2)
lines(x,y)
}
for(i in 0:10) cir4(i,10,5-i*0.5)
URLﾘﾝｸ(i.imgur.com)
888:コンブ薬屋
20/01/28 12:07:50 xazq0gTe.net
さあ、みさくら語に翻訳しよう
「何？街の薬剤師？大学教授の娘が、街の一薬剤師と結婚したいと言うのか？」
「いけませんの？ 」
「絶対に反対だよ。何代も続いている有名な製薬会社オーナーは別として、一般の薬剤師になるものは、会社に残りたくても残れず、仕方なく街の薬剤師になる場合が多いのだ」
「医学部教授と街の薬剤師。比較するだけでもおかしいよ」
889:卵の名無しさん
20/01/28 12:28:48.28 zoyO6+ty.net
社交性が無く、職場で干されているアスペルガー医者が建てたスレはここですか？
890:卵の名無しさん
20/01/28 13:27:34 G+B9uPjF.net
>>836
みさくら語ってド底辺頭脳の言語は知らんから、リンクも踏んでないよ。
891:卵の名無しさん
20/01/28 13:29:40 G+B9uPjF.net
新型コロナウイルスの統計処理が自分の計算と合わないのが不思議なんだよね。
2019-nCov（２０１９年新型コロナウイルス）の論文
URLﾘﾝｸ(www.thelancet.com)(20)30183-5/fulltext
は生データがついているので自分で統計処理が再現できるので遊べる。
2019-nCov感染者41例中ICU収容13例（うち現喫煙者0),28例非収容(うち現喫煙者3)である。
このことから、煙草を吸うと重症化が防げると結論できるか？
エントリーに０を含むときにはχ二乗検定は無効なのでFisher testで計算していると思うのだが
Fisher's Exact Test for Count Data
d
892:ata: cbind(hit, shot - hit) p-value = 0.539 alternative hypothesis: true odds ratio is not equal to 1 95 percent confidence interval: 0.000000 5.257728 sample estimates: odds ratio となって、論文のp値 0.31とは合致しないなぁ。
893:卵の名無しさん
20/01/28 13:47:56 G+B9uPjF.net
I hate the so-called Fisher exact test: a pointer to one of my favorite articles
URLﾘﾝｸ(statmodeling.stat.columbia.edu)
ベイズで計算してみるか。
# contingency stan for reproducibility of p.value
library(rstan)
options(mc.cores = parallel::detectCores())
rstan_options(auto_write = TRUE)
contingency2stan <- function(mat,alpha=0.05){
N=apply(mat,1,sum)
data <- list(mat=mat,N=N)
stanModel=readRDS('contingency.rds')
fit=sampling(stanModel,data=data, chain=1, iter=10000)
pars=c('y1','n1_y1','y0','n0_y0','d','RR','OR','NNT')
print(fit, pars=pars,digits=4,prob=c(0.025,0.5,0.975))
ms=rstan::extract(fit)
p.sim=NULL ;attach(ms) ; for(i in 1:length(p11)){
MAT=round(matrix(c(y1[i],n1_y1[i],y0[i],n0_y0[i]),ncol=2,byrow=TRUE))
p.sim=append(p.sim,chisq.test(MAT)$p.value)
}; detach(ms)
BEST::plotPost(p.sim,compVal = 0.05, xlab='p.value')
cat('Pr(p.value < ',alpha,') = ', mean(p.sim<alpha),'\n')
cat('Original p.value = ',chisq.test(mat)$p.value,'\n')
}
894:卵の名無しさん
20/01/28 13:57:41.14 G+B9uPjF.net
>>836
ド底辺スレでの宿題はまだできてないのか？
ド底辺シリツ医大卒の裏口バカは、空室なしの10人5部屋割付の数をまだ出せないのか？
部屋の名前を '底辺','私立','ガチ','裏口','バカ'と命名すると
このように
ガチ ガチ ガチ ガチ ガチ ガチ ガチ ガチ ガチ ガチ
ガチ ガチ ガチ ガチ ガチ ガチ ガチ ガチ ガチ バカ
ガチ ガチ ガチ ガチ ガチ ガチ ガチ ガチ ガチ 私立
ガチ ガチ ガチ ガチ ガチ ガチ ガチ ガチ ガチ 底辺
ガチ ガチ ガチ ガチ ガチ ガチ ガチ ガチ ガチ 裏口
ガチ ガチ ガチ ガチ ガチ ガチ ガチ ガチ バカ ガチ
から
裏口 裏口 裏口 裏口 裏口 裏口 裏口 裏口 底辺 裏口
裏口 裏口 裏口 裏口 裏口 裏口 裏口 裏口 裏口 ガチ
裏口 裏口 裏口 裏口 裏口 裏口 裏口 裏口 裏口 バカ
裏口 裏口 裏口 裏口 裏口 裏口 裏口 裏口 裏口 私立
裏口 裏口 裏口 裏口 裏口 裏口 裏口 裏口 裏口 底辺
裏口 裏口 裏口 裏口 裏口 裏口 裏口 裏口 裏口 裏口
まで並べて
裏口 底辺 裏口 バカ ガチ 底辺 底辺 私立 私立 バカ
バカ 裏口 私立 私立 私立 バカ ガチ 裏口 ガチ 底辺
バカ 私立 バカ バカ ガチ 裏口 底辺 バカ 裏口 裏口
裏口 裏口 ガチ 裏口 私立 底辺 裏口 裏口 バカ ガチ
底辺 バカ ガチ バカ ガチ 私立 裏口 バカ バカ 裏口
部屋が５種類あるのを数えるだけだぞ。
ド底辺シリツ医大卒の裏口バカって数もろくに数えられないのか？
895:コンブ薬屋
20/01/28 14:21:54 xazq0gTe.net
ほら、んほおぉぉぉっ、って言ってみな
んほおぉぉぉっ、ってW
896:卵の名無しさん
20/01/28 15:35:57.65 G+B9uPjF.net
最新のLancetに掲載された論文：新型コロナウイルスの治療に少量～中等量のステロイド（メチルプレドニゾロンで４０～１２０ｍｇ/日)使用の有無による死亡例の数字は
Death 4/9 (44.4%) 2/32 (6.3%)
なぜか、他の項目にあったｐ値が掲載されていない。
Fisher's Exact Test for Count Data
data: mat
p-value = 0.01481456
alternative hypothesis: true odds ratio is not equal to 1
95 percent confidence interval:
1.217274433 151.944072111
sample estimates:
odds ratio
10.93370175
なんでも確率変数にするのがベイズ統計。
ここで問題
ｐ値の９５％信頼区間とp値が0.05未満になる確率を求めよ。
897:卵の名無しさん
20/01/28 15:37:06.96 G+B9uPjF.net
0.05未満になる確率は
これを持って統計的に有意だというどの程度
898:の信頼度を持っていえるか、ってことなんだな。
899:卵の名無しさん
20/01/28 15:48:52.17 G+B9uPjF.net
## 2x2 tableからsensitivity,specificity,PPV,NPVの信頼区間を計算
cont2ci<- function(TP,FP,FN,TN){ # contingency table 2 confidence interval
library(BEST)
sensitivity=TP/(TP+FN)
specificity=TN/(TN+FP)
pLR=TP/FP # sensitivity/(1-specificity)
nLR=FN/TN # (1-TP)/(1-FP) = (1-sensitivity)/specificity
PPV=TP/(TP+FP)
NPV=TN/(TN+FN)
par(mfrow=c(3,2))
col=sample(colours(),1)
plotPost(sensitivity,showMode = F,col=col)
plotPost(specificity,showMode = F,col=col)
plotPost(PPV,showMode = F,col=col)
plotPost(NPV,showMode = F,col=col)
plotPost(pLR,showMode = F,col=col)
plotPost(nLR,showMode = F,col=col)
sn=hdi(sensitivity)[1:2]
sp=hdi(specificity)[1:2]
ppv=hdi(PPV)[1:2]
npv=hdi(NPV)[1:2]
plr=hdi(pLR)[1:2]
nlr=hdi(nLR)[1:2]
par(mfrow=c(1,1))
rbind(sn,sp,ppv,npv,plr,nlr)
}
900:卵の名無しさん
20/01/28 15:50:33.03 G+B9uPjF.net
# ポアソンモデル
k=1e5
TP=rpois(k,43)
FP=rpois(k,2)
FN=rpois(k,35)
TN=rpois(k,25)
cont2ci(TP,FP,FN,TN)
# 二項分布モデル(周辺度数固定)
k=1e5
TP=rbinom(k,43+35,43/(43+35))
FP=rbinom(k, 2+25, 2/( 2+25))
FN=rbinom(k,35+43,35/(35+43))
TN=rbinom(k,25+ 2,25/(25+ 2))
cont2ci(TP,FP,FN,TN)
# 二項分布モデル(総数固定)
k=1e5
mat=matrix(c(43,2,35,25),ncol=2)
colnames(mat)=c("逆転陽性","逆転陰性")
rownames(mat)=c("除菌成功","除菌失敗")
mat
N=apply(mat,1,sum) # 78 vs 27
N.total=sum(mat)
Succ=rbinom(k,N.total,78/N.total)
Fail=rbinom(k,N.total,27/N.total)
FP=rbinom(k,Fail, 2/Fail)
FN=rbinom(k,Succ,35/Succ)
TN=Fail - FP
TP=Succ - FN
cont2ci(TP,FP,FN,TN)
901:コンブ薬屋
20/01/28 17:55:30 xazq0gTe.net
では次の問題です
この病棟における入院期間の平均値と治癒率を求めよ
URLﾘﾝｸ(togetter.com)
902:卵の名無しさん
20/01/28 19:19:47.40 vtD6D8ed.net
０がエントリーにあると扱いが難しいな。
rm(list=ls())
mat.smoke=cbind(c(0,3),c(13,28)-c(0,3))
mat.smoke # 喫煙のICU入院のクロス表
Epi::twoby2(mat.smoke)
fmsb::oddsratio(mat.smoke)
library(BayesFactor)
library(BEST)
k=1e3
bfp=contingencyTableBF(mat.smoke,'poisson',post=T,iter=k)
matp=as.data.frame(bfp) # ポアソン分布で乱数化
apply(matp,2,mean)
bfi=contingencyTableBF(mat.smoke,'indepMulti','rows',post=T,iter=k)
bfi=as.data.frame(bfi) # 行の周辺度数固定の二項分布で乱数化
mati=bfi[,1:4]*sum(mat.smoke)
apply(mati,2,mean)
bfj=contingencyTableBF(mat.smoke,'jointMulti',post=T,iter=k)
matj=as.matrix(bfj)*sum(mat.smoke) # 総数数固定の二項分布で乱数化
apply(matj,2,mean)
903:卵の名無しさん
20/01/28 19:19:58.78 vtD6D8ed.net
mat2pf <- function(x){
fisher.test(matrix(round(x),nrow=2))$p.value # roundで整数化してFisherテスト
}
mat2pc <- function(x){
chisq.test(matrix(x,nrow=2))$p.value # χ二乗検定
}
ppf=apply(matp,1,mat2pf)
plotPost(ppf,compVal=0.05)
ppc=apply(matp,1,mat2pc)
plotPost(ppc,compVal=0.05)
pif=apply(mati,1,mat2pf)
plotPost(pif,compVal=0.05)
pic=apply(mati,1,mat2pc)
plotPost(pic,compVal=0.05)
pjf=apply(matj,1,mat2pf)
plotPost(pjf,compVal=0.05)
pjc=apply(matj,1,mat2pc)
plotPost(pjc,compVal=0.05)
source('tools.R')
jags4prop(0,3,13,28)
904:卵の名無しさん
20/01/28 19:22:05.70 vtD6D8ed.net
リンクは踏まないから、こういう風に問題を書けよ。
"某女子大には決して嘘をつかない学生と必ず嘘をつく学生がいることがわかっている。
この女子大の学生（嘘つきかどうかは不明）から
「あなたのいうことが正しければ手コキかフェラをしてあげる、間違っていれば何もしてあげない」と言われた。
この女子大生にフェラをしてもらうには何と言えばいいか？
905:卵の名無しさん
20/01/28 19:47:42.30 xazq0gTe.net
ここは下品ないんたねっとですね
まさにこのスレはド底辺と呼ぶにふさわしいですねW
906:卵の名無しさん
20/01/28 20:34:26 vtD6D8ed.net
これで、
# 2x2のクロス表からstanを使ってMCMCサンプリングしてFisherテストしてp.valueの配列からその分布を表示させる
library(rstan)
options(mc.cores = parallel::detectCores())
rstan_options(auto_write = TRUE)
contingency2stan2pf <- function(mat,alpha=0.05,print=T){
N=apply(mat,1,sum)
data <- list(mat=mat,N=N)
stanModel=readRDS('contingency.rds')
fit=sampling(stanModel,data=data, chain=1, iter=25000,warmup=5000)
pars=c('y1','n1_y1','y0','n0_y0','d','RR','OR','NNT')
if(print) print(fit, pars=pars,digits=4,prob=c(0.025,0.5,0.975))
ms=rstan::extract(fit)
p.sim=NULL
for(i in 1:length(p11)){
MAT=round(with(ms,matrix(c(y1[i],n1_y1[i],y0[i],n0_y0[i]),ncol=2,byrow=TRUE)))
p.sim=append(p.sim,fisher.test(MAT)$p.value)
}
BEST::plotPost(p.sim,compVal = 0.05, xlab='p.value')
cat('Pr(p.value < ',alpha,') = ', mean(p.sim<alpha),'\n')
cat('Original p.value = ',fisher.test(mat)$p.value,'\n')
print(summary(p.sim))
invisible(p.sim)
}
サンプル数を増やしてｐを小さくして信用させるという手法を見破れるかな。
907:卵の名無しさん
20/01/28 20:35:30 vtD6D8ed.net
>>854
ちゃんと答えないとフェラしてもらえないよ。
908:卵の名無しさん
20/01/28 20:38:30 vtD6D8ed.net
> (mat.steroid=cbind(c(4,2),c(9,32)-c(4,2))) # ステロイドと死亡数
[,1] [,2]
[1,] 4 5
[2,] 2 30
> contingency2stan2pf(mat.steroid,p=F)
Pr(p.value < 0.05 ) = 0.7024
Original p.value = 0.01481455782
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.0000121509 0.0055484536 0.0148145578 0.0481260381 0.0611450791 1.0000000000
909:卵の名無しさん
20/01/28 21:01:49.02 vtD6D8ed.net
>>855
フィッシャー検定は非負整数しか許さないから、こんな感じで
URLﾘﾝｸ(i.imgur.com)
グラフにすると途切れができてしまうのが美しくない。
χ二乗検定版を作ってみよう。
910:卵の名無しさん
20/01/28 21:03:52.80 IbKEPttD.net
国試浪人の素人童貞のスレはここですか？
911:卵の名無しさん
20/01/28 21:07:29.78 vtD6D8ed.net
>>858
# χ二乗テスト版
contingency2stan2pc <- function(mat,alpha=0.05,print=T){
N=apply(mat,1,sum)
data <- list(mat=mat,N=N)
stanModel=readRDS('contingency.rds')
fit=sampling(stanModel,data=data, chain=1, iter=25000,warmup=5000)
pars=c('y1','n1_y1','y0','n0_y0','d','RR','OR','NNT')
if(print) print(fit, pars=pars,digits=4,prob=c(0.025,0.5,0.975))
ms=rstan::extract(fit)
p.sim=NULL
for(i in 1:length(p11)){
MAT=with(ms,matrix(c(y1[i],n1_y1[i],y0[i],n0_y0[i]),ncol=2,byrow=TRUE))
p.sim=append(p.sim,chisq.test(MAT)$p.value)
}
BEST::plotPost(p.sim,compVal = 0.05, xlab='p.value')
cat('Pr(p.value < ',alpha,') = ', mean(p.sim<alpha),'\n')
cat('Original p.value = ',chisq.test(mat)$p.value,'\n')
print(summary(p.sim))
invisible(p.sim)
}
(mat.smoke=cbind(c(0,3),c(13,28)-c(0,3))) # 喫煙とICU入院のクロス表
contingency2stan2pc(mat.smoke,p=F)
(mat.steroid=cbind(c(4,2),c(9,32)-c(4,2))) # ステロイドと死亡数
p.sim=contingency2stan2pc(mat.steroid,p=F)
sapply(c(0.0001,0.001,0.01,0.05),function(x) mean(p.sim<x))
lines(density(p.sim))
この方が美しいな。
URLﾘﾝｸ(i.imgur.com)
912:卵の名無しさん
20/01/28 21:08:19.29 vtD6D8ed.net
>>859
これできないとフェラしてもらえないよ。
"某女子大には決して嘘をつかない学生と必ず嘘をつく学生がいることがわかっている。
この女子大の学生（嘘つきかどうかは不明）から
「あなたのいうことが正しければ手コキかフェラをしてあげる、間違っていれば何もしてあげない」と言われた。
この女子大生にフェラをしてもらうには何と言えばいいか？
913:卵の名無しさん
20/01/28 21:10:44.10 vtD6D8ed.net
人類の三大発明は
言語・貨幣・フェラチオである。
ｂｙ 吉田戦車
914:卵の名無しさん
20/01/28 21:52:45.89 cUmejCEO.net
お兄ちゃん、イサキは、イサキは獲れたのっ？
915:卵の名無しさん
20/01/28 21:57:08.65 vtD6D8ed.net
比率の比較でなくて平均値の比較のときはｐ値の信頼区間算出はブートストラップを使えば簡単にだせるな。
知能指数が平均101と99の各々500人からなる集団の比較
# t検定(等分散)
T.test=function(n1,n2,m1,m2,sd1,sd2){
SE12=sqrt((1/n1+1/n2)*((n1-1)*sd1^2+(n2-1)*sd2^2)/((n1-1)+(n2-1)))
T=(m1-m2)/SE12
p.value=2*pt(abs(T),n1-1+n2-1,lower.tail = FALSE)
return(c(p.value=p.value,t.statistic=T))
}
n=500
a=99
b=101
A=a+scale(rnorm(n))*15
B=b+scale(rnorm(n))*15
T.test(n,n,a,b,15,15)
t.test(A,B,var=T)$p.value
sim <- function(x=n) t.test(sample(A,x,rep=T),sample(B,x,rep=T))$p.value
k=1e4
pv=replicate(k,sim())
summary(pv)
quantile(pv,prob=c(0.005,0.025,0.5,0.975,0.995))
BEST::plotPost(pv,compVal = 0.05,showMode = T,breaks=50)
916:卵の名無しさん
20/01/28 22:00:00.87 vtD6D8ed.net
>>863
釣り人（嘘つきかどうかは不明）から
「あなたのいうことが正しければイサキかイサキの餌のどちらかをあげる、間違っていれば何もあげない」と言われた。
この釣り人からイサキをもらうには何と言えばいいか？
917:卵の名無しさん
20/01/29 05:29:27 qmYv2xxr.net
少し一般化、復元抽出するか否
918:かや個数が不均等でも計算できるように改変。 # 母集団Na,Nb個 平均Ma,Mb 標準偏差SDa,SDbからna,nbを抽出しt検定比較ときのp値の分布 pCI <- function(Na=500,Nb=Na,Ma=99,Mb=101,SDa=15,SDb=SDa,na=50,nb=na, rep=FALSE,k=1e3,print=TRUE){ A=Ma+scale(rnorm(Na))*SDa B=Mb+scale(rnorm(Nb))*SDb sim <- function(m) t.test(sample(A,na,replace=rep),sample(B,nb,replace=rep))$p.value pv=replicate(k,sim()) if(print){ print(T.test(n,n,a,b,15,15)) print(summary(pv)) BEST::plotPost(pv,compVal = 0.05,showMode = T,breaks=50)} invisible(pv) }
919:卵の名無しさん
20/01/29 05:46:52.67 qmYv2xxr.net
デフォルト設定では母集団から1/10の標本を抽出するように変更
# 母集団Na,Nb個 平均Ma,Mb 標準偏差SDa,SDbからna,nbを抽出しt検定比較ときのp値の分布
pCI <- function(Na=500,Nb=Na,Ma=99,Mb=101,SDa=15,SDb=SDa,na=Na/10,nb=Nb/10,
rep=FALSE,k=1e3,print=TRUE){
A=Ma+scale(rnorm(Na))*SDa
B=Mb+scale(rnorm(Nb))*SDb
sim <- function(m) t.test(sample(A,na,replace=rep),sample(B,nb,replace=rep))$p.value
pv=replicate(k,sim())
if(print){
print(t.test(A,B))
print(summary(pv))
BEST::plotPost(pv,compVal = 0.05,showMode = T,breaks=50)}
invisible(pv)
}
920:卵の名無しさん
20/01/29 05:53:26.38 qmYv2xxr.net
標準大学の新入生の知能指数の平均を100とする。
ド底辺シリツ医大の新入生の知能指数の平均が85であった。
各々から１/10を無作為抽出して知能指数をｔ検定したときのｐ値の期待値、中央値を求めよ。
また、ｐ値が0.05以上になってド底辺シリツ医の新入生の知能指数は統計的に有意差はないと主張できる確率はいくらか？
# 母集団Na,Nb個 平均Ma,Mb 標準偏差SDa,SDbからna,nbを抽出しt検定比較ときのp値の分布
pCI <- function(Na=500,Nb=Na,Ma=99,Mb=101,SDa=15,SDb=SDa,na=Na/10,nb=Nb/10,
rep=FALSE,k=1e3,print=TRUE){ # IQ:m=100,sd=15
A=Ma+scale(rnorm(Na))*SDa
B=Mb+scale(rnorm(Nb))*SDb
sim <- function(m) t.test(sample(A,na,replace=rep),sample(B,nb,replace=rep))$p.value
pv=replicate(k,sim())
if(print){
print(t.test(A,B))
print(summary(pv))
BEST::plotPost(pv,compVal = 0.05,showMode = T,breaks=50)}
invisible(pv)
}
pCI(Na=100,Ma=100,Mb=85)
921:卵の名無しさん
20/01/29 05:56:48.55 qmYv2xxr.net
>>868
この手法で新薬の血中濃度は老人でも統計的に有意差はないので、安心して老人にも処方できます
という商用パンプを目にする。
底辺シリツ医＝裏口容疑者はこういうのに騙されるんだよなぁ。
922:卵の名無しさん
20/01/29 07:45:56.70 qmYv2xxr.net
>>868
平均値100 標準偏差15で定義される知能指数で
標準大学の新入生の知能指数の平均が100
裏口シリツ医大の新入生の知能指数の平均が85であったとする。
1学年は100人として
各大学から1/10を無作為抽出して知能指数をｔ検定したときのｐ値の期待値、中央値を求めよ。
また、ｐ値が0.05以上になって裏口シリツ医大の新入生の知能指数は統計的に有意差はないと主張できる確率はいくらか？
923:卵の名無しさん
20/01/29 10:58:46.40 wNUvDIoi.net
>>863
大漁っ！！イサキぃぃ！！おにいちゃんかっこいいいいぃぃぃい ぃくううううう！
924:卵の名無しさん
20/01/29 12:49:44 NU71nZE3.net
私立出身医が一番多いスレであろう「開業医が集うスレ」のほうには書き込みせんのな？
925:卵の名無しさん
20/01/29 15:37:18.14 qmYv2xxr.net
>>872
書いているよ。
開業医スレでもシリツは馬鹿にされているじゃん。
ｽﾚﾘﾝｸ(hosp板:54番)
>>23
生データもある程度掲載されているから自分で統計処理が再現でいるな。
こういうのを考えてみると面白い。
Fisher testの結果が俺の計算と異なって気持ちが悪いままだが。
2019-nCov感染者41例中ICU収容13例（うち現喫煙者0),28例非収容(うち現喫煙者3)である。
このことから、煙草を吸うと重症化が防げると結論できるか？
926:卵の名無しさん
20/01/29 15:37:56.49 qmYv2xxr.net
>>871
そんなにレスじゃ、フェラもイサキも手に入れることはできんぞ。
927:卵の名無しさん
20/01/29 17:41:57.38 wNUvDIoi.net
では次の問題です
タバコは脳の病気です
この事実を証明せよ
URLﾘﾝｸ(togetter.com)
928:卵の名無しさん
20/01/29 17:46:37.77 qmYv2xxr.net
>853
929:卵の名無しさん
20/01/29 19:50:54 qmYv2xxr.net
>>872
俺が書かなくてもこういう投稿を目にする。
まあ、シリツ医という表記は俺の影響かもしれん。
自治医大とかを除外する意味で私立と書いていないわけだが。
>>
優秀な女子はぜひ国立医に入ってほしい
シリツ医は開業医のアホ子弟専用の専門学校
930:卵の名無しさん
20/01/29 21:34:38.87 wNUvDIoi.net
まあたしかに私立は学費が高いのは認める
わしも経済的理由で
合格した慶応は学費が高いので蹴って
地元イナカの国立にいった
931:卵の名無しさん
20/01/30 05:17:39 m/EdOA
932:9B.net
933:卵の名無しさん
20/01/30 08:20:09.67 b9qEKmQr.net
>>879
イナカモノにしてみれば
地元以外のところに大学が存在しているという時点で
学費以外にもアパート代や食費などの参入障壁があるんだな
934:卵の名無しさん
20/01/30 08:31:17 17lSJR5q.net
>>880
シリツ医大専門受験予備校のバイトは一コマ90分で手取り2万円だったな。大学の授業料が年間14万4千円の頃の話。
地方にはそんな割のいいバイトはない。
935:卵の名無しさん
20/01/30 20:25:51.32 m/EdOA9B.net
Rで困ったらstackoveflow.comだな。
library(Rmpfr)
vect = rep(0, 5)
for(i in 1:5){
vect[i] = mpfr(x=10^-(10*i), precBits=100)
}
# My vector has just turned to a list
vect
# Sum of list is an error
sum(vect)
# Turn it to a vector
converted_vect = Rmpfr::mpfr2array(vect, dim = length(vect))
converted_vect
# Now my sum and prod work fine and the precision is not lost
sum(converted_vect)
prod(converted_vect)
936:卵の名無しさん
20/01/30 21:10:39 m/EdOA9B.net
パッケージの内部関数を使って解決してくれた。
The function mpfr2array is not suposed to be called by the user, it is an internal tool for the package. However it's one way to solve the problem.
937:卵の名無しさん
20/01/30 21:21:54.99 m/EdOA9B.net
Rはデフォルトでは不定長さ整数を扱えないので工夫がいる。
その点はpythonやHaskellの方が有利。
"
H(n) = Σ[k=1,2,...,n] 1/k
とする。H(n)を既約分数で表したときの分子の整数をf(n)と表す。
(1）lim[n→∞] H(n) を求めよ、答えのみで良い。
(2）n=1,2,...に対して、f(n)に現れる1桁の整数を全て求めよ
"
rm(list=ls())
H <- function(n,prec=1000){ # Σ 1/kを既約分数表示する
GCD <- function(a,b){ # ユークリッドの互除法
r = a%%b # a=bq+r ⇒ a%%b=b%%rで最大公約数表示
while(r!=0){a = b ; b = r ; r = a%%b}
b }
library(Rmpfr)
one = mpfr(1, prec) # 1(one)を1000桁精度に設定
nn = 1:n # nn : 1 2 3 ... n
nume=numeric(n) # 分子の容器
for(i in nn) nume[i] = prod(nn[-i])*one # nnからi番目の要素を除いて乗算し精度アップ
nume = mpfr2array(nume, dim = length(nume)) # mpfr2arrayで加算を可能にする
Nume = sum(nume) # numeの総和を計算して分子に
Deno = factorial(n)*one # 分母のn!を精度アップ
gcd = GCD(Nume,Deno) # Numerator/Denominator約分するため最大公約数を計算
res=list(nume = Nume/gcd,deno=Deno/gcd,ratio=as.numeric(Nume/Deno)) # 最大公約数で除算して
cat(as.numeric(res$nume),'/',as.numeric(res$deno),'\n') # 数値化して分数表示
invisible(res)
}
938:卵の名無しさん
20/01/30 23:51:43 m/EdOA9B.net
#Reed-Frost Model
ReedFrost=function(
p=0.04,
N=100,
T=40)
{
q=1-p
I=numeric(T)
S=numeric(T)
I[1]=1
S[1]=N-I[1]
for(t in 1:(T-1)){
I[t+1]=S[t]*(1-q^I[t])
S[t+1]=S[t]-I[t+1]
}
plot(1:T,I,type="l",lwd=2, ylim=c(0,N
939:),xlab="time",ylab="persons",main=paste("Reed-Frost Model p= ",p)) lines(S,lty=2,col=2,lwd=2) lines(N-S,lty=3,col=3,lwd=2) legend("topright",bty="n",legend=c("Infected","Susceptible","Immunized"),lty=1:3,col=1:3,lwd=2) } par(mfrow=c(1,2)) ReedFrost(0.04) ReedFrost(0.015) ReedFrost(0.30) ReedFrost(0.03) ReedFrost(0.003)
940:卵の名無しさん
20/01/30 23:53:06 m/EdOA9B.net
こっちは、モンテカルロによるシミュレーション
# Reed-Frost and Greenwood epidemic models
# written by Dennis Chao (1/2009)
# reedfrost - the Reed-Frost epidemic model
# p = probability of transmission
# I0 = initial number of infecteds
# S0 = initial number of susceptibles
# n = number of trials
# greenwood = set to TRUE for the Greenwood model, otherwise run Reed-Frost
# outputs the number of infected and susceptibles over time (as I and S)
reedfrost <- function(p, I0, S0, n, greenwood=FALSE) {
S <- St <- rep(S0, n)
I <- It <- rep(I0, n)
q <- 1-p
time <- 0
while (sum(It)>0) {
if (greenwood)
It <- rbinom(n, St, ifelse(It>0, p, 0))
else
It <- rbinom(n, St, 1-(q^It))
St <- St-It
I <- rbind(I, It)
S <- rbind(S, St)
time <- time+1
}
I <- as.matrix(I)
S <- as.matrix(S)
list(I0=I0,S0=S0,p=p,n=n,I=I,S=S)
}
941:卵の名無しさん
20/01/31 00:59:52 24QiZJ0Y.net
f <- function(n,prec=1000){ # Σ 1/kを既約分数表示する
if(n==1){
cat(1,'\n')
invisible(1)
}else{
GCD <- function(a,b){ # ユークリッドの互除法
r = a%%b # a=bq+r ⇒ a%%b=b%%rで最大公約数表示
while(r!=0){a = b ; b = r ; r = a%%b}
b }
library(Rmpfr)
one = mpfr(1, prec) # 1(one)を1000桁精度に設定
nn = 1:n # nn : 1 2 3 ... n
nume=numeric(n) # 分子の容器
for(i in nn) nume[i] = prod(nn[-i])*one # nnからi番目の要素を除いて乗算し精度アップ
nume = mpfr2array(nume, dim = length(nume)) # mpfr2arrayで加算を可能にする
Nume = sum(nume) # numeの総和を計算して分子に
Deno=factorialZ(n) # 分母 n! = factorial(n*one)
gcd = GCD(Nume,Deno) # Numerator/Denominator約分するため最大公約数を計算
res=list(nume = Nume/gcd,deno=Deno/gcd,ratio=as.numeric(Nume/Deno)) # 最大公約数で除算して
# 分数表示 give.head=FALSEでheader除去,digits.dで桁数を指定
# capture.outputで変数に取り込み
nm = capture.output(str(res$nume, give.head=FALSE,digits.d = prec))
dn = capture.output(str(res$deno, give.head=FALSE,digits.d = prec))
cat(paste0(nm,'/',dn,'\n'))
invisible(res)
}}
for(i in 1:30) f(i)
942:卵の名無しさん
20/01/31 01:00:30 24QiZJ0Y.net
> for(i in 1:30) f(i)
1
3 / 2
11 / 6
25 / 12
137 / 60
49 / 20
363 / 140
761 / 280
7129 / 2520
7381 / 2520
83711 / 27720
86021 / 27720
1145993 / 360360
1171733 / 360360
1195757 / 360360
2436559 / 720720
42142223 / 12252240
14274301 / 4084080
275295799 / 77597520
55835135 / 15519504
18858053 / 5173168
19093197 / 5173168
122755644038509457 / 32872539188238750
186187999757029099 / 49308808782358125
14112026408124257248 / 3698160658676859375
185305423634953775872 / 48076088562799171875
5051322526706550956032 / 1298054391195577640625
11894590428248250515456 / 3028793579456347828125
1043915747995966839455744 / 263505041412702261046875
2255784105806550548873216 / 564653660170076273671875
943:卵の名無しさん
20/01/31 01:12:13.74 24QiZJ0Y.net
>>884
H(n) = Σ[k=1,2,...,n] 1/k
とする。H(n)を既約分数で表したときの分子の整数をf(n)と表す。
f(77)を求めよ。
> f(77)
17610982920730618962802441030965952272844514966520010106103127939813509744122599441432576
/ 35740194818708235
944:59764745233429885438685864430565417716720215849457565210956573486328125
945:卵の名無しさん
20/01/31 07:13:25 24QiZJ0Y.net
>>884
n=100と大きいとNAが混ざるな
946:卵の名無しさん
20/01/31 07:15:08 24QiZJ0Y.net
ｆ ＜－ ｆｕｎｃｔｉｏｎ（ｎ，ｐｒｅｃ＝１００００）｛ ＃ Σ １／ｋを既約分数表示する
ｉｆ（ｎ＝＝１）｛
ｃａｔ（? ?，１，?\ｎ?）
ｉｎｖｉｓｉｂｌｅ（１）
｝ｅｌｓｅ｛
ＧＣＤ ＜－ ｆｕｎｃｔｉｏｎ（ａ，ｂ）｛ ＃ ユークリッドの互除法
ｒ ＝ ａ％％ｂ ＃ ａ＝ｂｑ＋ｒ ⇒ ａ％％ｂ＝ｂ％％ｒで最大公約数表示
ｗｈｉｌｅ（ｒ！＝０）｛ａ ＝ ｂ ； ｂ ＝ ｒ ； ｒ ＝ ａ％％ｂ｝
ｂ ｝
947:卵の名無しさん
20/01/31 07:15:14 24QiZJ0Y.net
ｌｉｂｒａｒｙ（Ｒｍｐｆｒ）
ｏｎｅ ＝ ｍｐｆｒ（１， ｐｒｅｃ） ＃ １（ｏｎｅ）を１００００桁精度に設定
ｎｎ ＝ １：ｎ ＃ ｎｎ ： １ ２ ３ ．．． ｎ
ｎｕｍｅ＝ｎｕｍｅｒｉｃ（ｎ） ＃ 分子の容器
ｆｏｒ（ｉ ｉｎ ｎｎ） ｎｕｍｅ［ｉ］ ＝ ｐｒｏｄ（ｎｎ［－ｉ］）＊ｏｎｅ ＃ ｎｎからｉ番目の要素を除いて乗算し精度アップ
ｎｕｍｅ ＝ ｍｐｆｒ２ａｒｒａｙ（ｎｕｍｅ， ｄｉｍ ＝ ｌｅｎｇｔｈ（ｎｕｍｅ）） ＃ ｍｐｆｒ２ａｒｒａｙで加算を可能にする
Ｎｕｍｅ ＝ ｓｕｍ（ｎｕｍｅ） ＃ ｎｕｍｅの総和を計算して分子に
Ｄｅｎｏ＝ｆａｃｔｏｒｉａｌＺ（ｎ） ＃ 分母 ｎ！ ＝ ｆａｃｔｏｒｉａｌ（ｎ＊ｏｎｅ）
ｇｃｄ ＝ ＧＣＤ（Ｎｕｍｅ，Ｄｅｎｏ） ＃ Ｎｕｍｅｒａｔｏｒ／Ｄｅｎｏｍｉｎａｔｏｒ約分するため最大公約数を計算
ｒｅｓ＝ｌｉｓｔ（ｎｕｍｅ ＝ Ｎｕｍｅ／ｇｃｄ，ｄｅｎｏ＝Ｄｅｎｏ／ｇｃｄ，ｒａｔｉｏ＝ａｓ．ｎｕｍｅｒｉｃ（Ｎｕｍｅ／Ｄｅｎｏ）） ＃ 最大公約数で除算して
＃ 分数表示 ｇｉｖｅ．ｈｅａｄ＝ＦＡＬＳＥでｈｅａｄｅｒ除去，ｄｉｇｉｔｓ．ｄで桁数を指定
＃ ｃａｐｔｕｒｅ．ｏｕｔｐｕｔで切り取りｓｕｂｓｔｒで［１］を除去
＃ ｎｍ ＝ ｃａｐｔｕｒｅ．ｏｕｔｐｕｔ（ｓｔｒ（ｒｅｓ＄ｎｕｍｅ， ｇｉｖｅ．ｈｅａｄ＝ＦＡＬＳＥ，ｄｉｇｉｔｓ．ｄ ＝ ｐｒｅｃ）） ＮＡが混ざるバグあり
ｎｍ ＝ｓｕｂｓｔｒ（ｃａｐｔｕｒｅ．ｏｕｔｐｕｔ（ｒｅｓ＄ｎｕｍｅ）［２］，５，ｎｃｈａｒ（ｒｅｓ＄ｎｕｍｅ）） ＃ ?［１］ １２３４．．? 文字列の５文字目から最後まで
ｄｎ ＝ｓｕｂｓｔｒ（ｃａｐｔｕｒｅ．ｏｕｔｐｕｔ（ｒｅｓ＄ｄｅｎｏ）［２］，５，ｎｃｈａｒ（ｒｅｓ＄ｄｅｎｏ））
ｃａｔ（ｐａｓｔｅ０（ｎｍ，?／?，ｄｎ，?\ｎ?））
ｉｎｖｉｓｉｂｌｅ（ｒｅｓ）
｝｝
948:卵の名無しさん
20/01/31 07:15:42 24QiZJ0Y.net
> f(100)
3055216077446868329553816926933899676639525195878807877583434152044192757431459126874725081455196840519615954410565802448075620352
/ 588971222367687651371627846346807888288472382883312574253249804256440585603406374176100610302040933304083276457607746124267578125
949:卵の名無しさん
20/01/31 07:45:49.87 24QiZJ0Y.net
Reed-Frost モデル
(1) 集団内の感染者と感受性のあるものとの接触はランダムに起こる
(2) 感染者と感受性のあるものが接触して伝播する確率は一定である
(3) 感染のあと必ず免疫が起こる（再感染はしない）
(4) その集団は他の集団から隔離されている
(5) 上記の条件は各時間経過中一定である
ReedFrost=function(
p=0.04, # 1期間内での伝播確率
N=100, # 集団の人数
T=40) # 全期間
{
q=1-p # 伝播させない確率
I=numeric(T) # 感染者数 Infecteds
S=numeric(T) # 感受性のある人数 Susceptibles
I[1]=1 #一人から始まったとする
S[1]=N-I[1]
for(t in 1:(T-1)){
I[t+1]=S[t]*(1-q^I[t]) # Reed-Frost Model
S[t+1]=S[t]-I[t+1]
}
return(list(I,T))
}
950:卵の名無しさん
20/01/31 11:25:54.82 24QiZJ0Y.net
# simulation model using binominal random number
rm(list=ls())
reedfrost <- function(p, I0, S0, n, greenwood=FALSE) {
S <- St <- rep(S0, n) #
951: St : Suscepibles @ time t, S: I <- It <- rep(I0, n) # It : Infected @ time t q <- 1-p # probability of non-transmission time <- 0 while (sum(It)>0) { # until no new transmission if (greenwood) It <- rbinom(n, St, ifelse(It>0, p, 0)) else It <- rbinom(n, St, 1-(q^It)) # how many ppl newly trannsmitted among St St <- St-It # how many ppl left non-infected I <- rbind(I, It) S <- rbind(S, St) time <- time+1 } Inames=NULL for(i in 0:(nrow(I)-1)) Inames[i+1]=paste0('I',i) rownames(I)=Inames Snames=NULL for(i in 0:(nrow(I)-1)) Snames[i+1]=paste0('S',i) rownames(S)=Snames list(I0=I0,S0=S0,p=p,n=n,I=I,S=S) } re=reedfrost(0.05,1,99,100) nr=nrow(re$I) plot(0:(nr-1),re$I[,1],bty='l',type='l',lwd=2, ylim=c(0,max(re$I)),xlab='time',ylab='Infecteds') for(i in 2:n) lines(0:(nr-1),re$I[,i],lwd=2,col=sample(colours()))
952:卵の名無しさん
20/01/31 15:00:17 24QiZJ0Y.net
# SEIR MODEL
"
dS(t)/dt=-bS(t)I(t),
dE(t)/dt=bS(t)I(t)-aE(t) ,
dI(t)/dt=aE(t)-gI(t) ,
dR(t)/dt=gI(t)
a:発症率,b:感染率,g:回復率
"
remove (list = objects() )
graphics.off()
SEIR <- function(
# Parameters
contact_rate = 10, # number of contacts per day
transmission_probability = 0.07, # transmission probability
infectious_period = 5, # infectious period
latent_period = 2, # latent period
# Initial values for sub-populations.
s = 9990, # susceptible hosts
e = 9, # exposed hosts
i = 1, # infectious hosts
r = 0, # recovered hosts
# Output timepoints.
timepoints = seq (0, 100, by=0.5)
){
953:卵の名無しさん
20/01/31 15:00:39 24QiZJ0Y.net
library(deSolve)
# Function to compute derivatives of the differential equations.
seir_model = function (current_timepoint, state_values, parameters)
{
# create state variables (local variables)
S = state_values [1] # susceptibles
E = state_values [2] # exposed
I = state_values [3] # infectious
R = state_values [4] # recovered
with (
as.list (parameters), # variable names within parameters can be used
{
# compute derivatives
dS = (-beta * S * I)
dE = (beta * S * I) - (delta * E)
dI = (delta * E) - (gamma * I)
dR = (gamma * I)
# combine results
results = c (dS, dE, dI, dR)
list (results)
}
)
}
# Compute values of beta (tranmission rate) and gamma (recovery rate).
beta_value = contact_rate * transmission_probability
gamma_value = 1 / infectious_period
delta_value = 1 / latent_period
954:卵の名無しさん
20/01/31 15:00:45 24QiZJ0Y.net
# Compute Ro - Reproductive number.
Ro = beta_value / gamma_value
# Disease dynamics parameters.
parameter_list = c (beta = beta_value, gamma = gamma_value, delta = delta_value)
# Compute total population.
N = s + i + r + e
# Initial state values for the differential equations.
initial_values = c (S = s/N, E = e/N, I = i/N, R = r/N)
# Simulate the SEIR epidemic.
# ?lsoda # Solver for Ordinary Differential Equations (ODE), Switching Automatically Between Stiff and Non-stiff Methods
output = lsoda (initial_values, timepoints, seir_model, parameter_list)
# Plot dynamics of Susceptibles sub-population.
plot (S ~ time, data = output, pch='S', type='b', col = 'blue', bty='l', cex=0.75,
ylab='SEIR')
points(E ~ time, data = output, pch='E', type='b', col = 'yellow', cex=0.75)
points(I ~ time, data = output, pch='I', type='b', col = 'red', cex=0.75)
points(R ~ time, data = output, pch='R', type='b', col = 'green', cex=0.75)
}
SEIR()
args(SEIR)
SEIR(contact_rate=0.5,transmission_probability=0.05,infectious_period=365,latent_period=7,s=99,e=0,i=1,r=0,timepoints = 0:365*6)
955:卵の名無しさん
20/01/31 17:42:44 24QiZJ0Y.net
エンデミックな定常状態を(S?, I?)とおけば、S?N=1R0,I?N=μμ+γ(1?1R0)(15)である。すなわちエンデミックな状態における感受性人口比率と基本再生産数は逆数関係にあり、有病率(prevalence)I?/Nは1?1/R0に比例していて、その比例係数は、感染状態における平均滞在時間1/(μ+γ)とホストの寿命1/μの比である。これらの式はエンデミックな感染症におけるR0の推定式ともみなせる
956:卵の名無しさん
20/01/31 17:45:07.01 24QiZJ0Y.net
Rv?1という条件はワクチン接種率の条件として書き直せばv?1?1R0=H
957:卵の名無しさん
20/02/01 15:12:17 hIisy8jC.net
H(n) = Σ[k=1,2,...,n] 1/k
とする。H(n)を既約分数で表したときの分子の整数をf(n)と表す。
(1）lim[n→∞] H(n) を求めよ、答えのみで良い。
(2）n=1,2,...に対して、f(n)に現れる1桁の整数を全て求めよ
"
H <- function(n) sum(1/(1:n))
plot(sapply(1:1000,H),bty='l')
f <- function(n,prec=10000){ # Σ 1/kを既約分数表示する n>>=23で誤計算
if(n==1){
cat(n, ':' ,1,'\n')
invisible(1)
}else{
GCD <- function(a,b){ # ユークリッドの互除法
r = a%%b # a=bq+r ⇒ a%%b=b%%rで最大公約数表示
while(r!=0){a = b ; b = r ; r = a%%b}
b }
nn = 1:n # nn : 1 2 3 ... n
nume=list() # 分子の容器
length(nume)=n
Nume=0
library(gmp)
for(i in nn) Nume = Nume + prod.bigz(nn[-i]) # nnからi番目の要素を除いて乗算して総和を分子に
Deno=factorialZ(n) # 分母
gcd = GCD(Nume,Deno) # Numerator/Denominator約分するため最大公約数を計算
NUME=Nume/gcd
DENO=Deno/gcd
ratio=as.numeric(Nume/Deno)
RE=list(NUME,DENO,ratio) # 最大公約数で除算して約分
cat(n, ':' ,as.character(NUME),'/',as.character(DENO),'\n')
invisible(RE)
}}
958:卵の名無しさん
20/02/01 22:05:31.53 hIisy8jC.net
H(n) = Σ[k=1,2,...,n] 1/k
とする。H(n)を既約分数で表したときの分子の整数をf(n)と表す。
f <- function(n){ # Σ 1/kを既約分数表示する
if(n==1){
cat(n, ':' ,1,'\n')
invisible(1)
}else{
GCD <- function(a,b){ # ユークリッドの互除法
r = a%%b # a=bq+r ⇒ a%%b=b%%rで最大公約数表示
while(r!=0){a = b ; b = r ; r = a%%b}
b }
nn = 1:n # nn : 1 2 3 ... n
nume=list() # 分子の容器となるlist
length(nume)=n # を設定
Nume=0
library(gmp)
for(i in nn) Nume = Nume + prod.bigz(nn[-i]) # nnからi番目の要素を除いて乗算して総和を分子に
Deno=factorialZ(n) # 分母
gcd = GCD(Nume,Deno) # 約分するため最大公約数を計算
NUME=Nume/gcd
DENO=Deno/gcd
ratio=as.numeric(Nume/Deno)
RE=list(NUME,DENO,ratio)
cat(n, ':' ,as.character(NUME),'/',as.character(DENO),'\n')
invisible(RE)
}}
959:卵の名無しさん
20/02/02 13:53:23 +5dNqMpE.net
tbl <- function(x,v){ # vの要素がxにいくつあるか集計する
n=length(v)
hme=numeric(n) # how many entries?
for(i in 1:n) hme[i]=sum(x==v[i])
rbind(v,hme)
}
tbl(sample(10,rep=T),1:10)
960:卵の名無しさん
20/02/04 16:35:57 b64IHQrg.net
"URLﾘﾝｸ(this.kiji.is)
中国湖北省武漢市からチャーター機で日本へ帰国した邦人の新型コロナウイルス感染率が高いと、
中国で驚きの声が上がっている。中国当局が発表した同市の感染者の割合に比べ「39倍も高い」というのだ。
現地は医療現場が混乱しているため、実際には発表よりかなり多くの感染者がいる可能性がある。
日本政府はチャーター機計3便を武漢市に派遣し、邦人565人が帰国した。厚生労働省によると、
チャーター機に乗っていた感染者は、症状のない人も含め計8人。感染率は1.416％だ。
一方、1月31日現在、武漢市の感染者数は3215人で、感染率は0.036％にとどまった。"
r1=8
r2=3215
n1=565
n2=round(3215/(0.036/100))
poisson.test(c(r1,r2),c(n1,n2))
prop.test(c(r1,r2),c(n1,n2))
library(BayesFactor)
mat=cbind(c(r1,r2),c(n1,n2)-c(r1,r2))
fisher.test(mat)
chisq.test(mat)$p.value
bf=BayesFactor::contingencyTableBF(mat,sampleType = 'poisson',posterior = TRUE,iteration=1e5)
rbf=round(bf)
(r1/n1)/(r2/n2)
frp <- function(x) (x[1]/(x[1]+x[2]))/ (x[3]/(x[3]+x[4]))
rrp=apply(rbf,1,frp)
hist(rrp)
BEST::plotPost(rrp)
(r1/(n1-r1))/(r2/(n2-r2))
fop <- function(x) (x[1]/x[2])/ (x[3]/x[4])
orp=apply(rbf,1,fop)
BEST::plotPost(orp)
quantile(orp,c(0.005,0.025,0.5,0.975,0.995))
library(ex
961:act2x2) exact2x2::fisher.exact(mat) fisher.test(mat)
962:卵の名無しさん
20/02/04 16:57:56 b64IHQrg.net
# Reed-Frost Model 時刻t+1での感染者数=時刻ｔでの感受性人数*(少なくとも一人の感染者がでる確率)
# p=0.04;N=100;I0=1;T=8
ReedFrost=function(
p=8/565, # 1期間内での伝播確率
N=565, # 集団の人数
I0=8, # 最初の感染者数
T=30, # 全期間
print=TRUE # グラフ表示
)
{
q=1-p # 伝播させない確率
I=numeric(T) # 感染者数 newly Infecteds
S=numeric(T) # 感受性のある人数 Susceptibles
I[1]=I0 # I0人から始まったとする
S[1]=N-I[1]
for(t in 1:(T-1)){
I[t+1]=S[t]*(1-q^I[t]) # Reed-Frost Model
S[t+1]=S[t]-I[t+1]
}
if(print){
plot(1:T,I,type='b',pch=19,lwd=2, ylim=c(0,N),bty='l',
xlab="time",ylab="persons",main=paste("Reed-Frost Model p= ",round(p,3)))
points(S,lty=2,col=2,lwd=2,type='b',pch=19)
points(N-S,lty=3,col=3,lwd=2,type='b',pch=19)
legend("right",bty="n",legend=c("New Infection","Susceptible","Immunized"),lty=1:3,col=1:3,lwd=2)}
list(newly_infecteds=round(I),susceptibles=round(S))
}
ReedFrost(T=7)