臨床統計もおもしろいですよ、その3

**卵の名無しさん** · 2020/03/05(木) 20:17:05.03

　
　内科認定医受験の最低限の知識、
　製薬会社の示してくる臨床データ、
　論文の考察、
　論文を書くときの正当性、
　というのが、臨床統計の今までの目的の大きい部分でしたが、
　
　AI=機械学習の基本も、結局は統計学と確率に支配されます。
　そういう雑多な話をするスレです。
　
※前スレ
臨床統計もおもしろいですよ、その１
https://egg.2ch.net/test/read.cgi/hosp/1493809494/
臨床統計もおもしろいですよ、その2
https://egg.5ch.net/test/read.cgi/hosp/1540905566/

**卵の名無しさん** · 2020/05/11(月) 22:43:11.83

>>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)

**卵の名無しさん** · 2020/05/11(月) 22:43:15.93

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')

**卵の名無しさん** · 2020/05/11(月) 22:44:15.44

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()

**卵の名無しさん** · 2020/05/11(月) 22:45:23.11

reshape2 が　tidyr::pivot_longer(df4,col=1:10)　になっていた。
使っていないと、ggplot2も忘れてしまうなぁ。

**卵の名無しさん** · 2020/05/12(火) 07:14:59.85

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

**卵の名無しさん** · 2020/05/12(火) 07:58:45.50

えぇぇえぇぇぇっ

Ｆラン統計ジジイ　あらため　コロナ薬屋さん
および　コンプ薬屋さんって
裏口容疑者なんですかぁぁあ？

だと思いましたぁぁあぁぁ
早く出頭するんですよ
あなたたちのレスの全てから
頭の悪さを感じますからねえぇぇぇ

by 都内Ｓラン女子高生

**卵の名無しさん** · 2020/05/12(火) 08:31:27.05

交絡因子を考慮しての記憶減衰の階層ベイズモデル
デバッグできたので、事前分布をいろいろ変えてみて再現性を確認。
内視鏡休診で纏まった勉強する時間が取れて( ・∀・)ｲｲ!!

##### 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')

**卵の名無しさん** · 2020/05/12(火) 12:01:40.73

開業医ではないあなたが
先に開業医板を荒らすから
私がここに降臨して差し上げてるんですよ

それすら判らないくらい知能が低いんですかねぇ

女子高生にブチ切れる異常者のジジイって
どんなホラー映画の地獄絵図ですかＷ

by 都内Ｓラン女子高生

**卵の名無しさん** · 2020/05/12(火) 17:39:16.61

>>143
half cauchy にしても一様分布にしても事後分布はほとんど変化はないな。
モデルが優れているからだろう。

**卵の名無しさん** · 2020/05/12(火) 18:12:20.47

さて今日も学校は終わり。
お家に帰ってお勉強しましょう。

あ～あ、Ｆラン統計ジジイは
やっぱりくるくるぱーだったんですねぇ

by 都内Ｓラン女子高生

**卵の名無しさん** · 2020/05/12(火) 23:05:47.83

# 相関係数の差の検定
n1=150
r1=0.29
n2=120
r2=0.5
pnorm(-abs(atanh(r1)-atanh(r2))/sqrt(1/(n1-3)+1/(n2-3)))

**卵の名無しさん** · 2020/05/13(水) 05:47:04.95

離散量のプロットは重なってしまうと相関がわかりにくくなる。
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))

**卵の名無しさん** · 2020/05/13(水) 06:15:19.70

色の濃淡に頻度を反映させるという方法もあるな。

https://i.imgur.com/TnXbBKY.png

**卵の名無しさん** · 2020/05/13(水) 06:16:02.77

とりあえず、完成。多分、もっと便利なパッケージがあるんだろうな。

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)

**卵の名無しさん** · 2020/05/13(水) 12:26:32.87

学校の先生に習いましたが、統計的にはPCR検査はむやみにやらないのが正解なんですね

by 都内Ｓラン女子高生

**卵の名無しさん** · 2020/05/13(水) 14:17:48.98

【国民が誤解】安倍晋三「相談・受診の目安＝PCR実施の基準ではない。国民への周知が足りなかった」 [スタス★]
https://asahi.5ch.net/test/read.cgi/newsplus/1589284722/

**卵の名無しさん** · 2020/05/13(水) 14:18:56.03

世界の頭脳は実態把握には大量検査が必要とわかっているなぁ。
東大の児玉名誉教授も同じ意見だな。

https://www.covid19-yamanaka.com/cont4/14.html

