臨床統計もおもしろいですよ、その3
- 1 ：卵の名無しさん：2020/03/05(木) 20:17:05.03 ID:naSB8128.net
-
内科認定医受験の最低限の知識、
製薬会社の示してくる臨床データ、
論文の考察、
論文を書くときの正当性、
というのが、臨床統計の今までの目的の大きい部分でしたが、
AI=機械学習の基本も、結局は統計学と確率に支配されます。
そういう雑多な話をするスレです。
※前スレ
臨床統計もおもしろいですよ、その１
https://egg.2ch.net/test/read.cgi/hosp/1493809494/
臨床統計もおもしろいですよ、その2
https://egg.5ch.net/test/read.cgi/hosp/1540905566/
- 2 ：卵の名無しさん：2020/03/05(木) 20:17:46.76 ID:naSB8128.net
- 100人の集団（クラスター）で感染者は１人として１日に延べ１０回・人と接触し、
１人１回あたりの感染確率を１％、感染期間30日、潜伏期5日として
SEIRモデルで計算すると、感染のピークは110日めでとても１～２週間が山場とは言えない。
> SEIR2(contact_rate=10,transmission_probability=0.01
+ ,infectious_period=30,latent_period=5,mu=0,
+ nu=0, s=99,e=0,i=1,r=0,timepoints = seq(0,365,by=0.5),axes=TRUE)
Ro = 3 peak time I = 109.5 peak time E = 89
グラフにすると
https://i.imgur.com/EiG2HzI.jpg
- 3 ：卵の名無しさん：2020/03/05(木) 22:34:57.14 ID:vGY9zXsV.net
- https://i.imgur.com/obIkECO.jpg
- 4 ：卵の名無しさん：2020/03/06(金) 09:29:34.46 ID:H1NcHDj6.net
- 病的な虚言癖と妄想癖の精神科医 古根高の病名を診断するスレ
https://egg.5ch.net/test/read.cgi/hosp/1529634250/
古根高（サッポロファクトリーメンタルクリニック院長）（五浪で底辺医大に進学した知恵遅れ）のご尊顔
https://i.imgur.com/u8rBskz.jpg
トラブルメーカーゆえ市立札幌病院を解雇された過去あり
https://i.imgur.com/upWjhiA.jpg
- 5 ：卵の名無しさん：2020/03/07(土) 09:19:47.77 ID:jrmaVufM.net
- 週末はこれで遊べる。
nCov2019 for studying COVID-19 coronavirus outbreak
https://guangchuangyu.github.io/nCov2019/
- 6 ：卵の名無しさん：2020/03/07(土) 10:32:32.16 ID:jrmaVufM.net
- >>5
まず、Rを最新版に更新して
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("BiocStyle")
remotes::install_github("GuangchuangYu/nCov2019", dependencies = TRUE)
が必要だった。
- 7 ：卵の名無しさん：2020/03/07(土) 10:47:18.70 ID:jrmaVufM.net
- BiocStyle-defunct {BiocStyle} R Documentation
Defunct functions in package ‘BiocStyle’
Description
These functions are defunct and no longer available.
Details
The following functions are no longer available; use the replacement indicated below:
latex_old, latex2: latex
pdf_document_old, pdf_document2: pdf_document
html_document_old, html_document2: html_document
- 8 ：卵の名無しさん：2020/03/07(土) 15:19:10.32 ID:KnbHZezX.net
- >>1
受験を控えている高校生です
他のスレでアナタはコンプ性Ｆラン症候群という不治の病であるとお聞きしました
高学歴の私と違ってＦランであるアナタは今はどのようなお仕事をされてるんですか？
- 9 ：卵の名無しさん：2020/03/07(土) 19:47:32.16 ID:jrmaVufM.net
- Coronavirus COVID-19 outbreak statistics and forecast
http://www.bcloud.org/e/
これは使えるな
http://www.bcloud.org/e/session/a9b2a4a3e4ef04951a3e45743536a8f3/download/dataDownloadWorld?w=
- 10 ：卵の名無しさん：2020/03/08(日) 14:49:32.83 ID:WlKU3B2d.net
- SEIRモデルで有病率を１％に固定して、集団のサイズを変化させてシミュレーションしてみたけどピークは変わらないな。
このモデルでは集会規模の大小には影響されないということになるな。
https://i.imgur.com/343K91V.png
有病率を変化させて流行の変遷をグラフにすると、
https://i.imgur.com/SZ15LKT.png
https://i.imgur.com/gnJVFnd.png
有病率を４０％くらいに引き上げるとオリンピックのときには流行が収束していることになるｗｗｗｗ
- 11 ：卵の名無しさん：2020/03/08(日) 19:19:22.20 ID:WlKU3B2d.net
- rm(list=ls())
graphics.off()
f <- function(n,p) 1-(1-p)^n
vf=Vectorize(f)
n=seq(1,1000,by=1)
plot(n,vf(n,0.001),bty='l')
abline(h=0.5,lty=3)
"
藤井聡が 集会の参加人数と感染拡大確率の表を公開している。
://i.imgur.com/f3Tes5g.jpg
同じアルゴリズムを使って有病率が１/10000のときに感染拡大確率を0.05未満にしたい。
何人以上の集会を禁じることによってそれが達成できるか計算せよ。
"
F <- function(p=0.001,n=1000,p0=0.5){ # p=有病率,n:最大参加人数,p0:感染拡大確率
sub <- function(n,p) 1-(1-p)^n
uniroot(function(x) sub(x,p)-p0,c(0,n))$root
}
F()
F(p=1e-4,n=1e4,p0=0.05)
- 12 ：卵の名無しさん：2020/03/08(日) 19:19:38.25 ID:WlKU3B2d.net
- "
p0= 1-(1-p)^n
(1-p)^n=1-p0
n*log(1-p)=log(1-p0)
n=log(1-p0)/log(1-p)
"
pp02n <- function(p,p0) log(1-p0)/log(1-p) # p=有病率,p0:感染拡大確率
pp02n(1/1000,0.05)
p=1/10^seq(1,5,by=0.01)
p=rev(p)
plot(p,pp02n(p,0.05),log='x',type='l',bty='l',ylab='参加者数',xlab='有病率',
main='感染拡大確率別許容参加人数')
lines(p,pp02n(p,0.10),lty=2,lwd=2)
lines(p,pp02n(p,0.25),lty=3,lwd=3)
lines(p,pp02n(p,0.50),lty=4,lwd=4)
legend('top',bty='n',legend=c(0.05,0.10,0.25,0.5),lty=1:4,lwd=1:4)
"
p0= 1-(1-p)^n
(1-p)^n=1-p0
log(1-p)=log(1-p0)/n
1-p=exp(log(1-p0)/n)
p=1-exp(log(1-p0)/n)
"
n=465
p0=0.05
np02p <- function(n,p0) 1-exp(log(1-p0)/n)
np02p(465,0.05)
- 13 ：卵の名無しさん：2020/03/08(日) 20:01:02.96 ID:WlKU3B2d.net
- "
東京都は２１日、新型コロナウイルスの感染拡大を防ぐため、２２日から３月１５日までの３週間、都が主催する５００人以上の大規模な屋内イベントは、原則延期か中止にする方針を発表した.
"
by=1e-5
p0=seq(by,0.5,by)
plot(np02p(500,p0),p0,log='x',bty='l',type='l',lwd=2,main='500人参加集会',
xlab='有病率(log)',ylab='感染拡大確率')
- 14 ：卵の名無しさん：2020/03/09(月) 08:18:17.40 ID:cZTNWJLE.net
- Same Stats, Different Graphs: Generating Datasets with Varied Appearance and Identical Statistics through Simulated Annealing Justin Matejka and George Fitzmaurice Autodesk Research, Toronto Ontario Canada {first.last}@autodesk.com
- 15 ：卵の名無しさん：2020/03/09(月) 08:21:02.64 ID:cZTNWJLE.net
- https://blogs.sas.com/content/iml/2019/04/17/create-version-of-anscombes-quartet.html
- 16 ：卵の名無しさん：2020/03/09(月) 08:23:47.32 ID:cZTNWJLE.net
- >>14
https://www.researchgate.net/profile/Justin_Matejka/publication/
316652618_Same_Stats_Different_Graphs_Generating_Datasets_with_Varied_Appearance_and_Identical_Statistics_through_Simulated_Annealing/links/
599aef990f7e9b892bacff30/Same-Stats-Different-Graphs-Generating-Datasets-with-Varied-Appearance-and-Identical-Statistics-through-Simulated-Annealing.pdf
- 17 ：卵の名無しさん：2020/03/09(月) 08:24:55.87 ID:cZTNWJLE.net
- https://youtu.be/It4UA75z_KQ
- 18 ：卵の名無しさん：2020/03/09(月) 19:01:23.07 ID:Ost2Ton5.net
- Sys.setenv(lang='en')
rm(list=ls())
library(gtools)
n=4
a=permutations(n,n) # ｎの順列
r=nrow(a)
f<-function(x){ # x=c(1,2,3) -> rbind(a[1],a[2],a[3])
n=length(x)
ans=NULL
for(i in 1:n){
ans=rbind(ans,a[i,])
}
return(ans)
}
(b=combn(r,n,f)) # すべての行列の配列 b[,,i]
(c=dim(b)[3]) # その個数
diag2 <- function(x){ # 右上から左下への対角線の配列を返す
n=nrow(x)
ans=numeric(n)
for(i in 1:n) ans[i]=x[i,n+1-i]
return(ans)
}
- 19 ：卵の名無しさん：2020/03/09(月) 19:01:30.17 ID:Ost2Ton5.net
- g <- function(x){ # 列、対角線の要素がすべて異なればTRUEを返す
n=nrow(x)
if(length(unique(diag (x))) < n) return(FALSE)
if(length(unique(diag2(x))) < n) return(FALSE)
flag=TRUE
for(i in 1:n){
if(length(unique(x[,i])) < n){
flag=FALSE
break
}
}
return(flag)
}
count=0
for(i in 1:c){
count=count+g(b[,,i])
}
count
- 20 ：卵の名無しさん：2020/03/09(月) 20:15:13.42 ID:Ost2Ton5.net
- >>18
これバグがあるな
- 21 ：卵の名無しさん：2020/03/09(月) 21:27:07.44 ID:Ost2Ton5.net
- >>20
バグ修正版
https://rio2016.5ch.net/test/read.cgi/math/1501755792/789
- 22 ：卵の名無しさん：2020/03/10(火) 17:40:25.16 ID:Fd2PNckm.net
- n=1e3
p0=runif(n,0,1)
oz0=p0/(1-p0)
pLR=0.7/(1-0.9) #TP/FP
nLR=0.3/0.9 # FN/TN
oz1=oz0*nLR^2*pLR^2
p1=oz1/(1+oz1)
quantile(p1,c(.025,0.5,.95))
- 23 ：卵の名無しさん：2020/03/10(火) 18:13:44.48 ID:Fd2PNckm.net
- sn=0.7
sp=0.9
n=1e7
p0=runif(n,0,1)
oz0=p0/(1-p0)
pLR=sn/(1-sp) #TP/FP
nLR=(1-sn)/sp # FN/TN
oz1=oz0*nLR^2*pLR^2
p1=oz1/(1+oz1)
BEST::plotPost(p1,showMode = T)
HDInterval::hdi(p1)
quantile(p1,c(.025,0.5,.975))
summary(p1)
MAP <- function(x) {
dens <- density(x)
mode_i <- which.max(dens$y)
mode_y <- dens$y[mode_i]
c(x=mode_x, y=mode_y)
}
MAP(p1)
- 24 ：卵の名無しさん：2020/03/11(水) 07:17:27.96 ID:w+Rtdm4+.net
- PCR2prob <- function(
sn=0.7, # sensitivity
sp=0.9, # specificity
plus=2, # how many positive result?
minus=2, # how many negative result?
n=1e7,
print=TRUE) # how large the simulation
{
p0=runif(n,0,1)
oz0=p0/(1-p0) # prob -> odds
pLR=sn/(1-sp) # TP/FP
nLR=(1-sn)/sp # FN/TN
oz1=oz0*pLR^plus*nLR^minus # Bayesian formula
p1=oz1/(1+oz1) # odds -> prob
if(print){
BEST::plotPost(p1,showMode =T) # show mode instead of mean
print(HDInterval::hdi(p1)) # Highest Density Interval
print(quantile(p1,c(.025,0.5,.975))) # 95%CI by quantile
print(summary(p1)) # mean, median
}
invisible(p1)
}
PCR2prob(n=1e5,minus=3)
PCR2prob(n=1e5,minus=4)
- 25 ：卵の名無しさん：2020/03/11(水) 07:17:32.52 ID:w+Rtdm4+.net
- MAP <- function(x) {
dens <- density(x)
mode_i <- which.max(dens$y)
mode_x <- dens$x[mode_i]
mode_y <- dens$y[mode_i]
c(x=mode_x, y=mode_y)
}
MAP(p1)['x'] # show mode
f <- function(x) mean(PCR2prob(minus=x,n=1e4,print=F))
plot(sapply(1:10,f),bty='l',pch=19)
abline(h=c(0.05,0.1,0.5),col='gray',lty=3)
- 26 ：卵の名無しさん：2020/03/11(水) 10:07:10.64 ID:NmX2BxF2.net
- PCR2prob <- function(
sn=0.7, # sensitivity
sp=0.9, # specificity
plus=2, # how many positive result?
minus=2, # how many negative result?
n=1e7, # how large the simulation
p0=runif(n,0,1),
print=FALSE)
{
oz0=p0/(1-p0) # prob -> odds
pLR=sn/(1-sp) # TP/FP
nLR=(1-sn)/sp # FN/TN
oz1=oz0*pLR^plus*nLR^minus # Bayesian formula
p1=oz1/(1+oz1) # odds -> prob
if(print){
BEST::plotPost(p1,showMode =T) # show mode instead of mean
print(HDInterval::hdi(p1)) # Highest Density Interval
print(quantile(p1,c(.025,0.5,.975))) # 95%CI by quantile
print(summary(p1)) # mean, median
}
invisible(p1)
}
mean(PCR2prob(n=1e5,minus=3))
PCR2prob(n=1e5,minus=4)
- 27 ：卵の名無しさん：2020/03/11(水) 10:55:05.81 ID:NmX2BxF2.net
- PCR2prob <- function(
sn=0.7, # sensitivity
sp=0.9, # specificity
plus=2, # how many positive result?
minus=2, # how many negative result?
n=1e7, # how large the simulation
p0=rbeta(n,1,1), # prior distribution
print=TRUE)
{
oz0=p0/(1-p0) # prob -> odds
pLR=sn/(1-sp) # TP/FP
nLR=(1-sn)/sp # FN/TN
oz1=oz0*pLR^plus*nLR^minus # Bayesian formula
p1=oz1/(1+oz1) # odds -> prob
if(print){
BEST::plotPost(p1,showMode =T) # show mode instead of mean
print(HDInterval::hdi(p1)) # Highest Density Interval
print(quantile(p1,c(.025,0.5,.975))) # 95%CI by quantile
print(summary(p1)) # mean, median
}
invisible(p1) # return posterior distribution
}
- 28 ：卵の名無しさん：2020/03/11(水) 15:05:50.80 ID:w+Rtdm4+.net
- sim <- function(m=2){
count=0
i=0
while(count<m){
i=i+1
count = count + sample(6,1)==1
}
return(i)
}
k=1e5
re=replicate(k,sim())
tbl=table(re) ; tbl
which.max(tbl)
plot(tbl/k,bty='l')
- 29 ：卵の名無しさん：2020/03/11(水) 17:01:09.37 ID:w+Rtdm4+.net
- "1が累計m(=2)回出るまでサイコロを振って、振った回数を当てるギャンブルがある。
何回目に賭けるのがベストか？"
sim <- function(m=2){
pip1=0 # １の目の出た回数
i=0 # サイコロを振った回数
while(pip1 < m){
i=i+1
pip1 = pip1 + (sample(6,1)==1)
}
return(i)
}
k=1e5
re=replicate(k,sim(3))
tbl=table(re) ; tbl
which.max(tbl)
plot(tbl/k,bty='l')
- 30 ：卵の名無しさん：2020/03/11(水) 17:01:26.29 ID:w+Rtdm4+.net
- #
bg <- function(x,print=FALSE){ # big gambling
f <- function(n,m=x,p=1/6) choose(n-1,m-1)*p^(m-1)*(1-p)^(n-m)*p
nn=1:(10*x)
y=optimize(function(n) f(n),nn,maximum=TRUE)$maximum
if(print){
plot(nn,sapply(nn,f),bty='l',pch=19)
yy=c(floor(y),ceiling(y))
cat(c(f(yy[1]),f(yy[2])),'\n')
}
return(floor(y))
}
x=2:100
y=sapply(x,bg)
plot(x,y,bty='l')
(lm=lm(y~x))
summary(lm)
abline(lm)
y <- function(x) 6*x-6
y(1000)
- 31 ：卵の名無しさん：2020/03/12(木) 20:47:51.16 ID:Hoqw8S6P.net
- Rに移植
allocate.rooms <- function(m,n){ # m:rooms n:people
if(m==n) return(factorial(m))
else if(m==1) return(1)
else m*Recall(m,n-1) + m*Recall(m-1,n-1)
}
Cに移植
#include<stdio.h>
long factorial(long n) {
long re = 1;
long k;
for(k=1;k <=n;k++) {re *= k;}
return re;
}
long rooms(int m, int n){
if(m==n) { return factorial(m);}
else if(m==1){ return 1;}
else{
return m * rooms(m,n-1) + m * rooms(m-1,n-1);
}
}
void main( int argc, char *argv[] ){
int m,n;
long ways;
m=atoi(argv[1]);
n=atoi(argv[2]);
ways=rooms(m,n);
printf("%d\n",ways);
}
- 32 ：卵の名無しさん：2020/03/12(木) 20:48:19.73 ID:Hoqw8S6P.net
- m部屋n人空部屋なしの場合の数で
a[m,n] = a[m,n-1]×m + a[m-1,n-1]×m
*Main> let s m n = if m == n then (product [1..m]) else if m==1 then 1 else m*(s (m) (n-1)) + m*(s (m-1) (n-1))
*Main> s 6 12
953029440
- 33 ：卵の名無しさん：2020/03/14(土) 12:34:11.45 ID:E6fZCvx4.net
- ド底辺シリツ医大の新入生100人のうち10人は学力考査で合格した正規入学生で残りは裏口入学生であるとする。
誰が正規入学かを確かめる必要がある。審議官の女医（嘘つきかどうかは不明）が無作為に新入生を呼び出して
「あなたのいうことが正しければ手コキかフェラをしてあげる、そうでなければセンズリを命じる」という試問をする。
論理的思考ができる正規入学生はこの質問に回答してフェラをしてもらっているが、裏口入学生は正解がだせずにセンズリを命じられている。
正規入学生10人を同定できたら一連の試問は終了とする。
センズリを命じられる裏口入学生の期待値を求めよ。
rm(list=ls())
students=rep(1:2,c(10,90))
picked=NULL
flag=FALSE
sim <- function(){
while(flag==FALSE){
i=sample(length(students),1)
picked=c(picked,students[i])
students=students[-i]
flag=sum(picked==1)==10
}
sum(picked==2)
}
k=1e5
mean(replicate(k,sim()))
900/11
- 34 ：卵の名無しさん：2020/03/14(土) 14:06:57.29 ID:E6fZCvx4.net
- "
検査陽性 検査陰性
カレー頻食 a(=18) b(=30)
カレー稀食 c(=32) d(=20)
Usage
riskratio(X, Y, m1, m2, conf.level=0.95, p.calc.by.independence=TRUE)
Arguments
X The number of disease occurence among exposed cohort. # a
Y The number of disease occurence among non-exposed cohort. # c
m1 The number of individuals in exposed cohort group. # a+b
m2 The number of individuals in non-exposed cohort group. # c+d
"
curry <- function(
prev=0.10, # prevalence
N=1, # population
sn=0.7, # sensitivity
sp=0.9, # specificity
a=18,b=30,
c=32,d=20)
{
PPV=N*prev*sn/(N*prev*sn+N*(1-prev)*(1-sp))
NPV=N*(1-prev)*sp/(N*(1-prev)*sp+N*prev*(1-sn))
library(fmsb)
p0=riskratio(a,c,a+b,c+d)$p.value
p1=riskratio(a*PPV,c*PPV,a*PPV+b*NPV,c*PPV+d*NPV)$p.value
c(p0,p1)
}
curry()
- 35 ：卵の名無しさん：2020/03/14(土) 14:27:13.93 ID:E6fZCvx4.net
- 某国の新型コロナ感染症の有病率を0.1、PCR検査の感度を0.7 特異度を0.9とする
検査陽性陰性の人を無作為に50人ずつ集めて、カレーを頻回に食べているかを調査した結果が
以下の通りであった。
検査陽性 検査陰性
カレー頻食 a(=18) b(=30)
カレー稀食 c(=32) d(=20)
カレーを頻食すると新型コロナ感染に罹りにくいと結論できるか？ 危険率0.05で検定せよ。
- 36 ：卵の名無しさん：2020/03/14(土) 14:33:35.06 ID:E6fZCvx4.net
- >>35
有病率によって有意かどうか判断が変わるのは興味ぶかい。
https://i.imgur.com/qjjIjnm.png
- 37 ：卵の名無しさん：2020/03/14(土) 20:12:22.41 ID:E6fZCvx4.net
- 症数例から断定的な判断をするのは危険なので
周辺度数を固定してMCMCしてｐ値の信頼区間を出してみた。
#
prev=0.10 # prevalence
sn=0.7 # sensitivity
sp=0.9 # specificity
a=18
b=30
c=32
d=20
(x=matrix(c(a,b,c,d),2, byrow=TRUE))
N=sum(x)
PPV=N*prev*sn/(N*prev*sn+N*(1-prev)*(1-sp))
NPV=N*(1-prev)*sp/(N*(1-prev)*sp+N*prev*(1-sn))
library(BayesFactor)
bf=contingencyTableBF(x,sampleType = 'indepMulti',fixedMargin = 'cols',
posterior = TRUE, iteration=1e4)
head(bf)
N=sum(x)
a=N*bf[,1]
b=N*bf[,2]
c=N*bf[,3]
d=N*bf[,4]
library(fmsb)
p.value=riskratio(a*PPV,c*PPV,a*PPV+b*NPV,c*PPV+d*NPV)$p.value
BEST::plotPost(p.value,compVal = 0.05)
https://i.imgur.com/w2RzDLB.png
ｐ値の期待値は0.152
有意水準0.05で有意差ありとされる確率は約4割という結果がえられた。
- 38 ：卵の名無しさん：2020/03/14(土) 20:16:36.23 ID:E6fZCvx4.net
- >>37
タイプミス修正
症数例から
↓
少数例
- 39 ：卵の名無しさん：2020/03/15(日) 08:45:40.74 ID:LJXDSnIu.net
- https://id.fnshr.info/2019/07/20/r3-6-0/
で知った
例
rep("1:6",5)
paste(rep("1:6",5),collapse=',')
str2lang(paste("expand.grid(",paste(rep("1:6",5),collapse=','),')'))
eval(str2lang(paste("expand.grid(",paste(rep("1:6",5),collapse=','),')')))
(s=rep("1:6",5))
(str=paste(s,collapse=','))
(lang=str2lang(paste("expand.grid(",str,')')))
eval(lang)
- 40 ：卵の名無しさん：2020/03/15(日) 08:53:22.23 ID:LJXDSnIu.net
- ｆｎ ＜－ ｆｕｎｃｔｉｏｎ（ｎ）｛
ｌｉｂｒａｒｙ（ｇｍｐ）
ｇｒ＝ｅｖａｌ（ｓｔｒ２ｌａｎｇ（ｐａｓｔｅ（＂ｅｘｐａｎｄ．ｇｒｉｄ（＂，ｐａｓｔｅ（ｒｅｐ（＂１：６＂，ｎ），ｃｏｌｌａｐｓｅ＝＇，＇），＂）＂）））
ｒ＝ａｓ．ｂｉｇｑ（ｓｕｍ（ａｐｐｌｙ（ｇｒ，１，ｆｕｎｃｔｉｏｎ（ｘ） ｐｒｏｄ．ｂｉｇｚ（ｘ）％％６＝＝０）））／ａｓ．ｂｉｇｑ（ｎｒｏｗ（ｇｒ））
ｒ２＝ｃａｐｔｕｒｅ．ｏｕｔｐｕｔ（ｒ）［２］
ｓｕｂｓｔｒ（ｒ２，５，ｎｃｈａｒ（ｒ２））
｝
ｆｏｒ（ｉ ｉｎ １：９） ｃａｔ（ｉ，＇：＇，ｆｎ（ｉ），＇\ｎ＇）
＃ ｓｔｒ２ｌａｎｇ ｓｔｒ２ｅｘｐｒｅｓｓｉｏｎ
？ｓｔｒ２ｌａｎｇ
ｒｅｐ（＂１：６＂，５）
ｐａｓｔｅ（ｒｅｐ（＂１：６＂，５），ｃｏｌｌａｐｓｅ＝＇，＇）
ｓｔｒ２ｌａｎｇ（ｐａｓｔｅ（＂ｅｘｐａｎｄ．ｇｒｉｄ（＂，ｐａｓｔｅ（ｒｅｐ（＂１：６＂，５），ｃｏｌｌａｐｓｅ＝＇，＇），＇）＇））
ｅｖａｌ（ｓｔｒ２ｌａｎｇ（ｐａｓｔｅ（＂ｅｘｐａｎｄ．ｇｒｉｄ（＂，ｐａｓｔｅ（ｒｅｐ（＂１：６＂，５），ｃｏｌｌａｐｓｅ＝＇，＇），＇）＇）））
（ｓ＝ｒｅｐ（＂１：６＂，５））
（ｓｔｒ＝ｐａｓｔｅ（ｓ，ｃｏｌｌａｐｓｅ＝＇，＇））
（ｌａｎｇ＝ｓｔｒ２ｌａｎｇ（ｐａｓｔｅ（＂ｅｘｐａｎｄ．ｇｒｉｄ（＂，ｓｔｒ，＇）＇）））
ｅｖａｌ（ｌａｎｇ）
- 41 ：卵の名無しさん：2020/03/15(日) 16:22:13.79 ID:5ItQZoKJ.net
- >>1
受験を控えている女子高校生です
他のスレでアナタはコンプ性Ｆラン症候群という不治の病であるとお聞きしました
高学歴の私と違ってＦランであるアナタは今はどのようなお仕事をされてるんですか？
- 42 ：卵の名無しさん：2020/03/15(日) 18:08:53.16 ID:LJXDSnIu.net
- >41
高学歴ならこれに答えてみ！
2020/03/13 11:00時点で東京都で1524人検査して87人陽性と報告されている
# https://data-science.gr.jp/theory/tpd_negative_binomial_distribution.html
発見率を87/1524は一定仮定して
新型コロナ陽性患者を10人集めたいとする。
必要な被検者数の期待値と95%信頼区間を求めよ。
# Monte Carlo
p=87/1524
failure=0:500
success=10
fail=rnbinom(1e6,success,p)
hist(fail,freq=F,col='seashell2',main='',xlab='failures',ylab='probability',breaks=30)
lines(failure,dnbinom(failure,success,p),lwd=2)
mean(fail) + success
success*(1-p)/p + success # mean(fail)=kq/p
HDInterval::hdi(x)
quantile(x,c(0.025,0.5,0.975))
qnbinom(c(0.025,0.975),success,p)
- 43 ：卵の名無しさん：2020/03/15(日) 18:09:31.68 ID:LJXDSnIu.net
- シミュレーション解の方が回数上限なしでプログラムできるな。
# Simulation
sim <- function(){
i=0
s10=0
while(s10!=10){
i=i+1
s10 = s10 + sample(1:0,1,prob=c(p,1-p))
s10 == 10
}
i
}
k=1e4
re=replicate(k,sim())
mean(re)
HDInterval::hdi(re)
- 44 ：卵の名無しさん：2020/03/15(日) 18:33:54.47 ID:LJXDSnIu.net
- 2020/03/13 11:00時点で東京都で1524人検査して87人陽性と報告されている
# https://data-science.gr.jp/theory/tpd_negative_binomial_distribution.html
発見率は平均値87/1524、0.01-0.10の区間の二項分布に従うと仮定する。
新型コロナ陽性患者を10人集めたい。
必要な被検者数の期待値と95%信頼区間を求めよ。
library(BayesFactor)
bf=BayesFactor::proportionBF(87,1524, p, posterior = TRUE,iter=1e4,nullInterval = c(0.01,0.10))
head(bf)
pp=(bf[,'p'])
plot(pp,bty='l')
mean(pp)
BEST::plotPost(pp,showMode = T)
qq=1-p
r=10
trials=r*qq/pp + r
BEST::plotPost(trials)
HDInterval::hdi(trials)
- 45 ：卵の名無しさん：2020/03/15(日) 18:40:27.27 ID:LJXDSnIu.net
- >>41
ド底辺頭脳でも数くらい数えられるんだろ？
1つのサイコロを10000回投げ、出たの目の積が6の倍数になる確率を分数でもとめよ。
- 46 ：卵の名無しさん：2020/03/15(日) 19:19:43.78 ID:5ItQZoKJ.net
- >>45
どうでもいいから早くコンプ薬屋を呼んで来なよ
そんなだからモテないんだお
by 都内高学歴ＪＫ
- 47 ：卵の名無しさん：2020/03/15(日) 19:31:13.08 ID:LJXDSnIu.net
- >>46
数も数えられんの？ガイジ？
- 48 ：卵の名無しさん：2020/03/15(日) 21:05:46.74 ID:LJXDSnIu.net
- >>35
curry <- function(
prev=0.10, # prevalence
sn=0.7, # sensitivity
sp=0.9, # specificity
a=18,b=30,
c=32,d=20)
{
N=1 # population
PPV=N*prev*sn/(N*prev*sn+N*(1-prev)*(1-sp))
NPV=N*(1-prev)*sp/(N*(1-prev)*sp+N*prev*(1-sn))
library(fmsb)
p0=riskratio(a,c,a+b,c+d)$p.value
a=a*PPV+b*(1-NPV)
b=b*NPV+a*(1-PPV)
c=c*PPV+d*(1-NPV)
d=d*NPV+c*(1-PPV)
p1=riskratio(a,c,a+b,c+d)$p.value
c(p0,p1)
}
curry()
prevs=seq(0.01,0.99,by=0.01)
plot(prevs,sapply(prevs,function(x) curry(prev=x)[2]),bty='l',type='l',lwd=2,
ylim=c(0,1),ylab='p-value',xlab='prevalence')
curry(a=36,c=64,b=60,d=40)
- 49 ：卵の名無しさん：2020/03/15(日) 21:06:08.12 ID:LJXDSnIu.net
- prev=0.10 # prevalence
sn=0.7 # sensitivity
sp=0.9 # specificity
a=18
b=30
c=32
d=20
N=sum(x)
PPV=N*prev*sn/(N*prev*sn+N*(1-prev)*(1-sp))
NPV=N*(1-prev)*sp/(N*(1-prev)*sp+N*prev*(1-sn))
a=a*PPV+b*(1-NPV)
b=b*NPV+a*(1-PPV)
c=c*PPV+d*(1-NPV)
d=d*NPV+c*(1-PPV)
(x=matrix(c(a,b,c,d),2, byrow=TRUE))
library(BayesFactor)
bf=contingencyTableBF(x,sampleType = 'indepMulti',fixedMargin = 'cols',
posterior = TRUE, iteration=1e4)
head(bf)
N=sum(x)
a=N*bf[,1]
b=N*bf[,2]
c=N*bf[,3]
d=N*bf[,4]
library(fmsb)
a=a*PPV+b*(1-NPV)
b=b*NPV+a*(1-PPV)
c=c*PPV+d*(1-NPV)
d=d*NPV+c*(1-PPV)
p=riskratio(a,c,a+b,c+d)$p.value
p.value=as.vector(p)
BEST::plotPost(p.value,compVal = 0.05)
- 50 ：卵の名無しさん：2020/03/15(日) 21:10:18.08 ID:LJXDSnIu.net
- F-value = 2/(1/sensitivity + 1/PPV ) 感度とPPVの調和平均
- 51 ：卵の名無しさん：2020/03/16(月) 00:14:20.50 ID:RRM1UH+L.net
- # n個からr個を選んで得られる順列の総数をP(n, r)とする. 任意のr>1に対して, P(n, r)は平方数でないことを示せ.
rm(list=ls())
library(gmp)
library(Rmpfr)
fn <- function(n){
r=2:n
a=chooseZ(n,r)*factorialZ(r) # nPr
b=mpfr(a,1e5)
re=is.whole(sqrt(b))
r[which(re)]
}
fn(1000)
- 52 ：卵の名無しさん：2020/03/16(月) 06:35:35.55 ID:7zdTC761.net
- >>51
早くＦラン仲間のコンプ薬屋を呼んで来なよ
お勉強してくれるって言ってるお
by 都内高学歴ＪＫ
- 53 ：卵の名無しさん：2020/03/16(月) 06:40:26.18 ID:RRM1UH+L.net
- >>52
ド底辺頭脳でも数くらい数えられるんだろ？
1つのサイコロを10000回投げ、出たの目の積が6の倍数になる確率を分数でもとめよ。
- 54 ：卵の名無しさん：2020/03/16(月) 07:32:37.13 ID:RRM1UH+L.net
- 国内で患者数が大幅に増えたときに備えた医療提供体制の確保について
今後、国内で新型コロナウイルス感染症患者数が大幅に増えたときに備え、各都道府県、保健所設置市及び特別区（以下「都道府県等」という。）における外来を受診する患者数等について、
以下の数式を用いて計算いただき、ピーク時の医療需要の目安としてご活用の上、必要な医療提供体制を確保していただくようお願いいたします。
（１）（ピーク時において1 日あたり新たに新型コロナウイルス感染症を疑って外来を受診する患者数）＝
（0-14 歳人口）×0.18／100＋（15-64 歳人口） ×0.29／100＋（65 歳以上人口） ×0.51／100
（２）（ピーク時において1 日あたり新型コロナウイルス感染症で入院治療が必要な患者数）＝
（0-14 歳人口）×0.05／100＋（15-64 歳人口）×0.02／ 100＋（65 歳以上人口） ×0.56／100
（３）（ピーク時において1 日あたり新型コロナウイルス感染症で重症者として治療が必要な患者数）＝
（0-14 歳人口）×0.002／100＋（15-64 歳人口） ×0.001／100＋（65 歳以上人口） ×0.018／100
注１）ピーク時は、各都道府県等において疫学的関連性が把握できない程度に感染が拡大した時点から概ね３か月後に到来すると推計されている。ただし、公衆衛生上の対策を行うことにより、ピークが下がるとともに後ろ倒しされる。
注２）重症者とは、集中治療や人工呼吸器を要する管理が必要な患者を指す。
注３）当該計算式は、都道府県等の単位以下における医療提供体制を確保するためのものであるとともに、各都道府県等によってピークを迎える時期が異なるため、
全国の人口を用いて計算することや単純に各自治体が算出するピークの数値を足し合わせることは、不適切な取扱いとなることに留意いただきたい。なお、当該計算式については、
今後新たな知見等により変更される可能性がある。
注４）実際には、ピーク時に至るまでの日々の患者数の増加はばらつきがあり、増加曲線は推計通りの形にならない可能性が高いため、
現実の患者の発生動向も踏まえて適切に体制を確保することが必要。
注５）当該計算式については、今後新たな知見等により変更される可能性がある。
- 55 ：卵の名無しさん：2020/03/16(月) 07:32:47.99 ID:RRM1UH+L.net
- v=c(0.18,0.29,0.51,
0.05,0.02,0.56,
0.002,0.001,0.018)
(mat=matrix(v*100,3,byrow=T))
- 56 ：卵の名無しさん：2020/03/16(月) 07:35:19.73 ID:7zdTC761.net
- 早くＦラン仲間のコンプ薬屋を呼んで来なよ
童貞仲間でしょぉお
by 都内高学歴ＪＫ
- 57 ：卵の名無しさん：2020/03/16(月) 07:54:48.43 ID:RRM1UH+L.net
- v=c(0.18,0.29,0.51,
0.05,0.02,0.56,
0.002,0.001,0.018)
(mat=matrix(v*100,3,byrow=T))
# https://www.toukei.metro.tokyo.lg.jp/juukiy/2019/jy19qa0200.xls
(x=matrix(c(1601348,9035668,3103714)/1e4,ncol=1))
round(mat%*%x)
> round(mat%*%x)
[,1]
[1,] 44915
[2,] 19989
[3,] 681
重症者の5割が死亡するとすると、都内で340人の死者予想。
- 58 ：卵の名無しさん：2020/03/16(月) 07:56:56.15 ID:RRM1UH+L.net
- >>56
コンプはこれに答がだせなかったんだが、あんたは出せる？
解答できなきゃ、コンプと同レベルの頭脳と認定してあげよう。
24時間待ってあげる。
問題
某国の新型コロナ感染症の有病率を0.1、PCR検査の感度を0.7 特異度を0.9とする
検査陽性陰性の人を無作為に50人ずつ集めて、カレーを頻回に食べているかを調査した結果が
以下の通りであった。
検査陽性 検査陰性
カレー頻食 a(=18) b(=30)
カレー稀食 c(=32) d(=20)
カレーを頻食すると新型コロナ感染に罹りにくいと結論できるか？ 危険率0.05で検定せよ。
- 59 ：卵の名無しさん：2020/03/16(月) 08:12:57.28 ID:RRM1UH+L.net
- 国内で患者数が大幅に増えたときに備えた医療提供体制の確保について
今後、国内で新型コロナウイルス感染症患者数が大幅に増えたときに備え、各都道府県、保健所設置市及び特別区（以下「都道府県等」という。）における外来を受診する患者数等について、
以下の数式を用いて計算いただき、ピーク時の医療需要の目安としてご活用の上、必要な医療提供体制を確保していただくようお願いいたします。
（１）（ピーク時において1 日あたり新たに新型コロナウイルス感染症を疑って外来を受診する患者数）＝
（0-14 歳人口）×0.18／100＋（15-64 歳人口） ×0.29／100＋（65 歳以上人口） ×0.51／100
（２）（ピーク時において1 日あたり新型コロナウイルス感染症で入院治療が必要な患者数）＝
（0-14 歳人口）×0.05／100＋（15-64 歳人口）×0.02／ 100＋（65 歳以上人口） ×0.56／100
（３）（ピーク時において1 日あたり新型コロナウイルス感染症で重症者として治療が必要な患者数）＝
（0-14 歳人口）×0.002／100＋（15-64 歳人口） ×0.001／100＋（65 歳以上人口） ×0.018／100
問題 ： ある都市でピーク時に外来1000人、入院600人、重症20人であったとすると、この都市の14歳以下の人口は何人と推測されるか？
- 60 ：卵の名無しさん：2020/03/16(月) 13:47:58.64 ID:ulqLVmOp.net
- >>58
評価するには症例数が２桁少ないですわね
それはさておき昔の教授選では積み上げたペーパーの高さを比べたと聞きますが
Ｆランのアナタは採択されたペーパーの１本もあるのかしら？Ｗ
by 都内Ｓラン女子高生
- 61 ：卵の名無しさん：2020/03/16(月) 14:14:20.21 ID:ncInpF91.net
- やったー 国試無事合格！
- 62 ：卵の名無しさん：2020/03/16(月) 23:11:32.64 ID:5fhQT99v.net
- >>60
足りないのは標本数でなくてオマエの頭だよ。
- 63 ：卵の名無しさん：2020/03/16(月) 23:31:54.26 ID:5fhQT99v.net
- 数日前、科学技術省の共同予防および制御メカニズムの科学研究グループからの良いニュースがありました。
国家緊急予防および制御薬物工学技術研究センターと深セン第三人民病院LeレイおよびLi英Yチームによって完成された「ファピラビル」。
「新しいコロナウイルス肺炎（COVID-19）患者の安全性と有効性に関する臨床研究」は、予備的な臨床結果に達しています（登録番号：
ChiCTR2000029600）。
研究により、ファビピラビルは、ウイルスクリアランスを促進することにより、新しいコロナウイルス肺炎の進行を緩和することが示されています。研究結果は、中国工学院のジ
ャーナルに提出されました
合計80人の患者がこの臨床試験に登録され、35人の患者がファピラビルで治療され、対照群は年齢、性別、および疾患の重症度が治療群と一致した新冠動脈肺炎の患者45人
でした。 iタブレット治療。薬物投与からウイルス除去までの時間の中央値、治療の14日目の胸部画像の改善率、および安全性を2つのグループ間で比較しました。
結果は、ウイルス除去の観点から、治療後のファピラビル試験群の患者のウイルス除去時間の中央値（ウイルス核酸陰性）は、それぞれ4日と11日であった対照群のそれよりも
有意に短いことを示しました。別の重要な指標として、患者の胸部画像の改善、治療群と対照群の改善率はそれぞれ91.43％と62.22％でした。同時に、対照群と比較し
て、フェラビル治療群は副作用が少なく、忍容性が良好でした。 2つの患者グループのベースライン特性はすべて同等でした。
- 64 ：卵の名無しさん：2020/03/16(月) 23:35:45.65 ID:5fhQT99v.net
- >>63
PCRの方は中央値だけでデータがないので統計処理は辿れないが、比率の方が計算できる。
> round(35*0.9143)
[1] 32
> round(45*0.6222)
[1] 28
>
> library(fmsb)
> (mat=cbind(c(32,28),c(35,45)-c(32,28)))
[,1] [,2]
[1,] 32 3
[2,] 28 17
> chisq.test(mat)
Pearson's Chi-squared test with Yates' continuity correction
data: mat
X-squared = 7.4667, df = 1, p-value = 0.006285
> fisher.test(mat)
Fisher's Exact Test for Count Data
data: mat
p-value = 0.003669
alternative hypothesis: true odds ratio is not equal to 1
95 percent confidence interval:
1.59236 37.25281
sample estimates:
odds ratio
6.33637
- 65 ：卵の名無しさん：2020/03/16(月) 23:36:27.82 ID:5fhQT99v.net
- MCMCして信頼区間で比較すると、
https://i.imgur.com/3uZNIwI.png
- 66 ：卵の名無しさん：2020/03/17(火) 00:11:08.71 ID:+SeXrvpf.net
- >>63
PCR陰転までの日数の中央値が
アビガン群３５例で４日、対照群で１１日としても
次のような分布なら有効とは言い難い。
https://i.imgur.com/TjJshMu.png
> summary(A)
Min. 1st Qu. Median Mean 3rd Qu. Max.
1.00 3.00 4.00 13.86 24.50 47.00
> summary(C)
Min. 1st Qu. Median Mean 3rd Qu. Max.
1.00 5.00 11.00 11.09 16.00 20.00
グラフ化なしで統計処理をするな、という例だな。
- 67 ：卵の名無しさん：2020/03/17(火) 00:15:40.58 ID:+SeXrvpf.net
- 中央値比較での統計悪用シミュレーションプログラム
sim <- function(){
A=c(sort(sample(1:4,17,rep=T)),4,sort(sample(5:50,17,rep=T)))
C=c(sort(sample(1:11,23,rep=T)),11,sort(sample(11:20,23,rep=T)))
print(summary(A))
print(summary(C))
par(mfrow=c(2,1))
plot(table(A),bty='l',xlim=c(1,50),col='red')
plot(table(C),bty='l',xlim=c(1,50),col='blue')
}
sim()
- 68 ：卵の名無しさん：2020/03/17(火) 00:19:03.10 ID:+SeXrvpf.net
- こんな感じ
https://i.imgur.com/kL3pnn3.png
https://i.imgur.com/3HYkd5O.png
- 69 ：卵の名無しさん：2020/03/17(火) 11:29:17.69 ID:qm/RkVMU.net
- >>64
ド底辺頭脳はサイコロを６回投げたら１の目が１回でているとう風な計算しかできんみたいだな。
信頼区間とか確率分布という概念での思考ができなんよね。
周辺度数固定して二項分布で乱数発生させてχ二乗検定したときのｐ値の分布は以下のようになる。
https://i.imgur.com/5djScD1.png
https://i.imgur.com/RsFr7o4.png
危険率１％で有効といえる確率はほぼ５割であることがわかる。
- 70 ：卵の名無しさん：2020/03/17(火) 11:38:17.66 ID:qm/RkVMU.net
- >>60
エントリーに５以下の数値があるとχ二乗検定だとYatesの補正を使うけど、この程度の標本数は臨床では普通だよ。
中国のアビガンの予備試験も標本数が実薬群ととコントロール群を合わせて80だぞ。
足りないのはオマエの学力だよ。
- 71 ：卵の名無しさん：2020/03/17(火) 12:02:20.58 ID:RqUjOCLx.net
- >>70
まあ 学園カースト最下位のＦランさんが必死ですわね
昨日は先輩方からも国試合格のご報告がたくさんありましたが
Ｆランさんのお知り合いの方たちはいかがですかぁ？Ｗ
by 都内Ｓラン女子高生
- 72 ：卵の名無しさん：2020/03/17(火) 12:41:00.33 ID:+SeXrvpf.net
- >>71
答えられない低学力を隠そうと必死だな。
どうやって解くかわからんの？
- 73 ：卵の名無しさん：2020/03/17(火) 19:55:56.01 ID:RqUjOCLx.net
- >>72
わかりました
犯人はアナタですね Ｆランさん
- 74 ：卵の名無しさん：2020/03/17(火) 21:23:49.50 ID:qm/RkVMU.net
- ようやく完成した
# 球面上の一様分布
vertex <- function(r=1){
x=runif(1,-1,1) # x ~ 一様分布[-1,1]
phi=runif(1,-pi,pi) # φ ~ 一様分布[-π,π]
y=sqrt(1-x^2)*cos(phi) # √(1-x^2)*cos(φ)
z=sqrt(1-x^2)*sin(phi) # √(1-x^2)*sin(φ)
r*c(x,y,z)
}
https://i.imgur.com/xX0mTim.png
https://i.imgur.com/H7hs9w8.png
# 凸包四面体の体積
sim <- function(r=1,print=F){
v4=replicate(4,vertex()) # ４点の直交座標
if(print) print(v4)
abs(det(v4[,2:4]-v4[,1]))/6 # 四面体の体積
}
sim(print=T)
k=1e5
tetra=replicate(k,sim()) # k回のシミュレーション
mean(tetra)
summary(tetra)
BEST::plotPost(tetra)
- 75 ：卵の名無しさん：2020/03/17(火) 21:46:10.60 ID:qm/RkVMU.net
- # 単位球体内部一様分布の４点で四面体をつくるとき
rm(list=ls())
vertex <- function(){
x=runif(1,-1,1) # x ~ 一様分布[-1,1]
phi=runif(1,-pi,pi) # φ ~ 一様分布[-π,π]
r=runif(1,0,1) # r ~ 一様分布[0,1]
x=r^(1/3)*x # r^(1/3)*x
y=r^(1/3)*sqrt(1-x^2)*cos(phi) # r^(1/3)√(1-x^2)*cos(φ)
z=r^(1/3)*sqrt(1-x^2)*sin(phi) # r^(1/3)√(1-x^2)*sin(φ)
c(x,y,z)
}
vtx=replicate(5000,vertex())
par(mfrow=c(3,1))
x=vtx[1,] ; hist(x,col='pink')
y=vtx[2,] ; hist(y,col='orange')
z=vtx[3,] ; hist(z,col='darkgreen')
rgl::plot3d(x,y,z, col="skyblue")
par(mfrow=c(1,1))
# 四面体の体積
sim <- function(r=1,print=F){
v4=replicate(4,vertex()) # ４点の直交座標
if(print) print(v4)
abs(det(v4[,2:4]-v4[,1]))/6 # 四面体の体積
}
sim(print=T)
k=1e5
tetra=replicate(k,sim()) # k回のシミュレーション
mean(tetra)
summary(tetra)
BEST::plotPost(tetra)
- 76 ：卵の名無しさん：2020/03/17(火) 21:47:20.28 ID:qm/RkVMU.net
- core i7 Memory 16G はシミュレーションが捗って( ・∀・)ｲｲ!!
- 77 ：卵の名無しさん：2020/03/20(金) 07:04:24.48 ID:LvYpiZ17.net
- 臨床医学は確率事象を扱う稼業なので統計処理は必須。
これができない医師はアホといえる。
某国の新型コロナ感染症の有病率を0.1、
PCR検査の感度を0.7 特異度を0.9とする
検査陽性陰性の人を無作為に100人ずつ集めて、カレーを頻回に食べているかを調査した結果が
以下の通りであった。
検査陽性 検査陰性
カレー頻食 a(=36) b(=60)
カレー稀食 c(=64) d(=40)
カレーを頻食すると新型コロナ感染に罹りにくいと結論できるか？
- 78 ：卵の名無しさん：2020/03/21(土) 10:28:11.17 ID:eZCI10Rh.net
- Parameter Total (n=83)
Severe/critica l Group (n=25)
Ordinary Group (n=58)
P Value
Number of lobes
5 (4-5) 5 (5-5) 5 (3-5) 0.003
- 79 ：卵の名無しさん：2020/03/21(土) 10:29:01.74 ID:eZCI10Rh.net
- median (interquartile range, IQR)
- 80 ：卵の名無しさん：2020/03/21(土) 17:41:24.46 ID:o3rnHhhm.net
- library(coin)
# coin::wilcox_test(c(Ea,La) ~ factor(c(rep('Ea',n1),rep('La',n2))))
f4 <- function(x){ # filer p.value=0.003
# pv=t.test(sc[x[1],],or[x[2],])$p.value
SC=sc[x[1],]
OR=or[x[2],]
w=coin::wilcox_test(c(SC,OR) ~ factor(c(rep('SC',25),rep('OR',58))),distribution="exact")
pv=coin::pvalue(w)
if(round(pv,3)==0.003) return(c(x[1],x[2]))
else(return(c(0,0)))
}
- 81 ：卵の名無しさん：2020/03/22(日) 09:18:45.02 ID:/FOXgWJQ.net
- PCR <- function(n=1000,r=10,sen=0.7,spc=0.9,k=1e6,print=TRUE){
# n:検査件数 r:検査陽性数 sen=感度 spc=特異度 k:発生乱数の数,print:グラフ描画
prev=rbeta(k,1+r,1+n-r) # 有病率の事前分布を一様分布と仮定
PPV=prev*sen/( prev*sen + (1-prev)*(1-spc) ) # 検査特性を加味してPPV計算
adjtd.prev = PPV*(1-spc)/(sen - PPV*(sen+spc-1)) # そのPPVから有病率を計算
# NPVでも同じ結果
# NPV=(1-prev)*spc/( (1-prev)*spc + prev*(1-sen)) # 検査特性を加味してNPV計算
# adjtd.prev = (1-NPV)*spc/(spc - NPV*(sen+spc-1)) # そのNPVから有病率を計算
mean=mean(adjtd.prev) # 期待値
median=median(adjtd.prev) # 中央値
mode=density(adjtd.prev)$x[which.max(density(adjtd.prev)$y)] # 最頻値
LU=unlist(HDInterval::hdi(adjtd.prev))[1:2] # 95%信頼区間
re=c(mean=mean,median=median,mode=mode,LU) # 事後有病率の代表値
if(print){ # ヒストグラム描画
par(mfrow=c(2,2))
hist(prev,freq=F,xlim=c(0,0.5), main='prevalence') ; lines(density(prev))
hist(PPV,freq=F,xlim=c(0,0.5),main='PPV') ; lines(density(PPV))
hist(adjtd.prev,freq=F,xlim=c(0,0.5),main='adjusted prevalence ') ; lines(density(adjtd.prev))
BEST::plotPost(adjtd.prev) ; lines(density(adjtd.prev),col='skyblue')
par(mfrow=c(1,1))
}
print(round(re,4)) # 概算値表示
invisible(list(re,adjtd.prev)) # 代表値と事後有病率を非表示で返す
}
PCR(100,1)
PCR(1000,10)
PCR(10000,100)
PCR(100000,1000)
- 82 ：卵の名無しさん：2020/03/22(日) 20:13:20.91 ID:/FOXgWJQ.net
- # stanでMCMCさせることにした。
data{ // corona.stan
int n; // 1000
int x; // 10
real<lower=0,upper=1> sen; // 0.7
real<lower=0,upper=1> spc; // 0.9
}
parameters{
real<lower=0,upper=1> prev; // prevalence
}
transformed parameters{
real<lower=0,upper=1> p;
p = prev*sen + (1-prev)*(1-spc) ; // probability of positive test result
}
model{
x ~ binomial(n,p);
prev ~ beta(1,1); // prev ~ uniform(0,1)
}
- 83 ：卵の名無しさん：2020/03/22(日) 20:14:04.17 ID:/FOXgWJQ.net
- stanはコンパイルに時間がかかって面倒。
jagsの方が早くて( ・∀・)ｲｲ!!
rm(list=ls())
graphics.off()
options(scipen = 4)
options(digits=5)
library(rjags)
N=1000 ; X=10
dataList=list(n=N,x=X,sen=0.7,spc=0.9)
modelstring=paste0('
model
{
prev ~ dbeta(1,1)
p <- prev*sen + (1-prev)*(1-spc)
x ~ dbin(p,n)
}
')
writeLines(modelstring,'TEMPmodel.txt')
jagsModel = jags.model( file="TEMPmodel.txt" ,data=dataList)
update(jagsModel)
codaSamples = coda.samples( jagsModel ,
variable=c("prev","p"), n.iter=50000 )
js=as.matrix(codaSamples)
BEST::plotPost(js[,'prev'])
round(HDInterval::hdi(js[,'prev'])[1:2],7)
- 84 ：卵の名無しさん：2020/03/23(月) 12:03:19.86 ID:a9wYdNxH.net
- 何の責任もとれないド底辺ほど
他人の責任にうるさく自分の権利ばかりしつこく要求しがちですわね
コンプ薬屋のことですけど
Fランのアナタも同類ですわよね
早く呼んでくればぁ
by 都内Ｓラン女子高生
- 85 ：卵の名無しさん：2020/03/23(月) 13:38:33.35 ID:a9wYdNxH.net
- >>1
よそのスレで騒いでるみたいですけど
Fランのアナタには知性があるとは
医師の方々は考えてはいないようですよ
by 都内Ｓラン女子高生
- 86 ：卵の名無しさん：2020/03/23(月) 14:45:03.80 ID:a9wYdNxH.net
- >>1
国試に行き詰まって統計遊びですか
こういうときハゲワロスって言うんでしたっけ？
by ドＳなＳラン女子高生
- 87 ：卵の名無しさん：2020/03/23(月) 14:48:48.41 ID:a9wYdNxH.net
- タバコを吸わないものは自室に禁煙の張り紙など貼らない訳ですが
裏口裏口と盛んにワメキ回るアナタには
身に覚えがあると言うことでよろしいですねＷ
容疑者はオマエだぁっ（コナ○口調で）
by ドＳなＳラン女子高生
- 88 ：卵の名無しさん：2020/03/23(月) 15:00:20.99 ID:a9wYdNxH.net
- アナタの期待値は100発0中だぁっ
ち○ぴょろすぽ○ん
by ドＳなＳラン女子高生
- 89 ：卵の名無しさん：2020/03/23(月) 18:19:32.21 ID:a9wYdNxH.net
- 自分で作った統計ごっこのスレがここにあるのに
なーんでよそのスレに迷惑をかけに行くのかなん
やっぱりＦラン脳のゆえんなのかしらねんのねん
by ドＳなＳラン女子高生
- 90 ：卵の名無しさん：2020/03/25(水) 09:16:02.76 ID:htH1Knlx.net
- >>88
これやってみ、ド底辺頭脳のオマエには無理だろうな。
24時間待ってやるぞ。
ドツボ13は100発0中
ドツボ14は10発0中
ドツボ15は1発0中
とする。
各々10000発撃ったときドツボの命中数の期待値はいくらか？
- 91 ：卵の名無しさん：2020/03/25(水) 09:16:47.51 ID:htH1Knlx.net
- 正確な感度も特異度もわからないから有病率は計算できないという主張が数学板でされたから、
感度・特異度を確率変数としてベイズ階層モデルを組めば計算できること投稿しておいた。
最頻値と標準偏差からベータ分布のパラメータを計算するのに連立方程式を解く方が面倒だった。
- 92 ：卵の名無しさん：2020/03/25(水) 14:59:53.46 ID:Jsfh+Odg.net
- きゃー、統計ごっこのクルクルパー来たーっW
（キモワロス）
ド底辺仲間のコンプ薬屋はまだですか？
by ドＳなＳラン女子高生
- 93 ：卵の名無しさん：2020/03/25(水) 15:55:21.10 ID:C3cAic+l.net
- >>92
>>88
これやってみ、ド底辺頭脳のオマエには無理だろうな。
24時間待ってやるぞ。
ドツボ13は100発0中
ドツボ14は10発0中
ドツボ15は1発0中
とする。
各々10000発撃ったときドツボの命中数の期待値はいくらか？
- 94 ：卵の名無しさん：2020/03/26(木) 10:22:15.80 ID:F34mrbtj.net
- >>93
国試に合格できなかった還暦ジジイって本当ですかぁ？
アチコチで反日レスしてるみたいですねー
by ドＳなＳラン女子高生
- 95 ：卵の名無しさん：2020/03/26(木) 10:47:49.20 ID:Z3UzEFJB.net
- >>94
ドツボの宿題まだぁ？
- 96 ：卵の名無しさん：2020/03/26(木) 11:18:27.82 ID:F34mrbtj.net
- >>95
Ｆラン高の宿題はＦランのアナタ自身でされてくださいね
Ｓランのアタシはアナタと違ってヒマではございませんのでＷ
by ドＳなＳラン女子高生
- 97 ：卵の名無しさん：2020/03/26(木) 19:27:21.07 ID:O1C2WieZ.net
- >>96
ドツボ問題もできないドアホ。
これやってみ、ド底辺頭脳のオマエには無理だろうな。
24時間待ってやるぞ。
ドツボ13は100発0中
ドツボ14は10発0中
ドツボ15は1発0中
とする。
各々10000発撃ったときドツボの命中数の期待値はいくらか？
- 98 ：卵の名無しさん：2020/03/27(金) 07:08:38.46 ID:BC/LnmbB.net
- 95人調べて47人陽性。感度50-70% 特異度90%前後として
この集団の有病率はstanでMCMCして計算すると
mean lower upper
0.7558614 0.5214634 0.9999897
感度60% 特異度90%に固定すると
> pn2ip(95,47)
infected prevalence
75.0000000 0.7894737
- 99 ：卵の名無しさん：2020/04/04(土) 05:13:53.27 ID:ApKEo+4Y.net
- 麻薬取締法違反 札幌の薬剤師と法人を略式起訴
病院で管理する医療用麻薬の数量について道に虚偽の届け出をしたとして、
札幌区検は１９日までに麻薬取締法違反の罪で、札幌ひばりが丘病院（札幌
市厚別区）に勤務していた３０代の薬剤師の男と、同病院を運営する医療法人
潤和会（同区）を略式起訴した。札幌簡裁がそれぞれ２０万円以下の罰金刑
を言い渡す見通し。略式命令請求によると、男は２０１５年１１月、道に対し、
院内で使用した麻薬の使用量を偽って届け出るなどしたとされる。潤和会には
法人も罰する同法の両罰規定を適用した。
出典：北海道新聞 平成30年10月19日付
- 100 ：卵の名無しさん：2020/04/07(火) 08:43:20.38 ID:je8pKp8V.net
- 晋型コロナ肺炎に感度0.9，特異度0.9の迅速検査が開発されたと仮定する。
日本人1億2595万人からX人を無作為抽出して有病率を推定したい。
有病率の99%信頼区間幅を1%以内で検定したい。
何人を抽出すれば十分といえるか？
ss2sbj <- function(SEN,SPC,range=0.01,CONF=0.99,print=FALSE){
pr2sbj <- function(p0,sen=SEN,spc=SPC,R=range,conf.level=CONF){
z = qnorm(1-(1-conf.level)/2)
p = p0*sen+(1-p0)*spc
n = 4*z^2*p*(1-p)/R^2
return(n)
}
if(print){
prevalence=seq(0,1,by=0.001)
plot(prevalence,sapply(prevalence,pr2sbj),
ylab='subjects',bty='l',type='l',lwd=2)}
optimize(function(x) pr2sbj(x,SEN,SPC), c(0,1) ,maximum = TRUE)
}
- 101 ：卵の名無しさん：2020/04/22(水) 06:24:37.56 ID:XcybEEeM.net
- "n(=100)人の中から無作為にm(=10)人選んだらその中に少なくとも一人の感染者がいた。
全体で何人の感染者がいるかの期待値と95%信頼区間を求めよ。"
rm(list=ls())
fn <- function(n,m,print=FALSE){
library(gmp)
# 離散量で計算
x=0:n # 感染者数:x, 非感染数:n-x
pmf=1- chooseZ(n-x,m)/chooseZ(n,m) # 1 - (ｍ人全員非感染の確率)
pdf=pmf/sum(pmf) # 確率密度関数化して
Ev=sum(x*pdf) # 期待値を計算
if(print) print(Ev)
E=as.numeric(Ev)
# 信頼区間計算計算のために連続量とする
y=seq(0,n,by=0.05) # infected
p=1- choose(n-y,m)/choose(n,m) # Pr[at least one infected among m]
d=p/sum(p) # pmf -> pdf
(e=sum(y*d)) # mean of probability
if(print){plot(y,d,bty='l',xlab='infected',ylab='density',type='l')
print(quantile(n*p,prob=c(0.025,0.05,0.5,0.95,0.975)))
print(HDInterval::hdi(n*p)[1:2])}
return(E)
}
fn(100,5,print=T)
n=100
m=1:n
plot(m,sapply(m,function(x) fn(n,x)), type='p',pch=19,bty='l',ylim=c(0,100),
xlab='有名人',ylab='期待値',main='総人口:100人')
- 102 ：卵の名無しさん：2020/04/22(水) 06:25:30.78 ID:XcybEEeM.net
- "
n(=10)人の中から無作為にm(=2)人選んだらその中に少なくとも一人の感染者がいた。
全体で何人の感染者がいるかの期待値を求めよ。
"
f <- function(p){
# p:感染確率
p1=2*p*(1-p) # 一人だけ感染確率
p2=p^2 # 二人とも感染確率
(1*p1+2*p2)/(p1+p2) # 感染人数の期待値
}
curve(f(x),bty='l',lwd=2,xlab='p')
f(2/3)
# 有名人が感染
" 長時間必要なのて実行非推奨,メモリ16G以下はクラッシュする
library(gmp)
m=18200 # 有名人の数(桜を見る会参加人数)
n=1.268e5 # 日本の人口
x=0:n # 感染者数:x, 非感染数:n-x
pmf=try(1- chooseZ(n-x,m)/chooseZ(n,m)) # 1 - (ｍ人全員非感染の確率)
pdf=try(pmf/sum(pmf)) # 確率密度関数化して
(E=try(sum(x*pdf))) # 期待値を計算
as.numeric(E) # E=63463.27 (m=1000) , E=1154070201/18202=63403.48(m=1.268e5)
"
fn(100,20,print=T)
data.frame(有名人=1:10,期待値=sapply(1:10,function(x) fn(100,x)$mean))
data.frame(有名人=1:10*10,期待値=sapply(1:10*10,function(x) fn(100,x)$mean))
- 103 ：卵の名無しさん：2020/04/22(水) 06:26:27.83 ID:XcybEEeM.net
- fn4 <- function(n,m,k,quick=TRUE,print=FALSE){
# Pr[B|Ax] n人中x人感染者がいるときにm人中k人感染している確率
library(gmp)
x=0:n
if(quick){ Ev=sum(x*chooseZ(x,k))/sum(chooseZ(x,k))
}else{
pmf=chooseZ(x,k)/chooseZ(n,m)
pdf=pmf/sum(pmf)
Ev=sum(x*pdf)}
# Σ(x*(xCk/nCm))/Σ(xCk/nCm) = Σ(x*(xCk))/Σ(xCk)
if(print) print(Ev)
E=as.numeric(Ev)
return(E)
}
fn4=Vectorize(fn4)
fn4(100,10,0:10,quick=F)
- 104 ：卵の名無しさん：2020/04/24(金) 13:08:57.65 ID:YTSFapOE.net
- 休校続きでお勉強くらいしかすることがございません
今思うのは反日野党が政権をとってたら
世界はモノスゴクひどいことになっているであろうこと
by 都内Ｓラン女子高生
- 105 ：卵の名無しさん：2020/04/26(日) 12:36:58.52 ID:7Vo1jkoW.net
- 今にして思えば菅直人は実に優秀だったな。
事故当日に原子力緊急事態宣言を発令、命をかけて原発に突入して手動ベントを約束させる、
自衛隊１０万に派遣決定、
全原発停止、
東電の撤退阻止
と、戦力の逐次投入でなくてもてる力をすべて出し切って日本を守ろうとした。
原発事故の最悪シナリオも専門家に依頼して作成して、それをパワーポイントでレクチャーさせて閣僚で危機感を共有していたという。
一方、安倍は菅直人が海水注入を止めたとデマを流しつづけた。名誉毀損裁判途中でデマメルマガをこっそりと削除。
- 106 ：卵の名無しさん：2020/04/27(月) 08:26:59.50 ID:Sri1jAWu.net
- 今日も朝からクルクルパーのFラン統計ジジイが
他スレでお医者さんごっこだぁW
今思うのは反日野党が政権をとってたら
日本はとっくに壊滅していたであろうこと
だって反日なんだもの
by 都内Ｓラン女子高生
- 107 ：卵の名無しさん：2020/04/30(木) 07:03:38.11 ID:yBA/Llsz.net
- 推定陽性者数：
> summary(js$n.pos.sum)
Min. 1st Qu. Median Mean 3rd Qu. Max.
3888 5094 5951 31450 8608 421156
> hdi(js$n.pos.sum)[1:2]
lower upper
3888 172503
- 108 ：卵の名無しさん：2020/04/30(木) 23:55:05.47 ID:jagIaU9/.net
- 統計ジジイとコンプ薬屋はＦラン仲間
今思うのは反日野党が政権をとってたら
日本はとっくに壊滅していたであろうこと
だって反日なんだもの
by 都内Ｓラン女子高生
- 109 ：卵の名無しさん：2020/05/01(金) 19:26:10.83 ID:8qgl7OUN.net
- いまだに安倍が日本人の敵と気づかないドアホかよ。
- 110 ：卵の名無しさん：2020/05/01(金) 20:58:05.28 ID:jwZPeRWl.net
- >>109
あぁあ Ｆランジジイだぁ
コテハンを命名してあげましょう
コンプにならってコロナ薬屋とでも名乗るが良いですよ
by 都内Ｓラン女子高生
- 111 ：卵の名無しさん：2020/05/02(土) 11:49:24.12 ID:1l93SIuE.net
- 今日も朝からクルクルパーのFラン統計ジジイが
他スレでお医者さんごっこだぁW
今思うのは反日野党が政権をとってたら
日本はとっくに壊滅していたであろうこと
だって反日なんだもの
by 都内Ｓラン女子高生
- 112 ：卵の名無しさん：2020/05/05(火) 12:24:58.63 ID:JWXeCObj.net
- 今日も朝からクルクルパーの
Fラン統計ジジイ あらため コロナ薬屋が
他スレでお医者さんごっこだぁW
今思うのは反日野党が政権をとってたら
日本はとっくに壊滅していたであろうこと
だって反日なんだもの
by 都内Ｓラン女子高生
- 113 ：卵の名無しさん：2020/05/05(火) 15:39:04.77 ID:Yw3eOPQR.net
- https://toyokeizai.net/sp/visual/tko/covid19/
のデータを使って
P=14895/153581 # 2020/05/04
PCR検査の感度30-70％のモデルM1と感度70-90％のモデルM2のどちらが信憑性があるか、ベイズファクターで計算してみる。
M1は感度が最頻値60%標準偏差10%、M2は最頻値80%標準偏差10%のベータ分布に設定
特異度はいずれも最頻値95%標準偏差2.5%に設定し、有病率は一様分布を仮定
陽性数は、陽性率(P)＝真陽性率＋偽陽性率=有病率=有病率*感度 + (1-有病率)*(1-特異度)の二項分布に従うとする。
事後確率分布は
://i.imgur.com/81bK4KE.png
陽性率P=14895/153581=0.09698465での事後確率分布の密度比(Savage-Dickey density ratio)でベイズファクターを出すと
> d1/d2 # Savage-Dickey densiti ratio = BF12
[1] 1.007722
まあ、ちょっぴり、感度30-70％のモデルの方がいいかも、という結果。
陽性数/検査数の時系列データでもあればもう少し差がでるかもしれん。
- 114 ：卵の名無しさん：2020/05/06(水) 13:43:13.77 ID:Pv7HZN2f.net
- 味覚障害がある患者をPCR検査に回したら陰性であった。
PCR検査の感度は30-70%特異度は95-100%として、この患者が非感染者である確率の期待値と信頼区間を計算せよ。
(P.AA_COVID=1-72/202)
# https://jamanetwork.com/journals/jama/fullarticle/2765183
sn=Mv2ab(0.5,0.2^2)
sp=Mv2ab(0.95,0.05^2)
data=list(sn=sn,sp=sp, P.AA_COVID=P.AA_COVID)
modelString='
model{
# P(COVID | anosmia and/or ageusia)
P.COVID_AA = P.AA_COVID*P.COVID/
(P.AA_COVID*P.COVID + P.AA_nonCOVID*(1-P.COVID))
# P(PCR- | aosmia and/or ageusia)
P.nPCR = (1-P.COVID_AA)*spc + P.COVID_AA*(1-sen)
# = TN + FN
P.COVID_nPCR= (1-sen)*P.COVID/((1-sen)*P.COVID + spc*(1-P.COVID))
# P(COVID | PCR-) = P(PCR-|COVID)P(COVID)/P(PCR-)
# = P(PCR-|COVID)P(COVID) /
# (P(PCR-|COVID)P(COVID)+P(PCR-|nonCOVID)P(nonCOVID))
# prior
P.AA_nonCOVID ~ dunif(0.056,0.127)
# https://jamanetwork.com/journals/jama/fullarticle/2765183
P.COVID ~ dbeta(1,1) # prevalence
sen ~ dbeta(sn[1],sn[2]) # sensitivity 0.3-0.7
spc ~ dbeta(sp[1],sp[2]) # specificity 0.9-1.0
}
'
- 115 ：卵の名無しさん：2020/05/06(水) 14:24:40.20 ID:vlDdeBIj.net
- 今日も朝からクルクルパーの
Fラン統計ジジイ あらため コロナ薬屋が
他スレでお医者さんごっこだぁW
ちゃんとコテハン付けないとダメでちよぉぉぉ
今思うのは反日野党が政権をとってたら
日本はとっくに壊滅していたであろうこと
だって反日なんだもの
by 都内Ｓラン女子高生
- 116 ：卵の名無しさん：2020/05/06(水) 14:44:32.79 ID:C9G8mFSP.net
- ボクちゃん
ちゃんとコテハン付けるんでしよぉぉぉお
漢字はまだ読めないんでしたっけねえぇ
by 都内Ｓラン女子高生
- 117 ：卵の名無しさん：2020/05/06(水) 18:04:44.54 ID:Pv7HZN2f.net
- # s1,s2,n1,n2 -> BF10
BF01 <- function(s1,s2,n1,n2){
library(gmp)
bf=chooseZ(n1,s1)*chooseZ(n2,s2)/chooseZ(n1+n2,s1+s2)*(n1+1)*(n2+1)/(n1+n2+1)
BF01=as.numeric(bf)
return(c('in favor of null-hyposisis : H0'=BF01))
}
(bf01=BF01(s1=424,s2=5416,n1=777,n2=9072))
c('in favor of alternative hyposis : H1' = 1/as.numeric(bf01))
- 118 ：卵の名無しさん：2020/05/08(金) 12:35:41.21 ID:h0xW2bYs.net
- 今日も朝から人生が暇しかないの
Fラン統計ジジイ あらため コロナ薬屋が
他スレでお医者さんごっこだぁW
ちゃんとコテハン付けないとダメでちよぉぉぉ
今思うのは反日野党が政権をとってたら
日本はとっくに壊滅していたであろうこと
だって反日なんだもの
by 都内Ｓラン女子高生
- 119 ：卵の名無しさん：2020/05/08(金) 12:44:01.53 ID:h0xW2bYs.net
- Fラン統計ジジイ あらため コロナ薬屋さん
および コンプ薬屋さん って
やっぱりバカなんですねぇぇえぇぇぇ！
ちゃんとコテハン付けないとダメでちよぉぉぉ
今思うのは反日野党が政権をとってたら
日本はとっくに壊滅していたであろうこと
だって反日なんだもの
by 都内Ｓラン女子高生
- 120 ：卵の名無しさん：2020/05/08(金) 13:28:12.61 ID:h0xW2bYs.net
- Fラン統計ジジイ あらため コロナ薬屋さん
漢字は読めないんでしゅかぁぁあぁぁW
ちゃんとコテハン付けないとダメでちよぉぉぉ
今思うのは反日野党が政権をとってたら
日本はとっくに壊滅していたであろうこと
だって反日なんだもの
by 都内Ｓラン女子高生
- 121 ：卵の名無しさん：2020/05/08(金) 15:02:48.93 ID:h0xW2bYs.net
- ところで コンプ薬屋という方は
転売生活でもしてらっしゃるのでしょうか？
http://itest.5ch.net/egg/test/read.cgi/hosp/1588253497/242
レスを遡ると品不足をアッピールして転売を煽るような書き込みをされてましたね
ということは お仲間の Ｆラン統計ジジイあらため コロナ薬屋さんも ご職業は転売ヤーですか？
by 都内Ｓラン女子高生
- 122 ：卵の名無しさん：2020/05/08(金) 18:33:14.28 ID:6cl7ipwe.net
- それで 転売ヤーの コロナ薬屋さん コンプ薬屋さんは
次はどんな品を転売予定ですか？
それはそうと
今思うのは反日野党が政権をとってたら
日本はとっくに壊滅していたであろうこと
だって反日なんだもの
by 都内Ｓラン女子高生
- 123 ：卵の名無しさん：2020/05/08(金) 19:48:15.39 ID:p2nXiTw0.net
- library(rjags)
s1=22;s2=10;n1=158;n2=78
c(s1/n1,s2/n2,s1/n1-s2/n2)
shape1=shape2=0.5
data=list(s1=s1,s2=s2,n1=n1,n2=n2,shape1=shape1,shape2=shape2)
modelString='
model{
# data
s1 ~ dbin(p1,n1)
s2 ~ dbin(p2,n2)
# model
p1 ~ dbeta(shape1,shape2)
p2 ~ dbeta(shape1,shape2)
delta = p1 - p2
# priors
pri.delta= pri.p1 - pri.p2
pri.p1 ~ dbeta(shape1,shape2)
pri.p2 ~ dbeta(shape1,shape2)
}
'
writeLines(modelString,'tmp.txt')
n.iter=1e4 ; thin=1
jagsModel=jags.model('tmp.txt',data=data,n.chains=4,n.adapt=n.iter/5)
update(jagsModel)
codaSamples=coda.samples(jagsModel,c('delta','pri.delta'),n.iter,thin)
plot(codaSamples)
js=as.data.frame(as.matrix(codaSamples))
BEST::plotPost(js$delta,compVal = 0)
- 124 ：卵の名無しさん：2020/05/08(金) 19:48:21.41 ID:p2nXiTw0.net
- library(polspline)
fit=logspline(js$delta)
pri.fit=logspline(js$pri.delta)
xlim=c(-1,1) ; ylim=c(0,9)
plot(fit,xlim=xlim, ylim=ylim,bty='l',lty=2)
par(new=T)
plot(pri.fit,xlim=xlim,ylim=ylim)
dpost=dlogspline(0,fit)
dprio=dlogspline(0,pri.fit)
points(0,dpost)
points(0,dprio,pch=19)
dpost/dprio
- 125 ：卵の名無しさん：2020/05/08(金) 20:57:51.65 ID:6cl7ipwe.net
- 転売ヤーの コロナ薬屋さん コンプ薬屋さん
それが次の転売品ですか？
高音ぼったくりで通販サイトを荒らして
恥ずかしくないのですか？
転売品仕入れの時は夜中から並んで野糞をしてるって本当ですか？
by 都内Ｓラン女子高生
- 126 ：卵の名無しさん：2020/05/08(金) 23:58:03.61 ID:6cl7ipwe.net
- http://itest.5ch.net/egg/test/read.cgi/hosp/1588253497/296
コロナ薬屋さん
コンプ薬屋さんに共感されてますよ
やっぱりお仲間だったんですね
キンモー
by 都内Ｓラン女子高生
- 127 ：卵の名無しさん：2020/05/09(土) 07:25:33.37 ID:OE35Ej5j.net
- 今日も朝から人生が暇だけの
コロナ薬屋 コンプ薬屋が
他スレでお医者さんごっこだぁW
ちゃんとコテハン付けないとダメでちよぉぉぉ
今思うのは反日野党が政権をとってたら
日本はとっくに壊滅していたであろうこと
だって反日なんだもの
by 都内Ｓラン女子高生
- 128 ：卵の名無しさん：2020/05/09(土) 08:18:45.43 ID:OE35Ej5j.net
- それで 私立薬学部卒の Ｆラン統計ジジイ
あらため コロナ薬屋さん
あなたと違って大学は旧帝受験予定ですが
それはそうと昨夜は何を転売するため並んでたんですか？
野糞しながらぼったくり転売は恥ずかしくないのですか？
by 都内Ｓラン女子高生
- 129 ：卵の名無しさん：2020/05/10(日) 00:00:00.36 ID:gWKVru57.net
- 全てのパチンコ屋を営業停止にして
コロナ患者の収容施設にしよう
- 130 ：卵の名無しさん：2020/05/10(日) 05:08:19.48 ID:1fOrGJlo.net
- 3人の女子大生に18種類の体位を教えてtヶ月後に覚えていた体位の数は以下の通りである.
?は欠測値。
記憶保持率をexp(-α*t)+βのモデルに当て嵌めて
4人めの冬子を含めて欠測値を推定せよ。
t=1 2 4 7 12 21 35 59 99 200
治子 18 18 16 13 9 6 4 4 4 ?
奈津子 17 13 9 6 4 4 4 4 4 ?
亜希子 14 10 6 4 4 4 4 4 4 ?
冬子 ? ? ? ? ? ? ? ? ? ?
- 131 ：卵の名無しさん：2020/05/10(日) 08:38:11.59 ID:hThV3HgC.net
- それで 私立薬学部卒の Ｆラン統計ジジイ
あらため コロナ薬屋さん
お仕事は転売ヤーで間違いないですね
頭が悪いのに転売のお仕事の為に統計の勉強をされたんですね
by 都内Ｓラン女子高生
- 132 ：卵の名無しさん：2020/05/10(日) 08:55:40.13 ID:hThV3HgC.net
- 統計上は PCR件数と 死亡数は
正の相関関係がある と
なるほど PCR信者はバカなんですね
by 都内Ｓラン女子高生
- 133 ：卵の名無しさん：2020/05/11(月) 14:17:46.13 ID:8Fkcxw3I.net
- 今日も朝から人生が暇しかない
Ｆラン統計ジジイ あらため コロナ薬屋 が
他スレでお医者さんごっこだぁW
ちゃんとコテハン付けないとダメでちよぉぉぉ
今思うのは反日野党が政権をとってたら
日本はとっくに壊滅していたであろうこと
だって反日なんだもの
by 都内Ｓラン女子高生
- 134 ：卵の名無しさん：2020/05/11(月) 15:25:59.72 ID:8Fkcxw3I.net
- ネットを見ると、大学の９月入学を支持する高３生の泣き言、みたいなニュースがありますね
私の周りの高３生にはそんなことを支持する者は誰もいないのですが
どんなＦラン高の生徒の話なんでしょうね
受験は生モノ
勉強時間が足りないのや試験日に体調が悪くて不合格になるのは、本人の実力のうちでしょうに
by 都内Ｓラン女子高生
- 135 ：卵の名無しさん：2020/05/11(月) 19:32:29.89 ID:pfUv3Eic.net
- >>130
４人目の冬子の確率を面積比例させてみると
https://i.imgur.com/CHLo8FN.png
- 136 ：卵の名無しさん：2020/05/11(月) 20:28:09.06 ID:8Fkcxw3I.net
- えぇぇえぇぇぇっ
Ｆラン統計ジジイ あらため コロナ薬屋さんって
反日パヨクなんですかぁぁあ？
サイテーですね
下水のゾウリムシよりもはるかにサイッテーですね
ダンゴムシのごとくゴロゴロ転がりながら
海より深く反省しなさい
by 都内Ｓラン女子高生
- 137 ：卵の名無しさん：2020/05/11(月) 22:43:11.83 ID:pfUv3Eic.net
- >>130
わりと、面倒。
library(rjags)
t <- c(1, 2, 4, 7, 12, 21, 35, 59, 99, 200)
nt <- length(t)
slist <- 1:4
ns <- length(slist)
k <- matrix(c(18, 18, 16, 13, 9, 6, 4, 4, 4, NA,
17, 13, 9, 6, 4, 4, 4, 4, 4, NA,
14, 10, 6, 4, 4, 4, 4, 4, 4, NA,
NA, NA, NA, NA,NA,NA,NA,NA,NA, NA), nrow=ns, ncol=nt, byrow=T)
colnames(k)=paste0('t',as.character(t))
rownames(k)=paste0('subject',1:4)
n <- 18
dataList=list(t=t,k=k,nt=nt,ns=ns,n=n)
- 138 ：卵の名無しさん：2020/05/11(月) 22:43:15.93 ID:pfUv3Eic.net
- cat("model
{
for(i in 1:ns){
for(j in 1:nt){
k[i,j] ~ dbin(theta[i,j],n)
predk[i,j] ~ dbin(theta[i,j],n)
}
}
for(i in 1:ns){
for(j in 1:nt){
theta[i,j] <- min(1 - 1e-16, exp(-alpha[i]*t[j]) + beta[i])
} # min(1, ....) => Node inconsistent with parents -- ERROR
}
for(i in 1:ns){
alpha[i] ~ dbeta(shape1a,shape2a) # dnorm(mu.a,lambda.a)T(0,1)
beta[i] ~ dbeta(shape1b,shape2b) # dnorm(mu.b,lambda.b)T(0,1)
}
shape1a ~ dt(0,pow(2.5,-2),1)T(0,) # mu.a ~ dbeta(1,1)
shape2a ~ dt(0,pow(2.5,-2),1)T(0,) # mu.b ~ dbeta(1,1)
shape1b ~ dt(0,pow(2.5,-2),1)T(0,) # lambda.a ~ dgamma(0.001,0.001)T(0.001,)
shape2b ~ dt(0,pow(2.5,-2),1)T(0,) # lambda.b ~ dgamma(0.001,0.001)T(0.001,)
}
", file='tmp.txt')
- 139 ：卵の名無しさん：2020/05/11(月) 22:44:15.44 ID:pfUv3Eic.net
- pk4=NULL
for(i in 1:10){
pk4=cbind(pk4,eval(str2lang(paste0("js$\'predk[4," ,i, "]\'"))))
}
colnames(pk4)=c(paste0('t0',1:9),'t10')
boxplot(pk4)
head(pk4)
df4=as.data.frame(pk4)
head(df4)
predk4=tidyr::pivot_longer(df4,col=1:10)
head(predk4)
colnames(predk4)=c('time','items')
library(ggplot2)
g <- ggplot(predk4,aes(time,items))
g + geom_violin()
g + geom_boxplot()
- 140 ：卵の名無しさん：2020/05/11(月) 22:45:23.11 ID:pfUv3Eic.net
- reshape2 が tidyr::pivot_longer(df4,col=1:10) になっていた。
使っていないと、ggplot2も忘れてしまうなぁ。
- 141 ：卵の名無しさん：2020/05/12(火) 07:14:59.85 ID:oz4s+cD5.net
- g <- ggplot(predk4,aes(time,items))
g + geom_boxplot()
https://i.imgur.com/krliw0d.png
バイオリンプロットというよりナメクジプロットだな。
g + geom_violin()
https://i.imgur.com/fYZHA53.png
histの$densityを使って面積比で反映させてみた
https://i.imgur.com/bIigw6R.png
標準のboxplotでも十分な気がする。
https://i.imgur.com/L210Sf8.png
- 142 ：卵の名無しさん：2020/05/12(火) 07:58:45.50 ID:InIDJ33x.net
- えぇぇえぇぇぇっ
Ｆラン統計ジジイ あらため コロナ薬屋さん
および コンプ薬屋さんって
裏口容疑者なんですかぁぁあ？
だと思いましたぁぁあぁぁ
早く出頭するんですよ
あなたたちのレスの全てから
頭の悪さを感じますからねえぇぇぇ
by 都内Ｓラン女子高生
- 143 ：卵の名無しさん：2020/05/12(火) 08:31:27.05 ID:oz4s+cD5.net
- 交絡因子を考慮しての記憶減衰の階層ベイズモデル
デバッグできたので、事前分布をいろいろ変えてみて再現性を確認。
内視鏡休診で纏まった勉強する時間が取れて( ・∀・)ｲｲ!!
##### psychophysical function with confounder
cat('
model{
for(i in 1:nsubjs){
for(j in 1:nstim[i]){
z[i,j] ~ dbern(phi[i])
pi[i,j] ~ dunif(0,1)
r[i,j] ~ dbin(theta[i,j],n[i,j])
theta[i,j] <-
ifelse(z[i,j]==1,ilogit(alpha[i] + beta[i]*(x[i,j]-xmean[i])),pi[i,j])
}
}
for(i in 1:nsubjs){
phi[i] <- phi(probit)
alpha[i] ~ dnorm(mua,lambdaa)
beta[i] ~ dnorm(mub,lambdab)
}
probit ~ dnorm(muphi,lambdaphi)
muphi ~ dnorm(0,0.001)
mua ~ dnorm(0,0.001)
mub ~ dnorm(0,0.001)
lambdaphi ~ dt(0,pow(2.5,-2),1)T(0,) # dunif(0.3)
lambdaa ~ dt(0,pow(2.5,-2),1)T(0,) # dunif(0,1000)
lambdab ~ dt(0,pow(2.5,-2),1)T(0,) # dunif(0,1000)
}
',file='tmp.txt')
- 144 ：卵の名無しさん：2020/05/12(火) 12:01:40.73 ID:InIDJ33x.net
- 開業医ではないあなたが
先に開業医板を荒らすから
私がここに降臨して差し上げてるんですよ
それすら判らないくらい知能が低いんですかねぇ
女子高生にブチ切れる異常者のジジイって
どんなホラー映画の地獄絵図ですかＷ
by 都内Ｓラン女子高生
- 145 ：卵の名無しさん：2020/05/12(火) 17:39:16.61 ID:oz4s+cD5.net
- >>143
half cauchy にしても一様分布にしても事後分布はほとんど変化はないな。
モデルが優れているからだろう。
- 146 ：卵の名無しさん：2020/05/12(火) 18:12:20.47 ID:InIDJ33x.net
- さて今日も学校は終わり。
お家に帰ってお勉強しましょう。
あ～あ、Ｆラン統計ジジイは
やっぱりくるくるぱーだったんですねぇ
by 都内Ｓラン女子高生
- 147 ：卵の名無しさん：2020/05/12(火) 23:05:47.83 ID:oz4s+cD5.net
- # 相関係数の差の検定
n1=150
r1=0.29
n2=120
r2=0.5
pnorm(-abs(atanh(r1)-atanh(r2))/sqrt(1/(n1-3)+1/(n2-3)))
- 148 ：卵の名無しさん：2020/05/13(水) 05:47:04.95 ID:7qkgP5oh.net
- 離散量のプロットは重なってしまうと相関がわかりにくくなる。
jitterつけると重なりは防げるな。
頻度に応じて面積比でプロット
https://i.imgur.com/QVrj8hW.png
n=10 ; r=500
x=sample(n,r,rep=T,prob=rbeta(n,0.5,0.5))
y=sample(n,r,rep=T,prob=rbeta(n,0.5,0.5))
plot(rep(0,24),pch=1:24,cex=2,axes=F,ann=F,ylim=c(0,0.1)) ;axis(1,at=1:24)
xy=paste0(x,y)
cex=numeric()
for(i in 1:length(x)){
cex[i]=sum(xy[i]==xy)
}
layout(matrix(1:4,2,b=T))
par(mar=c(1,1,1,1))
plot(x,y,pch=16)
plot(x,jitter(y),pch=16)
plot(jitter(x),jitter(y),pch=16)
plot(x,y,pch=16,cex=sqrt(cex))
- 149 ：卵の名無しさん：2020/05/13(水) 06:15:19.70 ID:7qkgP5oh.net
- 色の濃淡に頻度を反映させるという方法もあるな。
https://i.imgur.com/TnXbBKY.png
- 150 ：卵の名無しさん：2020/05/13(水) 06:16:02.77 ID:7qkgP5oh.net
- とりあえず、完成。多分、もっと便利なパッケージがあるんだろうな。
graphics.off()
oldpar=par()
par(mar=c(1,1,1,1),oma=c(3,3,3,3),bty='l',ann=F,pch=16)
plot(rep(0,24),pch=1:24,cex=2,axes=F,ann=F,ylim=c(0,0.1)) ;axis(1,at=1:24)
n=10 ; r=250
x=sample(n,r,rep=T,prob=abs(rnorm(n)))
y=sample(n,r,rep=T,prob=rbeta(n,0.5,0.5))
plot(x,y)
xy=paste0(x,y)
cex=numeric()
for(i in 1:length(x)){
cex[i]=sum(xy[i]==xy)
}
layout(matrix(1:4,2,b=T))
plot(jitter(x),jitter(y))
plot(x,y,cex=sqrt(cex))
plot(x,y,cex=3,col=rgb(0.1,0.1,0.1,0.25))
plot(x,y,cex=sqrt(cex),col=rgb(0.1,0.1,0.1,0.5))
layout(1)
par(oldpar)
- 151 ：卵の名無しさん：2020/05/13(水) 12:26:32.87 ID:GtsY/MQT.net
- 学校の先生に習いましたが、統計的にはPCR検査はむやみにやらないのが正解なんですね
by 都内Ｓラン女子高生
- 152 ：卵の名無しさん：2020/05/13(水) 14:17:48.98 ID:7qkgP5oh.net
- 【国民が誤解】安倍晋三「相談・受診の目安＝PCR実施の基準ではない。国民への周知が足りなかった」 [スタス★]
https://asahi.5ch.net/test/read.cgi/newsplus/1589284722/
- 153 ：卵の名無しさん：2020/05/13(水) 14:18:56.03 ID:7qkgP5oh.net
- 世界の頭脳は実態把握には大量検査が必要とわかっているなぁ。
東大の児玉名誉教授も同じ意見だな。
https://www.covid19-yamanaka.com/cont4/14.html
- 154 ：卵の名無しさん：2020/05/13(水) 14:55:45.98 ID:7qkgP5oh.net
- >>150
色を変えて散布図を作成
https://i.imgur.com/l1eyD6o.png
- 155 ：卵の名無しさん：2020/05/13(水) 14:57:54.70 ID:7qkgP5oh.net
- 検診がメインだったから、内視鏡バイトは休診。
俺は６割の休業補償だが、銭ゲバ開業医はパート医を首にしているんだろうなぁ。
- 156 ：卵の名無しさん：2020/05/13(水) 15:27:10.62 ID:GtsY/MQT.net
- 今日も朝から人生が暇しかない
Ｆラン統計ジジイ あらため コロナ薬屋 が
発狂して連投だぁぁぁW
ちゃんとコテハン付けないとダメでちよぉぉぉ
今思うのは反日野党が政権をとってたら
日本はとっくに壊滅していたであろうこと
だって反日なんだもの
by 都内Ｓラン女子高生
- 157 ：卵の名無しさん：2020/05/13(水) 15:48:20.06 ID:7qkgP5oh.net
- >>156
内視鏡バイトは休診なので、６割の休業補償はあるけど。
あんたのところは、感染拡大阻止に休診にして職員には休業補償してんの？
まとまった時間がとれるので、https://bayesmodels.com/の英文テキストで独学中。
行き詰まったら数学版で相談。
- 158 ：卵の名無しさん：2020/05/13(水) 18:18:49.08 ID:7qkgP5oh.net
- 相関係数->FisherのZ変換->ベイズファクター
r=js$r
#The Fisher z-transformation:
r.F <- atanh(r)
fit.posterior <- logspline(r.F)
plot(fit.posterior)
posterior <- dlogspline(atanh(0), fit.posterior)
r.F.prior <- atanh(runif(n.iter,-1,1))
fit.prior <- logspline(r.F.prior)
prior <- dlogspline(atanh(0), fit.prior)
BF10 <- prior/posterior
BF10
- 159 ：卵の名無しさん：2020/05/14(木) 09:44:43.41 ID:coRQ6//C.net
- まあ、ベイズファクターの計算結果はさほど変わらんな。
> BEST::plotPost(js$r,compVal = 0,showCurve = F,col='lightgreen')
> lines(density(js$r),col='green')
> abline(h=0.5,col=4)
> fit.post=logspline(js$r)
> points(0,dlogspline(0,fit.post),pch=19,col='white',cex=1.5)
> points(0,dlogspline(0,fit.post),cex=1.5)
> points(0,dunif(0,-1,1),pch=19,cex=1.5)
> curve(dlogspline(x,fit.post),add=T)
> fit.F=logspline(atanh(r))
> curve(dlogspline(x,fit.F), col=2, add=T)
> d.prio=dunif(0,-1,1)
> d.post=dlogspline(0,fit.post)
> d.post/d.prio
[1] 2.917773
> dlogspline(atanh(0),fit.F)/d.prio
[1] 2.911407
- 160 ：卵の名無しさん：2020/05/14(木) 09:45:47.27 ID:coRQ6//C.net
- >>159
グラフにしてみた。
https://i.imgur.com/OyTi5Jb.png
- 161 ：卵の名無しさん：2020/05/14(木) 10:03:25.31 ID:coRQ6//C.net
- あんまり相関係数の信頼区間など気にしていなかった。
必要ならmcmcかbootstrapすれば出せるし、
corr.testのソースを
sink('cor.text.txt')
stats:::cor.test.default
sink()
で覗いてみた。
if (n > 3) {
if (!missing(conf.level) && (length(conf.level) !=
1 || !is.finite(conf.level) || conf.level < 0 ||
conf.level > 1))
stop("'conf.level' must be a single number between 0 and 1")
conf.int <- TRUE
z <- atanh(r)
sigma <- 1/sqrt(n - 3)
cint <- switch(alternative, less = c(-Inf, z + sigma *
qnorm(conf.level)), greater = c(z - sigma * qnorm(conf.level),
Inf), two.sided = z + c(-1, 1) * sigma * qnorm((1 +
conf.level)/2))
cint <- tanh(cint)
FisherのZ変換して正規分布近似して信頼区間を出して、それをZ変換の逆関数で出しているな。
- 162 ：卵の名無しさん：2020/05/14(木) 10:39:23.42 ID:coRQ6//C.net
- >>161
sample関数を使ってのbootstrap
https://i.imgur.com/gAIJTEQ.png
- 163 ：卵の名無しさん：2020/05/14(木) 11:21:55.13 ID:coRQ6//C.net
- cat('
"名","本人","家族含む"
"鳩山由紀夫首相",144269,144269
"菅直人副総理兼国家戦略担当相",905,2232
"原口一博総務相",914,1220
"千葉景子法相",3523,3523
"岡田克也外相",3273,8641
"藤井裕久財務相",14356,20214
"川端達夫文部科学相",4024,5583
"長妻昭厚生労働相",0,891
"赤松広隆農相",4864,5934
"直嶋正行経済産業相",3333,3333
"前原誠司国土交通相",741,1441
"小沢鋭仁環境相",2089,4014
"北沢俊美防衛相",309,609
"平野博文官房長官",1195,1875
"中井洽国家公安委員長",1296,1296
"亀井静香金融・郵政担当相",9427,18745
"福島瑞穂消費者・少子化担当相",12734,25000
"仙谷由人行政刷新担当相",1968,3987
', file='shisan.csv')
X=read.csv("shisan.csv")
X
library(boot)
md = boot(X, function(df,idx) median(df[idx,2]), R=1e4)
plot(md)
mn = boot(X, function(df,idx) mean(df[idx,2]), R=1e4)
plot(md)
boxplot(list(median=md$t,mean=mn$t),horizontal = TRUE, las=1)
- 164 ：卵の名無しさん：2020/05/14(木) 11:33:51.60 ID:coRQ6//C.net
- >>162
パッケージ bootをつかうと
library(boot)
df=data.frame(x,k)
bt=boot(df,function(df,idx) cor(df[idx,1],df[idx,2]),R=1e5)
plot(bt)
HDInterval::hdi(bt$t)[1:2]
BEST::plotPost(as.vector(bt$t),xlab=bquote(cor),compVal = 0)
> HDInterval::hdi(bt$t)[1:2]
[1] -0.06736828 0.31237556
https://i.imgur.com/vCeo1Yn.png
- 165 ：卵の名無しさん：2020/05/14(木) 11:37:20.35 ID:coRQ6//C.net
- >>156
アベノミクスでの経済成長率は震災があった民主の頃より低いという現実を認められないんだなぁ。
日本を破壊したいなら、安倍一択。
安倍の長期政権で中国も韓国も笑いが止まらないはず、
日本が衰退する政策ばかりを選択しているから。
野党が反感をもたれて安倍信者が増えることこそ日本破壊が進んで
パヨクの思う壺
- 166 ：卵の名無しさん：2020/05/17(日) 15:29:50.58 ID:H2bLOh/W.net
- 結局、これが核心部分
Rt =I t (number of new infections generated at time t) / Σ[s=1,t] I t-s * Ws ( = total infectiousness of infected individuals at time t)
Ws : an infectivity profile given by a probability distribution ws, dependent on time since infection of the case, s, but independent of calendar time, t.
E[I t] = Rt Σ[s=1,t] I t-s * Ws
Σ[s=1,t] I t-s * Wsの部分が畳み込み積分で
Ws ∝ serial interval
- 167 ：卵の名無しさん：2020/05/18(月) 07:52:43.74 ID:Tv0fQD3N.net
- # https://www.nejm.org/doi/full/10.1056/NEJMoa2001316
#--- incubation period ---
# from Li et al NEJM 2020
# lognormal mean = 5.2
ln.par1 = 1.434065
ln.par2 = 0.6612
x=rlnorm(1e6,ln.par1,ln.par2)
BEST::plotPost(x)
curve(dlnorm(x,ln.par1,ln.par2),add=T)
quantile(x,c(0.025,0.5,0.975))
"
c.onset = c.infected + c.incubation
d.onset = d.infected + d.incubation
c.onset - d.onset = onset.delay
c.infected - d.infected + c.incubation - d.incubation = delay
c.infected - d.infected = onset.delay + d.incubation - c.incubation
"
onset.delay = 2
d.incubation = rlnorm(1e6,ln.par1,ln.par2)
c.incubation = rlnorm(1e6,ln.par1,ln.par2)
infection.delay = onset.delay + d.incubation - c.incubation
BEST::plotPost(infection.delay,compVal = 0)
- 168 ：卵の名無しさん：2020/05/18(月) 13:03:08.05 ID:t88GB6TH.net
- stancode=
"
data{
real onset_delay;
real ln_par1;
real ln_par2;
}
parameters{
real <lower=0> d_incubation;
real <lower=0> c_incubation;
}
transformed parameters{
real infection_delay = onset_delay + d_incubation - c_incubation;
}
model{
d_incubation ~ lognormal(ln_par1,ln_par2);
c_incubation ~ lognormal(ln_par1,ln_par2);
}
"
model=stan_model(model_code = stancode)
fn.stan <- function(delay){
dataList=list(onset_delay=delay,ln_par1=ln.par1,ln_par2=ln.par2)
fit=sampling(model,data=dataList)
ms=rstan::extract(fit)
mean(ms$infection_delay < 0)
}
fn.stan(2)
onset_delays=0:20
y=sapply(onset_delays,fn.stan)
plot(onset_delays,y, ylab='Pr[ Infected Later ])',axes=F) ; axis(1)
- 169 ：卵の名無しさん：2020/05/18(月) 22:00:40.40 ID:mZhU0UjE.net
- 今日も朝から人生が暇しかない
Ｆラン統計ジジイ あらため コロナ薬屋 が
発狂して連投だぁぁぁW
ちゃんとコテハン付けないとダメでちよぉぉぉ
今思うのは反日野党が政権をとってたら
日本はとっくに壊滅していたであろうこと
だって反日なんだもの
by 都内Ｓラン女子高生
- 170 ：卵の名無しさん：2020/05/19(火) 07:38:54.76 ID:YFl3mfOu.net
- runJags=run.jags('TEMPmodel.txt',monitor=c('p1','p2','diff'),
data=dataList,n.chains=4,sample=10000,burnin=4000)
coda::gelman.plot(runJags)
codaSamples = as.mcmc.list(runJags)
- 171 ：卵の名無しさん：2020/05/19(火) 07:59:02.34 ID:Y2l7bcjw.net
- 今日も朝から人生が暇しかない
Ｆラン統計ジジイ あらため コロナ薬屋 が
発狂して連投だぁぁぁW
相変わらず臨床経験や社会経験がゼロなのが
丸わかりのレスですねぇぇぇぇぇ
ちゃんとコテハン付けないとダメでちよぉぉぉ
今思うのは反日野党が政権をとってたら
日本はとっくに壊滅していたであろうこと
だって反日なんだもの
by 都内Ｓラン女子高生
- 172 ：卵の名無しさん：2020/05/21(木) 12:57:01.59 ID:TgNVM+u3.net
- 今日のくるくるぱーのIDだああぁぁあ
id:3+Rla4zY
くるくるぱーが反日野党に肩入れして
自滅するのは自業自得だけど
善良な日本人は巻き込まれたくないですねえ
by 都内Ｓラン女子高生
- 173 ：卵の名無しさん：2020/05/21(木) 17:33:16.75 ID:TgNVM+u3.net
- 今日も朝から人生が暇しかない
Ｆラン統計ジジイ あらため コロナ薬屋 が
発狂して連投だぁぁぁW
相変わらず臨床経験や社会経験がゼロなのが
丸わかりのレスですねぇぇぇぇぇ
今思うのは反日野党が政権をとってたら
日本はとっくに壊滅していたであろうこと
だって反日なんだもの
by 都内Ｓラン女子高生
- 174 ：卵の名無しさん：2020/05/22(金) 17:43:59.44 ID:HsXVoBS6.net
- すいません、すごい基本的なことなのですが教えてください。
現在いる実験室で、様々細胞と遺伝子変異させた細胞に抗がん剤をかけて死亡率をみてます。
抗がん剤の量も比較してるのです
1: 細胞aとbに対して同量の抗がん剤を使用し死亡率を見る場合
2: 同種の細胞に違う濃度の抗がん剤をかけて死亡率を比較する場合。
両者とも何検体かとって平均の比較をする場合、対応するt検定ですか？、それともf値をみた上で対応しないt検定になるのですか？
- 175 ：卵の名無しさん：2020/05/22(金) 19:39:38.80 ID:lAN4YhgA.net
- >>174
比率の検定ならΧ（カイ）二乗検定。
- 176 ：卵の名無しさん：2020/05/23(土) 07:42:32.91 ID:Myh/vXaP.net
- 今日も朝から人生が暇しかない
Ｆラン統計ジジイ あらため コロナ薬屋 が
発狂して連投だぁぁぁW
身の程を弁えない謎の上から目線で政治談批評
老害アルアルですねぇぇぇぇぇ
今思うのは反日野党が政権をとってたら
日本はとっくに壊滅していたであろうこと
だって反日なんだもの
by 都内Ｓラン女子高生
- 177 ：卵の名無しさん：2020/05/23(土) 07:57:58.92 ID:PaqrAdk5.net
- 日本衰退を願う勢力は安倍一択。
観光立国という亡国政策で衰退途上国。
観光でしか生きていけない国になったときに中国からの出国禁止すれば日本は中国に平伏すしかなくなる。
- 178 ：卵の名無しさん：2020/05/23(土) 07:58:44.46 ID:PaqrAdk5.net
- >>174
対数変換して検定した方がいい場合もある。
- 179 ：卵の名無しさん：2020/05/23(土) 22:14:18.97 ID:/303gUSa.net
- テレビを見てるんですが
やっぱりマスコミはクズですね
問題のあった新聞社の社長は直ちに謝罪して
新聞社を解体するべきですね
反日野党は日本の足を引っ張ることしかしてませんね
存在意義がないとしか言えません
だって反日なんだもの
by 都内Ｓラン女子高生
- 180 ：卵の名無しさん：2020/05/24(日) 15:16:08.53 ID:3OldU+Ns.net
- >>178
対数変換した場合、検定方法は何使えばいいですか？
- 181 ：卵の名無しさん：2020/05/24(日) 18:15:55.34 ID:bhuCUyQ5.net
- >>180
平均の検定ならt検定
- 182 ：卵の名無しさん：2020/05/25(月) 12:16:59.94 ID:m/x5AhgL.net
- ある人物Dが新型コロナ肺炎に罹患したとする。行動調査によって発症前にキャバクラに行っており
接客したキャバ嬢Cが人物D発症の２日後に発症していたことがわかった。
キャバ嬢Cは人物Dから移されたと主張して１億円の賠償を求めている。
潜伏期間には幅がありキャバ嬢Cから移された可能性もあると主張してその確率を計算して賠償金を値切りたい。
いくら値切れるか計算せよ。
#--- incubation period ---
# from Li et al NEJM 2020
# lognormal mean = 5.2
ln_par1 = 1.434065
ln_par2 = 0.6612
Gt <- function(delay){
C=rlnorm(1e6,ln_par1,ln_par2)
D=rlnorm(1e6,ln_par1,ln_par2)
mean(C-D > delay)
}
Gt(2)
library(polspline)
delay=2
c=rlnorm(1e5,ln_par1,ln_par2)
d=rlnorm(1e5,ln_par1,ln_par2)
hist(c-d, freq=F,col='white',breaks=100,ylim=c(0,0.11),
xlim=c(-30,30),ann=F,axes=F) ; axis(1)
fit=logspline(c-d)
curve(dlogspline(x,fit),-30,30,ann=F,bty='n',add=T)
segments(delay,0,delay,dlogspline(delay,fit),pch=19,col=2)
curve(dlogspline(x,fit),delay,30,add=T,type='h',col=2)
1-plogspline(delay,fit)
fn <- function(delay) 1- plogspline(delay,fit)
curve(fn(x),0,14, bty='n' ,xlab='Delay', ylab='Probability')
- 183 ：卵の名無しさん：2020/05/25(月) 12:21:54.11 ID:FKqwK9d6.net
- >>181
t検定は二つのことを比べることしかできないですよね？
- 184 ：卵の名無しさん：2020/05/25(月) 16:55:34.95 ID:m/x5AhgL.net
- >>183
多重検定やりたいならTukeyとか色々ある。
- 185 ：卵の名無しさん：2020/05/25(月) 17:16:47.59 ID:FKqwK9d6.net
- >>184
図の点線が対照のmockで実線が遺伝子変化させた細胞、
横軸が加えた薬剤の濃度です。
ルシフェラーゼの比を見ています
同じ濃度で、モックと遺伝子変化させた群を比較させたい時は統計は何使えばいいですか？
https://i.imgur.com/VXfJ8TD.jpg
- 186 ：卵の名無しさん：2020/05/25(月) 17:49:48.35 ID:m/x5AhgL.net
- 好きなの使えば。
厳しいbonferroni補正したpairwise.t.testとか
- 187 ：卵の名無しさん：2020/05/25(月) 18:00:37.69 ID:m/x5AhgL.net
- 回帰係数の比較ならコクランアーミテッジだったかな。
Rだとprop.trend.testでできたはず！
- 188 ：卵の名無しさん：2020/05/25(月) 19:07:25.19 ID:I4VYZbvs.net
- 皆さんありがとうございます
すいませんjmpかエクセル統計で計算できるのだと嬉しいです。
- 189 ：卵の名無しさん：2020/05/25(月) 20:15:07.84 ID:ZmKINRjj.net
- やったね たえちゃん
- 190 ：卵の名無しさん：2020/05/25(月) 21:59:58.02 ID:m/x5AhgL.net
- Rだと
kruskal.testで多重比較で有意差確認してから
平均比較ならpairwise.t.test
比率比較ならpairwise.prop.test
補正法は指定できる。
エクセルのマクロは売られている。
使ったこともないけど。
https://bellcurve.jp/ex/function/kruskal.html
- 191 ：卵の名無しさん：2020/05/26(火) 07:48:20.04 ID:ZksetU7j.net
- 今日も朝から人生が暇しかない
Ｆラン統計ジジイ あらため コロナ薬屋 が
発狂して連投だぁぁぁW
それにしても反日野党や反日マスコミは日本の足を引っ張ることしかしてませんね
存在意義がないとしか言えません
だって反日なんだもの
by 都内Ｓラン女子高生
- 192 ：卵の名無しさん：2020/05/26(火) 12:50:37.25 ID:5lPxAo99.net
- "
>勤務医が院長のおれの悪口を妻に話してたそうだ
>腹が立ったからクビにする
を登場人物としてみた。
(問題)
新型コロナに勤務医が罹患。
勤務医が発症した翌日に院長が発症、その2日後に妻が発症した。
(1)妻が感染源である(最初に感染していた)確率を求めよ。
(2)感染順が妻→勤務医→院長の順である確率を求めよ。
#--- incubation period ---
# from Li et al NEJM 2020
# lognormal mean = 5.21 sd=3.86
# ln.par1 = 1.434065
# ln.par2 = 0.6612
rm(list=ls())
ln_par1 = 1.434065
ln_par2 = 0.6612
Aincub=rlnorm(1e6,ln_par1,ln_par2) # 勤務医
Bincub=rlnorm(1e6,ln_par1,ln_par2) # 院長
Cincub=rlnorm(1e6,ln_par1,ln_par2) # 妻
"
Ainfected=Aonset-Aincub
Binfected=Bonset-Bincub=Aonset+1-Bincub
Cinfected=Conset-Cincub=Bonset+2-Cincub=Aonset+3-Cincub
"
Aonset=0
Ainfected=Aonset-Aincub
Binfected=Aonset+1-Bincub
Cinfected=Aonset+3-Cincub
ABCinfected = as.data.frame(cbind(Ainfected,Binfected,Cinfected))
boxplot(ABCinfected)
- 193 ：卵の名無しさん：2020/05/26(火) 12:50:43.88 ID:5lPxAo99.net
- fn1 <- function(x) min(x)==x['Cinfected']
mean(apply(ABCinfected,1,fn1))
fn2 <- function(x){
x['Cinfected'] < x['Ainfected']
& x['Ainfected'] < x['Binfected']
}
mean(apply(ABCinfected,1,fn2))
- 194 ：卵の名無しさん：2020/05/26(火) 23:37:54.50 ID:+nMfC17p.net
- まず全裸になり
( ： )
（ ゜∀゜）ﾉ彡
<( )
ﾉωヽ
自分の尻を両手でバンバン叩きながら白目をむき
从
Дﾟ ） て
( ヾ) )ヾ て
< <
人__人__人__人__人__人__人__人__人__人__人
Σ て
Σ びっくりするほどユートピア！ て人__人_
Σ びっくりするほどユートピア！ て
⌒Ｙ⌒Ｙ⌒Ｙ) て
Ｙ⌒Ｙ⌒Ｙ⌒Ｙ⌒Ｙ⌒Ｙ⌒Ｙ⌒Ｙ⌒Ｙ⌒Ｙ⌒Ｙ⌒
＿_＿＿＿＿＿
|＿＿ ヽ(゜∀゜）ﾉ
|\＿〃´￣￣ ヽ..ヘ( )ミ
| |\,.-～´￣￣ ω > （∀゜ ）ﾉ
\|∫\ _,. - ､_,. - ､ \ ( ﾍ)
\ \＿＿＿＿＿＿ _\<
\ ||￣￣￣￣￣￣￣ |
\||＿＿＿＿＿＿＿ |
これを10分程続けると妙な脱力感に襲われ、解脱気分に浸れる。
- 195 ：卵の名無しさん：2020/05/27(水) 04:43:18.93 ID:uX2xSEBS.net
- x1,x2,...,xnの順に発症
その間隔はt1,t2,..,tn-1としたときにxiが感染源であった確率を算出するプログラムを作れ。
- 196 ：卵の名無しさん：2020/05/27(水) 05:15:26.64 ID:uX2xSEBS.net
- 安倍は新コロナは中国発と認めたから春節ウェルカムは外患誘致。
早くアメリカと足並み揃えて中国に損害賠償請求して国民に赦しを乞うべき。
観光立国という亡国政策で日本は衰退途上国。
次の世代の日本人は実習生として中国に出稼ぎがで立場が逆転。
日本人なら過労死自殺しても経営者殺害はしないから酷使される。
観光でしか生きていけない国になったときに中国からの訪日禁止すれば日本は中国に平伏すしかなくなる。
尖閣どころか沖縄を寄越せとか言われても差し出すことになりそうだな。
- 197 ：卵の名無しさん：2020/05/27(水) 06:40:49.96 ID:uX2xSEBS.net
- >>195
モデルは簡単なのでStanやJagsを使わずに完成
WhoInfectedFirst <- function( # 発症間隔から感染源の確率を計算
t=c(1,2), # 発症間隔
k=1e5 # 乱数発生数
){
ln_par1 = 1.434065 # 潜伏期間対数正規分布パラメータ
ln_par2 = 0.6612
n=length(t)+1 # 発症人数
# 潜伏期間
x_incub = matrix(rep(NA,n*k),ncol=n)
for(i in 1:n) x_incub[,i] = rlnorm(k,ln_par1,ln_par2)
# 感染日（一人目の発症日：０）
x_infected = matrix(rep(NA,n*k),ncol=n)
x_infected[,1]= -x_incub[,1]
for(i in 2:n) x_infected[,i] = sum(t[1:(i-1)]) - x_incub[,i]
# i番目の発症者が感染源の確率
fi <- function(i){
mean(apply(x_infected,1,function(x) which.min(x)==i))
}
data.frame(p=round(sapply(1:n,fi),2))
}
WhoInfectedFirst(c(1,2))
WhoInfectedFirst(c(1,0,1,0,0)) # 翌日2人発症 翌々日3人発症
- 198 ：卵の名無しさん：2020/05/27(水) 07:57:44.10 ID:KMcVqry5.net
- やったね たえちゃん
Ｆラン統計ジジイ あらため コロナ薬屋は
今日もクルクルパーだよ
- 199 ：卵の名無しさん：2020/05/27(水) 17:40:48.61 ID:uX2xSEBS.net
- >>197
忘れないようにStanで書いてみた。
################# MCMC by stan ############
stancode='
data {
int<lower=0> n;
vector[n-1] t;
real ln_par1;
real ln_par2;
}
parameters {
real<lower=0> x_incub[n];
}
model {
target += lognormal_lpdf(x_incub|ln_par1,ln_par2);
}
generated quantities{
real x_infected[n];
x_infected[1]= -x_incub[1];
for(i in 2:n){
x_infected[i] = sum(t[1:(i-1)]) - x_incub[i];
}
}
'
- 200 ：卵の名無しさん：2020/05/27(水) 17:40:55.22 ID:uX2xSEBS.net
- stanmodel=stan_model(model_code = stancode)
saveRDS(stanmodel,'cavaret.rds')
t=c(1,0,1,0,0)
n=length(t)+1
ln_par1 = 1.434065
ln_par2 = 0.6612
data=list(t=t,n=n,ln_par1=ln_par1,ln_par2=ln_par2)
fit=sampling(stanmodel,data=data,iter=1e5,thin=1,chains=4,
control=list(adapt_delta=0.95))
pars=c('x_infected')
print(fit,pars=pars)
stan_dens(fit,pars=pars,separate_chains = T)
ms=rstan::extract(fit)
head(ms$x_infected)
x_infected=ms$x_infected
fi <- function(i){ # i番目の発症者が感染源の確率
mean(apply(x_infected,1,function(x) which.min(x)==i))
}
data.frame(p=round(sapply(1:n,fi),2))
- 201 ：卵の名無しさん：2020/05/27(水) 17:41:13.27 ID:uX2xSEBS.net
- > data.frame(p=round(sapply(1:n,fi),2))
p
1 0.26
2 0.18
3 0.18
4 0.13
5 0.13
6 0.13
- 202 ：卵の名無しさん：2020/05/27(水) 17:41:42.49 ID:uX2xSEBS.net
- JAGSだと
library(rjags)
cat('
model{
ln_par1 = 1.434065 # 潜伏期間対数正規分布パラメータ
ln_par2 = 0.6612
# 潜伏期間
for(i in 1:n){
x_incub[i] ~ dlnorm(ln_par1,ln_par2)
}
# 感染日（一人目の発症日：０）
x_infected[1]= -x_incub[1]
for(i in 2:n){
x_infected[i] = sum(t[1:(i-1)]) - x_incub[i]
}
}',file='tmp.txt')
t=c(1,0,1,0,0)
n=length(t)+1
dataList=list(t=t,n=n)
jagsModel=jags.model('tmp.txt',data=dataList, n.chains=4,n.adapt=5e3)
update(jagsModel)
codaSamples=coda.samples(jagsModel,n.iter=1e5,thin=1,var=c('x_infected'))
x_infected=as.matrix(codaSamples)
fi <- function(i){ # i番目の発症者が感染源の確率
mean(apply(x_infected,1,function(x) which.min(x)==i))
}
data.frame(p=round(sapply(1:n,fi),2))
- 203 ：卵の名無しさん：2020/05/27(水) 17:44:39.93 ID:uX2xSEBS.net
- > data.frame(p=round(sapply(1:n,fi),2))
p
1 0.20
2 0.17
3 0.17
4 0.15
5 0.15
6 0.15
> WhoInfectedFirst(c(1,0,1,0,0)) # 翌日2人発症 翌々日3人発症
p
1 0.26
2 0.18
3 0.18
4 0.13
5 0.13
6 0.13
多数決でStanの値を採用しよう。
- 204 ：卵の名無しさん：2020/05/27(水) 17:46:07.20 ID:uX2xSEBS.net
- 発症した順が必ずしも感染順ではないので、
職員が発症しても責めたりせずに悪いのは春節ウェルカムした安倍のせいと心をまとめよう。
- 205 ：卵の名無しさん：2020/05/28(木) 17:57:18.03 ID:CdrIPYow.net
- # 武漢大量検査
fit=PCRs5a(6500000,218,SEN=0.6,SD1=0.2,SPC=0.999,SD2=0.001,
iter=1e5,warmup=2e4,thin=10,chain=1)$fit
print(fit,digits=7,prob=c(0.025,0.5,0.975),pars=c('spc','sen','prev','p'))
stan_dens(fit, separate_chains = T)
stan_plot(fit)
stan_ac(fit)
stan_trace(fit)
ms=rstan::extract(fit)
density2D(ms$sen,ms$spc)
summary(ms$p*1e7) ; MODE(ms$p*1e7)[1] ; hdi(ms$p*1e7)[1:2]
- 206 ：卵の名無しさん：2020/05/28(木) 17:58:09.17 ID:CdrIPYow.net
- > print(fit,digits=7,prob=c(0.025,0.5,0.975),pars=c('spc','sen','prev','p'))
Inference for Stan model: 451fbbe263cf71a87ea567fe0852d645.
1 chains, each with iter=100000; warmup=20000; thin=10;
post-warmup draws per chain=8000, total post-warmup draws=8000.
mean se_mean sd 2.5% 50% 97.5% n_eff Rhat
spc 0.9999755 0.0000001 0.0000074 0.9999644 0.9999742 0.9999920 8082 1.00007
sen 0.4561154 0.0025058 0.2177783 0.0776906 0.4454970 0.8752057 7553 1.00002
prev 0.0000336 0.0000009 0.0000767 0.0000007 0.0000186 0.0001547 8059 1.00005
p 0.0000341 0.0000000 0.0000023 0.0000296 0.0000341 0.0000388 7840 0.99988
> summary(ms$p*1e7) ; MODE(ms$p*1e7)[1] ; hdi(ms$p*1e7)[1:2]
Min. 1st Qu. Median Mean 3rd Qu. Max.
258 325 341 341 356 433
x
342.03
lower upper
296.90 388.22
- 207 ：卵の名無しさん：2020/05/28(木) 18:02:34.43 ID:CdrIPYow.net
- 検査は900万人らしいから最終的な陽性者数の95％信頼区間は267-349と出てきたな。
> summary(ms$p*9e6) ; MODE(ms$p*9e6)[1] ; hdi(ms$p*9e6)[1:2]
Min. 1st Qu. Median Mean 3rd Qu. Max.
233 293 307 307 320 390
x
307.83
lower upper
267.21 349.40
- 208 ：卵の名無しさん：2020/05/30(土) 10:34:55.69 ID:eqmXwXuQ.net
- >>182
乱数発生の代わりに数値積分で計算
gt = 2 # generation time (infection interval)
mu = 1.434065
sg = 0.6612
# 差のpdfの公式
# pdf1 <- function(x) dlnorm(x,mu,sg)
# pdf2 <- function(y) dlnorm(y,mu,sg)
# ∫[-∞,∞] pdf1(x+y)*pdf2(y) dy
f <- function(x,y) dlnorm(x+y,mu,sg)*dlnorm(y,mu,sg)
vf=Vectorize(f,vectorize.args = 'y')
pdf <- function(x) integrate(function(y) vf(x,y),
rel.tol = 1e-14,-Inf,Inf)$value
pdf=Vectorize(pdf)
curve(pdf(x),-30,30)
cdf <- function(x) integrate(pdf,x,Inf,rel.tol = 1e-14)$value
cdf(gt)
# 乱数発生との比較
c=rlnorm(1e6,mu,sg)
d=rlnorm(1e6,mu,sg)
hist(c-d, freq=F,breaks=200,ylim=c(0,0.11),col='skyblue',
xlim=c(-30,30),ann=F,axes=F) ; axis(1)
curve(pdf(x),add=T)
mean(c-d > 2)
- 209 ：卵の名無しさん：2020/05/30(土) 10:50:01.68 ID:eqmXwXuQ.net
- gt = 2 # generation time (infection interval)
mu = 1.434065
sg = 0.6612
# 差のpdfの公式
# pdf1 <- function(x) dlnorm(x,mu,sg)
# pdf2 <- function(y) dlnorm(y,mu,sg)
# ∫[-∞,∞] pdf1(x+y)*pdf2(y) dy
f <- function(x,y) dlnorm(x+y,mu,sg)*dlnorm(y,mu,sg)
vf=Vectorize(f,vectorize.args = 'y')
pdf <- function(x) integrate(function(y) vf(x,y),
rel.tol = 1e-14,-Inf,Inf)$value
pdf=Vectorize(pdf)
curve(pdf(x),-30,30)
cdf <- function(x) integrate(pdf,x,Inf,rel.tol = 1e-14)$value
cdf(gt)
# 発症間隔を変化させてグラフ化
gts=seq(0,14,by=0.5)
plot(gts,sapply(gts,cdf),type='l',xlab='発症間隔(日)',
ylab='逆順感染確率',las=1)
#
c=rlnorm(1e6,mu,sg)
d=rlnorm(1e6,mu,sg)
hist(c-d, freq=F,breaks=200,ylim=c(0,0.11),col='skyblue',
xlim=c(-30,30),ann=F,axes=F) ; axis(1)
curve(pdf(x),add=T)
mean(c-d > 2)
- 210 ：卵の名無しさん：2020/06/01(月) 17:57:09.41 ID:OULqunq+.net
- 今日のID:nu4+jSn4
- 211 ：卵の名無しさん：2020/06/01(月) 20:31:14.88 ID:nu4+jSn4.net
- >>210
ゴルゴ問題の答が出せるように勉強したのか？
- 212 ：卵の名無しさん：2020/06/26(金) 17:08:21.56 ID:vSOmQn7C.net
- "
拳銃を一発撃ったときに、
狙った相手を撃ち殺す確率は、
Ａは 1/3、Ｂは 1/2（50%）、Ｃは 1/1（100%）とします。
なお、この確率は、全員が知っているものとします。
拳銃を撃つ順番は、Ａ、Ｂ、Ｃの順番で、
以降は最後の一人が生き残るまでこの順番を繰り返すものとします。
Ａ、Ｂ、Ｃは拳銃を撃つときに誰を狙っても良いこととします。
ただし、一発で二人を狙うことはできません。
Ａ、Ｂ、Ｃの生き残る確率を求めなさい。
"
- 213 ：卵の名無しさん：2020/06/26(金) 17:08:35.96 ID:vSOmQn7C.net
- Duel <- function(
pa=1/3,
pb=1/2,
pc=1/1,
k=1e5
){
qa=1-pa # probability of dead ; alive
qb=1-pb
qc=1-pc
# simulation of 3 alive to 2 alive
# each shoots at the superior sniper
# Hence,A shoots C -> B shoots C -> C shoots B
f3 <- function(){
a=b=c=1
bshoot=FALSE # is next shooter B?
while(a+b+c==3){
c=rbinom(1,1,qa) # A shoots C
if(c==0) bshoot=TRUE # if A kills C, next shooter is B
c=rbinom(1,1,qb) # B shoots C
if(c==1) b=rbinom(1,1,qc) # if C alive, C shoots B
}
return(list(abc=c(a,b,c),bshoot=bshoot))
}
- 214 ：卵の名無しさん：2020/06/26(金) 17:08:51.89 ID:vSOmQn7C.net
- # simulation of 2 alive to 1 alive
# Notice who has the right to shoot
f2 <- function(){
re3=f3() # simulation to 2 suvivors
a=re3$abc[1] # 1 1
b=re3$abc[2] # 1 0
c=re3$abc[3] # 0 1
while(a+b+c==2){ # while 2 alive
if(c==0){ # when C dead
if(re3$bshoot){ # when A killed C, B can shoot A
a=rbinom(1,1,qb)
if(a==0) return(c(a,b,c))
}
else{ # when B killed C
b=rbinom(1,1,qa) # A can shoot B
if(b==0) return(c(a,b,c))
}
}
else{ # when C alive (C killed B, C shot at superior B)
c=rbinom(1,1,qa) # A can shoot C
if(c==0) return(c(a,b,c)) # if A killed C
a=rbinom(1,1,qc) # when A missed C,C can shoot A
if(a==0) return(c(a,b,c))
}
}
}
re2=replicate(k,f2())
return(apply(re2,1,mean))
}
Duel(1/3,1/2,1)
Duel(0.3,0.4,0.5)
- 215 ：卵の名無しさん：2020/06/28(日) 08:05:03.04 ID:uvFNh4/q.net
- # Debugged verision
rm(list=ls()) # clear workspace
p=c(1/3,1/2,1) # killing probablity
q=1-p # survival probability
# Three survivors to two survivors
f32 <- function(init=1){
abc=c(1,1,1) # dead or alive
n=length(abc) # initial survivors 3
s=sum(abc) # current survivors
sh=init # initial shooter
while(s==n){ # while 3 survivors
tmp=p
tmp[sh]=0 # set shooter for 0 probability (could be negative)
target=which.max(tmp) # target : the superior sniper
abc[target]=rbinom(1,1,q[sh]) # target dead(0) or alive(1)
# next shooter : if next sniper alive next index, otherwise survived sniper
sh=ifelse(abc[sh%%n+1]==1,sh%%n+1,(1:n)[-c(sh,target)])
s=sum(abc) # how many surviors left
}
list(abc=abc,sh=sh)
}
# demo
f32()
apply(replicate(1e5,f32()$abc),1,mean)
- 216 ：卵の名無しさん：2020/06/28(日) 08:05:24.51 ID:uvFNh4/q.net
- # Two survivors to one survivor
f21 <- function(){ # Three survivors to 2 survivor
re32=f32() # two surviors & next shooter
abc=re32$abc # two survivor
sh=re32$sh # shooter
sv=which(abc==1) # index of two suvivor
target=sv[sv!=sh] # index of target
s=sum(abc)
while(s==2){ # while two survivors, mutual shooting
abc[target]=rbinom(1,1,q[sh]) # target dead(0) or alive(1)
s=sum(abc) # how many surviors left
tmp=target # exhange shooter for target
target=sh
sh=tmp
}
abc
}
f21()
apply(replicate(1e5,f21()),1,mean)
- 217 ：卵の名無しさん：2020/07/05(日) 18:42:40.47 ID:T+5emEsp.net
- rm(list=ls())
mu=0.58 # mean of reproductive number
size=0.45 # dispersion parameter
(prob = size/(size+mu)) # its probability
Rt=rnbinom(1e5,size=size,mu=mu) # random numbers of negative binomial distribution
hist(Rt,breaks = 'scott',freq=F,ann=F) # show its histgram
sim <- function(n=10){ # simulation
infectee=0 # initial value
while(infectee!=n){ # while infectee is unequal to n
infector=sample(n,1) # prior discrete uniform distribution of infector number
infectee=sum(sample(Rt,infector)) # total number of infectee
}
return(infector) # when n infectee, return infector number
}
spreader=replicate(1e5,sim()) # simulation & calculation
hist(spreader,freq=F,ylab='',axes=F,breaks='scott') ; axis(1)
HDInterval::hdi(spreader)[1:2] # 95% credibility interval
BEST::plotPost(spreader) # graph with 95%CI & mean
quantile(spreader,c(0.025,0.5,0.975))
summary(spreader)
sum(spreader==1)/length(spreader) # the probability of single super-spreader
- 218 ：卵の名無しさん：2020/07/05(日) 18:43:03.70 ID:T+5emEsp.net
- > quantile(spreader,c(0.025,0.5,0.975))
2.5% 50% 97.5%
4 8 10
> summary(spreader)
Min. 1st Qu. Median Mean 3rd Qu. Max.
1.000 7.000 8.000 8.039 10.000 10.000
> sum(spreader==1)/length(spreader) # the probability of single super-spreader
[1] 0.0016
- 219 ：卵の名無しさん：2020/07/05(日) 20:04:10.52 ID:T+5emEsp.net
- デバッグ版
R=0.58 # mean of reproductive number
k=0.45 # dispersion parameter
(prob = k/(k+R)) # its probability
Rt=rnbinom(1e5,k,mu=R) # random numbers of negative binomial distribution
hist(Rt,breaks = 'scott',freq=F,ann=F) # show its histgram
sim <- function(n=10){ # simulation
infected=0 # initial value
while(infected!=n){ # while infected is unequal to n
infector=sample(n,1) # prior discrete uniform distribution of infector number
infectee=sum(sample(Rt,infector)) # number of infectee
infected=infectee+infector # number of infected
}
return(infector) # when n infected, return infector number
}
spreader=replicate(1e5,sim()) # simulation & calculation
hist(spreader,freq=F,ylab='',axes=F,breaks='scott') ; axis(1)
HDInterval::hdi(spreader)[1:2] # 95% credibility interval
BEST::plotPost(spreader) # graph with 95%CI & mean
summary(spreader)
sum(spreader==1)/length(spreader) # the probability of single super-spreader
- 220 ：卵の名無しさん：2020/07/05(日) 20:04:14.42 ID:T+5emEsp.net
- まあ、パーティー中に感染させる人数に再生産数を流用していいかは疑問ではあるが、
実効結果は
https://i.imgur.com/ApDSW1J.png
95% CI
> HDInterval::hdi(spreader)[1:2] # 95% credibility interval
lower upper
3 9
中央値、平均値、四分位値
> summary(spreader)
Min. 1st Qu. Median Mean 3rd Qu. Max.
1.000 6.000 7.000 6.844 8.000 10.000
1人のスーパースプレッダーの確率
> sum(spreader==1)/length(spreader) # the probability of single super-spreader
[1] 0.00104
- 221 ：卵の名無しさん：2020/07/07(火) 16:20:43.89 ID:vMNjk2cp.net
- # n種類のガチャからm種類を集めるまでの期待値
collector <- function(n=100,m=80,print=TRUE){
library(gmp)
x=(n-m+1):n
x=as.bigq(x)
y=sum(n/x)
if(print) print(y)
return(asNumeric(y))
}
collector(5,4)
collector(100,80)
- 222 ：卵の名無しさん：2020/07/09(木) 21:20:16.36 ID:N2v3KqqG.net
- library(randtests)
DEL <- permut(1:4,3,function(x) print(paste0(x,collapse = ''),quote=F))
- 223 ：卵の名無しさん：2020/07/10(金) 23:07:16.57 ID:1U5aVO+X.net
- まじで、エクセル、できればLibreOffice で出来る
医療統計の本を出して欲しい、この先生なら出来る。
- 224 ：卵の名無しさん：2020/07/10(金) 23:13:07.04 ID:Fs8YjPFQ.net
- 自分で、この先生とか言うなよ
恥ずかしい奴だな
- 225 ：卵の名無しさん：2020/07/13(月) 21:47:22.50 ID:brmpmlYq.net
- リアル223からです、数式を記載されている方とは別人です。
- 226 ：卵の名無しさん：2020/07/14(火) 17:33:19.37 ID:V6RZOqwd.net
- >>223
エクセルのソルバーなんて使えたものじゃないね。
エクセルはグラフも貧弱だし。
医療統計の本なら朝倉から丹波某が沢山書いた本が沢山でているよ。
ベイズはWINBUGの話で古いけど。
- 227 ：卵の名無しさん：2020/07/18(土) 23:13:50.16 ID:4YqdABa6.net
- "
容量8Lの袋と容量5Lの袋を使って池の水を丁度4L集めたい。
袋に目盛りはついていません。
袋から袋への移し替えは全量で行います。
池からとる水の量や池に捨てる水の量には制限はありません。
最初に片方に満たした作業を１回目として
丁度4Lを集めるのに最低何回の移動が必要か？
"
- 228 ：卵の名無しさん：2020/07/18(土) 23:14:05.95 ID:4YqdABa6.net
- abura <- function(
a7=8,
b3=5,
c5=4){
# starting from the bigger pitcher
movea7 <- function(xy){ # start from c(a7,0)
x=xy[1] ; y=xy[2]
# x==a7
if(x==a7) re=c(a7-(b3-y),b3)
# x==0
if(x==0) re=c(a7,y)
# y==b3
if(y==b3) re=c(x,0)
# y==0
if(y==0 & x!=a7){
if(x>=b3) re=c(x-b3,b3)
else re=c(0,x)
}
return(re)
}
STATUS=status=c(a7,0)
i=1
while(!identical(status,c(0,0))){#
i=i+1
status=movea7(status)
STATUS=rbind(STATUS,status)
}
rownames(STATUS)=1:nrow(STATUS)
colnames(STATUS)=c(paste0(a7,'L'),paste0(b3,'L'))
(Bigger=STATUS)
- 229 ：卵の名無しさん：2020/07/18(土) 23:14:23.30 ID:4YqdABa6.net
- # starting from the smaller pitcher
moveb3 <- function(xy){ # start from c(0,b3)
x=xy[1] ; y=xy[2]
if(y==b3){
if(x<=(a7-b3)) re=c(x+b3,0)
else re=c(a7, b3-(a7-x))
}
if(y==0) re=c(x,b3)
if(x==a7) re=c(0,y)
if(x==0) re=c(y,0)
return(re)
}
STATUS=status=c(0,b3)
i=1
while(!identical(status,c(0,0))){ # stop at c(5,0) for solution
i=i+1
status=moveb3(status)
STATUS=rbind(STATUS,status)
}
rownames(STATUS)=1:nrow(STATUS)
colnames(STATUS)=c(paste0(a7,'L'),paste0(b3,'L'))
(Smaller=STATUS)
min_Bigger=min(which(apply(Bigger,1,function(x) c5 %in% x)))
min_Smaller=min(which(apply(Smaller,1,function(x) c5 %in% x)))
list(Bigger=Bigger,Smaller=Smaller,
min_Bigger=min_Bigger,min_Smaller=min_Smaller)
}
- 230 ：卵の名無しさん：2020/07/26(日) 10:46:35.07 ID:sQTZfadO.net
- '%&%' <- function(x,y) paste0(x,y)
> '√2 =' %&% sqrt(2) %&% ', π = ' %&% pi
[1] "√2 =1.4142135623731, π = 3.14159265358979"
random <- function(n,char=c(LETTERS,letters,0:9)){
re=paste0(sample(char,n),collapse='')
cat(re,'\n')
invisible(re)
}
> random()
ujbDVK8pEtsNFTAo1BlkZP72cXRHrQvI0WdMmGyxOn3qz5awLhCg6S9eYJU4if
- 231 ：卵の名無しさん：2020/08/04(火) 21:12:57.26 ID:SiNiiw+2.net
- onCone <- function(p,q){ # (p,q) any figure on cone
x=sqrt(1-q^2)/pi*sin(pi*p/sqrt(1-q^2))
y=q
z=sqrt(1-q^2)/pi - sqrt(1-q^2)/pi *cos(pi*p/sqrt(1-q^2)) + q
list(x=x,y=y,z=z)
}
- 232 ：卵の名無しさん：2020/08/05(水) 16:36:41.12 ID:vZG88Q7w.net
- oncone <- function(p,q,α=pi/6){ # (p,q) 展開図上の座標、頂点の角度=２α
PQ=sqrt(p^2+q^2)
β=pi*sin(α)
θ=atan(p/q)
rdash=PQ*β/pi
γ=PQ*θ/rdash
x=rdash*sin(γ)
y=PQ/(tan(α)*sqrt(1+tan(α)^-2))
z=rdash*(1-cos(γ))
c(x=x,y=y,z=z)
}
https://i.imgur.com/qOojo9N.png
- 233 ：卵の名無しさん：2020/08/06(木) 07:04:19.79 ID:1+RzuAJZ.net
- numlockキーをbackspaceキーに変更するレジストリ（要再起動）
Windows Registry Editor Version 5.00
[HKEY_LOCAL_MACHINE\SYSTEM\CurrentControlSet\Control\Keyboard Layout]
"Scancode Map"=hex:00,00,00,00,00,00,00,00,02,00,00,00,0e,00,45,00,00,00,00,00
ソースは
https://www.souichi.club/technology/numlock-on/
- 234 ：卵の名無しさん：2020/08/09(日) 22:02:52.31 ID:nZS47/5e.net
- oncone <- function(p,q,α=40*pi/180){ # (p,q) 展開図上の座標、頂点の角度=２α
OA=sqrt(p^2+q^2)
θ=Arg(p+1i*q) # == atan2(p+1i*q)
β=pi*sin(α) # π*r=R*β ; r=R*sin(α)
rdash=OA*sin(α) # r'== OA*β/π, r'*π == OA*β : hemi-circle==arc(B-c)
δ=(pi/2-θ)/sin(α) # r'*δ == OA*(π/2-θ) ; δ == OA*(π/2-θ)/r'
x=rdash*sin(δ) # Re(E)
y=rdash/tan(α) # tan(α)=r'/y
z=rdash-rdash*cos(δ) # IM(E-B) ; rdash+rdash(cos(π-δ))
c(x=x,y=y,z=z)
}
https://i.imgur.com/Nq1xRAc.png
- 235 ：卵の名無しさん：2020/08/17(月) 22:48:07.01 ID:v42vnrN1.net
- #原子数=放射能/崩壊定数=放射能/(log2/半減期)=放射能*半減期*log2
#放射能が同じなら原子数∝半減期
cesium_now <- function(Date=NULL,RCs134=1,RCs137=1,TCs134 = 2.0652,TCs137 = 30.16171){
t=ifelse(is.null(Date),as.numeric((Sys.Date()-as.Date("2011/3/11"))/365.2425)
,as.numeric(as.Date(Date)-as.Date("2011/3/11"))/365.2425)
# mol ratio
cs <- RCs134*TCs134*(1/2)^(t/TCs134) + RCs134*TCs137*(1/2)^(t/TCs137)
cs0 <- RCs134*TCs134 + RCs134*TCs137
ratio=cs/cs0
# radioavtivity ratio
# decay constant = log(2)/half-life
# ratioactivity ∝ decay constant * mol
CS <- (1/2)^(t/TCs134) + (1/2)^(t/TCs137)
CS0 <- TCs134 + TCs137
Ratio=CS/CS0
list(year=t,mol_ratio=ratio,radioactivity_ratio=Ratio)
}
cesium_now()
cesium_now(NULL,1,1,2,30)
# 放射能比=1:1
TCs134 <- 2.0652 # 半減期(年）
TCs137 <- 30.16171
CS <- function(t) (1/2)^(t/TCs134) + (1/2)^(t/TCs137)
curve(CS(x),0,30)
uniroot(function(t,u0=1/2) CS(t)/CS(0)- u0, c(0,30))$root
- 236 ：卵の名無しさん：2020/08/22(土) 07:22:17.76 ID:d8M92xSQ.net
- # 問題「4/5より大きく5/6より小さい分数で、分母がいちばん小さい分数はなに？」
fn <- function(lo=4/5,up=5/6){
i=1
flg=FALSE
while(flg==FALSE){
for(j in 1:ceiling(i*up)){
flg = lo<j/i & j/i<up
if(flg==TRUE){
ans=paste0(j,'/',i)
break
}
}
i=i+1
}
cat(ans)
invisible(c(j,i))
}
fn()
fn(15/17,177/200)
# πの近似分数
fn(355/113,22/7)
- 237 ：卵の名無しさん：2020/09/27(日) 00:37:32.25 ID:ZaZXs1n5.net
- 電話番号0531454900の詳細情報 - 電話番号検索
http://web.archive.org/web/20191110142958/https://www.jpnumber.com/numberinfo_0531_45_4900.html
- 238 ：卵の名無しさん：2020/12/10(木) 21:17:39.06 ID:bHgC/hNU.net
- https://i.imgur.com/A2YWEIJ.png
VGB=integral(function(x) pi*green(x)^2, -r,b0) # volume green&blue
VR=integral(function(x) pi*red(x)^2, b0,a+R) # volume red
AG=integral(function(x) 2*pi*green(x), -r,b0) # area green
AB=integral(function(x) 2*pi*blue(x),b0,c+s) # area blue
AR=integral(function(x) 2*pi*red(x)^2, c+s,a+R) # area red
- 239 ：卵の名無しさん：2020/12/29(火) 00:31:19.09 ID:RkgAFAj8.net
- エセ医者がようやるわ。
- 240 ：卵の名無しさん：2020/12/29(火) 07:35:14.53 ID:v/EM6r6h.net
- >>239
では、同業者の意見を聞いてみましょう。
https://egg.5ch.net/test/read.cgi/hosp/1592662437/73
73 卵の名無しさん sage 2020/06/23(火) 13:24:47.79 ID:riQXI/fH
宮廷卒だけど、一括りに医師免許と言ってるが、私大卒など医者とは思ってへんよ
私大入学というインチキを経由したイシャモドキが、あんま調子のんなや
- 241 ：卵の名無しさん：2020/12/29(火) 09:55:53.10 ID:RkgAFAj8.net
- 自演w
- 242 ：卵の名無しさん：2020/12/29(火) 09:56:56.13 ID:RkgAFAj8.net
- 面白いねたしかに
AIみたいに定型文かコピペしか送ってこないんだもん
- 243 ：卵の名無しさん：2020/12/29(火) 10:03:17.22 ID:RkgAFAj8.net
- じゃあそんな君にこれをプレゼント！
うりゅうひろゆきの本スレにようこそ！
このウンコスレは基本的にさらしアゲて使用します！
ココは医師真似事務員ジイさんが、医師妬みの挙句、テメェで勝手に自称で医科歯科大を卒業した医師って設定で
妄想願望日記を発表するスレです
お医者さんのマトモなレスなど全く不要！
ゴミ箱まがいの公衆便所掃き溜め痰壷下品スレとして使用しましょうね
そんなわけで、いまだ懲りない恥さらし（自称）医科歯科事務員ジジイ登場～～
│ 彡川川川三三三ミ～ プウゥ～ン
｜ 川｜川／ ＼｜～ ポワ～ン ＿＿＿＿＿＿＿＿
｜ ∥｜∥ ◎---◎｜～ ／
｜ 川川∥ 3 ヽ～ ＜ 僕はうりゅうひろゆき
｜ 川川 ∴）д（∴）～ ＼＿＿＿＿＿＿＿＿
｜ 川川 ～ /～ カタカタカタ
｜ 川川∥ ～ /∥ ＿＿＿＿＿
｜ 川川川川＿＿＿／∥ | |￣￣＼ ＼
｜ ／ ＼＿＿| | |￣￣|
｜ / ＼医科歯科 | | |＿＿|
｜ | ＼ |つ |＿_|＿＿／ ／
｜ / 医師うらやましいぜ・・ |￣￣￣￣| 〔￣￣〕
よし、医科歯科の医者と称して
私立のお医者様を攻撃だ！なんたって、医科歯科の看板借りてんだからな
え？証明？？？自称だから医科歯科の学生証、卒業証書あるわけねえよ
同期の医師あげてみろ、ったて、そもそも医師じゃねえし居ねえよバカ
でも病院でもらった、医科歯科の封筒はもってま～す！これを証拠にしてやるぜ
専門医？欲しいけど・・いいや！そんなの必要ねえ！ 俺は事務員だし
- 244 ：卵の名無しさん：2021/01/28(木) 15:59:37.22 ID:RbXTChKG.net
- 某院でn人の患者感染。死亡率pとすると
予想される死亡者数とその95%信頼区間はいくらになるか？
"
p=0.20
n=100
(CI=HDInterval::hdi(qbinom,size=n,prob=p))
y=dbinom(0:n,n,p)
plot(0:n,y,type='h',bty='l')
sum(dbinom(CI[1]:CI[2],n,p))
- 245 ：卵の名無しさん：2021/02/01(月) 18:02:03.90 ID:e24XNIzt.net
- 「Googleドライブ」フォルダをDドライブにする
http://win10faq.blog.jp/archives/2208311.html
- 246 ：卵の名無しさん：2021/02/11(木) 08:40:29.92 ID:Sfff+JMX.net
- ウリュウのジジイ、5chでも社会でもまるで相手にされず！
- 247 ：卵の名無しさん：2021/05/03(月) 20:55:00.69 ID:ukTrFNeY.net
- 1 meter * 1 meter * 1 milli meter (water) = 1kg
1 mmHg柱/m2 = 13.6kg重/m2
13.6*9.8 N/m2 = 133.3 Pa =1.33hPa
1 hPa = 1/1.33 mmHg柱 = 0.752 mmHg柱
1 mmH20 = 1 kg重/m2
1*9.8 N/m2 = 9.8 Pa = 0.98 hPa
1 hPa = 1/0.98 = 1.02 mmH2O
- 248 ：卵の名無しさん：2021/05/17(月) 17:24:47.71 ID:vLTAoLk6.net
- https://news.mynavi.jp/article/20160317-a313/
about:config
privatebrowsing.autostartと入力して設定対象を絞り込み
- 249 ：卵の名無しさん：2021/06/16(水) 09:36:16.08 ID:Iaydx6qX.net
- ほんと庭なしバカ安達はくそ犬はうるせーわ道路でばかでかい声で会話するわやりたい放題
家族みんなデブだし
- 250 ：卵の名無しさん：2021/06/25(金) 06:15:51.74 ID:9PeLNEbV.net
- 【HKEY_CURRENT_USER > Control Panel > Personazation > Desktop Slideshow】
- 251 ：卵の名無しさん：2021/06/25(金) 17:47:33.56 ID:OGI6KuDR.net
- 尿瓶ジジイ＝トケジ＝ウリュウ＝ニセ内視鏡医 笑
- 252 ：卵の名無しさん：2021/07/03(土) 07:35:22.04 ID:C912OJj6.net
- tryCatch has a slightly complex syntax structure. However, once we understand the 4 parts which constitute a complete tryCatch call as shown below, it becomes easy to remember:
expr: [Required] R code(s) to be evaluated
error : [Optional] What should run if an error occured while evaluating the codes in expr
warning : [Optional] What should run if a warning occured while evaluating the codes in expr
finally : [Optional] What should run just before quitting the tryCatch call, irrespective of if expr ran successfully, with an error, or with a warning
tryCatch(
expr = {
# Your code...
# goes here...
# ...
},
error = function(e){
# (Optional)
# Do this if an error is caught...
},
warning = function(w){
# (Optional)
# Do this if an warning is caught...
},
finally = {
# (Optional)
# Do this at the end before quitting the tryCatch structure...
}
)
- 253 ：卵の名無しさん：2021/07/03(土) 07:36:30.21 ID:C912OJj6.net
- log_calculator <- function(x){
tryCatch(
expr = {
message(log(x))
message("Successfully executed the log(x) call.")
},
error = function(e){
message('Caught an error!')
print(e)
},
warning = function(w){
message('Caught an warning!')
print(w)
},
finally = {
message('All done, quitting.')
}
)
}
- 254 ：卵の名無しさん：2021/07/04(日) 06:05:02.42 ID:MvYVZ7Z5.net
- Windows Registry Editor Version 5.00
[HKEY_CURRENT_USER\Control Panel\Desktop\WindowMetrics]
- 255 ：卵の名無しさん：2021/07/06(火) 19:26:40.83 ID:7oMfliWM.net
- 古い低スペックPCに、Insider Programの正規版でWindows 11を入れてみました。
↓ Windows Insider Program公式サイト
https://insider.windows.com/ja-jp/
＜レジストリの変更＞
HKEY_LOCAL_MACHINE\SOFTWARE\Microsoft\WindowsSelfHost\UI\Selection
UIBranch ← Dev
UIContentType ←Mainline
UIRing ← External
HKEY_LOCAL_MACHINE\SOFTWARE\Microsoft\WindowsSelfHost\Applicability
BranchName ← Dev
ContentType ←Mainline
Ring ← External
HKEY_LOCAL_MACHINE\SYSTEM\Setup
新しいキー LabConfig を追加
LabConfigに以下のDWORD(32ビット値)を追加
BypassTPMCheck ← 1
BypassSecureBootCheck ← 1
Windows Updateでエラー0cx190012eが発生した場合、
C:\Users\Default\AppData\Local\Microsoft\Windows\WSUS
からsetupconfig.iniを削除します。
AppraiserRes.dllのダウンロードサイト
https://github.com/CodeProf14/Fix-TPM
C:\$WINDOWS.~BT\SourcesのDLLと置き換える
- 256 ：卵の名無しさん：2021/08/08(日) 22:23:50.74 ID:DBdOU9fb.net
- '%IN%' <- function(v,M) any(rowSums(M == v[col(M)]) == ncol(M))
c(4,5,6) %IN% matrix(c(1,2,3,
4,5,6,
7,8,9),
nrow=3,ncol=3,byrow=TRUE)
- 257 ：卵の名無しさん：2021/08/26(木) 18:32:55.73 ID:qkTk8DPf.net
- >>1
【全スレに注意連絡】
医療空間で連携して悪質な行為をする創.価.学.会.員.に注意してください。
【創.価.学.会.員.が連携して特定の患者をいじめたり、嫌がらせ、ガスライティング、不審な診察や治療、病院に受診できなくさせたり、たらい回し】
そのようなことが絶対にないように、みんなで不審な学.会.員の医療従事者の行動を見守りましょう。
そしておかしなことがあれば、そのまま放置せず、各所に報告や通報、ネットなどにも声をあげていきましょう。
隠蔽されやすい連携した嫌がらせは、そのように表沙汰にすることで弱体化していきます。
みんなで声をあげて悪質な連携を阻止していきましょう。
関連スレッド↓
医療法人アイクレセント いなだ小児科・アレルギー科 【佐賀】２
https://itest.5ch.net/egg/test/read.cgi/hosp/1624113085
集団ストーカー自称おおせき
https://lavender.5ch.net/test/read.cgi/minor/1622387124/l50
佐賀大学医学部附属病院
https://egg.5ch.net/test/read.cgi/hosp/1622027754/l50
情報交換などはこのスレでもできます↓
◆◆◆◆◆◆◆ 医療空間での創.価.学.会 ◆◆◆◆◆◆◆
https://egg.5ch.net/test/read.cgi/hosp/1628326388/l50
情報のリークはこのスレにも出来ます↓
誰にも言えない事を吐くスレ
https://egg.5ch.net/test/read.cgi/hosp/1627796592/l50 hdkdl
- 258 ：卵の名無しさん：2021/08/29(日) 23:14:57.19 ID:RCoITSyi.net
- [R] retrieving object names passed indirectly to a function
https://stat.ethz.ch/pipermail/r-help/2010-September/252840.html
foo <- function(p1, p2) {
p1name <- deparse(substitute(p1))
p2name <- deparse(substitute(p2))
# The real function runs a simulation, but
# here we just confirm the object names
# and values passed as args
cat("passed", p1name, "with value", p1, "\n")
cat("passed", p2name, "with value", p2, "\n")
}
- 259 ：卵の名無しさん：2021/10/30(土) 16:47:02.51 ID:k7cvh2/f.net
- "
レミフェンタニルの添付文書には，BMIが25を超える肥満患者には，標準体重を用いた投与量の調整を行うことが推奨されている．
高齢者:血圧低下などの副作用があらわれやすいため，開始用量の半減など，
患者の全身状態を観察しながら投与量を決定
"
# 2mg/20mL lwr=0.25μg/kg/min upr=0.05/kg/min ideal BMI=22
Ultiva<-function(cm,kg,age,lwr=0.25,upr=0.50,conc=2000/20,bmi=22){
BMI=kg/(cm/100)^2
LBM=ifelse(BMI>25,bmi*(cm/100)^2,kg) # lean body mass
L=lwr*LBM*60/conc
U=upr*LBM*60/conc
Lower=ifelse(age<75,L,L/2) # mL/h
Upper=ifelse(age<75,U,U/2)
round(c(BMI=BMI,LBM=LBM,'L(mL/h)'=Lower,'U(mL/h)'=Upper),2)
}
Ultiva(167.5,77.5,33)
Ultiva(137.2,50.8,87)
"
- 260 ：卵の名無しさん：2021/10/30(土) 16:47:07.75 ID:k7cvh2/f.net
- 持続投与量の目安は、セボフルランでは3-4 μg/kg/min、プロポフォールでは7.5 μg/kg/minと言われる。
添付文書:持続注入により投与する場合は、7μg/kg/分の投与速度で持続注入を開始する。
"
Eslax<-function(kg,sevo=FALSE,lwr=0.6,upr=0.9,conc=50e3/5){ # 50mg/5mL
precurarization=kg*0.03*1e3/conc
l=lwr*1e3*kg/conc # mL bolus
u=upr*1e3*kg/conc
L=ifelse(sevo,3*kg*60/conc,7*kg*60/conc) # on label
U=ifelse(sevo,4*kg*60/conc,7.5*kg*60/conc) # off label
cat('precararization(mL) =',precurarization,'\n')
cat('bolus(mL) =',l,'-',u,'\n')
cat('continuous(mL/h) =',L,'-',U,'\n')
}
Eslax(77.5,sevo=FALSE)
Eslax(50.8,sevo=TRUE)
- 261 ：卵の名無しさん：2021/10/31(日) 21:39:49.32 ID:asgKs/eN.net
- # Practically loading kg(mL/h) for 10 min and maintenance kg/10(ml/h)
Precedex<-function(kg,conc=200/50,lwr=0.2,upr=0.7,load=6){
I=kg*load/conc # 初期(10分)負荷(添付文書)
L=kg*lwr/conc # 維持下限(添付文書)
U=kg*upr/conc # 維持上限(添付文書)
i=kg*4/conc # 推奨初期負荷
m=kg*0.4/conc # 推奨維持
d=kg*0.1/conc # 譫妄予防
cat('sedation(mL/h) =',round(L,2),'-',round(U,2),'\n')
cat('delirium prevention(mL/h) =',round(d,2),'\n\n')
invisible(list(初期負荷=I,維持下限=L,維持上限=U,推奨初期負荷=i,推奨維持=m,譫妄予防=d))
}
Precedex(50)
> Precedex(50)
sedation(mL/h) = 2.5 - 8.75
delirium prevention(mL/h) = 1.25
- 262 ：卵の名無しさん：2021/10/31(日) 21:40:12.06 ID:asgKs/eN.net
- # pdfからcdfの逆関数を作ってhdiを表示させて逆関数を返す
# 両端での演算を回避 ∫[0,1]は∫[1/nxxx,1-1/nxx]で計算
pdf2hdi <- function(pdf,xMIN=0,xMAX=1,cred=0.95,Print=TRUE,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)
print(c(hdi[1],hdi[2]),digits=5)
if(Print){
par(mfrow=c(3,1))
plot(xx,sapply(xx,PDF),main='pdf',type='h',xlab='x',ylab='Density',col='lightgreen',bty='l')
legend('top',bty='n',legend=paste('HDI:',round(hdi,3)))
plot(xx,sapply(xx,cdf),main='cdf',type='h',xlab='x',ylab='Probability',col='lightblue',bty='l')
pp=seq(0,1,length=nxx)
pp=pp[-nxx]
pp=pp[-1]
plot(pp,sapply(pp,ICDF),type='l',xlab='p',ylab='x',main='ICDF',bty='l')
par(mfrow=c(1,1))
}
invisible(ICDF)
}
- 263 ：卵の名無しさん：2021/10/31(日) 21:41:36.08 ID:asgKs/eN.net
- # ベータ分布の差の分布
beta.diff<-function(r1,r2,n1,n2,a=1,b=1){
f <- function(x,y) dbeta(x+y,a+r1,b+n1-r1)*dbeta(y,a+r2,b+n2-r2)
vf=Vectorize(f,vectorize.args = 'y')
pdf <- function(x) integrate(function(y) vf(x,y),-1,1)$value
pdf=Vectorize(pdf)
invisible(pdf)
}
pdf=beta.diff(r1,r2,n1,n2)
- 264 ：卵の名無しさん：2021/11/03(水) 13:49:42.44 ID:WyXoRruB.net
- MAC2Sevo=function(MAC,age,N2O=0,MAC40sevo=1.8){
GOS=function(
age=85, # 年齢
FE_sevo=0.5, # 呼気セボフルラン濃度(%)
FE_N20=66.6){ # 呼気笑気濃度(%)
MAC40_sevo=MAC40sevo # 40歳セボフルランMAC(1.71,1.8,2.1と諸説あり)
MAC40_N2O=104 # 笑気MAC
MACage_sevo=MAC40_sevo*10^(-0.00269*(age-40))
MACage_N2O=MAC40_N2O*10^(-0.00269*(age-40))
MACage_total = FE_sevo/MACage_sevo + FE_N20/MACage_N2O
return(MACage_total) # 混合ガスのMAC
}
f=function(x) GOS(age,x,N2O) - MAC
f=Vectorize(f)
c('FEsevo(%)'=round(uniroot(f,c(0,5))$root,2))
}
MAC2Sevo(1.5,57)
- 265 ：卵の名無しさん：2021/11/04(木) 06:15:33.44 ID:6I+CHuG0.net
- # 年齢を設定してMACを笑気とセボフルレン濃度でMACをグラフ化
Age=87
GOS=function(
age=85, # 年齢
FE_sevo=0.5, # 呼気セボフルラン濃度(%)
FE_N20=66.6, # 呼気笑気濃度(%)
MAC40sevo=1.8){ # 40歳セボフルランMAC(1.71,1.8,2.1と諸説あり)
MAC40_sevo=MAC40sevo
MAC40_N2O=104 # 笑気MAC
MACage_sevo=MAC40_sevo*10^(-0.00269*(age-40))
MACage_N2O=MAC40_N2O*10^(-0.00269*(age-40))
MACage_total = FE_sevo/MACage_sevo + FE_N20/MACage_N2O
return(MACage_total) # 混合ガスのMAC
}
GOS=Vectorize(GOS)
sev=seq(0,3,by=0.1)
N2O=seq(0,70,by=1)
MAC=outer(sev,N2O,function(x,y)GOS(Age,x,y))
contour(sev,N2O,MAC,nlevels=50,
bty='l',xlab='Sevoflurane',ylab='nitrous oxide',main='MAC for age 87')
contour(sev,N2O,MAC,level=0.75,lwd=2,add=T)
- 266 ：卵の名無しさん：2022/04/06(水) 10:23:07.52 ID:HkRaCX9n.net
- f=\(n){
dat=c(0,rep(100,n-1))
m=mean(dat)
s=sd(dat)
50+(0-m)/s
}
nn=1000:5000
plot(nn,sapply(nn,f))
uniroot(f,c(1000,10000))
f(2501)
f(2502)
- 267 ：卵の名無しさん：2022/04/07(木) 06:36:23.98 ID:Ou/qCTfb.net
- calc = \(L,m,n){
f1 = \(x){
y = as.character(x)
strsplit(y,'') |> unlist() |> as.numeric()
}
f2 = \(x){
y = f1(x)
x%%n==0 & all(y==sort(y)) & sum(y)==m
}
x=((10^L-1)/9):(10^L-1)
x[sapply(x,f2)]
}
calc(3,16,2)
calc(5,24,3)
- 268 ：卵の名無しさん：2022/04/09(土) 06:30:01.59 ID:mZjVxAsM.net
- FF 83 F8 03 72
FF 83 F8 03 EB
- 269 ：卵の名無しさん：2022/04/09(土) 16:45:12.95 ID:wXi20OuD.net
- pH=7.32
pCO2=23
HCO3=16
pH_l=7.36
pH_u=7.44
pCO2_l=36
pCO2_u=44
HCO3_l=21
HCO3_u=27
if(pH<pH_l & pCO2<pCO2_l){
primary="metabolic acidosis"
expected_pCO2=40 + 1.2*(HCO3-24)
if(pCO2<expected_pCO2) additional= "secondary respiratory alkalosis")
if(pCO2>expected_pCO2) additional= "secondary respiratory acidosis")
}
if(pH>pH_u & pCO2>pCO2_u){
primary="metabolic alkalosis"
expected_pCO2=40 + 0.7*(HCO3-24)
if(pCO2<expected_pCO2) additional= "secondary respiratory alkalosis")
if(pCO2>expected_pCO2) additional= "secondary respiratory acidosis")
}
- 270 ：卵の名無しさん：2022/04/09(土) 16:45:22.36 ID:wXi20OuD.net
- if(pH<pH_l & pCO2>pCO2_u){
primary="respiratory aciosis"
status=NULL
if(HCO3_l<HCO3 & HCO3<HCO3_u) status="acute" else status="chronic"
if(status=="acute") expected_HCO3=24+0.1*(pCO2-40)
if(status=="chronic") expected_HCO3=24+0.4*(pCO2-40)
if(HCO3<expectedHCO3) additional="secondary metabolic acidosis")
if(HCO3>expectedHCO3) additional="secondary metabolic alkalosis")
}
if(pH>pH_u & pCO2<pCO2_l){
primary="respiratory alkalosis"
status=NULL
if(HCO3_l<HCO3 & HCO3<HCO3_u) status="acute" else status="chronic"
if(status=="acute") expected_HCO3=24+0.2*(40-pCO2)
if(status=="chronic") expected_HCO3=24+0.4*(40-pCO2)
if(HCO3<expectedHCO3) additional="secondary metabolic acidosis")
if(HCO3>expectedHCO3) additional="secondary metabolic alkalosis")
cat(paste(status,primary),'\n')
cat(additional)
- 271 ：卵の名無しさん：2022/04/09(土) 16:45:58.04 ID:wXi20OuD.net
- if(pH<pH_l & pCO2>pCO2_u){
primary="respiratory aciosis"
status=NULL
if(HCO3_l<HCO3 & HCO3<HCO3_u) status="acute" else status="chronic"
if(status=="acute") expected_HCO3=24+0.1*(pCO2-40)
if(status=="chronic") expected_HCO3=24+0.4*(pCO2-40)
if(HCO3<expectedHCO3) additional="secondary metabolic acidosis")
if(HCO3>expectedHCO3) additional="secondary metabolic alkalosis")
}
if(pH>pH_u & pCO2<pCO2_l){
primary="respiratory alkalosis"
status=NULL
if(HCO3_l<HCO3 & HCO3<HCO3_u) status="acute" else status="chronic"
if(status=="acute") expected_HCO3=24+0.2*(40-pCO2)
if(status=="chronic") expected_HCO3=24+0.4*(40-pCO2)
if(HCO3<expectedHCO3) additional="secondary metabolic acidosis")
if(HCO3>expectedHCO3) additional="secondary metabolic alkalosis")
cat(paste(status,primary),'\n')
cat(additional)
- 272 ：卵の名無しさん：2022/04/09(土) 17:37:32.78 ID:wXi20OuD.net
- Astrap <- \(pH,pCO2,HCO3){
pH_l=7.36
pH_u=7.44
pCO2_l=36
pCO2_u=44
HCO3_l=22
HCO3_u=26
if(pH_l<pH & pH<pH_u){
if(pCO2>pCO2_u){
primary="mixed respiratory acidosis"
additional="and metabolic alkalosis"
}
if(pCO2<pCO2_l){
primary="mixed respiratory alkalosis"
additional="and metabolic acidosis"
}
}
if(pCO2_l<pCO2 & pCO2<pCO2_u){
if(pH>pH_u){
primary="mixed metaboic alkalosis"
additional="and respiratory acidosis"
}
if(pH<pH_l){
primary="mixed metabolic acidosis"
additional="and respiratory alkalosis"
}
}
- 273 ：卵の名無しさん：2022/04/09(土) 17:38:10.02 ID:wXi20OuD.net
- if(pH<pH_l & pCO2<pCO2_l){
primary="primary metabolic acidosis"
expected_pCO2=40 + 1.2*(HCO3-24)
if(pCO2<expected_pCO2) additional= "secondary respiratory alkalosis"
if(pCO2>expected_pCO2) additional= "secondary respiratory acidosis"
}
if(pH>pH_u & pCO2>pCO2_u){
primary="primary metabolic alkalosis"
expected_pCO2=40 + 0.7*(HCO3-24)
if(pCO2<expected_pCO2) additional= "secondary respiratory alkalosis"
if(pCO2>expected_pCO2) additional= "secondary respiratory acidosis"
}
status=NULL
if(pH<pH_l & pCO2>pCO2_u){
primary="primary respiratory aciosis"
if(HCO3_l<HCO3 & HCO3<HCO3_u) status="acute" else status="chronic"
if(status=="acute") expected_HCO3=24+0.1*(pCO2-40)
if(status=="chronic") expected_HCO3=24+0.4*(pCO2-40)
if(HCO3<expected_HCO3) additional="secondary metabolic acidosis"
if(HCO3>expected_HCO3) additional="secondary metabolic alkalosis"
}
if(pH>pH_u & pCO2<pCO2_l){
primary="primary respiratory alkalosis"
if(HCO3_l<HCO3 & HCO3<HCO3_u) status="acute" else status="chronic"
if(status=="acute") expected_HCO3=24+0.2*(40-pCO2)
if(status=="chronic") expected_HCO3=24+0.4*(40-pCO2)
if(HCO3<expected_HCO3) additional="secondary metabolic acidosis"
if(HCO3>expected_HCO3) additional="secondary metabolic alkalosis"
}
cat(status,primary,'\n')
cat(additional)
}
- 274 ：卵の名無しさん：2022/04/09(土) 17:38:59.63 ID:wXi20OuD.net
- > Astrap(pH=7.23,pCO2=23,HCO3=16)
primary metabolic acidosis
secondary respiratory alkalosis
> Astrap(pH=7.54,pCO2=23,HCO3=17)
chronic primary respiratory alkalosis
secondary metabolic acidosis
> Astrap(pH=7.40, pCO2=23, HCO3=20)
mixed respiratory alkalosis
and metabolic acidosis
> Astrap(pH=7.40, pCO2=23, HCO3=35)
mixed respiratory alkalosis
and metabolic acidosis
> Astrap(pH=7.40, pCO2=23, HCO3=17)
mixed respiratory alkalosis
and metabolic acidosis
- 275 ：卵の名無しさん：2022/04/26(火) 16:51:59.52 ID:aUc6q9el.net
- >>1
【全スレ注意連絡】※告発内容はｺﾋﾟｰ自由ｶｸｻﾝ希望
【裏で集団ストーカー犯罪を指示してる医師、稲田成安に気をつけて】
http://www.inada-kids.com/wordpress/wp-content/uploads/greeting.jpg
佐賀県のいなだ小児科・アレルギー科 の小児科医、稲田成安（いなだ しげやす）に気をつけてください
自分に批判的な医療従事者や患者、その他の人たちを、探偵やソウカガッカイの集団ストーカーを使って調べ上げて、嫌がらせをやらせて仕返し、弾圧してきます
あらゆる嫌がらせで自殺や自滅に追い込む工作を仕掛けてきます
（※ ５ちゃんねるでは「ソウカガッカイ」という言葉を漢字四字で投稿するとアクセス禁止になる事があるのでカタカナ表記にしてあります）
もし稲田による集団ストーカー犯罪にあったら、証拠を取り、周囲や各機関に声を上げていきましょう。
そうする事で集団ストーカー組織は弱体化していきます
詳細は、↓のスレ
【ソウカ在日】稲田成安（いなだ しげやす）
https://egg.5ch.net/test/read.cgi/hosp/1644311065/
いなだ小児科・アレルギー科 【佐賀】
https://egg.5ch.net/test/read.cgi/hosp/1622365146/-100
黒猫世（くろねこよ）
https://rio2016.5ch.net/test/read.cgi/police/1648800949/l50 gyhjoooo
「 いなだ小児科・アレルギー科 集団ストーカー犯罪 」で検索
- 276 ：卵の名無しさん：2022/06/12(日) 19:47:04.91 ID:48n0Max4.net
- rm(list=ls())
library(RcppAlgos)
pm=permuteGeneral(v=1:3,freqs=c(3,3,4))
n=nrow(pm)
f=\(x){
t1=(1:10)[x==1]
t2=(1:10)[x==2]
t3=(1:10)[x==3]
if(t1[1]>t2[1]){
temp=t1
t1=t2
t2=temp
}
list(t1,t2,t3)
c(paste(LETTERS[t1],collapse=''),paste(LETTERS[t2],collapse=''),
paste(LETTERS[t3],collapse=''))
}
ans=unique(t(apply(pm,1,f)))
print(head(ans,15),quote=FALSE)
print(tail(ans,15),quote=FALSE)
- 277 ：卵の名無しさん：2022/06/12(日) 20:19:02.97 ID:48n0Max4.net
- "
10人の子供(A,B,C,....,J)を３チームに割り当てる。
各チームは１人以上の子供が割り当てられる。
何通りの分け方があるか？
例.
(A) (J) (B,C,D,E,F,G,H,I)
(A,C) (H,I,J) (B,D,E,F,G)
"
rm(list=ls())
gr=expand.grid(1:8,1:8)
gr1=gr[rowSums(gr)<10,]
f=\(x) 10-sum(x)
gr2=cbind(gr1,apply(gr1,1,f))
gr3=unique(t(apply(gr2,1,sort)))
gr3
ans=apply(gr3,1,\(x) choose(10,x[1])*choose(10-x[1],x[2])/ length(unique(x)) )
sum(ans)
ans
- 278 ：卵の名無しさん：2022/06/12(日) 20:25:08.35 ID:48n0Max4.net
- "10人の子供(A,B,C,....,J)を３人、３人、４人のチームに割り当てる。
何通りの分け方があるか、列挙せよ"
rm(list=ls())
library(RcppAlgos)
pm=permuteGeneral(v=1:3,freqs=c(3,3,4))
f=\(x){
t1=(1:10)[x==1]
t2=(1:10)[x==2]
t3=(1:10)[x==3]
if(t1[1]>t2[1]){
temp=t1
t1=t2
t2=temp
}
list(t1,t2,t3)
c(paste(LETTERS[t1],collapse=''),paste(LETTERS[t2],collapse=''),
paste(LETTERS[t3],collapse=''))
}
ans=unique(t(apply(pm,1,f)))
print(head(ans,20),quote=FALSE)
print(tail(ans,20),quote=FALSE)
- 279 ：卵の名無しさん：2022/06/12(日) 20:25:12.57 ID:48n0Max4.net
- "
10人の子供(A,B,C,....,J)を３チームに割り当てる。
各チームは１人以上の子供が割り当てられる。
何通りの分け方があるか？
例.
(A) (J) (B,C,D,E,F,G,H,I)
(A,C) (H,I,J) (B,D,E,F,G)
"
rm(list=ls())
gr=expand.grid(1:8,1:8)
gr1=gr[rowSums(gr)<10,]
f=\(x) 10-sum(x)
gr2=cbind(gr1,apply(gr1,1,f))
gr3=unique(t(apply(gr2,1,sort)))
ans=apply(gr3,1,\(x) choose(10,x[1])*choose(10-x[1],x[2])/ length(unique(x)) )
sum(ans)
res=cbind(gr3,ans)
colnames(res)=NULL
rownames(res)=NULL
res
- 280 ：卵の名無しさん：2022/08/15(月) 21:02:59.07 ID:WCPJCNZ4.net
- 多摩南部地域病院に入院した事あるけど看護師がクソ過ぎ
患者のことはどうでもよく自分達の仕事さえまわせればそれでいいって感じ
底意地の悪いクソ女看護師の清○桃○は要注意
そのうえ入院環境もクソ過ぎ
狭いトイレ、備品等もボロいものばかりで配慮が足りなさすぎる
病院自体のシステムも不親切だし二度と行きたくない
こんな病院は無いほうがいい
開業医もこんな病院紹介しないでほしい
- 281 ：卵の名無しさん：2022/10/07(金) 07:26:02.86 ID:i93G7UfB.net
- ganrikinto<-function(A=NULL,I=NULL,N=NULL,x=NULL){ # A:Principal I:Annual Interest N:months of repayment x:monthly repayment
if(is.null(I)){
I=uniroot(\(m) x-A*m/(1-(1+m)^(-N)),c(1e-12,1-1e-12))$root*12
return(c('Annual Interest'=I))
}
m=I/12
if(is.null(x)){
x=A*m/(1-(1+m)^(-N))
return(c('monthly repayment'=round(x,2)))
}
if(is.null(N)){
N=ifelse(x>A*m,log(x/(x-A*m))/log(m+1),Inf)
return(c('repayment period in months'=round(N,2)))
}
if(is.null(A)){
A=x/m*(1-(1+m)^(-N))
return(c('Principal'=round(A)))
}
}
- 282 ：卵の名無しさん：2024/01/04(木) 05:47:47.73 ID:s7ncVxDB.net
- f=\(n){
n1=n%%10
i=1
flg <- (n1*n1^i)%%10==n1
if(flg){ p=i
}else{
while(!flg & i<2024){
i=i+1
flg <- (n1*n1^i)%%10==n1
}
p=i
}
j=ifelse(2024%%p,2024%%p,p)
(n1^p)%%10
}
sapply(1:2024,f) |> table()
- 283 ：卵の名無しさん：2024/01/04(木) 11:13:27.39 ID:mx4lFG6x.net
- 2024＾2024
π^2024
の先頭の数字を算出
臨床医の嗜み：R言語での一般解を算出するコード
calc=function(m,n) floor(10^(n*log10(m) -floor(n*log10(m))))
calc(2024,2024)
calc(pi,2024)
- 284 ：卵の名無しさん：2024/01/05(金) 07:22:34.13 ID:lKTSD30o.net
- "
地球の公転周期は秒単位だと365日5時間48分46秒であるという
5時間48分46秒＝20926秒
閏年は400年に97回
1日＝86400秒
"
p=20926/86400 ; p
97/400 - p
8/33 - p
fn=\(x){
xp=x*p
fl=floor(xp)
ce=ceiling(xp)
dfl=abs(fl/x-p)
dce=abs(ce/x-p)
c(min(dfl,dce),ifelse(dfl<dce,fl,ce),x)
}
dat=t(sapply(1:40,fn))
rbind(dat[which.min(dat[,1]),])
dat=t(sapply(1:10000,fn))
rbind(dat[which.min(dat[,1]),])
n=seq(4,400,4)
dat=t(sapply(n,fn))
rbind(dat[which.min(dat[,1]),])
- 285 ：卵の名無しさん：2024/01/05(金) 19:36:46.56 ID:Nm2/P+r6.net
- "
短径20cm 長径24cmの楕円のピザを楕円弧の長さが等しいように３つの楕円扇形に分割する。
楕円扇形の面積が最大になるようにするにはどのように切ればよいか？
そのときの面積はいくつか？（答は整数値でよい）。
あらゆるリソースを用いてよい。
"
a=12
b=10
n=3
S=\(t1,t2) a*b/2*(t2-t1)
L=\(t1,t2) integrate(\(x) sqrt((-a*sin(x))^2+(b*cos(x))^2),t1,t2,abs.tol = 1e-24)$value
Lt=L(-pi,pi) ; Lt # 全周長
t12t2=\(t1) uniroot(\(x) L(t1,x)-Lt/n,c(t1,t1+pi),tol=1e-24)$root
t12S=\(t1) S(t1,t12t2(t1))
t12S=Vectorize(t12S)
curve(t12S,0,pi)
opt=optimize(t12S,c(-pi/2,pi/2),maximum = TRUE,tol=1e-24) ; opt
source('toolmini.R')
Plot(-a,a,zero=FALSE)
Ellip(0i,a,b,lty=3)
t1=opt$maximum ; t2=t12t2(t1)
tt=seq(t1,t2,length=100)
lines(a*cos(tt),b*sin(tt),lwd=2)
T1=a*cos(t1)+1i*b*sin(t1)
T2=a*cos(t2)+1i*b*sin(t2)
seg(0i,T1,lty=2,lwd=2)
seg(0i,T2,lty=2,lwd=2)
Angle(T1,0i,T2)
(pi*a*b - opt$obj)/2
S(t2,pi)
- 286 ：卵の名無しさん：2024/01/06(土) 16:06:08.99 ID:YqxPaf+q.net
- rm(list=ls())
par(bty='l')
calc=function(m=2024,n=2024) floor(10^(n*log10(m) -floor(n*log10(m))))
calc(2024,2024)
calc(pi,2024)
calc(exp(1),2024)
"
nを1以上2024以下の整数とする。
n^2024の先頭に現れる数字はなにが一番多いか？
頻度順に並べよ。
"
calc.top =\(n) sapply(1:n,calc) |> table() |> sort() |> rev()
re=calc.top(2024)
re
p=re/sum(re)
calc.top(1000)
sapply(1:1000,\(x) as.numeric(names(calc.top(x)[1])))
k=1:9
plot(k,p,pch=16)
lines(k,log10(1+1/k),type='h',col=2) # ベンフォードの法則
- 287 ：卵の名無しさん：2024/01/12(金) 16:58:50.71 ID:NVcGRVeI.net
- 医者板では全く相手にされないから数学板で高校生相手にイキれると思ってる尿瓶チンパンジジイ滑稽だね
- 288 ：卵の名無しさん：2024/04/08(月) 05:59:31.45 ID:0nrZk2Or.net
- これが現実
１次方程式もできないド底辺特殊シリツ医大卒の記録
http://imagizer.imag...g923/2715/RosCsf.jpg
何度読んでも馬鹿すぎる。
男女別の割合と全体での割合から男女比が計算できるとも思わないとは。
なんでこんなのが大学に入れるわけよ？
裏口入学以外に説明がつく？
中学生でも解ける一次方程式の問題だろ。
それすらできない馬鹿が自信を持って発言。
>患者の男女比が必要なのもわからないのか？
だとさ。
http://imagizer.imag...g923/9687/zNivZW.jpg
145 KB
新着レスの表示
掲示板に戻る
全部
前100
次100
最新50
read.cgi ver.24052200