**卵の名無しさん** · 2020/05/13(水) 14:55:45.98

>>150
色を変えて散布図を作成

https://i.imgur.com/l1eyD6o.png

**卵の名無しさん** · 2020/05/13(水) 14:57:54.70

検診がメインだったから、内視鏡バイトは休診。
俺は６割の休業補償だが、銭ゲバ開業医はパート医を首にしているんだろうなぁ。

**卵の名無しさん** · 2020/05/13(水) 15:27:10.62

今日も朝から人生が暇しかない
Ｆラン統計ジジイ　あらため　コロナ薬屋　が
発狂して連投だぁぁぁW

ちゃんとコテハン付けないとダメでちよぉぉぉ

今思うのは反日野党が政権をとってたら
日本はとっくに壊滅していたであろうこと
だって反日なんだもの

by 都内Ｓラン女子高生

**卵の名無しさん** · 2020/05/13(水) 15:48:20.06

>>156
内視鏡バイトは休診なので、６割の休業補償はあるけど。
あんたのところは、感染拡大阻止に休診にして職員には休業補償してんの？

まとまった時間がとれるので、https://bayesmodels.com/の英文テキストで独学中。
行き詰まったら数学版で相談。

**卵の名無しさん** · 2020/05/13(水) 18:18:49.08

相関係数->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

**卵の名無しさん** · 2020/05/14(木) 09:44:43.41

まあ、ベイズファクターの計算結果はさほど変わらんな。

> 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

**卵の名無しさん** · 2020/05/14(木) 09:45:47.27

>>159
グラフにしてみた。

https://i.imgur.com/OyTi5Jb.png

**卵の名無しさん** · 2020/05/14(木) 10:03:25.31

あんまり相関係数の信頼区間など気にしていなかった。
必要なら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変換の逆関数で出しているな。

**卵の名無しさん** · 2020/05/14(木) 10:39:23.42

>>161
sample関数を使ってのbootstrap

https://i.imgur.com/gAIJTEQ.png

**卵の名無しさん** · 2020/05/14(木) 11:21:55.13

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)

**卵の名無しさん** · 2020/05/14(木) 11:33:51.60

>>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

**卵の名無しさん** · 2020/05/14(木) 11:37:20.35

>>156
アベノミクスでの経済成長率は震災があった民主の頃より低いという現実を認められないんだなぁ。

日本を破壊したいなら、安倍一択。
安倍の長期政権で中国も韓国も笑いが止まらないはず、
日本が衰退する政策ばかりを選択しているから。

野党が反感をもたれて安倍信者が増えることこそ日本破壊が進んで
パヨクの思う壺

**卵の名無しさん** · 2020/05/17(日) 15:29:50.58

結局、これが核心部分

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

**卵の名無しさん** · 2020/05/18(月) 07:52:43.74

# 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)

**卵の名無しさん** · 2020/05/18(月) 13:03:08.05

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)

**卵の名無しさん** · 2020/05/18(月) 22:00:40.40

今日も朝から人生が暇しかない
Ｆラン統計ジジイ　あらため　コロナ薬屋　が
発狂して連投だぁぁぁW

ちゃんとコテハン付けないとダメでちよぉぉぉ

今思うのは反日野党が政権をとってたら
日本はとっくに壊滅していたであろうこと
だって反日なんだもの

by 都内Ｓラン女子高生

**卵の名無しさん** · 2020/05/19(火) 07:38:54.76

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)

**卵の名無しさん** · 2020/05/19(火) 07:59:02.34

今日も朝から人生が暇しかない
Ｆラン統計ジジイ　あらため　コロナ薬屋　が
発狂して連投だぁぁぁW

相変わらず臨床経験や社会経験がゼロなのが
丸わかりのレスですねぇぇぇぇぇ

ちゃんとコテハン付けないとダメでちよぉぉぉ

今思うのは反日野党が政権をとってたら
日本はとっくに壊滅していたであろうこと
だって反日なんだもの

by 都内Ｓラン女子高生

**卵の名無しさん** · 2020/05/21(木) 12:57:01.59

今日のくるくるぱーのIDだああぁぁあ
id:3+Rla4zY

くるくるぱーが反日野党に肩入れして
自滅するのは自業自得だけど
善良な日本人は巻き込まれたくないですねえ
by 都内Ｓラン女子高生

**卵の名無しさん** · 2020/05/21(木) 17:33:16.75

今日も朝から人生が暇しかない
Ｆラン統計ジジイ　あらため　コロナ薬屋　が
発狂して連投だぁぁぁW

相変わらず臨床経験や社会経験がゼロなのが
丸わかりのレスですねぇぇぇぇぇ

今思うのは反日野党が政権をとってたら
日本はとっくに壊滅していたであろうこと
だって反日なんだもの

by 都内Ｓラン女子高生

**卵の名無しさん** · 2020/05/22(金) 17:43:59.44

すいません、すごい基本的なことなのですが教えてください。

現在いる実験室で、様々細胞と遺伝子変異させた細胞に抗がん剤をかけて死亡率をみてます。
抗がん剤の量も比較してるのです
1: 細胞aとbに対して同量の抗がん剤を使用し死亡率を見る場合
2: 同種の細胞に違う濃度の抗がん剤をかけて死亡率を比較する場合。
両者とも何検体かとって平均の比較をする場合、対応するt検定ですか？、それともf値をみた上で対応しないt検定になるのですか？

**卵の名無しさん** · 2020/05/22(金) 19:39:38.80

>>174
比率の検定ならΧ（カイ）二乗検定。

**卵の名無しさん** · 2020/05/23(土) 07:42:32.91

今日も朝から人生が暇しかない
Ｆラン統計ジジイ　あらため　コロナ薬屋　が
発狂して連投だぁぁぁW

身の程を弁えない謎の上から目線で政治談批評
老害アルアルですねぇぇぇぇぇ

今思うのは反日野党が政権をとってたら
日本はとっくに壊滅していたであろうこと
だって反日なんだもの

by 都内Ｓラン女子高生

**卵の名無しさん** · 2020/05/23(土) 07:57:58.92

日本衰退を願う勢力は安倍一択。
観光立国という亡国政策で衰退途上国。
観光でしか生きていけない国になったときに中国からの出国禁止すれば日本は中国に平伏すしかなくなる。

**卵の名無しさん** · 2020/05/23(土) 07:58:44.46

>>174
対数変換して検定した方がいい場合もある。

**卵の名無しさん** · 2020/05/23(土) 22:14:18.97

テレビを見てるんですが
やっぱりマスコミはクズですね
問題のあった新聞社の社長は直ちに謝罪して
新聞社を解体するべきですね

反日野党は日本の足を引っ張ることしかしてませんね
存在意義がないとしか言えません
だって反日なんだもの

by 都内Ｓラン女子高生

**卵の名無しさん** · 2020/05/24(日) 15:16:08.53

>>178
対数変換した場合、検定方法は何使えばいいですか？

**卵の名無しさん** · 2020/05/24(日) 18:15:55.34

>>180
平均の検定ならt検定

**卵の名無しさん** · 2020/05/25(月) 12:16:59.94

ある人物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')

**卵の名無しさん** · 2020/05/25(月) 12:21:54.11

>>181
t検定は二つのことを比べることしかできないですよね？

**卵の名無しさん** · 2020/05/25(月) 16:55:34.95

>>183
多重検定やりたいならTukeyとか色々ある。

**卵の名無しさん** · 2020/05/25(月) 17:16:47.59

>>184
図の点線が対照のmockで実線が遺伝子変化させた細胞、
横軸が加えた薬剤の濃度です。
ルシフェラーゼの比を見ています
同じ濃度で、モックと遺伝子変化させた群を比較させたい時は統計は何使えばいいですか？
https://i.imgur.com/VXfJ8TD.jpg

**卵の名無しさん** · 2020/05/25(月) 17:49:48.35

好きなの使えば。
厳しいbonferroni補正したpairwise.t.testとか

**卵の名無しさん** · 2020/05/25(月) 18:00:37.69

回帰係数の比較ならコクランアーミテッジだったかな。
Rだとprop.trend.testでできたはず！

**卵の名無しさん** · 2020/05/25(月) 19:07:25.19

皆さんありがとうございます
すいませんjmpかエクセル統計で計算できるのだと嬉しいです。

**卵の名無しさん** · 2020/05/25(月) 20:15:07.84

やったね　たえちゃん

**卵の名無しさん** · 2020/05/25(月) 21:59:58.02

Rだと
kruskal.testで多重比較で有意差確認してから
平均比較ならpairwise.t.test
比率比較ならpairwise.prop.test
補正法は指定できる。

エクセルのマクロは売られている。
使ったこともないけど。
https://bellcurve.jp/ex/function/kruskal.html

**卵の名無しさん** · 2020/05/26(火) 07:48:20.04

今日も朝から人生が暇しかない
Ｆラン統計ジジイ　あらため　コロナ薬屋　が
発狂して連投だぁぁぁW

それにしても反日野党や反日マスコミは日本の足を引っ張ることしかしてませんね
存在意義がないとしか言えません
だって反日なんだもの

by 都内Ｓラン女子高生

**卵の名無しさん** · 2020/05/26(火) 12:50:37.25

"
>勤務医が院長のおれの悪口を妻に話してたそうだ
>腹が立ったからクビにする
を登場人物としてみた。
(問題)
新型コロナに勤務医が罹患。
勤務医が発症した翌日に院長が発症、その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)

**卵の名無しさん** · 2020/05/26(火) 12:50:43.88

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))

**卵の名無しさん** · 2020/05/26(火) 23:37:54.50

まず全裸になり
　　　　　　　　　　　　　( 　：　)
　　　　　　　（゜∀゜）ﾉ彡
　　　　　　 <(　　　)
　　　　　　　　ﾉωヽ

　自分の尻を両手でバンバン叩きながら白目をむき
　　　　　　　　　　　从
　　　　　　　Дﾟ　）　　て
　　　　　　　 (　ヾ) )ヾ　て
　　　　　　　 < <

　　　　　人__人__人__人__人__人__人__人__人__人__人
　　　　Σ　　　　　　　　　　　　　　　　　　　　　　　　　て
　　　　Σ　　びっくりするほどユートピア！　　　　　　　て人__人_
　　　　Σ　　　　　　　　　びっくりするほどユートピア！　　　　　　て
　　　　　⌒Ｙ⌒Ｙ⌒Ｙ)　　　　　　　　　　　　　　　　　　　　　　　　　て
　　　　　　　　　　　　Ｙ⌒Ｙ⌒Ｙ⌒Ｙ⌒Ｙ⌒Ｙ⌒Ｙ⌒Ｙ⌒Ｙ⌒Ｙ⌒Ｙ⌒
　＿_＿＿＿＿＿
　|＿＿　　　　　　ヽ(゜∀゜）ﾉ
　|\＿〃´￣￣ヽ..ヘ(　　 )ミ
　|　|\,.-～´￣￣　　 ω　>　（∀゜）ﾉ
　\|∫\　　 _,. - ､_,. - ､ \　(　　ﾍ)
　　 \ 　 \＿＿＿＿＿＿ _\<
　　　　\ 　||￣￣￣￣￣￣￣ |
　　　　　\||＿＿＿＿＿＿＿ |

これを10分程続けると妙な脱力感に襲われ、解脱気分に浸れる。

**卵の名無しさん** · 2020/05/27(水) 04:43:18.93

x1,x2,...,xnの順に発症
その間隔はt1,t2,..,tn-1としたときにxiが感染源であった確率を算出するプログラムを作れ。

**卵の名無しさん** · 2020/05/27(水) 05:15:26.64

安倍は新コロナは中国発と認めたから春節ウェルカムは外患誘致。
早くアメリカと足並み揃えて中国に損害賠償請求して国民に赦しを乞うべき。

観光立国という亡国政策で日本は衰退途上国。
次の世代の日本人は実習生として中国に出稼ぎがで立場が逆転。
日本人なら過労死自殺しても経営者殺害はしないから酷使される。

観光でしか生きていけない国になったときに中国からの訪日禁止すれば日本は中国に平伏すしかなくなる。
尖閣どころか沖縄を寄越せとか言われても差し出すことになりそうだな。

**卵の名無しさん** · 2020/05/27(水) 06:40:49.96

>>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人発症

**卵の名無しさん** · 2020/05/27(水) 07:57:44.10

やったね　たえちゃん

Ｆラン統計ジジイ　あらため　コロナ薬屋は
今日もクルクルパーだよ

**卵の名無しさん** · 2020/05/27(水) 17:40:48.61

>>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];
}
}

'

**卵の名無しさん** · 2020/05/27(水) 17:40:55.22

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))

**卵の名無しさん** · 2020/05/27(水) 17:41:13.27

> 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

**卵の名無しさん** · 2020/05/27(水) 17:41:42.49

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))

**卵の名無しさん** · 2020/05/27(水) 17:44:39.93

> 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の値を採用しよう。

**卵の名無しさん** · 2020/05/27(水) 17:46:07.20

発症した順が必ずしも感染順ではないので、
職員が発症しても責めたりせずに悪いのは春節ウェルカムした安倍のせいと心をまとめよう。

**卵の名無しさん** · 2020/05/28(木) 17:57:18.03

# 武漢大量検査
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]

**卵の名無しさん** · 2020/05/28(木) 17:58:09.17

> 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

**卵の名無しさん** · 2020/05/28(木) 18:02:34.43

検査は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

**卵の名無しさん** · 2020/05/30(土) 10:34:55.69

>>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)

**卵の名無しさん** · 2020/05/30(土) 10:50:01.68

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)

**卵の名無しさん** · 2020/06/01(月) 17:57:09.41

今日のID:nu4+jSn4

**卵の名無しさん** · 2020/06/01(月) 20:31:14.88

>>210
ゴルゴ問題の答が出せるように勉強したのか？

**卵の名無しさん** · 2020/06/26(金) 17:08:21.56

"
拳銃を一発撃ったときに、
狙った相手を撃ち殺す確率は、
Ａは 1/3、Ｂは 1/2（50%）、Ｃは 1/1（100%）とします。
なお、この確率は、全員が知っているものとします。
拳銃を撃つ順番は、Ａ、Ｂ、Ｃの順番で、
以降は最後の一人が生き残るまでこの順番を繰り返すものとします。
Ａ、Ｂ、Ｃは拳銃を撃つときに誰を狙っても良いこととします。
ただし、一発で二人を狙うことはできません。
Ａ、Ｂ、Ｃの生き残る確率を求めなさい。
"

**卵の名無しさん** · 2020/06/26(金) 17:08:35.96

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))
}

**卵の名無しさん** · 2020/06/26(金) 17:08:51.89

# 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)

**卵の名無しさん** · 2020/06/28(日) 08:05:03.04

# 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)

**卵の名無しさん** · 2020/06/28(日) 08:05:24.51

# 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)

**卵の名無しさん** · 2020/07/05(日) 18:42:40.47

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

**卵の名無しさん** · 2020/07/05(日) 18:43:03.70

> 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

**卵の名無しさん** · 2020/07/05(日) 20:04:10.52

デバッグ版

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

**卵の名無しさん** · 2020/07/05(日) 20:04:14.42

まあ、パーティー中に感染させる人数に再生産数を流用していいかは疑問ではあるが、

実効結果は

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

**卵の名無しさん** · 2020/07/07(火) 16:20:43.89

# 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)

**卵の名無しさん** · 2020/07/09(木) 21:20:16.36

library(randtests)
DEL <- permut(1:4,3,function(x) print(paste0(x,collapse = ''),quote=F))

**卵の名無しさん** · 2020/07/10(金) 23:07:16.57

まじで、エクセル、できればLibreOffice　で出来る
医療統計の本を出して欲しい、この先生なら出来る。

**卵の名無しさん** · 2020/07/10(金) 23:13:07.04

自分で、この先生とか言うなよ
恥ずかしい奴だな

**卵の名無しさん** · 2020/07/13(月) 21:47:22.50

リアル223からです、数式を記載されている方とは別人です。

**卵の名無しさん** · 2020/07/14(火) 17:33:19.37

>>223
エクセルのソルバーなんて使えたものじゃないね。
エクセルはグラフも貧弱だし。
医療統計の本なら朝倉から丹波某が沢山書いた本が沢山でているよ。
ベイズはWINBUGの話で古いけど。

**卵の名無しさん** · 2020/07/18(土) 23:13:50.16

"
容量8Lの袋と容量5Lの袋を使って池の水を丁度4L集めたい。
袋に目盛りはついていません。
袋から袋への移し替えは全量で行います。
池からとる水の量や池に捨てる水の量には制限はありません。
最初に片方に満たした作業を１回目として
丁度4Lを集めるのに最低何回の移動が必要か？
"

**卵の名無しさん** · 2020/07/18(土) 23:14:05.95

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)

**卵の名無しさん** · 2020/07/18(土) 23:14:23.30

# 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)
}

**卵の名無しさん** · 2020/07/26(日) 10:46:35.07

'%&%' <- 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

**卵の名無しさん** · 2020/08/04(火) 21:12:57.26

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)
}

**卵の名無しさん** · 2020/08/05(水) 16:36:41.12

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

**卵の名無しさん** · 2020/08/06(木) 07:04:19.79

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/

**卵の名無しさん** · 2020/08/09(日) 22:02:52.31

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

**卵の名無しさん** · 2020/08/17(月) 22:48:07.01

#原子数=放射能/崩壊定数=放射能/(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

**卵の名無しさん** · 2020/08/22(土) 07:22:17.76

# 問題「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)